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

determinant.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  *          Compute the determinant of a general square matrix
00008  *
00009  * Synopsis
00010  *      Matrix A;
00011  *      double A.determinant();
00012  * The matrix is assumed to be square. It is not altered.
00013  *
00014  * Method
00015  *      Gauss-Jordan transformations of the matrix with a slight
00016  *      modification to take advantage of the *column*-wise arrangement
00017  *      of Matrix elements. Thus we eliminate matrix's columns rather than
00018  *      rows in the Gauss-Jordan transformations. Note that determinant
00019  *      is invariant to matrix transpositions.
00020  *      The matrix is copied to a special object of type MatrixPivoting,
00021  *      where all Gauss-Jordan eliminations with full pivoting are to
00022  *      take place.
00023  *
00024  * $Id: determinant.cpp,v 1.1 2004/05/21 21:02:52 maxx Exp $
00025  *
00026  ************************************************************************
00027  */
00028 
00029 #include "LinAlg.h"
00030 #include <math.h>
00031 
00032 namespace linalg 
00033 {
00034     using namespace linalg;
00035     
00036 /*
00037  *------------------------------------------------------------------------
00038  *                      Class MatrixPivoting
00039  *
00040  * It is a descendant of a Matrix which keeps additional information
00041  * about what is being/has been pivoted 
00042  */
00043 
00044 class MatrixPivoting : public Matrix
00045 {
00046   typedef REAL * INDEX;                 // wanted to have typeof(index[0])
00047   INDEX * const row_index;              // row_index[i] = ptr to the i-th
00048                                         // matrix row, or 0 if the row
00049                                         // has been pivoted. Note,
00050   INDEX * const col_index;              // col_index[j] = ptr to the j-th
00051                                         // matrix col, or 0 if the col
00052                                         // has been pivoted.
00053 
00054                                 // Information about the pivot that was
00055                                 // just picked up
00056   double pivot_value;                   // Value of the pivoting element
00057   INDEX pivot_row;                      // pivot's location (ptrs)
00058   INDEX pivot_col;
00059   int pivot_odd;                        // parity of the pivot
00060                                         // (0 for even, 1 for odd)
00061 
00062   void pick_up_pivot(void);             // Pick up a pivot from
00063                                         // not-pivoted rows and cols
00064 
00065   MatrixPivoting(const MatrixPivoting&);    // Deliberately unimplemented:
00066   void operator = (const MatrixPivoting&);  // no copying/cloning allowed!
00067 
00068 public:
00069   MatrixPivoting(const Matrix& m);      // Construct an object 
00070                                         // for a given matrix
00071   ~MatrixPivoting(void);
00072 
00073   double pivoting_and_elimination(void);// Perform the pivoting, return
00074                                         // the pivot value times (-1)^(pi+pj)
00075                                         // (pi,pj - pivot el row & col)
00076 };
00077 
00078 
00079 /*
00080  *------------------------------------------------------------------------
00081  *              Constructing and decomissioning of MatrixPivoting
00082  */
00083 
00084 MatrixPivoting::MatrixPivoting(const Matrix& m)
00085   : Matrix(m), row_index(new INDEX[nrows]), 
00086     col_index(new INDEX[ncols]), pivot_value(0),
00087     pivot_row(0), pivot_col(0), pivot_odd(false)
00088 {
00089   assert( row_index != 0 && col_index != 0);
00090   register INDEX rp = elements;         // Fill in the row_index
00091   for(register INDEX * rip = row_index; rip<row_index+nrows;)
00092     *rip++ = rp++;
00093   register INDEX cp = elements;         // Fill in the col_index
00094   for(register INDEX * cip = col_index; cip<col_index+ncols; cp += nrows)
00095     *cip++ = cp;
00096 }
00097 
00098 MatrixPivoting::~MatrixPivoting(void)
00099 {
00100   is_valid();
00101   delete row_index;
00102   delete col_index;
00103 }
00104 
00105 /*
00106  *------------------------------------------------------------------------
00107  *                              Pivoting itself
00108  */
00109 
00110                         // Pick up a pivot, an element with the largest
00111                         // abs value from yet not-pivoted rows and cols
00112 void MatrixPivoting::pick_up_pivot(void)
00113 {
00114   register REAL max_elem = -1;          // Abs value of the largest element
00115   INDEX * prpp = 0;                     // Position of the pivot in row_index
00116   INDEX * pcpp = 0;                     // Position of the pivot in index
00117 
00118   int col_odd = 0;                      // Parity of the current column
00119   
00120   for(register INDEX * cpp = col_index; cpp < col_index + ncols; cpp++)
00121   {
00122     register const REAL * cp = *cpp;    // Column pointer for the curr col
00123     if( cp == 0 )                       // skip over already pivoted col
00124       continue;
00125     int row_odd = 0;                    // Parity of the current row
00126     for(register INDEX * rip = row_index; rip < row_index + nrows; rip++,cp++)
00127       if( *rip )
00128       {                                 // only if the row hasn't been pivoted
00129         const REAL v = *cp;
00130         if( fabs(v) > max_elem )
00131         {
00132           max_elem = fabs(v);           // Note the local max of col elements
00133           pivot_value = v;
00134           prpp = rip;
00135           pcpp = cpp;
00136           pivot_odd = row_odd ^ col_odd;
00137         }
00138         row_odd ^= 1;                   // Toggle parity for the next row
00139       }
00140     col_odd ^= 1;
00141   }
00142 
00143   assure( max_elem >= 0 && prpp != 0 && pcpp != 0,
00144          "All the rows and columns have been already pivoted and eliminated");
00145                                 // Note the position of the pivot and mark
00146                                 // the corresponding rows/cols as pivoted
00147   pivot_row = *prpp, *prpp = 0;
00148   pivot_col = *pcpp, *pcpp = 0;
00149 }
00150 
00151                         // Perform pivoting and gaussian elemination,
00152                         // return the pivot value times pivot_parity
00153                         // The procedure places zeros to the pivot_row of
00154                         // all not yet pivoted columns
00155                         // A[i,j] -= A[i,pivot_col]/pivot*A[pivot_row,j]
00156 double MatrixPivoting::pivoting_and_elimination(void)
00157 {
00158   pick_up_pivot();
00159   if( pivot_value == 0 )
00160     return 0;
00161 
00162   assert( pivot_row != 0 && pivot_col != 0 );
00163   
00164   register REAL * pcp;                          // Pivot column pointer
00165   register const INDEX * rip;                   // Current ptr in row_index
00166                                         // Divide the pivoted column by pivot
00167   for(pcp=pivot_col,rip=row_index; rip<row_index+nrows; pcp++,rip++)
00168     if( *rip )                          // Skip already pivoted rows
00169       *pcp /= pivot_value;
00170 
00171                                         // Eliminate all the elements from
00172                                         // the pivot_row in all not pivoted
00173                                         // columns
00174   const REAL * prp = pivot_row;         // Pivot row ptr
00175   for(register const INDEX * cpp = col_index; cpp < col_index + ncols; cpp++,prp+=nrows)
00176   {
00177     register REAL * cp = *cpp;          // A[*,j]
00178     if( cp == 0 )                       // skip over already pivoted col
00179       continue;
00180     const double fac = *prp;            // fac = A[pivot_row,j]
00181                                         // Do elimination stepping over pivoted rows
00182     for(pcp=pivot_col,rip=row_index; rip<row_index+nrows; pcp++,cp++,rip++)
00183       if( *rip )
00184         *cp -= *pcp * fac;
00185   }
00186         
00187   return pivot_odd ? -pivot_value : pivot_value;
00188 }
00189 
00190 
00191 /*
00192  *------------------------------------------------------------------------
00193  *                              Root module
00194  */
00195 
00196 double Matrix::determinant(void) const
00197 {
00198   is_valid();
00199 
00200   if( nrows != ncols )
00201     info(), _error("Can't obtain determinant of a non-square matrix");
00202 
00203   if( row_lwb != col_lwb )
00204     info(), _error("Row and col lower bounds are inconsistent");
00205 
00206   MatrixPivoting mp(*this);
00207 
00208   register double det = 1;
00209   register int k;
00210 
00211   for(k=0; k<ncols && det != 0; k++)
00212     det *= mp.pivoting_and_elimination();
00213 
00214   return det;
00215 }
00216 }

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