matrix multiplication with libatlas gives wrong result for big matrices

[Petr Tobiska]

Binary package hint: libatlas3gf-sse2

Complex matrix multiplication via zgemm gives a wrong result when using /usr/lib/sse2/atlas/liblapack.so.3gf. When the same binary uses /usr/lib/liblapack.so.3gf the result is ok.

sample code:
--------- atlas-test.cpp ---------------------------
/* compile:
   g++ -Wall -o atlas-test atlas-test.cpp -lblitz -llapack
   run:
   LD_LIBRARY_PATH=/usr/lib LD_PRELOAD=/usr/lib/liblapack.so.3gf \
or
   LD_LIBRARY_PATH=/usr/lib/sse2/atlas \
   LD_PRELOAD=/usr/lib/sse2/atlas/liblapack.so.3gf \
with
     ./atlas-test
*/

#include <iostream>
#include <blitz/array.h>

typedef std::complex<double> Complex;
typedef blitz::Array<Complex,2> cMatrix;
#define F77NAME(x) x ## _

/* SUBROUTINE ZGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB,
   BETA, C, LDC ) */
extern "C" void F77NAME(zgemm)
  (const char*,const char*,const int&,const int&,const int&,const Complex&,
   const Complex*,const int&,const Complex*,const int&,const Complex&,
   const Complex*,const int&);

int main() {
  blitz::firstIndex ii;
  blitz::secondIndex jj;

  const int M = 139, N = 128;
  cMatrix prod( M, N );
  prod = 0.1* ii + 0.4* jj;

  cMatrix C( M, M );
  char A_op = 'T', B_op = 'N';
  int M1 = M, N1 = N;

  F77NAME(zgemm)( &A_op, &B_op, M1, M1, N1, 1.0, prod.data(), N, prod.data(), N, 0.0, C.data(), M );

  std::cout << C( 0, 0 ) << " correct result = " << 0.4*0.4*N1*(N1-1)*(2*N1-1)/6 << std::endl;
  return 0;
}
----------------- end of atlas-test.cpp -------------

For M1 = 72 the program gives a correct result, for M1 > 72 not. Tested for other flavor of libatlas3gf (-base, -sse), similar behavior, but different results.

Ubuntu 8.04.2, package versions:
libatlas3gf-sse2_3.6.0-21.1ubuntu3_i386.deb Source: atlas
libblitz0ldbl_0.9-7_i386.deb Source: blitz++
g++-4.2_4.2.4-1ubuntu3_i386.deb Source: gcc-4.2
HP 6510b, cat /proc/cpuinfo
processor : 0
vendor_id : GenuineIntel
cpu family : 6
model : 23
model name : Intel(R) Core(TM)2 Duo CPU T9300 @ 2.50GHz
stepping : 6
cpu MHz : 800.000
cache size : 6144 KB
physical id : 0
siblings : 2
core id : 0
cpu cores : 2
fdiv_bug : no
hlt_bug : no
f00f_bug : no
coma_bug : no
fpu : yes
fpu_exception : yes
cpuid level : 10
wp : yes
flags : fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush dts acpi mmx fxsr sse sse2 ss ht tm pbe nx lm constant_tsc arch_perfmon pebs bts pni monitor ds_cpl vmx est tm2 ssse3 cx16 xtpr sse4_1 lahf_lm ida
bogomips : 4992.42
clflush size : 64

processor : 1 - the same as processor 0

[Dimitrios Symeonidis (azimout)]

On Karmic the results of the two experiments are the same:
LD_LIBRARY_PATH=/usr/lib;LD_PRELOAD=/usr/lib/liblapack.so.3gf;./atlas-test
gives me
(8968.96,26706.4) correct result = 110541
just like
LD_LIBRARY_PATH=/usr/lib/sse2/atlas;LD_PRELOAD=/usr/lib/sse2/atlas/liblapack.so.3gf;./atlas-test

Is this what you see, too?

[FP (fabrice-pardo)]

In response to https://bugs.launchpad.net/ubuntu/+source/atlas/+bug/376739/comments/2
I suspect /usr/lib/sse2/atlas/liblapack.so.3gf has nothing to do with that bug.
The faulty lib is /usr/lib/sse2/atlas/libblas.so.3gf
- If libblas3gf is not installed, then the command
    (export LD_LIBRARY_PATH=/usr/lib/sse2; ./atlas-test)
gives wrong result, because wrong /usr/lib/sse2/atlas/libblas.so.3gf is used.
- If libblas3gf is installed, then the commands
    (export LD_LIBRARY_PATH=/usr/lib; ./atlas-test)
    (export LD_PRELOAD=/usr/lib/libblas.so.3gf;./atlas-test)
give good result, because /usr/lib/libblas.so.3gf is used.

There is a regression
libatlas3gf-sse2_3.6.0-22ubuntu1_i386 -> good
libatlas3gf-sse2_3.6.0-22ubuntu2_i386 -> wrong
Cf. https://bugs.launchpad.net/debian/+source/atlas/+bug/406520/comments/3

Atlas library should be compiled using the processor which wil be used for.
Hence distributing such a precompiled package is not optimal.

Also, atlas is a very delicate piece of code. All atlas packages
should depend on package libatlas-test, and all tests should be passed
before installation can be finished.
Imagine the result of a an engineer computation of the structure of a bridge
or of a composite plane using a faulty atlas library!

A very small piece of code is gived at comment
https://bugs.launchpad.net/debian/+source/atlas/+bug/406520/comments/1
and attached here.

[Sylvestre Ledru (sylvestre)]

FYI, this bug is fixed with Debian Experimental ATLAS packages.

[Pauli Virtanen (pauli-virtanen)]

Maybe this is fixed in the 3.8.3-22ubuntu2 packages in Maverick?