ABACUS/src/MATRIX/lndet_LU_CX.cc

#include "ABACUS.h"
using namespace std;

complex<DP> ABACUS::lndet_LU_CX (SQMat_CX a)
{
  // Returns the ln of determinant of matrix a, through LU decomposition

  Vect_INT indx(a.size());

  SQMat_CX mat = a;

  DP d;

  complex<DP> logdet = 0.0;

  try {

    ABACUS::ludcmp_CX (mat, indx, d);

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

    for (int j = 0; j < mat.size(); j++) {
      logdet += log(mat[j][j]);
    }
  }

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

  return(logdet);
}