1

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.

1
  • Assuming you mean pointer assignment, you can't. The real part of vftmp is of course real(vftmp,16). You cannot do, for example, re => real(vftmp,16). BTW, 16 is a nonstandard kind type parameter. Commented Feb 22, 2019 at 0:19

1 Answer 1

1

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."

Sign up to request clarification or add additional context in comments.

3 Comments

I don't understand the notations buf(1::2) and buf(2::2) in your code. The value 2 after the colon, is it the stride? Could you explain please? Thanks.
It is the same as writing buf( 1 : 2*n1*n2 : 2 ) and buf( 2 : 2*n1*n2 : 2 ). We can omit any of a, b, inc in the notation a:b:inc, then they are assumed automatically to be the lower bound (for a), upper bound (for b), and 1 (for inc). (For example, if A is declared as A( 6 ), then A( 1::2 ) (= A( 1:6:2 )) means A(1), A(3), A(5), and A( 2::2 ) (= A( 2:6:2 )) means A(2), A(4), A(6).)
Just a late addition... The above method creates an array pointer for non-contiguous memory, which I think should not be passed to a subroutine that expects an array of contiguous memory (e.g. legacy library routines with explicit-shape or assumed-size dummy arguments). So please be careful...

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.