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