Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.
6881 Discussions

## matrix matrix multiplication using BLAS3 routine gemm

Beginner
213 Views

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_core.a \
!        -Wl,--end-group \
!        -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?