Stream-lined determinant
Computing the determinant of a general square matrix
slightly modified Gauss-Jordan eliminations using only sequential access to Matrix elements
From the LinAlg package; some code (and comments) have been snipped for clarity
Store the matrix being swept, along with what has been pivoted
class MatrixPivoting : public Matrix
{
typedef REAL * INDEX;
INDEX * const row_index;
// Information about the pivot that was just picked up
double pivot_value;
INDEX pivot_row;
INDEX pivot_col;
int pivot_odd;
void pick_up_pivot(void);
public:
MatrixPivoting(const Matrix& m);
~MatrixPivoting(void);
double pivoting_and_elimination(void);
};
Pick up a pivot, an element with the largest abs value from yet not-pivoted rows and cols
void MatrixPivoting::pick_up_pivot(void)
{
register REAL max_elem = -1; // Abs value of the largest element
INDEX * prpp = 0; // Position of the pivot in row_index
INDEX * pcpp = 0; // Position of the pivot in index
int col_odd = 0; // Parity of the current column
for(register INDEX * cpp = index; cpp < index + ncols; cpp++)
{
register const REAL * cp = *cpp;
if( cp == 0 ) // skip over already pivoted col
continue;
int row_odd = 0;
for(register INDEX * rip = row_index; rip < row_index + nrows; rip++,cp++)
if( *rip )
{
const REAL v = *cp;
if( ::abs(v) > max_elem )
{
max_elem = ::abs(v);
pivot_value = v;
prpp = rip;
pcpp = cpp;
pivot_odd = row_odd ^ col_odd;
}
row_odd ^= 1; // Toggle parity for the next row
}
col_odd ^= 1;
}
assure( max_elem >= 0 && prpp != 0 && pcpp != 0,
"All the rows and columns have been already pivoted and eliminated");
pivot_row = *prpp, *prpp = 0;
pivot_col = *pcpp, *pcpp = 0;
}
Gaussian elimination
double MatrixPivoting::pivoting_and_elimination(void)
{
pick_up_pivot();
if( pivot_value == 0 )
return 0;
assert( pivot_row != 0 && pivot_col != 0 );
register REAL * pcp; // Pivot column pointer
register const INDEX * rip; // Current ptr in row_index
// Divide the pivoted column by pivot
for(pcp=pivot_col,rip=row_index; rip
Next | Table of contents