74 lines
2.1 KiB
C
74 lines
2.1 KiB
C
#include "det.h"
|
|
|
|
|
|
// For reference see: `det_ge_real` in "R-4.2.1/src/modules/lapack/Lapack.c"
|
|
// of the R source code.
|
|
double det(
|
|
/* dim */ const int dimA,
|
|
/* matrix */ const double* A, const int ldA,
|
|
double* work_mem, int* info
|
|
) {
|
|
// if working memory size query, return immediately
|
|
if (work_mem == NULL) {
|
|
*info = dimA * (dimA + 1);
|
|
return 0.0;
|
|
}
|
|
|
|
// determinant of "zero size" matrix is 1 (by definition)
|
|
if (dimA == 0) {
|
|
return 1.0;
|
|
}
|
|
|
|
// copy `A` to (continuous) `work_mem` cause `dgetrf` works "in place"
|
|
for (int i = 0; i < dimA; ++i) {
|
|
memcpy(work_mem + i * dimA, A + i * ldA, dimA * sizeof(double));
|
|
}
|
|
|
|
// L U factorization of `A`
|
|
int error = 0;
|
|
int* ipvt = (int*)(work_mem + dimA * dimA);
|
|
F77_CALL(dgetrf)(&dimA, &dimA, work_mem, &dimA, ipvt, &error);
|
|
|
|
// check if an error occured which is the case iff `dgetrf` gives a negative error
|
|
*info |= (error < 0) * error;
|
|
// in both cases we return zero (ether the determinant is zero or error occured)
|
|
if (error) {
|
|
return 0.0;
|
|
}
|
|
|
|
// res <- det(A) = sign(P) * prod(diag(L)) where P is the pivoting permutation
|
|
double res = 1.0;
|
|
for (int i = 0; i < dimA; ++i) {
|
|
res *= (ipvt[i] != (i + 1) ? -1.0 : 1.0) * work_mem[i * (dimA + 1)];
|
|
}
|
|
|
|
return res;
|
|
}
|
|
|
|
/**
|
|
* R bindong to `det`
|
|
*/
|
|
extern SEXP R_det(SEXP A) {
|
|
// check if A is a real valued square matrix
|
|
if (!Rf_isReal(A) || !Rf_isMatrix(A) || Rf_nrows(A) != Rf_ncols(A)) {
|
|
Rf_error("`A` must be a real valued squae matrix");
|
|
}
|
|
|
|
// allocate working memory
|
|
int work_size;
|
|
(void)det(Rf_nrows(A), NULL, Rf_nrows(A), NULL, &work_size);
|
|
double* work_mem = (double*)R_alloc(work_size, sizeof(double));
|
|
|
|
// compute determinant (followed only by if-statement, no protection required)
|
|
int error = 0;
|
|
SEXP res = Rf_ScalarReal(
|
|
det(Rf_nrows(A), REAL(A), Rf_nrows(A), work_mem, &error)
|
|
);
|
|
|
|
if (error) {
|
|
Rf_error("Encountered error code %d in `det`", error);
|
|
}
|
|
|
|
return res;
|
|
}
|