Error : The CALL statement is invoking a function subprogram as a subroutine.

[vahid_a_1]

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

[mecej4 · Accepted solution]

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.

[vahid_a_1]

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

[mecej4]

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^-1f

   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.

[vahid_a_1]

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 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^-1f

   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.

[mecej4]

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.