Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages

matrix_inv.cpp

Go to the documentation of this file.
00001 // This may look like C code, but it is really -*- C++ -*-
00002 /*
00003  ************************************************************************
00004  *
00005  *                      Linear Algebra Package
00006  *
00007  *                      Find the matrix inverse
00008  *              for matrices of general and special forms
00009  *
00010  * $Id: matrix_inv.cpp,v 1.1 2004/05/21 21:02:52 maxx Exp $
00011  *
00012  ************************************************************************
00013  */
00014 
00015 #include "LinAlg.h"
00016 #include <math.h>
00017 #if defined(WIN32)
00018 #include <malloc.h>
00019 #else
00020 #include <alloca.h>
00021 #endif
00022 
00023 namespace linalg 
00024 {
00025     using namespace linalg;
00026     
00027 /*
00028  *------------------------------------------------------------------------
00029  *
00030  *              The most general (Gauss-Jordan) matrix inverse
00031  *
00032  * This method works for any matrix (which of course must be square and
00033  * non-singular). Use this method only if none of specialized algorithms
00034  * (for symmetric, tridiagonal, etc) matrices is applicable/available.
00035  * Also, the matrix to invert has to be _well_ conditioned:
00036  * Gauss-Jordan eliminations (even with pivoting) perform poorly for
00037  * near-singular matrices (e.g., Hilbert matrices).
00038  *
00039  * The method inverts matrix inplace and returns the determinant if
00040  * determ_ptr was specified as not nil. determinant will be exactly zero
00041  * if the matrix turns out to be (numerically) singular. If determ_ptr is
00042  * nil and matrix happens to be singular, throw up.
00043  *
00044  * The algorithm perform inplace Gauss-Jordan eliminations with
00045  * full pivoting. It was adapted from my Algol-68 "translation" (ca 1986)
00046  * of a FORTRAN code described in
00047  * Johnson, K. Jeffrey, "Numerical methods in chemistry", New York,
00048  * N.Y.: Dekker, c1980, 503 pp, p.221
00049  *
00050  * Note, since it's much more efficient to perform operations on matrix
00051  * columns rather than matrix rows (due to the layout of elements in the
00052  * matrix), the present method implements a "transposed" algorithm.
00053  * Also, this algorithm was modified to avoid random access to matrix
00054  * elements.
00055  */
00056 
00057 double Matrix::dummy_determinant_ref = 0;
00058 
00059 static bool is_dummy_det_ref(double &determ_ref)
00060 { return &determ_ref == &Matrix::dummy_determinant_ref; }
00061 
00062 Matrix& Matrix::invert(double &determ_ref)
00063 {
00064   is_valid();
00065   if( nrows != ncols )
00066     info(),
00067     _error("Matrix to invert must be square");
00068 
00069   double determinant = 1;
00070   const double singularity_tolerance = 1e-35;
00071 
00072                                 // Locations of pivots (indices start with 0)
00073 #ifdef WIN32
00074   struct Pivot { int row, col; } * const pivots = 
00075                         (Pivot*)malloc(ncols*sizeof(Pivot));
00076   bool * const was_pivoted = (bool*)malloc(nrows*sizeof(bool)); 
00077 #else
00078   struct Pivot { int row, col; } * const pivots = 
00079                         (Pivot*)alloca(ncols*sizeof(Pivot));
00080   bool * const was_pivoted = (bool*)alloca(nrows*sizeof(bool)); 
00081 #endif
00082   memset(was_pivoted,false,nrows*sizeof(bool));
00083 
00084   for(register Pivot * pivotp = &pivots[0]; pivotp < &pivots[ncols]; pivotp++)
00085   {
00086     const REAL * old_ppos = 0;          // Location of a pivot to be
00087     int prow = 0, pcol = 0;             // Pivot's row and column indices
00088     {                                   // Look through all non-pivoted cols
00089       REAL max_value = 0;               // (and rows) for a pivot (max elem)
00090       register const REAL * cp = elements;  // column pointer
00091       for(register int j=0; j<ncols; j++)
00092         if( !was_pivoted[j] )
00093         {                               // the following loop would increment
00094           REAL curr_value = 0;          // cp by nrows
00095           for(register int k=0; k<nrows; k++,cp++)
00096             if( !was_pivoted[k] && (curr_value = fabs(*cp)) > max_value )
00097               max_value = curr_value, prow = k, pcol = j, old_ppos = cp;
00098         }
00099         else
00100           cp += nrows;                  // and this branch would too
00101       if( max_value < singularity_tolerance )
00102         if( !is_dummy_det_ref(determ_ref) )
00103         {
00104           determ_ref = 0;
00105           return *this;
00106         }
00107         else
00108           _error("Matrix turns out to be singular: can't invert");
00109       pivotp->row = prow;
00110       pivotp->col = pcol;
00111     }
00112 
00113     REAL * const old_colp = const_cast<REAL*>(old_ppos) - prow;
00114     REAL * const new_colp = ( prow == pcol ? old_colp : elements + prow*nrows );
00115     if( prow != pcol )                  // Swap prow-th and pcol-th columns to
00116     {                                   // bring the pivot to the diagonal
00117       register REAL * cr = new_colp;
00118       register REAL * cc = old_colp;
00119       for(register int k=0; k<nrows; k++)
00120       {
00121         REAL temp = *cr; *cr++ = *cc; *cc++ = temp;
00122       }
00123     }
00124     was_pivoted[prow] = true;
00125 
00126     {                                   // Normalize the pivot column and
00127       register REAL * pivot_cp = new_colp;
00128       const double pivot_val = pivot_cp[prow];  // pivot is at the diagonal
00129       determinant *= pivot_val;         // correct the determinant
00130       pivot_cp[prow] = 1;
00131       for(register int k=0; k<nrows; k++)
00132         *pivot_cp++ /= pivot_val;
00133     }
00134 
00135     {                                   // Perform eliminations
00136       register REAL * pivot_rp = elements + prow;       // pivot row
00137       register REAL * ep = elements;            // sweeps all matrix' columns
00138       for(register int k=0; k<ncols; k++, pivot_rp += nrows)
00139         if( k != prow )
00140         {
00141           const double temp = *pivot_rp;
00142           *pivot_rp = 0;
00143           register const REAL * pivot_cp = new_colp;    // pivot column
00144           for(register int l=0; l<nrows; l++)
00145             *ep++ -= temp * *pivot_cp++;
00146         }
00147         else
00148           ep += nrows;
00149     }
00150   }
00151 
00152   int no_swaps = 0;             // Swap exchanged *rows* back in place
00153   for(register const Pivot * pvp = &pivots[ncols-1];
00154       pvp >= &pivots[0]; pvp--)
00155     if( pvp->row != pvp->col )
00156     {
00157       no_swaps++;
00158       register REAL * rp = elements + pvp->row;
00159       register REAL * cp = elements + pvp->col;
00160       for(register int k=0; k<ncols; k++, rp += nrows, cp += nrows)
00161       {
00162         const REAL temp = *rp; *rp = *cp; *cp = temp;
00163       }
00164     }
00165 
00166   if( !is_dummy_det_ref(determ_ref) )
00167     determ_ref = ( no_swaps & 1 ? -determinant : determinant );
00168 
00169   return *this;
00170 }
00171  
00172 }

Generated on Wed Dec 15 21:20:29 2004 for vuVolume by  doxygen 1.3.9.1