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

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

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

1 Solution Black Belt
356 Views

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.

5 Replies Black Belt
357 Views

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. Beginner
356 Views

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
``` Black Belt
356 Views

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. Beginner
356 Views

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. Black Belt
356 Views

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. 