4

I am trying to write a program where I want the allocatable array A to be of either rank 1, 2, or 3, depending on my input at run-time. I want to do this since the subsequent operations on A are similar, and I have defined in a module an interface work with module procedures that when acted on A, gives the desired result.

What I am doing currently is this:

program main
implicit none
integer :: rank,n=10
real*8, allocatable :: A1(:)
real*8, allocatable :: A2(:,:)
read (*,*) rank

if (rank.eq.1) then
    allocate (A1(n))
else if (rank.eq.2) then
    allocate (A2(n,n))
end if

! operate on the array
if (rank.eq.1) then
    call work(A1)
else if (rank.eq.2) then
    call work(A2)
end if

end program

Things would be much easier if somehow I could choose the rank of A, as then the if statements are not needed. Maybe this is not possible, but all help are appreciated.

Lun
  • 43
  • 3
  • 4
    It is not possible to allocate the rank of an array at run-time, as you seem to know. You've shown us a snippet indicating one of the possible kludges to work around this. There are others, such as judicious use of `include`; you can even wrap a rank-1 array inside a derived type and write operations to manipulate it as if it were of rank-n (choosing n at run-time). But if you tell us *why* you want this we might be able to tell you how to write your program without the facility. – High Performance Mark Jun 27 '16 at 15:45
  • Thank you very much for your quick reply Mark. The reason I want to do this is because I think that if it was possible, my code would look much cleaner (without all the `if` statements) and easier to modify. But perhaps it is overkill for what I am doing. I want my program to be able to solve the 1d time-dependent schrodinger equation (TDSE), 2D TDSE, and the array `A` could be the dimension of the wave function. Do you have any references to the solution you suggested with the derived type? – Lun Jun 27 '16 at 15:59

2 Answers2

5

Declare the array to be rank three. If a lower rank array is required, allocate the relevant trailing dimensions to be of size one.

real, allocatable :: array(:,:,:)
...
select case (desired_rank)
case (1) ; allocate(array(n,1,1))
case (2) ; allocate(array(n,n,1))
case (3) ; allocate(array(n,n,n))
case default ; error stop 'bad desired rank'
end select

You can then use an array section to get a contiguous slice of array that is consistent with your desired rank. Alternatively, write the relevant procedures that operate on array to take a rank three argument, and make them aware of the meaning of a size one extent for the higher dimensions.

IanH
  • 21,026
  • 2
  • 37
  • 59
  • Thanks, but not sure this completely solves my problem. I want to have a procedure `work` (contains subroutines that depends on the rank of input array) that when acted on `array`, gives out the desired result. If I use an array section to get a contiguous slice of array, I need to save the array slice before I can use the `work`? If I follow your alternative suggestion, I would have to define `work` such that it recognizes the three different cases, basically shifts the `if` statements of my original post into the subroutine. But maybe this is the only way without resorting to `select rank`. – Lun Jun 27 '16 at 22:07
  • You currently have three procedures behind one generic interface. I don't know what "save the array slice" means. You can take a contiguous rank one array slice using the syntax `array(:,i,j)`. In some cases processing of a higher rank array simply involves iteration over the last rank - the equivalent of the 'if' statements in that case become implicit in the single iteration if this sort of scheme is used for embedding lower rank arrays. If you were not invoking your procedures through a generic interface you could also use sequence association. – IanH Jun 27 '16 at 23:04
  • Thanks again. But suppose I have two procedures (depending on the rank of `array`) behind one generic interface, which outputs an array containing each `array`element multiplied by a number. To do what you suggest I either need to make `if` statements to select the slice, or I need to loop over the outer rank index from the main program and call the procedure (which is not desirable since I want to parallelize the outer loop). If only `array` in `case (1)` of your answer had rank one, then everything would work easily. :) – Lun Jun 28 '16 at 10:22
4

The next Fortran standard (2015) has the select rank construct similar to select case. My example uses the select case construct on the rank intrinsic of an assumed-rank dummy variable. 

    module my_type

  use, intrinsic :: iso_fortran_env, &
       ip => INT32, &
       wp => REAL64

  implicit none
  private
  public :: MyType

  type MyType
     real (wp)              :: rank0
     real (wp), allocatable :: rank1(:)
     real (wp), allocatable :: rank2(:,:)
     real (wp), allocatable :: rank3(:,:,:)
   contains
     procedure :: create => create_my_type
     procedure :: destroy => destroy_my_type
  end type MyType

contains

  subroutine create_my_type(this, array)
    ! calling arguments
    class (MyType), intent (in out) :: this
    real (wp),      intent (in)     :: array(..) !! Assumed-rank dummy variable

    ! local variables
    integer (ip), allocatable :: r(:)

    select case(rank(array))
    case (0)
       return
    case (1)
       r = shape(array)
       allocate( this%rank1(r(1)) )
    case (2)
       r = shape(array)
       allocate( this%rank2(r(1), r(2)) )
    case (3)
       r = shape(array)
       allocate( this%rank3(r(1), r(2), r(3)) )
    case default
       error stop 'array must have rank 0,1,2, or 3'
    end select

    ! Release memory
    if (allocated(r)) deallocate( r )

  end subroutine create_my_type


  subroutine destroy_my_type(this)
    ! calling arguments
    class (MyType), intent (in out) :: this

    if (allocated(this%rank1)) deallocate( this%rank1 )
    if (allocated(this%rank2)) deallocate( this%rank2 )
    if (allocated(this%rank3)) deallocate( this%rank3 )

  end subroutine destroy_my_type

end module my_type

program main

  use, intrinsic :: iso_fortran_env, only: &
       ip => INT32, &
       wp => REAL64

  use my_type, only: &
       MyType

  implicit none

  type (MyType) :: foo
  real (wp)     :: a0, a1(42), a2(42,42), a3(42,42,42)

  print *, rank(a0)
  print *, rank(a1)
  print *, rank(a2)
  print *, rank(a3)

  ! Allocate array of rank 3
  call foo%create(a3)

  print *, rank(foo%rank3)
  print *, shape(foo%rank3)
  print *, size(foo%rank3)

  ! Release memory
  call foo%destroy()

end program main
jlokimlin
  • 593
  • 4
  • 9
  • Yes, one can use assumed-rank arrays in F2015. But as their uses in an entirely Fortran program are very constrained they aren't too good for controlling program flow. That is, the value of `array` in `create_my_type` isn't freely available so one gets much the same control as just passing the literal constant `3` around. – francescalus Jun 27 '16 at 19:08
  • 1
    @francescalus F2015 adds capability beyond the further C interoperability TS. `select rank` has a similar model to `select type` (though it has some differences related to object attributes) - it makes available a name in a child block that has the rank corresponding to the select rank case statement that leads the block. You can then reference/define that object with its specific rank as if it was a normal array. But use of this F201X feature is pretty hypothetical today anyway. – IanH Jun 27 '16 at 20:17
  • Thanks for the extra detail, @IanH. I'll have to make sure my copy of the draft/proposal is up to date. Hopefully my first comment still makes sense as it refers to the example in this answer, rather than the more general `select rank` way. – francescalus Jun 27 '16 at 20:30
  • Thank you for the answer and all the comments. I am quite a newb at Fortran, but since F2015 is quite new, is it true to assume that it will contain bugs an that it won't be available on all systems and clusters? – Lun Jun 28 '16 at 09:32