102 lines
3.6 KiB
C
102 lines
3.6 KiB
C
// Suppress stripping R API prefixes, all API functions have the form `Rf_<name>`
|
|
// and public variables are called `R_<name>`.
|
|
#define R_NO_REMAP
|
|
#include <R.h>
|
|
#include <Rinternals.h>
|
|
|
|
/**
|
|
* Matrix Times a Vectorized Kronecker product
|
|
*
|
|
* C = A vec(B_1 %x% ... %x% B_r)
|
|
*
|
|
* Note the reverse order of the `B_k` in the Kronecker product.
|
|
*
|
|
* @param A matrix of dimensions `n` by `pp * qq`
|
|
* @param Bs list of matrices such that the product there Kronecker product has
|
|
* dimensions `pp` by `qq`.
|
|
*/
|
|
extern SEXP mtvk(SEXP A, SEXP Bs) {
|
|
|
|
if ((TYPEOF(A) != REALSXP) || !Rf_isMatrix(A)) {
|
|
Rf_error("First argument must be a numeric matrices");
|
|
}
|
|
|
|
// extract dimensions from `A` and `Bs`
|
|
size_t nrow = Rf_nrows(A);
|
|
size_t ncol = Rf_ncols(A);
|
|
size_t r = Rf_length(Bs);
|
|
size_t pp = 1; // `nrow(Reduce("%x%", Bs))` see below
|
|
size_t qq = 1; // `ncol(Reduce("%x%", Bs))` see below
|
|
|
|
// retrieve dimensions and direct memory access for each component in
|
|
// reverse order
|
|
// (Allocated memory freed at function exit and errors)
|
|
double** bs = (double**)R_alloc(r, sizeof(double*));
|
|
// dimensions of reshaped `A` columns, we treat the matrix `A` as an
|
|
// `nrow` by `p[1]` by ... by `p[2 r]` tensor where the `p[1:r]` dimensions
|
|
// as the row and `p[(r + 1):(2 r)]` the column dimensions of the Kronecker
|
|
// product components.
|
|
size_t* p = (size_t*)R_alloc(2 * r + 1, sizeof(size_t));
|
|
// Multi-Index into the reshaped `A`:
|
|
// `A[i, j] == do.call(array(A, dim = c(nrow, p), c(i, J)))`
|
|
size_t* J = (size_t*)R_alloc(2 * r + 1, sizeof(size_t));
|
|
|
|
// reversed order of Kronecker components
|
|
for (size_t k = 0; k < r; ++k) {
|
|
SEXP Bk = VECTOR_ELT(Bs, k);
|
|
if ((TYPEOF(Bk) == REALSXP) && Rf_isMatrix(Bk)) {
|
|
bs[r - 1 - k] = REAL(Bk);
|
|
pp *= (p[ r - 1 - k] = Rf_nrows(Bk));
|
|
qq *= (p[2 * r - 1 - k] = Rf_ncols(Bk));
|
|
// set all zero
|
|
J[k] = J[k + r] = 0;
|
|
} else {
|
|
Rf_error("Second argument must be a list of numeric matrices");
|
|
}
|
|
}
|
|
|
|
// Additional last element such that the condition `p[2 * r] <= (++J[2 * r])`
|
|
// (evaluated exactly once) is false. Allows to skip the check `k < 2 * r`
|
|
// in the loop for updating the Multi-Index `J`.
|
|
p[2 * r] = 127; // biggest value of smallest signed type
|
|
J[2 * r] = 0;
|
|
|
|
// check if dimensions match
|
|
if ((ncol != pp * qq) || (r < 1)) {
|
|
Rf_error("Dimension Missmatch");
|
|
}
|
|
|
|
// Create new R result vector
|
|
SEXP C = PROTECT(Rf_allocVector(REALSXP, nrow));
|
|
double* c = REAL(C);
|
|
memset(c, 0, nrow * sizeof(double));
|
|
|
|
// main operation (single iteration over `A`)
|
|
double* a = REAL(A);
|
|
for (size_t j = 0; j < ncol; ++j) {
|
|
// Compute `prod_k (B_k)_{J_k, J_k+r}` (identical for each row)
|
|
double prod_BJ = 1.0;
|
|
for (size_t k = 0; k < r; ++k) {
|
|
prod_BJ *= bs[k][J[k] + p[k] * J[k + r]];
|
|
}
|
|
|
|
// Add `j`th summand `A_i,J prod_k (B_k)_{J_k, J_k+r}` to each component
|
|
for (size_t i = 0; i < nrow; ++i) {
|
|
c[i] += a[i + j * nrow] * prod_BJ;
|
|
}
|
|
|
|
// Compute matching Multi-Index `J` for reshaped `A` columns
|
|
// In `R` the relation between `j` and `J` is
|
|
// `A[i, j] == do.call(array(A, dim = c(nrow, p), c(i, J)))`
|
|
// Note: the check `k < 2 * r` is skipped as described above.
|
|
for (size_t k = 0; p[k] <= (++J[k]); ++k) {
|
|
J[k] = 0;
|
|
}
|
|
}
|
|
|
|
// Onle the result object `C` needs to be unprotected
|
|
UNPROTECT(1);
|
|
|
|
return C;
|
|
}
|