This is the mail archive of the
fortran@gcc.gnu.org
mailing list for the GNU Fortran project.
Is it a bug: data aligment of matrix returned by gemm
- From: Maciej Zwierzycki <Maciej dot Zwierzycki at ifmpan dot poznan dot pl>
- To: fortran at gcc dot gnu dot org
- Date: Thu, 25 Jun 2009 14:38:31 +0200
- Subject: Is it a bug: data aligment of matrix returned by gemm
- Reply-to: Maciej dot Zwierzycki at ifmpan dot poznan dot pl
This is my first post. Greetings everyone!
I am not sure if the following behaviour is a bug or simply the matter
of my poor understanding of F90 standard so I am not submitting this a
a formal bug report.
Consider following code:
function matmul_func(a,b)
real(kind=8),intent(in) :: a(:,:), b(:,:)
real(kind=8) :: matmul_func(size(a,1),size(b,2))
call dgemm('N','N',size(a,1),size(b,2),size(a,2),1.0d0,a,size(a,1),&
&b,size(b,1),0.0d0,matmul_func,size(a,1))
end function
program test
integer, paramater :: range=2
real(kind=8) :: a(rank,rank), b(rank,rank), c(2,rank,rank)
interface
function matmul_func....
....
end function matmul_func
end interface
a(1,:)=1.0d0
a(2,:)=2.0d0
b(1,:)=(/1.0d0,2.0d0/)
b(2,:)=(/3.0d0,4.0d0/)
c(1,:,:)=matmul_func(a,b)
end program test
The result is incorrect in a sense the the resulting
2x2 matrix is put into c coulmnwise as:
c(1,1,1), c(2,1,1), c(1,2,1), c(2,2,1)
whereas the syntax suggests that the c(1,:,:) block should be filled,
i.e.:
c(1,1,1), c(1,2,1), c(1,1,2), c(1,2,2).
Is it a bug? The result can be understood in the above fashion and
corrected for, but it would be nice for gfortran to
honour the block syntax. Intel Fortran and g95 are doing just that.
If I pass the rank of the matrices explicitely, i.e.
function matmul_func(n,a,b)
real... :: a(n,n), b(n,n), matmul_func(n,n)
with the rest of the body of the function as previously the results are
correct.
Replacing geem with: matmul_func=matmul(a,b) yields correct
results. Also when GEMM is called in the main body of the program:
call dgemm(..... c(1,:,:),rank)
the result is correct.
I have used 4.3,4.4 and the 4.5-20090618 versions linking agains
reference netlib BLAS (also recompiled with 4.5) and intel's mkl and
ATLAS. The results are all the same.
I'd appreciate your comments
Regards
MZ