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 3257, Sat Mar 17 06:10:05 2018 UTC revision 3258, Sat Mar 17 06:18:00 2018 UTC
# Line 349  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 395  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 674  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 813  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 1170  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 1179  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 1232  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 1269  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 1281  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(MAKE_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.3257  
changed lines
  Added in v.3258

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