The R Project SVN
R-packages
(root)/branches/Matrix-APIchange/src/Mutils.c – Rev 4558
Rev 4418 | Go to most recent revision | Blame | Compare with Previous | Last modification | View Log | Download | RSS feed
#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;}SEXPset_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;}#if 0/** The following two csc_ functions are identically usable for rcs__** Check for unsorted columns in the row indices** @param ncol number of columns* @param p column pointers* @param i row indices** @return 0 if all columns are sorted, otherwise 1*/int csc_unsorted_columns(int ncol, const int p[], const int i[]){int j;for (j = 0; j < ncol; j++) {int ind, lst = p[j+1] - 1;for (ind = p[j]; ind < lst; ind++) {if (i[ind] > i[ind+1]) return 1;}}return 0;}/*** Sort the columns in a sparse column-oriented matrix so that each* column is in increasing order of row index.** @param ncol number of columns* @param p column pointers* @param i row indices* @param x values of nonzero elements*/void csc_sort_columns(int ncol, const int p[], int i[], double x[]){int j, maxdiff, *ord;double *dd = (double *) NULL;maxdiff = -1;for (j = 0; j < ncol; j++) {int diff = p[j+1] - p[j];if (diff > maxdiff) maxdiff = diff;}ord = Calloc(maxdiff, int);if (x) dd = Calloc(maxdiff, double);for (j = 0; j < ncol; j++) {int cLen = p[j+1] - p[j];if (cLen > 1) {int k, offset = p[j];for (k = 0; k < cLen; k++) ord[k] = k;R_qsort_int_I(i + offset, ord, 1, cLen);if (x) {for (k = 0; k < cLen; k++) dd[k] = x[ord[k] + offset];Memcpy(x + offset, dd, cLen);}}}Free(ord);if (x) Free(dd);}#endif#if 0/*** Check for sorted columns in an object that inherits from the* dgCMatrix class. Resort the columns if necessary.** @param m pointer to an object that inherits from the dgCMatrix class** @return m with the columns sorted by increasing row index*/SEXP csc_check_column_sorting(SEXP m){int *mp = INTEGER(GET_SLOT(m, Matrix_pSym)),*mi = INTEGER(GET_SLOT(m, Matrix_iSym)),ncol = INTEGER(GET_SLOT(m, Matrix_DimSym))[1];if (csc_unsorted_columns(ncol, mp, mi))csc_sort_columns(ncol, mp, mi, REAL(GET_SLOT(m, Matrix_xSym)));return m;}#endif/* 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*/SEXPMatrix_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 4096SEXP 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 = Calloc(Matrix_Error_Bufsiz, char); sprintfif (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);Free(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** @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]; \} \returnEND_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_getDiagvoid 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*/SEXPMatrix_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 src diagonal elements in an array of length nc* @param nc number of columns (and rows) in the matrix** @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;elseerror(_("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 = 0DUP_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_CASESUNPROTECT(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;}