Matrix r5059
Loading...
Searching...
No Matches
expand.c
Go to the documentation of this file.
1/* C implementation of methods for expand */
2
3#include "Mdefines.h"
4#include "M5.h"
5
7{
8 SEXP x = PROTECT(GET_SLOT(s_trf, Matrix_xSym));
9
10 SEXP P_ = PROTECT(newObject("pMatrix"));
11 char cl[] = "..CMatrix";
12 cl[0] = (TYPEOF(x) == CPLXSXP) ? 'z' : 'd';
13 cl[1] = 't';
14 SEXP T_ = PROTECT(newObject(cl));
15 cl[1] = 's';
16 SEXP D_ = PROTECT(newObject(cl));
17
18 int n = DIM(s_trf)[1], packed = XLENGTH(x) != (int_fast64_t) n * n;
19 SET_DIM(P_, n, n);
20 SET_DIM(T_, n, n);
21 SET_DIM(D_, n, n);
22
23 char ul = UPLO(s_trf),
24 ct = (TYPEOF(x) == CPLXSXP) ? TRANS(s_trf) : 'C';
25 if (ul != 'U') {
26 SET_UPLO(T_);
27 SET_UPLO(D_);
28 }
29 if (ct != 'C')
30 SET_TRANS(D_);
31 SET_DIAG(T_);
32
33 int i, j, s;
34 R_xlen_t n1a = (R_xlen_t) n + 1;
35
36 SEXP pivot = PROTECT(GET_SLOT(s_trf, Matrix_permSym)),
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;
40 D_pp[0] = 0;
41 j = 0;
42 while (j < n) {
43 if (ppivot[j] > 0) {
44 D_pp[j+1] = D_pp[j] + 1;
45 j += 1;
46 } else {
47 D_pp[j+1] = D_pp[j] + dp;
48 D_pp[j+2] = D_pp[j] + 3;
49 j += 2;
50 b -= 1;
51 }
52 }
53 SET_SLOT(D_, Matrix_pSym, D_p);
54 UNPROTECT(1); /* D_p */
55
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);
60
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));
63
64#define TEMPLATE(c) \
65 do { \
66 c##TYPE *T_px, *D_px = c##PTR(D_x), *px = c##PTR(x); \
67 \
68 j = 0; \
69 while (b--) { \
70 s = (ppivot[j] > 0) ? 1 : 2; \
71 dp = (ul == 'U') ? j : n - j - s; \
72 \
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)); \
79 \
80 P_pperm = INTEGER(P_perm); \
81 T_pp = INTEGER(T_p); \
82 T_pi = INTEGER(T_i); \
83 T_px = c##PTR(T_x); \
84 T_pp[0] = 0; \
85 \
86 for (i = 0; i < j; ++i) { \
87 T_pp[i+1] = 0; \
88 P_pperm[i] = i + 1; \
89 } \
90 for (i = j; i < j+s; ++i) { \
91 T_pp[i+1] = T_pp[i] + dp; \
92 P_pperm[i] = i + 1; \
93 } \
94 for (i = j+s; i < n; ++i) { \
95 T_pp[i+1] = T_pp[i]; \
96 P_pperm[i] = i + 1; \
97 } \
98 \
99 if (s == 1) { \
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; \
105 } else { \
106 P_pperm[j+1] = -ppivot[j]; \
107 P_pperm[-ppivot[j]-1] = j + 2; \
108 } \
109 \
110 if (ul == 'U') { \
111 for (i = 0; i < j; ++i) { \
112 *(T_pi++) = i; \
113 *(T_px++) = *(px++); \
114 } \
115 *(D_pi++) = j; \
116 *(D_px++) = *(px++); \
117 ++j; \
118 if (!packed) \
119 px += n - j; \
120 if (s == 2) { \
121 for (i = 0; i < j-1; ++i) { \
122 *(T_pi++) = i; \
123 *(T_px++) = *(px++); \
124 } \
125 *(D_pi++) = j - 1; \
126 *(D_pi++) = j; \
127 *(D_px++) = *(px++); \
128 *(D_px++) = *(px++); \
129 ++j; \
130 if (!packed) \
131 px += n - j; \
132 } \
133 } else { \
134 if (s == 2) { \
135 *(D_pi++) = j; \
136 *(D_pi++) = j + 1; \
137 *(D_px++) = *(px++); \
138 *(D_px++) = *(px++); \
139 for (i = j+2; i < n; ++i) { \
140 *(T_pi++) = i; \
141 *(T_px++) = *(px++); \
142 } \
143 ++j; \
144 if (!packed) \
145 px += j; \
146 } \
147 *(D_pi++) = j; \
148 *(D_px++) = *(px++); \
149 for (i = j+1; i < n; ++i) { \
150 *(T_pi++) = i; \
151 *(T_px++) = *(px++); \
152 } \
153 ++j; \
154 if (!packed) \
155 px += j; \
156 } \
157 \
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); \
162 \
163 if (ul == 'U') { \
164 SET_VECTOR_ELT(ans, k-1, P); \
165 SET_VECTOR_ELT(ans, k , T); \
166 k -= 2; \
167 } else { \
168 SET_VECTOR_ELT(ans, k , P); \
169 SET_VECTOR_ELT(ans, k+1, T); \
170 k += 2; \
171 } \
172 UNPROTECT(6); /* T_x, T_i, T_p, T, P_perm, P */ \
173 } \
174 } while (0)
175
176 if (TYPEOF(x) == CPLXSXP)
177 TEMPLATE(z);
178 else
179 TEMPLATE(d);
180
181 SET_SLOT(D_, Matrix_iSym, D_i);
182 SET_SLOT(D_, Matrix_xSym, D_x);
183 SET_VECTOR_ELT(ans, len-1, D_);
184
185 UNPROTECT(8); /* ans, D_x, D_i, pivot, D_, T_, P_, x */
186 return ans;
187}
#define SET_DIM(x, m, n)
Definition Mdefines.h:87
#define SET_UPLO(x)
Definition Mdefines.h:103
#define UPLO(x)
Definition Mdefines.h:101
#define SET_TRANS(x)
Definition Mdefines.h:108
#define TRANS(x)
Definition Mdefines.h:106
#define SET_DIAG(x)
Definition Mdefines.h:113
#define DIM(x)
Definition Mdefines.h:85
#define GET_SLOT(x, name)
Definition Mdefines.h:72
SEXP newObject(const char *)
Definition objects.c:13
#define SET_SLOT(x, name, value)
Definition Mdefines.h:73
#define TYPEOF(s)
Definition Mdefines.h:123
cholmod_common cl
Definition cholmod-etc.c:6
#define TEMPLATE(c)
SEXP denseBunchKaufman_expand(SEXP s_trf)
Definition expand.c:6
SEXP Matrix_permSym
Definition init.c:623
SEXP Matrix_xSym
Definition init.c:635
SEXP Matrix_iSym
Definition init.c:607
SEXP Matrix_pSym
Definition init.c:622