2 fail with tip chol.cc

[John W. Eaton]

On 21-Jan-2011, Doug Stewart wrote:

|  src/DLD-FUNCTIONS/bsxfun.cc ............................ PASS   72/72  
|   src/DLD-FUNCTIONS/cellfun.cc ........................... PASS   77/77  
|   src/DLD-FUNCTIONS/chol.cc .............................. PASS   24/26  
| FAIL 2
|   src/DLD-FUNCTIONS/conv2.cc ............................. PASS    4/4   
|   src/D

I traced this back to products like X'*X, which are giving me
uninitialized data errors when X is complex.  I can generate the
errors with something as simple as

  X = [ 0.5585528 + 0.0000000i  -0.1662088 - 0.0315341i;
       -0.1662088 + 0.0315341i   0.6760061 + 0.0000000i];

  X'*X

  Y = single (X);

  Y'*Y

It does not seem to be a problem for real-valued matrices.

The following patch avoids the problem for me, but looking at CHERK in
the reference BLAS, I doesn't appear that it should be necessary to
initialize the output matrix C if BETA is zero, so this looks like it
might be a bug in ATLAS.

I'm also attaching a Fortran program that will show the problem on my
system if I link with ATLAS (I have the libblas3gf package on a Debian
amd64 system).  When I use the reference BLAS version of CHERK, I
don't see a problem.

Jaroslav, can you duplicate this problem?  Would you prefer a
different fix?  For now, I'm checking in this patch.

jwe

diff --git a/liboctave/CMatrix.cc b/liboctave/CMatrix.cc
--- a/liboctave/CMatrix.cc
+++ b/liboctave/CMatrix.cc
@@ -3754,7 +3754,7 @@
         {
           octave_idx_type lda = a.rows ();
 
-          retval = ComplexMatrix (a_nr, b_nc);
+          retval = ComplexMatrix (a_nr, b_nc, 0.0);
           Complex *c = retval.fortran_vec ();
 
           const char ctra = get_blas_trans_arg (tra, cja);
diff --git a/liboctave/fCMatrix.cc b/liboctave/fCMatrix.cc
--- a/liboctave/fCMatrix.cc
+++ b/liboctave/fCMatrix.cc
@@ -3750,7 +3750,7 @@
         {
           octave_idx_type lda = a.rows ();
 
-          retval = FloatComplexMatrix (a_nr, b_nc);
+          retval = FloatComplexMatrix (a_nr, b_nc, 0.0);
           FloatComplex *c = retval.fortran_vec ();
 
           const char ctra = get_blas_trans_arg (tra, cja);

      program main

      real x, y, inf
      complex a(2,2), c(2,2)

      a(1,1) = cmplx (0.5585528, 0.0000000)
      a(1,2) = cmplx (-0.1662088, -0.0315341)
      a(2,1) = cmplx (-0.1662088, 0.0315341)
      a(2,2) = cmplx (0.6760061, 0.0000000)

      x = 1.0
      y = 0.0

      inf = x / y

      c(1,1) = inf
      c(1,2) = inf
      c(2,1) = inf
      c(2,1) = inf

      call cherk ('U', 'C', 2, 2, 1.0, a, 2, 0.0, c, 2)

      print *, a
      print *, c

      end