matrix matrix multiplication using BLAS3 routine gemm

eide — Mon, 12 Nov 2012 11:27:28 GMT

Hi,

I have some problems using the gemm routine, please see:

[fortran]

    program mm_test
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! test matrix matrix multiplication
!
! compile
!
!ifort   -w mm_test.f90 \
!        "/opt/intel/composer_xe_2011_sp1.9.293/mkl/lib/intel64"/libmkl_intel_lp64.a \
!        -Wl,--start-group \
!               "/opt/intel/composer_xe_2011_sp1.9.293/mkl/lib/intel64"/libmkl_intel_thread.a \
!               "/opt/intel/composer_xe_2011_sp1.9.293/mkl/lib/intel64"/libmkl_core.a \
!        -Wl,--end-group \
!        -L"/opt/intel/composer_xe_2011_sp1.9.293/mkl/../compiler/lib/intel64" -liomp5 -lpthread \
!        -o x_mm
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   implicit none
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   integer   (kind = 4)               :: i,j
   real      (kind = 4)               :: a(2,3),b(3,2),c(2,2),alpha, beta
   complex   (kind = 8)               :: ca(2,3),cb(3,2),cc(2,2),calpha, cbeta
   character (len = 1)                :: transa,transb
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! real example
   transa = 'N'
   transb = 'N'
   alpha = 1.0
   beta   = 0.0

   a(1,1) = 1.0
   a(1,2) = 2.0
   a(1,3) = 3.0
   a(2,1) = 4.0
   a(2,2) = 5.0
   a(2,3) = 6.0

   b = transpose(a)

   call sgemm(transa, transb, 2, 2, 3, alpha, a, 2, b, 3, beta, c, 2)
   print*,'result for real numbers'
   do i = 1,2
       do j = 1,2
           print*,i,j,c(i,j)
       enddo ! j
   enddo ! i
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! double complex numbers
   transa = 'N'
   transb = 'N'
   calpha = DCMPLX(1.0,0.0)
   cbeta   = DCMPLX(0.0,0.0)

   ca(1,1) = DCMPLX(1.0,0.0)
   ca(1,2) = DCMPLX(2.0,0.0)
   ca(1,3) = DCMPLX(3.0,0.0)
   ca(2,1) = DCMPLX(4.0,0.0)
   ca(2,2) = DCMPLX(5.0,0.0)
   ca(2,3) = DCMPLX(6.0,0.0)

   cb = dconjg(transpose(ca))

   call dzgemm(transa, transb, 2, 2, 3, calpha, ca, 2, cb, 3, cbeta, cc, 2)
   print*,'result for complex numbers'
   do i = 1,2
       do j = 1,2
           print*,i,j,cc(i,j)
       enddo ! j
   enddo ! i
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! double complex numbers test with internal conjugate transpose
   transa = 'N'
   transb = 'C'
   calpha = DCMPLX(1.0,0.0)
   cbeta   = DCMPLX(0.0,0.0)

   ca(1,1) = DCMPLX(1.0,0.0)
   ca(1,2) = DCMPLX(2.0,0.0)
   ca(1,3) = DCMPLX(3.0,0.0)
   ca(2,1) = DCMPLX(4.0,0.0)
   ca(2,2) = DCMPLX(5.0,0.0)
   ca(2,3) = DCMPLX(6.0,0.0)

   call dzgemm(transa, transb, 2, 2, 3, calpha, ca, 2, ca, 2, cbeta, cc, 2)
   print*,'result for complex numbers test'
   do i = 1,2
       do j = 1,2
           print*,i,j,cc(i,j)
       enddo ! j
   enddo ! i

   end program mm_test

[/fortran]

it produces the following screen output

[plain]

result for real numbers
           1           1   14.00000
           1           2   32.00000
           2           1   32.00000
           2           2   77.00000
result for complex numbers
           1           1 (15.0000000000000,0.000000000000000E+000)
           1           2 (36.0000000000000,0.000000000000000E+000)
           2           1 (0.000000000000000E+000,0.000000000000000E+000)
           2           2 (0.000000000000000E+000,0.000000000000000E+000)
result for complex numbers test
           1           1 (15.0000000000000,0.000000000000000E+000)
           1           2 (36.0000000000000,0.000000000000000E+000)
           2           1 (0.000000000000000E+000,0.000000000000000E+000)
           2           2 (0.000000000000000E+000,0.000000000000000E+000)

[/plain]

I wonder why the resulting complex matrix is not symmetric (like in the real number example), although it should be? Can anybody help me?

Thank you in advance.

Eide

You are not calling dzgemm

mecej4 — Mon, 12 Nov 2012 21:11:19 GMT

You are not calling dzgemm with arguments of the specified types. To catch this kind of error, add the following statement to your source code and the compiler will tell you what the errors are:
[fortran] include 'mkl.fi'[/fortran]

Hi mecej,

eide — Mon, 12 Nov 2012 21:40:33 GMT

Hi mecej, thank you very much. The first matrix has to be double precision - now it is correct. Thanks for the tip with the include file. Cheers Eide

topic matrix matrix multiplication using BLAS3 routine gemm in Intel® oneAPI Math Kernel Library

matrix matrix multiplication using BLAS3 routine gemm

You are not calling dzgemm

Hi mecej,