SCM

SCM Repository

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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : mmaechler 3186 #ifdef __GLIBC__
2 : mmaechler 3182 // to get strdup declared in glibc (when strict -std=c11 or -stdc99):
3 :     #define _POSIX_C_SOURCE 200809L
4 : mmaechler 3186 #endif
5 :    
6 : mmaechler 3182 #include <string.h>
7 :    
8 : bates 478 #include "dgCMatrix.h"
9 : bates 10
10 : maechler 1747 /* for Csparse_transpose() : */
11 :     #include "Csparse.h"
12 : bates 922 #include "chm_common.h"
13 : mmaechler 2327 /* -> Mutils.h / SPQR ... */
14 : bates 922
15 : maechler 956 /* FIXME -- we "forget" about dimnames almost everywhere : */
16 :    
17 : maechler 1661 /* for dgCMatrix _and_ lgCMatrix and others (but *not* ngC...) : */
18 :     SEXP xCMatrix_validate(SEXP x)
19 : bates 10 {
20 : maechler 1660 /* Almost everything now in Csparse_validate ( ./Csparse.c )
21 :     * *but* the checking of the 'x' slot : */
22 : mmaechler 3268 if (xlength(GET_SLOT(x, Matrix_iSym)) !=
23 :     xlength(GET_SLOT(x, Matrix_xSym)))
24 : maechler 1660 return mkString(_("lengths of slots 'i' and 'x' must match"));
25 : maechler 534
26 : bates 10 return ScalarLogical(1);
27 :     }
28 : maechler 534
29 : maechler 1968 /* for dgRMatrix _and_ lgRMatrix and others (but *not* ngC...) : */
30 :     SEXP xRMatrix_validate(SEXP x)
31 :     {
32 :     /* Almost everything now in Rsparse_validate ( ./Csparse.c )
33 :     * *but* the checking of the 'x' slot : */
34 : mmaechler 3268 if (xlength(GET_SLOT(x, Matrix_jSym)) !=
35 :     xlength(GET_SLOT(x, Matrix_xSym)))
36 : maechler 1968 return mkString(_("lengths of slots 'j' and 'x' must match"));
37 :    
38 :     return ScalarLogical(1);
39 :     }
40 :    
41 : maechler 1760 /* This and the following R_to_CMatrix() lead to memory-not-mapped seg.faults
42 :     * only with {32bit + R-devel + enable-R-shlib} -- no idea why */
43 : maechler 1747 SEXP compressed_to_TMatrix(SEXP x, SEXP colP)
44 : bates 10 {
45 : maechler 890 int col = asLogical(colP); /* 1 if "C"olumn compressed; 0 if "R"ow */
46 : maechler 1760 /* however, for Csparse, we now effectively use the cholmod-based
47 :     * Csparse_to_Tsparse() in ./Csparse.c ; maybe should simply write
48 :     * an as_cholmod_Rsparse() function and then do "as there" ...*/
49 : mmaechler 3213 SEXP indSym = col ? Matrix_iSym : Matrix_jSym, ans,
50 :     indP = PROTECT(GET_SLOT(x, indSym)),
51 :     pP = PROTECT(GET_SLOT(x, Matrix_pSym));
52 : bates 679 int npt = length(pP) - 1;
53 : bates 1867 char *ncl = strdup(class_P(x));
54 : mmaechler 2713 static const char *valid[] = { MATRIX_VALID_Csparse, MATRIX_VALID_Rsparse, ""};
55 : mmaechler 3147 int ctype = R_check_class_etc(x, valid);
56 : maechler 1747
57 :     if (ctype < 0)
58 : bates 1867 error(_("invalid class(x) '%s' in compressed_to_TMatrix(x)"), ncl);
59 : maechler 1747
60 : maechler 1766 /* replace 'C' or 'R' with 'T' :*/
61 : bates 1763 ncl[2] = 'T';
62 : mmaechler 3270 ans = PROTECT(NEW_OBJECT_OF_CLASS(ncl));
63 : maechler 1747
64 : maechler 1968 slot_dup(ans, x, Matrix_DimSym);
65 : maechler 1747 if((ctype / 3) % 4 != 2) /* not n..Matrix */
66 : maechler 1968 slot_dup(ans, x, Matrix_xSym);
67 : maechler 1747 if(ctype % 3) { /* s(ymmetric) or t(riangular) : */
68 : maechler 1968 slot_dup(ans, x, Matrix_uploSym);
69 : maechler 1747 if(ctype % 3 == 2) /* t(riangular) : */
70 : maechler 1968 slot_dup(ans, x, Matrix_diagSym);
71 : maechler 1747 }
72 :     SET_DimNames(ans, x);
73 : mmaechler 3055 // possibly asymmetric for symmetricMatrix is ok
74 : bates 677 SET_SLOT(ans, indSym, duplicate(indP));
75 : bates 679 expand_cmprPt(npt, INTEGER(pP),
76 :     INTEGER(ALLOC_SLOT(ans, col ? Matrix_jSym : Matrix_iSym,
77 :     INTSXP, length(indP))));
78 : bates 1763 free(ncl);
79 : mmaechler 3213 UNPROTECT(3);
80 : bates 10 return ans;
81 :     }
82 :    
83 : maechler 1747 SEXP R_to_CMatrix(SEXP x)
84 :     {
85 :     SEXP ans, tri = PROTECT(allocVector(LGLSXP, 1));
86 : bates 1867 char *ncl = strdup(class_P(x));
87 : mmaechler 2713 static const char *valid[] = { MATRIX_VALID_Rsparse, ""};
88 : mmaechler 3147 int ctype = R_check_class_etc(x, valid);
89 : maechler 1747 int *x_dims = INTEGER(GET_SLOT(x, Matrix_DimSym)), *a_dims;
90 : mmaechler 2455 PROTECT_INDEX ipx;
91 : maechler 1747
92 :     if (ctype < 0)
93 : bates 1867 error(_("invalid class(x) '%s' in R_to_CMatrix(x)"), ncl);
94 : maechler 1747
95 : maechler 1766 /* replace 'R' with 'C' : */
96 :     ncl[2] = 'C';
97 : mmaechler 3270 PROTECT_WITH_INDEX(ans = NEW_OBJECT_OF_CLASS(ncl), &ipx);
98 : maechler 1766
99 : maechler 1747 a_dims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));
100 :     /* reversed dim() since we will transpose: */
101 :     a_dims[0] = x_dims[1];
102 :     a_dims[1] = x_dims[0];
103 :    
104 :     /* triangular: */ LOGICAL(tri)[0] = 0;
105 :     if((ctype / 3) != 2) /* not n..Matrix */
106 : maechler 1968 slot_dup(ans, x, Matrix_xSym);
107 : maechler 1747 if(ctype % 3) { /* s(ymmetric) or t(riangular) : */
108 : maechler 1968 SET_SLOT(ans, Matrix_uploSym,
109 :     mkString((*uplo_P(x) == 'U') ? "L" : "U"));
110 : maechler 1747 if(ctype % 3 == 2) { /* t(riangular) : */
111 :     LOGICAL(tri)[0] = 1;
112 : maechler 1968 slot_dup(ans, x, Matrix_diagSym);
113 : maechler 1747 }
114 :     }
115 :     SET_SLOT(ans, Matrix_iSym, duplicate(GET_SLOT(x, Matrix_jSym)));
116 : maechler 1968 slot_dup(ans, x, Matrix_pSym);
117 : mmaechler 2455 REPROTECT(ans = Csparse_transpose(ans, tri), ipx);
118 : maechler 1747 SET_DimNames(ans, x);
119 : mmaechler 3055 // possibly asymmetric for symmetricMatrix is ok
120 : maechler 1771 free(ncl);
121 : maechler 1747 UNPROTECT(2);
122 :     return ans;
123 :     }
124 :    
125 : mmaechler 2814 /** Return a 2 column matrix '' cbind(i, j) '' of 0-origin index vectors (i,j)
126 :     * which entirely correspond to the (i,j) slots of
127 :     * as(x, "TsparseMatrix") :
128 :     */
129 : maechler 919 SEXP compressed_non_0_ij(SEXP x, SEXP colP)
130 :     {
131 :     int col = asLogical(colP); /* 1 if "C"olumn compressed; 0 if "R"ow */
132 :     SEXP ans, indSym = col ? Matrix_iSym : Matrix_jSym;
133 : mmaechler 3213 SEXP indP = PROTECT(GET_SLOT(x, indSym)),
134 :     pP = PROTECT(GET_SLOT(x, Matrix_pSym));
135 : mmaechler 2236 int i, *ij;
136 : dmbates 2759 int nouter = INTEGER(GET_SLOT(x, Matrix_DimSym))[col ? 1 : 0],
137 :     n_el = INTEGER(pP)[nouter]; /* is only == length(indP), if the
138 :     inner slot is not over-allocated */
139 : maechler 919
140 :     ij = INTEGER(ans = PROTECT(allocMatrix(INTSXP, n_el, 2)));
141 :     /* expand the compressed margin to 'i' or 'j' : */
142 : dmbates 2759 expand_cmprPt(nouter, INTEGER(pP), &ij[col ? n_el : 0]);
143 : maechler 919 /* and copy the other one: */
144 :     if (col)
145 :     for(i = 0; i < n_el; i++)
146 :     ij[i] = INTEGER(indP)[i];
147 :     else /* row compressed */
148 :     for(i = 0; i < n_el; i++)
149 :     ij[i + n_el] = INTEGER(indP)[i];
150 :    
151 : mmaechler 3213 UNPROTECT(3);
152 : maechler 919 return ans;
153 :     }
154 :    
155 : dmbates 2401 #if 0 /* unused */
156 : bates 1380 SEXP dgCMatrix_lusol(SEXP x, SEXP y)
157 : bates 332 {
158 : mmaechler 2237 SEXP ycp = PROTECT((TYPEOF(y) == REALSXP) ?
159 :     duplicate(y) : coerceVector(y, REALSXP));
160 : mmaechler 2227 CSP xc = AS_CSP__(x);
161 : maechler 1960 R_CheckStack();
162 : bates 310
163 : bates 1380 if (xc->m != xc->n || xc->m <= 0)
164 :     error(_("dgCMatrix_lusol requires a square, non-empty matrix"));
165 : mmaechler 2237 if (LENGTH(ycp) != xc->m)
166 : bates 1380 error(_("Dimensions of system to be solved are inconsistent"));
167 : bates 1383 if (!cs_lusol(/*order*/ 1, xc, REAL(ycp), /*tol*/ 1e-7))
168 : bates 1380 error(_("cs_lusol failed"));
169 : maechler 1960
170 : bates 1380 UNPROTECT(1);
171 :     return ycp;
172 :     }
173 : dmbates 2401 #endif
174 : bates 332
175 : mmaechler 2978 // called from package MatrixModels's R code
176 : mmaechler 2239 SEXP dgCMatrix_qrsol(SEXP x, SEXP y, SEXP ord)
177 : bates 1380 {
178 : mmaechler 2529 /* FIXME: extend this to work in multivariate case, i.e. y a matrix with > 1 column ! */
179 : mmaechler 2237 SEXP ycp = PROTECT((TYPEOF(y) == REALSXP) ?
180 :     duplicate(y) : coerceVector(y, REALSXP));
181 :     CSP xc = AS_CSP(x); /* <--> x may be dgC* or dtC* */
182 : mmaechler 3033 int order = asInteger(ord);
183 : mmaechler 2529 #ifdef _not_yet_do_FIXME__
184 :     const char *nms[] = {"L", "coef", "Xty", "resid", ""};
185 : mmaechler 2645 SEXP ans = PROTECT(Rf_mkNamed(VECSXP, nms));
186 : mmaechler 2529 #endif
187 : maechler 1960 R_CheckStack();
188 : bates 1380
189 : mmaechler 2239 if (order < 0 || order > 3)
190 : mmaechler 2529 error(_("dgCMatrix_qrsol(., order) needs order in {0,..,3}"));
191 : mmaechler 2239 /* --> cs_amd() --- order 0: natural, 1: Chol, 2: LU, 3: QR */
192 :     if (LENGTH(ycp) != xc->m)
193 :     error(_("Dimensions of system to be solved are inconsistent"));
194 :     /* FIXME? Note that qr_sol() would allow *under-determined systems;
195 : mmaechler 2529 * In general, we'd need LENGTH(ycp) = max(n,m)
196 :     * FIXME also: multivariate y (see above)
197 : mmaechler 2239 */
198 : bates 1380 if (xc->m < xc->n || xc->n <= 0)
199 : mmaechler 2237 error(_("dgCMatrix_qrsol(<%d x %d>-matrix) requires a 'tall' rectangular matrix"),
200 :     xc->m, xc->n);
201 : mmaechler 2239
202 :     /* cs_qrsol(): Tim Davis (2006) .. "8.2 Using a QR factorization", p.136f , calling
203 :     * ------- cs_sqr(order, ..), see p.76 */
204 : mmaechler 2529 /* MM: FIXME: write our *OWN* version of - the first case (m >= n) - of cs_qrsol()
205 :     * --------- which will (1) work with a *multivariate* y
206 :     * (2) compute coefficients properly, not overwriting RHS
207 :     */
208 : mmaechler 2239 if (!cs_qrsol(order, xc, REAL(ycp)))
209 :     /* return value really is 0 or 1 - no more info there */
210 : mmaechler 2529 error(_("cs_qrsol() failed inside dgCMatrix_qrsol()"));
211 : maechler 1960
212 : mmaechler 2239 /* Solution is only in the first part of ycp -- cut its length back to n : */
213 : mmaechler 3161 ycp = lengthgets(ycp, (R_xlen_t) xc->n);
214 : mmaechler 2879
215 : bates 332 UNPROTECT(1);
216 : bates 1380 return ycp;
217 : bates 332 }
218 : bates 1380
219 : mmaechler 2854 // Modified version of Tim Davis's cs_qr_mex.c file for MATLAB (in CSparse)
220 :     // Usage: [V,beta,p,R,q] = cs_qr(A) ;
221 : mmaechler 3115 SEXP dgCMatrix_QR(SEXP Ap, SEXP order, SEXP keep_dimnames)
222 : bates 1380 {
223 : mmaechler 2227 CSP A = AS_CSP__(Ap), D;
224 : mmaechler 2854 int io = INTEGER(order)[0];
225 : mmaechler 2933 Rboolean verbose = (io < 0);// verbose=TRUE, encoded with negative 'order'
226 : mmaechler 3115 int m0 = A->m, m = m0, n = A->n, ord = asLogical(order) ? 3 : 0, *p;
227 : maechler 1960 R_CheckStack();
228 : bates 1380
229 : mmaechler 2387 if (m < n) error(_("A must have #{rows} >= #{columns}")) ;
230 : mmaechler 3273 SEXP ans = PROTECT(NEW_OBJECT_OF_CLASS("sparseQR"));
231 : mmaechler 2854 int *dims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));
232 : maechler 2106 dims[0] = m; dims[1] = n;
233 : mmaechler 2854 css *S = cs_sqr(ord, A, 1); /* symbolic QR ordering & analysis*/
234 : mmaechler 2387 if (!S) error(_("cs_sqr failed"));
235 : mmaechler 3115 int keep_dimnms = asLogical(keep_dimnames);
236 :     if(keep_dimnms == NA_LOGICAL) { keep_dimnms = TRUE;
237 :     warning(_("dgcMatrix_QR(*, keep_dimnames = NA): NA taken as TRUE"));
238 :     }
239 : mmaechler 2854 if(verbose && S->m2 > m) // in ./cs.h , m2 := # of rows for QR, after adding fictitious rows
240 :     Rprintf("Symbolic QR(): Matrix structurally rank deficient (m2-m = %d)\n",
241 :     S->m2 - m);
242 :     csn *N = cs_qr(A, S); /* numeric QR factorization */
243 : mmaechler 2387 if (!N) error(_("cs_qr failed")) ;
244 : bates 1381 cs_dropzeros(N->L); /* drop zeros from V and sort */
245 : mmaechler 2854 D = cs_transpose(N->L, 1); cs_spfree(N->L);
246 :     N->L = cs_transpose(D, 1); cs_spfree(D);
247 : bates 1381 cs_dropzeros(N->U); /* drop zeros from R and sort */
248 : mmaechler 2854 D = cs_transpose(N->U, 1); cs_spfree(N->U) ;
249 :     N->U = cs_transpose(D, 1); cs_spfree(D);
250 : bates 1381 m = N->L->m; /* m may be larger now */
251 : mmaechler 2854 // MM: m := S->m2 also counting the ficticious rows (Tim Davis, p.72, 74f)
252 : bates 1383 p = cs_pinv(S->pinv, m); /* p = pinv' */
253 : mmaechler 3115 SEXP dn = R_NilValue; Rboolean do_dn = FALSE;
254 :     if(keep_dimnms) {
255 :     dn = GET_SLOT(Ap, Matrix_DimNamesSym);
256 :     do_dn = !isNull(VECTOR_ELT(dn, 0)) && m == m0;
257 :     // FIXME? also deal with case m > m0 ?
258 :     if(do_dn) { // keep rownames
259 :     dn = PROTECT(duplicate(dn));
260 :     SET_VECTOR_ELT(dn, 1, R_NilValue);
261 :     } else dn = R_NilValue;
262 :     }
263 :     SET_SLOT(ans, Matrix_VSym, Matrix_cs_to_SEXP(N->L, "dgCMatrix", 0, dn)); // "V"
264 :     Memcpy(REAL(ALLOC_SLOT(ans, Matrix_betaSym, REALSXP, n)), N->B, n);
265 :     Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, m)), p, m);
266 :     if(do_dn) {
267 :     UNPROTECT(1); // dn
268 :     dn = R_NilValue; do_dn = FALSE;
269 :     }
270 :     if (ord) {
271 :     Memcpy(INTEGER(ALLOC_SLOT(ans, install("q"), INTSXP, n)), S->q, n);
272 :     if(keep_dimnms) {
273 :     dn = GET_SLOT(Ap, Matrix_DimNamesSym);
274 :     do_dn = !isNull(VECTOR_ELT(dn, 1));
275 :     if(do_dn) {
276 :     dn = PROTECT(duplicate(dn));
277 :     // permute colnames by S->q : cn <- cn[ S->q ] :
278 :     SEXP cns = PROTECT(duplicate(VECTOR_ELT(dn, 1)));
279 :     for(int j=0; j < n; j++)
280 :     SET_STRING_ELT(VECTOR_ELT(dn, 1), j, STRING_ELT(cns, S->q[j]));
281 :     UNPROTECT(1);
282 :     SET_VECTOR_ELT(dn, 0, R_NilValue);
283 :     } else dn = R_NilValue;
284 :     }
285 :     } else
286 : bates 1397 ALLOC_SLOT(ans, install("q"), INTSXP, 0);
287 : mmaechler 3115 SET_SLOT(ans, install("R"), Matrix_cs_to_SEXP(N->U, "dgCMatrix", 0, dn));
288 :     if(do_dn) UNPROTECT(1); // dn
289 : bates 1380 cs_nfree(N);
290 :     cs_sfree(S);
291 :     cs_free(p);
292 :     UNPROTECT(1);
293 :     return ans;
294 :     }
295 :    
296 : mmaechler 2336 #ifdef Matrix_with_SPQR
297 : dmbates 2276 /**
298 :     * Return a SuiteSparse QR factorization of the sparse matrix A
299 :     *
300 : mmaechler 2322 * @param Ap (pointer to) a [m x n] dgCMatrix
301 :     * @param ordering integer SEXP specifying the ordering strategy to be used
302 :     * see SPQR/Include/SuiteSparseQR_definitions.h
303 :     * @param econ integer SEXP ("economy"): number of rows of R and columns of Q
304 :     * to return. The default is m. Using n gives the standard economy form.
305 :     * A value less than the estimated rank r is set to r, so econ=0 gives the
306 :     * "rank-sized" factorization, where nrow(R)==nnz(diag(R))==r.
307 : mmaechler 2327 * @param tol double SEXP: if tol <= -2 use SPQR's default,
308 :     * if -2 < tol < 0, then no tol is used; otherwise,
309 :     * tol > 0, use as tolerance: columns with 2-norm <= tol treated as 0
310 : dmbates 2276 *
311 : mmaechler 2327 *
312 :     * @return SEXP "SPQR" object with slots (Q, R, p, rank, Dim):
313 : mmaechler 2322 * Q: dgCMatrix; R: dgCMatrix [subject to change to dtCMatrix FIXME ?]
314 : mmaechler 2327 * p: integer: 0-based permutation (or length 0 <=> identity);
315 :     * rank: integer, the "revealed" rank Dim: integer, original matrix dim.
316 : dmbates 2276 */
317 : mmaechler 2322 SEXP dgCMatrix_SPQR(SEXP Ap, SEXP ordering, SEXP econ, SEXP tol)
318 : dmbates 2276 {
319 : mmaechler 2327 /* SEXP ans = PROTECT(allocVector(VECSXP, 4)); */
320 : mmaechler 3273 SEXP ans = PROTECT(NEW_OBJECT_OF_CLASS("SPQR"));
321 : mmaechler 2327
322 : dmbates 2299 CHM_SP A = AS_CHM_SP(Ap), Q, R;
323 : mmaechler 3001 SuiteSparse_long *E, rank;/* not always = int FIXME (Windows_64 ?) */
324 : mmaechler 2322
325 :     if ((rank = SuiteSparseQR_C_QR(asInteger(ordering),
326 :     asReal(tol),/* originally had SPQR_DEFAULT_TOL */
327 : mmaechler 3001 (SuiteSparse_long)asInteger(econ),/* originally had 0 */
328 : mmaechler 2661 A, &Q, &R, &E, &cl)) == -1)
329 : dmbates 2276 error(_("SuiteSparseQR_C_QR returned an error code"));
330 : mmaechler 2327
331 : mmaechler 3011 slot_dup(ans, Ap, Matrix_DimSym);
332 : mmaechler 2327 /* SET_VECTOR_ELT(ans, 0, */
333 :     /* chm_sparse_to_SEXP(Q, 0, 0, 0, "", R_NilValue)); */
334 :     SET_SLOT(ans, install("Q"),
335 :     chm_sparse_to_SEXP(Q, 0, 0, 0, "", R_NilValue));
336 :    
337 : mmaechler 2322 /* Also gives a dgCMatrix (not a dtC* *triangular*) :
338 :     * may make sense if to be used in the "spqr_solve" routines .. ?? */
339 : mmaechler 2327 /* SET_VECTOR_ELT(ans, 1, */
340 :     /* chm_sparse_to_SEXP(R, 0, 0, 0, "", R_NilValue)); */
341 :     SET_SLOT(ans, install("R"),
342 :     chm_sparse_to_SEXP(R, 0, 0, 0, "", R_NilValue));
343 : mmaechler 2661 cholmod_free_sparse(&Al, &cl);
344 :     cholmod_free_sparse(&R, &cl);
345 :     cholmod_free_sparse(&Q, &cl);
346 : dmbates 2276 if (E) {
347 : mmaechler 2661 int *Er;
348 :     SET_VECTOR_ELT(ans, 2, allocVector(INTSXP, A->ncol));
349 :     Er = INTEGER(VECTOR_ELT(ans, 2));
350 :     for (int i = 0; i < A->ncol; i++) Er[i] = (int) E[i];
351 : dmbates 2276 Free(E);
352 : mmaechler 2661 } else SET_VECTOR_ELT(ans, 2, allocVector(INTSXP, 0));
353 :     SET_VECTOR_ELT(ans, 3, ScalarInteger((int)rank));
354 : dmbates 2276 UNPROTECT(1);
355 :     return ans;
356 :     }
357 : mmaechler 2338 #endif
358 :     /* Matrix_with_SPQR */
359 : dmbates 2276
360 : bates 1463 /* Modified version of Tim Davis's cs_lu_mex.c file for MATLAB */
361 : mmaechler 3115 void install_lu(SEXP Ap, int order, double tol, Rboolean err_sing, Rboolean keep_dimnms)
362 : bates 1383 {
363 : mmaechler 2713 // (order, tol) == (1, 1) by default, when called from R.
364 : dmbates 2401 SEXP ans;
365 : bates 1383 css *S;
366 :     csn *N;
367 : dmbates 2401 int n, *p, *dims;
368 :     CSP A = AS_CSP__(Ap), D;
369 : maechler 1960 R_CheckStack();
370 : bates 1383
371 : bates 1385 n = A->n;
372 : bates 1383 if (A->m != n)
373 : mmaechler 2387 error(_("LU decomposition applies only to square matrices"));
374 : bates 1383 if (order) { /* not using natural order */
375 :     order = (tol == 1) ? 2 /* amd(S'*S) w/dense rows or I */
376 :     : 1; /* amd (A+A'), or natural */
377 :     }
378 : mmaechler 2713 S = cs_sqr(order, A, /*qr = */ 0); /* symbolic ordering */
379 : dmbates 2401 N = cs_lu(A, S, tol); /* numeric factorization */
380 : mmaechler 2488 if (!N) {
381 :     if(err_sing)
382 :     error(_("cs_lu(A) failed: near-singular A (or out of memory)"));
383 :     else {
384 :     /* No warning: The useR should be careful :
385 : mmaechler 2692 * Put NA into "LU" factor */
386 : mmaechler 2488 set_factors(Ap, ScalarLogical(NA_LOGICAL), "LU");
387 :     return;
388 :     }
389 :     }
390 : dmbates 2401 cs_dropzeros(N->L); /* drop zeros from L and sort it */
391 :     D = cs_transpose(N->L, 1);
392 :     cs_spfree(N->L);
393 :     N->L = cs_transpose(D, 1);
394 :     cs_spfree(D);
395 :     cs_dropzeros(N->U); /* drop zeros from U and sort it */
396 :     D = cs_transpose(N->U, 1);
397 :     cs_spfree(N->U);
398 :     N->U = cs_transpose(D, 1);
399 :     cs_spfree(D);
400 :     p = cs_pinv(N->pinv, n); /* p=pinv' */
401 : mmaechler 3273 ans = PROTECT(NEW_OBJECT_OF_CLASS("sparseLU"));
402 : mmaechler 2268 dims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));
403 :     dims[0] = n; dims[1] = n;
404 : mmaechler 3115 SEXP dn; Rboolean do_dn = FALSE;
405 :     if(keep_dimnms) {
406 :     dn = GET_SLOT(Ap, Matrix_DimNamesSym);
407 :     do_dn = !isNull(VECTOR_ELT(dn, 0));
408 :     if(do_dn) {
409 :     dn = PROTECT(duplicate(dn));
410 :     // permute rownames by p : rn <- rn[ p ] :
411 :     SEXP rn = PROTECT(duplicate(VECTOR_ELT(dn, 0)));
412 :     for(int i=0; i < n; i++)
413 :     SET_STRING_ELT(VECTOR_ELT(dn, 0), i, STRING_ELT(rn, p[i]));
414 :     UNPROTECT(1); // rn
415 :     SET_VECTOR_ELT(dn, 1, R_NilValue); // colnames(.) := NULL
416 :     }
417 :     }
418 : bates 1383 SET_SLOT(ans, install("L"),
419 : mmaechler 3115 Matrix_cs_to_SEXP(N->L, "dtCMatrix", 0, do_dn ? dn : R_NilValue));
420 :    
421 :     if(keep_dimnms) {
422 :     if(do_dn) {
423 :     UNPROTECT(1); // dn
424 :     dn = GET_SLOT(Ap, Matrix_DimNamesSym);
425 :     }
426 :     do_dn = !isNull(VECTOR_ELT(dn, 1));
427 :     if(do_dn) {
428 :     dn = PROTECT(duplicate(dn));
429 :     if(order) { // permute colnames by S->q : cn <- cn[ S->q ] :
430 :     SEXP cn = PROTECT(duplicate(VECTOR_ELT(dn, 1)));
431 :     for(int j=0; j < n; j++)
432 :     SET_STRING_ELT(VECTOR_ELT(dn, 1), j, STRING_ELT(cn, S->q[j]));
433 :     UNPROTECT(1); // cn
434 :     }
435 :     SET_VECTOR_ELT(dn, 0, R_NilValue); // rownames(.) := NULL
436 :     }
437 :     }
438 : bates 1383 SET_SLOT(ans, install("U"),
439 : mmaechler 3115 Matrix_cs_to_SEXP(N->U, "dtCMatrix", 0, do_dn ? dn : R_NilValue));
440 :     if(do_dn) UNPROTECT(1); // dn
441 : mmaechler 2247 Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_pSym, /* "p" */
442 : bates 1383 INTSXP, n)), p, n);
443 :     if (order)
444 :     Memcpy(INTEGER(ALLOC_SLOT(ans, install("q"),
445 :     INTSXP, n)), S->q, n);
446 :     cs_nfree(N);
447 :     cs_sfree(S);
448 :     cs_free(p);
449 :     UNPROTECT(1);
450 : dmbates 2401 set_factors(Ap, ans, "LU");
451 : bates 1385 }
452 :    
453 : mmaechler 3115 SEXP dgCMatrix_LU(SEXP Ap, SEXP orderp, SEXP tolp, SEXP error_on_sing, SEXP keep_dimnames)
454 : dmbates 2401 {
455 :     SEXP ans;
456 : mmaechler 2488 Rboolean err_sing = asLogical(error_on_sing);
457 : dmbates 2401 /* FIXME: dgCMatrix_LU should check ans for consistency in
458 :     * permutation type with the requested value - Should have two
459 :     * classes or two different names in the factors list for LU with
460 :     * permuted columns or not.
461 :     * OTOH, currently (order, tol) === (1, 1) always.
462 : dmbates 2403 * It is true that length(LU@q) does flag the order argument.
463 : dmbates 2401 */
464 :     if (!isNull(ans = get_factors(Ap, "LU")))
465 :     return ans;
466 : mmaechler 3115 int keep_dimnms = asLogical(keep_dimnames);
467 :     if(keep_dimnms == NA_LOGICAL) { keep_dimnms = TRUE;
468 :     warning(_("dgcMatrix_LU(*, keep_dimnames = NA): NA taken as TRUE"));
469 :     }
470 :     install_lu(Ap, asInteger(orderp), asReal(tolp), err_sing, keep_dimnms);
471 : dmbates 2401 return get_factors(Ap, "LU");
472 :     }
473 :    
474 : mmaechler 2889 SEXP dgCMatrix_matrix_solve(SEXP Ap, SEXP b, SEXP give_sparse)
475 : mmaechler 3115 // FIXME: add 'keep_dimnames' as argument
476 : bates 1385 {
477 : mmaechler 2889 Rboolean sparse = asLogical(give_sparse);
478 :     if(sparse) {
479 :     // FIXME: implement this
480 :     error(_("dgCMatrix_matrix_solve(.., sparse=TRUE) not yet implemented"));
481 :    
482 :     /* Idea: in the for(j = 0; j < nrhs ..) loop below, build the *sparse* result matrix
483 :     * ----- *column* wise -- which is perfect for dgCMatrix
484 :     * --> build (i,p,x) slots "increasingly" [well, allocate in batches ..]
485 :     *
486 :     * --> maybe first a protoype in R
487 :     */
488 :    
489 :     }
490 : dmbates 2401 SEXP ans = PROTECT(dup_mMatrix_as_dgeMatrix(b)),
491 :     lu, qslot;
492 :     CSP L, U;
493 :     int *bdims = INTEGER(GET_SLOT(ans, Matrix_DimSym)), *p, *q;
494 : bates 1385 int j, n = bdims[0], nrhs = bdims[1];
495 : mmaechler 3105 double *x, *ax = REAL(GET_SLOT(ans, Matrix_xSym));
496 :     C_or_Alloca_TO(x, n, double);
497 : bates 1385
498 : dmbates 2401 if (isNull(lu = get_factors(Ap, "LU"))) {
499 : mmaechler 3115 install_lu(Ap, /* order = */ 1, /* tol = */ 1.0, /* err_sing = */ TRUE,
500 :     /* keep_dimnames = */ TRUE);
501 : dmbates 2401 lu = get_factors(Ap, "LU");
502 :     }
503 :     qslot = GET_SLOT(lu, install("q"));
504 :     L = AS_CSP__(GET_SLOT(lu, install("L")));
505 :     U = AS_CSP__(GET_SLOT(lu, install("U")));
506 :     R_CheckStack();
507 : mmaechler 3011 if (U->n != n)
508 :     error(_("Dimensions of system to be solved are inconsistent"));
509 :     if(nrhs >= 1 && n >= 1) {
510 :     p = INTEGER(GET_SLOT(lu, Matrix_pSym));
511 :     q = LENGTH(qslot) ? INTEGER(qslot) : (int *) NULL;
512 : dmbates 2401
513 : mmaechler 3011 for (j = 0; j < nrhs; j++) {
514 :     cs_pvec(p, ax + j * n, x, n); /* x = b(p) */
515 :     cs_lsolve(L, x); /* x = L\x */
516 :     cs_usolve(U, x); /* x = U\x */
517 : mmaechler 3023 if (q) /* r(q) = x , hence
518 :     r = Q' U{^-1} L{^-1} P b = A^{-1} b */
519 : mmaechler 3011 cs_ipvec(q, x, ax + j * n, n);
520 :     else
521 :     Memcpy(ax + j * n, x, n);
522 :     }
523 : bates 1385 }
524 : mmaechler 3105 if(n >= SMALL_4_Alloca) Free(x);
525 : dmbates 2401 UNPROTECT(1);
526 : bates 1383 return ans;
527 :     }
528 : bates 1385
529 : mmaechler 2978 // called from package MatrixModels's R code:
530 : bates 1895 SEXP dgCMatrix_cholsol(SEXP x, SEXP y)
531 :     {
532 : mmaechler 2239 /* Solve Sparse Least Squares X %*% beta ~= y with dense RHS y,
533 :     * where X = t(x) i.e. we pass x = t(X) as argument,
534 :     * via "Cholesky(X'X)" .. well not really:
535 :     * cholmod_factorize("x", ..) finds L in X'X = L'L directly */
536 : mmaechler 2237 CHM_SP cx = AS_CHM_SP(x);
537 : mmaechler 2529 /* FIXME: extend this to work in multivariate case, i.e. y a matrix with > 1 column ! */
538 : mmaechler 3276 SEXP y_ = PROTECT(coerceVector(y, REALSXP));
539 :     CHM_DN cy = AS_CHM_DN(y_), rhs, cAns, resid;
540 : maechler 1960 CHM_FR L;
541 : mmaechler 2529 int n = cx->ncol;/* #{obs.} {x = t(X) !} */
542 :     double one[] = {1,0}, zero[] = {0,0}, neg1[] = {-1,0};
543 :     const char *nms[] = {"L", "coef", "Xty", "resid", ""};
544 : mmaechler 2645 SEXP ans = PROTECT(Rf_mkNamed(VECSXP, nms));
545 : maechler 1960 R_CheckStack();
546 : bates 1895
547 : mmaechler 2529 if (n < cx->nrow || n <= 0)
548 : bates 1895 error(_("dgCMatrix_cholsol requires a 'short, wide' rectangular matrix"));
549 : mmaechler 2529 if (cy->nrow != n)
550 : bates 1895 error(_("Dimensions of system to be solved are inconsistent"));
551 : mmaechler 2661 rhs = cholmod_allocate_dense(cx->nrow, 1, cx->nrow, CHOLMOD_REAL, &c);
552 : mmaechler 2529 /* cholmod_sdmult(A, transp, alpha, beta, X, Y, &c):
553 :     * Y := alpha*(A*X) + beta*Y or alpha*(A'*X) + beta*Y ;
554 :     * here: rhs := 1 * x %*% y + 0 = x %*% y = X'y */
555 : mmaechler 2661 if (!(cholmod_sdmult(cx, 0 /* trans */, one, zero, cy, rhs, &c)))
556 :     error(_("cholmod_sdmult error (rhs)"));
557 :     L = cholmod_analyze(cx, &c);
558 :     if (!cholmod_factorize(cx, L, &c))
559 :     error(_("cholmod_factorize failed: status %d, minor %d from ncol %d"),
560 : bates 1895 c.status, L->minor, L->n);
561 :     /* FIXME: Do this in stages so an "effects" vector can be calculated */
562 : mmaechler 2661 if (!(cAns = cholmod_solve(CHOLMOD_A, L, rhs, &c)))
563 :     error(_("cholmod_solve (CHOLMOD_A) failed: status %d, minor %d from ncol %d"),
564 : bates 1895 c.status, L->minor, L->n);
565 : mmaechler 2529 /* L : */
566 : bates 1895 SET_VECTOR_ELT(ans, 0, chm_factor_to_SEXP(L, 0));
567 : mmaechler 2529 /* coef : */
568 : bates 1895 SET_VECTOR_ELT(ans, 1, allocVector(REALSXP, cx->nrow));
569 :     Memcpy(REAL(VECTOR_ELT(ans, 1)), (double*)(cAns->x), cx->nrow);
570 : mmaechler 2529 /* X'y : */
571 : bates 1895 /* FIXME: Change this when the "effects" vector is available */
572 :     SET_VECTOR_ELT(ans, 2, allocVector(REALSXP, cx->nrow));
573 : mmaechler 2527 Memcpy(REAL(VECTOR_ELT(ans, 2)), (double*)(rhs->x), cx->nrow);
574 : mmaechler 2529 /* resid := y */
575 : mmaechler 2661 resid = cholmod_copy_dense(cy, &c);
576 : mmaechler 2529 /* cholmod_sdmult(A, transp, alp, bet, X, Y, &c):
577 :     * Y := alp*(A*X) + bet*Y or alp*(A'*X) + beta*Y ;
578 :     * here: resid := -1 * x' %*% coef + 1 * y = y - X %*% coef */
579 : mmaechler 2661 if (!(cholmod_sdmult(cx, 1/* trans */, neg1, one, cAns, resid, &c)))
580 :     error(_("cholmod_sdmult error (resid)"));
581 : mmaechler 2529 /* FIXME: for multivariate case, i.e. resid *matrix* with > 1 column ! */
582 :     SET_VECTOR_ELT(ans, 3, allocVector(REALSXP, n));
583 :     Memcpy(REAL(VECTOR_ELT(ans, 3)), (double*)(resid->x), n);
584 : bates 1895
585 : mmaechler 2661 cholmod_free_factor(&L, &c);
586 :     cholmod_free_dense(&rhs, &c);
587 :     cholmod_free_dense(&cAns, &c);
588 : mmaechler 3276 UNPROTECT(2);
589 : bates 1895 return ans;
590 :     }
591 : bates 1901
592 :    
593 : maechler 2015 /* Define all of
594 : maechler 1911 * dgCMatrix_colSums(....)
595 :     * igCMatrix_colSums(....)
596 : maechler 2015 * lgCMatrix_colSums_d(....)
597 :     * lgCMatrix_colSums_i(....)
598 :     * ngCMatrix_colSums_d(....)
599 :     * ngCMatrix_colSums_i(....)
600 : maechler 1911 */
601 : maechler 1920 #define _dgC_
602 :     #include "t_gCMatrix_colSums.c"
603 : bates 1901
604 : maechler 1920 #define _igC_
605 :     #include "t_gCMatrix_colSums.c"
606 : bates 1901
607 : maechler 1920 #define _lgC_
608 :     #include "t_gCMatrix_colSums.c"
609 : maechler 1911
610 : maechler 1920 #define _ngC_
611 :     #include "t_gCMatrix_colSums.c"
612 : maechler 1911
613 : maechler 1920 #define _lgC_mn
614 :     #include "t_gCMatrix_colSums.c"
615 : maechler 1911
616 : maechler 1920 #define _ngC_mn
617 :     #include "t_gCMatrix_colSums.c"
618 : maechler 1911
619 :    
620 : maechler 1920 SEXP lgCMatrix_colSums(SEXP x, SEXP NArm, SEXP spRes, SEXP trans, SEXP means)
621 :     {
622 :     if(asLogical(means)) /* ==> result will be "double" / "dsparseVector" */
623 :     return lgCMatrix_colSums_d(x, NArm, spRes, trans, means);
624 :     else
625 :     return lgCMatrix_colSums_i(x, NArm, spRes, trans, means);
626 :     }
627 : maechler 1911
628 : maechler 1920 SEXP ngCMatrix_colSums(SEXP x, SEXP NArm, SEXP spRes, SEXP trans, SEXP means)
629 :     {
630 :     if(asLogical(means)) /* ==> result will be "double" / "dsparseVector" */
631 :     return ngCMatrix_colSums_d(x, NArm, spRes, trans, means);
632 :     else
633 :     return ngCMatrix_colSums_i(x, NArm, spRes, trans, means);
634 :     }

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