How to make a matrix where its elements are functions, operate with them and the result still be a function?

Question

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?

Answer 1

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:

• have the actual argument a procedure pointer;
• have the dummy argument a procedure (non-pointer);
• allow argument association between a pointer and non-pointer.

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).