Как создать матрицу, где ее элементы являются функциями, работать с ними, и результат все еще является функцией?

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 

Проблема заключается в функции деривата, что я не знаю, как сказать, что аргумент 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 процедура (функция), определенная в модуле. Это, собственно, не процедура указатель. Так, действительно, компилятор жалуется.

Ну, как мы можем это исправить? Существует три очевидных подхода:

• имеет фактический аргумент указателя процедуры;
• содержит аргумент фиктивной процедуры (не указатель);
• Разрешить ассоциацию аргументов между указателем и не указателем.

Для первых, основная программа может иметь

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 становится указателем, связанным с этой фактической процедурой аргумента (и из-за намерения этот статус не может быть изменен в функции).