- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

I am trying to learn to use MKL routines. I have written a simple code for calculating the jacobian matrix. However, it ends up with the below error:

The CALL statement is invoking a function subprogram as a subroutine. [DJACOBI]

Here is my code and I am compiling it with: *ifort** -**mkl** newton2.f90*

INCLUDE '/home/vahid/intel/composer_xe_2015.3.187/mkl/include/mkl_rci.f90' program newton2 implicit none real*8, dimension(3) :: x real*8, dimension(3,3) :: fjac real*8 :: phi1,phi2,phi3 integer :: i,m,n integer, parameter :: maxit=100 real*8, parameter :: eps=0.0001 INCLUDE '/home/vahid/intel/composer_xe_2015.3.187/mkl/include/mkl_rci.fi' external fcn x(1)=0 x(2)=0 x(3)=0 n=3;m=3 !Call fcn (m, n, x, f) call djacobi(fcn , n, m, fjac, x, eps) write(*,*) fjac end program newton2 subroutine fcn (m, n, x, f) real*8, dimension(3) :: x,f real*8 :: phi1,phi2,phi3 phi1=0 ; phi2 = 0.1; phi3= 0.9; f(1) = 2*(800) *(x(1)-0.7) - 2*(12000)*(x(2)-0.45); f(2) = 2*(12000)*(x(2)-0.45)- 2*(1200)*(x(3)-0.9) ; f(3) = x(1)*phi1 + x(2)*phi2 + x(3)*phi3 - 0.42 ; end subroutine fcn

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

Please read the MKL documentation for "Nonlinear Optimization Solvers - Jacobian Matrix Calculation Routines". The proper invocation is

iret = djacobi(fcn , n, m, fjac, x, eps)

where **iret **is an integer status code. The include file mkl_rci.fi contains the interface of this routine, and your calling it as a subroutine is in error.

Link Copied

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

Please read the MKL documentation for "Nonlinear Optimization Solvers - Jacobian Matrix Calculation Routines". The proper invocation is

iret = djacobi(fcn , n, m, fjac, x, eps)

where **iret **is an integer status code. The include file mkl_rci.fi contains the interface of this routine, and your calling it as a subroutine is in error.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

Thanks @mesej4. Your comment solved the issue.

Now, I have updated the code accordingly and I am trying to finalize the code by writing a simple newton-raphson solver for a system of nonlinear equations. I always have difficulty in the codes which use "interface".

Here is the updated code. I would like to pass f values from the fcn "subroutine" to the main program and perform a "matmul" operation. However, the elements of array "f" are zero and no information passes through the subroutine to the main program.

program newton2 implicit none real*8, dimension(3) :: x,f real*8, dimension(3,3) :: fjac integer :: i,m,n integer, parameter :: maxit=100 !! ***** djacobi variables ****** real*8, parameter :: eps=0.0001 integer :: res !! ***** dgetri variables ****** integer :: lda,lwork,info integer, dimension(3) :: ipiv real*8 :: work(3) !! ***** dgetri variables ****** !INCLUDE '/home/vahid/intel/composer_xe_2015.3.187/mkl/include/mkl_rci.fi' INCLUDE 'mkl_rci.fi' INTERFACE fcn SUBROUTINE fcn REAL*8, DIMENSION(3) :: x real*8, dimension(3) :: f END SUBROUTINE fcn END INTERFACE !! ************* Input arguments x(1)=0 ; x(2)=0 ; x(3)=0 n=3;m=3;lda = 3; !! ************* Jacobi matrix of the user's objective function using the central difference res = djacobi(fcn , n, m, fjac, x, eps) write(*,*) fjac write(*,*) fjac(2,1) !! ************* The LU factorization of a general m-by-n matrix. ipiv = max(1,n) call dgetrf( m, n, fjac, lda, ipiv, info ) write(*,*) fjac write(*,*) info !! ************* The inverse inv(A) of a general matrix A. lwork=3*1000 call dgetri( n, fjac, lda, ipiv, work, lwork, info ) write(*,*) fjac write(*,*) info !DO i=1,maxit write(*,*) x write(*,*) f write(*,*) matmul(fjac,f) x = x - matmul(fjac,f) write(*,*) x end program newton2 subroutine fcn (m, n, x, f) real*8, dimension(3) :: x real*8, dimension(3), intent(out) :: f real*8 :: phi1,phi2,phi3 phi1=0.0 ; phi2 = 0.1D0; phi3= 0.9D0; f(1) = 2*(800) *(x(1)-0.7) - 2*(12000)*(x(2)-0.45); f(2) = 2*(12000)*(x(2)-0.45)- 2*(1200)*(x(3)-0.9) ; f(3) = x(1)*phi1 + x(2)*phi2 + x(3)*phi3 - 0.42 ; return end subroutine fcn

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

The include file 'mkl_rci.fi' provides the interfaces needed, so you need not provide one for **fcn** again. Here is a sketch of what you need:

...

do iter = 1,maxit

call fcn(...) ! to compute f() at the current x() <<=== you left this out in your code

call djacobi(...) ! to compute fjac() at the current x()

call dgesv() ! to solve for the Newton Raphson correction dx = J^{-1}f

compute ||dx|| or ||f||, decide whether to terminate iterations

x = x - dx

end do

If your objective is just to obtain the solution, all that you have to do is to take the example code **ex_nlsqp.f**, which is included in the MKL installation, change m and n, put in your functions, build and run.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

Dear mecej4,

My objective is to learn to use Intel MKL whenever I need it. Thanks for the suggestions. Everything works now and it solves my main equations. However, I have another question for you. Why using intel MKL makes the code very slow?! I mean very slow...

To b more precise, I am solving a PDE and at each iteration, the code calls newton-raphson solver to solve a system of nonlinear equations.

I have two implementation of the code now. Iwhere in the 1st one, the newton-raphson solver is written by myslef and the inverse of jacobian is calculated analytically and just placed in the code manually. However, in the 2nd implementation, I have modified the newton-raphson solver by your help and used the MKL routines to do everything in Fortran. (calculating the jacobian and its inverse)

The 2nd implementation is very slow. The ratio of simulation time for 2 million iterations is : (6minutes / ~3 hours) ---> (my newton-raphson/Intel MKL)

Both of the implementations work in serial mode.

mecej4 wrote:

The include file 'mkl_rci.fi' provides the interfaces needed, so you need not provide one for

fcnagain. Here is a sketch of what you need:...

do iter = 1,maxit

call fcn(...) ! to compute f() at the current x() <<=== you left this out in your code

call djacobi(...) ! to compute fjac() at the current x()

call dgesv() ! to solve for the Newton Raphson correction dx = J

^{-1}fcompute ||dx|| or ||f||, decide whether to terminate iterations

x = x - dx

end do

If your objective is just to obtain the solution, all that you have to do is to take the example code

ex_nlsqp.f, which is included in the MKL installation, change m and n, put in your functions, build and run.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

What is the size of the system of nonlinear equations? If it is large, calculating the Jacobian by calling DJacobi() will be very wasteful and expensive. At the very least, use a variation of Newton-Raphson in which the Jacobian is updated or calculated sparingly, instead of calculating every time -- straight Newton-Raphson is not a good method for solving large systems of nonlinear equations.

To solve a linear system A x = b, it is not a good approach to form inv(A) and multiply b with that. Instead, call a routine such as DGESV, which will solve the equations twice as fast.

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page