Matrix r4655
Loading...
Searching...
No Matches
factor.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 "factor.h"
6
7/* defined in ./attrib.c : */
8SEXP get_factor(SEXP, const char *);
9void set_factor(SEXP, const char *, SEXP);
10
11static
12SEXP dgeMatrix_trf_(SEXP obj, int warn)
13{
14 SEXP val = PROTECT(newObject("denseLU")),
15 dim = PROTECT(GET_SLOT(obj, Matrix_DimSym)),
16 dimnames = PROTECT(GET_SLOT(obj, Matrix_DimNamesSym));
17 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1], r = (m < n) ? m : n;
18 SET_SLOT(val, Matrix_DimSym, dim);
19 SET_SLOT(val, Matrix_DimNamesSym, dimnames);
20 if (r > 0) {
21 SEXP perm = PROTECT(allocVector(INTSXP, r)),
22 x = PROTECT(GET_SLOT(obj, Matrix_xSym)),
23 y = PROTECT(allocVector(TYPEOF(x), XLENGTH(x)));
24 int *pperm = INTEGER(perm), info;
25#ifdef MATRIX_ENABLE_ZMATRIX
26 if (TYPEOF(x) == CPLXSXP) {
27 Rcomplex *px = COMPLEX(x), *py = COMPLEX(y);
28 Matrix_memcpy(py, px, XLENGTH(y), sizeof(Rcomplex));
29 F77_CALL(zgetrf)(&m, &n, py, &m, pperm, &info);
30 ERROR_LAPACK_2(zgetrf, info, warn, U);
31 } else {
32#endif
33 double *px = REAL(x), *py = REAL(y);
34 Matrix_memcpy(py, px, XLENGTH(y), sizeof(double));
35 F77_CALL(dgetrf)(&m, &n, py, &m, pperm, &info);
36 ERROR_LAPACK_2(dgetrf, info, warn, U);
37#ifdef MATRIX_ENABLE_ZMATRIX
38 }
39#endif
40 SET_SLOT(val, Matrix_permSym, perm);
41 SET_SLOT(val, Matrix_xSym, y);
42 UNPROTECT(3); /* y, x, perm */
43 }
44 UNPROTECT(3); /* dimnames, dim, val */
45 return val;
46}
47
48static
49SEXP dsyMatrix_trf_(SEXP obj, int warn)
50{
51 SEXP val = PROTECT(newObject("BunchKaufman")),
52 dim = PROTECT(GET_SLOT(obj, Matrix_DimSym)),
53 dimnames = PROTECT(GET_SLOT(obj, Matrix_DimNamesSym)),
54 uplo = PROTECT(GET_SLOT(obj, Matrix_uploSym));
55 int n = INTEGER(dim)[1];
56 char ul = *CHAR(STRING_ELT(uplo, 0));
57 SET_SLOT(val, Matrix_DimSym, dim);
58 set_symmetrized_DimNames(val, dimnames, -1);
59 SET_SLOT(val, Matrix_uploSym, uplo);
60 if (n > 0) {
61 SEXP perm = PROTECT(allocVector(INTSXP, n)),
62 x = PROTECT(GET_SLOT(obj, Matrix_xSym)),
63 y = PROTECT(allocVector(TYPEOF(x), XLENGTH(x)));
64 int *pperm = INTEGER(perm), info, lwork = -1;
65#ifdef MATRIX_ENABLE_ZMATRIX
66 if (TYPEOF(x) == CPLXSXP) {
67 Rcomplex *px = COMPLEX(x), *py = COMPLEX(y), tmp, *work;
68 Matrix_memset(py, 0, XLENGTH(y), sizeof(Rcomplex));
69 F77_CALL(zlacpy)(&ul, &n, &n, px, &n, py, &n FCONE);
70 F77_CALL(zsytrf)(&ul, &n, py, &n, pperm, &tmp, &lwork, &info FCONE);
71 lwork = (int) tmp.r;
72 Matrix_Calloc(work, lwork, Rcomplex);
73 F77_CALL(zsytrf)(&ul, &n, py, &n, pperm, work, &lwork, &info FCONE);
74 Matrix_Free(work, lwork);
75 ERROR_LAPACK_2(zsytrf, info, warn, D);
76 } else {
77#endif
78 double *px = REAL(x), *py = REAL(y), tmp, *work;
79 Matrix_memset(py, 0, XLENGTH(y), sizeof(double));
80 F77_CALL(dlacpy)(&ul, &n, &n, px, &n, py, &n FCONE);
81 F77_CALL(dsytrf)(&ul, &n, py, &n, pperm, &tmp, &lwork, &info FCONE);
82 lwork = (int) tmp;
83 Matrix_Calloc(work, lwork, double);
84 F77_CALL(dsytrf)(&ul, &n, py, &n, pperm, work, &lwork, &info FCONE);
85 Matrix_Free(work, lwork);
86 ERROR_LAPACK_2(dsytrf, info, warn, D);
87#ifdef MATRIX_ENABLE_ZMATRIX
88 }
89#endif
90 SET_SLOT(val, Matrix_permSym, perm);
91 SET_SLOT(val, Matrix_xSym, y);
92 UNPROTECT(3); /* y, x, perm */
93 }
94 UNPROTECT(4); /* uplo, dimnames, dim, val */
95 return val;
96}
97
98static
99SEXP dspMatrix_trf_(SEXP obj, int warn)
100{
101 SEXP val = PROTECT(newObject("pBunchKaufman")),
102 dim = PROTECT(GET_SLOT(obj, Matrix_DimSym)),
103 dimnames = PROTECT(GET_SLOT(obj, Matrix_DimNamesSym)),
104 uplo = PROTECT(GET_SLOT(obj, Matrix_uploSym));
105 int n = INTEGER(dim)[1];
106 char ul = *CHAR(STRING_ELT(uplo, 0));
107 SET_SLOT(val, Matrix_DimSym, dim);
108 set_symmetrized_DimNames(val, dimnames, -1);
109 SET_SLOT(val, Matrix_uploSym, uplo);
110 if (n > 0) {
111 SEXP perm = PROTECT(allocVector(INTSXP, n)),
112 x = PROTECT(GET_SLOT(obj, Matrix_xSym)),
113 y = PROTECT(allocVector(TYPEOF(x), XLENGTH(x)));
114 int *pperm = INTEGER(perm), info;
115#ifdef MATRIX_ENABLE_ZMATRIX
116 if (TYPEOF(x) == CPLXSXP) {
117 Rcomplex *px = COMPLEX(x), *py = COMPLEX(y);
118 Matrix_memcpy(py, px, XLENGTH(y), sizeof(Rcomplex));
119 F77_CALL(zsptrf)(&ul, &n, py, pperm, &info FCONE);
120 ERROR_LAPACK_2(zsptrf, info, warn, D);
121 } else {
122#endif
123 double *px = REAL(x), *py = REAL(y);
124 Matrix_memcpy(py, px, XLENGTH(y), sizeof(double));
125 F77_CALL(dsptrf)(&ul, &n, py, pperm, &info FCONE);
126 ERROR_LAPACK_2(dsptrf, info, warn, D);
127#ifdef MATRIX_ENABLE_ZMATRIX
128 }
129#endif
130 SET_SLOT(val, Matrix_permSym, perm);
131 SET_SLOT(val, Matrix_xSym, y);
132 UNPROTECT(3); /* y, x, perm */
133 }
134 UNPROTECT(4); /* uplo, dimnames, dim, val */
135 return val;
136}
137
138static
139SEXP dpoMatrix_trf_(SEXP obj, int warn, int pivot, double tol)
140{
141 SEXP val = PROTECT(newObject("Cholesky")),
142 dim = PROTECT(GET_SLOT(obj, Matrix_DimSym)),
143 dimnames = PROTECT(GET_SLOT(obj, Matrix_DimNamesSym)),
144 uplo = PROTECT(GET_SLOT(obj, Matrix_uploSym));
145 int n = INTEGER(dim)[1];
146 char ul = *CHAR(STRING_ELT(uplo, 0));
147 SET_SLOT(val, Matrix_DimSym, dim);
148 set_symmetrized_DimNames(val, dimnames, -1);
149 SET_SLOT(val, Matrix_uploSym, uplo);
150 if (n > 0) {
151 SEXP x = PROTECT(GET_SLOT(obj, Matrix_xSym)),
152 y = PROTECT(allocVector(TYPEOF(x), XLENGTH(x)));
153 int info;
154#ifdef MATRIX_ENABLE_ZMATRIX
155 if (TYPEOF(x) == CPLXSXP) {
156 Rcomplex *px = COMPLEX(x), *py = COMPLEX(y);
157 Matrix_memset(py, 0, XLENGTH(y), sizeof(Rcomplex));
158 F77_CALL(zlacpy)(&ul, &n, &n, px, &n, py, &n FCONE);
159 if (!pivot) {
160 F77_CALL(zpotrf)(&ul, &n, py, &n, &info FCONE);
161 ERROR_LAPACK_3(zpotrf, info, warn, 6);
162 } else {
163 SEXP perm = PROTECT(allocVector(INTSXP, n));
164 int *pperm = INTEGER(perm), rank;
165 Rcomplex *work = (Rcomplex *) R_alloc((size_t) 2 * n, sizeof(Rcomplex));
166 F77_CALL(zpstrf)(&ul, &n, py, &n, pperm, &rank, &tol, work, &info FCONE);
167 ERROR_LAPACK_4(zpstrf, info, rank, warn);
168 if (info > 0) {
169 int j, d = n - rank;
170 py += (R_xlen_t) rank * n + rank;
171 for (j = rank; j < n; ++j) {
172 Matrix_memset(py, 0, d, sizeof(Rcomplex));
173 py += n;
174 }
175 }
176 SET_SLOT(val, Matrix_permSym, perm);
177 UNPROTECT(1); /* perm */
178 }
179 } else {
180#endif
181 double *px = REAL(x), *py = REAL(y);
182 Matrix_memset(py, 0, XLENGTH(y), sizeof(double));
183 F77_CALL(dlacpy)(&ul, &n, &n, px, &n, py, &n FCONE);
184 if (!pivot) {
185 F77_CALL(dpotrf)(&ul, &n, py, &n, &info FCONE);
186 ERROR_LAPACK_3(dpotrf, info, warn, 6);
187 } else {
188 SEXP perm = PROTECT(allocVector(INTSXP, n));
189 int *pperm = INTEGER(perm), rank;
190 double *work = (double *) R_alloc((size_t) 2 * n, sizeof(double));
191 F77_CALL(dpstrf)(&ul, &n, py, &n, pperm, &rank, &tol, work, &info FCONE);
192 ERROR_LAPACK_4(dpstrf, info, rank, warn);
193 if (info > 0) {
194 int j, d = n - rank;
195 py += (R_xlen_t) rank * n + rank;
196 for (j = rank; j < n; ++j) {
197 Matrix_memset(py, 0, d, sizeof(double));
198 py += n;
199 }
200 }
201 SET_SLOT(val, Matrix_permSym, perm);
202 UNPROTECT(1); /* perm */
203 }
204#ifdef MATRIX_ENABLE_ZMATRIX
205 }
206#endif
207 SET_SLOT(val, Matrix_xSym, y);
208 UNPROTECT(2); /* y, x */
209 }
210 UNPROTECT(4); /* uplo, dimnames, dim, val */
211 return val;
212}
213
214static
215SEXP dppMatrix_trf_(SEXP obj, int warn)
216{
217 SEXP val = PROTECT(newObject("pCholesky")),
218 dim = PROTECT(GET_SLOT(obj, Matrix_DimSym)),
219 dimnames = PROTECT(GET_SLOT(obj, Matrix_DimNamesSym)),
220 uplo = PROTECT(GET_SLOT(obj, Matrix_uploSym));
221 int n = INTEGER(dim)[1];
222 char ul = *CHAR(STRING_ELT(uplo, 0));
223 SET_SLOT(val, Matrix_DimSym, dim);
224 set_symmetrized_DimNames(val, dimnames, -1);
225 SET_SLOT(val, Matrix_uploSym, uplo);
226 if (n > 0) {
227 SEXP x = PROTECT(GET_SLOT(obj, Matrix_xSym)),
228 y = PROTECT(allocVector(TYPEOF(x), XLENGTH(x)));
229 int info;
230#ifdef MATRIX_ENABLE_ZMATRIX
231 if (TYPEOF(x) == CPLXSXP) {
232 Rcomplex *px = COMPLEX(x), *py = COMPLEX(y);
233 Matrix_memcpy(py, px, XLENGTH(y), sizeof(Rcomplex));
234 F77_CALL(zpptrf)(&ul, &n, py, &info FCONE);
235 ERROR_LAPACK_3(zpptrf, info, warn, 6);
236 } else {
237#endif
238 double *px = REAL(x), *py = REAL(y);
239 Matrix_memcpy(py, px, XLENGTH(y), sizeof(double));
240 F77_CALL(dpptrf)(&ul, &n, py, &info FCONE);
241 ERROR_LAPACK_3(dpptrf, info, warn, 6);
242#ifdef MATRIX_ENABLE_ZMATRIX
243 }
244#endif
245 SET_SLOT(val, Matrix_xSym, y);
246 UNPROTECT(2); /* y, x */
247 }
248 UNPROTECT(4); /* uplo, dimnames, dim, val */
249 return val;
250}
251
252SEXP dgeMatrix_trf(SEXP obj, SEXP warn)
253{
254 SEXP val = get_factor(obj, "denseLU");
255 if (isNull(val)) {
256 PROTECT(val = dgeMatrix_trf_(obj, asInteger(warn)));
257 set_factor(obj, "denseLU", val);
258 UNPROTECT(1);
259 }
260 return val;
261}
262
263SEXP dsyMatrix_trf(SEXP obj, SEXP warn)
264{
265 SEXP val = get_factor(obj, "BunchKaufman");
266 if (isNull(val)) {
267 PROTECT(val = dsyMatrix_trf_(obj, asInteger(warn)));
268 set_factor(obj, "BunchKaufman", val);
269 UNPROTECT(1);
270 }
271 return val;
272}
273
274SEXP dspMatrix_trf(SEXP obj, SEXP warn)
275{
276 SEXP val = get_factor(obj, "pBunchKaufman");
277 if (isNull(val)) {
278 PROTECT(val = dspMatrix_trf_(obj, asInteger(warn)));
279 set_factor(obj, "pBunchKaufman", val);
280 UNPROTECT(1);
281 }
282 return val;
283}
284
285SEXP dpoMatrix_trf(SEXP obj, SEXP warn, SEXP pivot, SEXP tol)
286{
287 int pivot_ = asLogical(pivot);
288 SEXP val = get_factor(obj, (pivot_) ? "Cholesky~" : "Cholesky");
289 if (isNull(val)) {
290 double tol_ = asReal(tol);
291 PROTECT(val = dpoMatrix_trf_(obj, asInteger(warn), pivot_, tol_));
292 set_factor(obj, (pivot_) ? "Cholesky~" : "Cholesky", val);
293 UNPROTECT(1);
294 }
295 return val;
296}
297
298SEXP dppMatrix_trf(SEXP obj, SEXP warn)
299{
300 SEXP val = get_factor(obj, "pCholesky");
301 if (isNull(val)) {
302 PROTECT(val = dppMatrix_trf_(obj, asInteger(warn)));
303 set_factor(obj, "pCholesky", val);
304 UNPROTECT(1);
305 }
306 return val;
307}
308
309SEXP dgeMatrix_sch(SEXP x, SEXP vectors, SEXP isDGE)
310{
311// 'x' is either a traditional matrix or a dgeMatrix, as indicated by isDGE.
312 int *dims, n, vecs = asLogical(vectors), is_dge = asLogical(isDGE),
313 info, izero = 0, lwork = -1, nprot = 1;
314
315 if(is_dge) {
316 dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
317 } else { // traditional matrix
318 dims = INTEGER(getAttrib(x, R_DimSymbol));
319 if(!isReal(x)) { // may not be "numeric" ..
320 x = PROTECT(coerceVector(x, REALSXP)); // -> maybe error
321 nprot++;
322 }
323 }
324 double *work, tmp;
325 const char *nms[] = {"WR", "WI", "T", "Z", ""};
326 SEXP val = PROTECT(Rf_mkNamed(VECSXP, nms));
327
328 n = dims[0];
329 if (n != dims[1] || n < 1)
330 error(_("dgeMatrix_Schur: argument x must be a non-null square matrix"));
331 const R_xlen_t n2 = ((R_xlen_t)n) * n; // = n^2
332
333 SET_VECTOR_ELT(val, 0, allocVector(REALSXP, n));
334 SET_VECTOR_ELT(val, 1, allocVector(REALSXP, n));
335 SET_VECTOR_ELT(val, 2, allocMatrix(REALSXP, n, n));
336 Memcpy(REAL(VECTOR_ELT(val, 2)),
337 REAL(is_dge ? GET_SLOT(x, Matrix_xSym) : x),
338 n2);
339 SET_VECTOR_ELT(val, 3, allocMatrix(REALSXP, vecs ? n : 0, vecs ? n : 0));
340 F77_CALL(dgees)(vecs ? "V" : "N", "N", NULL, dims, (double *) NULL, dims, &izero,
341 (double *) NULL, (double *) NULL, (double *) NULL, dims,
342 &tmp, &lwork, (int *) NULL, &info FCONE FCONE);
343 if (info) error(_("dgeMatrix_Schur: first call to dgees failed"));
344 lwork = (int) tmp;
345 Matrix_Calloc(work, lwork, double);
346
347 F77_CALL(dgees)(vecs ? "V" : "N", "N", NULL, dims, REAL(VECTOR_ELT(val, 2)), dims,
348 &izero, REAL(VECTOR_ELT(val, 0)), REAL(VECTOR_ELT(val, 1)),
349 REAL(VECTOR_ELT(val, 3)), dims, work, &lwork,
350 (int *) NULL, &info FCONE FCONE);
351 Matrix_Free(work, lwork);
352 if (info) error(_("dgeMatrix_Schur: dgees returned code %d"), info);
353 UNPROTECT(nprot);
354 return val;
355}
356
357#define DO_FREE(_A_, _S_, _N_) \
358do { \
359 if (!(_A_)) \
360 _A_ = Matrix_cs_spfree(_A_); \
361 if (!(_S_)) \
362 _S_ = Matrix_cs_sfree (_S_); \
363 if (!(_N_)) \
364 _N_ = Matrix_cs_nfree (_N_); \
365} while (0)
366
367#define DO_SORT(_A_) \
368do { \
369 Matrix_cs_dropzeros(_A_); \
370 T = Matrix_cs_transpose(_A_, 1); \
371 if (!T) { \
372 DO_FREE(T, *S, *N); \
373 return 0; \
374 } \
375 _A_ = Matrix_cs_spfree(_A_); \
376 _A_ = Matrix_cs_transpose(T, 1); \
377 if (!(_A_)) { \
378 DO_FREE(T, *S, *N); \
379 return 0; \
380 } \
381 T = Matrix_cs_spfree(T); \
382} while (0)
383
384static
385int dgCMatrix_trf_(const Matrix_cs *A, Matrix_css **S, Matrix_csn **N,
386 int order, double tol)
387{
388 Matrix_cs *T = NULL;
389 if (!(*S = Matrix_cs_sqr(order, A, 0)) ||
390 !(*N = Matrix_cs_lu(A, *S, tol))) {
391 DO_FREE(T, *S, *N);
392 return 0;
393 }
394 DO_SORT((*N)->L);
395 DO_SORT((*N)->U);
396 return 1;
397}
398
399SEXP dgCMatrix_trf(SEXP obj, SEXP order, SEXP tol, SEXP doError)
400{
401 double tol_ = asReal(tol);
402 if (ISNAN(tol_))
403 error(_("'%s' is not a number"), "tol");
404
405 int order_ = asInteger(order);
406 if (order_ == NA_INTEGER)
407 order_ = (tol_ == 1.0) ? 2 : 1;
408 else if (order_ < 0 || order_ > 3)
409 order_ = 0;
410
411 SEXP val = get_factor(obj, (order_) ? "sparseLU~" : "sparseLU");
412 if (!isNull(val))
413 return val;
414 PROTECT(val = newObject("sparseLU"));
415
416 Matrix_cs *A = M2CXS(obj, 1);
417 MCS_XTYPE_SET(A->xtype);
418
419 Matrix_css *S = NULL;
420 Matrix_csn *N = NULL;
421 int *P = NULL;
422
423 if (A->m != A->n)
424 error(_("LU factorization of m-by-n %s requires m == n"),
425 ".gCMatrix");
426 if (!dgCMatrix_trf_(A, &S, &N, order_, tol_) ||
427 !(P = Matrix_cs_pinv(N->pinv, A->m))) {
428 if (!P) {
429 S = Matrix_cs_sfree(S);
430 N = Matrix_cs_nfree(N);
431 }
432 if (asLogical(doError))
433 error(_("LU factorization of %s failed: out of memory or near-singular"),
434 ".gCMatrix");
435 /* Defensive code will check with is(., "sparseLU") : */
436 UNPROTECT(1); /* val */
437 return ScalarLogical(NA_LOGICAL);
438 }
439
440 SEXP dim = PROTECT(GET_SLOT(obj, Matrix_DimSym));
441 SET_SLOT(val, Matrix_DimSym, dim);
442 UNPROTECT(1); /* dim */
443
444 SEXP dimnames = PROTECT(GET_SLOT(obj, Matrix_DimNamesSym));
445 SET_SLOT(val, Matrix_DimNamesSym, dimnames);
446 UNPROTECT(1); /* dimnames */
447
448 SEXP L = PROTECT(CXS2M(N->L, 1, 't')),
449 U = PROTECT(CXS2M(N->U, 1, 't')),
450 uplo = PROTECT(mkString("L"));
451 SET_SLOT(L, Matrix_uploSym, uplo);
452 SET_SLOT(val, Matrix_LSym, L);
453 SET_SLOT(val, Matrix_USym, U);
454 UNPROTECT(3); /* uplo, U, L */
455
456 SEXP p = PROTECT(allocVector(INTSXP, A->m));
457 Matrix_memcpy(INTEGER(p), P, A->m, sizeof(int));
458 SET_SLOT(val, Matrix_pSym, p);
459 UNPROTECT(1); /* p */
460 if (order_ > 0) {
461 SEXP q = PROTECT(allocVector(INTSXP, A->n));
462 Matrix_memcpy(INTEGER(q), S->q, A->n, sizeof(int));
463 SET_SLOT(val, Matrix_qSym, q);
464 UNPROTECT(1); /* q */
465 }
466
467 S = Matrix_cs_sfree(S);
468 N = Matrix_cs_nfree(N);
469 P = Matrix_cs_free(P);
470
471 set_factor(obj, (order_) ? "sparseLU~" : "sparseLU", val);
472 UNPROTECT(1); /* val */
473 return val;
474}
475
476static
477int dgCMatrix_orf_(const Matrix_cs *A, Matrix_css **S, Matrix_csn **N,
478 int order)
479{
480 Matrix_cs *T = NULL;
481 if (!(*S = Matrix_cs_sqr(order, A, 1)) ||
482 !(*N = Matrix_cs_qr(A, *S))) {
483 DO_FREE(T, *S, *N);
484 return 0;
485 }
486 DO_SORT((*N)->L);
487 DO_SORT((*N)->U);
488 return 1;
489}
490
491SEXP dgCMatrix_orf(SEXP obj, SEXP order, SEXP doError)
492{
493 int order_ = asInteger(order);
494 if (order_ < 0 || order_ > 3)
495 order_ = 0;
496
497 SEXP val = get_factor(obj, (order_) ? "sparseQR~" : "sparseQR");
498 if (!isNull(val))
499 return val;
500 PROTECT(val = newObject("sparseQR"));
501
502 Matrix_cs *A = M2CXS(obj, 1);
503 MCS_XTYPE_SET(A->xtype);
504
505 Matrix_css *S = NULL;
506 Matrix_csn *N = NULL;
507 int *P = NULL;
508
509 if (A->m < A->n)
510 error(_("QR factorization of m-by-n %s requires m >= n"),
511 ".gCMatrix");
512 if (!dgCMatrix_orf_(A, &S, &N, order_) ||
513 !(P = Matrix_cs_pinv(S->pinv, S->m2))) {
514 if (!P) {
515 S = Matrix_cs_sfree(S);
516 N = Matrix_cs_nfree(N);
517 }
518 if (asLogical(doError))
519 error(_("QR factorization of %s failed: out of memory"),
520 ".gCMatrix");
521 /* Defensive code will check with is(., "sparseQR") : */
522 UNPROTECT(1); /* val */
523 return ScalarLogical(NA_LOGICAL);
524 }
525
526 SEXP dim = PROTECT(GET_SLOT(obj, Matrix_DimSym));
527 SET_SLOT(val, Matrix_DimSym, dim);
528 UNPROTECT(1); /* dim */
529
530 SEXP dimnames = PROTECT(GET_SLOT(obj, Matrix_DimNamesSym));
531 SET_SLOT(val, Matrix_DimNamesSym, dimnames);
532 UNPROTECT(1); /* dimnames */
533
534 SEXP V = PROTECT(CXS2M(N->L, 1, 'g')),
535 R = PROTECT(CXS2M(N->U, 1, 'g'));
536 SET_SLOT(val, Matrix_VSym, V);
537 SET_SLOT(val, Matrix_RSym, R);
538 UNPROTECT(2); /* R, V */
539
540 SEXP beta = PROTECT(allocVector(REALSXP, A->n));
541 Matrix_memcpy(REAL(beta), N->B, A->n, sizeof(double));
542 SET_SLOT(val, Matrix_betaSym, beta);
543 UNPROTECT(1); /* beta */
544
545 SEXP p = PROTECT(allocVector(INTSXP, S->m2));
546 Matrix_memcpy(INTEGER(p), P, S->m2, sizeof(int));
547 SET_SLOT(val, Matrix_pSym, p);
548 UNPROTECT(1); /* p */
549 if (order_ > 0) {
550 SEXP q = PROTECT(allocVector(INTSXP, A->n));
551 Matrix_memcpy(INTEGER(q), S->q, A->n, sizeof(int));
552 SET_SLOT(val, Matrix_qSym, q);
553 UNPROTECT(1); /* q */
554 }
555
556 S = Matrix_cs_sfree(S);
557 N = Matrix_cs_nfree(N);
558 P = Matrix_cs_free(P);
559
560 set_factor(obj, (order_) ? "sparseQR~" : "sparseQR", val);
561 UNPROTECT(1); /* val */
562 return val;
563}
564
565#undef DO_FREE
566#undef DO_SORT
567
568static
569int dpCMatrix_trf_(cholmod_sparse *A, cholmod_factor **L,
570 int perm, int ldl, int super, double mult)
571{
572 /* defined in ./cholmod-common.c : */
573 void R_cholmod_common_envget(void);
574 void R_cholmod_common_envset(void);
575
576 R_cholmod_common_envset();
577
578 if (*L == NULL) {
579 if (perm == 0) {
580 c.nmethods = 1;
581 c.method[0].ordering = CHOLMOD_NATURAL;
582 c.postorder = 0;
583 }
584
585 c.supernodal = (super == NA_LOGICAL) ? CHOLMOD_AUTO :
586 ((super != 0) ? CHOLMOD_SUPERNODAL : CHOLMOD_SIMPLICIAL);
587
588 *L = cholmod_analyze(A, &c);
589 }
590
591 if (super == NA_LOGICAL)
592 super = (*L)->is_super;
593 if (super != 0)
594 ldl = 0;
595
596 c.final_asis = 0;
597 c.final_super = super != 0;
598 c.final_ll = ldl == 0;
599 c.final_pack = 1;
600 c.final_monotonic = 1;
601
602 double beta[2];
603 beta[0] = mult;
604 beta[1] = 0.0;
605 int res = cholmod_factorize_p(A, beta, NULL, 0, *L, &c);
606
607 R_cholmod_common_envget();
608
609 return res;
610}
611
612SEXP dpCMatrix_trf(SEXP obj,
613 SEXP perm, SEXP ldl, SEXP super, SEXP mult)
614{
615 int perm_ = asLogical(perm), ldl_ = asLogical(ldl),
616 super_ = asLogical(super);
617 double mult_ = asReal(mult);
618 if (!R_FINITE(mult_))
619 error(_("'%s' is not a number or not finite"), "mult");
620
621 SEXP trf = R_NilValue;
622 char nm[] = "spdCholesky";
623 if (perm_)
624 nm[1] = 'P';
625 if (super_ != NA_LOGICAL && super_ != 0)
626 ldl_ = 0;
627 if (super_ == NA_LOGICAL || super_ == 0) {
628 if (ldl_)
629 nm[2] = 'D';
630 trf = get_factor(obj, nm);
631 }
632 if (isNull(trf) && (super_ == NA_LOGICAL || super_ != 0)) {
633 nm[0] = 'S';
634 nm[2] = 'd';
635 trf = get_factor(obj, nm);
636 }
637
638 int cached = !isNull(trf);
639 if (cached && mult_ == 0.0)
640 return trf;
641
642 PROTECT_INDEX pid;
643 PROTECT_WITH_INDEX(trf, &pid);
644 cholmod_sparse *A = M2CHS(obj, 1);
645 cholmod_factor *L = NULL;
646
647 SEXP uplo = GET_SLOT(obj, Matrix_uploSym);
648 char ul = *CHAR(STRING_ELT(uplo, 0));
649 A->stype = (ul == 'U') ? 1 : -1;
650
651 if (cached) {
652 L = M2CHF(trf, 1);
653 L = cholmod_copy_factor(L, &c);
654 dpCMatrix_trf_(A, &L, perm_, ldl_, super_, mult_);
655 } else {
656 dpCMatrix_trf_(A, &L, perm_, ldl_, super_, mult_);
657 if (super_ == NA_LOGICAL) {
658 nm[0] = (L->is_super) ? 'S' : 's';
659 nm[2] = (L->is_ll ) ? 'd' : 'D';
660 }
661 }
662 REPROTECT(trf = CHF2M(L, 1), pid);
663 cholmod_free_factor(&L, &c);
664
665 SEXP dimnames = PROTECT(GET_SLOT(obj, Matrix_DimNamesSym));
666 set_symmetrized_DimNames(trf, dimnames, -1);
667 UNPROTECT(1); /* dimnames */
668
669 if (!cached && mult_ == 0.0)
670 set_factor(obj, nm, trf);
671 UNPROTECT(1); /* trf */
672 return trf;
673}
674
675SEXP BunchKaufman_expand(SEXP obj, SEXP packed)
676{
677 SEXP P_ = PROTECT(newObject("pMatrix")),
678 T_ = PROTECT(newObject("dtCMatrix")),
679 D_ = PROTECT(newObject("dsCMatrix")),
680 dim = PROTECT(GET_SLOT(obj, Matrix_DimSym));
681 int i, j, s, n = INTEGER(dim)[0];
682 R_xlen_t n1a = (R_xlen_t) n + 1;
683 if (n > 0) {
684 SET_SLOT(P_, Matrix_DimSym, dim);
685 SET_SLOT(T_, Matrix_DimSym, dim);
686 SET_SLOT(D_, Matrix_DimSym, dim);
687 }
688 UNPROTECT(1); /* dim */
689
690 SEXP uplo = PROTECT(GET_SLOT(obj, Matrix_uploSym));
691 int upper = *CHAR(STRING_ELT(uplo, 0)) == 'U';
692 if (!upper) {
693 SET_SLOT(T_, Matrix_uploSym, uplo);
694 SET_SLOT(D_, Matrix_uploSym, uplo);
695 }
696 UNPROTECT(1); /* uplo */
697
698 SEXP diag = PROTECT(mkString("U"));
699 SET_SLOT(T_, Matrix_diagSym, diag);
700 UNPROTECT(1); /* diag */
701
702 SEXP pivot = PROTECT(GET_SLOT(obj, Matrix_permSym)),
703 D_p = PROTECT(allocVector(INTSXP, n1a));
704 int *ppivot = INTEGER(pivot), *D_pp = INTEGER(D_p),
705 b = n, dp = (upper) ? 1 : 2;
706 D_pp[0] = 0;
707 j = 0;
708 while (j < n) {
709 if (ppivot[j] > 0) {
710 D_pp[j+1] = D_pp[j] + 1;
711 j += 1;
712 } else {
713 D_pp[j+1] = D_pp[j] + dp;
714 D_pp[j+2] = D_pp[j] + 3;
715 j += 2;
716 --b;
717 }
718 }
719 SET_SLOT(D_, Matrix_pSym, D_p);
720 UNPROTECT(1); /* D_p */
721
722 SEXP P, P_perm, T, T_p, T_i, T_x,
723 D_i = PROTECT(allocVector(INTSXP, D_pp[n])),
724 D_x = PROTECT(allocVector(REALSXP, D_pp[n])),
725 x = PROTECT(GET_SLOT(obj, Matrix_xSym));
726 int *P_pperm, *T_pp, *T_pi, *D_pi = INTEGER(D_i);
727 double *T_px, *D_px = REAL(D_x), *px = REAL(x);
728
729 int unpacked = !asLogical(packed);
730
731 R_xlen_t len = (R_xlen_t) 2 * b + 1, k = (upper) ? len - 2 : 0;
732 SEXP res = PROTECT(allocVector(VECSXP, len));
733
734 j = 0;
735 while (b--) {
736 s = (ppivot[j] > 0) ? 1 : 2;
737 dp = (upper) ? j : n - j - s;
738
739 PROTECT(P = duplicate(P_));
740 PROTECT(P_perm = allocVector(INTSXP, n));
741 PROTECT(T = duplicate(T_));
742 PROTECT(T_p = allocVector(INTSXP, n1a));
743 PROTECT(T_i = allocVector(INTSXP, (R_xlen_t) s * dp));
744 PROTECT(T_x = allocVector(REALSXP, (R_xlen_t) s * dp));
745
746 P_pperm = INTEGER(P_perm);
747 T_pp = INTEGER(T_p);
748 T_pi = INTEGER(T_i);
749 T_px = REAL(T_x);
750 T_pp[0] = 0;
751
752 for (i = 0; i < j; ++i) {
753 T_pp[i+1] = 0;
754 P_pperm[i] = i + 1;
755 }
756 for (i = j; i < j+s; ++i) {
757 T_pp[i+1] = T_pp[i] + dp;
758 P_pperm[i] = i + 1;
759 }
760 for (i = j+s; i < n; ++i) {
761 T_pp[i+1] = T_pp[i];
762 P_pperm[i] = i + 1;
763 }
764
765 if (s == 1) {
766 P_pperm[j] = ppivot[j];
767 P_pperm[ppivot[j]-1] = j + 1;
768 } else if (upper) {
769 P_pperm[j] = -ppivot[j];
770 P_pperm[-ppivot[j]-1] = j + 1;
771 } else {
772 P_pperm[j+1] = -ppivot[j];
773 P_pperm[-ppivot[j]-1] = j + 2;
774 }
775
776 if (upper) {
777 for (i = 0; i < j; ++i) {
778 *(T_pi++) = i;
779 *(T_px++) = *(px++);
780 }
781 *(D_pi++) = j;
782 *(D_px++) = *(px++);
783 ++j;
784 if (unpacked)
785 px += n - j;
786 if (s == 2) {
787 for (i = 0; i < j-1; ++i) {
788 *(T_pi++) = i;
789 *(T_px++) = *(px++);
790 }
791 *(D_pi++) = j - 1;
792 *(D_pi++) = j;
793 *(D_px++) = *(px++);
794 *(D_px++) = *(px++);
795 ++j;
796 if (unpacked)
797 px += n - j;
798 }
799 } else {
800 if (s == 2) {
801 *(D_pi++) = j;
802 *(D_pi++) = j + 1;
803 *(D_px++) = *(px++);
804 *(D_px++) = *(px++);
805 for (i = j+2; i < n; ++i) {
806 *(T_pi++) = i;
807 *(T_px++) = *(px++);
808 }
809 ++j;
810 if (unpacked)
811 px += j;
812 }
813 *(D_pi++) = j;
814 *(D_px++) = *(px++);
815 for (i = j+1; i < n; ++i) {
816 *(T_pi++) = i;
817 *(T_px++) = *(px++);
818 }
819 ++j;
820 if (unpacked)
821 px += j;
822 }
823
824 SET_SLOT(P, Matrix_permSym, P_perm);
825 SET_SLOT(T, Matrix_pSym, T_p);
826 SET_SLOT(T, Matrix_iSym, T_i);
827 SET_SLOT(T, Matrix_xSym, T_x);
828
829 if (upper) {
830 SET_VECTOR_ELT(res, k-1, P);
831 SET_VECTOR_ELT(res, k , T);
832 k -= 2;
833 } else {
834 SET_VECTOR_ELT(res, k , P);
835 SET_VECTOR_ELT(res, k+1, T);
836 k += 2;
837 }
838 UNPROTECT(6); /* T_x, T_i, T_p, T, P_perm, P */
839 }
840
841 SET_SLOT(D_, Matrix_iSym, D_i);
842 SET_SLOT(D_, Matrix_xSym, D_x);
843 SET_VECTOR_ELT(res, len-1, D_);
844
845 UNPROTECT(8); /* res, x, D_x, D_i, pivot, D_, T_, P_ */
846 return res;
847}
848
849SEXP CHMfactor_diag_get(SEXP obj, SEXP square)
850{
851 cholmod_factor *L = M2CHF(obj, 1);
852 int n = (int) L->n, square_ = asLogical(square);
853 SEXP y = PROTECT(allocVector(REALSXP, n));
854 double *py = REAL(y);
855 if (L->is_super) {
856 int k, j, nc,
857 nsuper = (int) L->nsuper,
858 *psuper = (int *) L->super,
859 *ppi = (int *) L->pi,
860 *ppx = (int *) L->px;
861 double *px = (double *) L->x, *px_;
862 R_xlen_t nr1a;
863 for (k = 0; k < nsuper; ++k) {
864 nc = psuper[k+1] - psuper[k];
865 nr1a = (R_xlen_t) (ppi[k+1] - ppi[k]) + 1;
866 px_ = px + ppx[k];
867 for (j = 0; j < nc; ++j) {
868 *py = *px_;
869 if (square_)
870 *py *= *py;
871 ++py;
872 px_ += nr1a;
873 }
874 }
875 } else {
876 square_ = square_ && L->is_ll;
877 int j, *pp = (int *) L->p;
878 double *px = (double *) L->x;
879 for (j = 0; j < n; ++j) {
880 *py = px[pp[j]];
881 if (square_)
882 *py *= *py;
883 ++py;
884 }
885 }
886 UNPROTECT(1);
887 return y;
888}
889
890SEXP CHMfactor_update(SEXP obj, SEXP parent, SEXP mult)
891{
892 /* defined in ./objects.c : */
893 char Matrix_shape(SEXP);
894
895 double mult_ = asReal(mult);
896 if (!R_FINITE(mult_))
897 error(_("'%s' is not a number or not finite"), "mult");
898
899 cholmod_factor *L = cholmod_copy_factor(M2CHF(obj, 1), &c);
900 cholmod_sparse *A = M2CHS(parent, 1);
901 if (Matrix_shape(parent) == 's') {
902 SEXP uplo = GET_SLOT(parent, Matrix_uploSym);
903 char ul = *CHAR(STRING_ELT(uplo, 0));
904 A->stype = (ul == 'U') ? 1 : -1;
905 }
906
907 dpCMatrix_trf_(A, &L, 0, !L->is_ll, L->is_super, mult_);
908
909 SEXP res = PROTECT(CHF2M(L, 1));
910 cholmod_free_factor(&L, &c);
911
912 SEXP dimnames = PROTECT(GET_SLOT(obj, Matrix_DimNamesSym));
913 SET_SLOT(res, Matrix_DimNamesSym, dimnames);
914 UNPROTECT(1);
915
916 UNPROTECT(1);
917 return res;
918}
919
920SEXP CHMfactor_updown(SEXP obj, SEXP parent, SEXP update)
921{
922 /* defined in ./objects.c : */
923 char Matrix_shape(SEXP);
924
925 cholmod_factor *L = cholmod_copy_factor(M2CHF(obj, 1), &c);
926 cholmod_sparse *A = M2CHS(parent, 1);
927 if (Matrix_shape(parent) == 's') {
928 SEXP uplo = GET_SLOT(parent, Matrix_uploSym);
929 char ul = *CHAR(STRING_ELT(uplo, 0));
930 A->stype = (ul == 'U') ? 1 : -1;
931 }
932
933 cholmod_updown(asLogical(update) != 0, A, L, &c);
934
935 SEXP res = PROTECT(CHF2M(L, 1));
936 cholmod_free_factor(&L, &c);
937
938 SEXP dimnames = PROTECT(GET_SLOT(obj, Matrix_DimNamesSym));
939 SET_SLOT(res, Matrix_DimNamesSym, dimnames);
940 UNPROTECT(1);
941
942 UNPROTECT(1);
943 return res;
944}
ERROR_LAPACK_4
#define ERROR_LAPACK_4(_ROUTINE_, _INFO_, _RANK_, _WARN_)
Definition Lapack-etc.h:57
FCONE
#define FCONE
Definition Lapack-etc.h:18
ERROR_LAPACK_3
#define ERROR_LAPACK_3(_ROUTINE_, _INFO_, _WARN_, _NPROTECT_)
Definition Lapack-etc.h:41
ERROR_LAPACK_2
#define ERROR_LAPACK_2(_ROUTINE_, _INFO_, _WARN_, _LETTER_)
Definition Lapack-etc.h:28
_
#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_RSym
SEXP Matrix_RSym
Definition Msymbols.h:6
Matrix_xSym
SEXP Matrix_xSym
Definition Msymbols.h:22
Matrix_iSym
SEXP Matrix_iSym
Definition Msymbols.h:13
Matrix_DimNamesSym
SEXP Matrix_DimNamesSym
Definition Msymbols.h:2
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
set_symmetrized_DimNames
void set_symmetrized_DimNames(SEXP obj, SEXP dn, int J)
Definition attrib.c:132
R_cholmod_common_envset
void R_cholmod_common_envset(void)
Definition cholmod-common.c:1181
R_cholmod_common_envget
void R_cholmod_common_envget(void)
Definition cholmod-common.c:1217
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
CHF2M
SEXP CHF2M(cholmod_factor *L, int values)
Definition cholmod-etc.c:174
c
cholmod_common c
Definition cholmod-etc.c:5
Matrix_cs_lu
Matrix_csn * Matrix_cs_lu(const Matrix_cs *A, const Matrix_css *S, double tol)
Definition cs-etc.c:162
Matrix_cs_free
void * Matrix_cs_free(void *p)
Definition cs-etc.c:114
Matrix_cs_qr
Matrix_csn * Matrix_cs_qr(const Matrix_cs *A, const Matrix_css *S)
Definition cs-etc.c:247
CXS2M
SEXP CXS2M(Matrix_cs *A, int values, char shape)
Definition cs-etc.c:40
Matrix_cs_sfree
Matrix_css * Matrix_cs_sfree(Matrix_css *S)
Definition cs-etc.c:271
Matrix_cs_pinv
int * Matrix_cs_pinv(const int *p, int n)
Definition cs-etc.c:223
Matrix_cs_sqr
Matrix_css * Matrix_cs_sqr(int order, const Matrix_cs *A, int qr)
Definition cs-etc.c:374
M2CXS
Matrix_cs * M2CXS(SEXP obj, int values)
Definition cs-etc.c:6
Matrix_cs_nfree
Matrix_csn * Matrix_cs_nfree(Matrix_csn *N)
Definition cs-etc.c:186
MCS_XTYPE_SET
#define MCS_XTYPE_SET(_VALUE_)
Definition cs-etc.h:12
BunchKaufman_expand
SEXP BunchKaufman_expand(SEXP obj, SEXP packed)
Definition factor.c:675
dsyMatrix_trf_
static SEXP dsyMatrix_trf_(SEXP obj, int warn)
Definition factor.c:49
dgCMatrix_orf
SEXP dgCMatrix_orf(SEXP obj, SEXP order, SEXP doError)
Definition factor.c:491
dgeMatrix_sch
SEXP dgeMatrix_sch(SEXP x, SEXP vectors, SEXP isDGE)
Definition factor.c:309
dgCMatrix_orf_
static int dgCMatrix_orf_(const Matrix_cs *A, Matrix_css **S, Matrix_csn **N, int order)
Definition factor.c:477
DO_FREE
#define DO_FREE(_A_, _S_, _N_)
Definition factor.c:357
CHMfactor_updown
SEXP CHMfactor_updown(SEXP obj, SEXP parent, SEXP update)
Definition factor.c:920
dppMatrix_trf
SEXP dppMatrix_trf(SEXP obj, SEXP warn)
Definition factor.c:298
dpoMatrix_trf
SEXP dpoMatrix_trf(SEXP obj, SEXP warn, SEXP pivot, SEXP tol)
Definition factor.c:285
DO_SORT
#define DO_SORT(_A_)
Definition factor.c:367
set_factor
void set_factor(SEXP, const char *, SEXP)
Definition attrib.c:203
dpCMatrix_trf_
static int dpCMatrix_trf_(cholmod_sparse *A, cholmod_factor **L, int perm, int ldl, int super, double mult)
Definition factor.c:569
dspMatrix_trf_
static SEXP dspMatrix_trf_(SEXP obj, int warn)
Definition factor.c:99
dppMatrix_trf_
static SEXP dppMatrix_trf_(SEXP obj, int warn)
Definition factor.c:215
dpCMatrix_trf
SEXP dpCMatrix_trf(SEXP obj, SEXP perm, SEXP ldl, SEXP super, SEXP mult)
Definition factor.c:612
CHMfactor_update
SEXP CHMfactor_update(SEXP obj, SEXP parent, SEXP mult)
Definition factor.c:890
dpoMatrix_trf_
static SEXP dpoMatrix_trf_(SEXP obj, int warn, int pivot, double tol)
Definition factor.c:139
dgeMatrix_trf_
static SEXP dgeMatrix_trf_(SEXP obj, int warn)
Definition factor.c:12
dsyMatrix_trf
SEXP dsyMatrix_trf(SEXP obj, SEXP warn)
Definition factor.c:263
get_factor
SEXP get_factor(SEXP, const char *)
Definition attrib.c:189
dgeMatrix_trf
SEXP dgeMatrix_trf(SEXP obj, SEXP warn)
Definition factor.c:252
dgCMatrix_trf
SEXP dgCMatrix_trf(SEXP obj, SEXP order, SEXP tol, SEXP doError)
Definition factor.c:399
dspMatrix_trf
SEXP dspMatrix_trf(SEXP obj, SEXP warn)
Definition factor.c:274
dgCMatrix_trf_
static int dgCMatrix_trf_(const Matrix_cs *A, Matrix_css **S, Matrix_csn **N, int order, double tol)
Definition factor.c:385
CHMfactor_diag_get
SEXP CHMfactor_diag_get(SEXP obj, SEXP square)
Definition factor.c:849
Matrix_shape
char Matrix_shape(SEXP obj)
Definition objects.c:102
Matrix_cs_numeric
Definition cs-etc.h:41
Matrix_cs_numeric::B
double * B
Definition cs-etc.h:45
Matrix_cs_numeric::L
Matrix_cs * L
Definition cs-etc.h:42
Matrix_cs_numeric::pinv
int * pinv
Definition cs-etc.h:44
Matrix_cs_numeric::U
Matrix_cs * U
Definition cs-etc.h:43
Matrix_cs_sparse
Definition cs-etc.h:17
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
Matrix_cs_symbolic
Definition cs-etc.h:29
Matrix_cs_symbolic::q
int * q
Definition cs-etc.h:31
Matrix_cs_symbolic::pinv
int * pinv
Definition cs-etc.h:30
Matrix_cs_symbolic::m2
int m2
Definition cs-etc.h:35
Matrix_memset
void * Matrix_memset(void *dest, int ch, R_xlen_t length, size_t size)
Definition utils.c:8
Matrix_memcpy
void * Matrix_memcpy(void *dest, const void *src, R_xlen_t length, size_t size)
Definition utils.c:70
Generated by doxygen 1.9.8