SCM

SCM Repository

[matrix] Diff of /pkg/Matrix/src/Csparse.c
ViewVC logotype

Diff of /pkg/Matrix/src/Csparse.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3147, Thu Oct 29 16:56:10 2015 UTC revision 3270, Fri Mar 23 08:50:48 2018 UTC
# Line 151  Line 151 
151       * BUT, much worse (FIXME!), it also transforms CHOLMOD_PATTERN ("n") matrices       * BUT, much worse (FIXME!), it also transforms CHOLMOD_PATTERN ("n") matrices
152       * to numeric (CHOLMOD_REAL) ones {and we "revert" via chm_dense_to_SEXP()}: */       * to numeric (CHOLMOD_REAL) ones {and we "revert" via chm_dense_to_SEXP()}: */
153      CHM_DN chxd = cholmod_sparse_to_dense(chxs, &c);      CHM_DN chxd = cholmod_sparse_to_dense(chxs, &c);
154        /* FIXME: The above FAILS for prod(dim(.)) > INT_MAX
155         * ----
156         * TODO: use cholmod_l_* but also the 'cl' global ==> many changes in chm_common.[ch]
157         * >>>>>>>>>>> TODO <<<<<<<<<<<<
158         * CHM_DN chxd = cholmod_l_sparse_to_dense(chxs, &cl); */
159        //                   ^^^ important when prod(dim(.)) > INT_MAX
160      int Rkind = (chxs->xtype == CHOLMOD_PATTERN)? -1 : Real_kind(x);      int Rkind = (chxs->xtype == CHOLMOD_PATTERN)? -1 : Real_kind(x);
161    
162      SEXP ans = chm_dense_to_SEXP(chxd, 1, Rkind, GET_SLOT(x, Matrix_DimNamesSym),      SEXP ans = chm_dense_to_SEXP(chxd, 1, Rkind, GET_SLOT(x, Matrix_DimNamesSym),
# Line 240  Line 246 
246      ncl[0] = (r_kind == x_double ? 'd' :      ncl[0] = (r_kind == x_double ? 'd' :
247                (r_kind == x_logical ? 'l' :                (r_kind == x_logical ? 'l' :
248                 /* else (for now):  r_kind == x_integer : */ 'i'));                 /* else (for now):  r_kind == x_integer : */ 'i'));
249      PROTECT(ans = NEW_OBJECT(MAKE_CLASS(ncl)));      PROTECT(ans = NEW_OBJECT_OF_CLASS(ncl));
250      // create a correct 'x' slot:      // create a correct 'x' slot:
251      switch(r_kind) {      switch(r_kind) {
252          int i;          int i;
# Line 343  Line 349 
349      if (!(chx->stype))      if (!(chx->stype))
350          error(_("Nonsymmetric matrix in Csparse_symmetric_to_general"));          error(_("Nonsymmetric matrix in Csparse_symmetric_to_general"));
351      chgx = cholmod_copy(chx, /* stype: */ 0, chx->xtype, &c);      chgx = cholmod_copy(chx, /* stype: */ 0, chx->xtype, &c);
     /* xtype: pattern, "real", complex or .. */  
352      return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "",      return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "",
353                                symmetric_DimNames(GET_SLOT(x, Matrix_DimNamesSym)));                                symmetric_DimNames(GET_SLOT(x, Matrix_DimNamesSym)));
354  }  }
# Line 389  Line 394 
394          }          }
395          UNPROTECT(1);          UNPROTECT(1);
396      }      }
397      /* xtype: pattern, "real", complex or .. */      /* Rkind: pattern, "real", complex or .. */
398      return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "", dns);      return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "", dns);
399  }  }
400    
# Line 407  Line 412 
412      tmp = VECTOR_ELT(dn, 0);    /* swap the dimnames */      tmp = VECTOR_ELT(dn, 0);    /* swap the dimnames */
413      SET_VECTOR_ELT(dn, 0, VECTOR_ELT(dn, 1));      SET_VECTOR_ELT(dn, 0, VECTOR_ELT(dn, 1));
414      SET_VECTOR_ELT(dn, 1, tmp);      SET_VECTOR_ELT(dn, 1, tmp);
415      if(!isNull(tmp = getAttrib(dn, R_NamesSymbol))) { // swap names(dimnames(.)):      tmp = PROTECT(getAttrib(dn, R_NamesSymbol));
416        if(!isNull(tmp)) { // swap names(dimnames(.)):
417          SEXP nms_dns = PROTECT(allocVector(VECSXP, 2));          SEXP nms_dns = PROTECT(allocVector(VECSXP, 2));
418          SET_VECTOR_ELT(nms_dns, 1, STRING_ELT(tmp, 0));          SET_VECTOR_ELT(nms_dns, 1, STRING_ELT(tmp, 0));
419          SET_VECTOR_ELT(nms_dns, 0, STRING_ELT(tmp, 1));          SET_VECTOR_ELT(nms_dns, 0, STRING_ELT(tmp, 1));
420          setAttrib(dn, R_NamesSymbol, nms_dns);          setAttrib(dn, R_NamesSymbol, nms_dns);
421          UNPROTECT(1);          UNPROTECT(1);
422      }      }
423      UNPROTECT(1);  
424      return chm_sparse_to_SEXP(chxt, 1, /* SWAP 'uplo' for triangular */      SEXP ans = chm_sparse_to_SEXP(chxt, 1, /* SWAP 'uplo' for triangular */
425                                tr ? ((*uplo_P(x) == 'U') ? -1 : 1) : 0,                                tr ? ((*uplo_P(x) == 'U') ? -1 : 1) : 0,
426                                Rkind, tr ? diag_P(x) : "", dn);                                Rkind, tr ? diag_P(x) : "", dn);
427        UNPROTECT(2);
428        return ans;
429  }  }
430    
431  /** @brief  A %*% B  - for matrices of class CsparseMatrix (R package "Matrix")  /** @brief  A %*% B  - for matrices of class CsparseMatrix (R package "Matrix")
# Line 665  Line 673 
673      if(cha->xtype == CHOLMOD_PATTERN) {      if(cha->xtype == CHOLMOD_PATTERN) {
674          /* warning(_("Csparse_dense_prod(): cholmod_sdmult() not yet implemented for pattern./ ngCMatrix" */          /* warning(_("Csparse_dense_prod(): cholmod_sdmult() not yet implemented for pattern./ ngCMatrix" */
675          /*        " --> slightly inefficient coercion")); */          /*        " --> slightly inefficient coercion")); */
   
