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...