- Generated on Thu Oct 6 2011 for Matrix by
1.7.4
|
Matrix $Rev: 2718 $ at $LastChangedDate: 2011-10-06 11:45:17 +0200 (Thu, 06 Oct 2011) $
|
00001 #include <limits.h> 00002 00003 #include "Mutils.h" 00004 #include <R_ext/Lapack.h> 00005 00006 #if 0 /* defined(R_VERSION) && R_VERSION >= R_Version(2, 7, 0) * 00007 * La_norm_type() & La_rcond_type() are now in R_ext/Lapack.h 00008 * but because of the 'module-mess' that's not sufficient */ 00009 #else 00010 char La_norm_type(const char *typstr) 00011 { 00012 char typup; 00013 00014 if (strlen(typstr) != 1) 00015 error( 00016 _("argument type[1]='%s' must be a one-letter character string"), 00017 typstr); 00018 typup = toupper(*typstr); 00019 if (typup == '1') 00020 typup = 'O'; /* alias */ 00021 else if (typup == 'E') 00022 typup = 'F'; 00023 else if (typup != 'M' && typup != 'O' && typup != 'I' && typup != 'F') 00024 error(_("argument type[1]='%s' must be one of 'M','1','O','I','F' or 'E'"), 00025 typstr); 00026 return typup; 00027 } 00028 00029 char La_rcond_type(const char *typstr) 00030 { 00031 char typup; 00032 00033 if (strlen(typstr) != 1) 00034 error( 00035 _("argument type[1]='%s' must be a one-letter character string"), 00036 typstr); 00037 typup = toupper(*typstr); 00038 if (typup == '1') 00039 typup = 'O'; /* alias */ 00040 else if (typup != 'O' && typup != 'I') 00041 error(_("argument type[1]='%s' must be one of '1','O', or 'I'"), 00042 typstr); 00043 return typup; 00044 } 00045 #endif 00046 00047 double get_double_by_name(SEXP obj, char *nm) 00048 { 00049 SEXP nms = getAttrib(obj, R_NamesSymbol); 00050 int i, len = length(obj); 00051 00052 if ((!isReal(obj)) || (length(obj) > 0 && nms == R_NilValue)) 00053 error(_("object must be a named, numeric vector")); 00054 for (i = 0; i < len; i++) { 00055 if (!strcmp(nm, CHAR(STRING_ELT(nms, i)))) { 00056 return REAL(obj)[i]; 00057 } 00058 } 00059 return R_NaReal; 00060 } 00061 00062 SEXP 00063 set_double_by_name(SEXP obj, double val, char *nm) 00064 { 00065 SEXP nms = getAttrib(obj, R_NamesSymbol); 00066 int i, len = length(obj); 00067 00068 if ((!isReal(obj)) || (length(obj) > 0 && nms == R_NilValue)) 00069 error(_("object must be a named, numeric vector")); 00070 for (i = 0; i < len; i++) { 00071 if (!strcmp(nm, CHAR(STRING_ELT(nms, i)))) { 00072 REAL(obj)[i] = val; 00073 return obj; 00074 } 00075 } 00076 { 00077 SEXP nx = PROTECT(allocVector(REALSXP, len + 1)), 00078 nnms = allocVector(STRSXP, len + 1); 00079 00080 setAttrib(nx, R_NamesSymbol, nnms); 00081 for (i = 0; i < len; i++) { 00082 REAL(nx)[i] = REAL(obj)[i]; 00083 SET_STRING_ELT(nnms, i, duplicate(STRING_ELT(nms, i))); 00084 } 00085 REAL(nx)[len] = val; 00086 SET_STRING_ELT(nnms, len, mkChar(nm)); 00087 UNPROTECT(1); 00088 return nx; 00089 } 00090 } 00091 00092 SEXP as_det_obj(double val, int log, int sign) 00093 { 00094 SEXP det = PROTECT(allocVector(VECSXP, 2)), 00095 nms = PROTECT(allocVector(STRSXP, 2)), 00096 vv = PROTECT(ScalarReal(val)); 00097 00098 setAttrib(det, R_NamesSymbol, nms); 00099 SET_STRING_ELT(nms, 0, mkChar("modulus")); 00100 SET_STRING_ELT(nms, 1, mkChar("sign")); 00101 setAttrib(vv, install("logarithm"), ScalarLogical(log)); 00102 SET_VECTOR_ELT(det, 0, vv); 00103 SET_VECTOR_ELT(det, 1, ScalarInteger(sign)); 00104 setAttrib(det, R_ClassSymbol, mkString("det")); 00105 UNPROTECT(3); 00106 return det; 00107 } 00108 00109 SEXP get_factors(SEXP obj, char *nm) 00110 { 00111 SEXP fac = GET_SLOT(obj, Matrix_factorSym), 00112 nms = getAttrib(fac, R_NamesSymbol); 00113 int i, len = length(fac); 00114 00115 if ((!isNewList(fac)) || (length(fac) > 0 && nms == R_NilValue)) 00116 error(_("'factors' slot must be a named list")); 00117 for (i = 0; i < len; i++) { 00118 if (!strcmp(nm, CHAR(STRING_ELT(nms, i)))) { 00119 return VECTOR_ELT(fac, i); 00120 } 00121 } 00122 return R_NilValue; 00123 } 00124 00125 /* In the past this function installed a duplicate of val in the 00126 * factors slot for obj then returned the (typically unprotected) 00127 * val. This is now changed to return the duplicate, which will be 00128 * protected if obj is protected. */ 00129 SEXP set_factors(SEXP obj, SEXP val, char *nm) 00130 { 00131 SEXP fac = GET_SLOT(obj, Matrix_factorSym), 00132 nms = getAttrib(fac, R_NamesSymbol), nfac, nnms; 00133 int i, len = length(fac); 00134 00135 if ((!isNewList(fac)) || (length(fac) > 0 && nms == R_NilValue)) 00136 error(_("'factors' slot must be a named list")); 00137 PROTECT(val); /* set_factors(..) may be called as "finalizer" after UNPROTECT()*/ 00138 for (i = 0; i < len; i++) { 00139 if (!strcmp(nm, CHAR(STRING_ELT(nms, i)))) { 00140 SET_VECTOR_ELT(fac, i, duplicate(val)); 00141 return val; 00142 } 00143 } 00144 nfac = PROTECT(allocVector(VECSXP, len + 1)); 00145 nnms = PROTECT(allocVector(STRSXP, len + 1)); 00146 setAttrib(nfac, R_NamesSymbol, nnms); 00147 for (i = 0; i < len; i++) { 00148 SET_VECTOR_ELT(nfac, i, VECTOR_ELT(fac, i)); 00149 SET_STRING_ELT(nnms, i, duplicate(STRING_ELT(nms, i))); 00150 } 00151 SET_VECTOR_ELT(nfac, len, duplicate(val)); 00152 SET_STRING_ELT(nnms, len, mkChar(nm)); 00153 SET_SLOT(obj, Matrix_factorSym, nfac); 00154 UNPROTECT(3); 00155 return VECTOR_ELT(nfac, len); 00156 } 00157 00158 #if 0 /* unused */ 00159 /* useful for all the ..CMatrix classes (and ..R by [0] <-> [1]); but unused */ 00160 SEXP CMatrix_set_Dim(SEXP x, int nrow) 00161 { 00162 int *dims = INTEGER(GET_SLOT(x, Matrix_DimSym)); 00163 00164 dims[0] = nrow; 00165 dims[1] = length(GET_SLOT(x, Matrix_pSym)) - 1; 00166 return x; 00167 } 00168 #endif /* unused */ 00169 00170 /* Fill in the "trivial remainder" in n*m array ; 00171 * typically the 'x' slot of a "dtrMatrix" : 00172 * But should be usable for double/logical/int/complex : */ 00173 00174 #define MAKE_TRIANGULAR_BODY(_TO_, _FROM_, _ZERO_, _ONE_) \ 00175 { \ 00176 int i, j, *dims = INTEGER(GET_SLOT(_FROM_, Matrix_DimSym)); \ 00177 int n = dims[0], m = dims[1]; \ 00178 \ 00179 if (*uplo_P(_FROM_) == 'U') { \ 00180 for (j = 0; j < n; j++) \ 00181 for (i = j+1; i < m; i++) \ 00182 _TO_[i + j*m] = _ZERO_; \ 00183 } else { \ 00184 for (j = 1; j < n; j++) \ 00185 for (i = 0; i < j && i < m; i++) \ 00186 _TO_[i + j*m] = _ZERO_; \ 00187 } \ 00188 if (*diag_P(_FROM_) == 'U') { \ 00189 j = (n < m) ? n : m; \ 00190 for (i = 0; i < j; i++) \ 00191 _TO_[i * (m + 1)] = _ONE_; \ 00192 } \ 00193 } 00194 00195 void make_d_matrix_triangular(double *to, SEXP from) 00196 MAKE_TRIANGULAR_BODY(to, from, 0., 1.) 00197 void make_i_matrix_triangular(int *to, SEXP from) 00198 MAKE_TRIANGULAR_BODY(to, from, 0, 1) 00199 00200 00201 /* Should work for double/logical/int/complex : */ 00202 #define MAKE_SYMMETRIC_BODY(_TO_, _FROM_) \ 00203 { \ 00204 int i, j, n = INTEGER(GET_SLOT(_FROM_, Matrix_DimSym))[0]; \ 00205 \ 00206 if (*uplo_P(_FROM_) == 'U') { \ 00207 for (j = 0; j < n; j++) \ 00208 for (i = j+1; i < n; i++) \ 00209 _TO_[i + j*n] = _TO_[j + i*n]; \ 00210 } else { \ 00211 for (j = 1; j < n; j++) \ 00212 for (i = 0; i < j && i < n; i++) \ 00213 _TO_[i + j*n] = _TO_[j + i*n]; \ 00214 } \ 00215 } 00216 00217 void make_d_matrix_symmetric(double *to, SEXP from) 00218 MAKE_SYMMETRIC_BODY(to, from) 00219 00220 void make_i_matrix_symmetric(int *to, SEXP from) 00221 MAKE_SYMMETRIC_BODY(to, from) 00222 00223 00224 #define Matrix_Error_Bufsiz 4096 00225 00236 SEXP check_scalar_string(SEXP sP, char *vals, char *nm) 00237 { 00238 SEXP val = ScalarLogical(1); 00239 char *buf; 00240 /* only allocate when needed: in good case, none is needed */ 00241 #define SPRINTF buf = Alloca(Matrix_Error_Bufsiz, char); R_CheckStack(); sprintf 00242 00243 if (length(sP) != 1) { 00244 SPRINTF(buf, _("'%s' slot must have length 1"), nm); 00245 } else { 00246 const char *str = CHAR(STRING_ELT(sP, 0)); 00247 if (strlen(str) != 1) { 00248 SPRINTF(buf, _("'%s' must have string length 1"), nm); 00249 } else { 00250 int i, len; 00251 for (i = 0, len = strlen(vals); i < len; i++) { 00252 if (str[0] == vals[i]) 00253 return R_NilValue; 00254 } 00255 SPRINTF(buf, _("'%s' must be in '%s'"), nm, vals); 00256 } 00257 } 00258 /* 'error' returns : */ 00259 val = mkString(buf); 00260 return val; 00261 #undef SPRINTF 00262 } 00263 00264 /* FIXME? Something like this should be part of the R API ? 00265 * But then, R has the more general compute_identical() 00266 * in src/main/identical.c: Rboolean compute_identical(SEXP x, SEXP y); 00267 */ 00268 Rboolean equal_string_vectors(SEXP s1, SEXP s2) 00269 { 00270 Rboolean n1 = isNull(s1), n2 = isNull(s2); 00271 if (n1 || n2) 00272 return (n1 == n2) ? TRUE : FALSE; 00273 else if (TYPEOF(s1) != STRSXP || TYPEOF(s2) != STRSXP) { 00274 error(_("'s1' and 's2' must be \"character\" vectors")); 00275 return FALSE; /* -Wall */ 00276 } else { 00277 int n = LENGTH(s1), i; 00278 if (n != LENGTH(s2)) 00279 return FALSE; 00280 for(i = 0; i < n; i++) { 00281 /* note that compute_identical() code for STRSXP 00282 is careful about NA's which we don't need */ 00283 if(strcmp(CHAR(STRING_ELT(s1, i)), 00284 CHAR(STRING_ELT(s2, i)))) 00285 return FALSE; 00286 } 00287 return TRUE; /* they *are* equal */ 00288 } 00289 } 00290 00291 00292 SEXP dense_nonpacked_validate(SEXP obj) 00293 { 00294 int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym)); 00295 if ((dims[0] * dims[1]) != length(GET_SLOT(obj, Matrix_xSym))) 00296 return mkString(_("length of x slot != prod(Dim)")); 00297 return ScalarLogical(1); 00298 } 00299 00300 00301 #define PACKED_TO_FULL(TYPE) \ 00302 TYPE *packed_to_full_ ## TYPE(TYPE *dest, const TYPE *src, \ 00303 int n, enum CBLAS_UPLO uplo) \ 00304 { \ 00305 int i, j, pos = 0; \ 00306 \ 00307 AZERO(dest, n*n); \ 00308 for (j = 0; j < n; j++) { \ 00309 switch(uplo) { \ 00310 case UPP: \ 00311 for (i = 0; i <= j; i++) dest[i + j * n] = src[pos++]; \ 00312 break; \ 00313 case LOW: \ 00314 for (i = j; i < n; i++) dest[i + j * n] = src[pos++]; \ 00315 break; \ 00316 default: \ 00317 error(_("'uplo' must be UPP or LOW")); \ 00318 } \ 00319 } \ 00320 return dest; \ 00321 } 00322 00323 PACKED_TO_FULL(double) 00324 PACKED_TO_FULL(int) 00325 00326 #define FULL_TO_PACKED(TYPE) \ 00327 TYPE *full_to_packed_ ## TYPE(TYPE *dest, const TYPE *src, int n, \ 00328 enum CBLAS_UPLO uplo, enum CBLAS_DIAG diag) \ 00329 { \ 00330 int i, j, pos = 0; \ 00331 \ 00332 for (j = 0; j < n; j++) { \ 00333 switch(uplo) { \ 00334 case UPP: \ 00335 for (i = 0; i <= j; i++) \ 00336 dest[pos++] = (i == j && diag== UNT) ? 1 : src[i + j*n]; \ 00337 break; \ 00338 case LOW: \ 00339 for (i = j; i < n; i++) \ 00340 dest[pos++] = (i == j && diag== UNT) ? 1 : src[i + j*n]; \ 00341 break; \ 00342 default: \ 00343 error(_("'uplo' must be UPP or LOW")); \ 00344 } \ 00345 } \ 00346 return dest; \ 00347 } 00348 00349 FULL_TO_PACKED(double) 00350 FULL_TO_PACKED(int) 00351 00352 00353 00363 void d_packed_getDiag(double *dest, SEXP x, int n) 00364 { 00365 double *xx = REAL(GET_SLOT(x, Matrix_xSym)); 00366 00367 #define END_packed_getDiag \ 00368 int j, pos = 0; \ 00369 \ 00370 if (*uplo_P(x) == 'U') { \ 00371 for(pos= 0, j=0; j < n; pos += 1+(++j)) dest[j] = xx[pos]; \ 00372 } else { \ 00373 for(pos= 0, j=0; j < n; pos += (n - j), j++) dest[j] = xx[pos]; \ 00374 } \ 00375 return 00376 00377 END_packed_getDiag; 00378 } 00379 00380 void l_packed_getDiag(int *dest, SEXP x, int n) 00381 { 00382 int *xx = LOGICAL(GET_SLOT(x, Matrix_xSym)); 00383 00384 END_packed_getDiag; 00385 } 00386 00387 #undef END_packed_getDiag 00388 00389 void tr_d_packed_getDiag(double *dest, SEXP x) 00390 { 00391 int n = INTEGER(GET_SLOT(x, Matrix_DimSym))[0]; 00392 SEXP val = PROTECT(allocVector(REALSXP, n)); 00393 double *v = REAL(val); 00394 00395 if (*diag_P(x) == 'U') { 00396 int j; 00397 for (j = 0; j < n; j++) v[j] = 1.; 00398 } else { 00399 d_packed_getDiag(v, x, n); 00400 } 00401 UNPROTECT(1); 00402 return; 00403 } 00404 00405 void tr_l_packed_getDiag( int *dest, SEXP x) 00406 { 00407 int n = INTEGER(GET_SLOT(x, Matrix_DimSym))[0]; 00408 SEXP val = PROTECT(allocVector(LGLSXP, n)); 00409 int *v = LOGICAL(val); 00410 00411 if (*diag_P(x) == 'U') { 00412 int j; 00413 for (j = 0; j < n; j++) v[j] = 1; 00414 } else { 00415 l_packed_getDiag(v, x, n); 00416 } 00417 UNPROTECT(1); 00418 return; 00419 } 00420 00421 SEXP Matrix_expand_pointers(SEXP pP) 00422 { 00423 int n = length(pP) - 1; 00424 int *p = INTEGER(pP); 00425 SEXP ans = PROTECT(allocVector(INTSXP, p[n])); 00426 00427 expand_cmprPt(n, p, INTEGER(ans)); 00428 UNPROTECT(1); 00429 return ans; 00430 } 00431 00432 00441 SEXP 00442 Matrix_getElement(SEXP list, char *nm) { 00443 SEXP names = getAttrib(list, R_NamesSymbol); 00444 int i; 00445 00446 for (i = 0; i < LENGTH(names); i++) 00447 if (!strcmp(CHAR(STRING_ELT(names, i)), nm)) 00448 return(VECTOR_ELT(list, i)); 00449 return R_NilValue; 00450 } 00451 00460 static double * 00461 install_diagonal(double *dest, SEXP A) 00462 { 00463 int nc = INTEGER(GET_SLOT(A, Matrix_DimSym))[0]; 00464 int i, ncp1 = nc + 1, unit = *diag_P(A) == 'U'; 00465 double *ax = REAL(GET_SLOT(A, Matrix_xSym)); 00466 00467 AZERO(dest, nc * nc); 00468 for (i = 0; i < nc; i++) 00469 dest[i * ncp1] = (unit) ? 1. : ax[i]; 00470 return dest; 00471 } 00472 00473 static int * 00474 install_diagonal_int(int *dest, SEXP A) 00475 { 00476 int nc = INTEGER(GET_SLOT(A, Matrix_DimSym))[0]; 00477 int i, ncp1 = nc + 1, unit = *diag_P(A) == 'U'; 00478 int *ax = INTEGER(GET_SLOT(A, Matrix_xSym)); 00479 00480 AZERO(dest, nc * nc); 00481 for (i = 0; i < nc; i++) 00482 dest[i * ncp1] = (unit) ? 1 : ax[i]; 00483 return dest; 00484 } 00485 00486 00494 /* NOTA BENE: If you enlarge this list, do change '14' and '6' below ! 00495 * --------- ddiMatrix & ldiMatrix are no longer ddense or ldense on the R level, 00496 * --- --- but are still dealt with here. 00497 */ 00498 00499 #define ddense_CLASSES \ 00500 "dgeMatrix", "dtrMatrix", \ 00501 "dsyMatrix", "dpoMatrix", "ddiMatrix", \ 00502 "dtpMatrix", "dspMatrix", "dppMatrix", \ 00503 /* sub classes of those above:*/ \ 00504 /* dtr */ "Cholesky", "LDL", "BunchKaufman", \ 00505 /* dtp */ "pCholesky", "pBunchKaufman", \ 00506 /* dpo */ "corMatrix" 00507 00508 #define ldense_CLASSES \ 00509 "lgeMatrix", "ltrMatrix", \ 00510 "lsyMatrix", "ldiMatrix", \ 00511 "ltpMatrix", "lspMatrix" 00512 00513 #define ndense_CLASSES \ 00514 "ngeMatrix", "ntrMatrix", \ 00515 "nsyMatrix", \ 00516 "ntpMatrix", "nspMatrix" 00517 00518 /* Generalized -- "geMatrix" -- dispatch where needed : */ 00519 SEXP dup_mMatrix_as_geMatrix(SEXP A) 00520 { 00521 SEXP ans, ad = R_NilValue, an = R_NilValue; /* -Wall */ 00522 static const char *valid[] = { 00523 "_NOT_A_CLASS_",/* *_CLASSES defined in ./Mutils.h */ 00524 ddense_CLASSES, /* 14 */ 00525 ldense_CLASSES, /* 6 */ 00526 ndense_CLASSES, /* 5 */ 00527 ""}; 00528 int sz, ctype = Matrix_check_class_etc(A, valid), 00529 nprot = 1; 00530 enum dense_enum { ddense, ldense, ndense } M_type = ddense /* -Wall */; 00531 00532 if (ctype > 0) { /* a [nld]denseMatrix object */ 00533 ad = GET_SLOT(A, Matrix_DimSym); 00534 an = GET_SLOT(A, Matrix_DimNamesSym); 00535 M_type = (ctype <= 14) ? ddense : 00536 ((ctype <= 14+6) ? ldense : ndense); 00537 } 00538 else if (ctype < 0) { /* not a (recognized) classed matrix */ 00539 00540 if (isReal(A)) 00541 M_type = ddense; 00542 else if (isInteger(A)) { 00543 A = PROTECT(coerceVector(A, REALSXP)); 00544 nprot++; 00545 M_type = ddense; 00546 } 00547 else if (isLogical(A)) 00548 M_type = ldense; 00549 else 00550 error(_("invalid class '%s' to dup_mMatrix_as_geMatrix"), 00551 class_P(A)); 00552 00553 #define DUP_MMATRIX_NON_CLASS \ 00554 if (isMatrix(A)) { /* "matrix" */ \ 00555 ad = getAttrib(A, R_DimSymbol); \ 00556 an = getAttrib(A, R_DimNamesSymbol); \ 00557 } else {/* maybe "numeric" (incl integer,logical) --> (n x 1) */\ 00558 int* dd = INTEGER(ad = PROTECT(allocVector(INTSXP, 2))); \ 00559 nprot++; \ 00560 dd[0] = LENGTH(A); \ 00561 dd[1] = 1; \ 00562 an = R_NilValue; \ 00563 } \ 00564 ctype = 0 00565 00566 DUP_MMATRIX_NON_CLASS; 00567 } 00568 00569 ans = PROTECT(NEW_OBJECT(MAKE_CLASS(M_type == ddense ? "dgeMatrix" : 00570 (M_type == ldense ? "lgeMatrix" : 00571 "ngeMatrix")))); 00572 #define DUP_MMATRIX_SET_1 \ 00573 SET_SLOT(ans, Matrix_DimSym, duplicate(ad)); \ 00574 SET_SLOT(ans, Matrix_DimNamesSym, (LENGTH(an) == 2) ? \ 00575 duplicate(an): allocVector(VECSXP, 2)); \ 00576 sz = INTEGER(ad)[0] * INTEGER(ad)[1] 00577 00578 DUP_MMATRIX_SET_1; 00579 00580 if(M_type == ddense) { /* ddense -> dge */ 00581 00582 double *ansx; 00583 00584 #define DUP_MMATRIX_ddense_CASES \ 00585 ansx = REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, sz)); \ 00586 switch(ctype) { \ 00587 case 0: /* unclassed real matrix */ \ 00588 Memcpy(ansx, REAL(A), sz); \ 00589 break; \ 00590 case 1: /* dgeMatrix */ \ 00591 Memcpy(ansx, REAL(GET_SLOT(A, Matrix_xSym)), sz); \ 00592 break; \ 00593 case 2: /* dtrMatrix and subclasses */ \ 00594 case 9: case 10: case 11: /* --- Cholesky, LDL, BunchKaufman */ \ 00595 Memcpy(ansx, REAL(GET_SLOT(A, Matrix_xSym)), sz); \ 00596 make_d_matrix_triangular(ansx, A); \ 00597 break; \ 00598 case 3: /* dsyMatrix */ \ 00599 case 4: /* dpoMatrix + subclass */ \ 00600 case 14: /* --- corMatrix */ \ 00601 Memcpy(ansx, REAL(GET_SLOT(A, Matrix_xSym)), sz); \ 00602 make_d_matrix_symmetric(ansx, A); \ 00603 break; \ 00604 case 5: /* ddiMatrix */ \ 00605 install_diagonal(ansx, A); \ 00606 break; \ 00607 case 6: /* dtpMatrix + subclasses */ \ 00608 case 12: case 13: /* --- pCholesky, pBunchKaufman */ \ 00609 packed_to_full_double(ansx, REAL(GET_SLOT(A, Matrix_xSym)), \ 00610 INTEGER(ad)[0], \ 00611 *uplo_P(A) == 'U' ? UPP : LOW); \ 00612 make_d_matrix_triangular(ansx, A); \ 00613 break; \ 00614 case 7: /* dspMatrix */ \ 00615 case 8: /* dppMatrix */ \ 00616 packed_to_full_double(ansx, REAL(GET_SLOT(A, Matrix_xSym)), \ 00617 INTEGER(ad)[0], \ 00618 *uplo_P(A) == 'U' ? UPP : LOW); \ 00619 make_d_matrix_symmetric(ansx, A); \ 00620 break; \ 00621 } /* switch(ctype) */ 00622 00623 DUP_MMATRIX_ddense_CASES 00624 00625 } 00626 else { /* M_type == ldense || M_type = ndense */ 00627 /* ldense -> lge */ 00628 /* ndense -> nge */ 00629 int *ansx = LOGICAL(ALLOC_SLOT(ans, Matrix_xSym, LGLSXP, sz)); 00630 00631 switch(ctype) { 00632 case 0: /* unclassed logical matrix */ 00633 Memcpy(ansx, LOGICAL(A), sz); 00634 break; 00635 00636 case 1+14: /* lgeMatrix */ 00637 case 1+14+6: /* ngeMatrix */ 00638 Memcpy(ansx, LOGICAL(GET_SLOT(A, Matrix_xSym)), sz); 00639 break; 00640 case 2+14: /* ltrMatrix */ 00641 case 2+14+6: /* ntrMatrix */ 00642 Memcpy(ansx, LOGICAL(GET_SLOT(A, Matrix_xSym)), sz); 00643 make_i_matrix_triangular(ansx, A); 00644 break; 00645 case 3+14: /* lsyMatrix */ 00646 case 3+14+6: /* nsyMatrix */ 00647 Memcpy(ansx, LOGICAL(GET_SLOT(A, Matrix_xSym)), sz); 00648 make_i_matrix_symmetric(ansx, A); 00649 break; 00650 case 4+14: /* ldiMatrix */ 00651 case 4+14+6: /* ndiMatrix */ 00652 install_diagonal_int(ansx, A); 00653 break; 00654 case 5+14: /* ltpMatrix */ 00655 case 5+14+6: /* ntpMatrix */ 00656 packed_to_full_int(ansx, LOGICAL(GET_SLOT(A, Matrix_xSym)), 00657 INTEGER(ad)[0], 00658 *uplo_P(A) == 'U' ? UPP : LOW); 00659 make_i_matrix_triangular(ansx, A); 00660 break; 00661 case 6+14: /* lspMatrix */ 00662 case 6+14+6: /* nspMatrix */ 00663 packed_to_full_int(ansx, LOGICAL(GET_SLOT(A, Matrix_xSym)), 00664 INTEGER(ad)[0], 00665 *uplo_P(A) == 'U' ? UPP : LOW); 00666 make_i_matrix_symmetric(ansx, A); 00667 break; 00668 00669 default: 00670 error(_("unexpected ctype = %d in dup_mMatrix_as_geMatrix"), ctype); 00671 } /* switch(ctype) */ 00672 00673 } /* if(M_type == .) */ 00674 00675 UNPROTECT(nprot); 00676 return ans; 00677 } 00678 00679 SEXP dup_mMatrix_as_dgeMatrix(SEXP A) 00680 { 00681 SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix"))), 00682 ad = R_NilValue , an = R_NilValue; /* -Wall */ 00683 static const char *valid[] = {"_NOT_A_CLASS_", ddense_CLASSES, ""}; 00684 int ctype = Matrix_check_class_etc(A, valid), 00685 nprot = 1, sz; 00686 double *ansx; 00687 00688 if (ctype > 0) { /* a ddenseMatrix object */ 00689 ad = GET_SLOT(A, Matrix_DimSym); 00690 an = GET_SLOT(A, Matrix_DimNamesSym); 00691 } 00692 else if (ctype < 0) { /* not a (recognized) classed matrix */ 00693 00694 DUP_MMATRIX_NON_CLASS; 00695 00696 if (isInteger(A) || isLogical(A)) { 00697 A = PROTECT(coerceVector(A, REALSXP)); 00698 nprot++; 00699 } 00700 if (!isReal(A)) 00701 error(_("invalid class '%s' to dup_mMatrix_as_dgeMatrix"), 00702 class_P(A)); 00703 } 00704 00705 DUP_MMATRIX_SET_1; 00706 00707 DUP_MMATRIX_ddense_CASES 00708 00709 UNPROTECT(nprot); 00710 return ans; 00711 } 00712 00713 00714 SEXP new_dgeMatrix(int nrow, int ncol) 00715 { 00716 SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix"))), 00717 ad = PROTECT(allocVector(INTSXP, 2)); 00718 00719 INTEGER(ad)[0] = nrow; 00720 INTEGER(ad)[1] = ncol; 00721 SET_SLOT(ans, Matrix_DimSym, ad); 00722 SET_SLOT(ans, Matrix_DimNamesSym, allocVector(VECSXP, 2)); 00723 ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nrow * ncol); 00724 00725 UNPROTECT(2); 00726 return ans; 00727 } 00728 00738 SEXP m_encodeInd(SEXP ij, SEXP di, SEXP chk_bnds) 00739 { 00740 SEXP ans; 00741 int *ij_di = NULL, n, *Di = INTEGER(di), *IJ, *j_; 00742 Rboolean check_bounds = asLogical(chk_bnds); 00743 00744 ij = PROTECT(coerceVector(ij, INTSXP)); 00745 if(!isMatrix(ij) || 00746 (ij_di = INTEGER(getAttrib(ij, R_DimSymbol)))[1] != 2) 00747 error(_("Argument ij must be 2-column integer matrix")); 00748 n = ij_di[0]; 00749 IJ = INTEGER(ij); 00750 j_ = IJ+n;/* pointer offset! */ 00751 00752 if((Di[0] * (double) Di[1]) >= 1 + (double)INT_MAX) { /* need double */ 00753 ans = PROTECT(allocVector(REALSXP, n)); 00754 double *ii = REAL(ans), nr = (double) Di[0]; 00755 #define do_ii_FILL(_i_, _j_) \ 00756 int i; \ 00757 if(check_bounds) { \ 00758 for(i=0; i < n; i++) { \ 00759 if(_i_[i] == NA_INTEGER || _j_[i] == NA_INTEGER) \ 00760 ii[i] = NA_INTEGER; \ 00761 else { \ 00762 if(_i_[i] < 0 || _i_[i] >= Di[0]) \ 00763 error(_("subscript 'i' out of bounds in M[ij]")); \ 00764 if(_j_[i] < 0 || _j_[i] >= Di[1]) \ 00765 error(_("subscript 'j' out of bounds in M[ij]")); \ 00766 ii[i] = _i_[i] + _j_[i] * nr; \ 00767 } \ 00768 } \ 00769 } else { \ 00770 for(i=0; i < n; i++) \ 00771 ii[i] = (_i_[i] == NA_INTEGER || _j_[i] == NA_INTEGER) \ 00772 ? NA_INTEGER : _i_[i] + _j_[i] * nr; \ 00773 } 00774 00775 do_ii_FILL(IJ, j_); 00776 } else { 00777 ans = PROTECT(allocVector(INTSXP, n)); 00778 int *ii = INTEGER(ans), nr = Di[0]; 00779 00780 do_ii_FILL(IJ, j_); 00781 } 00782 UNPROTECT(2); 00783 return ans; 00784 } 00785 00796 SEXP m_encodeInd2(SEXP i, SEXP j, SEXP di, SEXP chk_bnds) 00797 { 00798 SEXP ans; 00799 int n = LENGTH(i); 00800 int *Di = INTEGER(di), *i_, *j_; 00801 Rboolean check_bounds = asLogical(chk_bnds); 00802 00803 if(LENGTH(j) != n || !isInteger(i) || !isInteger(j)) 00804 error(_("i and j must be integer vectors of the same length")); 00805 i_ = INTEGER(i); 00806 j_ = INTEGER(j); 00807 00808 if((Di[0] * (double) Di[1]) >= 1 + (double)INT_MAX) { /* need double */ 00809 ans = PROTECT(allocVector(REALSXP, n)); 00810 double *ii = REAL(ans), nr = (double) Di[0]; 00811 00812 do_ii_FILL(i_, j_); 00813 } else { 00814 ans = PROTECT(allocVector(INTSXP, n)); 00815 int *ii = INTEGER(ans), nr = Di[0]; 00816 00817 do_ii_FILL(i_, j_); 00818 } 00819 UNPROTECT(1); 00820 return ans; 00821 } 00822 #undef do_ii_FILL 00823 00824 00825 #if R_VERSION < R_Version(2, 13, 0) 00826 00836 int Matrix_check_class_and_super(SEXP x, const char **valid, SEXP rho) 00837 { 00838 int ans; 00839 SEXP cl = getAttrib(x, R_ClassSymbol); 00840 const char *class = CHAR(asChar(cl)); 00841 for (ans = 0; ; ans++) { 00842 if (!strlen(valid[ans])) // empty string 00843 break; 00844 if (!strcmp(class, valid[ans])) return ans; 00845 } 00846 /* if not found directly, now search the non-virtual super classes :*/ 00847 if(IS_S4_OBJECT(x)) { 00848 /* now try the superclasses, i.e., try is(x, "...."); superCl := 00849 .selectSuperClasses(getClass("...")@contains, dropVirtual=TRUE) */ 00850 SEXP classExts, superCl, _call; 00851 int i; 00852 PROTECT(_call = lang2(install("getClassDef"), cl)); 00853 classExts = GET_SLOT(eval(_call, rho), 00854 install("contains")); 00855 UNPROTECT(1); 00856 PROTECT(classExts); 00857 PROTECT(_call = lang3(install(".selectSuperClasses"), classExts, 00858 /* dropVirtual = */ ScalarLogical(1))); 00859 superCl = eval(_call, rho); 00860 UNPROTECT(2); 00861 PROTECT(superCl); 00862 for(i=0; i < length(superCl); i++) { 00863 const char *s_class = CHAR(STRING_ELT(superCl, i)); 00864 for (ans = 0; ; ans++) { 00865 if (!strlen(valid[ans])) 00866 break; 00867 if (!strcmp(s_class, valid[ans])) { 00868 UNPROTECT(1); 00869 return ans; 00870 } 00871 } 00872 } 00873 UNPROTECT(1); 00874 } 00875 return -1; 00876 } 00877 #endif 00878 00889 int Matrix_check_class_etc(SEXP x, const char **valid) 00890 { 00891 static SEXP s_M_classEnv = NULL; 00892 SEXP cl = getAttrib(x, R_ClassSymbol), rho = R_GlobalEnv, pkg; 00893 /* PROTECT(cl); */ 00894 if(!s_M_classEnv) 00895 s_M_classEnv = install(".M.classEnv"); 00896 00897 pkg = getAttrib(cl, R_PackageSymbol); /* ==R== packageSlot(class(x)) */ 00898 if(!isNull(pkg)) { /* find rho := correct class Environment */ 00899 SEXP clEnvCall; 00900 /* need to make sure we find ".M.classEnv" even if Matrix is not 00901 attached, but just namespace-loaded: */ 00902 00903 /* Matrix::: does not work here either ... : 00904 * rho = eval(lang2(install("Matrix:::.M.classEnv"), cl), */ 00905 00906 /* Now make this work via .onLoad() hack in ../R/zzz.R : */ 00907 PROTECT(clEnvCall = lang2(s_M_classEnv, cl)); 00908 rho = eval(clEnvCall, R_GlobalEnv); 00909 UNPROTECT(1); 00910 00911 if(!isEnvironment(rho)) 00912 error(_("could not find correct environment; please report!")); 00913 } 00914 /* UNPROTECT(1); */ 00915 return Matrix_check_class_and_super(x, valid, rho); 00916 }
1.7.4