From bug-octave-request at che dot utexas dot edu  Wed Jun 29 09:09:26 1994
Subject: Re: octave for linux 
From: John Eaton <jwe at che dot utexas dot edu>
To: Janne dot Sinkkonen at helsinki dot fi
cc: matthew at redbird dot umsl dot edu (Matthew Feldt), bug-octave@che.utexas.edu
Date: Wed, 29 Jun 94 09:09:05 EDT

janne at avocado dot pc dot helsinki dot fi (Janne Sinkkonen) wrote:

:  > am confused about some results I have seen.  If you could offer insight
:  > to the following results or need more information please contact me at
:  > the e-mail address below.
: 
:  > octave:1> a = [ 1 2 3; 4 5 6; 7 8 9]
: 
: What's the problem? Your matrix is singular:
: 
: octave:262> cond(a)
: ans =  6.0262e+16
:
: octave:263> svd(a)
: ans =
: 
:    1.6848e+01
:    1.0684e+00
:    2.7958e-16
: 
: So, don't expect very accurate results from inverse().

True, but you should get something like

  octave:1> det ([1,2,3;4,5,6;7,8,9])
  warning: det: matrix singular to machine precision, rcond = 2.05597e-18
  octave:2> inv ([1,2,3;4,5,6;7,8,9])
  warning: inverse: matrix singular to machine precision, rcond = 2.05597e-18

instead.  This is what happens for me on a DECstation, SPARCstation,
and an RS/6000.

(Though I suppose a result should be returned, or the warnings should
become errors instead.)

I suspect that the problem is due to the following code in Matrix.cc
(which is repeated several times in various functions):

  if (rcond + 1.0 == rcond)
    ... matrix is singular to machine precision ...

and that gcc is placing the temporaries in the x87 coprocessor's
extended precision registers, which ends up making this statement
false, instead of true as it should be.

Declaring rcond as `volatile' will probably fix the problem.  Can
someone try this and tell me if it does?

: I guess determinant is most probably NOT calculated by the familiar
: sum of products of permutations rule, but instead by doing something
: like prod(eig(a)) but more clever.

Since the source for Octave is available, there's no need to guess.

The computations for inverse and determinant are computed using the
Linpack routines dgeco and dgedi.

Thanks,

jwe