如何制作一个矩阵，其中的元素是函数，对它们进行运算，结果仍然是一个函数？

Question

我正在使用 Fortran 我正在尝试创建其元素为函数的矩阵。我也想和他们一起操作，结果仍然是一个函数。所以这就是我的尝试

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

问题出在函数 Derivate 中，我不知道如何说参数 G 是矩阵函数，因此我收到一条错误消息

  DD(:,:) = Derivate(1d0,2d0,Inversa)
                         1
Error: Expected a procedure pointer for argument ‘g’ at (1)

这就是为什么我使用抽象接口，它应该说它是一个函数，但它并没有像我预期的那样工作

我也试过在模块部分做一个指针，就是

type(ptr_wrapper) :: DD(matrix_size,matrix_size)

但我收到一条错误消息

Error: Unexpected data declaration statement in CONTAINS section at (1)

我想在模块部分和程序中制作所有矩阵，只用感兴趣的值计算它们。

我做错了什么？

Answer 1

查看函数 Derivate 伪参数 G 声明为

procedure(f), pointer:: G

这是一个过程指针。错误消息证实了这一点。

在这种情况下，要传递给 Derivate 的实际参数也应该是一个过程指针。让我们看看参数是什么：

DD(:,:) = Derivate(...,Inversa)

Inversa是一个过程（函数），定义在模块中。至关重要的是，它不是过程 pointer。所以，确实，编译器会抱怨。

那么，我们该如何解决这个问题？共有三种明显的方法：

让实参成为过程指针；
将伪参数设置为过程（非指针）；
允许指针和非指针之间的参数关联。

首先，主程序可以有

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

对于 Derivate 的实现，它不使用参数 G 的指针性质：仅引用目标。这意味着其他两个选项可用。

我们可以使伪参数不是指针，

function Derivate(...,G)
   procedure(f) :: G
end function

像

一样使用

DD(:,:) = Derivate(...,Inversa)

我们的第三个选择来自将伪参数定义为

function Derivate(...,G)
   procedure(f), pointer, intent(in) :: G
end function

其中，引用与第二种情况相同。

当伪参数过程指针具有intent(in)属性时，允许与作为指针赋值的有效目标的非指针过程相关联。在这种情况下 G 成为与实际参数过程关联的指针（并且由于意图，该状态不能在函数中更改）。

如何制作一个矩阵，其中的元素是函数，对它们进行运算，结果仍然是一个函数？

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

module

fortran

function-pointers

derivative

matrix-inverse