113 lines
3.8 KiB
C
113 lines
3.8 KiB
C
// 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;
|
|
}
|