Community
cancel
Showing results for
Did you mean: Beginner
230 Views

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

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 :

```!===============================================================================
!
! 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```

1 Solution Employee
230 Views

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

5 Replies Employee
230 Views

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 Black Belt
230 Views

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 x0 = [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. Beginner
230 Views
```Dear Fiona,
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.39263e-17 1.74818e-07 1.74818e-07 3.09908e-10 Final value of F -4.26428e-28 -1.78746e-14 -4.10048e-17 -4.78886e-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)``` Beginner
230 Views

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 Employee
231 Views

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 