Rounding problems BLAS SGEMM??

[Jokin]

Hello all,

Yesterday I had problems using the BLAS SGEMM routine. Today my problems arise when I compare the results I get.

For A:
3 3
-0.584736229226734 -0.491227207020988 0.645584520657102
0.703072824765721 -0.703876587238667 0.101224270870166
0.404687713311944 0.513082431018259 0.756950641504942

B:
3 3
-0.584734597596954 0.703073681013485 0.404688583282598
-0.491232818974689 -0.703876705893774 0.513076895272347
0.645581728315873 0.101217498320586 0.756953928649831

I get, using BLAS SGEMM:
3 3
1.000000000000000 -0.000004813849046 0.000004339528459
0.000004803181582 1.000000000000000 0.000004858287411
-0.000004346621154 -0.000004867824373 1.000000000000000

The result should be (using Matlab):

0.999999999979023 -4.81464564645873e-006 4.33273549999003e-006
4.81462467211147e-006 0.999999999976693 4.84087351078355e-006
-4.33275880712403e-006 -4.8408526502064e-006 0.999999999978897

Any idea?

Thanks in advance,
J. Zurutuza

[buttari]

Hello all,

Yesterday I had problems using the BLAS SGEMM routine. Today my problems arise when I compare the results I get.

For A:
3 3
-0.584736229226734 -0.491227207020988 0.645584520657102
0.703072824765721 -0.703876587238667 0.101224270870166
0.404687713311944 0.513082431018259 0.756950641504942

B:
3 3
-0.584734597596954 0.703073681013485 0.404688583282598
-0.491232818974689 -0.703876705893774 0.513076895272347
0.645581728315873 0.101217498320586 0.756953928649831

I get, using BLAS SGEMM:
3 3
1.000000000000000 -0.000004813849046 0.000004339528459
0.000004803181582 1.000000000000000 0.000004858287411
-0.000004346621154 -0.000004867824373 1.000000000000000

The result should be (using Matlab):

0.999999999979023 -4.81464564645873e-006 4.33273549999003e-006
4.81462467211147e-006 0.999999999976693 4.84087351078355e-006
-4.33275880712403e-006 -4.8408526502064e-006 0.999999999978897

Any idea?

Thanks in advance,
J. Zurutuza

Zurutuza,
can you show us your code? there's no way we can help you with this little information.

Alfredo

[Jokin]

I attach the code. The program reads two data sets froma A.txt and B.txt and produces athird file (CAB.txt) with the solution.

Thankyou,
J. Zurutuza

c Programa en fortran77 que resuelve un producto de
c dos matrices, usando BLAS
program mul_mat
implicit none
c
c Declaracion de variables
real ALPHA, BETA
integer I,J
integer K0, LDA0,LDB0,LDC0, M0, N0, MaxTam
parameter (MaxTam=1500)

real A0(MaxTam,MaxTam), B0(MaxTam,MaxTam), C0(MaxTam,MaxTam)

c Reading A.txt.......
open (10, File='A.txt')
read (10,*) LDA0, K0
write (*,*) 'Filas y Columnas de A: ',LDA0,'*', K0
READ (10,*) ((A0(I,J),J=1,K0),I=1,LDA0)
c write (*,*) ((A0(I,J), J=1,K0), I=1,LDA0)
write (*,*) 'Matriz A.txt leida. Leyendo B.txt'
close (10)
c Leyendo los elementos de B:
open (10, File='B.txt')
read (10,*) LDB0, N0
write (*,*) 'Filas y Columnas de B: ',LDB0,'*', N0
READ (10,*) ((B0(I,J),J=1,N0),I=1,LDB0)
c write (*,*) ((B0(I,J), J=1,N0), I=1,LDB0)
write (*,*) 'Matriz B.txt leida. Generando C.txt'
write (*,*) 'Dimensiones de C: ', LDA0,'*', N0
close (10)

c Actualiza LDC0:
M0 = LDA0
LDC0 = LDA0
ALPHA = 1.0
BETA = 0.0

c Mi llamada, a SGEMM
c OJO con MaxDim
call SGEMM('No transpuesta de A','No transpuesta de B',
$ M0, N0, K0, ALPHA, A0, MaxTam, B0, MaxTam, BETA,C0, MaxTam)

c Nuevos formatos de salida:
print 1000
c Bucle para leer todas las filas de A:
do 100 I=1, LDA0
print 1010, (A0(I,J), J=1,K0)
100 continue
print 1020
c Para B:
do 110 I=1, LDB0
print 1010, (B0(I,J), J=1,N0)
110 continue
print 1040
c Para C:
do 120 I=1, LDC0
print 1010, (C0(I,J), J=1,N0)
120 continue

open (10, File='AB.txt')
write (10,*) M0, N0
do 130 I=1, LDC0
write (10,1010), (C0(I,J), J=1,N0)
130 continue

c write (10,1030) ((C0(I,J), J=1,N0), I=1,M0)
close (10)

1000 format (1X, 'A:')
1010 format (1500(F30.15))
1020 format (/1X, 'B:')
1030 format (1500(F25.15))
1040 format (/1X, 'AB:')

end

[Julien Langou]

Hello,

I did not look at your code but your output is correct. Keep in mind that
SGEMM performs operation in single precision arithmetic. So the result
you will obtain from SGEMM will be accurate up to 1e-7. As far as I can
tell there is nothing to complain about your result.

Julien.

[Jokin]

Yes. The code is ok, but the routine to be called was not. I have tried with the DGESVD instead and it works perfectly.

Sorry for these simple-novicer questions, and thankyou very much for your patience,
J. Zurutuza