I'm using Fortran I'm trying to create matrices where their elements are functions. Also I'd like to operate with them and the result still be a function. So here is what I try
module Greeninverse
use, intrinsic :: iso_fortran_env, only: dp => real64
implicit none
real(dp), public, parameter :: wl = 1d0
real(dp), public, parameter :: wr = 1d0
integer, public, parameter :: matrix_size = 5
type ptr_wrapper
procedure(f), nopass, pointer :: func
end type ptr_wrapper
abstract interface
function f(x1,x2)
import
real(dp), intent(in) :: x1
real(dp), intent(in) :: x2
complex (dp), dimension(matrix_size,matrix_size):: f
end function f
end interface
contains
function Sigma(x1) result(S)
real(dp),intent(in) :: x1
complex(dp), dimension(matrix_size,matrix_size) :: S
real(dp):: aux_wr1,aux_wl1
complex(dp) :: S11, Snn
integer :: i,j
aux_wr1 = 1-x1**2/(2d0*wr)
aux_wl1 = 1-x1**2/(2d0*wl)
S11 = dcmplx(.5*(x1**2-2d0*wl), 2.0*wL*dsqrt(1-aux_wL1**2))
Snn = dcmplx(.5*(x1**2-2d0*wr), 2.0*wr*dsqrt(1-aux_wr1**2))
do i = 1, matrix_size
do j=i,matrix_size
S(i,j) = 0d0
S(j,i) = 0d0
end do
end do
S(1,1) = S11
S(matrix_size,matrix_size) = Snn
end function Sigma
function Omega(x1) result(Om)
real(dp),intent(in) :: x1
real(dp),dimension(matrix_size, matrix_size) :: Om
integer :: i,j
do i=1,matrix_size
do j= i, matrix_size
Om(i,j) = 0d0
Om(j,i) = 0d0
end do
end do
do i = 1,matrix_size
Om(i,i) = x1**2
end do
end function Omega
! Now I'd like to add them and take the inverse of the sum and still be a function
function Inversa(x1,x2) result (G0inv)
real(dp), intent(in) :: x1
real(dp), intent(in) :: x2
complex(dp), dimension(matrix_size,matrix_size) :: G0inv
complex(dp),dimension(matrix_size,matrix_size) :: Gaux
! Down here all these variables are needed by ZGETRF and ZGETRI
DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: WORK
Integer:: LWORK = matrix_size*matrix_size
Integer, Allocatable, dimension(:) :: IPIV
Integer :: INFO, LDA = matrix_size, M = matrix_size, N = matrix_size
Integer DeAllocateStatus
external liblapack
allocate(work(Lwork))
allocate(IPIV(N))
Gaux = Omega(x1)+Sigma(x2)
CALL ZGETRF (M, N, Gaux, LDA, IPIV, INFO)
! This calculates LU descomposition of a matrix and overwrites it
CALL ZGETRI(N, Gaux, N, IPIV, WORK, LWORK, INFO)
! This calculates the inverse of a matrix knowing its LU descomposition and overwrites it
G0inv = Gaux
end function Inversa
! Now I'd like to derive it
function Derivate(x1,x2,G) result(d)
! This function is supposed to derivate a matrix which its elements are functions but of two variables; x1 and x2. And it only derives respect the first variable
implicit none
real(dp), intent(in) :: x1
real(dp), intent(in) :: x2
procedure(f),pointer:: G
complex(dp),dimension(matrix_size,matrix_size) :: d
real(dp) :: h = 1.0E-6
d = (1.0*G(x1-2*h,x2) - 8.0*G(x1-h,x2) + 8.0*G(x1+h,x2) - 1.0*G(x1+2*h,x2))/(12.0*h)
end function Derivate
end module Greeninverse
program Greentest3
use, intrinsic :: iso_fortran_env, only: dp => real64
use Greeninverse
implicit none
real(dp) :: W(matrix_size,matrix_size)
complex(dp) :: S(matrix_size,matrix_size)
complex(dp) :: G(matrix_size,matrix_size)
complex(dp) :: DD(matrix_size,matrix_size)
W(:,:) = Omega(1d0)
S(:,:) = Sigma(2d0)
G(:,:) = Inversa(1d0,2d0)
DD(:,:) = Derivate(1d0,2d0,Inversa)
print*, W
print*, S
print*, G
print*, DD
end program Greentest3
The problem is in the function Derivate that I don't know how to say that the argument G is a matrix function and because of that I get an error message
DD(:,:) = Derivate(1d0,2d0,Inversa)
1
Error: Expected a procedure pointer for argument ‘g’ at (1)
That's why I use the abstract interface that it's supposed to say that is a function but it doesn't work as I expected
I tried also to make a pointer in the module section, that is
type(ptr_wrapper) :: DD(matrix_size,matrix_size)
but I get an error message
Error: Unexpected data declaration statement in CONTAINS section at (1)
I'd like to make all the matrices in the module section and in the program just evaluate them in the values of interest.
What am I doing wrong?
Looking at the function Derivate the dummy argument G is declared like
procedure(f), pointer:: G
This is a procedure pointer. The error message confirms that.
The actual argument to be passed to Derivate is, in this case, expected also to be a procedure pointer. Let's look at what the argument is:
DD(:,:) = Derivate(...,Inversa)
Inversa is a procedure (function), defined in the module. It, crucially, isn't a procedure pointer. So, indeed, the compiler complains.
Well, how do we go about fixing this? There are three obvious approaches:
For the first, the main program could have
procedure(f), pointer :: Inversa_ptr ! We've a procedure pointer...
Inversa_ptr => Inversa ! ... which we point at our procedure...
DD(:,:) = Derivate(...,Inversa_ptr) ! ... and is then the argument
For the Derivate as it is implemented, it doesn't use the pointer nature of the argument G: just the target is referenced. This means that the other two options become available.
We can make the dummy argument not a pointer, having
function Derivate(...,G)
procedure(f) :: G
end function
used like
DD(:,:) = Derivate(...,Inversa)
The third of our choices comes from defining the dummy argument as
function Derivate(...,G)
procedure(f), pointer, intent(in) :: G
end function
where, again, the reference is as in the second case.
When the dummy argument procedure pointer has the intent(in) attribute, it is allowed to be associated with a non-pointer procedure which is a valid target in pointer assignment. In this case G becomes pointer associated with that actual argument procedure (and because of the intent, that status can't be changed in the function).