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

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)

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

我还尝试在模块部分中创建一个指针，即

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 是一个过程(函数)，在模块中定义。至关重要的是，它不是一个过程指针。所以，编译器确实会提示。

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

• 让实际参数成为过程指针；
• 让虚拟参数成为一个过程(非指针)；
• 允许指针和非指针之间的参数关联。

首先，主程序可以有

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 成为与该实际参数过程关联的指针(并且由于意图，该状态无法在函数中更改)。