CODEWITHSUNDEEP
fortrangfortranblas
Emil Andre Løvbak
Emil Andre Løvbak

Reputation: 41

Unexpected results from SGEMM

I'm trying to replace the following instance of matmul in the code calculating A*B':

print *, mat%A
print *, mat%B
print *, mat%dimsA
print *, mat%dimsB
C = matmul(mat%A,transpose(mat%B))
print *, C

with SGEMM:

print *, mat%A
print *, mat%B
print *, mat%dimsA
print *, mat%dimsB
call SGEMM('N', 'T', mat%dimsA(1), mat%dimsB(1), mat%dimsA(2), 1, mat%A, &
    mat%dimsA(1), mat%B, mat%dimsB(1), 0, C, mat%dimsA(1))
print *, C

Where mat is an type that contains the following relevant data:

integer :: dimsA(2), dimsB(2)
real, pointer :: A(:,:),B(:,:)

Running the code on a small test where A is:

1 1 1 
2 2 2

and B is:

1 1 1

I get the following outputs for matmul:

1.00000000       2.00000000       1.00000000       2.00000000       1.00000000       2.00000000    
1.00000000       1.00000000       1.00000000    
2           3
1           3
3.00000000       6.00000000

and SGEMM:

1.00000000       2.00000000       1.00000000       2.00000000       1.00000000       2.00000000    
1.00000000       1.00000000       1.00000000    
2           3
1           3
4.20389539E-45   8.40779079E-45

Am I wrong in assuming these two blocks of code should have the same result?

Upvotes: 1

Views: 158

Answers (1)

Alexander Vogt
Alexander Vogt

Reputation: 18118

This is because you are passing integers for alpha and beta to the lapack subroutine. The compiler does not have the interface of SGEMM available and cannot warn you of the resulting type confusion...

Just use the notation for (default precision) reals - 1. instead of 1, and 0. instead of 0:

call SGEMM('N', 'T', mat%dimsA(1), mat%dimsB(1), mat%dimsA(2), 1., mat%A, &
    mat%dimsA(1), mat%B, mat%dimsB(1), 0., C, mat%dimsA(1))

Upvotes: 1

Related Questions