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