SCM Repository
[matrix] / pkg / Matrix / src / dgCMatrix.c |
View of /pkg/Matrix/src/dgCMatrix.c
Parent Directory
|
Revision Log
Revision 3270 -
(download)
(as text)
(annotate)
Fri Mar 23 08:50:48 2018 UTC (10 months, 3 weeks ago) by mmaechler
File size: 22732 byte(s)
Fri Mar 23 08:50:48 2018 UTC (10 months, 3 weeks ago) by mmaechler
File size: 22732 byte(s)
provide NEW_OBJECT_OF_CLASS() wrapper. ==> Safe NEW_OBJECT(MAKE_CLASS())
#ifdef __GLIBC__ // to get strdup declared in glibc (when strict -std=c11 or -stdc99): #define _POSIX_C_SOURCE 200809L #endif #include <string.h> #include "dgCMatrix.h" /* for Csparse_transpose() : */ #include "Csparse.h" #include "chm_common.h" /* -> Mutils.h / SPQR ... */ /* FIXME -- we "forget" about dimnames almost everywhere : */ /* for dgCMatrix _and_ lgCMatrix and others (but *not* ngC...) : */ SEXP xCMatrix_validate(SEXP x) { /* Almost everything now in Csparse_validate ( ./Csparse.c ) * *but* the checking of the 'x' slot : */ if (xlength(GET_SLOT(x, Matrix_iSym)) != xlength(GET_SLOT(x, Matrix_xSym))) return mkString(_("lengths of slots 'i' and 'x' must match")); return ScalarLogical(1); } /* for dgRMatrix _and_ lgRMatrix and others (but *not* ngC...) : */ SEXP xRMatrix_validate(SEXP x) { /* Almost everything now in Rsparse_validate ( ./Csparse.c ) * *but* the checking of the 'x' slot : */ if (xlength(GET_SLOT(x, Matrix_jSym)) != xlength(GET_SLOT(x, Matrix_xSym))) return mkString(_("lengths of slots 'j' and 'x' must match")); return ScalarLogical(1); } /* This and the following R_to_CMatrix() lead to memory-not-mapped seg.faults * only with {32bit + R-devel + enable-R-shlib} -- no idea why */ SEXP compressed_to_TMatrix(SEXP x, SEXP colP) { int col = asLogical(colP); /* 1 if "C"olumn compressed; 0 if "R"ow */ /* however, for Csparse, we now effectively use the cholmod-based * Csparse_to_Tsparse() in ./Csparse.c ; maybe should simply write * an as_cholmod_Rsparse() function and then do "as there" ...*/ SEXP indSym = col ? Matrix_iSym : Matrix_jSym, ans, indP = PROTECT(GET_SLOT(x, indSym)), pP = PROTECT(GET_SLOT(x, Matrix_pSym)); int npt = length(pP) - 1; char *ncl = strdup(class_P(x)); static const char *valid[] = { MATRIX_VALID_Csparse, MATRIX_VALID_Rsparse, ""}; int ctype = R_check_class_etc(x, valid); if (ctype < 0) error(_("invalid class(x) '%s' in compressed_to_TMatrix(x)"), ncl); /* replace 'C' or 'R' with 'T' :*/ ncl[2] = 'T'; ans = PROTECT(NEW_OBJECT_OF_CLASS(ncl)); slot_dup(ans, x, Matrix_DimSym); if((ctype / 3) % 4 != 2) /* not n..Matrix */ slot_dup(ans, x, Matrix_xSym); if(ctype % 3) { /* s(ymmetric) or t(riangular) : */ slot_dup(ans, x, Matrix_uploSym); if(ctype % 3 == 2) /* t(riangular) : */ slot_dup(ans, x, Matrix_diagSym); } SET_DimNames(ans, x); // possibly asymmetric for symmetricMatrix is ok SET_SLOT(ans, indSym, duplicate(indP)); expand_cmprPt(npt, INTEGER(pP), INTEGER(ALLOC_SLOT(ans, col ? Matrix_jSym : Matrix_iSym, INTSXP, length(indP)))); free(ncl); UNPROTECT(3); return ans; } SEXP R_to_CMatrix(SEXP x) { SEXP ans, tri = PROTECT(allocVector(LGLSXP, 1)); char *ncl = strdup(class_P(x)); static const char *valid[] = { MATRIX_VALID_Rsparse, ""}; int ctype = R_check_class_etc(x, valid); int *x_dims = INTEGER(GET_SLOT(x, Matrix_DimSym)), *a_dims; PROTECT_INDEX ipx; if (ctype < 0) error(_("invalid class(x) '%s' in R_to_CMatrix(x)"), ncl); /* replace 'R' with 'C' : */ ncl[2] = 'C'; PROTECT_WITH_INDEX(ans = NEW_OBJECT_OF_CLASS(ncl), &ipx); a_dims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2)); /* reversed dim() since we will transpose: */ a_dims[0] = x_dims[1]; a_dims[1] = x_dims[0]; /* triangular: */ LOGICAL(tri)[0] = 0; if((ctype / 3) != 2) /* not n..Matrix */ slot_dup(ans, x, Matrix_xSym); if(ctype % 3) { /* s(ymmetric) or t(riangular) : */ SET_SLOT(ans, Matrix_uploSym, mkString((*uplo_P(x) == 'U') ? "L" : "U")); if(ctype % 3 == 2) { /* t(riangular) : */ LOGICAL(tri)[0] = 1; slot_dup(ans, x, Matrix_diagSym); } } SET_SLOT(ans, Matrix_iSym, duplicate(GET_SLOT(x, Matrix_jSym))); slot_dup(ans, x, Matrix_pSym); REPROTECT(ans = Csparse_transpose(ans, tri), ipx); SET_DimNames(ans, x); // possibly asymmetric for symmetricMatrix is ok free(ncl); UNPROTECT(2); return ans; } /** Return a 2 column matrix '' cbind(i, j) '' of 0-origin index vectors (i,j) * which entirely correspond to the (i,j) slots of * as(x, "TsparseMatrix") : */ SEXP compressed_non_0_ij(SEXP x, SEXP colP) { int col = asLogical(colP); /* 1 if "C"olumn compressed; 0 if "R"ow */ SEXP ans, indSym = col ? Matrix_iSym : Matrix_jSym; SEXP indP = PROTECT(GET_SLOT(x, indSym)), pP = PROTECT(GET_SLOT(x, Matrix_pSym)); int i, *ij; int nouter = INTEGER(GET_SLOT(x, Matrix_DimSym))[col ? 1 : 0], n_el = INTEGER(pP)[nouter]; /* is only == length(indP), if the inner slot is not over-allocated */ ij = INTEGER(ans = PROTECT(allocMatrix(INTSXP, n_el, 2))); /* expand the compressed margin to 'i' or 'j' : */ expand_cmprPt(nouter, INTEGER(pP), &ij[col ? n_el : 0]); /* and copy the other one: */ if (col) for(i = 0; i < n_el; i++) ij[i] = INTEGER(indP)[i]; else /* row compressed */ for(i = 0; i < n_el; i++) ij[i + n_el] = INTEGER(indP)[i]; UNPROTECT(3); return ans; } #if 0 /* unused */ SEXP dgCMatrix_lusol(SEXP x, SEXP y) { SEXP ycp = PROTECT((TYPEOF(y) == REALSXP) ? duplicate(y) : coerceVector(y, REALSXP)); CSP xc = AS_CSP__(x); R_CheckStack(); if (xc->m != xc->n || xc->m <= 0) error(_("dgCMatrix_lusol requires a square, non-empty matrix")); if (LENGTH(ycp) != xc->m) error(_("Dimensions of system to be solved are inconsistent")); if (!cs_lusol(/*order*/ 1, xc, REAL(ycp), /*tol*/ 1e-7)) error(_("cs_lusol failed")); UNPROTECT(1); return ycp; } #endif // called from package MatrixModels's R code SEXP dgCMatrix_qrsol(SEXP x, SEXP y, SEXP ord) { /* FIXME: extend this to work in multivariate case, i.e. y a matrix with > 1 column ! */ SEXP ycp = PROTECT((TYPEOF(y) == REALSXP) ? duplicate(y) : coerceVector(y, REALSXP)); CSP xc = AS_CSP(x); /* <--> x may be dgC* or dtC* */ int order = asInteger(ord); #ifdef _not_yet_do_FIXME__ const char *nms[] = {"L", "coef", "Xty", "resid", ""}; SEXP ans = PROTECT(Rf_mkNamed(VECSXP, nms)); #endif R_CheckStack(); if (order < 0 || order > 3) error(_("dgCMatrix_qrsol(., order) needs order in {0,..,3}")); /* --> cs_amd() --- order 0: natural, 1: Chol, 2: LU, 3: QR */ if (LENGTH(ycp) != xc->m) error(_("Dimensions of system to be solved are inconsistent")); /* FIXME? Note that qr_sol() would allow *under-determined systems; * In general, we'd need LENGTH(ycp) = max(n,m) * FIXME also: multivariate y (see above) */ if (xc->m < xc->n || xc->n <= 0) error(_("dgCMatrix_qrsol(<%d x %d>-matrix) requires a 'tall' rectangular matrix"), xc->m, xc->n); /* cs_qrsol(): Tim Davis (2006) .. "8.2 Using a QR factorization", p.136f , calling * ------- cs_sqr(order, ..), see p.76 */ /* MM: FIXME: write our *OWN* version of - the first case (m >= n) - of cs_qrsol() * --------- which will (1) work with a *multivariate* y * (2) compute coefficients properly, not overwriting RHS */ if (!cs_qrsol(order, xc, REAL(ycp))) /* return value really is 0 or 1 - no more info there */ error(_("cs_qrsol() failed inside dgCMatrix_qrsol()")); /* Solution is only in the first part of ycp -- cut its length back to n : */ ycp = lengthgets(ycp, (R_xlen_t) xc->n); UNPROTECT(1); return ycp; } // Modified version of Tim Davis's cs_qr_mex.c file for MATLAB (in CSparse) // Usage: [V,beta,p,R,q] = cs_qr(A) ; SEXP dgCMatrix_QR(SEXP Ap, SEXP order, SEXP keep_dimnames) { CSP A = AS_CSP__(Ap), D; int io = INTEGER(order)[0]; Rboolean verbose = (io < 0);// verbose=TRUE, encoded with negative 'order' int m0 = A->m, m = m0, n = A->n, ord = asLogical(order) ? 3 : 0, *p; R_CheckStack(); if (m < n) error(_("A must have #{rows} >= #{columns}")) ; SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("sparseQR"))); int *dims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2)); dims[0] = m; dims[1] = n; css *S = cs_sqr(ord, A, 1); /* symbolic QR ordering & analysis*/ if (!S) error(_("cs_sqr failed")); int keep_dimnms = asLogical(keep_dimnames); if(keep_dimnms == NA_LOGICAL) { keep_dimnms = TRUE; warning(_("dgcMatrix_QR(*, keep_dimnames = NA): NA taken as TRUE")); } if(verbose && S->m2 > m) // in ./cs.h , m2 := # of rows for QR, after adding fictitious rows Rprintf("Symbolic QR(): Matrix structurally rank deficient (m2-m = %d)\n", S->m2 - m); csn *N = cs_qr(A, S); /* numeric QR factorization */ if (!N) error(_("cs_qr failed")) ; cs_dropzeros(N->L); /* drop zeros from V and sort */ D = cs_transpose(N->L, 1); cs_spfree(N->L); N->L = cs_transpose(D, 1); cs_spfree(D); cs_dropzeros(N->U); /* drop zeros from R and sort */ D = cs_transpose(N->U, 1); cs_spfree(N->U) ; N->U = cs_transpose(D, 1); cs_spfree(D); m = N->L->m; /* m may be larger now */ // MM: m := S->m2 also counting the ficticious rows (Tim Davis, p.72, 74f) p = cs_pinv(S->pinv, m); /* p = pinv' */ SEXP dn = R_NilValue; Rboolean do_dn = FALSE; if(keep_dimnms) { dn = GET_SLOT(Ap, Matrix_DimNamesSym); do_dn = !isNull(VECTOR_ELT(dn, 0)) && m == m0; // FIXME? also deal with case m > m0 ? if(do_dn) { // keep rownames dn = PROTECT(duplicate(dn)); SET_VECTOR_ELT(dn, 1, R_NilValue); } else dn = R_NilValue; } SET_SLOT(ans, Matrix_VSym, Matrix_cs_to_SEXP(N->L, "dgCMatrix", 0, dn)); // "V" Memcpy(REAL(ALLOC_SLOT(ans, Matrix_betaSym, REALSXP, n)), N->B, n); Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, m)), p, m); if(do_dn) { UNPROTECT(1); // dn dn = R_NilValue; do_dn = FALSE; } if (ord) { Memcpy(INTEGER(ALLOC_SLOT(ans, install("q"), INTSXP, n)), S->q, n); if(keep_dimnms) { dn = GET_SLOT(Ap, Matrix_DimNamesSym); do_dn = !isNull(VECTOR_ELT(dn, 1)); if(do_dn) { dn = PROTECT(duplicate(dn)); // permute colnames by S->q : cn <- cn[ S->q ] : SEXP cns = PROTECT(duplicate(VECTOR_ELT(dn, 1))); for(int j=0; j < n; j++) SET_STRING_ELT(VECTOR_ELT(dn, 1), j, STRING_ELT(cns, S->q[j])); UNPROTECT(1); SET_VECTOR_ELT(dn, 0, R_NilValue); } else dn = R_NilValue; } } else ALLOC_SLOT(ans, install("q"), INTSXP, 0); SET_SLOT(ans, install("R"), Matrix_cs_to_SEXP(N->U, "dgCMatrix", 0, dn)); if(do_dn) UNPROTECT(1); // dn cs_nfree(N); cs_sfree(S); cs_free(p); UNPROTECT(1); return ans; } #ifdef Matrix_with_SPQR /** * Return a SuiteSparse QR factorization of the sparse matrix A * * @param Ap (pointer to) a [m x n] dgCMatrix * @param ordering integer SEXP specifying the ordering strategy to be used * see SPQR/Include/SuiteSparseQR_definitions.h * @param econ integer SEXP ("economy"): number of rows of R and columns of Q * to return. The default is m. Using n gives the standard economy form. * A value less than the estimated rank r is set to r, so econ=0 gives the * "rank-sized" factorization, where nrow(R)==nnz(diag(R))==r. * @param tol double SEXP: if tol <= -2 use SPQR's default, * if -2 < tol < 0, then no tol is used; otherwise, * tol > 0, use as tolerance: columns with 2-norm <= tol treated as 0 * * * @return SEXP "SPQR" object with slots (Q, R, p, rank, Dim): * Q: dgCMatrix; R: dgCMatrix [subject to change to dtCMatrix FIXME ?] * p: integer: 0-based permutation (or length 0 <=> identity); * rank: integer, the "revealed" rank Dim: integer, original matrix dim. */ SEXP dgCMatrix_SPQR(SEXP Ap, SEXP ordering, SEXP econ, SEXP tol) { /* SEXP ans = PROTECT(allocVector(VECSXP, 4)); */ SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("SPQR"))); CHM_SP A = AS_CHM_SP(Ap), Q, R; SuiteSparse_long *E, rank;/* not always = int FIXME (Windows_64 ?) */ if ((rank = SuiteSparseQR_C_QR(asInteger(ordering), asReal(tol),/* originally had SPQR_DEFAULT_TOL */ (SuiteSparse_long)asInteger(econ),/* originally had 0 */ A, &Q, &R, &E, &cl)) == -1) error(_("SuiteSparseQR_C_QR returned an error code")); slot_dup(ans, Ap, Matrix_DimSym); /* SET_VECTOR_ELT(ans, 0, */ /* chm_sparse_to_SEXP(Q, 0, 0, 0, "", R_NilValue)); */ SET_SLOT(ans, install("Q"), chm_sparse_to_SEXP(Q, 0, 0, 0, "", R_NilValue)); /* Also gives a dgCMatrix (not a dtC* *triangular*) : * may make sense if to be used in the "spqr_solve" routines .. ?? */ /* SET_VECTOR_ELT(ans, 1, */ /* chm_sparse_to_SEXP(R, 0, 0, 0, "", R_NilValue)); */ SET_SLOT(ans, install("R"), chm_sparse_to_SEXP(R, 0, 0, 0, "", R_NilValue)); cholmod_free_sparse(&Al, &cl); cholmod_free_sparse(&R, &cl); cholmod_free_sparse(&Q, &cl); if (E) { int *Er; SET_VECTOR_ELT(ans, 2, allocVector(INTSXP, A->ncol)); Er = INTEGER(VECTOR_ELT(ans, 2)); for (int i = 0; i < A->ncol; i++) Er[i] = (int) E[i]; Free(E); } else SET_VECTOR_ELT(ans, 2, allocVector(INTSXP, 0)); SET_VECTOR_ELT(ans, 3, ScalarInteger((int)rank)); UNPROTECT(1); return ans; } #endif /* Matrix_with_SPQR */ /* Modified version of Tim Davis's cs_lu_mex.c file for MATLAB */ void install_lu(SEXP Ap, int order, double tol, Rboolean err_sing, Rboolean keep_dimnms) { // (order, tol) == (1, 1) by default, when called from R. SEXP ans; css *S; csn *N; int n, *p, *dims; CSP A = AS_CSP__(Ap), D; R_CheckStack(); n = A->n; if (A->m != n) error(_("LU decomposition applies only to square matrices")); if (order) { /* not using natural order */ order = (tol == 1) ? 2 /* amd(S'*S) w/dense rows or I */ : 1; /* amd (A+A'), or natural */ } S = cs_sqr(order, A, /*qr = */ 0); /* symbolic ordering */ N = cs_lu(A, S, tol); /* numeric factorization */ if (!N) { if(err_sing) error(_("cs_lu(A) failed: near-singular A (or out of memory)")); else { /* No warning: The useR should be careful : * Put NA into "LU" factor */ set_factors(Ap, ScalarLogical(NA_LOGICAL), "LU"); return; } } cs_dropzeros(N->L); /* drop zeros from L and sort it */ D = cs_transpose(N->L, 1); cs_spfree(N->L); N->L = cs_transpose(D, 1); cs_spfree(D); cs_dropzeros(N->U); /* drop zeros from U and sort it */ D = cs_transpose(N->U, 1); cs_spfree(N->U); N->U = cs_transpose(D, 1); cs_spfree(D); p = cs_pinv(N->pinv, n); /* p=pinv' */ ans = PROTECT(NEW_OBJECT(MAKE_CLASS("sparseLU"))); dims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2)); dims[0] = n; dims[1] = n; SEXP dn; Rboolean do_dn = FALSE; if(keep_dimnms) { dn = GET_SLOT(Ap, Matrix_DimNamesSym); do_dn = !isNull(VECTOR_ELT(dn, 0)); if(do_dn) { dn = PROTECT(duplicate(dn)); // permute rownames by p : rn <- rn[ p ] : SEXP rn = PROTECT(duplicate(VECTOR_ELT(dn, 0))); for(int i=0; i < n; i++) SET_STRING_ELT(VECTOR_ELT(dn, 0), i, STRING_ELT(rn, p[i])); UNPROTECT(1); // rn SET_VECTOR_ELT(dn, 1, R_NilValue); // colnames(.) := NULL } } SET_SLOT(ans, install("L"), Matrix_cs_to_SEXP(N->L, "dtCMatrix", 0, do_dn ? dn : R_NilValue)); if(keep_dimnms) { if(do_dn) { UNPROTECT(1); // dn dn = GET_SLOT(Ap, Matrix_DimNamesSym); } do_dn = !isNull(VECTOR_ELT(dn, 1)); if(do_dn) { dn = PROTECT(duplicate(dn)); if(order) { // permute colnames by S->q : cn <- cn[ S->q ] : SEXP cn = PROTECT(duplicate(VECTOR_ELT(dn, 1))); for(int j=0; j < n; j++) SET_STRING_ELT(VECTOR_ELT(dn, 1), j, STRING_ELT(cn, S->q[j])); UNPROTECT(1); // cn } SET_VECTOR_ELT(dn, 0, R_NilValue); // rownames(.) := NULL } } SET_SLOT(ans, install("U"), Matrix_cs_to_SEXP(N->U, "dtCMatrix", 0, do_dn ? dn : R_NilValue)); if(do_dn) UNPROTECT(1); // dn Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_pSym, /* "p" */ INTSXP, n)), p, n); if (order) Memcpy(INTEGER(ALLOC_SLOT(ans, install("q"), INTSXP, n)), S->q, n); cs_nfree(N); cs_sfree(S); cs_free(p); UNPROTECT(1); set_factors(Ap, ans, "LU"); } SEXP dgCMatrix_LU(SEXP Ap, SEXP orderp, SEXP tolp, SEXP error_on_sing, SEXP keep_dimnames) { SEXP ans; Rboolean err_sing = asLogical(error_on_sing); /* FIXME: dgCMatrix_LU should check ans for consistency in * permutation type with the requested value - Should have two * classes or two different names in the factors list for LU with * permuted columns or not. * OTOH, currently (order, tol) === (1, 1) always. * It is true that length(LU@q) does flag the order argument. */ if (!isNull(ans = get_factors(Ap, "LU"))) return ans; int keep_dimnms = asLogical(keep_dimnames); if(keep_dimnms == NA_LOGICAL) { keep_dimnms = TRUE; warning(_("dgcMatrix_LU(*, keep_dimnames = NA): NA taken as TRUE")); } install_lu(Ap, asInteger(orderp), asReal(tolp), err_sing, keep_dimnms); return get_factors(Ap, "LU"); } SEXP dgCMatrix_matrix_solve(SEXP Ap, SEXP b, SEXP give_sparse) // FIXME: add 'keep_dimnames' as argument { Rboolean sparse = asLogical(give_sparse); if(sparse) { // FIXME: implement this error(_("dgCMatrix_matrix_solve(.., sparse=TRUE) not yet implemented")); /* Idea: in the for(j = 0; j < nrhs ..) loop below, build the *sparse* result matrix * ----- *column* wise -- which is perfect for dgCMatrix * --> build (i,p,x) slots "increasingly" [well, allocate in batches ..] * * --> maybe first a protoype in R */ } SEXP ans = PROTECT(dup_mMatrix_as_dgeMatrix(b)), lu, qslot; CSP L, U; int *bdims = INTEGER(GET_SLOT(ans, Matrix_DimSym)), *p, *q; int j, n = bdims[0], nrhs = bdims[1]; double *x, *ax = REAL(GET_SLOT(ans, Matrix_xSym)); C_or_Alloca_TO(x, n, double); if (isNull(lu = get_factors(Ap, "LU"))) { install_lu(Ap, /* order = */ 1, /* tol = */ 1.0, /* err_sing = */ TRUE, /* keep_dimnames = */ TRUE); lu = get_factors(Ap, "LU"); } qslot = GET_SLOT(lu, install("q")); L = AS_CSP__(GET_SLOT(lu, install("L"))); U = AS_CSP__(GET_SLOT(lu, install("U"))); R_CheckStack(); if (U->n != n) error(_("Dimensions of system to be solved are inconsistent")); if(nrhs >= 1 && n >= 1) { p = INTEGER(GET_SLOT(lu, Matrix_pSym)); q = LENGTH(qslot) ? INTEGER(qslot) : (int *) NULL; for (j = 0; j < nrhs; j++) { cs_pvec(p, ax + j * n, x, n); /* x = b(p) */ cs_lsolve(L, x); /* x = L\x */ cs_usolve(U, x); /* x = U\x */ if (q) /* r(q) = x , hence r = Q' U{^-1} L{^-1} P b = A^{-1} b */ cs_ipvec(q, x, ax + j * n, n); else Memcpy(ax + j * n, x, n); } } if(n >= SMALL_4_Alloca) Free(x); UNPROTECT(1); return ans; } // called from package MatrixModels's R code: SEXP dgCMatrix_cholsol(SEXP x, SEXP y) { /* Solve Sparse Least Squares X %*% beta ~= y with dense RHS y, * where X = t(x) i.e. we pass x = t(X) as argument, * via "Cholesky(X'X)" .. well not really: * cholmod_factorize("x", ..) finds L in X'X = L'L directly */ CHM_SP cx = AS_CHM_SP(x); /* FIXME: extend this to work in multivariate case, i.e. y a matrix with > 1 column ! */ CHM_DN cy = AS_CHM_DN(coerceVector(y, REALSXP)), rhs, cAns, resid; CHM_FR L; int n = cx->ncol;/* #{obs.} {x = t(X) !} */ double one[] = {1,0}, zero[] = {0,0}, neg1[] = {-1,0}; const char *nms[] = {"L", "coef", "Xty", "resid", ""}; SEXP ans = PROTECT(Rf_mkNamed(VECSXP, nms)); R_CheckStack(); if (n < cx->nrow || n <= 0) error(_("dgCMatrix_cholsol requires a 'short, wide' rectangular matrix")); if (cy->nrow != n) error(_("Dimensions of system to be solved are inconsistent")); rhs = cholmod_allocate_dense(cx->nrow, 1, cx->nrow, CHOLMOD_REAL, &c); /* cholmod_sdmult(A, transp, alpha, beta, X, Y, &c): * Y := alpha*(A*X) + beta*Y or alpha*(A'*X) + beta*Y ; * here: rhs := 1 * x %*% y + 0 = x %*% y = X'y */ if (!(cholmod_sdmult(cx, 0 /* trans */, one, zero, cy, rhs, &c))) error(_("cholmod_sdmult error (rhs)")); L = cholmod_analyze(cx, &c); if (!cholmod_factorize(cx, L, &c)) error(_("cholmod_factorize failed: status %d, minor %d from ncol %d"), c.status, L->minor, L->n); /* FIXME: Do this in stages so an "effects" vector can be calculated */ if (!(cAns = cholmod_solve(CHOLMOD_A, L, rhs, &c))) error(_("cholmod_solve (CHOLMOD_A) failed: status %d, minor %d from ncol %d"), c.status, L->minor, L->n); /* L : */ SET_VECTOR_ELT(ans, 0, chm_factor_to_SEXP(L, 0)); /* coef : */ SET_VECTOR_ELT(ans, 1, allocVector(REALSXP, cx->nrow)); Memcpy(REAL(VECTOR_ELT(ans, 1)), (double*)(cAns->x), cx->nrow); /* X'y : */ /* FIXME: Change this when the "effects" vector is available */ SET_VECTOR_ELT(ans, 2, allocVector(REALSXP, cx->nrow)); Memcpy(REAL(VECTOR_ELT(ans, 2)), (double*)(rhs->x), cx->nrow); /* resid := y */ resid = cholmod_copy_dense(cy, &c); /* cholmod_sdmult(A, transp, alp, bet, X, Y, &c): * Y := alp*(A*X) + bet*Y or alp*(A'*X) + beta*Y ; * here: resid := -1 * x' %*% coef + 1 * y = y - X %*% coef */ if (!(cholmod_sdmult(cx, 1/* trans */, neg1, one, cAns, resid, &c))) error(_("cholmod_sdmult error (resid)")); /* FIXME: for multivariate case, i.e. resid *matrix* with > 1 column ! */ SET_VECTOR_ELT(ans, 3, allocVector(REALSXP, n)); Memcpy(REAL(VECTOR_ELT(ans, 3)), (double*)(resid->x), n); cholmod_free_factor(&L, &c); cholmod_free_dense(&rhs, &c); cholmod_free_dense(&cAns, &c); UNPROTECT(1); return ans; } /* Define all of * dgCMatrix_colSums(....) * igCMatrix_colSums(....) * lgCMatrix_colSums_d(....) * lgCMatrix_colSums_i(....) * ngCMatrix_colSums_d(....) * ngCMatrix_colSums_i(....) */ #define _dgC_ #include "t_gCMatrix_colSums.c" #define _igC_ #include "t_gCMatrix_colSums.c" #define _lgC_ #include "t_gCMatrix_colSums.c" #define _ngC_ #include "t_gCMatrix_colSums.c" #define _lgC_mn #include "t_gCMatrix_colSums.c" #define _ngC_mn #include "t_gCMatrix_colSums.c" SEXP lgCMatrix_colSums(SEXP x, SEXP NArm, SEXP spRes, SEXP trans, SEXP means) { if(asLogical(means)) /* ==> result will be "double" / "dsparseVector" */ return lgCMatrix_colSums_d(x, NArm, spRes, trans, means); else return lgCMatrix_colSums_i(x, NArm, spRes, trans, means); } SEXP ngCMatrix_colSums(SEXP x, SEXP NArm, SEXP spRes, SEXP trans, SEXP means) { if(asLogical(means)) /* ==> result will be "double" / "dsparseVector" */ return ngCMatrix_colSums_d(x, NArm, spRes, trans, means); else return ngCMatrix_colSums_i(x, NArm, spRes, trans, means); }
R-Forge@R-project.org | ViewVC Help |
Powered by ViewVC 1.0.0 |