Solved: Problem in using MKL Trust Region dtrnlsp solving system of nonlinear equations

peng_z_ · ‎12-18-2016

Hello,

Currently, I am using the dtrnlsp_solve (trust region )to sovle a systems of 4 nonlinear equations. I am following the example provided in the https://software.intel.com/en-us/node/471544 Howerver , the powell method failured to get solution.

When I resorted to the fsolve function in Matlab, which uses the trust region method by default as well (pls refer to the link https://www.mathworks.com/help/optim/ug/fsolve.html, I can sucessful get the root of the equations.

Can you pls offer me some suggestions on how to solve it using Dtrnlsp_solver? Thank you in advance。

the matlab code is shown below:

function F = EQs(x)
k=190000/3/(1-2*0.3);
g=190000/(1+0.3)/2;
pe=-167.6485 ;
qe=500.0385  ;
p=pe+k*x(1);
q=qe-3*g*x(2);
peeq=x(3);
sm=500+1000*(peeq);
freal=x(4);
ds=sinh(-1.5*p/sm);
dc=cosh(-1.5*p/sm);
a=0.04*exp(-0.5*((peeq-0.3)/0.1)^2)/((2*3.14159)^0.5)/0.1;

F(1) = x(1)* 2*q/sm/sm+x(2)*(-3*1.5*freal*ds/sm);
F(2) = (q/sm)^2+2*freal*1.5*dc-(1+2.25*freal^2);
F(3)=x(3)-(-p*x(1)+q*x(2))/(1-freal)/sm
F(4)=x(4)-(1-freal)*x(1)-a*peeq

 >>options = optimoptions('fsolve','Display','iter','PlotFcn',@optimplotfirstorderopt);

>>x0=0;

>>[x,fval]  = fsolve(@EQs,x0,options)

the fortran code is shown below :

!===============================================================================
! Copyright 2004-2016 Intel Corporation All Rights Reserved.
!
! The source code,  information  and material  ("Material") contained  herein is
! owned by Intel Corporation or its  suppliers or licensors,  and  title to such
! Material remains with Intel  Corporation or its  suppliers or  licensors.  The
! Material  contains  proprietary  information  of  Intel or  its suppliers  and
! licensors.  The Material is protected by  worldwide copyright  laws and treaty
! provisions.  No part  of  the  Material   may  be  used,  copied,  reproduced,
! modified, published,  uploaded, posted, transmitted,  distributed or disclosed
! in any way without Intel's prior express written permission.  No license under
! any patent,  copyright or other  intellectual property rights  in the Material
! is granted to  or  conferred  upon  you,  either   expressly,  by implication,
! inducement,  estoppel  or  otherwise.  Any  license   under such  intellectual
! property rights must be express and approved by Intel in writing.
!
! Unless otherwise agreed by Intel in writing,  you may not remove or alter this
! notice or  any  other  notice   embedded  in  Materials  by  Intel  or Intel's
! suppliers or licensors in any way.
!===============================================================================

*   Content : TR Solver Fortran-77 example
*
********************************************************************************

!** NONLINEAR LEAST SQUARE PROBLEM WITHOUT BOUNDARY CONSTRAINTS
      PROGRAM EXAMPLE_DTRNLSP_POWELL
        IMPLICIT NONE
!** HEADER-FILE WITH DEFINITIONS (CONSTANTS, EXTERNALS)
        INCLUDE "mkl_rci.fi"
!** USER'S OBJECTIVE FUNCTION
        EXTERNAL            EXTENDED_POWELL
!** N - NUMBER OF FUNCTION VARIABLES
        INTEGER             N
        PARAMETER           (N = 4)
!** M - DIMENSION OF FUNCTION VALUE
        INTEGER             M
        PARAMETER           (M = 4)
        
         DOUBLE PRECISION    pe,qe,k,g  
         common pe,qe,k,g
!** SOLUTION VECTOR. CONTAINS VALUES X FOR F(X)
        DOUBLE PRECISION    X (N)
!** PRECISIONS FOR STOP-CRITERIA (SEE MANUAL FOR MORE DETAILS)
        DOUBLE PRECISION    EPS (6)
!** JACOBI CALCULATION PRECISION
        DOUBLE PRECISION    JAC_EPS
!** REVERSE COMMUNICATION INTERFACE PARAMETER
        INTEGER             RCI_REQUEST
!** FUNCTION (F(X)) VALUE VECTOR
        DOUBLE PRECISION    FVEC (M)
!** JACOBI MATRIX
        DOUBLE PRECISION    FJAC (M, N)
!** NUMBER OF ITERATIONS
        INTEGER             ITER
!** NUMBER OF STOP-CRITERION
        INTEGER             ST_CR
!** CONTROLS OF RCI CYCLE
        INTEGER             SUCCESSFUL
!** MAXIMUM NUMBER OF ITERATIONS
        INTEGER             ITER1
!** MAXIMUM NUMBER OF ITERATIONS OF CALCULATION OF TRIAL-STEP
        INTEGER             ITER2
!** INITIAL STEP BOUND
        DOUBLE PRECISION    RS
!** INITIAL AND FINAL RESIDUALS
        DOUBLE PRECISION    R1, R2
!** TR SOLVER HANDLE
        INTEGER*8           HANDLE
!** CYCLE'S COUNTERS
        INTEGER             I, J
!** RESULTS OF INPUT PARAMETER CHECKING
        INTEGER             INFO(6)
!** SET PRECISIONS FOR STOP-CRITERIA
        DO I = 1, 6
            EPS (I) = 1.D-5
        END DO
!** SET MAXIMUM NUMBER OF ITERATIONS
        ITER1 = 1000
!** SET MAXIMUM NUMBER OF ITERATIONS OF CALCULATION OF TRIAL-STEP
        ITER2 = 100
!** SET INITIAL STEP BOUND
        RS = 100.D0
!** PRECISIONS FOR JACOBI CALCULATION
        JAC_EPS = 1.D-8
!** SET THE INITIAL GUESS
        x=0
 
      
!** SET INITIAL VALUES
        DO I = 1, M
            FVEC (I) = 0.D0
            DO J = 1, N
                FJAC (I, J) = 0.D0
            END DO
        END DO
!** INITIALIZE SOLVER (ALLOCATE MEMORY, SET INITIAL VALUES)
!**   HANDLE    IN/OUT: TR SOLVER HANDLE
!**   N         IN:     NUMBER OF FUNCTION VARIABLES
!**   M         IN:     DIMENSION OF FUNCTION VALUE
!**   X         IN:     SOLUTION VECTOR. CONTAINS VALUES X FOR F(X)
!**   EPS       IN:     PRECISIONS FOR STOP-CRITERIA
!**   ITER1     IN:     MAXIMUM NUMBER OF ITERATIONS
!**   ITER2     IN:     MAXIMUM NUMBER OF ITERATIONS OF CALCULATION OF TRIAL-STEP
!**   RS        IN:     INITIAL STEP BOUND
        IF (DTRNLSP_INIT (HANDLE, N, M, X, EPS, ITER1, ITER2, RS)
     &      .NE. TR_SUCCESS) THEN
!** IF FUNCTION DOES NOT COMPLETE SUCCESSFULLY THEN PRINT ERROR MESSAGE
            PRINT *, '| ERROR IN DTRNLSP_INIT'
!** RELEASE INTERNAL Intel(R) MKL MEMORY THAT MIGHT BE USED FOR COMPUTATIONS.
!** NOTE: IT IS IMPORTANT TO CALL THE ROUTINE BELOW TO AVOID MEMORY LEAKS
!** UNLESS YOU DISABLE Intel(R) MKL MEMORY MANAGER
            CALL MKL_FREE_BUFFERS
!** AND STOP
            STOP 1
        END IF
!** CHECKS THE CORRECTNESS OF HANDLE AND ARRAYS CONTAINING JACOBIAN MATRIX, 
!** OBJECTIVE FUNCTION, LOWER AND UPPER BOUNDS, AND STOPPING CRITERIA.
        IF (DTRNLSP_CHECK (HANDLE, N, M, FJAC, FVEC, EPS, 
     &     INFO) .NE. TR_SUCCESS) THEN
!** IF FUNCTION DOES NOT COMPLETE SUCCESSFULLY THEN PRINT ERROR MESSAGE
            PRINT *, '| ERROR IN DTRNLSPBC_INIT'
!** RELEASE INTERNAL Intel(R) MKL MEMORY THAT MIGHT BE USED FOR COMPUTATIONS.
!** NOTE: IT IS IMPORTANT TO CALL THE ROUTINE BELOW TO AVOID MEMORY LEAKS
!** UNLESS YOU DISABLE Intel(R) MKL MEMORY MANAGER
            CALL MKL_FREE_BUFFERS
!** AND STOP
            STOP 1
        ELSE
!** THE HANDLE IS NOT VALID.
            IF( INFO(1) .NE. 0 .OR. 
!** THE FJAC ARRAY IS NOT VALID.
     &              INFO(2) .NE. 0 .OR. 
!** THE FVEC ARRAY IS NOT VALID.
     &              INFO(3) .NE. 0 .OR. 
!** THE EPS ARRAY IS NOT VALID.
     &              INFO(4) .NE. 0 ) THEN
                PRINT *, '| INPUT PARAMETERS ARE NOT VALID'
!** RELEASE INTERNAL Intel(R) MKL MEMORY THAT MIGHT BE USED FOR COMPUTATIONS.
!** NOTE: IT IS IMPORTANT TO CALL THE ROUTINE BELOW TO AVOID MEMORY LEAKS
!** UNLESS YOU DISABLE Intel(R) MKL MEMORY MANAGER
                CALL MKL_FREE_BUFFERS
!** AND STOP
                STOP 1
            END IF
        END IF
!** SET INITIAL RCI CYCLE VARIABLES
        RCI_REQUEST = 0
        SUCCESSFUL = 0
!** RCI CYCLE
        DO WHILE (SUCCESSFUL == 0)
!** CALL TR SOLVER
!**   HANDLE        IN/OUT: TR SOLVER HANDLE
!**   FVEC          IN:     VECTOR
!**   FJAC          IN:     JACOBI MATRIX
!**   RCI_REQUEST   IN/OUT: RETURN NUMBER WHICH DENOTE NEXT STEP FOR PERFORMING
            IF (DTRNLSP_SOLVE (HANDLE, FVEC, FJAC, RCI_REQUEST)
     &              .NE. TR_SUCCESS) THEN
!** IF FUNCTION DOES NOT COMPLETE SUCCESSFULLY THEN PRINT ERROR MESSAGE
                PRINT *, '| ERROR IN DTRNLSP_SOLVE'
!** RELEASE INTERNAL Intel(R) MKL MEMORY THAT MIGHT BE USED FOR COMPUTATIONS.
!** NOTE: IT IS IMPORTANT TO CALL THE ROUTINE BELOW TO AVOID MEMORY LEAKS
!** UNLESS YOU DISABLE Intel(R) MKL MEMORY MANAGER
                CALL MKL_FREE_BUFFERS
!** AND STOP
                STOP 1
            END IF
!** ACCORDING WITH RCI_REQUEST VALUE WE DO NEXT STEP
            SELECT CASE (RCI_REQUEST)
            CASE (-1, -2, -3, -4, -5, -6)
!**   STOP RCI CYCLE
                SUCCESSFUL = 1
            CASE (1)
!**   RECALCULATE FUNCTION VALUE
!**     M               IN:     DIMENSION OF FUNCTION VALUE
!**     N               IN:     NUMBER OF FUNCTION VARIABLES
!**     X               IN:     SOLUTION VECTOR
!**     FVEC            OUT:    FUNCTION VALUE F(X)
	   k=190000/3/(1-2*0.3);
         g=190000/(1+0.3)/2;
         pe= -167.6485
         qe= 500.0385;
                CALL EXTENDED_POWELL (M, N, X, FVEC)
            CASE (2)
!**   COMPUTE JACOBI MATRIX
!**     EXTENDED_POWELL IN:     EXTERNAL OBJECTIVE FUNCTION
!**     N               IN:     NUMBER OF FUNCTION VARIABLES
!**     M               IN:     DIMENSION OF FUNCTION VALUE
!**     FJAC            OUT:    JACOBI MATRIX
!**     X               IN:     SOLUTION VECTOR
!**     JAC_EPS         IN:     JACOBI CALCULATION PRECISION
                IF (DJACOBI (EXTENDED_POWELL, N, M, FJAC, X, JAC_EPS)
     &                  .NE. TR_SUCCESS) THEN
!** IF FUNCTION DOES NOT COMPLETE SUCCESSFULLY THEN PRINT ERROR MESSAGE
                    PRINT *, '| ERROR IN DJACOBI'
!** RELEASE INTERNAL Intel(R) MKL MEMORY THAT MIGHT BE USED FOR COMPUTATIONS.
!** NOTE: IT IS IMPORTANT TO CALL THE ROUTINE BELOW TO AVOID MEMORY LEAKS
!** UNLESS YOU DISABLE Intel(R) MKL MEMORY MANAGER
                    CALL MKL_FREE_BUFFERS
!** AND STOP
                    STOP 1
                END IF
            ENDSELECT
        END DO
!** GET SOLUTION STATUSES
!**   HANDLE            IN: TR SOLVER HANDLE
!**   ITER              OUT: NUMBER OF ITERATIONS
!**   ST_CR             OUT: NUMBER OF STOP CRITERION
!**   R1                OUT: INITIAL RESIDUALS
!**   R2                OUT: FINAL RESIDUALS
        IF (DTRNLSP_GET (HANDLE, ITER, ST_CR, R1, R2)
     &          .NE. TR_SUCCESS) THEN
!** IF FUNCTION DOES NOT COMPLETE SUCCESSFULLY THEN PRINT ERROR MESSAGE
            PRINT *, '| ERROR IN DTRNLSP_GET'
!** RELEASE INTERNAL Intel(R) MKL MEMORY THAT MIGHT BE USED FOR COMPUTATIONS.
!** NOTE: IT IS IMPORTANT TO CALL THE ROUTINE BELOW TO AVOID MEMORY LEAKS
!** UNLESS YOU DISABLE Intel(R) MKL MEMORY MANAGER
            CALL MKL_FREE_BUFFERS
!** AND STOP
            STOP 1
        END IF
!** FREE HANDLE MEMORY
        IF (DTRNLSP_DELETE (HANDLE) .NE. TR_SUCCESS) THEN
!** IF FUNCTION DOES NOT COMPLETE SUCCESSFULLY THEN PRINT ERROR MESSAGE
            PRINT *, '| ERROR IN DTRNLSP_DELETE'
!** RELEASE INTERNAL Intel(R) MKL MEMORY THAT MIGHT BE USED FOR COMPUTATIONS.
!** NOTE: IT IS IMPORTANT TO CALL THE ROUTINE BELOW TO AVOID MEMORY LEAKS
!** UNLESS YOU DISABLE Intel(R) MKL MEMORY MANAGER
            CALL MKL_FREE_BUFFERS
!** AND STOP
            STOP 1
        END IF

!** RELEASE INTERNAL Intel(R) MKL MEMORY THAT MIGHT BE USED FOR COMPUTATIONS.
!** NOTE: IT IS IMPORTANT TO CALL THE ROUTINE BELOW TO AVOID MEMORY LEAKS
!** UNLESS YOU DISABLE Intel(R) MKL MEMORY MANAGER
        CALL MKL_FREE_BUFFERS
!** IF FINAL RESIDUAL LESS THEN REQUIRED PRECISION THEN PRINT PASS
        IF (R2 .LT. 1.D-6) THEN
            PRINT *, '|         DTRNLSP POWELL............PASS'
            print *,"x=",x
            print *,"r2=",r2
            !pause
            !STOP 0
!** ELSE PRINT FAILED
        ELSE
            PRINT *, '|         DTRNLSP POWELL............FAILED'
            print *, 'r2=',r2
            print *, 'x=',x
            !STOP 1
		
        END IF
      END PROGRAM EXAMPLE_DTRNLSP_POWELL

!** ROUTINE FOR EXTENDED POWELL FUNCTION CALCULATION
!**   M     IN:     DIMENSION OF FUNCTION VALUE
!**   N     IN:     NUMBER OF FUNCTION VARIABLES
!**   X     IN:     VECTOR FOR FUNCTION CALCULATING
!**   F     OUT:    FUNCTION VALUE F(X)
      SUBROUTINE EXTENDED_POWELL (M, N, X, F)
        IMPLICIT NONE
        INTEGER M, N
        DOUBLE PRECISION X(*),F(*),K,G,PE,QE,SM
     & ,P,Q,PEEQ,FREAL,DS,DC,A
        INTEGER I
       COMMON PE,QE,K,G  
		  P=PE+K*X(1);
            Q=QE-3*G*X(2);
		  PEEQ=X(3);
		  SM=500+1000*(PEEQ);
            FREAL=X(4);
            DS=SINH(-1.5*P/SM);
            DC=COSH(-1.5*P/SM);
            A=0.04*EXP(-0.5*((PEEQ-0.3)/0.1)**2)/((2*3.14159)**0.5)/0.1;
          
          F(1) = X(1)* 2*Q/SM/SM+X(2)*(-3*1.5*FREAL*DS/SM);
          F(2) = (Q/SM)**2+2*FREAL*1.5*DC-(1+2.25*FREAL**2);
	    F(3) = X(3)-(-P*X(1)+Q*X(2))/(1-FREAL)/SM
	    F(4) = X(4)-(1-FREAL)*X(1)-A*PEEQ

      END SUBROUTINE EXTENDED_POWELL

Zhen_Z_Intel · ‎12-20-2016

Dear customer,

The calculation relay on EPS, you set the EPS as 1E-5, it is too large for your sample. You'd better set smaller than the smallest value of x you would like to get, for instance 1.D-17. The result of MKL should be same as Matlab. That means all result of x would not be smaller than 1.D-17, thus the function value will also be different. You could print the result of function "FVEC" to check with it. Thanks.

Best regards,
Fiona

View solution in original post

Zhen_Z_Intel · ‎12-19-2016

Dear customer,

In MKL the final result is 1.540586674360078E-004, failed because you set the limitation for PASS is 1.D-6. I tried with your matlab code, seems there's an error. Could you please provide the output for matlab program and the init value & result value of x(1)~x(4) and the final f(x) value. Thanks.

Best regards,
Fiona

mecej4 · ‎12-20-2016

A good rule to follow when solving nonlinear equations is not to use all zeros as the guess for the solution, and a second rule is not to give up when an arbitrary initial guess fails to yield a solution.

I used x₀ = [1d-7, 1d-7, 1d-7, 1d-7], and the solution found by TRNLSP was x = [3.3209d-12, 1.7481d-07, 1.7482d-07, 3.1324d-10] .

Perhaps you need to re-scale your variables x and F so that you can work with reasonable values, or use only a relative error in the termination criterion.

peng_z_ · ‎12-20-2016

Dear Fiona,  
Thank you for your prompt reply.
I just changed "x0=0" to "x0=[0,0,0,0]" in the Matlab code;
Here is a table showing the input (x0)j (j=1,2,3,4),
the solution x, and the residual F. The correct code is posted below the table. 
Seems like Matlab can get the maxinum residual to a order of 1e-14.

index j	1	2	3	4
initial guess x0	0	0	0	0
Final value of x	6.392630365754150e-17	1.74817817657479e-07	1.74817817600321e-07	3.09908247048591e-10
Final value of F	-4.26428382337290e-28	-1.78745906964650e-14	-4.10047855672419e-17	-4.78885875879918e-21

!*********save the function as EQs.m file*****

function F = EQs(x)

k=190000/3/(1-2*0.3);
g=190000/(1+0.3)/2;
pe=-167.6485 ;
qe=500.0385  ;
p=pe+k*x(1);
q=qe-3*g*x(2);
peeq=x(3);
sm=500+1000*(peeq);
freal=x(4);
ds=sinh(-1.5*p/sm);
dc=cosh(-1.5*p/sm);
a=0.04*exp(-0.5*((peeq-0.3)/0.1)^2)/((2*3.14159)^0.5)/0.1;

F(1) = x(1)* 2*q/sm/sm+x(2)*(-3*1.5*freal*ds/sm);
F(2) = (q/sm)^2+2*freal*1.5*dc-(1+2.25*freal^2);
F(3)=x(3)-(-p*x(1)+q*x(2))/(1-freal)/sm
F(4)=x(4)-(1-freal)*x(1)-a*peeq

!**************end of EQs.m file**************

>>options = optimoptions('fsolve','Display','iter','PlotFcn',@optimplotfirstorderopt);
>>x0=[0,0,0,0];
>>[x,fval]  = fsolve(@EQs,x0,options)

peng_z_ · ‎12-20-2016

Hi mecej4,

Thanks for your reply, your result is very close to what I have got via Matlab using initial guess x0=[0,0,0,0], pls see the result in my above comment.

However, I tried x0=[1d-7,1d-7,1d-7,1d-7,], failed still. Have you changed other line of Fortran code?
What do you mean by "rescale" x and F? Could you pls elaborate on how to do so in TRNLSP ?

Kind regards,

P

Zhen_Z_Intel · ‎12-20-2016

Dear customer,

The calculation relay on EPS, you set the EPS as 1E-5, it is too large for your sample. You'd better set smaller than the smallest value of x you would like to get, for instance 1.D-17. The result of MKL should be same as Matlab. That means all result of x would not be smaller than 1.D-17, thus the function value will also be different. You could print the result of function "FVEC" to check with it. Thanks.

Best regards,
Fiona