11 char cl[] =
"..CMatrix";
12 cl[0] = (
TYPEOF(x) == CPLXSXP) ?
'z' :
'd';
18 int n =
DIM(s_trf)[1], packed = XLENGTH(x) != (int_fast64_t) n * n;
23 char ul =
UPLO(s_trf),
34 R_xlen_t n1a = (R_xlen_t) n + 1;
37 D_p = PROTECT(Rf_allocVector(INTSXP, n1a));
38 int *ppivot = INTEGER(pivot), *D_pp = INTEGER(D_p),
39 b = n, dp = (ul ==
'U') ? 1 : 2;
44 D_pp[j+1] = D_pp[j] + 1;
47 D_pp[j+1] = D_pp[j] + dp;
48 D_pp[j+2] = D_pp[j] + 3;
56 SEXP P, P_perm, T, T_p, T_i, T_x,
57 D_i = PROTECT(Rf_allocVector(INTSXP, D_pp[n])),
58 D_x = PROTECT(Rf_allocVector(
TYPEOF(x), D_pp[n]));
59 int *P_pperm, *T_pp, *T_pi, *D_pi = INTEGER(D_i);
61 R_xlen_t len = (R_xlen_t) 2 * b + 1, k = (ul ==
'U') ? len - 2 : 0;
62 SEXP ans = PROTECT(Rf_allocVector(VECSXP, len));
66 c##TYPE *T_px, *D_px = c##PTR(D_x), *px = c##PTR(x); \
70 s = (ppivot[j] > 0) ? 1 : 2; \
71 dp = (ul == 'U') ? j : n - j - s; \
73 PROTECT(P = Rf_duplicate(P_)); \
74 PROTECT(P_perm = Rf_allocVector(INTSXP, n)); \
75 PROTECT(T = Rf_duplicate(T_)); \
76 PROTECT(T_p = Rf_allocVector(INTSXP, n1a)); \
77 PROTECT(T_i = Rf_allocVector(INTSXP, (R_xlen_t) s * dp)); \
78 PROTECT(T_x = Rf_allocVector(TYPEOF(x), (R_xlen_t) s * dp)); \
80 P_pperm = INTEGER(P_perm); \
81 T_pp = INTEGER(T_p); \
82 T_pi = INTEGER(T_i); \
86 for (i = 0; i < j; ++i) { \
90 for (i = j; i < j+s; ++i) { \
91 T_pp[i+1] = T_pp[i] + dp; \
94 for (i = j+s; i < n; ++i) { \
95 T_pp[i+1] = T_pp[i]; \
100 P_pperm[j] = ppivot[j]; \
101 P_pperm[ppivot[j]-1] = j + 1; \
102 } else if (ul == 'U') { \
103 P_pperm[j] = -ppivot[j]; \
104 P_pperm[-ppivot[j]-1] = j + 1; \
106 P_pperm[j+1] = -ppivot[j]; \
107 P_pperm[-ppivot[j]-1] = j + 2; \
111 for (i = 0; i < j; ++i) { \
113 *(T_px++) = *(px++); \
116 *(D_px++) = *(px++); \
121 for (i = 0; i < j-1; ++i) { \
123 *(T_px++) = *(px++); \
127 *(D_px++) = *(px++); \
128 *(D_px++) = *(px++); \
137 *(D_px++) = *(px++); \
138 *(D_px++) = *(px++); \
139 for (i = j+2; i < n; ++i) { \
141 *(T_px++) = *(px++); \
148 *(D_px++) = *(px++); \
149 for (i = j+1; i < n; ++i) { \
151 *(T_px++) = *(px++); \
158 SET_SLOT(P, Matrix_permSym, P_perm); \
159 SET_SLOT(T, Matrix_pSym, T_p); \
160 SET_SLOT(T, Matrix_iSym, T_i); \
161 SET_SLOT(T, Matrix_xSym, T_x); \
164 SET_VECTOR_ELT(ans, k-1, P); \
165 SET_VECTOR_ELT(ans, k , T); \
168 SET_VECTOR_ELT(ans, k , P); \
169 SET_VECTOR_ELT(ans, k+1, T); \
183 SET_VECTOR_ELT(ans, len-1, D_);