SCM

SCM Repository

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

Diff of /pkg/src/chm_common.c

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

revision 2298, Fri Oct 17 13:45:12 2008 UTC revision 2299, Fri Oct 17 16:07:38 2008 UTC
# Line 4  Line 4 
4  Rboolean isValid_Csparse(SEXP x); /* -> Csparse.c */  Rboolean isValid_Csparse(SEXP x); /* -> Csparse.c */
5    
6  cholmod_common c;  cholmod_common c;
7  cholmod_common cl;  /* cholmod_common cl; */
8    
9  static int stype(int ctype, SEXP x)  static int stype(int ctype, SEXP x)
10  {  {
# Line 64  Line 64 
64    
65      /* R_alloc the vector storage for dest and copy the contents from src */      /* R_alloc the vector storage for dest and copy the contents from src */
66      np1 = src->ncol + 1;      np1 = src->ncol + 1;
67      nnz = (int) cholmod_nnz(src, &c);      nnz = (int) cholmod_l_nnz(src, &c);
68      dest->p = (void*) Memcpy((   int*)R_alloc(sizeof(   int), np1),      dest->p = (void*) Memcpy((   int*)R_alloc(sizeof(   int), np1),
69                               (   int*)(src->p), np1);                               (   int*)(src->p), np1);
70      dest->i = (void*) Memcpy((   int*)R_alloc(sizeof(   int), nnz),      dest->i = (void*) Memcpy((   int*)R_alloc(sizeof(   int), nnz),
# Line 107  Line 107 
107          error("invalid object passed to as_cholmod_sparse");          error("invalid object passed to as_cholmod_sparse");
108      memset(ans, 0, sizeof(cholmod_sparse)); /* zero the struct */      memset(ans, 0, sizeof(cholmod_sparse)); /* zero the struct */
109    
110      ans->itype = CHOLMOD_INT;   /* characteristics of the system */      ans->itype = CHOLMOD_LONG;  /* characteristics of the system */
111      ans->dtype = CHOLMOD_DOUBLE;      ans->dtype = CHOLMOD_DOUBLE;
112      ans->packed = TRUE;      ans->packed = TRUE;
113                                  /* slots always present */                                  /* slots always present */
# Line 128  Line 128 
128      ans->sorted = check_sorted_chm(ans);      ans->sorted = check_sorted_chm(ans);
129      if (!(ans->sorted)) { /* sort columns */      if (!(ans->sorted)) { /* sort columns */
130          if(sort_in_place) {          if(sort_in_place) {
131              if (!cholmod_sort(ans, &c))              if (!cholmod_l_sort(ans, &c))
132                  error(_("in_place cholmod_sort returned an error code"));                  error(_("in_place cholmod_l_sort returned an error code"));
133              ans->sorted = 1;              ans->sorted = 1;
134          }          }
135          else {          else {
136              CHM_SP tmp = cholmod_copy_sparse(ans, &c);              CHM_SP tmp = cholmod_l_copy_sparse(ans, &c);
137              if (!cholmod_sort(tmp, &c))              if (!cholmod_l_sort(tmp, &c))
138                  error(_("cholmod_sort returned an error code"));                  error(_("cholmod_l_sort returned an error code"));
139    
140  #ifdef DEBUG_Matrix  #ifdef DEBUG_Matrix
141              /* This "triggers" exactly for return values of dtCMatrix_sparse_solve():*/              /* This "triggers" exactly for return values of dtCMatrix_sparse_solve():*/
142              /* Don't want to translate this: want it report */              /* Don't want to translate this: want it report */
143              Rprintf("Note: as_cholmod_sparse() needed cholmod_sort()ing\n");              Rprintf("Note: as_cholmod_l_sparse() needed cholmod_l_sort()ing\n");
144  #endif  #endif
145              chm2Ralloc(ans, tmp);              chm2Ralloc(ans, tmp);
146              cholmod_free_sparse(&tmp, &c);              cholmod_l_free_sparse(&tmp, &c);
147          }          }
148      }      }
149    
150      if (check_Udiag && ctype % 3 == 2 && (*diag_P(x) == 'U')) { /* diagU2N(.)  "in place" : */      if (check_Udiag && ctype % 3 == 2 && (*diag_P(x) == 'U')) { /* diagU2N(.)  "in place" : */
151          double one[] = {1, 0};          double one[] = {1, 0};
152          CHM_SP eye = cholmod_speye(ans->nrow, ans->ncol, ans->xtype, &c);          CHM_SP eye = cholmod_l_speye(ans->nrow, ans->ncol, ans->xtype, &c);
153          CHM_SP tmp = cholmod_add(ans, eye, one, one, TRUE, TRUE, &c);          CHM_SP tmp = cholmod_l_add(ans, eye, one, one, TRUE, TRUE, &c);
154    
155          chm2Ralloc(ans, tmp);          chm2Ralloc(ans, tmp);
156          cholmod_free_sparse(&tmp, &c);          cholmod_l_free_sparse(&tmp, &c);
157          cholmod_free_sparse(&eye, &c);          cholmod_l_free_sparse(&eye, &c);
158      } /* else :      } /* else :
159         * NOTE: if(*diag_P(x) == 'U'), the diagonal is lost (!);         * NOTE: if(*diag_P(x) == 'U'), the diagonal is lost (!);
160         * ---- that may be ok, e.g. if we are just converting from/to Tsparse,         * ---- that may be ok, e.g. if we are just converting from/to Tsparse,
# Line 168  Line 168 
168   * optionally, free a or free both a and its the pointers to its contents.   * optionally, free a or free both a and its the pointers to its contents.
169   *   *
170   * @param a matrix to be converted   * @param a matrix to be converted
171   * @param dofree 0 - don't free a; > 0 cholmod_free a; < 0 Free a   * @param dofree 0 - don't free a; > 0 cholmod_l_free a; < 0 Free a
172   * @param uploT 0 - not triangular; > 0 upper triangular; < 0 lower   * @param uploT 0 - not triangular; > 0 upper triangular; < 0 lower
173   * @param Rkind - vector type to store for a->xtype == CHOLMOD_REAL,   * @param Rkind - vector type to store for a->xtype == CHOLMOD_REAL,
174   *                0 - REAL; 1 - LOGICAL   *                0 - REAL; 1 - LOGICAL
# Line 183  Line 183 
183  {  {
184      SEXP ans;      SEXP ans;
185      char *cls = "";/* -Wall */      char *cls = "";/* -Wall */
186      int *dims, nnz, *ansp, *ansi, *aii = (int*)(a->i), *api = (int*)(a->p),      int *dims, nnz, *aii = (int*)(a->i), *api = (int*)(a->p);
         longi = (a->itype) == CHOLMOD_LONG;  
     UF_long *ail = (UF_long*)(a->i), *apl = (UF_long*)(a->p);  
187    
188      PROTECT(dn);  /* dn is usually UNPROTECTed before the call */      PROTECT(dn);  /* dn is usually UNPROTECTed before the call */
189    
190                                  /* ensure a is sorted and packed */                                  /* ensure a is sorted and packed */
191      if (!a->sorted || !a->packed)      if (!a->sorted || !a->packed) cholmod_l_sort(a, &c);
         longi ? cholmod_l_sort(a, &cl) : cholmod_sort(a, &c);  
192                                  /* determine the class of the result */                                  /* determine the class of the result */
193      switch(a->xtype){      switch(a->xtype){
194      case CHOLMOD_PATTERN:      case CHOLMOD_PATTERN:
# Line 214  Line 211 
211      }      }
212      ans = PROTECT(NEW_OBJECT(MAKE_CLASS(cls)));      ans = PROTECT(NEW_OBJECT(MAKE_CLASS(cls)));
213                                  /* allocate and copy common slots */                                  /* allocate and copy common slots */
214      nnz = longi ? cholmod_l_nnz(a, &cl) : cholmod_nnz(a, &c);      nnz = cholmod_l_nnz(a, &c);
215      dims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));      dims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));
216      dims[0] = a->nrow; dims[1] = a->ncol;      dims[0] = a->nrow; dims[1] = a->ncol;
217      ansp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, a->ncol + 1));      Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, a->ncol + 1)), api, a->ncol + 1);
218      ansi = INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nnz));      Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nnz)), aii, nnz);
     for (int j = 0; j <= a->ncol; j++) ansp[j] = longi ? (int)(apl[j]) : api[j];  
     for (int p = 0; p < nnz; p++) ansi[p] = longi ? (int)(ail[p]) : aii[p];  
