33 行
512 B
C++
33 行
512 B
C++
#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);
|
|
}
|