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 1547, Mon Sep 11 14:49:39 2006 UTC revision 1548, Mon Sep 11 22:13:07 2006 UTC
# Line 18  Line 18 
18      cholmod_sparse *ans = Calloc(1, cholmod_sparse);      cholmod_sparse *ans = Calloc(1, cholmod_sparse);
19      char *valid[] = {"dgCMatrix", "dsCMatrix", "dtCMatrix",      char *valid[] = {"dgCMatrix", "dsCMatrix", "dtCMatrix",
20                       "lgCMatrix", "lsCMatrix", "ltCMatrix",                       "lgCMatrix", "lsCMatrix", "ltCMatrix",
21                         "ngCMatrix", "nsCMatrix", "ntCMatrix",
22                       "zgCMatrix", "zsCMatrix", "ztCMatrix",                       "zgCMatrix", "zsCMatrix", "ztCMatrix",
23                       ""};                       ""};
24      int *dims, ctype = Matrix_check_class(class_P(x), valid);      int *dims, ctype = Matrix_check_class(class_P(x), valid);
# Line 39  Line 40 
40                                  /* slots always present */                                  /* slots always present */
41      ans->i = (void *) INTEGER(islot);      ans->i = (void *) INTEGER(islot);
42      ans->p = (void *) INTEGER(GET_SLOT(x, Matrix_pSym));      ans->p = (void *) INTEGER(GET_SLOT(x, Matrix_pSym));
                                 /* set the xtype and any elements */  
     switch(ctype) {  
     case 0:  
     case 1:  
     case 2:  
         ans->xtype = CHOLMOD_REAL;  
         ans->x = (void *) REAL(GET_SLOT(x, Matrix_xSym));  
         break;  
     case 3:  
     case 4:  
     case 5:  
         ans->xtype = CHOLMOD_PATTERN;  
         break;  
     case 6:  
     case 7:  
     case 8:  
         ans->xtype = CHOLMOD_COMPLEX;  
         ans->x = (void *) COMPLEX(GET_SLOT(x, Matrix_xSym));  
         break;  
     }  
                                 /* set the stype */  
     switch(ctype % 3) {  
     case 0: ans->stype = 0; break;  
     case 1:  
         ans->stype =  
             (!strcmp(CHAR(asChar(getAttrib(x, Matrix_uploSym))), "U")) ?  
             1 : -1;  
         break;  
     case 2: ans->stype = 0; break; /* Note that triangularity property is lost */  
 /*      error("triangular matrices not yet mapped to CHOLMOD"); */  
     }  
