CODEWITHSUNDEEP
fortran
Simon Ranjith
Simon Ranjith

Reputation: 37

How to assign pointers to a complex 3d array in Fortran

I would like to know how to assign two pointers, one to the real part of a complex 3d array and another to the imaginary part of the same array in Fortran.

Let's say I have defined a 3d array as such:

complex*16, dimension(:,:,:), allocatable, target :: vftmp

and I would like to assign a pointer to the real part of vftmp(2,1,1) and a pointer to the imaginary part of vftmp(2,1,1). Could someone help me with a snippet please? Thanks.

Upvotes: 1

Views: 872

Answers (1)

roygvib
roygvib

Reputation: 7395

I hope something like the following is possible

real, pointer :: re
complex, target :: z
re => z % re
! or
real, pointer :: re(:,:,:)
complex, target :: z(2,3,4)
re => z(:,:,:) % re

but it seems not (or possible with very new compilers...?) So a workaround approach below:

1) If the goal is to get (scalar) pointers for the Re and Im parts of a single element of a complex array, I guess we can use c_f_pointer such that

module testmod
contains
    subroutine getreim_ptr( z, re, im )
        use iso_c_binding
        implicit none
        complex, target, intent(in) :: z
        real, pointer :: re, im, buf(:)

        call c_f_pointer( c_loc( z ), buf, [ 2 ] )

        re => buf( 1 )
        im => buf( 2 )
    end subroutine
end module

program main
    use testmod
    implicit none
    complex :: z( 2, 3 )
    real, pointer :: re, im

    !! Test array.
    z = 0.0
    z( 1, 1 ) = ( 1.0, -1.0 )

    !! Get pointers for the Re/Im parts of z(1,1).
    call getreim_ptr( z( 1, 1 ), re, im )

    print *, "z(1,:) = ", z(1,:)
    print *, "z(2,:) = ", z(2,:)
    print *, "re = ", re
    print *, "im = ", im
end

Result (gfortran-8.2):

 z(1,:) =   (1.00000000,-1.00000000)  (0.00000000,0.00000000)  (0.00000000,0.00000000)
 z(2,:) =   (0.00000000,0.00000000)   (0.00000000,0.00000000)  (0.00000000,0.00000000)
 re =    1.00000000    
 im =   -1.00000000

2) If the goal is to get array pointers for the entire complex array, I guess we can use rank-remapping pointer assignments (to point to non-contiguous memory with constant gaps). For example, in the 2D case (for simplicity),

re( 1:n1, 1:n2 ) => buf( 1::2 )
im( 1:n1, 1:n2 ) => buf( 2::2 )

where re and im are 2D array pointers and buf is a real 1D array pointer that points to an allocatable 2D complex array (via c_f_pointer). A minimum example may look like this:

module testmod
contains
    subroutine getreim_ptr2d( zarr, re, im )
        use iso_c_binding
        implicit none
        complex, allocatable, target, intent(in) :: zarr(:,:)
        real, pointer :: re(:,:), im(:,:), buf(:)
        integer :: n1, n2

        n1 = size( zarr, 1 )
        n2 = size( zarr, 2 )
        call c_f_pointer( c_loc( zarr ), buf, [ size(zarr) * 2 ] )

        re( 1:n1, 1:n2 ) => buf( 1::2 )
        im( 1:n1, 1:n2 ) => buf( 2::2 )
    end subroutine
end module

program main
    use testmod
    implicit none
    complex, allocatable :: zarr(:,:)
    real, pointer :: re(:,:), im(:,:)
    integer i

    !! Prepare a test array (zarr).
    allocate( zarr( 2, 3 ) )
    zarr(1,:) = [( complex( 100 + i, -100 -i ), i=1,3 )]
    zarr(2,:) = [( complex( 200 + i, -200 -i ), i=1,3 )]
    print *, "shape( zarr ) = ", shape( zarr )
    print *, "zarr(1,:) = ", zarr(1,:)
    print *, "zarr(2,:) = ", zarr(2,:)

    call getreim_ptr2d( zarr, re, im )

    print *
    print *, "shape( re ) = ", shape( re )
    print *, "re(1,:) = ", re(1,:)
    print *, "re(2,:) = ", re(2,:)
    print *
    print *, "shape( im ) = ", shape( im )
    print *, "im(1,:) = ", im(1,:)
    print *, "im(2,:) = ", im(2,:)
end program

Result (gfortran 8.2):

 shape( zarr ) =            2           3
 zarr(1,:) =  (101.000000,-101.000000)  (102.000000,-102.000000)  (103.000000,-103.000000)
 zarr(2,:) =  (201.000000,-201.000000)  (202.000000,-202.000000)  (203.000000,-203.000000)

 shape( re ) =            2           3
 re(1,:) =    101.000000    102.000000    103.000000    
 re(2,:) =    201.000000    202.000000    203.000000    

 shape( im ) =            2           3
 im(1,:) =   -101.000000   -102.000000   -103.000000    
 im(2,:) =   -201.000000   -202.000000   -203.000000

Below are some materials we can find on the net:

The New Features of Fortran 2003 (N1597): 3.7 "Pointer assignment"

"...Remapping of the elements of a rank-one array is permitted:

p(1:m,1:2*m) => a(1:2*m*m)

The mapping is in array-element order and the target array must be large enough. The bounds may be any scalar integer expressions. The limitation to rank-one arrays is because pointer arrays need not occupy contiguous storage:

a => b(1:10:2)

but all the gaps have the same length in the rank-one case."

Fortran 2003 extensions: 5.4.3 Rank-remapping Pointer Assignment (this page)

"...This feature allows a multi-dimensional pointer to point to a single-dimensional object. For example:

REAL,POINTER :: diagonal(:),matrix(:,:),base(:)
...
ALLOCATE(base(n*n))
matrix(1:n,1:n) => base
diagonal => base(::n+1)
!
! DIAGONAL now points to the diagonal elements of MATRIX.
!

Note that when rank-remapping, the values for both the lower and upper bounds must be explicitly specified for all dimensions, there are no defaults."

Upvotes: 1

Related Questions