676          // This *fails* to produce a CHOLMOD_REAL ..          // This *fails* to produce a CHOLMOD_REAL ..
677          // CHM_SP chd = cholmod_l_copy(cha, cha->stype, CHOLMOD_REAL, &c);          // CHM_SP chd = cholmod_l_copy(cha, cha->stype, CHOLMOD_REAL, &cl);
678          // --> use our Matrix-classes  
679            // --> use our Matrix-classes: they work:
680          SEXP da = PROTECT(nz2Csparse(a, x_double)); nprot++;          SEXP da = PROTECT(nz2Csparse(a, x_double)); nprot++;
681          cha = AS_CHM_SP(da);          cha = AS_CHM_SP(da);
682      }      }
# Line 804  Line 812 
812  #define CSPARSE_CAT(_KIND_)                                             \  #define CSPARSE_CAT(_KIND_)                                             \
813      CHM_SP chx = AS_CHM_SP__(x), chy = AS_CHM_SP__(y);                  \      CHM_SP chx = AS_CHM_SP__(x), chy = AS_CHM_SP__(y);                  \
814      R_CheckStack();                                                     \      R_CheckStack();                                                     \
815      int Rk_x = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : -3,     \      int Rk_x = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : x_pattern, \
816          Rk_y = (chy->xtype != CHOLMOD_PATTERN) ? Real_kind(y) : -3, Rkind; \          Rk_y = (chy->xtype != CHOLMOD_PATTERN) ? Real_kind(y) : x_pattern, Rkind; \
817      if(Rk_x == -3 || Rk_y == -3) { /* at least one of them is patter"n" */ \      if(Rk_x == x_pattern || Rk_y == x_pattern) { /* at least one of them is patter"n" */ \
818          if(Rk_x == -3 && Rk_y == -3) { /* fine */                       \          if(Rk_x == x_pattern && Rk_y == x_pattern) { /* fine */         \
819          } else { /* only one is a patter"n"                             \          } else { /* only one is a patter"n"                             \
820                    * "Bug" in cholmod_horzcat()/vertcat(): returns patter"n" matrix if one of them is */ \                    * "Bug" in cholmod_horzcat()/vertcat():               \
821                      * returns patter"n" matrix if one of them is */       \
822              Rboolean ok;                                                \              Rboolean ok;                                                \
823              if(Rk_x == -3) {                                            \              if(Rk_x == x_pattern) {                                     \
824                  ok = chm_MOD_xtype(CHOLMOD_REAL, chx, &c); Rk_x = 0;    \                  ok = chm_MOD_xtype(CHOLMOD_REAL, chx, &c); Rk_x = 0;    \
825              } else if(Rk_y == -3) {                                     \              } else if(Rk_y == x_pattern) {                              \
826                  ok = chm_MOD_xtype(CHOLMOD_REAL, chy, &c); Rk_y = 0;    \                  ok = chm_MOD_xtype(CHOLMOD_REAL, chy, &c); Rk_y = 0;    \
827              } else                                                      \              } else                                                      \
828                  error(_("Impossible Rk_x/Rk_y in Csparse_%s(), please report"), _KIND_); \                  error(_("Impossible Rk_x/Rk_y in Csparse_%s(), please report"), _KIND_); \
# Line 926  Line 935 
935      if (csize >= 0 && !isInteger(j))      if (csize >= 0 && !isInteger(j))
936          error(_("Index j must be NULL or integer"));          error(_("Index j must be NULL or integer"));
937    
938        /* Must treat 'NA's in i[] and j[] here -- they are *not* treated by Cholmod!
939         * haveNA := ...
940           if(haveNA) {
941             a. i = removeNA(i); j =removeNA(j), and remember where they were
942             b. ans = CHM_SUB(.., i, j)
943             c. add NA rows and/or columns to 'ans' according to
944                place of NA's in i and/or j.
945           } else {
946             ans = CHM_SUB(.....)  // == current code
947           }
948         */
949  #define CHM_SUB(_M_, _i_, _j_)                                  \  #define CHM_SUB(_M_, _i_, _j_)                                  \
950      cholmod_submatrix(_M_,                                      \      cholmod_submatrix(_M_,                                      \
951                        (rsize < 0) ? NULL : INTEGER(_i_), rsize, \                        (rsize < 0) ? NULL : INTEGER(_i_), rsize, \
# Line 935  Line 955 
955      if (!chx->stype) {/* non-symmetric Matrix */      if (!chx->stype) {/* non-symmetric Matrix */
956          ans = CHM_SUB(chx, i, j);          ans = CHM_SUB(chx, i, j);
957      }      }
958      else {      else { /* symmetric : "dsCMatrix";
959          /* for now, cholmod_submatrix() only accepts "generalMatrix" */                currently, cholmod_submatrix() only accepts "generalMatrix" */
960          CHM_SP tmp = cholmod_copy(chx, /* stype: */ 0, chx->xtype, &c);          CHM_SP tmp = cholmod_copy(chx, /* stype: */ 0, chx->xtype, &c);
961          ans = CHM_SUB(tmp, i, j);          ans = CHM_SUB(tmp, i, j);
962          cholmod_free_sparse(&tmp, &c);          cholmod_free_sparse(&tmp, &c);
# Line 948  Line 968 
968      /*  dn = PROTECT(allocVector(VECSXP, 2)); */      /*  dn = PROTECT(allocVector(VECSXP, 2)); */
969      return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "", /* dimnames: */ R_NilValue);      return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "", /* dimnames: */ R_NilValue);
970  }  }
971    #undef CHM_SUB
972    
973  #define _d_Csp_  #define _d_Csp_
974  #include "t_Csparse_subassign.c"  #include "t_Csparse_subassign.c"
# Line 1149  Line 1170 
1170                      int index1)                      int index1)
1171  {  {
1172      SEXP ans;      SEXP ans;
1173      int *ij = (int*)NULL, *tri, *trj,      int *ij = (int*)NULL, *tri, *trj, nrow = -1, ncol = -1;
         mi, mj, mp, nrow = -1, ncol = -1;  
1174      int xtype = -1;             /* -Wall */      int xtype = -1;             /* -Wall */
1175      CHM_TR T;      CHM_TR T;
1176      CHM_SP A;      CHM_SP A;
# Line 1158  Line 1178 
1178      if (np < 0 || nnz < 0)      if (np < 0 || nnz < 0)
1179          error(_("negative vector lengths not allowed: np = %d, nnz = %d"),          error(_("negative vector lengths not allowed: np = %d, nnz = %d"),
1180                np, nnz);                np, nnz);
1181      if (1 != ((mi = (i == (int*)NULL)) +      int mi = (i == (int*)NULL), // := missing 'i'
1182                (mj = (j == (int*)NULL)) +          mj = (j == (int*)NULL), // := missing 'j'
1183                (mp = (p == (int*)NULL))))          mp = (p == (int*)NULL); // := missing 'p'
1184        if ((mi + mj + mp) != 1)
1185          error(_("exactly 1 of 'i', 'j' or 'p' must be NULL"));          error(_("exactly 1 of 'i', 'j' or 'p' must be NULL"));
1186      if (mp) {      if (mp) {
1187          if (np) error(_("np = %d, must be zero when p is NULL"), np);          if (np) error(_("np = %d, must be zero when p is NULL"), np);
# Line 1211  Line 1232 
1232      /* check the class name */      /* check the class name */
1233      if (strlen(cls) != 8)      if (strlen(cls) != 8)
1234          error(_("strlen of cls argument = %d, should be 8"), strlen(cls));          error(_("strlen of cls argument = %d, should be 8"), strlen(cls));
1235      if (!strcmp(cls + 2, "CMatrix"))      if (strcmp(cls + 2, "CMatrix"))
1236          error(_("cls = \"%s\" does not end in \"CMatrix\""), cls);          error(_("cls = \"%s\" does not end in \"CMatrix\""), cls);
1237      switch(cls[0]) {      switch(cls[0]) {
1238      case 'd':      case 'd':
# Line 1240  Line 1261 
1261      A = cholmod_triplet_to_sparse(T, nnz, &c);      A = cholmod_triplet_to_sparse(T, nnz, &c);
1262      cholmod_free_triplet(&T, &c);      cholmod_free_triplet(&T, &c);
1263      /* copy the information to the SEXP */      /* copy the information to the SEXP */
1264      ans = PROTECT(NEW_OBJECT(MAKE_CLASS(cls)));      ans = PROTECT(NEW_OBJECT_OF_CLASS(cls));
1265  // FIXME: This has been copied from chm_sparse_to_SEXP in  chm_common.c  // FIXME: This has been copied from chm_sparse_to_SEXP in  chm_common.c
1266      /* allocate and copy common slots */      /* allocate and copy common slots */
1267      nnz = cholmod_nnz(A, &c);      nnz = cholmod_nnz(A, &c);
# Line 1248  Line 1269 
1269      dims[0] = A->nrow; dims[1] = A->ncol;      dims[0] = A->nrow; dims[1] = A->ncol;
1270      Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, A->ncol + 1)), (int*)A->p, A->ncol + 1);      Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, A->ncol + 1)), (int*)A->p, A->ncol + 1);
1271      Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nnz)), (int*)A->i, nnz);      Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nnz)), (int*)A->i, nnz);
1272      switch(cls[1]) {      switch(cls[0]) {
1273      case 'd':      case 'd':
1274          Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nnz)), (double*)A->x, nnz);          Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nnz)), (double*)A->x, nnz);
1275          break;          break;
# Line 1260  Line 1281 
1281      UNPROTECT(1);      UNPROTECT(1);
1282      return ans;      return ans;
1283  }  }
1284    
1285    /**
1286     * Create a Csparse matrix object from a traditional R matrix
1287     *
1288     * @param x   traditional R matrix (numeric, logical, ...)
1289     * @param cls class (a string)
1290     *
1291     * @return an SEXP of a class inheriting from CsparseMatrix.
1292     */
1293    SEXP matrix_to_Csparse(SEXP x, SEXP cls)
1294    {
1295        if (!isMatrix(x))
1296            error(_("%s must be (traditional R) matrix"), "'x'");
1297        SEXP d_x  = getAttrib(x, R_DimSymbol),
1298            dn_x  = getAttrib(x, R_DimNamesSymbol);
1299        int nr = INTEGER(d_x)[0],
1300            nc = INTEGER(d_x)[1];
1301    
1302        if (!(isString(cls) && LENGTH(cls) == 1))
1303            error(_("%s must be character string"), "'cls'");
1304        R_xlen_t ii, n = XLENGTH(x);
1305        int xtype = -1;
1306        if (n != ((R_xlen_t) nr) * nc)
1307            error(_("nrow * ncol = %d * %d must equal length(x) = %ld"), nr, nc, n);
1308    
1309        const char *ccls = CHAR(STRING_ELT(cls, 0));
1310        if (strlen(ccls) != 9)
1311            error(_("strlen of cls argument = %d, should be 9"), strlen(ccls));
1312        if (strcmp(ccls + 2, "CMatrix"))
1313            error(_("cls = \"%s\" does not end in \"CMatrix\""), ccls);
1314        switch(ccls[0]) {
1315        case 'd':
1316        case 'l':
1317            xtype = CHOLMOD_REAL;
1318        break;
1319        case 'n':
1320            xtype = CHOLMOD_PATTERN;
1321            break;
1322        default:
1323            error(_("cls = \"%s\" must begin with 'd', 'l' or 'n' for now"), ccls);
1324        }
1325        /* if (ccls[1] != 'g') */
1326        /*  error(_("Only 'g'eneral sparse matrix types allowed")); */
1327    
1328        SEXP ans = PROTECT(NEW_OBJECT_OF_CLASS(ccls));
1329        SET_SLOT(ans, Matrix_DimSym, d_x);
1330        SET_SLOT(ans, Matrix_DimNamesSym, (!isNull(dn_x) && LENGTH(dn_x) == 2)
1331                 ? duplicate(dn_x)
1332                 : allocVector(VECSXP, 2));
1333    
1334        int nz = 0, // current number of nonzero entries
1335            nnz = imax2(256, imax2(nr,nc));/* nnz := final number of nonzero entries, yet unknown;
1336                                               -- must start with guess and then grow */
1337        int *rp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, nc + 1)),
1338            *ri = Calloc(nnz, int); // to become i slot -- of not-yet-known length nnz
1339    
1340        rp[0] = 0; // always
1341    
1342        switch(TYPEOF(x)) {
1343        case LGLSXP: {
1344            if(xtype == CHOLMOD_PATTERN) {
1345    #         define _PATTERN_x
1346    #         include "t_matrix_to_Csp.c"
1347            } else {
1348    #         define _LOGICAL_x
1349    #         include "t_matrix_to_Csp.c"
1350            }
1351            break;
1352        }
1353        case REALSXP: {
1354    #       define _DOUBLE_x
1355    #       include "t_matrix_to_Csp.c"
1356            break;
1357        }
1358    /* case INTSXP: we would have to use
1359            x = coerceVector(x, REALSXP));
1360       and then fall through to REALSXP case, but we must *not* modify 'x' here
1361       FIXME: use a macro or (inline?) function with argument (y), where
1362       -----  SEXP y = PROTECT(coerceVector(x, REALSXP))
1363    
1364       ==> give error in INTSXP case, so caller (in R) must set  storage.mode(x) <- "double"
1365    */
1366    #ifdef _USING_INTEGER_NOT_READY__
1367        case INTSXP: {
1368    #       define _INTEGER_x
1369    #       include "t_matrix_to_Csp.c"
1370            break;
1371        }
1372    #endif
1373    #ifdef _USING_COMPLEX_NOT_READY__
1374        case CPLXSXP: {
1375    #       define _COMPLEX_x
1376    #       include "t_matrix_to_Csp.c"
1377            break;
1378        }
1379    #endif
1380        default:
1381            error(_("%s must be a logical or double vector"), "'x'");
1382            break;
1383        }
1384    
1385        Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym,  INTSXP, nnz)), ri, nnz);
1386        Free(ri);
1387    
1388        UNPROTECT(1);
1389        return ans;
1390    }
1391    
1392    
1393    
1394    

Legend:
Removed from v.3147  
changed lines
  Added in v.3270

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business University of Wisconsin - Madison Powered By FusionForge