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 2222, Fri Jul 18 10:00:25 2008 UTC revision 2223, Fri Jul 18 23:04:48 2008 UTC
# Line 16  Line 16 
16   *   *
17   * @return ans containing pointers to the slots of x.   * @return ans containing pointers to the slots of x.
18   */   */
19  CHM_SP as_cholmod_sparse(CHM_SP ans, SEXP x)  CHM_SP as_cholmod_sparse(CHM_SP ans, SEXP x, Rboolean check_Udiag)
20  {  {
21      char *valid[] = {"dgCMatrix", "dsCMatrix", "dtCMatrix",      char *valid[] = {"dgCMatrix", "dsCMatrix", "dtCMatrix",
22                       "lgCMatrix", "lsCMatrix", "ltCMatrix",                       "lgCMatrix", "lsCMatrix", "ltCMatrix",
# Line 25  Line 25 
25                       ""};                       ""};
26      int *dims = INTEGER(GET_SLOT(x, Matrix_DimSym)),      int *dims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
27          ctype = Matrix_check_class(class_P(x), valid);          ctype = Matrix_check_class(class_P(x), valid);
28        Rboolean do_Udiag = FALSE;
29      SEXP islot = GET_SLOT(x, Matrix_iSym);      SEXP islot = GET_SLOT(x, Matrix_iSym);
30    
31      if (ctype < 0) error("invalid class of object to as_cholmod_sparse");      if (ctype < 0) error("invalid class of object to as_cholmod_sparse");
# Line 44  Line 45 
45      ans->i = (void *) INTEGER(islot);      ans->i = (void *) INTEGER(islot);
46      ans->p = (void *) INTEGER(GET_SLOT(x, Matrix_pSym));      ans->p = (void *) INTEGER(GET_SLOT(x, Matrix_pSym));
47    
48  #define AS_CHM_FINISH                                                           \  #define AS_CHM_COMMON                                                           \
49                                    /* set the stype */                             \
50        switch(ctype % 3) {                                                         \
51        case 0: /* g(eneral) */                                                     \
52            ans->stype = 0; break;                                                  \
53        case 1: /* s(ymmetric) */                                                   \
54            ans->stype = (*uplo_P(x) == 'U') ? 1 : -1;                              \
55            break;                                                                  \
56        case 2: /* t(riangular) -- Note that triangularity property is lost! */     \
57            ans->stype = 0;                                                         \
58            do_Udiag = (check_Udiag && (*diag_P(x) == 'U'));                        \
59            break;                                                                  \
60        }                                                                           \
61                                  /* set the xtype and any elements */            \                                  /* set the xtype and any elements */            \
62      switch(ctype / 3) {                                                         \      switch(ctype / 3) {                                                         \
63      case 0: /* "d" */                                                           \      case 0: /* "d" */                                                           \
# Line 62  Line 75 
75          ans->xtype = CHOLMOD_COMPLEX;                                           \          ans->xtype = CHOLMOD_COMPLEX;                                           \
76          ans->x = (void *) COMPLEX(GET_SLOT(x, Matrix_xSym));                    \          ans->x = (void *) COMPLEX(GET_SLOT(x, Matrix_xSym));                    \
77          break;                                                                  \          break;                                                                  \
78      }                                                                           \      }
                                 /* set the stype */                             \  
     switch(ctype % 3) {                                                         \  
     case 0: /* g(eneral) */                                                     \  
         ans->stype = 0; break;                                                  \  
     case 1: /* s(ymmetric) */                                                   \  
         ans->stype = (*uplo_P(x) == 'U') ? 1 : -1;                              \  
         break;                                                                  \  
     case 2: /* t(riangular) -- Note that triangularity property is lost! */     \  
         ans->stype = 0;                                                         \  
         /* NOTE: if(*diag_P(x) == 'U'), the diagonal is lost (!); */            \  
         /* ---- that may be ok, e.g. if we are just converting from/to Tsparse, */ \  
         /*      but is *not* at all ok, e.g. when used before matrix products */   \  
         break;                                                                  \  
     }                                                                           \  
     return ans  
79    
80      AS_CHM_FINISH;      AS_CHM_COMMON;
81    
82        if(do_Udiag) {  /* diagU2N(.)  "in place" : */
83            double one[] = {1, 0};
84            CHM_SP tmp = cholmod_copy (ans, ans->stype, ans->xtype, &c);
85            CHM_SP eye = cholmod_speye(ans->nrow, ans->ncol, ans->xtype, &c);
86    
87            ans = cholmod_add(tmp, eye, one, one, TRUE, TRUE, &c);
88    
89            cholmod_free_sparse(&tmp, &c);
90            cholmod_free_sparse(&eye, &c);
91        } /* else :
92           * NOTE: if(*diag_P(x) == 'U'), the diagonal is lost (!);
93           * ---- that may be ok, e.g. if we are just converting from/to Tsparse,
94           *      but is *not* at all ok, e.g. when used before matrix products */
95    
96        return ans;
97  }  }
98    
99  /**  /**
# Line 186  Line 200 
200   *   *
201   * @return ans containing pointers to the slots of x.   * @return ans containing pointers to the slots of x.
202   */   */
203  CHM_TR as_cholmod_triplet(CHM_TR ans, SEXP x)  CHM_TR as_cholmod_triplet(CHM_TR ans, SEXP x, Rboolean check_Udiag)
204  {  {
205      char *valid[] = {"dgTMatrix", "dsTMatrix", "dtTMatrix",      char *valid[] = {"dgTMatrix", "dsTMatrix", "dtTMatrix",
206                       "lgTMatrix", "lsTMatrix", "ltTMatrix",                       "lgTMatrix", "lsTMatrix", "ltTMatrix",
207                       "ngTMatrix", "nsTMatrix", "ntTMatrix",                       "ngTMatrix", "nsTMatrix", "ntTMatrix",
208                       "zgTMatrix", "zsTMatrix", "ztTMatrix",                       "zgTMatrix", "zsTMatrix", "ztTMatrix",
209                       ""};                       ""};
210      int *dims, ctype = Matrix_check_class(class_P(x), valid);      int *dims, ctype = Matrix_check_class(class_P(x), valid), m;
211      SEXP islot;      SEXP islot;
212        Rboolean do_Udiag;
213    
214      if (ctype < 0) error("invalid class of object to as_cholmod_triplet");      if (ctype < 0) error("invalid class of object to as_cholmod_triplet");
215      memset(ans, 0, sizeof(cholmod_triplet)); /* zero the struct */      memset(ans, 0, sizeof(cholmod_triplet)); /* zero the struct */
# Line 206  Line 221 
221      dims = INTEGER(GET_SLOT(x, Matrix_DimSym));      dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
222      ans->nrow = dims[0];      ans->nrow = dims[0];
223      ans->ncol = dims[1];      ans->ncol = dims[1];
224    
225        do_Udiag = (check_Udiag && ctype % 3 == 2 && (*diag_P(x) == 'U'));
226    
227      islot = GET_SLOT(x, Matrix_iSym);      islot = GET_SLOT(x, Matrix_iSym);
228      ans->nnz = ans->nzmax = LENGTH(islot);      m = LENGTH(islot);
229        ans->nnz = ans->nzmax = do_Udiag ? m + dims[0] : m;
230                                  /* slots always present */                                  /* slots always present */
231      ans->i = (void *) INTEGER(islot);      ans->i = (void *) INTEGER(islot);
232      ans->j = (void *) INTEGER(GET_SLOT(x, Matrix_jSym));      ans->j = (void *) INTEGER(GET_SLOT(x, Matrix_jSym));
233    
234      AS_CHM_FINISH;      AS_CHM_COMMON;
235    
236        if(do_Udiag) {  /* Tsparse_diagU2N(.)  [ ./Tsparse.c ]  "in place" : */
237            int k = m + dims[0];
238            int *a_i, *a_j;
239    
240            /* TODO? instead of reallocating, don't do the 2nd part of AS_CHM_COMMON()
241             * ---- above, and allocate to correct length + Memcpy() here, as in
242             * Tsparse_diagU2N() */
243            if(cholmod_reallocate_triplet((size_t) k, ans, &c))
244                error(_("as_cholmod_triplet(): could not reallocate for internal diagU2N()"
245                          ));
246            a_i = ans->i;
247            a_j = ans->j;
248            /* add (@i, @j)[k+m] = k, @x[k+m] = 1.   for k = 0,..,(n-1) */
249            for(k=0; k < dims[0]; k++) {
250                a_i[k+m] = k;
251                a_j[k+m] = k;
252    
253                switch(ctype / 3) {
254                case 0: { /* "d" */
255                    double *a_x = ans->x;
256                    a_x[k+m] = 1.;
257                    break;
258                }
259                case 1: { /* "l" */
260                    int *a_x = ans->x;
261                    a_x[k+m] = 1;
262                    break;
263                }
264                case 2: /* "n" */
265                    break;
266                case 3: { /* "z" */
267                    double *a_x = ans->x;
268                    a_x[2*(k+m)  ] = 1.;
269                    a_x[2*(k+m)+1] = 0.;
270                    break;
271                }
272                }
273            }
274    
275        } /* else :
276           * NOTE: if(*diag_P(x) == 'U'), the diagonal is lost (!);
277           * ---- that may be ok, e.g. if we are just converting from/to Tsparse,
278           *      but is *not* at all ok, e.g. when used before matrix products */
279    
280        return ans;
281  }  }
282    
283  /**  /**

Legend:
Removed from v.2222  
changed lines
  Added in v.2223

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