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 :

For dimension, I have not the same number of dimensions (3 or 4)
For allocate, the number of arguments is variable (3 or 4)
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.

The situation hasn't really changed since this previous Q&A - https://stackoverflow.com/questions/38328846/dynamic-dimensions-of-arrays/ - similar questions continue to pop up with some regularity. — High Performance Mark, May 28 '22 at 11:52
@HighPerformanceMark. So cpp is the solution ? Or is there anohter one ? — Stef1611, May 28 '22 at 12:06
If you follow the chain of answers I pointed you at and read them you'll see there are a number of Fortran approaches to solving your problem. I don't know what you mean by *cpp* so have no idea if it is a solution or not. And if you read my comment to the first Q I reference, you might understand that I am not persuaded of the *need* for variable rank arrays. And I'm commenting precisely because I am not planning to answer. — High Performance Mark, May 28 '22 at 12:09
@HighPerformanceMark. I posted an answer using cpp. I know it is ugly, ... — Stef1611, May 28 '22 at 13:26

Federico Perini · Accepted Answer · 2022-05-30T08:26:30.927

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

@federicp. Thanks for your answer. I am discovering polymorphism in fortran (abstract interface, abstract type ...). I still have a problem with your code. `call grid_init(myGrid,5` returns an error (myGrid is `class(grid2D), allocatable :: myGrid`)? Could you give me an complete example ? — Stef1611, May 30 '22 at 07:21
If you have `class(grid2D), allocatable :: myGrid` then your variable can only be `grid2D` or any inherited types, so you want to define it with the root class, i.e., `class(grid), allocatable :: myGrid`. Plus, when you have a polymorphic object, you always have to call routines through the type interface, e.g.: `call grid%init(...)` and not the non-polymorphic versions like `call init_2d(....)`. Only with the first method there is polymorphism — Federico Perini, May 30 '22 at 07:40
Not sure to well understand. Now, it is compiling but there is a segmentation fault. I suppose, if I want to work with 2D (or3D) arrays, I declare then as grid2D (or grid3D) in the main ? — Stef1611, May 30 '22 at 07:56
The main program I used is very simple : `class(grid2D), allocatable :: myGrid`, `call myGrid%init(5)` — Stef1611, May 30 '22 at 08:01
Thanks a lot. I discovered my problem just before you sent this example. I declared my variable with `class` and not `type`. I got the explanation here https://stackoverflow.com/a/29652905/7462275 — Stef1611, May 30 '22 at 08:38

score 0 · Answer 2 · answered May 28 '22 at 13:25

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

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

2 Answers2

Linked