- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

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

Link Copied

2 Replies

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

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.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

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

Topic Options

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page