Matrix r5059
Loading...
Searching...
No Matches
solve.c
Go to the documentation of this file.
1/* C implementation of methods for solve */
2
3#include "Lapack-etc.h"
4#include "cs-etc.h"
5#include "cholmod-etc.h"
6#include "Mdefines.h"
7#include "M5.h"
8#include "idz.h"
9
10static
11void solveDN(SEXP rdn, SEXP adn, SEXP bdn)
12{
13 SEXP s;
14 if ((s = VECTOR_ELT(adn, 1)) != R_NilValue)
15 SET_VECTOR_ELT(rdn, 0, s);
16 if ((s = VECTOR_ELT(bdn, 1)) != R_NilValue)
17 SET_VECTOR_ELT(rdn, 1, s);
18 PROTECT(adn = Rf_getAttrib(adn, R_NamesSymbol));
19 PROTECT(bdn = Rf_getAttrib(bdn, R_NamesSymbol));
20 if (adn != R_NilValue || bdn != R_NilValue) {
21 PROTECT(s = Rf_allocVector(STRSXP, 2));
22 if (adn != R_NilValue)
23 SET_STRING_ELT(s, 0, STRING_ELT(adn, 1));
24 if (bdn != R_NilValue)
25 SET_STRING_ELT(s, 1, STRING_ELT(bdn, 1));
26 Rf_setAttrib(rdn, R_NamesSymbol, s);
27 UNPROTECT(1);
28 }
29 UNPROTECT(2);
30 return;
31}
32
33SEXP denseLU_solve(SEXP s_a, SEXP s_b)
34{
35
36#define SOLVE_START \
37 int *padim = DIM(s_a), m = padim[0], n = padim[1]; \
38 if (m != n) \
39 Rf_error(_("'%s' is not square"), "a"); \
40 if (s_b != R_NilValue) { \
41 int *pbdim = DIM(s_b); \
42 if (pbdim[0] != m) \
43 Rf_error(_("dimensions of '%s' and '%s' are inconsistent"), \
44 "a", "b"); \
45 n = pbdim[1]; \
46 }
47
48#define SOLVE_FINISH \
49 SEXP rdimnames = PROTECT(DIMNAMES(r, 0)), \
50 adimnames = PROTECT(DIMNAMES(s_a, 0)); \
51 if (s_b == R_NilValue) \
52 cpyDN(rdimnames, adimnames, 1); \
53 else { \
54 SEXP bdimnames = PROTECT(DIMNAMES(s_b, 0)); \
55 solveDN(rdimnames, adimnames, bdimnames); \
56 UNPROTECT(1); /* bdimnames */ \
57 } \
58 UNPROTECT(2); /* adimnames, rdimnames */
59
60 SOLVE_START;
61
62 SEXP ax = GET_SLOT(s_a, Matrix_xSym), bx = R_NilValue, obx = bx;
63 if (s_b != R_NilValue) {
64 PROTECT(ax);
65 PROTECT(bx = obx = GET_SLOT(s_b, Matrix_xSym));
66 if (TYPEOF(ax) != TYPEOF(bx)) {
67 ax = Rf_coerceVector(ax, CPLXSXP);
68 bx = Rf_coerceVector(bx, CPLXSXP);
69 }
70 UNPROTECT(2);
71 }
72 PROTECT(ax);
73 PROTECT(bx);
74
75 char rcl[] = ".geMatrix";
76 rcl[0] = (TYPEOF(ax) == CPLXSXP) ? 'z' : 'd';
77 SEXP r = PROTECT(newObject(rcl));
78
79 SET_DIM(r, m, n);
80
81 if (m > 0) {
82 SEXP apivot = PROTECT(GET_SLOT(s_a, Matrix_permSym)), rx;
83 int info;
84 if (s_b == R_NilValue) {
85 PROTECT(rx = duplicateVector(ax));
86 int lwork = -1;
87 if (TYPEOF(ax) == CPLXSXP) {
88 Rcomplex work0, *work = &work0;
89 F77_CALL(zgetri)(&m, COMPLEX(rx), &m, INTEGER(apivot),
90 work, &lwork, &info);
91 ERROR_LAPACK_1(zgetri, info);
92 lwork = (int) work0.r;
93 work = (Rcomplex *) R_alloc((size_t) lwork, sizeof(Rcomplex));
94 F77_CALL(zgetri)(&m, COMPLEX(rx), &m, INTEGER(apivot),
95 work, &lwork, &info);
96 ERROR_LAPACK_2(zgetri, info, 2, U);
97 } else {
98 double work0, *work = &work0;
99 F77_CALL(dgetri)(&m, REAL(rx), &m, INTEGER(apivot),
100 work, &lwork, &info);
101 ERROR_LAPACK_1(dgetri, info);
102 lwork = (int) work0;
103 work = (double *) R_alloc((size_t) lwork, sizeof(double ));
104 F77_CALL(dgetri)(&m, REAL(rx), &m, INTEGER(apivot),
105 work, &lwork, &info);
106 ERROR_LAPACK_2(dgetri, info, 2, U);
107 }
108 } else {
109 PROTECT(rx = (bx == obx) ? duplicateVector(bx) : bx);
110 if (TYPEOF(ax) == CPLXSXP) {
111 F77_CALL(zgetrs)("N", &m, &n, COMPLEX(ax), &m, INTEGER(apivot),
112 COMPLEX(rx), &m, &info FCONE);
113 ERROR_LAPACK_1(zgetrs, info);
114 } else {
115 F77_CALL(dgetrs)("N", &m, &n, REAL(ax), &m, INTEGER(apivot),
116 REAL(rx), &m, &info FCONE);
117 ERROR_LAPACK_1(dgetrs, info);
118 }
119 }
120 SET_SLOT(r, Matrix_xSym, rx);
121 UNPROTECT(2); /* rx, apivot */
122 }
123
124 SOLVE_FINISH;
125
126 UNPROTECT(3); /* r, bx, ax */
127 return r;
128}
129
130SEXP denseBunchKaufman_solve(SEXP s_a, SEXP s_b)
131{
132 SOLVE_START;
133
134 SEXP ax = GET_SLOT(s_a, Matrix_xSym), oax = ax,
135 bx = R_NilValue, obx = bx;
136 if (s_b != R_NilValue) {
137 PROTECT(ax);
138 PROTECT(bx = obx = GET_SLOT(s_b, Matrix_xSym));
139 if (TYPEOF(ax) != TYPEOF(bx)) {
140 ax = Rf_coerceVector(ax, CPLXSXP);
141 bx = Rf_coerceVector(bx, CPLXSXP);
142 }
143 UNPROTECT(2);
144 }
145 PROTECT(ax);
146 PROTECT(bx);
147
148 int packed = XLENGTH(ax) != (int_fast64_t) m * m;
149
150 char rcl[] = "...Matrix";
151 rcl[0] = (TYPEOF(ax) == CPLXSXP) ? 'z' : 'd';
152 if (s_b == R_NilValue) {
153 rcl[1] = 's';
154 rcl[2] = (!packed) ? 'y' : 'p';
155 } else {
156 rcl[1] = 'g';
157 rcl[2] = 'e';
158 }
159 SEXP r = PROTECT(newObject(rcl));
160
161 SET_DIM(r, m, n);
162
163 char aul = UPLO(s_a);
164 if (s_b == R_NilValue && aul != 'U')
165 SET_UPLO(r);
166
167 char act = (TYPEOF(ax) == CPLXSXP && ax == oax) ? TRANS(s_a) : 'C';
168 if (s_b == R_NilValue && act != 'C')
169 SET_TRANS(r);
170
171 if (m > 0) {
172 SEXP apivot = PROTECT(GET_SLOT(s_a, Matrix_permSym)), rx;
173 int info;
174 if (s_b == R_NilValue) {
175 PROTECT(rx = duplicateVector(ax));
176 if (TYPEOF(ax) == CPLXSXP) {
177 Rcomplex *work = (Rcomplex *) R_alloc((size_t) m, sizeof(Rcomplex));
178 if (act == 'C') {
179 if (!packed) {
180 F77_CALL(zhetri)(&aul, &m, COMPLEX(rx), &m, INTEGER(apivot),
181 work, &info FCONE);
182 ERROR_LAPACK_2(zhetri, info, 2, D);
183 } else {
184 F77_CALL(zhptri)(&aul, &m, COMPLEX(rx), INTEGER(apivot),
185 work, &info FCONE);
186 ERROR_LAPACK_2(zhptri, info, 2, D);
187 }
188 } else {
189 if (!packed) {
190 F77_CALL(zsytri)(&aul, &m, COMPLEX(rx), &m, INTEGER(apivot),
191 work, &info FCONE);
192 ERROR_LAPACK_2(zsytri, info, 2, D);
193 } else {
194 F77_CALL(zsptri)(&aul, &m, COMPLEX(rx), INTEGER(apivot),
195 work, &info FCONE);
196 ERROR_LAPACK_2(zsptri, info, 2, D);
197 }
198 }
199 } else {
200 double *work = (double *) R_alloc((size_t) m, sizeof(double ));
201 if (!packed) {
202 F77_CALL(dsytri)(&aul, &m, REAL(rx), &m, INTEGER(apivot),
203 work, &info FCONE);
204 ERROR_LAPACK_2(dsytri, info, 2, D);
205 } else {
206 F77_CALL(dsptri)(&aul, &m, REAL(rx), INTEGER(apivot),
207 work, &info FCONE);
208 ERROR_LAPACK_2(dsptri, info, 2, D);
209 }
210 }
211 } else {
212 PROTECT(rx = (bx == obx) ? duplicateVector(bx) : bx);
213 if (TYPEOF(ax) == CPLXSXP) {
214 if (act == 'C') {
215 if (!packed) {
216 F77_CALL(zhetrs)(&aul, &m, &n, COMPLEX(ax), &m, INTEGER(apivot),
217 COMPLEX(rx), &m, &info FCONE);
218 ERROR_LAPACK_1(zhetrs, info);
219 } else {
220 F77_CALL(zhptrs)(&aul, &m, &n, COMPLEX(ax), INTEGER(apivot),
221 COMPLEX(rx), &m, &info FCONE);
222 ERROR_LAPACK_1(zhptrs, info);
223 }
224 } else {
225 if (!packed) {
226 F77_CALL(zsytrs)(&aul, &m, &n, COMPLEX(ax), &m, INTEGER(apivot),
227 COMPLEX(rx), &m, &info FCONE);
228 ERROR_LAPACK_1(zsytrs, info);
229 } else {
230 F77_CALL(zsptrs)(&aul, &m, &n, COMPLEX(ax), INTEGER(apivot),
231 COMPLEX(rx), &m, &info FCONE);
232 ERROR_LAPACK_1(zsptrs, info);
233 }
234 }
235 } else {
236 if (!packed) {
237 F77_CALL(dsytrs)(&aul, &m, &n, REAL(ax), &m, INTEGER(apivot),
238 REAL(rx), &m, &info FCONE);
239 ERROR_LAPACK_1(dsytrs, info);
240 } else {
241 F77_CALL(dsptrs)(&aul, &m, &n, REAL(ax), INTEGER(apivot),
242 REAL(rx), &m, &info FCONE);
243 ERROR_LAPACK_1(dsptrs, info);
244 }
245 }
246 }
247 SET_SLOT(r, Matrix_xSym, rx);
248 UNPROTECT(2); /* rx, apivot */
249 }
250
251 SOLVE_FINISH;
252
253 UNPROTECT(3); /* r, bx, ax */
254 return r;
255}
256
257SEXP denseCholesky_solve(SEXP s_a, SEXP s_b)
258{
259 SOLVE_START;
260
261 SEXP ax = GET_SLOT(s_a, Matrix_xSym), bx = R_NilValue;
262 if (s_b != R_NilValue) {
263 PROTECT(ax);
264 PROTECT(bx = GET_SLOT(s_b, Matrix_xSym));
265 if (TYPEOF(ax) != TYPEOF(bx)) {
266 ax = Rf_coerceVector(ax, CPLXSXP);
267 bx = Rf_coerceVector(bx, CPLXSXP);
268 }
269 UNPROTECT(2);
270 }
271 PROTECT(ax);
272 PROTECT(bx);
273
274 int packed = XLENGTH(ax) != (int_fast64_t) m * m;
275
276 char rcl[] = "...Matrix";
277 rcl[0] = (TYPEOF(ax) == CPLXSXP) ? 'z' : 'd';
278 if (s_b == R_NilValue) {
279 rcl[1] = 'p';
280 rcl[2] = (!packed) ? 'o' : 'p';
281 } else {
282 rcl[1] = 'g';
283 rcl[2] = 'e';
284 }
285 SEXP r = PROTECT(newObject(rcl));
286
287 SET_DIM(r, m, n);
288
289 char aul = UPLO(s_a);;
290 if (s_b == R_NilValue && aul != 'U')
291 SET_UPLO(r);
292
293 if (m > 0) {
294 SEXP aperm = PROTECT(Rf_getAttrib(s_a, Matrix_permSym)), rx;
295 int info, pivoted = TYPEOF(aperm) == INTSXP && LENGTH(aperm) > 0;
296 if (s_b == R_NilValue) {
297 PROTECT(rx = Rf_allocVector(TYPEOF(ax), XLENGTH(ax)));
298 if (TYPEOF(ax) == CPLXSXP) {
299 memcpy(COMPLEX(rx), COMPLEX(ax), sizeof(Rcomplex) * (size_t) XLENGTH(ax));
300 if (!packed) {
301 F77_CALL(zpotri)(&aul, &m, COMPLEX(rx), &m, &info FCONE);
302 ERROR_LAPACK_2(zpotri, info, 2, L);
303 zsymperm2(COMPLEX(rx), NULL,
304 m, aul, (pivoted) ? INTEGER(aperm) : NULL, 1, 1);
305 } else {
306 F77_CALL(zpptri)(&aul, &m, COMPLEX(rx), &info FCONE);
307 ERROR_LAPACK_2(zpptri, info, 2, L);
308 zsymperm1(COMPLEX(rx), NULL,
309 m, aul, (pivoted) ? INTEGER(aperm) : NULL, 1, 1);
310 }
311 } else {
312 memcpy( REAL(rx), REAL(ax), sizeof( double) * (size_t) XLENGTH(ax));
313 if (!packed) {
314 F77_CALL(dpotri)(&aul, &m, REAL(rx), &m, &info FCONE);
315 ERROR_LAPACK_2(dpotri, info, 2, L);
316 dsymperm2( REAL(rx), NULL,
317 m, aul, (pivoted) ? INTEGER(aperm) : NULL, 1, 1);
318 } else {
319 F77_CALL(dpptri)(&aul, &m, REAL(rx), &info FCONE);
320 ERROR_LAPACK_2(dpptri, info, 2, L);
321 dsymperm1( REAL(rx), NULL,
322 m, aul, (pivoted) ? INTEGER(aperm) : NULL, 1, 1);
323 }
324 }
325 } else {
326 PROTECT(rx = Rf_allocVector(TYPEOF(ax), XLENGTH(bx)));
327 if (TYPEOF(ax) == CPLXSXP) {
328 zrowperm2(COMPLEX(rx), COMPLEX(bx),
329 m, n, (pivoted) ? INTEGER(aperm) : NULL, 1, 0);
330 if (!packed) {
331 F77_CALL(zpotrs)(&aul, &m, &n, COMPLEX(ax), &m,
332 COMPLEX(rx), &m, &info FCONE);
333 ERROR_LAPACK_1(zpotrs, info);
334 } else {
335 F77_CALL(zpptrs)(&aul, &m, &n, COMPLEX(ax),
336 COMPLEX(rx), &m, &info FCONE);
337 ERROR_LAPACK_1(zpptrs, info);
338 }
339 zrowperm2(COMPLEX(rx), NULL,
340 m, n, (pivoted) ? INTEGER(aperm) : NULL, 1, 1);
341 } else {
342 drowperm2( REAL(rx), REAL(bx),
343 m, n, (pivoted) ? INTEGER(aperm) : NULL, 1, 0);
344 if (!packed) {
345 F77_CALL(dpotrs)(&aul, &m, &n, REAL(ax), &m,
346 REAL(rx), &m, &info FCONE);
347 ERROR_LAPACK_1(dpotrs, info);
348 } else {
349 F77_CALL(dpptrs)(&aul, &m, &n, REAL(ax),
350 REAL(rx), &m, &info FCONE);
351 ERROR_LAPACK_1(dpptrs, info);
352 }
353 drowperm2( REAL(rx), NULL,
354 m, n, (pivoted) ? INTEGER(aperm) : NULL, 1, 1);
355 }
356 }
357 SET_SLOT(r, Matrix_xSym, rx);
358 UNPROTECT(2); /* rx, aperm */
359 }
360
361 SOLVE_FINISH;
362
363 UNPROTECT(3); /* r, bx, ax */
364 return r;
365}
366
367SEXP trMatrix_solve(SEXP s_a, SEXP s_b)
368{
369 SOLVE_START;
370
371 SEXP ax = GET_SLOT(s_a, Matrix_xSym), bx = R_NilValue, obx = bx;
372 if (s_b != R_NilValue) {
373 PROTECT(ax);
374 PROTECT(bx = obx = GET_SLOT(s_b, Matrix_xSym));
375 if (TYPEOF(ax) != TYPEOF(bx)) {
376 ax = Rf_coerceVector(ax, CPLXSXP);
377 bx = Rf_coerceVector(bx, CPLXSXP);
378 }
379 UNPROTECT(2);
380 }
381 PROTECT(ax);
382 PROTECT(bx);
383
384 int packed = XLENGTH(ax) != (int_fast64_t) m * m;
385
386 char rcl[] = "...Matrix";
387 rcl[0] = (TYPEOF(ax) == CPLXSXP) ? 'z' : 'd';
388 if (s_b == R_NilValue) {
389 rcl[1] = 't';
390 rcl[2] = (!packed) ? 'r' : 'p';
391 } else {
392 rcl[1] = 'g';
393 rcl[2] = 'e';
394 }
395 SEXP r = PROTECT(newObject(rcl));
396
397 SET_DIM(r, m, n);
398
399 char aul = UPLO(s_a);
400 if (s_b == R_NilValue && aul != 'U')
401 SET_UPLO(r);
402
403 char anu = DIAG(s_a);
404 if (s_b == R_NilValue && anu != 'N')
405 SET_DIAG(r);
406
407 if (m > 0) {
408 SEXP rx;
409 int info;
410 if (s_b == R_NilValue) {
411 PROTECT(rx = duplicateVector(ax));
412 if (TYPEOF(ax) == CPLXSXP) {
413 if (!packed) {
414 F77_CALL(ztrtri)(&aul, &anu, &m, COMPLEX(rx), &m,
415 &info FCONE FCONE);
416 ERROR_LAPACK_2(ztrtri, info, 2, A);
417 } else {
418 F77_CALL(ztptri)(&aul, &anu, &m, COMPLEX(rx),
419 &info FCONE FCONE);
420 ERROR_LAPACK_2(ztptri, info, 2, A);
421 }
422 } else {
423 if (!packed) {
424 F77_CALL(dtrtri)(&aul, &anu, &m, REAL(rx), &m,
425 &info FCONE FCONE);
426 ERROR_LAPACK_2(dtrtri, info, 2, A);
427 } else {
428 F77_CALL(dtptri)(&aul, &anu, &m, REAL(rx),
429 &info FCONE FCONE);
430 ERROR_LAPACK_2(dtptri, info, 2, A);
431 }
432 }
433 } else {
434 PROTECT(rx = (obx == bx) ? duplicateVector(bx) : bx);
435 if (TYPEOF(ax) == CPLXSXP) {
436 if (!packed) {
437 F77_CALL(ztrtrs)(&aul, "N", &anu, &m, &n, COMPLEX(ax), &m,
438 COMPLEX(rx), &m, &info FCONE FCONE FCONE);
439 ERROR_LAPACK_1(ztrtrs, info);
440 } else {
441 F77_CALL(ztptrs)(&aul, "N", &anu, &m, &n, COMPLEX(ax),
442 COMPLEX(rx), &m, &info FCONE FCONE FCONE);
443 ERROR_LAPACK_1(ztptrs, info);
444 }
445 } else {
446 if (!packed) {
447 F77_CALL(dtrtrs)(&aul, "N", &anu, &m, &n, REAL(ax), &m,
448 REAL(rx), &m, &info FCONE FCONE FCONE);
449 ERROR_LAPACK_1(dtrtrs, info);
450 } else {
451 F77_CALL(dtptrs)(&aul, "N", &anu, &m, &n, REAL(ax),
452 REAL(rx), &m, &info FCONE FCONE FCONE);
453 ERROR_LAPACK_1(dtptrs, info);
454 }
455 }
456 }
457 SET_SLOT(r, Matrix_xSym, rx);
458 UNPROTECT(1); /* rx */
459 }
460
461 SOLVE_FINISH;
462
463 UNPROTECT(3); /* r, bx, ax */
464 return r;
465}
466
467SEXP sparseLU_solve(SEXP s_a, SEXP s_b, SEXP s_sparse)
468{
469 SOLVE_START;
470
471 SEXP r,
472 aL = PROTECT(GET_SLOT(s_a, Matrix_LSym)),
473 aU = PROTECT(GET_SLOT(s_a, Matrix_USym)),
474 ap = PROTECT(GET_SLOT(s_a, Matrix_pSym)),
475 aq = PROTECT(GET_SLOT(s_a, Matrix_qSym));
476 int i, j,
477 *pap = (LENGTH(ap)) ? INTEGER(ap) : NULL,
478 *paq = (LENGTH(aq)) ? INTEGER(aq) : NULL;
479 Matrix_cs *L = M2CXS(aL, 1), *U = M2CXS(aU, 1);
480
481#define PROMOTE(A) \
482 do { \
483 size_t n = (size_t) (A)->nzmax; \
484 double *x = (double *) (A)->x; \
485 Rcomplex *y = (Rcomplex *) R_alloc(n, sizeof(Rcomplex)); \
486 while (n-- > 0) { (*(y)).r = *(x++); (*(y++)).i = 0.0; } \
487 (A)->x = y; \
488 (A)->xtype = CXSPARSE_COMPLEX; \
489 } while (0)
490
491 if (L->xtype != CXSPARSE_COMPLEX && s_b != R_NilValue &&
492 TYPEOF(GET_SLOT(s_b, Matrix_xSym)) == CPLXSXP) {
493 PROMOTE(L);
494 PROMOTE(U);
495 }
496 CXSPARSE_XTYPE_SET(L->xtype);
497 if (!Rf_asLogical(s_sparse)) {
498 if ((int_fast64_t) m * n > R_XLEN_T_MAX)
499 Rf_error(_("attempt to allocate vector of length exceeding %s"),
500 "R_XLEN_T_MAX");
501 char rcl[] = ".geMatrix";
502 rcl[0] = (L->xtype == CXSPARSE_COMPLEX) ? 'z' : 'd';
503 PROTECT(r = newObject(rcl));
504 SET_DIM(r, m, n);
505 R_xlen_t mn = (R_xlen_t) m * n;
506 SEXP rx = PROTECT(Rf_allocVector((L->xtype == CXSPARSE_COMPLEX) ? CPLXSXP : REALSXP, mn));
507 if (s_b == R_NilValue) {
508
509#define SOLVE_DENSE(c) \
510 do { \
511 c##TYPE *prx = c##PTR(rx), \
512 *work = (c##TYPE *) R_alloc((size_t) m, sizeof(c##TYPE)); \
513 memset(prx, 0, sizeof(c##TYPE) * (size_t) mn); \
514 for (j = 0; j < n; ++j) { \
515 prx[j] = c##UNIT; \
516 Matrix_cs_pvec(pap, prx, work, m); \
517 Matrix_cs_lsolve(L, work); \
518 Matrix_cs_usolve(U, work); \
519 Matrix_cs_ipvec(paq, work, prx, m); \
520 prx += m; \
521 } \
522 } while (0)
523
524 if (L->xtype == CXSPARSE_COMPLEX)
525 SOLVE_DENSE(z);
526 else
527 SOLVE_DENSE(d);
528
529#undef SOLVE_DENSE
530
531 } else {
532 SEXP bx = GET_SLOT(s_b, Matrix_xSym);
533 if (L->xtype == CXSPARSE_COMPLEX && TYPEOF(bx) != CPLXSXP) {
534 PROTECT(bx);
535 bx = Rf_coerceVector(bx, CPLXSXP);
536 UNPROTECT(1);
537 }
538 PROTECT(bx);
539
540#define SOLVE_DENSE(c) \
541 do { \
542 c##TYPE *prx = c##PTR(rx), *pbx = c##PTR(bx), \
543 *work = (c##TYPE *) R_alloc((size_t) m, sizeof(c##TYPE)); \
544 for (j = 0; j < n; ++j) { \
545 Matrix_cs_pvec(pap, pbx, work, m); \
546 Matrix_cs_lsolve(L, work); \
547 Matrix_cs_usolve(U, work); \
548 Matrix_cs_ipvec(paq, work, prx, m); \
549 prx += m; \
550 pbx += m; \
551 } \
552 } while (0)
553
554 if (L->xtype == CXSPARSE_COMPLEX)
555 SOLVE_DENSE(z);
556 else
557 SOLVE_DENSE(d);
558
559#undef SOLVE_DENSE
560
561 UNPROTECT(1); /* bx */
562 }
563 SET_SLOT(r, Matrix_xSym, rx);
564 UNPROTECT(1); /* rx */
565 } else {
566 Matrix_cs *B = NULL, *X = NULL;
567 if (s_b == R_NilValue) {
568 B = Matrix_cs_speye(m, m, 1, 0);
569 if (B && pap)
570 for (i = 0; i < m; ++i)
571 B->i[pap[i]] = i;
572 } else {
573 B = M2CXS(s_b, 1);
574 if (L->xtype == CXSPARSE_COMPLEX && B->xtype != L->xtype)
575 PROMOTE(B);
576 if (B && pap) {
577 int *papinv = Matrix_cs_pinv(pap, m);
578 if (!papinv)
579 ERROR_OOM(__func__);
580 B = Matrix_cs_permute(B, papinv, NULL, 1);
581 papinv = Matrix_cs_free(papinv);
582 }
583 }
584 if (!B)
585 ERROR_OOM(__func__);
586 int Bfr = s_b == R_NilValue || pap;
587
588 int k, top, nz, nzmax,
589 *iwork = (int *) R_alloc((size_t) m * 2, sizeof(int));
590
591#define SOLVE_SPARSE_TRIANGULAR(c, _A_, _ALO_, _BFR_) \
592 do { \
593 X = Matrix_cs_spalloc(m, n, B->nzmax, 1, 0); \
594 if (!X) { \
595 if (_BFR_) \
596 B = Matrix_cs_spfree(B); \
597 ERROR_OOM(__func__); \
598 } \
599 c##TYPE *X__x = (c##TYPE *) X->x; \
600 X->p[0] = nz = 0; \
601 nzmax = X->nzmax; \
602 for (j = 0, k = 0; j < n; ++j) { \
603 top = Matrix_cs_spsolve(_A_, B, j, iwork, work, NULL, _ALO_); \
604 if (m - top > INT_MAX - nz) { \
605 if (_BFR_) \
606 B = Matrix_cs_spfree(B); \
607 X = Matrix_cs_spfree(X); \
608 Rf_error(_("attempt to construct %s with more than %s nonzero elements"), \
609 "sparseMatrix", "2^31-1"); \
610 } \
611 nz += m - top; \
612 if (nz > nzmax) { \
613 nzmax = (nz <= INT_MAX / 2) ? 2 * nz : INT_MAX; \
614 if (!Matrix_cs_sprealloc(X, nzmax)) { \
615 if (_BFR_) \
616 B = Matrix_cs_spfree(B); \
617 X = Matrix_cs_spfree(X); \
618 ERROR_OOM(__func__); \
619 } \
620 X__x = (c##TYPE *) X->x; \
621 } \
622 X->p[j + 1] = nz; \
623 if (_ALO_) { \
624 for (i = top; i < m; ++i) { \
625 X->i[k] = iwork[i] ; \
626 X__x[k] = work[iwork[i]]; \
627 ++k; \
628 } \
629 } else { \
630 for (i = m - 1; i >= top; --i) { \
631 X->i[k] = iwork[i] ; \
632 X__x[k] = work[iwork[i]]; \
633 ++k; \
634 } \
635 } \
636 } \
637 if (_BFR_) \
638 B = Matrix_cs_spfree(B); \
639 B = X; \
640 } while (0)
641
642#define SOLVE_SPARSE(c) \
643 do { \
644 c##TYPE *work = (c##TYPE *) R_alloc((size_t) m, sizeof(c##TYPE)); \
645 SOLVE_SPARSE_TRIANGULAR(c, L, 1, Bfr); \
646 SOLVE_SPARSE_TRIANGULAR(c, U, 0, 1); \
647 } while (0)
648
649 if (L->xtype == CXSPARSE_COMPLEX)
650 SOLVE_SPARSE(z);
651 else
652 SOLVE_SPARSE(d);
653
654#undef SOLVE_SPARSE
655
656 if (paq) {
657 X = Matrix_cs_permute(B, paq, NULL, 1);
658 B = Matrix_cs_spfree(B);
659 if (!X)
660 ERROR_OOM(__func__);
661 B = X;
662 }
663
664 /* Drop zeros from B and sort it : */
665 Matrix_cs_dropzeros(B);
666 X = Matrix_cs_transpose(B, 1);
667 B = Matrix_cs_spfree(B);
668 if (!X)
669 ERROR_OOM(__func__);
670 B = Matrix_cs_transpose(X, 1);
671 X = Matrix_cs_spfree(X);
672 if (!B)
673 ERROR_OOM(__func__);
674
675 PROTECT(r = CXS2M(B, 1, 'g'));
676 B = Matrix_cs_spfree(B);
677 }
678
679 SOLVE_FINISH;
680
681 UNPROTECT(5); /* r, aq, ap, aU, aL */
682 return r;
683}
684
685static
686int strmatch(const char *s, const char **nms)
687{
688 int i = 0;
689 while (nms[i][0] != '\0') {
690 if (strcmp(s, nms[i]) == 0)
691 return i;
692 ++i;
693 }
694 return -1;
695}
696
697SEXP sparseCholesky_solve(SEXP s_a, SEXP s_b, SEXP s_sparse,
698 SEXP s_system)
699{
700 static const char *valid[] = {
701 "A", "LDLt", "LD", "DLt", "L", "Lt", "D", "P", "Pt", "" };
702 int ivalid = -1;
703 if (TYPEOF(s_system) != STRSXP || LENGTH(s_system) < 1 ||
704 (s_system = STRING_ELT(s_system, 0)) == NA_STRING ||
705 (ivalid = strmatch(CHAR(s_system), valid)) < 0)
706 Rf_error(_("invalid '%s' to '%s'"), "system", __func__);
707
708 SOLVE_START;
709
710 SEXP r;
711 size_t m_ = (size_t) m, n_ = (size_t) n;
712 cholmod_factor *L = M2CHF(s_a, 1);
713 if (!Rf_asLogical(s_sparse)) {
714 if ((int_fast64_t) m * n > R_XLEN_T_MAX)
715 Rf_error(_("attempt to allocate vector of length exceeding %s"),
716 "R_XLEN_T_MAX");
717 cholmod_dense *B = NULL, *X = NULL;
718 if (s_b == R_NilValue) {
719 B = cholmod_eye(m_, n_, L->xtype + L->dtype, &c);
720 if (!B)
721 ERROR_OOM(__func__);
722 X = cholmod_solve(ivalid, L, B, &c);
723 cholmod_free_dense(&B, &c);
724 if (!X)
725 ERROR_OOM(__func__);
726 PROTECT(r = CHD2M(X, 'N', (ivalid < 2) ? ((L->is_ll) ? 'p' : 's') : ((ivalid < 7) ? 't' : 'g')));
727 } else {
728 B = M2CHD(s_b, 'N');
729 X = cholmod_solve(ivalid, L, B, &c);
730 if (!X)
731 ERROR_OOM(__func__);
732 PROTECT(r = CHD2M(X, 'N', 'g'));
733 }
734 cholmod_free_dense(&X, &c);
735 } else {
736 cholmod_sparse *B = NULL, *X = NULL;
737 if (s_b == R_NilValue) {
738 B = cholmod_speye(m_, n_, L->xtype + L->dtype, &c);
739 if (!B)
740 ERROR_OOM(__func__);
741 X = cholmod_spsolve(ivalid, L, B, &c);
742 cholmod_free_sparse(&B, &c);
743 if (X && ivalid < 7) {
744 if (!X->sorted)
745 cholmod_sort(X, &c);
746 B = cholmod_copy(X, (ivalid == 2 || ivalid == 4) ? -1 : 1, 1, &c);
747 cholmod_free_sparse(&X, &c);
748 X = B;
749 }
750 if (!X)
751 ERROR_OOM(__func__);
752 PROTECT(r = CHS2M(X, 1, (ivalid < 2) ? ((L->is_ll) ? 'p' : 's') : ((ivalid < 7) ? 't' : 'g')));
753 } else {
754 B = M2CHS(s_b, 1);
755 X = cholmod_spsolve(ivalid, L, B, &c);
756 if (!X)
757 ERROR_OOM(__func__);
758 PROTECT(r = CHS2M(X, 1, 'g'));
759 }
760 cholmod_free_sparse(&X, &c);
761 }
762 if (s_b == R_NilValue && (ivalid == 2 || ivalid == 4))
763 SET_UPLO(r);
764
765 SOLVE_FINISH;
766
767 UNPROTECT(1); /* r */
768 return r;
769}
770
771SEXP tCMatrix_solve(SEXP s_a, SEXP s_b, SEXP s_sparse)
772{
773 SOLVE_START;
774
775 SEXP r;
776 char aul = UPLO(s_a);
777 int i, j;
778 Matrix_cs *A = M2CXS(s_a, 1);
779 if (A->xtype != CXSPARSE_COMPLEX && s_b != R_NilValue &&
780 TYPEOF(GET_SLOT(s_b, Matrix_xSym)) == CPLXSXP)
781 PROMOTE(A);
782 CXSPARSE_XTYPE_SET(A->xtype);
783 if (!Rf_asLogical(s_sparse)) {
784 if ((int_fast64_t) m * n > R_XLEN_T_MAX)
785 Rf_error(_("attempt to allocate vector of length exceeding %s"),
786 "R_XLEN_T_MAX");
787 char rcl[] = "...Matrix";
788 rcl[0] = (A->xtype == CXSPARSE_COMPLEX) ? 'z' : 'd';
789 rcl[1] = (s_b == R_NilValue) ? 't' : 'g';
790 rcl[2] = (s_b == R_NilValue) ? 'r' : 'e';
791 PROTECT(r = newObject(rcl));
792 SET_DIM(r, m, n);
793 R_xlen_t mn = (R_xlen_t) m * n;
794 SEXP rx = PROTECT(Rf_allocVector((A->xtype == CXSPARSE_COMPLEX) ? CPLXSXP : REALSXP, mn));
795 if (s_b == R_NilValue) {
796
797#define SOLVE_DENSE(c) \
798 do { \
799 c##TYPE *prx = c##PTR(rx); \
800 memset(prx, 0, sizeof(c##TYPE) * (size_t) mn); \
801 for (j = 0; j < n; ++j) { \
802 prx[j] = c##UNIT; \
803 if (aul == 'U') \
804 Matrix_cs_usolve(A, prx); \
805 else \
806 Matrix_cs_lsolve(A, prx); \
807 prx += m; \
808 } \
809 } while (0)
810
811 if (A->xtype == CXSPARSE_COMPLEX)
812 SOLVE_DENSE(z);
813 else
814 SOLVE_DENSE(d);
815
816#undef SOLVE_DENSE
817
818 } else {
819 SEXP bx = GET_SLOT(s_b, Matrix_xSym);
820 if (A->xtype == CXSPARSE_COMPLEX && TYPEOF(bx) != CPLXSXP) {
821 PROTECT(bx);
822 bx = Rf_coerceVector(bx, CPLXSXP);
823 UNPROTECT(1);
824 }
825 PROTECT(bx);
826
827#define SOLVE_DENSE(c) \
828 do { \
829 c##TYPE *prx = c##PTR(rx), *pbx = c##PTR(bx); \
830 memcpy(prx, pbx, sizeof(c##TYPE) * (size_t) mn); \
831 for (j = 0; j < n; ++j) { \
832 if (aul == 'U') \
833 Matrix_cs_usolve(A, prx); \
834 else \
835 Matrix_cs_lsolve(A, prx); \
836 prx += m; \
837 } \
838 } while (0)
839
840 if (A->xtype == CXSPARSE_COMPLEX)
841 SOLVE_DENSE(z);
842 else
843 SOLVE_DENSE(d);
844
845#undef SOLVE_DENSE
846
847 UNPROTECT(1); /* bx */
848 }
849 SET_SLOT(r, Matrix_xSym, rx);
850 UNPROTECT(1); /* rx */
851 } else {
852 Matrix_cs *B = NULL, *X = NULL;
853 if (s_b == R_NilValue)
854 B = Matrix_cs_speye(m, m, 1, 0);
855 else {
856 B = M2CXS(s_b, 1);
857 if (A->xtype == CXSPARSE_COMPLEX && B->xtype != A->xtype)
858 PROMOTE(B);
859 }
860 if (!B)
861 ERROR_OOM(__func__);
862
863 int k, top, nz, nzmax,
864 *iwork = (int *) R_alloc((size_t) m * 2, sizeof(int));
865
866#define SOLVE_SPARSE(c) \
867 do { \
868 c##TYPE *work = (c##TYPE *) R_alloc((size_t) m, sizeof(c##TYPE)); \
869 SOLVE_SPARSE_TRIANGULAR(c, A, aul != 'U', s_b == R_NilValue); \
870 } while (0)
871
872 if (A->xtype == CXSPARSE_COMPLEX)
873 SOLVE_SPARSE(z);
874 else
875 SOLVE_SPARSE(d);
876
877#undef SOLVE_SPARSE
878
879 /* Drop zeros from B and sort it : */
880 Matrix_cs_dropzeros(B);
881 X = Matrix_cs_transpose(B, 1);
882 B = Matrix_cs_spfree(B);
883 if (!X)
884 ERROR_OOM(__func__);
885 B = Matrix_cs_transpose(X, 1);
886 X = Matrix_cs_spfree(X);
887 if (!B)
888 ERROR_OOM(__func__);
889
890 PROTECT(r = CXS2M(B, 1, (s_b == R_NilValue) ? 't' : 'g'));
891 B = Matrix_cs_spfree(B);
892 }
893 if (s_b == R_NilValue && aul != 'U')
894 SET_UPLO(r);
895
896 SOLVE_FINISH;
897
898 UNPROTECT(1); /* r */
899 return r;
900}
901
902SEXP sparseQR_matmult(SEXP s_qr, SEXP s_y, SEXP s_op,
903 SEXP s_complete, SEXP s_yxjj)
904{
905 SEXP V = PROTECT(GET_SLOT(s_qr, Matrix_VSym));
906 Matrix_cs *V_ = M2CXS(V, 1);
907 if (V_->xtype != CXSPARSE_COMPLEX &&
908 ((s_y != R_NilValue &&
909 TYPEOF(GET_SLOT(s_y, Matrix_xSym)) == CPLXSXP) ||
910 (s_yxjj != R_NilValue &&
911 TYPEOF(s_yxjj) == CPLXSXP)))
912 PROMOTE(V_);
913 CXSPARSE_XTYPE_SET(V_->xtype);
914
915 SEXP beta = PROTECT(GET_SLOT(s_qr, Matrix_betaSym));
916 double *pbeta = REAL(beta);
917
918 SEXP p = PROTECT(GET_SLOT(s_qr, Matrix_pSym));
919 int *pp = (LENGTH(p) > 0) ? INTEGER(p) : NULL;
920
921 int m = V_->m, r = V_->n, n, i, j, op = Rf_asInteger(s_op),
922 nprotect = 5;
923
924 SEXP yx;
925 if (s_y == R_NilValue) {
926 n = (Rf_asLogical(s_complete)) ? m : r;
927 if ((int_fast64_t) m * n > R_XLEN_T_MAX)
928 Rf_error(_("attempt to allocate vector of length exceeding %s"),
929 "R_XLEN_T_MAX");
930 R_xlen_t mn = (R_xlen_t) m * n;
931 PROTECT(yx = Rf_allocVector((V_->xtype == CXSPARSE_COMPLEX) ? CPLXSXP : REALSXP, mn));
932
933#define TEMPLATE(c) \
934 do { \
935 c##TYPE *pyx = c##PTR(yx); \
936 memset(pyx, 0, sizeof(c##TYPE) * (size_t) mn); \
937 if (s_yxjj == R_NilValue) { \
938 c##TYPE one = c##UNIT; \
939 c##NAME(copy2)((size_t) n, \
940 pyx, (size_t) m + 1, &one, 0); \
941 } else { \
942 s_yxjj = Rf_coerceVector(s_yxjj, c##TYPESXP); \
943 c##TYPE *pyxjj = c##PTR(s_yxjj); \
944 c##NAME(copy2)((size_t) n, \
945 pyx, (size_t) m + 1, pyxjj, 1); \
946 } \
947 } while (0)
948
949 if (V_->xtype == CXSPARSE_COMPLEX)
950 TEMPLATE(z);
951 else
952 TEMPLATE(d);
953
954#undef TEMPLATE
955
956 } else {
957 int *pydim = DIM(s_y);
958 if (pydim[0] != m)
959 Rf_error(_("dimensions of '%s' and '%s' are inconsistent"),
960 "qr", "y");
961 n = pydim[1];
962
963 yx = GET_SLOT(s_y, Matrix_xSym);
964 if (V_->xtype == CXSPARSE_COMPLEX && TYPEOF(yx) != CPLXSXP) {
965 PROTECT(yx);
966 yx = Rf_coerceVector(yx, CPLXSXP);
967 UNPROTECT(1);
968 }
969 PROTECT(yx);
970 }
971
972 char acl[] = ".geMatrix";
973 acl[0] = (V_->xtype == CXSPARSE_COMPLEX) ? 'z' : 'd';
974 SEXP a = PROTECT(newObject(acl));
975
976 int *padim = DIM(a);
977 padim[0] = (op != 0) ? m : r;
978 padim[1] = n;
979
980 SEXP ax;
981 if (s_y == R_NilValue && padim[0] == m)
982 ax = yx;
983 else {
984 R_xlen_t mn = (R_xlen_t) padim[0] * padim[1];
985 PROTECT(ax = Rf_allocVector((V_->xtype == CXSPARSE_COMPLEX) ? CPLXSXP : REALSXP, mn));
986 ++nprotect;
987 }
988 SET_SLOT(a, Matrix_xSym, ax);
989
990#define TEMPLATE(c) \
991 do { \
992 c##TYPE *pyx = c##PTR(yx), *pax = c##PTR(ax), *work = NULL; \
993 if (op < 5) \
994 work = (c##TYPE *) R_alloc((size_t) m, sizeof(c##TYPE)); \
995 switch (op) { \
996 case 0: /* qr.coef : A = P2 R1^{-1} Q1' P1 y */ \
997 { \
998 SEXP R = PROTECT(GET_SLOT(s_qr, Matrix_RSym)), \
999 q = PROTECT(GET_SLOT(s_qr, Matrix_qSym)); \
1000 Matrix_cs *R_ = M2CXS(R, 1); \
1001 if (V_->xtype == CXSPARSE_COMPLEX && R_->xtype != V_->xtype) \
1002 PROMOTE(R_); \
1003 int *pq = (LENGTH(q) > 0) ? INTEGER(q) : NULL; \
1004 for (j = 0; j < n; ++j) { \
1005 Matrix_cs_pvec(pp, pyx, work, m); \
1006 for (i = 0; i < r; ++i) \
1007 Matrix_cs_happly(V_, i, pbeta[i], work); \
1008 Matrix_cs_usolve(R_, work); \
1009 Matrix_cs_ipvec(pq, work, pax, r); \
1010 pyx += m; \
1011 pax += r; \
1012 } \
1013 UNPROTECT(2); /* q, R */ \
1014 break; \
1015 } \
1016 case 1: /* qr.fitted : A = P1' Q1 Q1' P1 y */ \
1017 for (j = 0; j < n; ++j) { \
1018 Matrix_cs_pvec(pp, pyx, work, m); \
1019 for (i = 0; i < r; ++i) \
1020 Matrix_cs_happly(V_, i, pbeta[i], work); \
1021 if (r < m) \
1022 memset(work + r, 0, sizeof(c##TYPE) * (size_t) (m - r)); \
1023 for (i = r - 1; i >= 0; --i) \
1024 Matrix_cs_happly(V_, i, pbeta[i], work); \
1025 Matrix_cs_ipvec(pp, work, pax, m); \
1026 pyx += m; \
1027 pax += m; \
1028 } \
1029 break; \
1030 case 2: /* qr.resid : A = P1' Q2 Q2' P1 y */ \
1031 for (j = 0; j < n; ++j) { \
1032 Matrix_cs_pvec(pp, pyx, work, m); \
1033 for (i = 0; i < r; ++i) \
1034 Matrix_cs_happly(V_, i, pbeta[i], work); \
1035 if (r > 0) \
1036 memset(work, 0, sizeof(c##TYPE) * (size_t) r); \
1037 for (i = r - 1; i >= 0; --i) \
1038 Matrix_cs_happly(V_, i, pbeta[i], work); \
1039 Matrix_cs_ipvec(pp, work, pax, m); \
1040 pyx += m; \
1041 pax += m; \
1042 } \
1043 break; \
1044 case 3: /* qr.qty {w/ perm.} : A = Q' P1 y */ \
1045 for (j = 0; j < n; ++j) { \
1046 Matrix_cs_pvec(pp, pyx, work, m); \
1047 memcpy(pax, work, sizeof(c##TYPE) * (size_t) m); \
1048 for (i = 0; i < r; ++i) \
1049 Matrix_cs_happly(V_, i, pbeta[i], pax); \
1050 pyx += m; \
1051 pax += m; \
1052 } \
1053 break; \
1054 case 4: /* qr.qy {w/ perm.} : A = P1' Q y */ \
1055 for (j = 0; j < n; ++j) { \
1056 memcpy(work, pyx, sizeof(c##TYPE) * (size_t) m); \
1057 for (i = r - 1; i >= 0; --i) \
1058 Matrix_cs_happly(V_, i, pbeta[i], work); \
1059 Matrix_cs_ipvec(pp, work, pax, m); \
1060 pyx += m; \
1061 pax += m; \
1062 } \
1063 break; \
1064 case 5: /* qr.qty {w/o perm.} : A = Q' y */ \
1065 if (ax != yx) \
1066 memcpy(pax, pyx, sizeof(c##TYPE) * (size_t) m * (size_t) n); \
1067 for (j = 0; j < n; ++j) { \
1068 for (i = 0; i < r; ++i) \
1069 Matrix_cs_happly(V_, i, pbeta[i], pax); \
1070 pax += m; \
1071 } \
1072 break; \
1073 case 6: /* qr.qy {w/o perm.} : A = Q y */ \
1074 if (ax != yx) \
1075 memcpy(pax, pyx, sizeof(c##TYPE) * (size_t) m * (size_t) n); \
1076 for (j = 0; j < n; ++j) { \
1077 for (i = r - 1; i >= 0; --i) \
1078 Matrix_cs_happly(V_, i, pbeta[i], pax); \
1079 pax += m; \
1080 } \
1081 break; \
1082 default: \
1083 Rf_error(_("invalid '%s' to '%s'"), "op", __func__); \
1084 break; \
1085 } \
1086 } while (0)
1087
1088 if (V_->xtype == CXSPARSE_COMPLEX)
1089 TEMPLATE(z);
1090 else
1091 TEMPLATE(d);
1092
1093#undef TEMPLATE
1094
1095 UNPROTECT(nprotect); /* ax, a, yx, p, beta, V */
1096 return a;
1097}
FCONE
#define FCONE
Definition Lapack-etc.h:13
ERROR_LAPACK_1
#define ERROR_LAPACK_1(_ROUTINE_, _INFO_)
Definition Lapack-etc.h:18
ERROR_LAPACK_2
#define ERROR_LAPACK_2(_ROUTINE_, _INFO_, _WARN_, _LETTER_)
Definition Lapack-etc.h:25
SET_DIM
#define SET_DIM(x, m, n)
Definition Mdefines.h:87
ERROR_OOM
#define ERROR_OOM(_FUNC_)
Definition Mdefines.h:141
_
#define _(String)
Definition Mdefines.h:66
SET_UPLO
#define SET_UPLO(x)
Definition Mdefines.h:103
DIAG
#define DIAG(x)
Definition Mdefines.h:111
UPLO
#define UPLO(x)
Definition Mdefines.h:101
duplicateVector
SEXP duplicateVector(SEXP)
Definition utils.c:42
SET_TRANS
#define SET_TRANS(x)
Definition Mdefines.h:108
TRANS
#define TRANS(x)
Definition Mdefines.h:106
SET_DIAG
#define SET_DIAG(x)
Definition Mdefines.h:113
DIM
#define DIM(x)
Definition Mdefines.h:85
GET_SLOT
#define GET_SLOT(x, name)
Definition Mdefines.h:72
newObject
SEXP newObject(const char *)
Definition objects.c:13
SET_SLOT
#define SET_SLOT(x, name, value)
Definition Mdefines.h:73
TYPEOF
#define TYPEOF(s)
Definition Mdefines.h:123
M2CHF
cholmod_factor * M2CHF(SEXP obj, int values)
Definition cholmod-etc.c:39
M2CHS
cholmod_sparse * M2CHS(SEXP obj, int values)
Definition cholmod-etc.c:132
c
cholmod_common c
Definition cholmod-etc.c:5
CHS2M
SEXP CHS2M(cholmod_sparse *A, int values, char shape)
Definition cholmod-etc.c:328
M2CHD
cholmod_dense * M2CHD(SEXP obj, char trans)
Definition cholmod-etc.c:176
CHD2M
SEXP CHD2M(cholmod_dense *A, char trans, char shape)
Definition cholmod-etc.c:375
Matrix_cs_transpose
Matrix_cs * Matrix_cs_transpose(const Matrix_cs *A, int values)
Definition cs-etc.c:385
Matrix_cs_spfree
Matrix_cs * Matrix_cs_spfree(Matrix_cs *A)
Definition cs-etc.c:337
Matrix_cs_free
void * Matrix_cs_free(void *p)
Definition cs-etc.c:113
Matrix_cs_permute
Matrix_cs * Matrix_cs_permute(const Matrix_cs *A, const int *pinv, const int *q, int values)
Definition cs-etc.c:197
CXS2M
SEXP CXS2M(Matrix_cs *A, int values, char shape)
Definition cs-etc.c:42
Matrix_cs_speye
Matrix_cs * Matrix_cs_speye(int m, int n, int values, int triplet)
Definition cs-etc.c:307
Matrix_cs_pinv
int * Matrix_cs_pinv(const int *p, int n)
Definition cs-etc.c:222
Matrix_cs_dropzeros
int Matrix_cs_dropzeros(Matrix_cs *A)
Definition cs-etc.c:101
M2CXS
Matrix_cs * M2CXS(SEXP obj, int values)
Definition cs-etc.c:6
CXSPARSE_XTYPE_SET
#define CXSPARSE_XTYPE_SET(_VALUE_)
Definition cs-etc.h:12
CXSPARSE_COMPLEX
#define CXSPARSE_COMPLEX
Definition cs-etc.h:9
dsymperm2
void dsymperm2(double *, const double *, int, char, const int *, int, int)
zsymperm2
void zsymperm2(Rcomplex *, const Rcomplex *, int, char, const int *, int, int)
zsymperm1
void zsymperm1(Rcomplex *, const Rcomplex *, int, char, const int *, int, int)
drowperm2
void drowperm2(double *, const double *, int, int, const int *, int, int)
dsymperm1
void dsymperm1(double *, const double *, int, char, const int *, int, int)
zrowperm2
void zrowperm2(Rcomplex *, const Rcomplex *, int, int, const int *, int, int)
Matrix_permSym
SEXP Matrix_permSym
Definition init.c:623
Matrix_betaSym
SEXP Matrix_betaSym
Definition init.c:603
Matrix_xSym
SEXP Matrix_xSym
Definition init.c:635
Matrix_VSym
SEXP Matrix_VSym
Definition init.c:602
Matrix_LSym
SEXP Matrix_LSym
Definition init.c:599
Matrix_qSym
SEXP Matrix_qSym
Definition init.c:627
Matrix_USym
SEXP Matrix_USym
Definition init.c:601
Matrix_pSym
SEXP Matrix_pSym
Definition init.c:622
sparseCholesky_solve
SEXP sparseCholesky_solve(SEXP s_a, SEXP s_b, SEXP s_sparse, SEXP s_system)
Definition solve.c:697
denseBunchKaufman_solve
SEXP denseBunchKaufman_solve(SEXP s_a, SEXP s_b)
Definition solve.c:130
denseLU_solve
SEXP denseLU_solve(SEXP s_a, SEXP s_b)
Definition solve.c:33
SOLVE_SPARSE
#define SOLVE_SPARSE(c)
SOLVE_FINISH
#define SOLVE_FINISH
PROMOTE
#define PROMOTE(A)
TEMPLATE
#define TEMPLATE(c)
strmatch
static int strmatch(const char *s, const char **nms)
Definition solve.c:686
sparseQR_matmult
SEXP sparseQR_matmult(SEXP s_qr, SEXP s_y, SEXP s_op, SEXP s_complete, SEXP s_yxjj)
Definition solve.c:902
tCMatrix_solve
SEXP tCMatrix_solve(SEXP s_a, SEXP s_b, SEXP s_sparse)
Definition solve.c:771
sparseLU_solve
SEXP sparseLU_solve(SEXP s_a, SEXP s_b, SEXP s_sparse)
Definition solve.c:467
solveDN
static void solveDN(SEXP rdn, SEXP adn, SEXP bdn)
Definition solve.c:11
SOLVE_START
#define SOLVE_START
SOLVE_DENSE
#define SOLVE_DENSE(c)
denseCholesky_solve
SEXP denseCholesky_solve(SEXP s_a, SEXP s_b)
Definition solve.c:257
trMatrix_solve
SEXP trMatrix_solve(SEXP s_a, SEXP s_b)
Definition solve.c:367
Matrix_cs_sparse
Definition cs-etc.h:17
Matrix_cs_sparse::i
int * i
Definition cs-etc.h:22
Matrix_cs_sparse::n
int n
Definition cs-etc.h:20
Matrix_cs_sparse::m
int m
Definition cs-etc.h:19
Matrix_cs_sparse::xtype
int xtype
Definition cs-etc.h:25
Generated by doxygen 1.12.0