Mutils.c
#include "Mutils.h"
#include <R_ext/Lapack.h>
char norm_type(const char *typstr)
{
char typup;
if (strlen(typstr) != 1)
error(
_("argument type[1]='%s' must be a character string of string length 1"),
typstr);
typup = toupper(*typstr);
if (typup == '1') typup = 'O'; /* aliases */
if (typup == 'E') typup = 'F';
if (typup != 'M' && typup != 'O' && typup != 'I' && typup != 'F')
error(_("argument type[1]='%s' must be one of 'M','1','O','I','F' or 'E'"),
typstr);
return typup;
}
char rcond_type(const char *typstr)
{
char typup;
if (strlen(typstr) != 1)
error(_("argument type[1]='%s' must be a character string of string length 1"),
typstr);
typup = toupper(*typstr);
if (typup == '1') typup = 'O'; /* alias */
if (typup != 'O' && typup != 'I')
error(_("argument type[1]='%s' must be one of '1','O', or 'I'"),
typstr);
return typup;
}
double get_double_by_name(SEXP obj, char *nm)
{
SEXP nms = getAttrib(obj, R_NamesSymbol);
int i, len = length(obj);
if ((!isReal(obj)) || (length(obj) > 0 && nms == R_NilValue))
error(_("object must be a named, numeric vector"));
for (i = 0; i < len; i++) {
if (!strcmp(nm, CHAR(STRING_ELT(nms, i)))) {
return REAL(obj)[i];
}
}
return R_NaReal;
}
SEXP
set_double_by_name(SEXP obj, double val, char *nm)
{
SEXP nms = getAttrib(obj, R_NamesSymbol);
int i, len = length(obj);
if ((!isReal(obj)) || (length(obj) > 0 && nms == R_NilValue))
error("object must be a named, numeric vector");
for (i = 0; i < len; i++) {
if (!strcmp(nm, CHAR(STRING_ELT(nms, i)))) {
REAL(obj)[i] = val;
return obj;
}
}
{
SEXP nx = PROTECT(allocVector(REALSXP, len + 1)),
nnms = allocVector(STRSXP, len + 1);
setAttrib(nx, R_NamesSymbol, nnms);
for (i = 0; i < len; i++) {
REAL(nx)[i] = REAL(obj)[i];
SET_STRING_ELT(nnms, i, duplicate(STRING_ELT(nms, i)));
}
REAL(nx)[len] = val;
SET_STRING_ELT(nnms, len, mkChar(nm));
UNPROTECT(1);
return nx;
}
}
SEXP as_det_obj(double val, int log, int sign)
{
SEXP det = PROTECT(allocVector(VECSXP, 2)),
nms = allocVector(STRSXP, 2),
vv = ScalarReal(val);
setAttrib(det, R_NamesSymbol, nms);
SET_STRING_ELT(nms, 0, mkChar("modulus"));
SET_STRING_ELT(nms, 1, mkChar("sign"));
setAttrib(vv, install("logarithm"), ScalarLogical(log));
SET_VECTOR_ELT(det, 0, vv);
SET_VECTOR_ELT(det, 1, ScalarInteger(sign));
setAttrib(det, R_ClassSymbol, mkString("det"));
UNPROTECT(1);
return det;
}
SEXP get_factors(SEXP obj, char *nm)
{
SEXP fac = GET_SLOT(obj, Matrix_factorSym),
nms = getAttrib(fac, R_NamesSymbol);
int i, len = length(fac);
if ((!isNewList(fac)) || (length(fac) > 0 && nms == R_NilValue))
error("factors slot must be a named list");
for (i = 0; i < len; i++) {
if (!strcmp(nm, CHAR(STRING_ELT(nms, i)))) {
return VECTOR_ELT(fac, i);
}
}
return R_NilValue;
}
SEXP set_factors(SEXP obj, SEXP val, char *nm)
{
SEXP fac = GET_SLOT(obj, Matrix_factorSym),
nms = getAttrib(fac, R_NamesSymbol), nfac, nnms;
int i, len = length(fac);
if ((!isNewList(fac)) || (length(fac) > 0 && nms == R_NilValue))
error("factors slot must be a named list");
for (i = 0; i < len; i++) {
if (!strcmp(nm, CHAR(STRING_ELT(nms, i)))) {
SET_VECTOR_ELT(fac, i, duplicate(val));
return val;
}
}
nfac = PROTECT(allocVector(VECSXP, len + 1));
nnms = PROTECT(allocVector(STRSXP, len + 1));
setAttrib(nfac, R_NamesSymbol, nnms);
for (i = 0; i < len; i++) {
SET_VECTOR_ELT(nfac, i, VECTOR_ELT(fac, i));
SET_STRING_ELT(nnms, i, duplicate(STRING_ELT(nms, i)));
}
SET_VECTOR_ELT(nfac, len, duplicate(val));
SET_STRING_ELT(nnms, len, mkChar(nm));
SET_SLOT(obj, Matrix_factorSym, nfac);
UNPROTECT(2);
return val;
}
/*MM: this is useful for all the ..CMatrix classes
(and ..R by [0] <-> [1]): */
SEXP dgCMatrix_set_Dim(SEXP x, int nrow)
{
int *dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
dims[0] = nrow;
dims[1] = length(GET_SLOT(x, Matrix_pSym)) - 1;
return x;
}
/* Fill in the "trivial remainder" in n*m array ;
* typically the 'x' slot of a "dtrMatrix" :
* But should be usable for double/logical/int/complex : */
#define MAKE_TRIANGULAR_BODY(_TO_, _FROM_, _ZERO_, _ONE_) \
{ \
int i, j, *dims = INTEGER(GET_SLOT(_FROM_, Matrix_DimSym)); \
int n = dims[0], m = dims[1]; \
\
if (*uplo_P(_FROM_) == 'U') { \
for (j = 0; j < n; j++) \
for (i = j+1; i < m; i++) \
_TO_[i + j*m] = _ZERO_; \
} else { \
for (j = 1; j < n; j++) \
for (i = 0; i < j && i < m; i++) \
_TO_[i + j*m] = _ZERO_; \
} \
if (*diag_P(_FROM_) == 'U') { \
j = (n < m) ? n : m; \
for (i = 0; i < j; i++) \
_TO_[i * (m + 1)] = _ONE_; \
} \
}
void make_d_matrix_triangular(double *to, SEXP from)
MAKE_TRIANGULAR_BODY(to, from, 0., 1.)
void make_i_matrix_triangular(int *to, SEXP from)
MAKE_TRIANGULAR_BODY(to, from, 0, 1)
/* Should work for double/logical/int/complex : */
#define MAKE_SYMMETRIC_BODY(_TO_, _FROM_) \
{ \
int i, j, n = INTEGER(GET_SLOT(_FROM_, Matrix_DimSym))[0]; \
\
if (*uplo_P(_FROM_) == 'U') { \
for (j = 0; j < n; j++) \
for (i = j+1; i < n; i++) \
_TO_[i + j*n] = _TO_[j + i*n]; \
} else { \
for (j = 1; j < n; j++) \
for (i = 0; i < j && i < n; i++) \
_TO_[i + j*n] = _TO_[j + i*n]; \
} \
}
void make_d_matrix_symmetric(double *to, SEXP from)
MAKE_SYMMETRIC_BODY(to, from)
void make_i_matrix_symmetric(int *to, SEXP from)
MAKE_SYMMETRIC_BODY(to, from)
/**
* Create a named vector of type TYP
*
* @param TYP a vector SEXP type (e.g. REALSXP)
* @param names names of list elements with null string appended
*
* @return pointer to a named vector of type TYP
*/
SEXP
Matrix_make_named(int TYP, char **names)
{
SEXP ans, nms;
int i, n;
for (n = 0; strlen(names[n]) > 0; n++) {}
ans = PROTECT(allocVector(TYP, n));
nms = PROTECT(allocVector(STRSXP, n));
for (i = 0; i < n; i++) SET_STRING_ELT(nms, i, mkChar(names[i]));
setAttrib(ans, R_NamesSymbol, nms);
UNPROTECT(2);
return ans;
}
#define Matrix_Error_Bufsiz 4096
SEXP check_scalar_string(SEXP sP, char *vals, char *nm)
{
SEXP val = ScalarLogical(1);
char *buf;
/* only allocate when needed: in good case, none is needed */
#define SPRINTF buf = Alloca(Matrix_Error_Bufsiz, char); R_CheckStack(); sprintf
if (length(sP) != 1) {
SPRINTF(buf, _("'%s' slot must have length 1"), nm);
} else {
const char *str = CHAR(STRING_ELT(sP, 0));
if (strlen(str) != 1) {
SPRINTF(buf, _("'%s' must have string length 1"), nm);
} else {
int i, len, match;
for (i = 0, len = strlen(vals), match = 0; i < len; i++) {
if (str[0] == vals[i])
return R_NilValue;
}
SPRINTF(buf, _("'%s' must be in '%s'"), nm, vals);
}
}
/* 'error' returns : */
val = mkString(buf);
return val;
#undef SPRINTF
}
SEXP dense_nonpacked_validate(SEXP obj)
{
int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym));
if ((dims[0] * dims[1]) != length(GET_SLOT(obj, Matrix_xSym)))
return mkString(_("length of x slot != prod(Dim)"));
return ScalarLogical(1);
}
#define PACKED_TO_FULL(TYPE) \
TYPE *packed_to_full_ ## TYPE(TYPE *dest, const TYPE *src, \
int n, enum CBLAS_UPLO uplo) \
{ \
int i, j, pos = 0; \
\
AZERO(dest, n*n); \
for (j = 0; j < n; j++) { \
switch(uplo) { \
case UPP: \
for (i = 0; i <= j; i++) dest[i + j * n] = src[pos++]; \
break; \
case LOW: \
for (i = j; i < n; i++) dest[i + j * n] = src[pos++]; \
break; \
default: \
error(_("'uplo' must be UPP or LOW")); \
} \
} \
return dest; \
}
PACKED_TO_FULL(double)
PACKED_TO_FULL(int)
#define FULL_TO_PACKED(TYPE) \
TYPE *full_to_packed_ ## TYPE(TYPE *dest, const TYPE *src, int n, \
enum CBLAS_UPLO uplo, enum CBLAS_DIAG diag) \
{ \
int i, j, pos = 0; \
\
for (j = 0; j < n; j++) { \
switch(uplo) { \
case UPP: \
for (i = 0; i <= j; i++) \
dest[pos++] = (i == j && diag== UNT) ? 1 : src[i + j*n]; \
break; \
case LOW: \
for (i = j; i < n; i++) \
dest[pos++] = (i == j && diag== UNT) ? 1 : src[i + j*n]; \
break; \
default: \
error(_("'uplo' must be UPP or LOW")); \
} \
} \
return dest; \
}
FULL_TO_PACKED(double)
FULL_TO_PACKED(int)
/**
* Copy the diagonal elements of the packed denseMatrix x to dest
*
* @param dest vector of length ncol(x)
* @param x pointer to an object representing a packed array
* @param n number of columns in the matrix represented by x
*
* @return dest
*/
void d_packed_getDiag(double *dest, SEXP x, int n)
{
double *xx = REAL(GET_SLOT(x, Matrix_xSym));
#define END_packed_getDiag \
int j, pos = 0; \
\
if (*uplo_P(x) == 'U') { \
for(pos= 0, j=0; j < n; pos += 1+(++j)) dest[j] = xx[pos]; \
} else { \
for(pos= 0, j=0; j < n; pos += (n - j), j++) dest[j] = xx[pos]; \
} \
return
END_packed_getDiag;
}
void l_packed_getDiag(int *dest, SEXP x, int n)
{
int *xx = LOGICAL(GET_SLOT(x, Matrix_xSym));
END_packed_getDiag;
}
#undef END_packed_getDiag
void tr_d_packed_getDiag(double *dest, SEXP x)
{
int n = INTEGER(GET_SLOT(x, Matrix_DimSym))[0];
SEXP val = PROTECT(allocVector(REALSXP, n));
double *v = REAL(val);
if (*diag_P(x) == 'U') {
int j;
for (j = 0; j < n; j++) v[j] = 1.;
} else {
d_packed_getDiag(v, x, n);
}
UNPROTECT(1);
return;
}
void tr_l_packed_getDiag( int *dest, SEXP x)
{
int n = INTEGER(GET_SLOT(x, Matrix_DimSym))[0];
SEXP val = PROTECT(allocVector(LGLSXP, n));
int *v = LOGICAL(val);
if (*diag_P(x) == 'U') {
int j;
for (j = 0; j < n; j++) v[j] = 1;
} else {
l_packed_getDiag(v, x, n);
}
UNPROTECT(1);
return;
}
SEXP Matrix_expand_pointers(SEXP pP)
{
int n = length(pP) - 1;
int *p = INTEGER(pP);
SEXP ans = PROTECT(allocVector(INTSXP, p[n]));
expand_cmprPt(n, p, INTEGER(ans));
UNPROTECT(1);
return ans;
}
/**
* Return the element of a given name from a named list
*
* @param list
* @param nm name of desired element
*
* @return element of list with name nm
*/
SEXP
Matrix_getElement(SEXP list, char *nm) {
SEXP names = getAttrib(list, R_NamesSymbol);
int i;
for (i = 0; i < LENGTH(names); i++)
if (!strcmp(CHAR(STRING_ELT(names, i)), nm))
return(VECTOR_ELT(list, i));
return R_NilValue;
}
/**
* Zero a square matrix of size nc then copy a vector to the diagonal
*
* @param dest destination array of length nc * nc
* @param A pointer to a square Matrix object
*
* @return dest
*/
static double *
install_diagonal(double *dest, SEXP A)
{
int nc = INTEGER(GET_SLOT(A, Matrix_DimSym))[0];
int i, ncp1 = nc + 1, unit = *diag_P(A) == 'U';
double *ax = REAL(GET_SLOT(A, Matrix_xSym));
AZERO(dest, nc * nc);
for (i = 0; i < nc; i++)
dest[i * ncp1] = (unit) ? 1. : ax[i];
return dest;
}
static int *
install_diagonal_int(int *dest, SEXP A)
{
int nc = INTEGER(GET_SLOT(A, Matrix_DimSym))[0];
int i, ncp1 = nc + 1, unit = *diag_P(A) == 'U';
int *ax = INTEGER(GET_SLOT(A, Matrix_xSym));
AZERO(dest, nc * nc);
for (i = 0; i < nc; i++)
dest[i * ncp1] = (unit) ? 1 : ax[i];
return dest;
}
/** Duplicate a [dln]denseMatrix or a numeric matrix or vector
* as a [dln]geMatrix.
* This is for the many *_matrix_{prod,crossprod,tcrossprod,etc.}
* functions that work with both classed and unclassed matrices.
*
* @param A either a denseMatrix object or a matrix object
*/
/* NOTA BENE: If you enlarge this list, do change '14' and '6' below !*/
#define ddense_CLASSES \
"dgeMatrix", "dtrMatrix", \
"dsyMatrix", "dpoMatrix", "ddiMatrix", \
"dtpMatrix", "dspMatrix", "dppMatrix", \
/* sub classes of those above:*/ \
/* dtr */ "Cholesky", "LDL", "BunchKaufman", \
/* dtp */ "pCholesky", "pBunchKaufman", \
/* dpo */ "corMatrix"
#define ldense_CLASSES \
"lgeMatrix", "ltrMatrix", \
"lsyMatrix", "ldiMatrix", \
"ltpMatrix", "lspMatrix"
#define ndense_CLASSES \
"ngeMatrix", "ntrMatrix", \
"nsyMatrix", \
"ntpMatrix", "nspMatrix"
/* Generalized -- "geMatrix" -- dispatch where needed : */
SEXP dup_mMatrix_as_geMatrix(SEXP A)
{
SEXP ans, ad = R_NilValue, an = R_NilValue; /* -Wall */
const char *cl = class_P(A);
char *valid[] = {"_NOT_A_CLASS_",/* *_CLASSES defined in ./Mutils.h */
ddense_CLASSES, /* 14 */
ldense_CLASSES, /* 6 */
ndense_CLASSES, /* 5 */
""};
int sz, ctype = Matrix_check_class(cl, valid), nprot = 1;
enum dense_enum { ddense, ldense, ndense } M_type = ddense /* -Wall */;
if (ctype > 0) { /* a [nld]denseMatrix object */
ad = GET_SLOT(A, Matrix_DimSym);
an = GET_SLOT(A, Matrix_DimNamesSym);
M_type = (ctype <= 14) ? ddense :
((ctype <= 14+6) ? ldense : ndense);
}
else if (ctype < 0) { /* not a (recognized) classed matrix */
if (isReal(A))
M_type = ddense;
else if (isLogical(A))
M_type = ldense;
else
error(_("invalid class `%s' to dup_mMatrix_as_geMatrix"), cl);
#define DUP_MMATRIX_NON_CLASS \
if (isMatrix(A)) { /* "matrix" */ \
ad = getAttrib(A, R_DimSymbol); \
an = getAttrib(A, R_DimNamesSymbol); \
} else {/* maybe "numeric" (incl integer,logical) --> (n x 1) */\
int* dd = INTEGER(ad = PROTECT(allocVector(INTSXP, 2))); \
nprot++; \
dd[0] = LENGTH(A); \
dd[1] = 1; \
an = R_NilValue; \
} \
ctype = 0
DUP_MMATRIX_NON_CLASS;
}
ans = PROTECT(NEW_OBJECT(MAKE_CLASS(M_type == ddense ? "dgeMatrix" :
(M_type == ldense ? "lgeMatrix" :
"ngeMatrix"))));
#define DUP_MMATRIX_SET_1 \
SET_SLOT(ans, Matrix_DimSym, duplicate(ad)); \
SET_SLOT(ans, Matrix_DimNamesSym, (LENGTH(an) == 2) ? \
duplicate(an): allocVector(VECSXP, 2)); \
sz = INTEGER(ad)[0] * INTEGER(ad)[1]
DUP_MMATRIX_SET_1;
if(M_type == ddense) { /* ddense -> dge */
double *ansx;
#define DUP_MMATRIX_ddense_CASES \
ansx = REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, sz)); \
switch(ctype) { \
case 0: /* unclassed real/logical matrix */ \
Memcpy(ansx, REAL(A), sz); \
break; \
case 1: /* dgeMatrix */ \
Memcpy(ansx, REAL(GET_SLOT(A, Matrix_xSym)), sz); \
break; \
case 2: /* dtrMatrix and subclasses */ \
case 9: case 10: case 11: /* --- Cholesky, LDL, BunchKaufman */ \
Memcpy(ansx, REAL(GET_SLOT(A, Matrix_xSym)), sz); \
make_d_matrix_triangular(ansx, A); \
break; \
case 3: /* dsyMatrix */ \
case 4: /* dpoMatrix + subclass */ \
case 14: /* --- corMatrix */ \
Memcpy(ansx, REAL(GET_SLOT(A, Matrix_xSym)), sz); \
make_d_matrix_symmetric(ansx, A); \
break; \
case 5: /* ddiMatrix */ \
install_diagonal(ansx, A); \
break; \
case 6: /* dtpMatrix + subclasses */ \
case 12: case 13: /* --- pCholesky, pBunchKaufman */ \
packed_to_full_double(ansx, REAL(GET_SLOT(A, Matrix_xSym)), \
INTEGER(ad)[0], \
*uplo_P(A) == 'U' ? UPP : LOW); \
make_d_matrix_triangular(ansx, A); \
break; \
case 7: /* dspMatrix */ \
case 8: /* dppMatrix */ \
packed_to_full_double(ansx, REAL(GET_SLOT(A, Matrix_xSym)), \
INTEGER(ad)[0], \
*uplo_P(A) == 'U' ? UPP : LOW); \
make_d_matrix_symmetric(ansx, A); \
break; \
} /* switch(ctype) */
DUP_MMATRIX_ddense_CASES
}
else { /* M_type == ldense || M_type = ndense */
/* ldense -> lge */
/* ndense -> nge */
int *ansx = LOGICAL(ALLOC_SLOT(ans, Matrix_xSym, LGLSXP, sz));
switch(ctype) {
case 0: /* unclassed real/logical matrix */
Memcpy(ansx, LOGICAL(A), sz);
break;
case 1+14: /* lgeMatrix */
case 1+14+6: /* ngeMatrix */
Memcpy(ansx, LOGICAL(GET_SLOT(A, Matrix_xSym)), sz);
break;
case 2+14: /* ltrMatrix */
case 2+14+6: /* ntrMatrix */
Memcpy(ansx, LOGICAL(GET_SLOT(A, Matrix_xSym)), sz);
make_i_matrix_triangular(ansx, A);
break;
case 3+14: /* lsyMatrix */
case 3+14+6: /* nsyMatrix */
Memcpy(ansx, LOGICAL(GET_SLOT(A, Matrix_xSym)), sz);
make_i_matrix_symmetric(ansx, A);
break;
case 4+14: /* ldiMatrix */
case 4+14+6: /* ndiMatrix */
install_diagonal_int(ansx, A);
break;
case 5+14: /* ltpMatrix */
case 5+14+6: /* ntpMatrix */
packed_to_full_int(ansx, LOGICAL(GET_SLOT(A, Matrix_xSym)),
INTEGER(ad)[0],
*uplo_P(A) == 'U' ? UPP : LOW);
make_i_matrix_triangular(ansx, A);
break;
case 6+14: /* lspMatrix */
case 6+14+6: /* nspMatrix */
packed_to_full_int(ansx, LOGICAL(GET_SLOT(A, Matrix_xSym)),
INTEGER(ad)[0],
*uplo_P(A) == 'U' ? UPP : LOW);
make_i_matrix_symmetric(ansx, A);
break;
default:
error(_("unexpected ctype = %d in dup_mMatrix_as_geMatrix"), ctype);
} /* switch(ctype) */
} /* if(M_type == .) */
UNPROTECT(nprot);
return ans;
}
SEXP dup_mMatrix_as_dgeMatrix(SEXP A)
{
SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix"))),
ad = R_NilValue , an = R_NilValue; /* -Wall */
const char *cl = class_P(A);
char *valid[] = {"_NOT_A_CLASS_", ddense_CLASSES, ""};
int ctype = Matrix_check_class(cl, valid), nprot = 1, sz;
double *ansx;
if (ctype > 0) { /* a ddenseMatrix object */
ad = GET_SLOT(A, Matrix_DimSym);
an = GET_SLOT(A, Matrix_DimNamesSym);
}
else if (ctype < 0) { /* not a (recognized) classed matrix */
DUP_MMATRIX_NON_CLASS;
if (isInteger(A) || isLogical(A)) {
A = PROTECT(coerceVector(A, REALSXP));
nprot++;
}
if (!isReal(A))
error(_("invalid class `%s' to dup_mMatrix_as_dgeMatrix"), cl);
}
DUP_MMATRIX_SET_1;
DUP_MMATRIX_ddense_CASES
UNPROTECT(nprot);
return ans;
}
SEXP new_dgeMatrix(int nrow, int ncol)
{
SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix"))),
ad = PROTECT(allocVector(INTSXP, 2));
double *ansx;
int sz = nrow * ncol;
INTEGER(ad)[0] = nrow;
INTEGER(ad)[1] = ncol;
SET_SLOT(ans, Matrix_DimSym, ad);
SET_SLOT(ans, Matrix_DimNamesSym, allocVector(VECSXP, 2));
ansx = REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, sz));
UNPROTECT(2);
return ans;
}
