#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;
}