liboctave matrix inverse function inaccurate

[John W. Eaton]

On 16-Oct-2007, Alex Roper wrote:

| I'm attempting to use liboctave to do matrix math for me in my CS171 computer 
| graphics class here at Caltech. Everything works smoothly except for matrix 
| inversion, which returns a matrix of correct dimension full of 0's. I believe 
| it to be tied to some sort of precision issue, but other matrix operations 
| seem to work perfectly.
| 
| Inverting the identity matrix succeeds, as does (occasionally) inverting 
| matrices that have close to integral results.
| 
| I understand this may not be an issue with Octave itself, given the same 
| matrix inversion succeeds when using Octave interactively, but since other 
| matrix operations succeed I'm not sure. I'd appreciate any advice people can 
| give.
| 
| Repeat-By:
| ---------
| 
| I am attaching full source code and a makefile (a total of about 155 lines of 
| code, but most of that is testing matrix routines that succeed). The failure 
| occurs at lines 59-69 of matrixtest.cpp (the end). oMatrix is a very thin 
| class inherited from Matrix

A quick fix that does not require patching Octave is to use

   octave_idx_type info;
   double rcond;
   oMatrix b = A.inverse(info, rcond);

in your code instead of

   oMatrix b = A.inverse();

I've checked in the following patch to my Octave sources so that the
calculation will work correctly even if you do not request INFO and
RCOND from the inverse method.

Thanks,

jwe


liboctave/ChangeLog:

2007-10-16  John W. Eaton  <address@hidden>

        * dMatrix.cc (Matrix::inverse): Only check rcond == 0 if the
        matrix is hermitian or calc_cond is true.
        * CMatrix.cc (ComplexMatrix::inverse): Likewise.


Index: liboctave/CMatrix.cc
===================================================================
RCS file: /cvs/octave/liboctave/CMatrix.cc,v
retrieving revision 1.136
diff -u -r1.136 CMatrix.cc
--- liboctave/CMatrix.cc        12 Oct 2007 21:27:13 -0000      1.136
+++ liboctave/CMatrix.cc        16 Oct 2007 17:46:38 -0000
@@ -1202,7 +1202,7 @@
       if (!mattype.is_hermitian ())
        ret = finverse(mattype, info, rcond, force, calc_cond);
 
-      if (rcond == 0.)
+      if ((mattype.is_hermitian () || calc_cond) && rcond == 0.)
        ret = ComplexMatrix (rows (), columns (), Complex (octave_Inf, 0.));
     }
 
Index: liboctave/dMatrix.cc
===================================================================
RCS file: /cvs/octave/liboctave/dMatrix.cc,v
retrieving revision 1.139
diff -u -r1.139 dMatrix.cc
--- liboctave/dMatrix.cc        12 Oct 2007 21:27:15 -0000      1.139
+++ liboctave/dMatrix.cc        16 Oct 2007 17:46:39 -0000
@@ -871,7 +871,7 @@
       if (!mattype.is_hermitian ())
        ret = finverse(mattype, info, rcond, force, calc_cond);
 
-      if (rcond == 0.)
+      if ((mattype.is_hermitian () || calc_cond) && rcond == 0.)
        ret = Matrix (rows (), columns (), octave_Inf);
     }