!integration matrix test ! this is a routine for integration functions with an example at the end that integrates a matrix whose elements are functions !It's written in fortran 2003 module intmat use, intrinsic :: iso_fortran_env, only: dp => real64 implicit none ! Wrapper for a procedure pointer so we can make an array of them type ptr_wrapper procedure(f), nopass, pointer :: func end type ptr_wrapper ! Use an interface for the function rather than "extern" functions ! Change the interface to whatever you want. You could take multiple reals, integers, etc abstract interface pure function f(x1) import real(dp), intent(in) :: x1 real(dp) :: f end function f end interface ! The intent(in attribute of argument x1 means that x1 cannot be changed inside the function. Note that the return type of func needs to be declared. If this is omitted, somme compilers will not compile !When declaring variables inside functions and subroutines that need to be passed in or out, intent may be added to the declaration. !intent(in) means that the variable value can enter, but not be changed !intent(out) means the variable is set inside the procedure and sent back to the main program with any initial values ignored. !intent(inout) means that the variable comes in with a value and leaves with a value (default). contains subroutine DGS0(x1,x2,fx,S) ! Integra por gausseanos. LLama a la función DGS8. ! x1 y x2 son los límites, fx es la función y S el resultado real(dp), intent(in) :: x1 real(dp), intent(in) :: x2 procedure(f), pointer :: fx real(dp):: u, v, c, sl, sg, sp, sa, sf real(dp), intent(out) :: S real(dp) :: eps = 1.0E-6 S = 0d0 U = x1 1 V = x2 IF ( U .LT. x2 ) THEN SF = DGS8 (u,v,fx) 2 C = (U+V)/2d0 SL = DGS8 (u,c,fx) SG = DGS8 (c,v,fx) SP = SL+SG SA = DABS(SP-SF)/EPS IF ( SA.GE.1.D0 ) THEN V = C SF = SL GOTO 2 END IF U = V S = S+SP GOTO 1 END IF end subroutine DGS0 function DGS8(x1,x2,fx) result (d) ! rutina que integra. Es llamada por DGS0. ! x1 y x2 son los límites, fx la función, d el resultado implicit none real(dp), intent(in) :: x1 real(dp), intent(in) :: x2 procedure(f), pointer :: fx ! real (dp), intent(out) :: d real(dp) :: d real(dp) :: c,h,x,s c = (x1+x2)/2 h = (x2-x1)/2 X = .960289856497536D0*H S = .101228536290376D0*(fx(C+X)+fx(C-X)) X = .796666477413627D0*H S = S + .222381034453374D0*(fx(C+X)+fx(C-X)) X = .525532409916329D0*H S = S + .313706645877887D0*(fx(C+X)+fx(C-X)) X = .183434642495650D0*H S = S + .362683783378362D0*(fx(C+X)+fx(C-X)) d = S * H end function DGS8 pure function test_func(x1) real(dp), intent(in) :: x1 real(dp) test_func test_func = 2d0 * x1 end function test_func pure function test_func2(x1) real(dp), intent(in) :: x1 real(dp) test_func2 test_func2 = x1**2 end function test_func2 end module intmat program integral_matrix use, intrinsic :: iso_fortran_env, only: dp => real64 use intmat implicit none integer i, j real(dp) :: dx, S ! real(dp) :: pi2 = acos(0d0) type(ptr_wrapper) :: w(2,2) !Point to the function that you want each index to point to w(1,1)%func => test_func w(1,2)%func => test_func2 ! Pointing to the generic function wasn't working so I had to point to the specific double value trig functions w(2,1)%func => dcos w(2,2)%func => dsin do i = 1, 2 do j = 1, 2 call DGS0(0d0,1.5d0, w(i,j)%func,S) dx = S print *, dx end do end do print *, 'Done :)' end program integral_matrix