Fast 3D array Multiplication with DGEMM.....help

[Zioivan]

Hello!

I wrote this routine that takes two 3D arrays and performs a multiplication.......

function matmull(a,b)
! This in-line function performs array muliplication of two 3-D arrays
double precision :: a,:,,b,:,
double precision :: matmull(size(a,1),size(b,3))
integer i

matmull=0.0
do i=1,size(b,2)
matmull=matmull + matmul(a,:,i),b,i,)
enddo
end function matmull

I'd like to use the BLAS routine DGEMM that should be faster than the previous one, but it doesn't works:

function matmull(a,b)
! This in-line function performs array multiplication of two 3-D arrays
double precision :: a,:,,b,:,
double precision :: matmull(size(a,1),size(b,3))
integer :: M,N,K
M=size(a,1)
N=size(b,3)
K=size(a,3)*size(a,3)
CALL DGEMM('N','N',M,N,K,1.0,a,M,b,K,0.0,matmull,M)
end function matmull

Where is the mistake? I think there is a problem with the index of a but I'm not shure......

Thank you for the help.
I

 

[salgerman]

The mistake, as far as I can tell, is that matrix multiplication is only defined between two 2-dimensional matrices; but you are passing 3-dimensional matrices to DGEMM.

If you want to multiply 3-dimensional matrices, you need to stick to your own custom multiplication function...or keep passing one plane at a time...but you need to make sure the plane you are passing truly behaves (and is stored in memory) in the exact same way as a 2-d matrix.

 

[Zioivan]

Thank you!

If I put the DGEMM function in a do loop:

do i=1,size(b,2)
CALL DGEMM('N','N',M,N,K,1.0,a,:,i),M,b,i,,K,0.0,atimesb,M)
matmull=matmull + atimesb
enddo

the results are the same of the "standard" implementation but is slower! Is it possible to reshape the matrixs before in order to obtain a 2d array?

Thanks

 

[salgerman]

Never tried it; it sounds nonsensical but you may onto something...you can either work out the algebra to prove it right or wrong, which would require you to know where the elements end up when you re-shape a 3-d array into a 2-d one or you could simple write a small fortran program give it a try.

 

[Zioivan]

Thanks. Just to clarify the following function performs array multiplication of two 3-D arrays. In tensor notation, it performs a_{ijk}b_{jkl}:

function matmull(a,b)
double precision :: a,:,,b,:,
double precision :: matmull(size(a,1),size(b,3))
matmull=matmul(reshape(a,(/size(a,1),size(a,2)*size(a,3)/)),reshape(b,(/size(b,1)*size(b,2),size(b,3)/)))
end function matmull

I'd like to use the DGEMM instead of matmul.........