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<row_index+nrows; pcp++,rip++)
    if( *rip )	
      *pcp /= pivot_value;

  const REAL * prp = pivot_row;	
  for(register const INDEX * cpp = index; cpp < index + ncols; cpp++,prp+=nrows)
  {
    register REAL * cp = *cpp;
    if( cp == 0 )
      continue;
    const double fac = *prp; 		// fac = A[pivot_row,j]
    					// Do elimination stepping over pivoted rows
    for(pcp=pivot_col,rip=row_index; rip<row_index+nrows; pcp++,cp++,rip++)
      if( *rip )
        *cp -= *pcp * fac;
  }
        
  return pivot_odd ? -pivot_value : pivot_value;
}



Next | Table of contents