43    
44      return ans;  #define AS_CHM_FINISH                                                           \
45                                    /* set the xtype and any elements */            \
46        switch(ctype / 3) {                                                         \
47        case 0: /* "d" */                                                           \
48            ans->xtype = CHOLMOD_REAL;                                              \
49            ans->x = (void *) REAL(GET_SLOT(x, Matrix_xSym));                       \
50            break;                                                                  \
51        case 1: /* "l" */                                                           \
52            ans->xtype = CHOLMOD_REAL;                                              \
53            ans->x = (void *) LOGICAL(GET_SLOT(x, Matrix_xSym));                    \
54            break;                                                                  \
55        case 2: /* "n" */                                                           \
56            ans->xtype = CHOLMOD_PATTERN;                                           \
57            break;                                                                  \
58        case 3: /* "z" */                                                           \
59            ans->xtype = CHOLMOD_COMPLEX;                                           \
60            ans->x = (void *) COMPLEX(GET_SLOT(x, Matrix_xSym));                    \
61            break;                                                                  \
62        }                                                                           \
63                                    /* set the stype */                             \
64        switch(ctype % 3) {                                                         \
65        case 0: /* g(eneral) */                                                     \
66            ans->stype = 0; break;                                                  \
67        case 1: /* s(ymmetric) */                                                   \
68            ans->stype =                                                            \
69                (!strcmp(CHAR(asChar(getAttrib(x, Matrix_uploSym))), "U")) ?        \
70                1 : -1;                                                             \
71            break;                                                                  \
72        case 2: /* t(riangular) */                                                  \
73            ans->stype = 0; break; /* Note that triangularity property is lost */   \
74    /*      error("triangular matrices not yet mapped to CHOLMOD"); */              \
75        }                                                                           \
76        return ans
77    
78        AS_CHM_FINISH;
79  }  }
80    
81  /**  /**
# Line 87  Line 91 
91   *   *
92   * @return SEXP containing a copy of a   * @return SEXP containing a copy of a
93   */   */
94  SEXP chm_sparse_to_SEXP(cholmod_sparse *a, int dofree, int uploT,  SEXP chm_sparse_to_SEXP(cholmod_sparse *a, int dofree, int uploT, int Rkind,
95                          char* diag, SEXP dn)                          char* diag, SEXP dn)
96  {  {
97      SEXP ans;      SEXP ans;
# Line 101  Line 105 
105                                  /* determine the class of the result */                                  /* determine the class of the result */
106      switch(a->xtype){      switch(a->xtype){
107      case CHOLMOD_PATTERN:      case CHOLMOD_PATTERN:
108          cl = uploT ? "ltCMatrix": ((a->stype) ? "lsCMatrix" : "lgCMatrix");          cl = uploT ? "ntCMatrix": ((a->stype) ? "nsCMatrix" : "ngCMatrix");
109          break;          break;
110      case CHOLMOD_REAL:      case CHOLMOD_REAL:
111            switch(Rkind) {
112            case 0:
113          cl = uploT ? "dtCMatrix": ((a->stype) ? "dsCMatrix" : "dgCMatrix");          cl = uploT ? "dtCMatrix": ((a->stype) ? "dsCMatrix" : "dgCMatrix");
114          break;          break;
115            case 1:
116                cl = uploT ? "ltCMatrix": ((a->stype) ? "lsCMatrix" : "lgCMatrix");
117                break;
118            }
119            break;
120      case CHOLMOD_COMPLEX:      case CHOLMOD_COMPLEX:
121          cl = uploT ? "ztCMatrix": ((a->stype) ? "zsCMatrix" : "zgCMatrix");          cl = uploT ? "ztCMatrix": ((a->stype) ? "zsCMatrix" : "zgCMatrix");
122          break;          break;
# Line 121  Line 132 
132             (int *) a->i, nnz);             (int *) a->i, nnz);
133                                  /* copy data slot if present */                                  /* copy data slot if present */
134      if (a->xtype == CHOLMOD_REAL)      if (a->xtype == CHOLMOD_REAL)
135            switch(Rkind) {
136            case 0:
137          Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nnz)),          Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nnz)),
138                 (double *) a->x, nnz);                 (double *) a->x, nnz);
139                break;
140            case 1:
141                Memcpy(LOGICAL(ALLOC_SLOT(ans, Matrix_xSym, LGLSXP, nnz)),
142                       (double *) a->x, nnz);
143                break;
144            }
145      if (a->xtype == CHOLMOD_COMPLEX)      if (a->xtype == CHOLMOD_COMPLEX)
146          error("complex sparse matrix code not yet written");          error("complex sparse matrix code not yet written");
147  /*      Memcpy(COMPLEX(ALLOC_SLOT(ans, Matrix_xSym, CPLXSXP, nnz)), */  /*      Memcpy(COMPLEX(ALLOC_SLOT(ans, Matrix_xSym, CPLXSXP, nnz)), */
# Line 160  Line 179 
179      cholmod_triplet *ans = Calloc(1, cholmod_triplet);      cholmod_triplet *ans = Calloc(1, cholmod_triplet);
180      char *valid[] = {"dgTMatrix", "dsTMatrix", "dtTMatrix",      char *valid[] = {"dgTMatrix", "dsTMatrix", "dtTMatrix",
181                       "lgTMatrix", "lsTMatrix", "ltTMatrix",                       "lgTMatrix", "lsTMatrix", "ltTMatrix",
182                         "ngTMatrix", "nsTMatrix", "ntTMatrix",
183                       "zgTMatrix", "zsTMatrix", "ztTMatrix",                       "zgTMatrix", "zsTMatrix", "ztTMatrix",
184                       ""};                       ""};
185      int *dims, ctype = Matrix_check_class(class_P(x), valid);      int *dims, ctype = Matrix_check_class(class_P(x), valid);
# Line 179  Line 199 
199                                  /* slots always present */                                  /* slots always present */
200      ans->i = (void *) INTEGER(islot);      ans->i = (void *) INTEGER(islot);
201      ans->j = (void *) INTEGER(GET_SLOT(x, Matrix_jSym));      ans->j = (void *) INTEGER(GET_SLOT(x, Matrix_jSym));
                                 /* set the xtype and any elements */  
     switch(ctype) {  
     case 0:  
     case 1:  
     case 2:  
         ans->xtype = CHOLMOD_REAL;  
         ans->x = (void *) REAL(GET_SLOT(x, Matrix_xSym));  
         break;  
     case 3:  
     case 4:  
     case 5:  
         ans->xtype = CHOLMOD_PATTERN;  
         break;  
     case 6:  
     case 7:  
     case 8:  
         ans->xtype = CHOLMOD_COMPLEX;  
         ans->x = (void *) COMPLEX(GET_SLOT(x, Matrix_xSym));  
         break;  
     }  
                                 /* set the stype */  
     switch(ctype % 3) {  
     case 0: ans->stype = 0; break;  
     case 1:  
         ans->stype =  
             (!strcmp(CHAR(asChar(getAttrib(x, Matrix_uploSym))), "U")) ?  
             1 : -1;  
         break;  
     case 2: ans->stype = 0; break; /* Note that triangularity property is lost */  
 /*      error("triangular matrices not yet mapped to CHOLMOD"); */  
     }  
