CODEWITHSUNDEEP
fortran
Monica
Monica

Reputation: 1070

Reducing 4 dimentional array to 2 Fortran

I am trying to reshape 4-dimensional array in Fortran to 2 dimensional one:

do k=1,nA
   do m=1,nA
      do l=1,nB
         do n=1,nB

           p = k * nA + m
           q = l * nB + n

           Lpq(p,q) = Aij(k,m,l,n) 

         end do
      end do
   end do
end do

What is the proper way of defining indexes p and q?

Upvotes: 2

Views: 100

Answers (2)

HanT
HanT

Reputation: 151

It depends of what you want to do with Aij and Lpq. If you just want to change the indexing to have easier access to the array elements and want to use Lpq as a read-only array, you might well use the EQUIVALENCE statement. Moreover, this will save memory space, because it only re-interprets the array. If you want to treat them as independ arrays you can make a copy of Lpq next, and use the copy.

  integer Aij(nA,nA,nB,nB)
  integer Lpq(nA*nA,nB*nB), Lpq_copy(nA*nA,nB*nB)
  equivalence (Aij, Lpq)

! fill Aij with test data
  do k=1,nA
  do m=1,nA
  do l=1,nB
  do n=1,nB
    Aij(k,m,l,n) = k*10000+m*1000+l*100+n
  end do
  end do
  end do
  end do

  print *, Lpq

! make a copy if needed
  Lpq_copy = Lpq
  Lpq(1,1) = 1 ! this will affect Aij, but not Lpq_copy

Upvotes: 3

kvantour
kvantour

Reputation: 26471

The following answer makes the following assumption:

  • Both Aij and Lpq have default start and end indices (i.e. start indices are ONE and end indices are the corresponding dimensions. I.e.

    TYPE, DIMENSION(nA,nA,nB,nB) :: Aij
TYPE, DIMENSION(nA*nA,nB*nB) :: Lpq

There are many ways you can map Aij on Lpq, but the one you defined is where the second and fourth index in Aij (referenced as m and n) are the fastest running indices in Lpq. This means:

Aij(1,1  ,1,1) => Lpq(1,1)
...
Aij(1,nA ,1,1) => Lpq(nA,1)
Aij(2,1  ,1,1) => Lpq(nA+1,1)
...
Aij(2,nA ,1,1) => Lpq(2*nA,1)
...
Aij(nA,nA,1,1) => Lpq(nA*nA,1)

A similar mapping exists for the other dimensions.

The current mapping you implemented (i.e. p = k * nA + m) sets p=nA+1 when k=m=1 and returns p=(nA+1)*nA for k=m=nA. So you are missing the lower range and overshoot the upper range.

you have to define the mapping as:

p = (k-1) * nA + m
q = (l-1) * nB + n

Note, however, as mentioned in the comments, Fortran has a intrinsic procedure reshape that could do these kind of mappings for you. Sadly, in your particular case it is not as handy due to the choice of the fastest index.

  • If you would choose the fastest indices to be k and l, then you can use a single reshape as the system is already in column-major-order:

    Lpq = reshape(Aij,[nA*nA,nB*nB])

  • If you choose the fastest indices to be m and n, then you need two calls to reshape. The first to swap two indices and bring the system in a column manjor order, and the second to do the real reshape.

    Bij = reshape(Aij,[nA,nA,nB,nB],order=[2,1,4,3])
Lpq = reshape(Bij,[nA*nA,nB*nB])

Upvotes: 3

Related Questions