ABACUS/src/MATRIX/lndet_LU_dstry.cc

#include "ABACUS.h"
using namespace std;

complex<DP> ABACUS::lndet_LU_dstry (SQMat_DP& a)
{
  // Returns the ln of determinant of matrix a, through LU decomposition
  // Does not preserve a, unlike lndet_LU

  Vect_INT indx(a.size());

  DP d;

  complex<DP> logdet = 0.0;

  try {

    ABACUS::ludcmp (a, indx, d);

    logdet = log(complex<DP>(d));

    for (int j = 0; j < a.size(); j++) {
      logdet += log(complex<DP>(a[j][j]));
    }
  }

  catch (Divide_by_zero) {
    logdet = log(-1.0);  // reset to nan
  }

  return(logdet);
}