Complex left-divide and SVD give ireproducible results

[John W. Eaton]

On 28-Jul-2008, Dirk Laurie wrote:

| --------
| Bug report for Octave 2.1.73 configured for i486-pc-linux-gnu    
| 
| Description:
| -----------
| 
| The following operations give irreproducible results:
|     z=A\b  when A is a nonsquare complex matrix
|     [U,S,V]=svd(A) when A is a complex matrix
| 
| By 'irreproducible' is meant that the results are not identical
| to the last bit when the operation is run more than once on
| identical input, although the results in both cases are numerically
| correct in the sense that b-A*z and A-U*S*V respectively are at
| roundoff error level.  The two bugs are reported as one because 
| it is thought probable that they share a common cause.
| 
| A possible cause might be that an uninitialized value is used in 
| a comparison of which the result would normally not critically 
| affect the output, e.g. whether columns in a full-rank orthogonal 
| factorization are to be interchanged. 
| 
| 
| Repeat-By:
| ---------
| 
| The left-divide operation is so capricious that it suffices to 
| generate a random test case.  The following code exposes the bug:
|  
| A=rand(200,40)+1i*rand(200,40); b=rand(200,1); x0=A\b;
| for iter=1:10, x=A\b; err(iter)=norm(x0-x,Inf); end
| err
| 
| The otput should be all zeros.  If I save the above code in an M-file,
| I get nonzero results most of the time, whether running the code in 
| an interactive shell or from the command line.
| 
| The SVD operation can also be exposed by a random test case, 
| albeit one in which the matrix is scaled down, as follows:
| 
| A=eps*(rand(40,8)+1i*rand(40,8)); [U0,S0,V0]=svd(A);
| for iter=1:10,
|     [U,S,V]=svd(A);
|     err(iter)=max([norm(U-U0,Inf),norm(S-S0,Inf),norm(V-V0,Inf)]);
|     end
| err
| 
| The same comments as before apply to the output.
| 
| The same errors are also present, in fact manifesting more often, on
| all more recent Octave releases that I have access to, including
| Octave 3.0.0 as supplied with Ubuntu Hardy.

Of the systems I have easy access to, I'm only able to reproduce this
on one with Octave 2.1.73.  The failure happens on a 32-bit AMD system
running Debian.  I'm unable to reproduce the problem with Octave 3.0.0
on the same system.  The difference appears to be that Octave 3.0 is
linked with

  liblapack.so.3gf => /usr/lib/liblapack.so.3gf (0xa633d000)
  libblas.so.3gf => /usr/lib/libblas.so.3gf (0xa627c000)

which are provided by

  ii  libblas3gf                 1.2-1.5
  ii  liblapack3gf               3.1.1-0.4

and Octave 2.1.73 is linked with

  liblapack.so.3 => /usr/lib/atlas/liblapack.so.3 (0xa6b2e000)
  libblas.so.3 => /usr/lib/atlas/libblas.so.3 (0xa67cf000)

which are provided by

  ii  blas                       1.1-14
  ii  lapack3                    3.0.20000531a-6.1

I'd guess that the problem is in the linear algebra package, not
Octave.

jwe