--- pkg/Matrix/src/Csparse.c 2015/03/26 10:00:49 3069 +++ pkg/Matrix/src/Csparse.c 2015/10/29 16:56:10 3147 @@ -1,5 +1,8 @@ - /* Sparse matrices in compressed column-oriented form */ - +/** @file Csparse.c + * The "CsparseMatrix" class from R package Matrix: + * + * Sparse matrices in compressed column-oriented form + */ #include "Csparse.h" #include "Tsparse.h" #include "chm_common.h" @@ -101,8 +104,8 @@ return ScalarLogical(1); } -/** - * From a CsparseMatrix, produce a dense one. +/** @brief From a CsparseMatrix, produce a dense one. + * * Directly deals with symmetric, triangular and general. * Called from ../R/Csparse.R's C2dense() * @@ -122,7 +125,7 @@ ctype = 0; // <- default = "dgC" static const char *valid[] = { MATRIX_VALID_Csparse, ""}; if(is_sym_or_tri == NA_INTEGER) { // find if is(x, "symmetricMatrix") : - ctype = Matrix_check_class_etc(x, valid); + ctype = R_check_class_etc(x, valid); is_sym = (ctype % 3 == 1); is_tri = (ctype % 3 == 2); } else { @@ -130,7 +133,7 @@ is_tri = is_sym_or_tri < 0; // => both are FALSE iff is_.. == 0 if(is_sym || is_tri) - ctype = Matrix_check_class_etc(x, valid); + ctype = R_check_class_etc(x, valid); } CHM_SP chxs = AS_CHM_SP__(x);// -> chxs->stype = +- 1 <==> symmetric R_CheckStack(); @@ -186,19 +189,27 @@ } // FIXME: do not go via CHM (should not be too hard, to just *drop* the x-slot, right? -SEXP Csparse_to_nz_pattern(SEXP x, SEXP tri) +SEXP Csparse2nz(SEXP x, Rboolean tri) { CHM_SP chxs = AS_CHM_SP__(x); CHM_SP chxcp = cholmod_copy(chxs, chxs->stype, CHOLMOD_PATTERN, &c); - int tr = asLogical(tri); R_CheckStack(); return chm_sparse_to_SEXP(chxcp, 1/*do_free*/, - tr ? ((*uplo_P(x) == 'U') ? 1 : -1) : 0, + tri ? ((*uplo_P(x) == 'U') ? 1 : -1) : 0, /* Rkind: pattern */ 0, - /* diag = */ tr ? diag_P(x) : "", + /* diag = */ tri ? diag_P(x) : "", GET_SLOT(x, Matrix_DimNamesSym)); } +SEXP Csparse_to_nz_pattern(SEXP x, SEXP tri) +{ + int tr_ = asLogical(tri); + if(tr_ == NA_LOGICAL) { + warning(_("Csparse_to_nz_pattern(x, tri = NA): 'tri' is taken as TRUE")); + tr_ = TRUE; + } + return Csparse2nz(x, (Rboolean) tr_); +} // n.CMatrix --> [dli].CMatrix (not going through CHM!) SEXP nz_pattern_to_Csparse(SEXP x, SEXP res_kind) @@ -211,8 +222,16 @@ SEXP nz2Csparse(SEXP x, enum x_slot_kind r_kind) { const char *cl_x = class_P(x); - if(cl_x[0] != 'n') error(_("not a 'n.CMatrix'")); - if(cl_x[2] != 'C') error(_("not a CsparseMatrix")); + // quick check - if ok, fast + if(cl_x[0] != 'n' || cl_x[2] != 'C') { + // e.g. class = "A", from setClass("A", contains = "ngCMatrix") + static const char *valid[] = { MATRIX_VALID_nCsparse, ""}; + int ctype = R_check_class_etc(x, valid); + if(ctype < 0) + error(_("not a 'n.CMatrix'")); + else // fine : get a valid cl_x class_P()-like string : + cl_x = valid[ctype]; + } int nnz = LENGTH(GET_SLOT(x, Matrix_iSym)); SEXP ans; char *ncl = alloca(strlen(cl_x) + 1); /* not much memory required */ @@ -261,7 +280,7 @@ int is_sym = asLogical(symm); if(is_sym == NA_LOGICAL) { // find if is(x, "symmetricMatrix") : static const char *valid[] = { MATRIX_VALID_Csparse, ""}; - int ctype = Matrix_check_class_etc(x, valid); + int ctype = R_check_class_etc(x, valid); is_sym = (ctype % 3 == 1); } return chm_dense_to_matrix( @@ -401,7 +420,15 @@ Rkind, tr ? diag_P(x) : "", dn); } -/* NOTA BENE: cholmod_ssmult(A,B, ...) -> ./CHOLMOD/MatrixOps/cholmod_ssmult.c +/** @brief A %*% B - for matrices of class CsparseMatrix (R package "Matrix") + * + * @param a + * @param b + * @param bool_arith + * + * @return + * + * NOTA BENE: cholmod_ssmult(A,B, ...) -> ./CHOLMOD/MatrixOps/cholmod_ssmult.c * --------- computes a patter*n* matrix __always_ when * *one* of A or B is pattern*n*, because of this (line 73-74): --------------------------------------------------------------------------- @@ -410,14 +437,12 @@ --------------------------------------------------------------------------- * ==> Often need to copy the patter*n* to a *l*ogical matrix first !!! */ - SEXP Csparse_Csparse_prod(SEXP a, SEXP b, SEXP bool_arith) { CHM_SP cha = AS_CHM_SP(a), chb = AS_CHM_SP(b), chc; R_CheckStack(); - // const char *cl_a = class_P(a), *cl_b = class_P(b); static const char *valid_tri[] = { MATRIX_VALID_tri_Csparse, "" }; char diag[] = {'\0', '\0'}; int uploT = 0, nprot = 1, @@ -455,8 +480,8 @@ * Note that in that case, the multiplication itself should happen * faster. But there's no support for that in CHOLMOD */ - if(Matrix_check_class_etc(a, valid_tri) >= 0 && - Matrix_check_class_etc(b, valid_tri) >= 0) + if(R_check_class_etc(a, valid_tri) >= 0 && + R_check_class_etc(b, valid_tri) >= 0) if(*uplo_P(a) == *uplo_P(b)) { /* both upper, or both lower tri. */ uploT = (*uplo_P(a) == 'U') ? 1 : -1; if(*diag_P(a) == 'U' && *diag_P(b) == 'U') { /* return UNIT-triag. */ @@ -476,8 +501,16 @@ return chm_sparse_to_SEXP(chc, 1, uploT, /*Rkind*/0, diag, dn); } -/* trans = FALSE: crossprod(a,b) - * trans = TRUE : tcrossprod(a,b) */ +/** @brief [t]crossprod (, ) + * + * @param a a "CsparseMatrix" object + * @param b a "CsparseMatrix" object + * @param trans trans = FALSE: crossprod(a,b) + * trans = TRUE : tcrossprod(a,b) + * @param bool_arith logical (TRUE / NA / FALSE): Should boolean arithmetic be used. + * + * @return a CsparseMatrix, the (t)cross product of a and b. + */ SEXP Csparse_Csparse_crossprod(SEXP a, SEXP b, SEXP trans, SEXP bool_arith) { int tr = asLogical(trans), nprot = 1, @@ -487,7 +520,6 @@ chb = AS_CHM_SP(b), chTr, chc; R_CheckStack(); - // const char *cl_a = class_P(a), *cl_b = class_P(b); static const char *valid_tri[] = { MATRIX_VALID_tri_Csparse, "" }; char diag[] = {'\0', '\0'}; int uploT = 0; @@ -502,16 +534,25 @@ SEXP da = PROTECT(nz2Csparse(a, x_double)); nprot++; cha = AS_CHM_SP(da); R_CheckStack(); - a_is_n = FALSE; + // a_is_n = FALSE; } else if(b_is_n && (force_num || (maybe_bool && !a_is_n))) { // coerce 'b' to double SEXP db = PROTECT(nz2Csparse(b, x_double)); nprot++; chb = AS_CHM_SP(db); R_CheckStack(); - b_is_n = FALSE; + // b_is_n = FALSE; + } + else if(do_bool == TRUE) { // Want boolean arithmetic: sufficient if *one* is pattern: + if(!a_is_n && !b_is_n) { + // coerce 'a' to pattern + SEXP da = PROTECT(Csparse2nz(a, /* tri = */ + R_check_class_etc(a, valid_tri) >= 0)); nprot++; + cha = AS_CHM_SP(da); + R_CheckStack(); + // a_is_n = TRUE; + } } - chTr = cholmod_transpose((tr) ? chb : cha, chb->xtype, &c); chc = cholmod_ssmult((tr) ? cha : chTr, (tr) ? chTr : chb, /*out_stype:*/ 0, /* values : */ do_bool != TRUE, @@ -520,8 +561,8 @@ /* Preserve triangularity and unit-triangularity if appropriate; * see Csparse_Csparse_prod() for comments */ - if(Matrix_check_class_etc(a, valid_tri) >= 0 && - Matrix_check_class_etc(b, valid_tri) >= 0) + if(R_check_class_etc(a, valid_tri) >= 0 && + R_check_class_etc(b, valid_tri) >= 0) if(*uplo_P(a) != *uplo_P(b)) { /* one 'U', the other 'L' */ uploT = (*uplo_P(b) == 'U') ? 1 : -1; if(*diag_P(a) == 'U' && *diag_P(b) == 'U') { /* return UNIT-triag. */ @@ -545,7 +586,7 @@ /** * All (dense * sparse) Matrix products and cross products * - * f( f() %*% f() ) for f in {t(), identity()} + * f( f() %*% f() ) where f () is either t () [tranpose] or the identity. * * @param a CsparseMatrix (n x m) * @param b numeric vector, matrix, or denseMatrix (m x k) or (k x m) if `transp` is '2' or 'B' @@ -589,9 +630,9 @@ /* repeating a "cheap part" of mMatrix_as_dgeMatrix2(b, .) to see if * we have a vector that we might 'transpose_if_vector' : */ static const char *valid[] = {"_NOT_A_CLASS_", MATRIX_VALID_ddense, ""}; - /* int ctype = Matrix_check_class_etc(b, valid); + /* int ctype = R_check_class_etc(b, valid); * if (ctype > 0) /.* a ddenseMatrix object */ - if (Matrix_check_class_etc(b, valid) < 0) { + if (R_check_class_etc(b, valid) < 0) { // not a ddenseM*: is.matrix() or vector: b_is_vector = !isMatrix(b); } @@ -670,8 +711,9 @@ } -/* Computes x'x or x x' -- *also* for Tsparse (triplet = TRUE) - see Csparse_Csparse_crossprod above for x'y and x y' */ +/** @brief Computes x'x or x x' -- *also* for Tsparse (triplet = TRUE) + see Csparse_Csparse_crossprod above for x'y and x y' +*/ SEXP Csparse_crossprod(SEXP x, SEXP trans, SEXP triplet, SEXP bool_arith) { int tripl = asLogical(triplet), @@ -683,7 +725,7 @@ SEXP xx = PROTECT(Tsparse_diagU2N(x)); CHM_TR cht = tripl ? AS_CHM_TR__(xx) : (CHM_TR) NULL; int nprot = 2; #endif - CHM_SP chcp, chxt, + CHM_SP chcp, chxt, chxc, chx = (tripl ? cholmod_triplet_to_sparse(cht, cht->nnz, &c) : AS_CHM_SP(x)); @@ -691,6 +733,7 @@ R_CheckStack(); Rboolean x_is_n = (chx->xtype == CHOLMOD_PATTERN), + x_is_sym = chx->stype != 0, force_num = (do_bool == FALSE); if(x_is_n && force_num) { @@ -699,11 +742,23 @@ chx = AS_CHM_SP(dx); R_CheckStack(); } + else if(do_bool == TRUE && !x_is_n) { // Want boolean arithmetic; need patter[n] + // coerce 'x' to pattern + static const char *valid_tri[] = { MATRIX_VALID_tri_Csparse, "" }; + SEXP dx = PROTECT(Csparse2nz(x, /* tri = */ + R_check_class_etc(x, valid_tri) >= 0)); nprot++; + chx = AS_CHM_SP(dx); + R_CheckStack(); + } if (!tr) chxt = cholmod_transpose(chx, chx->xtype, &c); - chcp = cholmod_aat((!tr) ? chxt : chx, (int *) NULL, 0, - /* mode: */ chx->xtype, &c); + if (x_is_sym) // cholmod_aat() does not like symmetric + chxc = cholmod_copy(tr ? chx : chxt, /* stype: */ 0, + chx->xtype, &c); + // CHOLMOD/Core/cholmod_aat.c : + chcp = cholmod_aat(x_is_sym ? chxc : (tr ? chx : chxt), + (int *) NULL, 0, /* mode: */ chx->xtype, &c); if(!chcp) { UNPROTECT(1); error(_("Csparse_crossprod(): error return from cholmod_aat()")); @@ -721,13 +776,13 @@ return chm_sparse_to_SEXP(chcp, 1, 0, 0, "", dn); } -/* Csparse_drop(x, tol): drop entries with absolute value < tol, i.e, -* at least all "explicit" zeros */ +/** @brief Csparse_drop(x, tol): drop entries with absolute value < tol, i.e, + * at least all "explicit" zeros. */ SEXP Csparse_drop(SEXP x, SEXP tol) { const char *cl = class_P(x); /* dtCMatrix, etc; [1] = the second character =?= 't' for triangular */ - int tr = (cl[1] == 't'); + int tr = (cl[1] == 't'); // FIXME - rather R_check_class_etc(..) CHM_SP chx = AS_CHM_SP__(x); CHM_SP ans = cholmod_copy(chx, chx->stype, chx->xtype, &c); double dtol = asReal(tol); @@ -742,6 +797,8 @@ GET_SLOT(x, Matrix_DimNamesSym)); } +/** @brief Horizontal Concatenation - cbind( , ) + */ SEXP Csparse_horzcat(SEXP x, SEXP y) { #define CSPARSE_CAT(_KIND_) \ @@ -774,6 +831,8 @@ 1, 0, Rkind, "", R_NilValue); } +/** @brief Vertical Concatenation - rbind( , ) + */ SEXP Csparse_vertcat(SEXP x, SEXP y) { CSPARSE_CAT("vertcat"); @@ -846,10 +905,10 @@ } /** - * "Indexing" aka subsetting : Compute x[i,j], also for vectors i and j + * Indexing aka subsetting : Compute x[i,j], also for vectors i and j * Working via CHOLMOD_submatrix, see ./CHOLMOD/MatrixOps/cholmod_submatrix.c * @param x CsparseMatrix - * @param i row indices (0-origin), or NULL (R's) + * @param i row indices (0-origin), or NULL (R, not C) * @param j columns indices (0-origin), or NULL * * @return x[i,j] still CsparseMatrix --- currently, this loses dimnames