9void solveDN(SEXP rdn, SEXP adn, SEXP bdn)
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));
21 SET_STRING_ELT(s, 0, STRING_ELT(adn, 1));
23 SET_STRING_ELT(s, 1, STRING_ELT(bdn, 1));
24 setAttrib(rdn, R_NamesSymbol, s);
35 SEXP adim = GET_SLOT(a, Matrix_DimSym); \
36 int *padim = INTEGER(adim), m = padim[0], n = padim[1]; \
38 error(_("'%s' is not square"), "a"); \
40 SEXP bdim = GET_SLOT(b, Matrix_DimSym); \
41 int *pbdim = INTEGER(bdim); \
43 error(_("dimensions of '%s' and '%s' are inconsistent"), \
49 SEXP rdimnames = PROTECT(GET_SLOT(r, Matrix_DimNamesSym)), \
50 adimnames = PROTECT(GET_SLOT(a, Matrix_DimNamesSym)); \
52 revDN(rdimnames, adimnames); \
54 SEXP bdimnames = PROTECT(GET_SLOT(b, Matrix_DimNamesSym)); \
55 solveDN(rdimnames, adimnames, bdimnames); \
64 char rcl[] =
".geMatrix";
65 rcl[0] = (TYPEOF(ax) == CPLXSXP) ?
'z' :
'd';
69 int *prdim = INTEGER(rdim);
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),
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),
93 double work0, *work = &work0;
94 F77_CALL(dgetri)(&m, REAL(rx), &m, INTEGER(apivot),
98 work = (
double *) R_alloc((
size_t) lwork,
sizeof(double ));
99 F77_CALL(dgetri)(&m, REAL(rx), &m, INTEGER(apivot),
100 work, &lwork, &info);
102#ifdef MATRIX_ENABLE_ZMATRIX
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);
117 F77_CALL(dgetrs)(
"N", &m, &n, REAL(ax), &m, INTEGER(apivot),
118 REAL(rx), &m, &info
FCONE);
120#ifdef MATRIX_ENABLE_ZMATRIX
140 XLENGTH(ax) == (R_xlen_t) m * m;
142 char rcl[] =
"...Matrix";
143 rcl[0] = (TYPEOF(ax) == CPLXSXP) ?
'z' :
'd';
149 rcl[2] = (unpacked) ?
'y' :
'p';
154 int *prdim = INTEGER(rdim);
159 char aul = CHAR(STRING_ELT(auplo, 0))[0];
160 if (isNull(b) && aul !=
'U') {
172#ifdef MATRIX_ENABLE_ZMATRIX
173 if (TYPEOF(ax) == CPLXSXP) {
174 Rcomplex *work = (Rcomplex *) R_alloc((
size_t) m,
sizeof(Rcomplex));
176 F77_CALL(zsytri)(&aul, &m, COMPLEX(rx), &m, INTEGER(apivot),
180 F77_CALL(zsptri)(&aul, &m, COMPLEX(rx), INTEGER(apivot),
186 double *work = (
double *) R_alloc((
size_t) m,
sizeof(double ));
188 F77_CALL(dsytri)(&aul, &m, REAL(rx), &m, INTEGER(apivot),
192 F77_CALL(dsptri)(&aul, &m, REAL(rx), INTEGER(apivot),
196#ifdef MATRIX_ENABLE_ZMATRIX
204#ifdef MATRIX_ENABLE_ZMATRIX
205 if (TYPEOF(ax) == CPLXSXP) {
207 F77_CALL(zsytrs)(&aul, &m, &n, COMPLEX(ax), &m, INTEGER(apivot),
208 COMPLEX(rx), &m, &info
FCONE);
211 F77_CALL(zsptrs)(&aul, &m, &n, COMPLEX(ax), INTEGER(apivot),
212 COMPLEX(rx), &m, &info
FCONE);
218 F77_CALL(dsytrs)(&aul, &m, &n, REAL(ax), &m, INTEGER(apivot),
219 REAL(rx), &m, &info
FCONE);
222 F77_CALL(dsptrs)(&aul, &m, &n, REAL(ax), INTEGER(apivot),
223 REAL(rx), &m, &info
FCONE);
226#ifdef MATRIX_ENABLE_ZMATRIX
246 XLENGTH(ax) == (R_xlen_t) m * m;
248 char rcl[] =
"...Matrix";
249 rcl[0] = (TYPEOF(ax) == CPLXSXP) ?
'z' :
'd';
255 rcl[2] = (unpacked) ?
'o' :
'p';
260 int *prdim = INTEGER(rdim);
265 char aul = CHAR(STRING_ELT(auplo, 0))[0];
266 if (isNull(b) && aul !=
'U') {
274 int info, pivoted = TYPEOF(aperm) == INTSXP && LENGTH(aperm) > 0;
278#ifdef MATRIX_ENABLE_ZMATRIX
279 if (TYPEOF(ax) == CPLXSXP) {
281 F77_CALL(zpotri)(&aul, &m, COMPLEX(rx), &m, &info
FCONE);
284 zsymperm2(COMPLEX(rx), n, aul, INTEGER(aperm), 1, 1);
286 F77_CALL(zpptri)(&aul, &m, COMPLEX(rx), &info
FCONE);
291 size_t lwork = (size_t) n * n;
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');
302 F77_CALL(dpotri)(&aul, &m, REAL(rx), &m, &info
FCONE);
305 dsymperm2( REAL(rx), n, aul, INTEGER(aperm), 1, 1);
307 F77_CALL(dpptri)(&aul, &m, REAL(rx), &info
FCONE);
312 size_t lwork = (size_t) n * n;
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');
320#ifdef MATRIX_ENABLE_ZMATRIX
328#ifdef MATRIX_ENABLE_ZMATRIX
329 if (TYPEOF(ax) == CPLXSXP) {
331 zrowperm2(COMPLEX(rx), m, n, INTEGER(aperm), 1, 0);
333 F77_CALL(zpotrs)(&aul, &m, &n, COMPLEX(ax), &m,
334 COMPLEX(rx), &m, &info
FCONE);
337 F77_CALL(zpptrs)(&aul, &m, &n, COMPLEX(ax),
338 COMPLEX(rx), &m, &info
FCONE);
342 zrowperm2(COMPLEX(rx), m, n, INTEGER(aperm), 1, 1);
346 drowperm2( REAL(rx), m, n, INTEGER(aperm), 1, 0);
348 F77_CALL(dpotrs)(&aul, &m, &n, REAL(ax), &m,
349 REAL(rx), &m, &info
FCONE);
352 F77_CALL(dpptrs)(&aul, &m, &n, REAL(ax),
353 REAL(rx), &m, &info
FCONE);
357 drowperm2( REAL(rx), m, n, INTEGER(aperm), 1, 1);
358#ifdef MATRIX_ENABLE_ZMATRIX
378 XLENGTH(ax) == (R_xlen_t) m * m;
380 char rcl[] =
"...Matrix";
381 rcl[0] = (TYPEOF(ax) == CPLXSXP) ?
'z' :
'd';
387 rcl[2] = (unpacked) ?
'r' :
'p';
392 int *prdim = INTEGER(rdim);
397 char aul = CHAR(STRING_ELT(auplo, 0))[0];
398 if (isNull(b) && aul !=
'U') {
405 char adi = CHAR(STRING_ELT(adiag, 0))[0];
406 if (isNull(b) && adi !=
'N') {
418#ifdef MATRIX_ENABLE_ZMATRIX
419 if (TYPEOF(ax) == CPLXSXP) {
421 F77_CALL(ztrtri)(&aul, &adi, &m, COMPLEX(rx), &m,
425 F77_CALL(ztptri)(&aul, &adi, &m, COMPLEX(rx),
432 F77_CALL(dtrtri)(&aul, &adi, &m, REAL(rx), &m,
436 F77_CALL(dtptri)(&aul, &adi, &m, REAL(rx),
440#ifdef MATRIX_ENABLE_ZMATRIX
448#ifdef MATRIX_ENABLE_ZMATRIX
449 if (TYPEOF(ax) == CPLXSXP) {
451 F77_CALL(ztrtrs)(&aul,
"N", &adi, &m, &n, COMPLEX(ax), &m,
455 F77_CALL(ztptrs)(&aul,
"N", &adi, &m, &n, COMPLEX(ax),
462 F77_CALL(dtrtrs)(&aul,
"N", &adi, &m, &n, REAL(ax), &m,
467 F77_CALL(dtptrs)(&aul,
"N", &adi, &m, &n, REAL(ax),
476#ifdef MATRIX_ENABLE_ZMATRIX
713 static const char *
valid[] = {
714 "A",
"LDLt",
"LD",
"DLt",
"L",
"Lt",
"D",
"P",
"Pt",
"" };
716 if (TYPEOF(system) != STRSXP || LENGTH(system) < 1 ||
717 (system = STRING_ELT(system, 0)) == NA_STRING ||
719 error(
_(
"invalid '%s' to '%s'"),
"system", __func__);
725 cholmod_factor *L =
M2CHF(a, 1);
726 if (!asLogical(sparse)) {
727 cholmod_dense *B = NULL, *X = NULL;
729 B = cholmod_allocate_dense(m, n, m, L->xtype, &
c);
732 R_xlen_t m1a = (R_xlen_t) m + 1;
734#define EYE(_CTYPE_, _ONE_) \
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) { \
744#ifdef MATRIX_ENABLE_ZMATRIX
745 if (L->xtype == CHOLMOD_COMPLEX)
753 X = cholmod_solve(ivalid, L, B, &
c);
754 cholmod_free_dense(&B, &
c);
757 PROTECT(r =
CHD2M(X, 0,
758 (ivalid < 2) ?
'p' : ((ivalid < 7) ?
't' :
'g')));
761 X = cholmod_solve(ivalid, L, B, &
c);
764 PROTECT(r =
CHD2M(X, 0,
'g'));
766 cholmod_free_dense(&X, &
c);
768 cholmod_sparse *B = NULL, *X = NULL;
770 B = cholmod_speye(m, n, L->xtype, &
c);
773 X = cholmod_spsolve(ivalid, L, B, &
c);
774 cholmod_free_sparse(&B, &
c);
775 if (X && ivalid < 7) {
779 (ivalid == 2 || ivalid == 4) ? -1 : 1, 1, &
c);
780 cholmod_free_sparse(&X, &
c);
785 PROTECT(r =
CHS2M(X, 1,
786 (ivalid < 2) ?
's' : ((ivalid < 7) ?
't' :
'g')));
789 X = cholmod_spsolve(ivalid, L, B, &
c);
792 PROTECT(r =
CHS2M(X, 1,
'g'));
794 cholmod_free_sparse(&X, &
c);
796 if (isNull(b) && (ivalid == 2 || ivalid == 4)) {
797 SEXP uplo = PROTECT(mkString(
"L"));
941 double *pbeta = REAL(beta);
944 int *pp = (LENGTH(p) > 0) ? INTEGER(p) : NULL;
946 int m = V_->
m, r = V_->
n, n, i, j, op_ = asInteger(op), nprotect = 5;
950 n = (asLogical(complete)) ? m : r;
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));
955#define EYE(_CTYPE_, _PTR_, _ONE_) \
957 _CTYPE_ *pyx = _PTR_(yx); \
958 Matrix_memset(pyx, 0, mn, sizeof(_CTYPE_)); \
959 if (isNull(yxjj)) { \
960 for (j = 0; j < n; ++j) { \
964 } else if (TYPEOF(yxjj) == TYPEOF(yx) && XLENGTH(yxjj) >= n) { \
965 _CTYPE_ *pyxjj = _PTR_(yxjj); \
966 for (j = 0; j < n; ++j) { \
972 error(_("invalid '%s' to '%s'"), "yxjj", __func__); \
975#ifdef MATRIX_ENABLE_ZMATRIX
980 EYE(
double, REAL, 1.0);
986 int *pydim = INTEGER(ydim);
988 error(
_(
"dimensions of '%s' and '%s' are inconsistent"),
995 char acl[] =
".geMatrix";
1000 int *padim = INTEGER(adim);
1001 padim[0] = (op_ != 0) ? m : r;
1005 if (isNull(y) && padim[0] == m)
1008 R_xlen_t mn = (R_xlen_t) padim[0] * padim[1];
1009 PROTECT(ax = allocVector(
IF_COMPLEX(CPLXSXP, REALSXP), mn));
1013#define MATMULT(_CTYPE_, _PTR_) \
1015 _CTYPE_ *pyx = _PTR_(yx), *pax = _PTR_(ax), *work = NULL; \
1017 work = (_CTYPE_ *) R_alloc((size_t) m, sizeof(_CTYPE_)); \
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); \
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); \
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); \
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); \
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); \
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); \
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); \
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); \
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); \
1104 error(_("invalid '%s' to '%s'"), "op", __func__); \
1109#ifdef MATRIX_ENABLE_ZMATRIX
1117 UNPROTECT(nprotect);