How to have an allocatable array of variable rank (number of dimensions)?

Question

I have two programs : 2D and 3D versions. For example, this is a simplified 2D version :

program main

    integer, parameter :: nDim = 2
    integer :: size
    real, dimension(:,:,:), allocatable :: myArray

    size=4 ! Normally read from a file
    allocate(myArray(1:size,1:size,1:nDim))
    call initArray(myArray)

contains
    
    subroutine initArray(myArray)
    
        real, dimension(1:size,1:size,1:nDim) :: myArray
    
        myArray(1:size,1:size,1:nDim) = 5.d0
     
        return  
    end subroutine initArray
    
end program main

And the 3D version,

program main

    integer, parameter :: nDim = 3
    integer :: size
    real, dimension(:,:,:,:), allocatable :: myArray

    size=4 ! Normally read from a file
    allocate(myArray(1:size,1:size,1:size,1:nDim))
    call initArray(myArray)

contains
    
    subroutine initArray(myArray)
    
        real, dimension(1:size,1:size,1:size,1:nDim) :: myArray
    
        myArray(1:size,1:size,1:size,1:nDim) = 5.d0
     
        return  
    end subroutine initArray
    
end program main

These programs are very similar and I would like to have only one program where the parameter nDim determine everything.

But I have problem with the following statements and instructions :

1. For dimension, I have not the same number of dimensions (3 or 4)
2. For allocate, the number of arguments is variable (3 or 4)
3. For initialisating myArray, the number of dimension is variable (3 or 4)

Is there a solution purely in fortran or should I use the C preprocessor ?

Thanks for answer.

Answer 1

This is a typical example where object-oriented programming helps. You have two programs that respond to the same API, but internal computations will be different (2D or 3D). You want to have 2d and 3d grids both extend a generic "grid type"

module grids
   implicit none

   type, abstract, public :: grid
      contains
      procedure (grid_init), deferred :: init
   end type grid

   type, public, extends(grid) :: grid2D
      real, allocatable :: myArray(:,:)
      contains
         procedure :: init=>init_2d
   end type grid2D
   
   type, public, extends(grid) :: grid3D
      real, allocatable :: myArray(:,:,:)
      contains
         procedure :: init=>init_3d
   end type grid3D

   ! Define procedure APIs 
   abstract interface
      subroutine grid_init(this,size)
         import grid
         class(grid), intent(inout) :: this
         integer, intent(in) :: size
      end subroutine grid_init
   end interface

   contains

   ! Define ACTUAL procedures for the 2d and 3d grids
   subroutine init_3d(this, size)
      class(grid3D), intent(inout) :: this
      integer, intent(in) :: size
      allocate(this%myArray(size, size,size))
   end subroutine init_3d
   subroutine init_2d(this, size)
      class(grid2D), intent(inout) :: this
      integer, intent(in) :: size
      allocate(this%myArray(size,size))
   end subroutine init_2d
end module grids

The actual implementation is up to you, but the key points are:

• Expose all common "operations" (on a grid, for example) through the polymorphic API; define them in the abstract class, so the compiler will check that all you've done in the child classes is right;

• Hide all dimension-dependent code and data inside the derived types;

• Ideally, you would try to define all routines of your public interface in a dimension-independent way. This is often not possible, in which case, you'll revert to use type-checking in those places where a shared call is not possible.

For example, you may need 2 or 3 input parameters to something: you can do

class(grid), allocatable :: myGrid

select type (myGrid)
   type is (grid2D); [...]
   type is (grid3D); [...]
   class default; [...]
end select

A simple test program will look like

program test_grids
   use grids
   implicit none

   type(grid2D) :: fixed_2d
   type(grid3D) :: fixed_3d
   class(grid), allocatable :: polymorphic

   ! Initialize non-polymorphic grids (fixed type)
   call fixed_2d%init(5)
   call fixed_3d%init(10)

   ! Initialize polymorphic grid as 2D
   allocate(grid2D :: polymorphic)
   call polymorphic%init(10)

   ! Initialize polymorphic grid as 3D
   deallocate(polymorphic)
   allocate(grid3D :: polymorphic)
   call polymorphic%init(10)

end program test_grids

Hope this helps, Federico

Answer 2

I know, it is ugly but it works and it is convenient to use

#define nDim 3

#if nDim==2
#  define xdimCpp :,:
#  define xboundsCpp(X) 1:X,1:X
#elif nDim = 3
#  define xdimCpp :,:,:
#  define xboundsCpp(X) 1:X,1:X,1:X
#endif
  
program main

    integer :: size
    real, dimension(xdimCpp,:), allocatable :: myArray

    size=4 ! Normally read from a file
    allocate(myArray(xboundsCpp(size),1:nDim))
    call initArray(myArray)
    print *,myArray
    
contains
    
    subroutine initArray(myArray)
    
        real, dimension(xboundsCpp(size),1:nDim) :: myArray
    
        myArray(xboundsCpp(size),1:nDim) = 5.d0
     
        return  
    end subroutine initArray
    
end program main