SCM

SCM Repository

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

Annotation of /pkg/src/cs_utils.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2234 - (view) (download) (as text)

1 : bates 1250 #include "cs_utils.h"
2 :    
3 :     /* Borrowed from one of Tim Davis' examples in the CSparse Demo directory */
4 :     /* 1 if A is square & upper tri., -1 if square & lower tri., 0 otherwise */
5 :     static int is_sym (cs *A)
6 :     {
7 :     int is_upper, is_lower, j, p, n = A->n, m = A->m, *Ap = A->p, *Ai = A->i ;
8 :     if (m != n) return (0) ;
9 :     is_upper = 1 ;
10 :     is_lower = 1 ;
11 :     for (j = 0 ; j < n ; j++)
12 :     {
13 :     for (p = Ap [j] ; p < Ap [j+1] ; p++)
14 :     {
15 :     if (Ai [p] > j) is_upper = 0 ;
16 :     if (Ai [p] < j) is_lower = 0 ;
17 :     }
18 :     }
19 :     return (is_upper ? 1 : (is_lower ? -1 : 0)) ;
20 :     }
21 :    
22 : dmbates 2216
23 : bates 1250 /**
24 : dmbates 2216 * Create an identity matrix of size n as a cs struct. The structure
25 :     * must be freed with cs_free by the caller.
26 :     *
27 :     * @param n size of identity matrix to construct.
28 :     *
29 :     * @return pointer to a cs object that contains the identity matrix.
30 :     */
31 :     static CSP csp_eye(int n)
32 :     {
33 :     CSP eye = cs_spalloc(n, n, n, 1, 0);
34 :     int *ep = eye->p, *ei = eye->i;
35 :     double *ex = eye->x;
36 : mmaechler 2227
37 : dmbates 2216 if (n <= 0) error("csp_eye argument n must be positive");
38 : dmbates 2228 eye->nz = -1; /* compressed column storage */
39 : dmbates 2216 for (int j = 0; j < n; j++) {
40 :     ep[j] = ei[j] = j;
41 :     ex[j] = 1;
42 :     }
43 : dmbates 2228 eye->nzmax = ep[n] = n;
44 : dmbates 2216 return eye;
45 :     }
46 :    
47 :     /**
48 : maechler 2115 * Create a cs object with the contents of x. Typically called via AS_CSP()
49 : bates 1250 *
50 : bates 2104 * @param ans pointer to a cs struct. This is allocated in the caller
51 :     * so it is easier to keep track of where it should be freed - in many
52 :     * applications the memory can be allocated with alloca and
53 :     * automatically freed on exit from the caller.
54 : bates 1250 * @param x pointer to an object that inherits from CsparseMatrix
55 :     *
56 :     * @return pointer to a cs object that contains pointers
57 :     * to the slots of x.
58 :     */
59 : mmaechler 2227 cs *Matrix_as_cs(cs *ans, SEXP x, Rboolean check_Udiag)
60 : bates 1250 {
61 : mmaechler 2227 char *valid[] = {"dgCMatrix", "dtCMatrix", ""};
62 :     /* had also "dsCMatrix", but that only stores one triangle */
63 : maechler 1945 int *dims, ctype = Matrix_check_class(class_P(x), valid);
64 : bates 1250 SEXP islot;
65 :    
66 : mmaechler 2234 if (ctype < 0) error("invalid class of 'x' in Matrix_as_cs(a, x)");
67 : bates 1250 /* dimensions and nzmax */
68 :     dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
69 :     ans->m = dims[0]; ans->n = dims[1];
70 :     islot = GET_SLOT(x, Matrix_iSym);
71 : bates 1256 ans->nz = -1; /* indicates compressed column storage */
72 :     ans->nzmax = LENGTH(islot);
73 : bates 1250 ans->i = INTEGER(islot);
74 :     ans->p = INTEGER(GET_SLOT(x, Matrix_pSym));
75 :     ans->x = REAL(GET_SLOT(x, Matrix_xSym));
76 :    
77 : dmbates 2228 if(check_Udiag && ctype == 1 && (*diag_P(x) == 'U') && ans->m == ans->n) { /* diagU2N(.) : */
78 : mmaechler 2227 int n = dims[0];
79 :     CSP I_n = csp_eye(n);
80 :     /* tmp := 1*ans + 1*eye -- result is newly allocated in cs_add(): */
81 :     CSP tmp = cs_add(ans, I_n, 1., 1.);
82 : dmbates 2228 int nz = (tmp->p)[n];
83 : mmaechler 2227
84 :     /* content(ans) := content(tmp) : */
85 :     ans->nzmax = nz;
86 : mmaechler 2230 /* The ans "slots" were pointers to x@ <slots>; all need new content now: */
87 : mmaechler 2234 ans->p = Memcpy(( int*) R_alloc(sizeof( int), n+1),
88 :     ( int*) tmp->p, n+1);
89 :     ans->i = Memcpy(( int*) R_alloc(sizeof( int), nz),
90 :     ( int*) tmp->i, nz);
91 :     ans->x = Memcpy((double*) R_alloc(sizeof(double), nz),
92 :     (double*) tmp->x, nz);
93 : mmaechler 2227
94 :     cs_spfree(I_n);
95 :     cs_spfree(tmp);
96 :     }
97 : bates 1250 return ans;
98 :     }
99 :    
100 : maechler 1960 /**
101 :     * Copy the contents of a to an appropriate CsparseMatrix object and,
102 :     * optionally, free a or free both a and the pointers to its contents.
103 :     *
104 :     * @param a matrix to be converted
105 :     * @param cl the name of the S4 class of the object to be generated
106 :     * @param dofree 0 - don't free a; > 0 cs_free a; < 0 Free a
107 :     *
108 :     * @return SEXP containing a copy of a
109 :     */
110 :     SEXP Matrix_cs_to_SEXP(cs *a, char *cl, int dofree)
111 :     {
112 :     SEXP ans;
113 :     char *valid[] = {"dgCMatrix", "dsCMatrix", "dtCMatrix", ""};
114 :     int *dims, ctype = Matrix_check_class(cl, valid), nz;
115 :    
116 :     if (ctype < 0)
117 : maechler 2115 error(_("invalid class of object to Matrix_cs_to_SEXP"));
118 : maechler 1960 ans = PROTECT(NEW_OBJECT(MAKE_CLASS(cl)));
119 :     /* allocate and copy common slots */
120 :     dims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));
121 :     dims[0] = a->m; dims[1] = a->n;
122 :     Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, a->n + 1)),
123 :     a->p, a->n + 1);
124 :     nz = a->p[a->n];
125 :     Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nz)), a->i, nz);
126 :     Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), a->x, nz);
127 : maechler 2115 if (ctype > 0) { /* dsC or dtC */
128 : maechler 1960 int uplo = is_sym(a);
129 : maechler 2115 if (!uplo)
130 :     error(_("cs matrix not compatible with class '%s'"), valid[ctype]);
131 :     if (ctype == 2) /* dtC* */
132 :     SET_SLOT(ans, Matrix_diagSym, mkString("N"));
133 : maechler 1960 SET_SLOT(ans, Matrix_uploSym, mkString(uplo < 0 ? "L" : "U"));
134 :     }
135 :     if (dofree > 0) cs_spfree(a);
136 :     if (dofree < 0) Free(a);
137 :     UNPROTECT(1);
138 :     return ans;
139 :     }
140 :    
141 : maechler 2115 #if 0 /* unused ------------------------------------*/
142 : mmaechler 2227 /* -------------------------------------*/
143 :    
144 : bates 1250 /**
145 : maechler 1960 * Populate a css object with the contents of x.
146 : bates 1380 *
147 : maechler 1960 * @param ans pointer to a csn object
148 : bates 1380 * @param x pointer to an object of class css_LU or css_QR.
149 :     *
150 :     * @return pointer to a cs object that contains pointers
151 :     * to the slots of x.
152 :     */
153 : maechler 1960 css *Matrix_as_css(css *ans, SEXP x)
154 : bates 1380 {
155 : maechler 1736 char *cl = class_P(x);
156 : bates 1380 *valid[] = {"css_LU", "css_QR", ""};
157 :     int *nz = INTEGER(GET_SLOT(x, install("nz"))),
158 : maechler 1945 ctype = Matrix_check_class(cl, valid);
159 : bates 1380
160 :     if (ctype < 0) error("invalid class of object to Matrix_as_css");
161 : bates 1383 ans->q = INTEGER(GET_SLOT(x, install("Q")));
162 : bates 1380 ans->m2 = nz[0]; ans->lnz = nz[1]; ans->unz = nz[2];
163 :     switch(ctype) {
164 :     case 0: /* css_LU */
165 : bates 1383 ans->pinv = (int *) NULL;
166 : bates 1380 ans->parent = (int *) NULL;
167 :     ans->cp = (int *) NULL;
168 :     break;
169 :     case 1: /* css_QR */
170 : bates 1383 ans->pinv = INTEGER(GET_SLOT(x, install("Pinv")));
171 : bates 1380 ans->parent = INTEGER(GET_SLOT(x, install("parent")));
172 :     ans->cp = INTEGER(GET_SLOT(x, install("cp")));
173 :     break;
174 :     default:
175 :     error("invalid class of object to Matrix_as_css");
176 :     }
177 :     return ans;
178 :     }
179 :    
180 :     /**
181 : maechler 1960 * Populate a csn object with the contents of x.
182 : bates 1380 *
183 : maechler 1960 * @param ans pointer to a csn object
184 : bates 1380 * @param x pointer to an object of class csn_LU or csn_QR.
185 :     *
186 :     * @return pointer to a cs object that contains pointers
187 :     * to the slots of x.
188 :     */
189 : maechler 1960 csn *Matrix_as_csn(csn *ans, SEXP x)
190 : bates 1380 {
191 : bates 1461 char *valid[] = {"csn_LU", "csn_QR", ""};
192 : maechler 1945 int ctype = Matrix_check_class(class_P(x), valid);
193 : bates 1380
194 :     if (ctype < 0) error("invalid class of object to Matrix_as_csn");
195 :     ans->U = Matrix_as_cs(GET_SLOT(x, install("U")));
196 :     ans->L = Matrix_as_cs(GET_SLOT(x, install("L")));
197 :     switch(ctype) {
198 :     case 0:
199 :     ans->B = (double*) NULL;
200 : bates 1383 ans->pinv = INTEGER(GET_SLOT(x, install("Pinv")));
201 : bates 1380 break;
202 :     case 1:
203 :     ans->B = REAL(GET_SLOT(x, install("beta")));
204 : bates 1383 ans->pinv = (int*) NULL;
205 : bates 1380 break;
206 :     default:
207 :     error("invalid class of object to Matrix_as_csn");
208 :     }
209 :     return ans;
210 :     }
211 :    
212 :     /**
213 :     * Copy the contents of S to a css_LU or css_QR object and,
214 :     * optionally, free S or free both S and the pointers to its contents.
215 :     *
216 :     * @param a css object to be converted
217 :     * @param cl the name of the S4 class of the object to be generated
218 :     * @param dofree 0 - don't free a; > 0 cs_free a; < 0 Free a
219 :     * @param m number of rows in original matrix
220 :     * @param n number of columns in original matrix
221 :     *
222 :     * @return SEXP containing a copy of S
223 :     */
224 : bates 1381 SEXP Matrix_css_to_SEXP(css *S, char *cl, int dofree, int m, int n)
225 : bates 1380 {
226 :     SEXP ans;
227 :     char *valid[] = {"css_LU", "css_QR", ""};
228 : maechler 1945 int *nz, ctype = Matrix_check_class(cl, valid);
229 : maechler 1736
230 : bates 1380 if (ctype < 0)
231 : maechler 1906 error("Inappropriate class '%s' for Matrix_css_to_SEXP", cl);
232 : bates 1380 ans = PROTECT(NEW_OBJECT(MAKE_CLASS(cl)));
233 :     /* allocate and copy common slots */
234 : bates 1383 Memcpy(INTEGER(ALLOC_SLOT(ans, install("Q"), INTSXP, n)), S->q, n);
235 : bates 1380 nz = INTEGER(ALLOC_SLOT(ans, install("nz"), INTSXP, 3));
236 :     nz[0] = S->m2; nz[1] = S->lnz; nz[2] = S->unz;
237 :     switch(ctype) {
238 :     case 0:
239 :     break;
240 :     case 1:
241 :     Memcpy(INTEGER(ALLOC_SLOT(ans, install("Pinv"), INTSXP, m)),
242 : bates 1383 S->pinv, m);
243 : bates 1380 Memcpy(INTEGER(ALLOC_SLOT(ans, install("parent"), INTSXP, n)),
244 :     S->parent, n);
245 :     Memcpy(INTEGER(ALLOC_SLOT(ans, install("cp"), INTSXP, n)),
246 :     S->cp, n);
247 :     break;
248 :     default:
249 : maechler 1906 error("Inappropriate class '%s' for Matrix_css_to_SEXP", cl);
250 : bates 1380 }
251 :     if (dofree > 0) cs_sfree(S);
252 :     if (dofree < 0) Free(S);
253 :     UNPROTECT(1);
254 :     return ans;
255 :     }
256 :    
257 :     /**
258 :     * Copy the contents of N to a csn_LU or csn_QR object and,
259 :     * optionally, free N or free both N and the pointers to its contents.
260 :     *
261 :     * @param a csn object to be converted
262 :     * @param cl the name of the S4 class of the object to be generated
263 :     * @param dofree 0 - don't free a; > 0 cs_free a; < 0 Free a
264 :     *
265 :     * @return SEXP containing a copy of S
266 :     */
267 : bates 1381 SEXP Matrix_csn_to_SEXP(csn *N, char *cl, int dofree)
268 : bates 1380 {
269 :     SEXP ans;
270 :     char *valid[] = {"csn_LU", "csn_QR", ""};
271 : maechler 1945 int ctype = Matrix_check_class(cl, valid), n = (N->U)->n;
272 : maechler 1736
273 : bates 1380 if (ctype < 0)
274 : maechler 1906 error("Inappropriate class '%s' for Matrix_csn_to_SEXP", cl);
275 : bates 1380 ans = PROTECT(NEW_OBJECT(MAKE_CLASS(cl)));
276 :     /* allocate and copy common slots */
277 :     /* FIXME: Use the triangular matrix classes for csn_LU */
278 :     SET_SLOT(ans, install("L"), /* these are free'd later if requested */
279 :     Matrix_cs_to_SEXP(N->L, "dgCMatrix", 0));
280 :     SET_SLOT(ans, install("U"),
281 :     Matrix_cs_to_SEXP(N->U, "dgCMatrix", 0));
282 :     switch(ctype) {
283 :     case 0:
284 :     Memcpy(INTEGER(ALLOC_SLOT(ans, install("Pinv"), INTSXP, n)),
285 : bates 1383 N->pinv, n);
286 : bates 1380 break;
287 :     case 1:
288 :     Memcpy(REAL(ALLOC_SLOT(ans, install("beta"), REALSXP, n)),
289 :     N->B, n);
290 :     break;
291 :     default:
292 : maechler 1906 error("Inappropriate class '%s' for Matrix_csn_to_SEXP", cl);
293 : bates 1380 }
294 :     if (dofree > 0) cs_nfree(N);
295 :     if (dofree < 0) {
296 :     Free(N->L); Free(N->U); Free(N);
297 :     }
298 :     UNPROTECT(1);
299 :     return ans;
300 :     }
301 :    
302 : bates 1383 #endif /* unused */

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