219                                  /* copy data slot if present */                                  /* copy data slot if present */
220      if (a->xtype == CHOLMOD_REAL) {      if (a->xtype == CHOLMOD_REAL) {
221          int i, *m_x;          int i, *m_x;
# Line 248  Line 243 
243      if (a->stype)               /* slot for symmetricMatrix */      if (a->stype)               /* slot for symmetricMatrix */
244          SET_SLOT(ans, Matrix_uploSym,          SET_SLOT(ans, Matrix_uploSym,
245                   mkString((a->stype > 0) ? "U" : "L"));                   mkString((a->stype > 0) ? "U" : "L"));
246      if (dofree > 0)      if (dofree > 0) cholmod_l_free_sparse(&a, &c);
         longi ? cholmod_l_free_sparse(&a, &cl) : cholmod_free_sparse(&a, &c);  
247      if (dofree < 0) Free(a);      if (dofree < 0) Free(a);
248      if (dn != R_NilValue)      if (dn != R_NilValue)
249          SET_SLOT(ans, Matrix_DimNamesSym, duplicate(dn));          SET_SLOT(ans, Matrix_DimNamesSym, duplicate(dn));
# Line 287  Line 281 
281      if (ctype < 0) error("invalid class of object to as_cholmod_triplet");      if (ctype < 0) error("invalid class of object to as_cholmod_triplet");
282      memset(ans, 0, sizeof(cholmod_triplet)); /* zero the struct */      memset(ans, 0, sizeof(cholmod_triplet)); /* zero the struct */
283    
284      ans->itype = CHOLMOD_INT;   /* characteristics of the system */      ans->itype = CHOLMOD_LONG;  /* characteristics of the system */
285      ans->dtype = CHOLMOD_DOUBLE;      ans->dtype = CHOLMOD_DOUBLE;
286      ans->x = ans->z = (void *) NULL;      ans->x = ans->z = (void *) NULL;
287                                  /* dimensions and nzmax */                                  /* dimensions and nzmax */
# Line 313  Line 307 
307          /* TODO? instead of reallocating, don't do the 2nd part of AS_CHM_COMMON()          /* TODO? instead of reallocating, don't do the 2nd part of AS_CHM_COMMON()
308           * ---- above, and allocate to correct length + Memcpy() here, as in           * ---- above, and allocate to correct length + Memcpy() here, as in
309           * Tsparse_diagU2N() */           * Tsparse_diagU2N() */
310          if(cholmod_reallocate_triplet((size_t) k, ans, &c))          if(cholmod_l_reallocate_triplet((size_t) k, ans, &c))
311              error(_("as_cholmod_triplet(): could not reallocate for internal diagU2N()"              error(_("as_cholmod_l_triplet(): could not reallocate for internal diagU2N()"
312                        ));                        ));
313          a_i = ans->i;          a_i = ans->i;
314          a_j = ans->j;          a_j = ans->j;
# Line 358  Line 352 
352   * optionally, free a or free both a and its the pointers to its contents.   * optionally, free a or free both a and its the pointers to its contents.
353   *   *
354   * @param a matrix to be converted   * @param a matrix to be converted
355   * @param dofree 0 - don't free a; > 0 cholmod_free a; < 0 Free a   * @param dofree 0 - don't free a; > 0 cholmod_l_free a; < 0 Free a
356   * @param uploT 0 - not triangular; > 0 upper triangular; < 0 lower   * @param uploT 0 - not triangular; > 0 upper triangular; < 0 lower
357   * @param Rkind - vector type to store for a->xtype == CHOLMOD_REAL,   * @param Rkind - vector type to store for a->xtype == CHOLMOD_REAL,
358   *                0 - REAL; 1 - LOGICAL   *                0 - REAL; 1 - LOGICAL
# Line 438  Line 432 
432      if (a->stype)      if (a->stype)
433          SET_SLOT(ans, Matrix_uploSym,          SET_SLOT(ans, Matrix_uploSym,
434                   mkString((a->stype > 0) ? "U" : "L"));                   mkString((a->stype > 0) ? "U" : "L"));
435      if (dofree > 0) cholmod_free_triplet(&a, &c);      if (dofree > 0) cholmod_l_free_triplet(&a, &c);
436      if (dofree < 0) Free(a);      if (dofree < 0) Free(a);
437      if (dn != R_NilValue)      if (dn != R_NilValue)
438          SET_SLOT(ans, Matrix_DimNamesSym, duplicate(dn));          SET_SLOT(ans, Matrix_DimNamesSym, duplicate(dn));
# Line 566  Line 560 
560   *   *
561   * @return CHOLMOD_OK if successful   * @return CHOLMOD_OK if successful
562   */   */
 int R_cholmod_start(CHM_CM c)  
 {  
     int res;  
     if (!(res = cholmod_start(c)))  
         error(_("Unable to initialize cholmod: error code %d"), res);  
     c->print_function = R_cholmod_printf; /* Rprintf gives warning */  
     /* Since we provide an error handler, it may not be a good idea to allow CHOLMOD printing,  
      * because that's not easily suppressed on the R level :  
      * Hence consider, at least temporarily *  c->print_function = NULL; */  
     c->error_handler = R_cholmod_error;  
     return TRUE;  
 }  
   
 /**  
  * Initialize the CHOLMOD library and replace the print and error functions  
  * by R-specific versions.  
  *  
  * @param c pointer to a cholmod_common structure to be initialized  
  *  
  * @return CHOLMOD_OK if successful  
  */  
563  int R_cholmod_l_start(CHM_CM cl)  int R_cholmod_l_start(CHM_CM cl)
564  {  {
565      int res;      int res;
# Line 605  Line 578 
578   * optionally, free a or free both a and its pointer to its contents.   * optionally, free a or free both a and its pointer to its contents.
579   *   *
580   * @param a matrix to be converted   * @param a matrix to be converted
581   * @param dofree 0 - don't free a; > 0 cholmod_free a; < 0 Free a   * @param dofree 0 - don't free a; > 0 cholmod_l_free a; < 0 Free a
582   * @param Rkind type of R matrix to be generated (special to this function)   * @param Rkind type of R matrix to be generated (special to this function)
583   * @param dn   -- dimnames [list(.,.) or NULL]   * @param dn   -- dimnames [list(.,.) or NULL]
584   *   *
# Line 670  Line 643 
643  /*             (complex *) a->x, ntot); */  /*             (complex *) a->x, ntot); */
644      } else error("code for cholmod_dense with holes not yet written");      } else error("code for cholmod_dense with holes not yet written");
645    
646      if (dofree > 0) cholmod_free_dense(&a, &c);      if (dofree > 0) cholmod_l_free_dense(&a, &c);
647      if (dofree < 0) Free(a);      if (dofree < 0) Free(a);
648      if (dn != R_NilValue)      if (dn != R_NilValue)
649          SET_SLOT(ans, Matrix_DimNamesSym, duplicate(dn));          SET_SLOT(ans, Matrix_DimNamesSym, duplicate(dn));
# Line 683  Line 656 
656   * or free both a and its pointer to its contents.   * or free both a and its pointer to its contents.
657   *   *
658   * @param a cholmod_dense structure to be converted   * @param a cholmod_dense structure to be converted
659   * @param dofree 0 - don't free a; > 0 cholmod_free a; < 0 Free a   * @param dofree 0 - don't free a; > 0 cholmod_l_free a; < 0 Free a
660   * @param dn either R_NilValue or an SEXP suitable for the Dimnames slot.   * @param dn either R_NilValue or an SEXP suitable for the Dimnames slot.
661   *   *
662   * @return SEXP containing a copy of a as a matrix object   * @return SEXP containing a copy of a as a matrix object
# Line 712  Line 685 
685  /*             (complex *) a->x, a->nz); */  /*             (complex *) a->x, a->nz); */
686      } else error("code for cholmod_dense with holes not yet written");      } else error("code for cholmod_dense with holes not yet written");
687    
688      if (dofree > 0) cholmod_free_dense(&a, &c);      if (dofree > 0) cholmod_l_free_dense(&a, &c);
689      if (dofree < 0) Free(a);      if (dofree < 0) Free(a);
690      if (dn != R_NilValue)      if (dn != R_NilValue)
691          setAttrib(ans, R_DimNamesSymbol, duplicate(dn));          setAttrib(ans, R_DimNamesSymbol, duplicate(dn));
# Line 754  Line 727 
727      if (ctype < 0) error("invalid class of object to as_cholmod_factor");      if (ctype < 0) error("invalid class of object to as_cholmod_factor");
728      memset(ans, 0, sizeof(cholmod_factor)); /* zero the struct */      memset(ans, 0, sizeof(cholmod_factor)); /* zero the struct */
729    
730      ans->itype = CHOLMOD_INT;   /* characteristics of the system */      ans->itype = CHOLMOD_LONG;  /* characteristics of the system */
731      ans->dtype = CHOLMOD_DOUBLE;      ans->dtype = CHOLMOD_DOUBLE;
732      ans->z = (void *) NULL;      ans->z = (void *) NULL;
733      ans->xtype = (ctype < 2) ? CHOLMOD_REAL : CHOLMOD_PATTERN;      ans->xtype = (ctype < 2) ? CHOLMOD_REAL : CHOLMOD_PATTERN;
# Line 804  Line 777 
777          ans->next = INTEGER(GET_SLOT(x, install("nxt")));          ans->next = INTEGER(GET_SLOT(x, install("nxt")));
778          ans->prev = INTEGER(GET_SLOT(x, install("prv")));          ans->prev = INTEGER(GET_SLOT(x, install("prv")));
779      }      }
780      if (!cholmod_check_factor(ans, &c))      if (!cholmod_l_check_factor(ans, &c))
781          error(_("failure in as_cholmod_factor"));          error(_("failure in as_cholmod_factor"));
782      return ans;      return ans;
783  }  }
# Line 875  Line 848 
848    
849      }      }
850      if(dofree) {      if(dofree) {
851          if (dofree > 0) cholmod_free_factor(&f, &c);          if (dofree > 0) cholmod_l_free_factor(&f, &c);
852          else /* dofree < 0 */ Free(f);          else /* dofree < 0 */ Free(f);
853      }      }
854      UNPROTECT(1);      UNPROTECT(1);
# Line 895  Line 868 
868   */   */
869  void chm_diagN2U(CHM_SP chx, int uploT, Rboolean do_realloc)  void chm_diagN2U(CHM_SP chx, int uploT, Rboolean do_realloc)
870  {  {
871      int i, n = chx->nrow, nnz = (int)cholmod_nnz(chx, &c),      int i, n = chx->nrow, nnz = (int)cholmod_l_nnz(chx, &c),
872          n_nnz = nnz - n, /* the new nnz : we will have removed  n entries */          n_nnz = nnz - n, /* the new nnz : we will have removed  n entries */
873          i_to = 0, i_from = 0;          i_to = 0, i_from = 0;
874    
# Line 903  Line 876 
876          error(_("chm_diagN2U(<non-square matrix>): nrow=%d, ncol=%d"),          error(_("chm_diagN2U(<non-square matrix>): nrow=%d, ncol=%d"),
877                n, chx->ncol);                n, chx->ncol);
878    
879      if (!chx->sorted || !chx->packed) cholmod_sort(chx, &c);      if (!chx->sorted || !chx->packed) cholmod_l_sort(chx, &c);
880                                  /* dimensions and nzmax */                                  /* dimensions and nzmax */
881    
882  #define _i(I) (   (int*) chx->i)[I]  #define _i(I) (   (int*) chx->i)[I]
# Line 955  Line 928 
928  #undef _p  #undef _p
929    
930      if(do_realloc) /* shorten (i- and x-slots from nnz to n_nnz */      if(do_realloc) /* shorten (i- and x-slots from nnz to n_nnz */
931          cholmod_reallocate_sparse(n_nnz, chx, &c);          cholmod_l_reallocate_sparse(n_nnz, chx, &c);
932      return;      return;
933  }  }
934    

Legend:
Removed from v.2298  
changed lines
  Added in v.2299

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