3D array data storage 2

[dibas]

Hi.
I'm working with 1D and 2D f90 arrays. Actually, for a fixed value of a
looping/iteration index, I have to form a 2D (dimension(6,6) )
array that holds some Jacobian matrix. Calling this iteration counter i,
where i ranges from 0 to n (say n=5000), this means I need to compute
the 6x6 Jacobian matrices (n+1) times.

In order to store ALL the Jacobian matrices, I have thought of defining
a new 3D array - dimension(1:6, 1:6, 0:n) . To fill out the 3D array, I
want to use one of the two ways shown in the following piece/extract of
code. What I then want to know is: (i) are these valid ways of doing what
I mean to do? (ii) if yes, which one would be more preferable?

In subsequent stages of my work, I will need the 6x6 2D jacobians per
slice of i. I would be happy if I could also get suggestions on how to
retrieve that. Would I need the "reshape" intrinsic function for that?

..........
 real, dimension(-1:n+1)  :: u1, u2, u3, u4, u5, u6
 real, dimension(-1:n+1)  :: fun1, fun2, fun3, fun4, fun5, fun6
 real, dimension(1:6, -1:n+1)  :: fun_tmp, u_tmp
 real, dimension(1:6, 0:n)     :: dfun, du
 real, dimension(1:6,1:6)       :: jac_slice
 real, dimension(1:6, 1:6, 0:n) :: jac
 integer  :: i, k  

 
! Fill fun_tmp by rows
   fun_tmp(1, -1:n+1) = fun1    ; fun_tmp(2, -1:n+1) = fun2
   fun_tmp(3, -1:n+1) = fun3    ; fun_tmp(4, -1:n+1) = fun4
   fun_tmp(5, -1:n+1) = fun5    ; fun_tmp(6, -1:n+1) = fun6

! Fill up u_tmp by rows
   u_tmp(1, -1:n+1) = u1    ; u_tmp(2, -1:n+1) = u2
   u_tmp(3, -1:n+1) = u3    ; u_tmp(4, -1:n+1) = u4
   u_tmp(5, -1:n+1) = u5    ; u_tmp(6, -1:n+1) = u6

! Compute central differences for fun_tmp and u_tmp
   dfun = fun_tmp(:, 1:) - fun_tmp(:, :n-1)  ! 6x(n+1) - 2D array
   du   = u_tmp(:, 1:) - u_tmp(:, :n-1)      ! 6x(n+1) - 2D array

!!! Now compute & store 6x6x(n+1) 3D array---The part I need HELP on!!!

! Option 1: Loop over all Jacobians for all "slices" of i in one stage
  do i = 0, n
     do k = 1, 6            ! loop over columns of "jac" matrix
       jac(:, k, i) = dfun(:, i) / du(k, i)   ! Fill matrix column-wise   
     end do
  end do

! Option 2: Loop over ONE 6x6 2D Jacobian per slice of i (internal loop),
!           then store it & proceed to next slice
  do i = 0, n
     do k = 1, 6              
      jac_slice(:, k) = dfun(:, i) / du(k, i)   
     end do
     jac(1:6, 1:6, i) = jac_slice(1:6, 1:6)
  end do
..............

Your help will be highly appreciated.

 

[gummibaer]

Sorry, where is the problem ?

You need n matrices dimensioned 6 by 6.

So to save them to an array declared as

real jac (6, 6, n)

is the perfectly natural thing to do. Personally I would write the loops differently for more clarity, but this is my own thing. My code would look like

do i = 1, n
    do k1 = 1, 6
        do k2 = 1, 6
            jac (k1, k2, i) = dfun (k2, i) / du (k1 , i)   ! if this is what you want to do
        enddo
    enddo
enddo

Just to keep in control of the flow of evaluation and for better debugging I would introduce this inner k2-loop.

To retrieve a certain matrix you set the third index of jac to the proper value and you can work like

l = (some integer)
do k1 = 1 : 6
    do k2 = 1 : 6
       result = some_maths (jac (k1, k2, l))
    enddo
enddo

Norbert


The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.

 

[dibas]

Thank you for the assurance and extra tip, gummibaer.
When I said "retrieve", I actually meant to say "access".
I need to access the 2D arrays (6x6 matrices) for each i,
and use it as the matrix "A" in Ax=b and solve a system
of equations. In the point you made concerning this, would
I be right to say the variable "result" is a 2D array (6x6
matrix), because that's how it appears to me?

 

[salgerman]

Hhhmmm...I am not entirely sure of what you are trying to do, but be careful how you use such a '2D' matrix...if you are about to pass such matrix to a subroutine that is not yours and that expects a 2D matrix...the individual elements may not quite be lined up as expected...

without testing, I would simply actually copy the 6x6 matrix out of the 3D matrix and THEN pass the actual 2D matrix to the 2D system of equation solver...

...maybe things will be o.k., just something to be aware and consider.

my 2 cents.

 

[gummibaer]

dibas,
I just tried to indicate that you could do anything in this loop with the matrix and produce any result you want. If you define 'result' to be a scalar in type declaration, then you can perform 'some_maths' that result in a scalar, if you declare result a matrix you can do 'some_maths that yields a matrix.

By the way arrays are saved in fortran, even passing the whole 6 x 6 x l matrix to a procedures expecting a 6 x 6 matrix should work - if you observe the way arrays are saved to memory: The first index is run through its complete range, then the second is increased. An after the second has run through all its range the third is increased.

So the l, as I named it must (!) be the last of the three indices and you should (I said 'should') be all right. But there is no harm in testing this.

But your linker might complain for the different ranks of matrices anyway, so as Salgerman said, copying the slice you need to a 6 by 6 matrix to pass it as actual argument is the safest to do.

Norbert


The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.