// The need for `USE_FC_LEN_T` and `FCONE` is due to a Fortran character string
// to C incompatibility. See: Writing R Extentions: 6.6.1 Fortran character strings
#define USE_FC_LEN_T
#include <R.h>
#include <Rinternals.h>
#include <R_ext/BLAS.h>
#include <R_ext/Lapack.h>
#ifndef FCONE
#define FCONE
#endif

// Singular Value Decomposition
// @note assumes passed arguments to be "clean", "properly checked"
int svd(
    const int p, const int q,
    double* A, const int ldA,
    double* d,
    double* U, const int ldU,
    double* Vt, const int ldVt,
    double* work_mem, int work_size, int* i_work    // if unknown set to 0, -1, 0
) {
    // LAPACK information return variable (in case of an error, LAPACK sets this
    // to a non-zero value)
    int info = 0;

    // check if work memory is supplied, if _not_ query for appropriate size and
    // allocate the memory
    if (!work_mem || work_size < 1 || !i_work) {
        // create (declare) work memory
        i_work = (int*)R_alloc(8 * (p < q ? p : q), sizeof(int));
        double tmp;         // variable to store work memory size query result
        work_mem = &tmp;    // point to `tmp` as it will contain the result of
                            // the work memory size query

        // request appropriate work memory size
        F77_CALL(dgesdd)(
            "S", &p, &q, A, &ldA, d, U, &ldU, Vt, &ldVt,
            work_mem, &work_size, i_work, &info
        FCONE);

        // allocate work memory
        work_size = (int)tmp;   // "read" work memory size query result
        work_mem = (double*)R_alloc(work_size, sizeof(double));

        // in case of an error, return error code stored in `info`
        if (info) { return info; }
    }

    // actual SVD computation
    F77_CALL(dgesdd)(
        "S", &p, &q, A, &ldA, d, U, &ldU, Vt, &ldVt,
        work_mem, &work_size, i_work, &info
    FCONE);

    return info;
}

// Singular Valued Decomposition (R binding)
// @note educational purpose, same as `La.svd`
extern SEXP R_svd(SEXP A) {
    // Check if we got a real-valued matrix
    if (!isMatrix(A) || !isReal(A)) {
        error("Require a real-valued matrix");
    }

    // extract matrix dimensions
    int* dims = INTEGER(coerceVector(getAttrib(A, R_DimSymbol), INTSXP));
    int p = dims[0];        // = nrow(A)
    int q = dims[1];        // = ncol(A)
    int m = p < q ? p : q;  // = min(p, q) = min(dim(A))

    // Check if dimensions are not degenerate
    if (p < 1 || q < 1) {
        error("Expected positive matrix dimensions");
    }

    // create R objects to store the SVD result `A = U diag(d) V^T`
    SEXP U = PROTECT(allocMatrix(REALSXP, p, m));
    SEXP d = PROTECT(allocVector(REALSXP, m));
    SEXP Vt = PROTECT(allocMatrix(REALSXP, m, q));

    // Call C SVD routine
    int info = svd(
        p, q,           // nrow(A), ncol(A)
        REAL(A), p,     // mem. addr. of A, ldA (leading dimension of A)
        REAL(d),        // mem. addr. of d
        REAL(U), p,     // mem. addr. of U, ldU (leading dimension of U)
        REAL(Vt), m,    // mem. addr. of V^T, ldVt (leading dimension of V^T)
        0, -1, 0        // work mem. pointer, work mem. size and int work mem
    );                  //  set to `0, -1, 0` to indicate "unknown"

    // Check LAPACK info
    if (info) {
        error("error code %d from LAPACK routine 'dgesdd'", info);
    }

    // Create R list containint SVD components
    SEXP result = PROTECT(allocVector(VECSXP, 3));
    SEXP names = PROTECT(allocVector(STRSXP, 3));
    SET_VECTOR_ELT(result, 0, d);
    SET_VECTOR_ELT(result, 1, U);
    SET_VECTOR_ELT(result, 2, Vt);
    SET_STRING_ELT(names, 0, mkChar("d"));
    SET_STRING_ELT(names, 1, mkChar("u"));
    SET_STRING_ELT(names, 2, mkChar("vt"));
    setAttrib(result, R_NamesSymbol, names);

    // Release created objects to the garbage collector
    UNPROTECT(5);

    return result;
}