This is a follow up of this question.

I have an array D(:,:,:) of size NxMxM. Typically, for the problem that I am considering now, it is M=400 and N=600000 (I reshaped the array in order to give the biggest size to the first entry).

Therefore, for each value l of the first entry, D(l,:,:) is an MxM matrix in a certain basis. I need to perform a change of components of this matrix using a basis set vec(:,:), of size MxM, so as to obtain the matrices D_mod(l,:,:).

I think that the easiest way to do it is with:

         D_mod=0.0d0
          
         do b=1,M
         do a=1,M

         do   nu=1,M
         do   mu=1,M
             D_mod(:,mu,nu)=D_mod(:,mu,nu) vec(mu,a)*vec(nu,b)*D(:,a,b)
         end do
         end do

         end do
         end do

Is there a way to improve the speed of this calculation (also using LAPACK/BLAS libraries)?

I was considering this approach: reshaping D into a N x M^2 matrix D_re; doing the tensor product vec(:,:) x vec(:,:) and reshaping it in order to obtain an M^2 x M^2 matrix vecsq_re(:,:) (this motivates this question); finally, doing the matrix product of these two matrices with zgemm. However, I am not sure this is a good strategy.

EDIT I am sorry, I wrote the question too fast and too late. The size can be up to 600000, yes, but I usually adopt strategies to reduce it by a factor 10 at least. The code is supposed to run on nodes with 100 Gb of memory

CodePudding user response：

As @IanBush has said, your D array is enormous, and you're likely to need some kind of high-memory machine or cluster of machines to process it all at once. However, you probably don't need to process it all at once.

Before we get to that, let's first imagine you don't have any memory issues. For the operation you have described, D looks like an array of N matrices, each of size M*M. You say you have "reshaped the array in order to give the biggest size to the first entry", but for this problem this is the exact opposite of what you want. Fortran is a column-major language, so iterating across the first index of an array is much faster than iterating across the last. In practice, this means that an example triple-loop like

do i=1,N
  do j=1,M
    do k=1,M
      D(i,j,k) = D(i,j,k)  1
    enddo
  enddo
enddo

will run much slower ¹ than the re-ordered triple-loop

do k=1,M
  do j=1,M
    do i=1,N
      D(i,j,k) = D(i,j,k)  1
    enddo
  enddo
enddo

and so you can immediately speed everything up by transposing D and D_mod from N*M*M arrays to M*M*N arrays and rearranging your loops. You can also speed everything up by replacing your manually-written matrix multiplication with matmul (or BLAS/LAPACK), to give

do i=1,N
  D_mod(:,:,i) = matmul(matmul(vec , D(:,:,i)),transpose(vec))
enddo

Now that you're doing matrix multiplication one matrix at a time, you can also find a solution for your memory usage problems: instead of loading everything into memory and trying to do everything at once, just load one D matrix at a time into an M*M array, calculate the corresponding entry of D_mod, and write it out to disk before loading the next matrix.

¹ if your compiler doesn't just optimise the loop order.