Dsfact.hh File Reference

Go to the source code of this file.

Functions
template<class Matrix , unsigned int n, unsigned int idim>
bool	Dsfact (Matrix &rhs, typename Matrix::value_type &det)

Function Documentation

Dsfact()

template<class Matrix , unsigned int n, unsigned int idim>

bool Dsfact	(	Matrix &	rhs,
		typename Matrix::value_type &	det
	)

Dsfact. Compute determinant of a symmetric, positive definite matrix of dimension $idim$ and order $n$.

Author

T. Glebe

                                                         {
 
#ifdef XXX
  /* Function Body */
  if (idim < n || n <= 0) {
    return false;
  }
#endif
 
  typename Matrix::value_type* a = rhs.Array();
 
#ifdef XXX
  const typename Matrix::value_type* A = rhs.Array();
  typename Matrix::value_type array[Matrix::size];
  typename Matrix::value_type* a = array;
 
  // copy contents of matrix to working place
  for(unsigned int i=0; i<Matrix::size; ++i) {
    array[i] = A[i];
  }
#endif
 
  /* Local variables */
  static unsigned int i, j, l;
 
  /* Parameter adjustments */
  a -= idim + 1;
 
  /* sfactd.inc */
  det = 1.;
  for (j = 1; j <= n; ++j) {
    const unsigned int ji = j * idim;
    const unsigned int jj = j + ji;
 
    if (a[jj] <= 0.) {
      det = 0.;
      return false;
    }
 
    const unsigned int jp1 = j + 1;
    const unsigned int jpi = jp1 * idim;
 
    det *= a[jj];
    a[jj] = 1. / a[jj];
 
    for (l = jp1; l <= n; ++l) {
      a[j + l * idim] = a[jj] * a[l + ji];
 
      const unsigned int lj = l + jpi;
 
      for (i = 1; i <= j; ++i) {
    a[lj] -= a[l + i * idim] * a[i + jpi];
      } // for i
    } // for l
  } // for j
 
  return true;
} // end of Dsfact