202    
203      return ans;      AS_CHM_FINISH;
204  }  }
205    
206  /**  /**
# Line 227  Line 216 
216   *   *
217   * @return SEXP containing a copy of a   * @return SEXP containing a copy of a
218   */   */
219  SEXP chm_triplet_to_SEXP(cholmod_triplet *a, int dofree, int uploT,  SEXP chm_triplet_to_SEXP(cholmod_triplet *a, int dofree, int uploT, int Rkind,
220                           char* diag, SEXP dn)                           char* diag, SEXP dn)
221  {  {
222      SEXP ans;      SEXP ans;
# Line 238  Line 227 
227                                  /* determine the class of the result */                                  /* determine the class of the result */
228      switch(a->xtype){      switch(a->xtype){
229      case CHOLMOD_PATTERN:      case CHOLMOD_PATTERN:
230          cl = uploT ? "ltTMatrix" :          cl = uploT ? "ntTMatrix" :
231              ((a->stype) ? "lsTMatrix" : "lgTMatrix");              ((a->stype) ? "nsTMatrix" : "ngTMatrix");
232          break;          break;
233      case CHOLMOD_REAL:      case CHOLMOD_REAL:
234            switch(Rkind) {
235            case 0:
236          cl = uploT ? "dtTMatrix" :          cl = uploT ? "dtTMatrix" :
237              ((a->stype) ? "dsTMatrix" : "dgTMatrix");              ((a->stype) ? "dsTMatrix" : "dgTMatrix");
238          break;          break;
239            case 1:
240                cl = uploT ? "ltTMatrix" :
241                    ((a->stype) ? "lsTMatrix" : "lgTMatrix");
242                break;
243            }
244            break;
245      case CHOLMOD_COMPLEX:      case CHOLMOD_COMPLEX:
246          cl = uploT ? "ztTMatrix" :          cl = uploT ? "ztTMatrix" :
247              ((a->stype) ? "zsTMatrix" : "zgTMatrix");              ((a->stype) ? "zsTMatrix" : "zgTMatrix");
# Line 261  Line 258 
258             (int *) a->j, a->nnz);             (int *) a->j, a->nnz);
259                                  /* copy data slot if present */                                  /* copy data slot if present */
260      if (a->xtype == CHOLMOD_REAL)      if (a->xtype == CHOLMOD_REAL)
261            switch(Rkind) {
262            case 0:
263          Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, a->nnz)),          Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, a->nnz)),
264                 (double *) a->x, a->nnz);                 (double *) a->x, a->nnz);
265      if (a->xtype == CHOLMOD_COMPLEX)              break;
266            case 1:
267                Memcpy(LOGICAL(ALLOC_SLOT(ans, Matrix_xSym, LGLSXP, a->nnz)),
268                       (double *) a->x, a->nnz);
269                break;
270            }
271        else if (a->xtype == CHOLMOD_COMPLEX)
272          error("complex sparse matrix code not yet written");          error("complex sparse matrix code not yet written");
273  /*      Memcpy(COMPLEX(ALLOC_SLOT(ans, Matrix_xSym, CPLXSXP, a->nnz)), */  /*      Memcpy(COMPLEX(ALLOC_SLOT(ans, Matrix_xSym, CPLXSXP, a->nnz)), */
274  /*             (complex *) a->x, a->nz); */  /*             (complex *) a->x, a->nz); */
# Line 299  Line 304 
304      cholmod_dense *ans = Calloc(1, cholmod_dense);      cholmod_dense *ans = Calloc(1, cholmod_dense);
305      char *valid[] = {"dmatrix", "dgeMatrix",      char *valid[] = {"dmatrix", "dgeMatrix",
306                       "lmatrix", "lgeMatrix",                       "lmatrix", "lgeMatrix",
307                         "nmatrix", "ngeMatrix",
308                       "zmatrix", "zgeMatrix", ""};                       "zmatrix", "zgeMatrix", ""};
309      int dims[2], ctype = Matrix_check_class(class_P(x), valid), nprot = 0;      int dims[2], ctype = Matrix_check_class(class_P(x), valid), nprot = 0;
310    
# Line 310  Line 316 
316              nprot++;              nprot++;
317          }          }
318          ctype = (isReal(x) ? 0 :          ctype = (isReal(x) ? 0 :
319                   (isLogical(x) ? 2 :                   (isLogical(x) ? 2 : /* logical -> default to "l", not "n" */
320                    (isComplex(x) ? 4 : -1)));                    (isComplex(x) ? 6 : -1)));
321      } else Memcpy(dims, INTEGER(GET_SLOT(x, Matrix_DimSym)), 2);      } else Memcpy(dims, INTEGER(GET_SLOT(x, Matrix_DimSym)), 2);
322      if (ctype < 0) error("invalid class of object to as_cholmod_dense");      if (ctype < 0) error("invalid class of object to as_cholmod_dense");
323                                  /* characteristics of the system */                                  /* characteristics of the system */
# Line 323  Line 329 
329      ans->nzmax = dims[0] * dims[1];      ans->nzmax = dims[0] * dims[1];
330                                  /* set the xtype and any elements */                                  /* set the xtype and any elements */
331      switch(ctype / 2) {      switch(ctype / 2) {
332      case 0:      case 0: /* "d" & "l" */
333        case 1:
334          ans->xtype = CHOLMOD_REAL;          ans->xtype = CHOLMOD_REAL;
335          ans->x = (void *) REAL((ctype % 2) ? GET_SLOT(x, Matrix_xSym) : x);          ans->x = (void *) REAL((ctype % 2) ? GET_SLOT(x, Matrix_xSym) : x);
336          break;          break;
337      case 1:      case 2: /* "n" */
338          ans->xtype = CHOLMOD_PATTERN;          ans->xtype = CHOLMOD_PATTERN;
339          ans->x = (void *) LOGICAL((ctype % 2) ? GET_SLOT(x, Matrix_xSym) : x);          ans->x = (void *) LOGICAL((ctype % 2) ? GET_SLOT(x, Matrix_xSym) : x);
340          break;          break;
341      case 3:      case 3: /* "z" */
342          ans->xtype = CHOLMOD_COMPLEX;          ans->xtype = CHOLMOD_COMPLEX;
343          ans->x = (void *) COMPLEX((ctype % 2) ? GET_SLOT(x, Matrix_xSym) : x);          ans->x = (void *) COMPLEX((ctype % 2) ? GET_SLOT(x, Matrix_xSym) : x);
344          break;          break;
# Line 349  Line 356 
356   *   *
357   * @return SEXP containing a copy of a   * @return SEXP containing a copy of a
358   */   */
359  SEXP chm_dense_to_SEXP(cholmod_dense *a, int dofree)  SEXP chm_dense_to_SEXP(cholmod_dense *a, int dofree, int Rkind)
360  {  {
361      SEXP ans;      SEXP ans;
362      char *cl;      char *cl = ""; /* -Wall */
363      int *dims, ntot;      int *dims, ntot;
364                                  /* determine the class of the result */                                  /* determine the class of the result */
365      cl = (a->xtype == CHOLMOD_PATTERN) ? "lgeMatrix" :      switch(a->xtype) {
366          ((a->xtype == CHOLMOD_REAL) ? "dgeMatrix" :      case CHOLMOD_PATTERN:
367           ((a->xtype == CHOLMOD_COMPLEX) ? "zgeMatrix" : ""));          cl = "ngeMatrix"; break;
368      if (!strlen(cl)) error("unknown xtype");      case CHOLMOD_REAL:
369            switch(Rkind) {
370            case 0:  cl = "dgeMatrix"; break;
371            case 1:  cl = "lgeMatrix"; break;
372            default: error("unknown 'Rkind'");
373            }
374            break;
375        case CHOLMOD_COMPLEX:
376            cl = "zgeMatrix"; break;
377        default:
378            error("unknown xtype");
379        }
380    
381      ans = PROTECT(NEW_OBJECT(MAKE_CLASS(cl)));      ans = PROTECT(NEW_OBJECT(MAKE_CLASS(cl)));
382                                  /* allocate and copy common slots */                                  /* allocate and copy common slots */
# Line 367  Line 385 
385      ntot = dims[0] * dims[1];      ntot = dims[0] * dims[1];
386      if (a->d == a->nrow) {      /* copy data slot if present */      if (a->d == a->nrow) {      /* copy data slot if present */
387          if (a->xtype == CHOLMOD_REAL)          if (a->xtype == CHOLMOD_REAL)
388                switch(Rkind) {
389                case 0:
390              Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, ntot)),              Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, ntot)),
391                     (double *) a->x, ntot);                     (double *) a->x, ntot);
392          if (a->xtype == CHOLMOD_COMPLEX)                  break;
393                case 1:
394                    Memcpy(LOGICAL(ALLOC_SLOT(ans, Matrix_xSym, LGLSXP, ntot)),
395                           (double *) a->x, ntot);
396                    break;
397                }
398            else if (a->xtype == CHOLMOD_COMPLEX)
399              error("complex sparse matrix code not yet written");              error("complex sparse matrix code not yet written");
400  /*      Memcpy(COMPLEX(ALLOC_SLOT(ans, Matrix_xSym, CPLXSXP, a->nnz)), */  /*      Memcpy(COMPLEX(ALLOC_SLOT(ans, Matrix_xSym, CPLXSXP, a->nnz)), */
401  /*             (complex *) a->x, a->nz); */  /*             (complex *) a->x, a->nz); */
# Line 448  Line 474 
474  cholmod_factor *as_cholmod_factor(SEXP x)  cholmod_factor *as_cholmod_factor(SEXP x)
475  {  {
476      cholmod_factor *ans = Calloc(1, cholmod_factor);      cholmod_factor *ans = Calloc(1, cholmod_factor);
477      char *valid[] = {"dCHMsuper", "dCHMsimpl", "lCHMsuper", "lCHMsimpl", ""};      char *valid[] = {"dCHMsuper", "dCHMsimpl", "nCHMsuper", "nCHMsimpl", ""};
478      int *type = INTEGER(GET_SLOT(x, install("type"))),      int *type = INTEGER(GET_SLOT(x, install("type"))),
479          ctype = Matrix_check_class(class_P(x), valid);          ctype = Matrix_check_class(class_P(x), valid);
480      SEXP tmp;      SEXP tmp;
# Line 531  Line 557 
557          class = f->is_super ? "dCHMsuper" : "dCHMsimpl";          class = f->is_super ? "dCHMsuper" : "dCHMsimpl";
558          break;          break;
559      case CHOLMOD_PATTERN:      case CHOLMOD_PATTERN:
560          class = f->is_super ? "lCHMsuper" : "lCHMsimpl";          class = f->is_super ? "nCHMsuper" : "nCHMsimpl";
561          break;          break;
562      default:      default:
563          error(_("f->xtype of %d not recognized"), f->xtype);          error(_("f->xtype of %d not recognized"), f->xtype);

Legend:
Removed from v.1547  
changed lines
  Added in v.1548

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