Matrix r5059
Loading...
Searching...
No Matches
forceSymmetric.c
Go to the documentation of this file.
1/* C implementation of methods for forceSymmetric */
2
3#include "Mdefines.h"
4#include "M5.h"
5#include "idz.h"
6
7SEXP dense_force_symmetric(SEXP from, const char *class,
8 char op_ul, char op_ct)
9{
10 if (class[0] != 'z')
11 op_ct = '\0';
12
13 char
14 ul0 = '\0', ul1 = 'U',
15 ct0 = '\0', ct1 = (class[0] == 'z') ? 'C' : '\0',
16 nu0 = '\0';
17 if (class[1] != 'g')
18 ul0 = ul1 = UPLO(from);
19 if (class[1] == 's' && class[0] == 'z')
20 ct0 = ct1 = TRANS(from);
21 if (class[1] == 't')
22 nu0 = DIAG(from);
23 if (op_ul != '\0')
24 ul1 = op_ul;
25 if (op_ct != '\0')
26 ct1 = op_ct;
27
28 if (class[1] == 's' && ul0 == ul1 && ct0 == ct1)
29 return from;
30
31 int packed = class[2] == 'p';
32
33 char cl[] = ".s.Matrix";
34 cl[0] = class[0];
35 cl[2] = (packed) ? 'p' : 'y';
36 SEXP to = PROTECT(newObject(cl));
37
38 int *pdim = DIM(from), n = pdim[1];
39 if (pdim[0] != n)
40 Rf_error(_("attempt to symmetrize a non-square matrix"));
41 SET_DIM(to, n, n);
42 SET_DIMNAMES(to, -(class[1] != 's'), DIMNAMES(from, 0));
43 if (ul1 != 'U')
44 SET_UPLO(to);
45 if (ct1 != 'C' && ct1 != '\0')
46 SET_TRANS(to);
47
48 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym));
49
50 if ((class[1] == 'g') ||
51 (class[1] == 's' && ul0 == ul1 && ct0 != 'C') ||
52 (class[1] == 't' && ul0 == ul1 && nu0 == 'N'))
53 SET_SLOT(to, Matrix_xSym, x0);
54 else {
55 SEXP x1 = PROTECT(Rf_allocVector(TYPEOF(x0), XLENGTH(x0)));
56 size_t n_ = (size_t) n;
57
58#define TEMPLATE(c) \
59 do { \
60 c##TYPE *px0 = c##PTR(x0), *px1 = c##PTR(x1); \
61 if (!packed) \
62 c##NAME(force2)(px1, px0, n_, ul0, ct0, nu0); \
63 else if (ul0 == ul1) \
64 c##NAME(force1)(px1, px0, n_, ul0, ct0, nu0); \
65 else { \
66 c##NAME(trans1)(px1, px0, n_, ul0, ct0); \
67 c##NAME(force1)(px1, NULL, n_, ul1, ct0, nu0); \
68 } \
69 } while (0)
70
71 SWITCH4(class[0], TEMPLATE);
72
73#undef TEMPLATE
74
75 SET_SLOT(to, Matrix_xSym, x1);
76 UNPROTECT(1); /* x1 */
77 }
78
79 UNPROTECT(2); /* x0, to */
80 return to;
81}
82
83SEXP sparse_force_symmetric(SEXP from, const char *class,
84 char op_ul, char op_ct)
85{
86 if (class[0] != 'z')
87 op_ct = '\0';
88
89 char
90 ul0 = '\0', ul1 = 'U',
91 ct0 = '\0', ct1 = (class[0] == 'z') ? 'C' : '\0',
92 nu0 = '\0';
93 if (class[1] != 'g')
94 ul0 = ul1 = UPLO(from);
95 if (class[1] == 's' && class[0] == 'z')
96 ct0 = ct1 = TRANS(from);
97 if (class[1] == 't')
98 nu0 = DIAG(from);
99 if (op_ul != '\0')
100 ul1 = op_ul;
101 if (op_ct != '\0')
102 ct1 = op_ct;
103
104 if (class[1] == 's' && ul0 == ul1 && ct0 == ct1)
105 return from;
106
107 char cl[] = ".s.Matrix";
108 cl[0] = class[0];
109 cl[2] = class[2];
110 SEXP to = PROTECT(newObject(cl));
111
112 int *pdim = DIM(from), n = pdim[1];
113 if (pdim[0] != n)
114 Rf_error(_("attempt to symmetrize a non-square matrix"));
115 SET_DIM(to, n, n);
116 SET_DIMNAMES(to, -(class[1] != 's'), DIMNAMES(from, 0));
117 if (ul1 != 'U')
118 SET_UPLO(to);
119 if (ct1 != 'C' && ct1 != '\0')
120 SET_TRANS(to);
121
122 int up = (class[0] != 'R') == (ul1 == 'U');
123
124 if (class[2] != 'T') {
125
126 SEXP iSym = (class[2] == 'C') ? Matrix_iSym : Matrix_jSym,
127 p0 = PROTECT(GET_SLOT(from, Matrix_pSym)),
128 i0 = PROTECT(GET_SLOT(from, iSym));
129 int *pp0 = INTEGER(p0), *pi0 = INTEGER(i0),
130 j, k, kend, nnz0 = pp0[n], nnz1 = 0;
131 pp0++;
132
133 if ((class[1] == 's' && ul0 == ul1) ||
134 (class[1] == 't' && ul0 == ul1 && nu0 == 'N')) {
135
136 SET_SLOT(to, Matrix_pSym, p0);
137 SET_SLOT(to, iSym, i0);
138
139 if (class[0] != 'n') {
140 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym));
141 if (class[1] != 's' || ct0 != 'C')
142 SET_SLOT(to, Matrix_xSym, x0);
143 else {
144 SEXP x1 = PROTECT(Rf_allocVector(CPLXSXP, nnz0));
145 Rcomplex *px0 = COMPLEX(x0), *px1 = COMPLEX(x1);
146 for (j = 0, k = 0; j < n; ++j) {
147 kend = pp0[j];
148 while (k < kend) {
149 if (pi0[k] == j)
150 zASSIGN_PROJ_REAL(px1[k], px0[k]);
151 else
152 zASSIGN_IDEN (px1[k], px0[k]);
153 ++k;
154 }
155 }
156 SET_SLOT(to, Matrix_xSym, x1);
157 UNPROTECT(1); /* x1 */
158 }
159 UNPROTECT(1); /* x0 */
160 }
161
162 } else {
163
164 SEXP p1 = PROTECT(Rf_allocVector(INTSXP, XLENGTH(p0)));
165 int *pp1 = INTEGER(p1);
166 *(pp1++) = 0;
167 SET_SLOT(to, Matrix_pSym, p1);
168
169 if (class[1] == 'g')
170 for (j = 0, k = 0; j < n; ++j) {
171 kend = pp0[j];
172 while (k < kend) {
173 if ((up) ? pi0[k] <= j : pi0[k] >= j)
174 ++nnz1;
175 ++k;
176 }
177 pp1[j] = nnz1;
178 }
179 else if (class[1] == 's')
180 nnz1 = nnz0;
181 else if (nu0 == 'N')
182 for (j = 0, k = 0; j < n; ++j) {
183 kend = pp0[j];
184 if (k < kend && pi0[(up) ? k : kend - 1] == j)
185 ++nnz1;
186 pp1[j] = nnz1;
187 k = kend;
188 }
189 else if (ul0 == ul1) {
190 for (j = 0; j < n; ++j)
191 pp1[j] = ++nnz1 + pp0[j];
192 nnz1 = nnz0 + n;
193 }
194 else
195 for (j = 0; j < n; ++j)
196 pp1[j] = ++nnz1;
197
198 SEXP i1 = PROTECT(Rf_allocVector(INTSXP, nnz1));
199 int *pi1 = INTEGER(i1);
200 SET_SLOT(to, iSym, i1);
201
202#define TEMPLATE(c) \
203 do { \
204 c##TYPE *px0 = NULL, *px1 = NULL; \
205 c##IF_NPATTERN( \
206 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym)), \
207 x1 = PROTECT(Rf_allocVector(c##TYPESXP, nnz1)); \
208 px0 = c##PTR(x0); \
209 px1 = c##PTR(x1); \
210 SET_SLOT(to, Matrix_xSym, x1); \
211 UNPROTECT(2); /* x1, x0 */ \
212 ); \
213 if (class[1] == 'g') \
214 for (j = 0, k = 0; j < n; ++j) { \
215 kend = pp0[j]; \
216 while (k < kend) { \
217 if ((up) ? pi0[k] <= j : pi0[k] >= j) { \
218 *(pi1++) = pi0[k]; \
219 c##IF_NPATTERN( \
220 *(px1++) = px0[k]; \
221 ); \
222 } \
223 ++k; \
224 } \
225 } \
226 else if (class[1] == 's') { \
227 int *iwork = NULL; \
228 Matrix_Calloc(iwork, n, int); \
229 c##csptrans(pp1, pi1, px1, pp0, pi0, px0, n, n, ct0, iwork); \
230 Matrix_Free(iwork, n); \
231 if (ct0 == 'C') \
232 for (j = 0, k = 0; j < n; ++j) { \
233 kend = pp1[j]; \
234 if (k < kend && pi1[(up) ? kend - 1 : k] == j) \
235 c##SET_PROJ_REAL(px1[(up) ? kend - 1 : k]); \
236 k = kend; \
237 } \
238 } \
239 else if (nu0 == 'N') \
240 for (j = 0, k = 0; j < n; ++j) { \
241 kend = pp0[j]; \
242 if (k < kend && pi0[(up) ? k : kend - 1] == j) { \
243 *(pi1++) = j; \
244 c##IF_NPATTERN( \
245 *(px1++) = px0[(up) ? k : kend - 1]; \
246 ); \
247 } \
248 k = kend; \
249 } \
250 else if (ul0 == ul1) \
251 for (j = 0, k = 0; j < n; ++j) { \
252 kend = pp0[j]; \
253 if (!up) { \
254 *(pi1++) = j; \
255 c##IF_NPATTERN( \
256 *(px1++) = c##UNIT; \
257 ); \
258 } \
259 while (k < kend) { \
260 *(pi1++) = pi0[k]; \
261 c##IF_NPATTERN( \
262 *(px1++) = px0[k]; \
263 ); \
264 ++k; \
265 } \
266 if (up) { \
267 *(pi1++) = j; \
268 c##IF_NPATTERN( \
269 *(px1++) = c##UNIT; \
270 ); \
271 } \
272 } \
273 else \
274 for (j = 0; j < n; ++j) { \
275 *(pi1++) = j; \
276 c##IF_NPATTERN( \
277 *(px1++) = c##UNIT; \
278 ); \
279 } \
280 } while (0)
281
282 SWITCH5(class[0], TEMPLATE);
283
284#undef TEMPLATE
285
286 UNPROTECT(2); /* i1, p1 */
287
288 }
289
290 UNPROTECT(2); /* i0, p0 */
291
292 } else {
293
294 SEXP i0 = PROTECT(GET_SLOT(from, Matrix_iSym)),
295 j0 = PROTECT(GET_SLOT(from, Matrix_jSym));
296 int *pi0 = INTEGER(i0), *pj0 = INTEGER(j0);
297 R_xlen_t k, nnz0 = XLENGTH(i0), nnz1 = 0;
298
299 if ((class[1] == 's') ||
300 (class[1] == 't' && ul0 == ul1 && nu0 == 'N')) {
301
302 if (class[1] != 's' || ul0 == ul1) {
303 SET_SLOT(to, Matrix_iSym, i0);
304 SET_SLOT(to, Matrix_jSym, j0);
305 } else {
306 SET_SLOT(to, Matrix_iSym, j0);
307 SET_SLOT(to, Matrix_jSym, i0);
308 }
309
310 if (class[0] != 'n') {
311 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym));
312 if (class[1] != 's' || ct0 != 'C')
313 SET_SLOT(to, Matrix_xSym, x0);
314 else {
315 SEXP x1 = PROTECT(Rf_allocVector(CPLXSXP, XLENGTH(x0)));
316 Rcomplex *px0 = COMPLEX(x0), *px1 = COMPLEX(x1);
317 for (k = 0; k < nnz0; ++k)
318 if (pi0[k] == pj0[k])
319 zASSIGN_PROJ_REAL(px1[k], px0[k]);
320 else if (ul0 == ul1)
321 zASSIGN_IDEN (px1[k], px0[k]);
322 else
323 zASSIGN_CONJ (px1[k], px0[k]);
324 SET_SLOT(to, Matrix_xSym, x1);
325 UNPROTECT(1); /* x1 */
326 }
327 UNPROTECT(1); /* x0 */
328 }
329
330 } else {
331
332 if (class[1] == 'g' || nu0 == 'N') {
333 for (k = 0; k < nnz0; ++k)
334 if ((up) ? pi0[k] <= pj0[k] : pi0[k] >= pj0[k])
335 ++nnz1;
336 }
337 else
338 nnz1 = (ul0 == ul1) ? n + nnz0 : n;
339
340 SEXP i1 = PROTECT(Rf_allocVector(INTSXP, nnz1)),
341 j1 = PROTECT(Rf_allocVector(INTSXP, nnz1));
342 int *pi1 = INTEGER(i1), *pj1 = INTEGER(j1), j;
343 SET_SLOT(to, Matrix_iSym, i1);
344 SET_SLOT(to, Matrix_jSym, j1);
345
346#define TEMPLATE(c) \
347 do { \
348 c##IF_NPATTERN( \
349 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym)), \
350 x1 = PROTECT(Rf_allocVector(c##TYPESXP, nnz1)); \
351 c##TYPE *px0 = c##PTR(x0), *px1 = c##PTR(x1); \
352 SET_SLOT(to, Matrix_xSym, x1); \
353 UNPROTECT(2); /* x1, x0 */ \
354 ); \
355 if (class[1] == 'g' || nu0 == 'N') \
356 for (k = 0; k < nnz0; ++k) { \
357 if ((up) ? pi0[k] <= pj0[k] : pi0[k] >= pj0[k]) { \
358 *(pi1++) = pi0[k]; \
359 *(pj1++) = pj0[k]; \
360 c##IF_NPATTERN( \
361 *(px1++) = px0[k]; \
362 ); \
363 } \
364 } \
365 else { \
366 if (ul0 == ul1) { \
367 memcpy(pi1, pi0, sizeof( int) * (size_t) nnz0); \
368 memcpy(pj1, pj0, sizeof( int) * (size_t) nnz0); \
369 pi1 += nnz0; \
370 pj1 += nnz0; \
371 c##IF_NPATTERN( \
372 memcpy(px1, px0, sizeof(c##TYPE) * (size_t) nnz0); \
373 px1 += nnz0; \
374 ); \
375 } \
376 for (j = 0; j < n; ++j) { \
377 *(pi1++) = *(pj1++) = j; \
378 c##IF_NPATTERN( \
379 *(px1++) = c##UNIT; \
380 ); \
381 } \
382 } \
383 } while (0)
384
385 SWITCH5(class[0], TEMPLATE);
386
387#undef TEMPLATE
388
389 UNPROTECT(2); /* j1, i1 */
390
391 }
392
393 UNPROTECT(2); /* j0, i0 */
394
395 }
396
397 UNPROTECT(1); /* to */
398 return to;
399}
400
401SEXP R_dense_force_symmetric(SEXP s_from, SEXP s_uplo, SEXP s_trans)
402{
403 const char *class = Matrix_class(s_from, valid_dense, 6, __func__);
404
405 char ul = '\0', ct = '\0';
406 if (s_uplo != R_NilValue)
407 VALID_UPLO (s_uplo , ul);
408 VALID_TRANS(s_trans, ct);
409
410 return dense_force_symmetric(s_from, class, ul, ct);
411}
412
413SEXP R_sparse_force_symmetric(SEXP s_from, SEXP s_uplo, SEXP s_trans)
414{
415 const char *class = Matrix_class(s_from, valid_sparse, 6, __func__);
416
417 char ul = '\0', ct = '\0';
418 if (s_uplo != R_NilValue)
419 VALID_UPLO (s_uplo , ul);
420 VALID_TRANS(s_trans, ct);
421
422 return sparse_force_symmetric(s_from, class, ul, ct);
423}
zASSIGN_IDEN
#define zASSIGN_IDEN(x, y)
Definition M5.h:166
SWITCH5
#define SWITCH5(c, template)
Definition M5.h:327
zASSIGN_CONJ
#define zASSIGN_CONJ(x, y)
Definition M5.h:172
zASSIGN_PROJ_REAL
#define zASSIGN_PROJ_REAL(x, y)
Definition M5.h:178
SWITCH4
#define SWITCH4(c, template)
Definition M5.h:315
SET_DIM
#define SET_DIM(x, m, n)
Definition Mdefines.h:87
valid_dense
const char * valid_dense[]
Definition objects.c:3
_
#define _(String)
Definition Mdefines.h:66
SET_UPLO
#define SET_UPLO(x)
Definition Mdefines.h:103
DIAG
#define DIAG(x)
Definition Mdefines.h:111
UPLO
#define UPLO(x)
Definition Mdefines.h:101
SET_DIMNAMES
#define SET_DIMNAMES(x, mode, value)
Definition Mdefines.h:98
VALID_UPLO
#define VALID_UPLO(s, c)
Definition Mdefines.h:158
Matrix_class
const char * Matrix_class(SEXP, const char **, int, const char *)
Definition objects.c:112
DIMNAMES
#define DIMNAMES(x, mode)
Definition Mdefines.h:96
SET_TRANS
#define SET_TRANS(x)
Definition Mdefines.h:108
TRANS
#define TRANS(x)
Definition Mdefines.h:106
DIM
#define DIM(x)
Definition Mdefines.h:85
GET_SLOT
#define GET_SLOT(x, name)
Definition Mdefines.h:72
VALID_TRANS
#define VALID_TRANS(s, c)
Definition Mdefines.h:166
newObject
SEXP newObject(const char *)
Definition objects.c:13
valid_sparse
const char * valid_sparse[]
Definition Mdefines.h:328
SET_SLOT
#define SET_SLOT(x, name, value)
Definition Mdefines.h:73
TYPEOF
#define TYPEOF(s)
Definition Mdefines.h:123
cl
cholmod_common cl
Definition cholmod-etc.c:6
sparse_force_symmetric
SEXP sparse_force_symmetric(SEXP from, const char *class, char op_ul, char op_ct)
Definition forceSymmetric.c:83
TEMPLATE
#define TEMPLATE(c)
dense_force_symmetric
SEXP dense_force_symmetric(SEXP from, const char *class, char op_ul, char op_ct)
Definition forceSymmetric.c:7
R_sparse_force_symmetric
SEXP R_sparse_force_symmetric(SEXP s_from, SEXP s_uplo, SEXP s_trans)
Definition forceSymmetric.c:413
R_dense_force_symmetric
SEXP R_dense_force_symmetric(SEXP s_from, SEXP s_uplo, SEXP s_trans)
Definition forceSymmetric.c:401
Matrix_xSym
SEXP Matrix_xSym
Definition init.c:635
Matrix_iSym
SEXP Matrix_iSym
Definition init.c:607
Matrix_jSym
SEXP Matrix_jSym
Definition init.c:610
Matrix_pSym
SEXP Matrix_pSym
Definition init.c:622
Generated by doxygen 1.12.0