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 2685, Fri Aug 5 19:52:10 2011 UTC revision 2889, Thu Aug 8 21:06:22 2013 UTC
# Line 37  Line 37 
37      return Csparse_validate_(x, FALSE);      return Csparse_validate_(x, FALSE);
38  }  }
39    
40    
41    #define _t_Csparse_validate
42    #include "t_Csparse_validate.c"
43    
44    #define _t_Csparse_sort
45    #include "t_Csparse_validate.c"
46    
47    // R: .validateCsparse(x, sort.if.needed = FALSE) :
48  SEXP Csparse_validate2(SEXP x, SEXP maybe_modify) {  SEXP Csparse_validate2(SEXP x, SEXP maybe_modify) {
49      return Csparse_validate_(x, asLogical(maybe_modify));      return Csparse_validate_(x, asLogical(maybe_modify));
50  }  }
51    
52  SEXP Csparse_validate_(SEXP x, Rboolean maybe_modify)  // R: Matrix:::.sortCsparse(x) :
53  {  SEXP Csparse_sort (SEXP x) {
54      /* NB: we do *NOT* check a potential 'x' slot here, at all */     int ok = Csparse_sort_2(x, TRUE); // modifying x directly
55      SEXP pslot = GET_SLOT(x, Matrix_pSym),     if(!ok) warning(_("Csparse_sort(x): x is not a valid (apart from sorting) CsparseMatrix"));
56          islot = GET_SLOT(x, Matrix_iSym);     return x;
     Rboolean sorted, strictly;  
     int j, k,  
         *dims = INTEGER(GET_SLOT(x, Matrix_DimSym)),  
         nrow = dims[0],  
         ncol = dims[1],  
         *xp = INTEGER(pslot),  
         *xi = INTEGER(islot);  
   
     if (length(pslot) != dims[1] + 1)  
         return mkString(_("slot p must have length = ncol(.) + 1"));  
     if (xp[0] != 0)  
         return mkString(_("first element of slot p must be zero"));  
     if (length(islot) < xp[ncol]) /* allow larger slots from over-allocation!*/  
         return  
             mkString(_("last element of slot p must match length of slots i and x"));  
     for (j = 0; j < xp[ncol]; j++) {  
         if (xi[j] < 0 || xi[j] >= nrow)  
             return mkString(_("all row indices must be between 0 and nrow-1"));  
     }  
     sorted = TRUE; strictly = TRUE;  
     for (j = 0; j < ncol; j++) {  
         if (xp[j] > xp[j + 1])  
             return mkString(_("slot p must be non-decreasing"));  
         if(sorted) /* only act if >= 2 entries in column j : */  
             for (k = xp[j] + 1; k < xp[j + 1]; k++) {  
                 if (xi[k] < xi[k - 1])  
                     sorted = FALSE;  
                 else if (xi[k] == xi[k - 1])  
                     strictly = FALSE;  
             }  
     }  
     if (!sorted) {  
         if(maybe_modify) {  
             CHM_SP chx = (CHM_SP) alloca(sizeof(cholmod_sparse));  
             R_CheckStack();  
             as_cholmod_sparse(chx, x, FALSE, TRUE);/*-> cholmod_l_sort() ! */  
             /* as chx = AS_CHM_SP__(x)  but  ^^^^ sorting x in_place !!! */  
   
             /* Now re-check that row indices are *strictly* increasing  
              * (and not just increasing) within each column : */  
             for (j = 0; j < ncol; j++) {  
                 for (k = xp[j] + 1; k < xp[j + 1]; k++)  
                     if (xi[k] == xi[k - 1])  
                         return mkString(_("slot i is not *strictly* increasing inside a column (even after cholmod_l_sort)"));  
             }  
         } else { /* no modifying sorting : */  
             return mkString(_("row indices are not sorted within columns"));  
         }  
     } else if(!strictly) {  /* sorted, but not strictly */  
         return mkString(_("slot i is not *strictly* increasing inside a column"));  
     }  
     return ScalarLogical(1);  
57  }  }
58    
59  SEXP Rsparse_validate(SEXP x)  SEXP Rsparse_validate(SEXP x)
# Line 189  Line 145 
145      if(cl_x[2] != 'C') error(_("not a CsparseMatrix"));      if(cl_x[2] != 'C') error(_("not a CsparseMatrix"));
146      int nnz = LENGTH(GET_SLOT(x, Matrix_iSym));      int nnz = LENGTH(GET_SLOT(x, Matrix_iSym));
147      SEXP ans;      SEXP ans;
148      char *ncl = strdup(cl_x);      char *ncl = alloca(strlen(cl_x) + 1); /* not much memory required */
149        strcpy(ncl, cl_x);
150      double *dx_x; int *ix_x;      double *dx_x; int *ix_x;
151      ncl[0] = (r_kind == x_double ? 'd' :      ncl[0] = (r_kind == x_double ? 'd' :
152                (r_kind == x_logical ? 'l' :                (r_kind == x_logical ? 'l' :
# Line 266  Line 223 
223    
224  SEXP Csparse_general_to_symmetric(SEXP x, SEXP uplo)  SEXP Csparse_general_to_symmetric(SEXP x, SEXP uplo)
225  {  {
226        int *adims = INTEGER(GET_SLOT(x, Matrix_DimSym)), n = adims[0];
227        if(n != adims[1]) {
228            error(_("Csparse_general_to_symmetric(): matrix is not square!"));
229            return R_NilValue; /* -Wall */
230        }
231      CHM_SP chx = AS_CHM_SP__(x), chgx;      CHM_SP chx = AS_CHM_SP__(x), chgx;
232      int uploT = (*CHAR(STRING_ELT(uplo,0)) == 'U') ? 1 : -1;      int uploT = (*CHAR(STRING_ELT(uplo,0)) == 'U') ? 1 : -1;
233      int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;      int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
234      R_CheckStack();      R_CheckStack();
   
235      chgx = cholmod_copy(chx, /* stype: */ uploT, chx->xtype, &c);      chgx = cholmod_copy(chx, /* stype: */ uploT, chx->xtype, &c);
236      /* xtype: pattern, "real", complex or .. */      /* xtype: pattern, "real", complex or .. */
237      return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "",      return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "",
# Line 566  Line 527 
527      }      }
528      else { /* triangular with diag='N'): now drop the diagonal */      else { /* triangular with diag='N'): now drop the diagonal */
529          /* duplicate, since chx will be modified: */          /* duplicate, since chx will be modified: */
530          CHM_SP chx = AS_CHM_SP__(duplicate(x));          SEXP xx = PROTECT(duplicate(x));
531            CHM_SP chx = AS_CHM_SP__(xx);
532          int uploT = (*uplo_P(x) == 'U') ? 1 : -1,          int uploT = (*uplo_P(x) == 'U') ? 1 : -1,
533              Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;              Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
534          R_CheckStack();          R_CheckStack();
535    
536          chm_diagN2U(chx, uploT, /* do_realloc */ FALSE);          chm_diagN2U(chx, uploT, /* do_realloc */ FALSE);
537    
538          return chm_sparse_to_SEXP(chx, /*dofree*/ 0/* or 1 ?? */,          SEXP ans = chm_sparse_to_SEXP(chx, /*dofree*/ 0/* or 1 ?? */,
539                                    uploT, Rkind, "U",                                    uploT, Rkind, "U",
540                                    GET_SLOT(x, Matrix_DimNamesSym));                                    GET_SLOT(x, Matrix_DimNamesSym));
541            UNPROTECT(1);// only now !
542            return ans;
543      }      }
544  }  }
545    
# Line 601  Line 565 
565      if (csize >= 0 && !isInteger(j))      if (csize >= 0 && !isInteger(j))
566          error(_("Index j must be NULL or integer"));          error(_("Index j must be NULL or integer"));
567    
568      if (chx->stype) /* symmetricMatrix */      if (!chx->stype) {/* non-symmetric Matrix */
         /* for now, cholmod_submatrix() only accepts "generalMatrix" */  
         chx = cholmod_copy(chx, /* stype: */ 0, chx->xtype, &c);  
   
569      return chm_sparse_to_SEXP(cholmod_submatrix(chx,      return chm_sparse_to_SEXP(cholmod_submatrix(chx,
570                                  (rsize < 0) ? NULL : INTEGER(i), rsize,                                  (rsize < 0) ? NULL : INTEGER(i), rsize,
571                                  (csize < 0) ? NULL : INTEGER(j), csize,                                  (csize < 0) ? NULL : INTEGER(j), csize,
# Line 612  Line 573 
573                                1, 0, Rkind, "",                                1, 0, Rkind, "",
574                                /* FIXME: drops dimnames */ R_NilValue);                                /* FIXME: drops dimnames */ R_NilValue);
575  }  }
576                                    /* for now, cholmod_submatrix() only accepts "generalMatrix" */
577        CHM_SP tmp = cholmod_copy(chx, /* stype: */ 0, chx->xtype, &c);
578        CHM_SP ans = cholmod_submatrix(tmp,
579                                       (rsize < 0) ? NULL : INTEGER(i), rsize,
580                                       (csize < 0) ? NULL : INTEGER(j), csize,
581                                       TRUE, TRUE, &c);
582        cholmod_free_sparse(&tmp, &c);
583        return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "", R_NilValue);
584    }
585    
586  #define _d_Csp_  #define _d_Csp_
587  #include "t_Csparse_subassign.c"  #include "t_Csparse_subassign.c"
# Line 712  Line 682 
682      case diag_backpermuted:      case diag_backpermuted:
683          for_DIAG(v[i] = x_x[i_from]);          for_DIAG(v[i] = x_x[i_from]);
684    
685          warning(_("resultKind = 'diagBack' (back-permuted) is experimental"));          warning(_("%s = '%s' (back-permuted) is experimental"),
686                    "resultKind", "diagBack");
687          /* now back_permute : */          /* now back_permute : */
688          for(i = 0; i < n; i++) {          for(i = 0; i < n; i++) {
689              double tmp = v[i]; v[i] = v[perm[i]]; v[perm[i]] = tmp;              double tmp = v[i]; v[i] = v[perm[i]]; v[perm[i]] = tmp;

Legend:
Removed from v.2685  
changed lines
  Added in v.2889

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