SCM

SCM Repository

[matrix] Annotation of /pkg/Matrix/src/Csparse.c
ViewVC logotype

Annotation of /pkg/Matrix/src/Csparse.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1029 - (view) (download) (as text)
Original Path: pkg/src/Csparse.c

1 : maechler 925 /* Sparse matrices in compress column-oriented form */
2 : bates 922 #include "Csparse.h"
3 :     #ifdef USE_CHOLMOD
4 :     #include "chm_common.h"
5 :     #endif /* USE_CHOLMOD */
6 :    
7 :     SEXP Csparse_validate(SEXP x)
8 :     {
9 :     SEXP pslot = GET_SLOT(x, Matrix_pSym),
10 :     islot = GET_SLOT(x, Matrix_iSym);
11 :     int j, ncol = length(pslot) - 1,
12 :     *dims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
13 :     nrow, *xp = INTEGER(pslot),
14 :     *xi = INTEGER(islot);
15 :    
16 :     nrow = dims[0];
17 :     if (length(pslot) <= 0)
18 :     return mkString(_("slot p must have length > 0"));
19 :     if (xp[0] != 0)
20 :     return mkString(_("first element of slot p must be zero"));
21 :     if (length(islot) != xp[ncol])
22 :     return mkString(_("last element of slot p must match length of slots i and x"));
23 :     for (j = 0; j < ncol; j++) {
24 :     if (xp[j] > xp[j+1])
25 :     return mkString(_("slot p must be non-decreasing"));
26 :     }
27 :     for (j = 0; j < length(islot); j++) {
28 :     if (xi[j] < 0 || xi[j] >= nrow)
29 :     return mkString(_("all row indices must be between 0 and nrow-1"));
30 :     }
31 :     return ScalarLogical(1);
32 :     }
33 :    
34 :     SEXP Csparse_to_Tsparse(SEXP x)
35 :     {
36 :     #ifdef USE_CHOLMOD
37 : maechler 925 cholmod_sparse *chxs = as_cholmod_sparse(x);
38 : bates 922 cholmod_triplet *chxt = cholmod_sparse_to_triplet(chxs, &c);
39 :    
40 : maechler 1029 free(chxs);
41 : bates 922 return chm_triplet_to_SEXP(chxt, 1);
42 :     #else
43 :     error("General conversion requires CHOLMOD");
44 :     return R_NilValue; /* -Wall */
45 :     #endif /* USE_CHOLMOD */
46 :     }
47 :    
48 :     SEXP Csparse_transpose(SEXP x)
49 :     {
50 :     #ifdef USE_CHOLMOD
51 :     cholmod_sparse *chx = as_cholmod_sparse(x);
52 :     cholmod_sparse *chxt = cholmod_transpose(chx, (int) chx->xtype, &c);
53 :    
54 : maechler 1029 free(chx);
55 : bates 922 return chm_sparse_to_SEXP(chxt, 1);
56 :     #else
57 :     error("General conversion requires CHOLMOD");
58 :     return R_NilValue; /* -Wall */
59 :     #endif /* USE_CHOLMOD */
60 :     }
61 :    
62 :    
63 :     SEXP Csparse_Csparse_prod(SEXP a, SEXP b)
64 :     {
65 :     #ifdef USE_CHOLMOD
66 :     cholmod_sparse *cha = as_cholmod_sparse(a), *chb = as_cholmod_sparse(b);
67 :     cholmod_sparse *chc = cholmod_ssmult(cha, chb, 0, (int) cha->xtype, 1, &c);
68 :    
69 : maechler 1029 free(cha); free(chb);
70 : bates 922 return chm_sparse_to_SEXP(chc, 1);
71 :     #else
72 :     error("General multiplication requires CHOLMOD");
73 :     return R_NilValue; /* -Wall */
74 :     #endif /* USE_CHOLMOD */
75 :     }
76 :    
77 :     SEXP Csparse_dense_prod(SEXP a, SEXP b)
78 :     {
79 :     #ifdef USE_CHOLMOD
80 :     cholmod_sparse *cha = as_cholmod_sparse(a);
81 :     cholmod_dense *chb = as_cholmod_dense(b);
82 :     cholmod_dense *chc = cholmod_allocate_dense(cha->nrow, chb->ncol,
83 :     cha->nrow, chb->xtype, &c);
84 :     double alpha = 1, beta = 0;
85 :    
86 :     cholmod_sdmult(cha, 0, &alpha, &beta, chb, chc, &c);
87 : maechler 1029 free(cha); free(chb);
88 : bates 923 return chm_dense_to_SEXP(chc, 1);
89 : bates 922 #else
90 :     error("General multiplication requires CHOLMOD");
91 :     return R_NilValue; /* -Wall */
92 :     #endif /* USE_CHOLMOD */
93 :     }
94 : maechler 925
95 : bates 928 SEXP Csparse_crossprod(SEXP x, SEXP trans, SEXP triplet)
96 : bates 922 {
97 :     #ifdef USE_CHOLMOD
98 : maechler 957 int trip = asLogical(triplet),
99 :     tr = asLogical(trans); /* gets reversed because _aat is tcrossprod */
100 : bates 928 cholmod_triplet
101 :     *cht = trip ? as_cholmod_triplet(x) : (cholmod_triplet*) NULL;
102 :     cholmod_sparse *chcp, *chxt,
103 :     *chx = trip ? cholmod_triplet_to_sparse(cht, cht->nnz, &c)
104 :     : as_cholmod_sparse(x);
105 : bates 922
106 : bates 923 if (!tr)
107 :     chxt = cholmod_transpose(chx, (int) chx->xtype, &c);
108 : bates 928 chcp = cholmod_aat((!tr) ? chxt : chx, (int *) NULL, 0, chx->xtype, &c);
109 : maechler 957 if(!chcp)
110 :     error("Csparse_crossprod(): error return from cholmod_aat()");
111 : bates 923
112 : bates 930 if (trip) {
113 :     cholmod_free_sparse(&chx, &c);
114 : maechler 1029 free(cht);
115 : bates 930 } else {
116 : maechler 1029 free(chx);
117 : bates 930 }
118 : bates 923 if (!tr) cholmod_free_sparse(&chxt, &c);
119 :     return chm_sparse_to_SEXP(chcp, 1);
120 : bates 922 #else
121 : bates 930 error("General crossproduct requires CHOLMOD");
122 : bates 922 return R_NilValue; /* -Wall */
123 :     #endif /* USE_CHOLMOD */
124 :     }
125 : bates 923

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