CODEWITHSUNDEEP
multidimensional-arrayfortranstorage
Enlico
Enlico

Reputation: 28416

store multi dimensional arrays to be operated on along all the dimensions

Foreword

The Fortran program that I'm writing should deal with 1D, 2D and 3D problems depending on ndims, which can be 1, 2 or 3 and is read from an input file.

In these cases the quantity/ies of interest can be stored in arrays (one could be named phi)

  1. of rank dims (ALLOCATABLE(:) or ALLOCATABLE(:,:) or ALLOCATABLE(:,:,:)),
  2. or in arrays of rank 3 (ALLOCATABLE(:,:,:) with third dimension to be set equal to 1 in 2D or both second and third dimensions equal to 1 in 1D);

both cases are explained well in this answer. The first approach seems more elegant to me, but in the following I assume the second one, which is definitely simpler.

These quantities have to be operated on by several subroutines (e.g. mysub) along the ndims dimensions (along "pencils" should give a graphic idea), so I should call something like

SELECT CASE (ndims)

! 3D case
CASE (3)
  DO j = ...
    DO k = ...
      CALL mysub(phi(:,j,k))
    END DO
  END DO
  DO i = ...
    DO k = ...
      CALL mysub(phi(i,:,k))
    END DO
  END DO
  DO i = ...
    DO j = ...
      CALL mysub(phi(i,j,:))
    END DO
  END DO

! 2D case
CASE (2)
  DO j = ...
    DO k = ...
      CALL mysub(phi(:,j,1))
    END DO
  END DO
  DO i = ...
    DO k = ...
      CALL mysub(phi(i,:,1))
    END DO
  END DO

! 1D case
CASE (1)
  DO j = ...
    DO k = ...
      CALL mysub(phi(:,1,1))
    END DO
  END DO
END SELECT

Actual question

Can anyone suggest me (or help me to to devise!) a different way of store phi (maybe involving derived data types?) so that I can collapse the preceding code as follows?

DO id = 1, ndims
  CALL mysub2(phi,id)
END DO

(Here mysub2 acts on mysub's place.)

So the question is how should I store phi, so that I can substitute the first code with the second one?

Maybe I could return to the foreword and decide to follow the point 1., in which case would be easier to write a generic interface. This would be, I think, just a way to "hide" exactly what the SELECT CASE would do. Which of the two (SELECT CASE/generic INTERFACE) would be more efficient?

Are these the only two ways to face this problem?

Upvotes: 0

Views: 156

Answers (2)

IanH
IanH

Reputation: 21431

Perhaps I have misunderstood, but I think the answer to the specific question is to not make any changes to the storage or declaration of phi at all.

In the original code, three dimensional data (differentiating rank of the data from rank of the array used to store the data) is processed in slices along the first dimension, then the second, then the third. Two dimensional data is processed along the first, then the second, and one dimensional data is processed along the first only.

So with id going from 1 to the number of dimensions in the data, consider the following implementation of mysub2:

SUBROUTINE mysub2(phi, id)
  TYPE(pink_elephant), INTENT(IN) :: phi(:,:,:)
  INTEGER, INTENT(IN) :: id

  INTEGER :: i, j, k

  SELECT CASE (id)
  CASE (1)
    DO j = ...
      DO k = ...
        CALL mysub(phi(:,j,k))
      END DO
    END DO

  CASE (2)
    DO i = ...
      DO k = ...
        CALL mysub(phi(i,:,k))
      END DO
    END DO

  CASE (3)
    DO i = ...
      DO j = ...
        CALL mysub(phi(i,j,:))
      END DO
    END DO

  END SELECT
END SUBROUTINE mysub2

~~

Generic interfaces can always be resolved at "compile time" - the specific procedure (not type bound) or binding (type bound) that will be invoked by a particular CALL statement or function reference can be determined just from looking at the declarations in the code.

If you have a situation where "runtime" information is going to affect the choice of procedure, then there has to be some other executable mechanism, other than or additional to the resolution of a generic, that comes into play - if statement, select case, dynamic dispatch, etc, etc, etc.

Asking whether generic resolution is more efficient than a executable decision is therefore not particularly meaningful - they are different things.

Upvotes: 1

PeMa
PeMa

Reputation: 1696

You probably want something like this:

program test

   integer :: j,ndims
   integer :: n ! rank of each dimension, could also be read from input an allocated separately

   type arr
      real(8) :: x(n) ! one array for each dimension
   end type

   type(arr),allocatable :: phi

   read(*,*) ndims
   allocate(phi(ndims))

   do j=1,ndims
      call mysub(phi(j)%x) ! acts on the array in dimension j
   end do

contains

   subroutine mysub(x)
   ...
   end subroutine

end program

Upvotes: 1

Related Questions