DGEMM is not working

[kon-konstantinov]

Dear all,

I amd trying to use DGEMM subroutine from BLAS to multiply two large matrices.
However my test code is not working. What am I doing wrong? Thanks in advance.

The code:
program test_blas
use mkl95_precision, only: wp => dp ! in program head
USE mkl95_LAPACK, ONLY: GETRF, GETRI, SYTRF, SYTRI, POTRF, POTRI
USE mkl95_BLAS, ONLY: GEMM

!implicit none
integer :: i, j, k, l, m, n, lda, ldb, ldc
real(8), DIMENSION(:,:), ALLOCATABLE :: A, B, C, S
character(len = 1) :: transa, transb
real(8)::alfa,beta


transa = 'N'; transb = 'N'
m = 2
n = 2
k = 2
alfa=1.d0; beta=1.d0

if ((transa.eq.'N').or.(transa.eq.'n')) then
lda = m
allocate(a(m, k))
else
lda = k
allocate(a(k, m))
end if
if ((transb.eq.'N').or.(transb.eq.'n')) then
ldb = k
allocate(b(k, n))
else
ldb = n
allocate(b(n, k))
end if
ldc = m
allocate(c(m, n))


print*, " Dimensions of A : ", ubound(a)
print*, " Dimensions of B : ", ubound(b)
print*, " Dimensions of C : ", ubound(c)

!A(1,1) = 1; A(1,2) = 2; A(1,3) = 3
!A(2,1) = 4; A(2,2) = 5; A(2,3) = 6
!
!B(1,1) = 2; B(1,2) = 3
!B(2,1) = 4; B(2,2) = 5
!B(3,1) = 6; B(3,2) = 7

A(1,1) = 1; A(1,2) = 2
A(2,1) = 4; A(2,2) = 5
!
B(1,1) = 2; B(1,2) = 3
B(2,1) = 4; B(2,2) = 5


print*, " Original A : "
do i = 1, m
write(*,fmt='(13f10.3)') a(i,:)
enddo
print*, " Original B : "
do i = 1, n
write(*,fmt='(13f10.3)') b(i,:)
enddo

! Call DGEMM subroutine
! call DGEMM(a, b, c, transa, transb, alpha, beta)

call dgemm(transa, transb, m, n, k, alfa, a, lda, b, ldb, beta, c, ldc)


Print*, " Result matrix C = AxB "
do i = 1, ldc
write(*,fmt='(i8,13f10.3)') i,c(i,:)
enddo


end program test_blas

The compilation:

/opt/intel/bin/ifort construct_G.f90 -L/opt/intel/composerxe-2011.2.137/mkl/lib/intel64 -I/opt/intel/composerxe-2011.2.137/mkl/include/intel64/ilp64 -lmkl_blas95_ilp64 -lmkl_lapack95_ilp64 -Wl,--start-group /opt/intel/composerxe-2011.2.137/mkl/lib/intel64/libmkl_intel_ilp64.a /opt/intel/composerxe-2011.2.137/mkl/lib/intel64/libmkl_intel_thread.a /opt/intel/composerxe-2011.2.137/mkl/lib/intel64/libmkl_core.a -Wl,--end-group -openmp -lpthread -o construct_G

The result:

Dimensions of A : 2 2
Dimensions of B : 2 2
Dimensions of C : 2 2
Original A :
1.000 2.000
4.000 5.000
Original B :
2.000 3.000
4.000 5.000

So I am getting nothing.

Kon

[mecej4]

You link with the ILP64 libraries, yet I do not see anything in your program or the compilation command to tell the compiler that default integers should be treated as eight-byte integers. Your USE statements did not do one of their intended functions since the MKL subroutine that you call is not declared in those modules.

You can do one of:

(i) Use LP64 libraries

(ii) Use ILP64 libraries and change the declarations of the integer variables involved in the calls to 8-byte integers.

(iii) Use ILP64 libraries and use a compiler option that promotes default integers to 8-byte integers.

(iv) Change the DGEMM call to a call to the generic GEMM subroutine.

[kon-konstantinov]

Thank you very much for your suggestions. I've change the dafault integers to 8-byte integers
and everything worked. As per using "USE mkl95_BLAS, ONLY: GEMM", changing to DGEMM is not working.
It requires the generic name gemm. Anyway thanks again.

Kon