Discrepancy between LAPACK Dgemm and Gonum's Dgemm

[Patricio Whittingslow]

I'm trying to understand why two programs that I believe should be doing the same thing are having different effect.

https://github.com/soypat/dgemm-bug

Here's the Go program: https://github.com/soypat/dgemm-bug/blob/main/gover/main.go

Fortran LAPACK program: https://github.com/soypat/dgemm-bug/blob/main/main.f90

[Vladimír Chalupecký]

It's the storage order difference between Fortran and Go. rhs is a 1d array [ 2  1  -1  3 ... ] given to this particular DGEMM call as a 2x2 matrix. Fortran goes column-wise, so the rhs matrix is [ 2 -1; 1 3 ], while Go goes row-wise and so the rhs matrix is [ 2 1; -1 3 ].

[Patricio Whittingslow]

Holy moly, thank you. This makes so much sense now. Right, so to recap, one must take great care when an 1D array is passed as a matrix in Fortran since this means elements will be accessed fundamentally differently. Luckily in this case it can easily be solved by setting the B matrix argument to transposed.

Thank you again Vladimir!

[Patricio Whittingslow]

I have deleted the repository but maintained a gist for posterity's sake.
https://gist.github.com/soypat/8dd5bca0a99a3defd858e58ea2af4946

[Brendan Tracey]

Note that in general it's not as simple as setting the transpose flag. In particular, if the stride is longer than the number of columns, the data will be read incorrectly

--
You received this message because you are subscribed to the Google Groups "gonum-dev" group.
To unsubscribe from this group and stop receiving emails from it, send an email to gonum-dev+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/gonum-dev/6fd604e1-f0a7-4a1c-950a-e686f71572efn%40googlegroups.com.

[jorgeneo560]

i had a similar problem and i dont know how to solve it, in fortran i do

input

 predgemm matrix a
            
       4.0000     1.0000     1.0000     2.0000
       0.0000     3.0000     4.0000     1.0000
       0.0000     1.0000     3.0000     1.0000
       0.0000     0.0000     0.0000     6.0000
 predgemm matrix c
           
      -4.0000     7.0000     1.0000    12.0000
      -9.0000     2.0000    -2.0000    -2.0000
      -4.0000     2.0000    -2.0000     8.0000
      -1.0000     1.0000    -1.0000    19.0000

is=4

nb=3

mb=1

ldc=lda=4

CALL dgemm2( 'N', 'N', is-1, nb, mb, -one,a( 1, is ), lda, c( is, js ), ldc, one,c( 1, js ), ldc) 

 dgemm matrix a
             1          2          3          4
 1      4.0000     1.0000     1.0000     2.0000
 2      0.0000     3.0000     4.0000     1.0000
 3      0.0000     1.0000     3.0000     1.0000
 4      0.0000     0.0000     0.0000     6.0000
 dgemm matrix c
             1          2          3          4
 1     -2.0000     5.0000     3.0000    12.0000
 2     -8.0000     1.0000    -1.0000    -2.0000
 3     -3.0000     1.0000    -1.0000     8.0000
 4     -1.0000     1.0000    -1.0000    19.0000

in go

same values

[]float64 len: 16, cap: 16, [4,0,0,0,1,3,1,0,1,4,3,0,2,1,1,6]

[]float64 len: 16, cap: 16, [-4,-9,-4,-0.9999999999999998,7,2,2,1,1,-2,-2,-0.9999999999999998,12,-2,8,19]

bi.Dgemm(blas.Trans, blas.NoTrans, is-1, nb, mb, -1, a[lda*(is-1):], lda, c[is-1+ldc*(js-1)+3:], ldc, 1, c[ldc*(js-1):], ldc)

and getting

[]float64 len: 16, cap: 16, [-8,-11,-6,-0.9999999999999998,5,1,1,1,0,-3,-2,-0.9999999999999998,12,-2,8,19]