Solved: Dear customer,

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

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