Matrix r5059
Loading...
Searching...
No Matches
skewpart.c
Go to the documentation of this file.
1/* C implementation of methods for skewpart */
2
3#include "Mdefines.h"
4#include "M5.h"
5#include "idz.h"
6
7/* defined in ./coerce.c : */
8SEXP dense_as_kind(SEXP, const char *, char, int);
9SEXP sparse_as_kind(SEXP, const char *, char);
10
11SEXP dense_skewpart(SEXP from, const char *class, char op_ct)
12{
13 PROTECT(from = dense_as_kind(from, class, ',', 0));
14
15 if (class[0] != 'z')
16 op_ct = '\0';
17
18 int packed = class[2] == 'p';
19
20 char cl[] = "...Matrix";
21 cl[0] = (class[0] == 'z') ? 'z' : 'd';
22 cl[1] = (class[1] == 's') ? 's' : 'g';
23 cl[2] = (class[1] == 's') ? class[2] : 'e';
24 SEXP to = PROTECT(newObject(cl));
25
26 int *pdim = DIM(from), n = pdim[1];
27 if (pdim[0] != n)
28 Rf_error((op_ct == 'C')
29 ? _("attempt to get skew-Hermitian part of non-square matrix")
30 : _("attempt to get skew-symmetric part of non-square matrix"));
31 SET_DIM(to, n, n);
32 SET_DIMNAMES(to, -(class[1] != 's'), DIMNAMES(from, 0));
33
34 char ul = '\0', ct = '\0', nu = '\0';
35 if (class[1] != 'g' && (ul = UPLO(from)) != 'U' && class[1] == 's')
36 SET_UPLO(to);
37 if (class[1] == 's' && class[0] == 'z' && (ct = TRANS(from)) != 'C')
38 SET_TRANS(to);
39 if (class[1] == 't' && (nu = DIAG(from)) != 'N')
40 ;
41
42 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym)), x1;
43
44 if (class[1] == 's') {
45
46 /* Skew-symmetric part of Hermitian matrix is imaginary part */
47 /* Skew-Hermitian part of symmetric matrix is imaginary part */
48 if (op_ct == ct)
49 PROTECT(x1 = allocZero(TYPEOF(x0), XLENGTH(x0)));
50 else {
51 PROTECT(x1 = Rf_allocVector(CPLXSXP, XLENGTH(x0)));
52 zvimag(COMPLEX(x1), COMPLEX(x0), (size_t) XLENGTH(x0));
53 }
54
55 } else {
56
57 if ((int_fast64_t) n * n > R_XLEN_T_MAX)
58 Rf_error(_("attempt to allocate vector of length exceeding %s"),
59 "R_XLEN_T_MAX");
60 PROTECT(x1 = Rf_allocVector(TYPEOF(x0), (R_xlen_t) n * n));
61
62 int i, j;
63
64#define TEMPLATE(c) \
65 do { \
66 c##TYPE *px0 = c##PTR(x0), *pu0 = px0, *pl0 = px0; \
67 c##TYPE *px1 = c##PTR(x1), *pu1 = px1, *pl1 = px1; \
68 if (class[1] == 'g') { \
69 if (op_ct == 'C') \
70 for (j = 0; j < n; ++j) { \
71 for (i = 0; i < j; ++i) { \
72 c##ASSIGN_IDEN(*pu1, *pu0); \
73 c##DECREMENT_CONJ(*pu1, *pl0); \
74 c##ASSIGN_CONJ(*pl1, *pu1); \
75 c##MULTIPLY(*pu1, 0.5); \
76 c##MULTIPLY(*pl1, -0.5); \
77 pu0 += 1; \
78 pu1 += 1; \
79 pl0 += n; \
80 pl1 += n; \
81 } \
82 c##ASSIGN_PROJ_IMAG(*pu1, *pu0); \
83 pu0 += n - j; \
84 pu1 += n - j; \
85 pl0 = px0 + j + 1; \
86 pl1 = px1 + j + 1; \
87 } \
88 else \
89 for (j = 0; j < n; ++j) { \
90 for (i = 0; i < j; ++i) { \
91 c##ASSIGN_IDEN(*pu1, *pu0); \
92 c##DECREMENT_IDEN(*pu1, *pl0); \
93 c##ASSIGN_IDEN(*pl1, *pu1); \
94 c##MULTIPLY(*pu1, 0.5); \
95 c##MULTIPLY(*pl1, -0.5); \
96 pu0 += 1; \
97 pu1 += 1; \
98 pl0 += n; \
99 pl1 += n; \
100 } \
101 c##SET_ZERO(*pu1); \
102 pu0 += n - j; \
103 pu1 += n - j; \
104 pl0 = px0 + j + 1; \
105 pl1 = px1 + j + 1; \
106 } \
107 } else { \
108 if (ul == 'U') \
109 for (j = 0; j < n; ++j) { \
110 for (i = 0; i < j; ++i) { \
111 c##ASSIGN_IDEN(*pu1, *pu0); \
112 c##ASSIGN_IDEN(*pl1, *pu0); \
113 c##MULTIPLY(*pu1, 0.5); \
114 c##MULTIPLY(*pl1, -0.5); \
115 pu0 += 1; \
116 pu1 += 1; \
117 pl1 += n; \
118 } \
119 if (nu != 'N' || op_ct != 'C') \
120 c##SET_ZERO(*pu1); \
121 else \
122 c##ASSIGN_PROJ_IMAG(*pu1, *pu0); \
123 pu0 += 1; \
124 pu1 += 1; \
125 pl1 += n; \
126 if (!packed) \
127 pu0 += n - j - 1; \
128 pu1 += n - j - 1; \
129 pl1 = px1 + j + 1; \
130 } \
131 else \
132 for (j = 0; j < n; ++j) { \
133 if (!packed) \
134 pl0 += j; \
135 pl1 += j; \
136 pu1 = pl1; \
137 if (nu != 'N' || op_ct != 'C') \
138 c##SET_ZERO(*pl1); \
139 else \
140 c##ASSIGN_PROJ_IMAG(*pl1, *pl0); \
141 pl0 += 1; \
142 pl1 += 1; \
143 pu1 += n; \
144 for (i = j + 1; i < n; ++i) { \
145 c##ASSIGN_IDEN(*pl1, *pl0); \
146 c##ASSIGN_IDEN(*pu1, *pl0); \
147 c##MULTIPLY(*pl1, 0.5); \
148 c##MULTIPLY(*pu1, -0.5); \
149 pl0 += 1; \
150 pl1 += 1; \
151 pu1 += n; \
152 } \
153 } \
154 } \
155 } while (0)
156
157 SWITCH2(class[0], TEMPLATE);
158
159#undef TEMPLATE
160
161 }
162
163 SET_SLOT(to, Matrix_xSym, x1);
164
165 UNPROTECT(4); /* x1, x0, to, from */
166 return to;
167}
168
169SEXP sparse_skewpart(SEXP from, const char *class, char op_ct)
170{
171 PROTECT(from = sparse_as_kind(from, class, ','));
172
173 if (class[0] != 'z')
174 op_ct = '\0';
175
176 char cl[] = "...Matrix";
177 cl[0] = (class[0] == 'z') ? 'z' : 'd';
178 cl[1] = (class[1] == 's') ? 's' : 'g';
179 cl[2] = class[2];
180 SEXP to = PROTECT(newObject(cl));
181
182 int *pdim = DIM(from), n = pdim[1];
183 if (pdim[0] != n)
184 Rf_error((op_ct == 'C')
185 ? _("attempt to get skew-Hermitian part of non-square matrix")
186 : _("attempt to get skew-symmetric part of non-square matrix"));
187 SET_DIM(to, n, n);
188 SET_DIMNAMES(to, -(class[1] != 's'), DIMNAMES(from, 0));
189
190 char ul = '\0', ct = '\0';
191 if (class[1] != 'g' && (ul = UPLO(from)) != 'U' && class[1] == 's')
192 SET_UPLO(to);
193 if (class[1] == 's' && class[0] == 'z' && (ct = TRANS(from)) != 'C')
194 SET_TRANS(to);
195
196 if (class[1] == 's' && op_ct == ct) {
197 if (class[2] != 'T') {
198 SEXP p1 = PROTECT(allocZero(INTSXP, (R_xlen_t) n + 1));
199 SET_SLOT(to, Matrix_pSym, p1);
200 UNPROTECT(1); /* p1 */
201 }
202 UNPROTECT(2); /* to, from */
203 return to;
204 }
205
206 if (class[2] != 'T') {
207
208 SEXP iSym = (class[2] == 'C') ? Matrix_iSym : Matrix_jSym,
209 p0 = PROTECT(GET_SLOT(from, Matrix_pSym)),
210 i0 = PROTECT(GET_SLOT(from, iSym)),
211 x0 = PROTECT(GET_SLOT(from, Matrix_xSym));
212 int *pp0 = INTEGER(p0), *pi0 = INTEGER(i0),
213 j, k, kend, nnz0 = pp0[n], nnz1;
214 pp0++;
215
216 if (class[1] == 's') {
217
218 SET_SLOT(to, Matrix_pSym, p0);
219 SET_SLOT(to, iSym, i0);
220
221 /* Skew-symmetric part of Hermitian matrix is imaginary part */
222 /* Skew-Hermitian part of symmetric matrix is imaginary part */
223 SEXP x1 = PROTECT(Rf_allocVector(CPLXSXP, nnz0));
224 zvimag(COMPLEX(x1), COMPLEX(x0), (size_t) nnz0);
225 SET_SLOT(to, Matrix_xSym, x1);
226 UNPROTECT(1); /* x1 */
227
228 } else {
229
230 int *iwork = NULL;
231 size_t liwork = (size_t) ((int_fast64_t) n + n + 1 + nnz0);
232 Matrix_Calloc(iwork, liwork, int);
233
234 int *pp0_ = iwork + n + 1, *pi0_ = iwork + n + n + 1;
235 ncsptrans(pp0_ - 1, pi0_, NULL, pp0 - 1, pi0, NULL, n, n, 'T', iwork);
236 memcpy(iwork, pp0 - 1, sizeof(int) * (size_t) n);
237
238 SEXP p1 = PROTECT(Rf_allocVector(INTSXP, XLENGTH(p0)));
239 int *pp1 = INTEGER(p1), k_, kend_;
240 *(pp1++) = 0;
241 SET_SLOT(to, Matrix_pSym, p1);
242
243 for (j = 0, k = 0, k_ = 0; j < n; ++j) {
244 kend = pp0 [j];
245 kend_ = pp0_[j];
246 pp1[j] = 0;
247 while (k < kend) {
248 if (pi0[k] >= j)
249 k = kend;
250 else {
251 while (k_ < kend_ && pi0_[k_] < pi0[k]) {
252 ++pp1[j];
253 ++pp1[pi0_[k_]];
254 ++k_;
255 }
256 ++pp1[j];
257 ++pp1[pi0[k]];
258 if (k_ < kend_ && pi0_[k_] == pi0[k])
259 ++k_;
260 ++k;
261 }
262 }
263 while (k_ < kend_) {
264 if (pi0_[k_] >= j)
265 k_ = kend_;
266 else {
267 ++pp1[j];
268 ++pp1[pi0_[k_]];
269 ++k_;
270 }
271 }
272 }
273 for (j = 0; j < n; ++j)
274 pp1[j] += pp1[j - 1];
275 nnz1 = pp1[n - 1];
276 for (j = n - 1; j >= 0; --j)
277 pp1[j] = pp1[j - 1];
278
279 SEXP i1 = PROTECT(Rf_allocVector(INTSXP, nnz1)),
280 x1 = PROTECT(Rf_allocVector(TYPEOF(x0), nnz1));
281 int *pi1 = INTEGER(i1), l;
282 SET_SLOT(to, iSym, i1);
283 SET_SLOT(to, Matrix_xSym, x1);
284
285#define TEMPLATE(c) \
286 do { \
287 c##TYPE *px0 = c##PTR(x0), *px1 = c##PTR(x1); \
288 if (op_ct == 'C') \
289 for (j = 0, k = 0, k_ = 0; j < n; ++j) { \
290 kend = pp0 [j]; \
291 kend_ = pp0_[j]; \
292 while (k < kend) { \
293 if (pi0[k] > j) \
294 k = kend; \
295 else { \
296 while (k_ < kend_ && pi0_[k_] < pi0[k]) { \
297 l = iwork[pi0_[k_]]++; \
298 c##ASSIGN_CONJ(px1[pp1[ j]], px0[l]); \
299 c##ASSIGN_IDEN(px1[pp1[pi0_[k_]]], px0[l]); \
300 c##MULTIPLY (px1[pp1[ j]], -0.5); \
301 c##MULTIPLY (px1[pp1[pi0_[k_]]], 0.5); \
302 pi1[pp1[ j]++] = pi0_[k_]; \
303 pi1[pp1[pi0_[k_]]++] = j; \
304 ++k_; \
305 } \
306 l = iwork[j]++; \
307 if (pi0[k] == j) \
308 ++k_; \
309 else { \
310 c##ASSIGN_IDEN(px1[pp1[ j]], px0[k]); \
311 c##ASSIGN_CONJ(px1[pp1[pi0[k]]], px0[k]); \
312 if (k_ < kend_ && pi0_[k_] == pi0[k]) { \
313 l = iwork[pi0[k]]++; \
314 c##DECREMENT_CONJ(px1[pp1[ j]], px0[l]); \
315 c##DECREMENT_IDEN(px1[pp1[pi0[k]]], px0[l]); \
316 ++k_; \
317 } \
318 c##MULTIPLY (px1[pp1[ j]], 0.5); \
319 c##MULTIPLY (px1[pp1[pi0[k]]], -0.5); \
320 pi1[pp1[ j]++] = pi0[k]; \
321 pi1[pp1[pi0[k]]++] = j; \
322 } \
323 ++k; \
324 } \
325 } \
326 while (k_ < kend_) { \
327 if (pi0_[k_] >= j) \
328 k_ = kend_; \
329 else { \
330 l = iwork[pi0_[k_]]++; \
331 c##ASSIGN_CONJ(px1[pp1[ j]], px0[l]); \
332 c##ASSIGN_IDEN(px1[pp1[pi0_[k_]]], px0[l]); \
333 c##MULTIPLY (px1[pp1[ j]], -0.5); \
334 c##MULTIPLY (px1[pp1[pi0_[k_]]], 0.5); \
335 pi1[pp1[ j]++] = pi0_[k_]; \
336 pi1[pp1[pi0_[k_]]++] = j; \
337 ++k_; \
338 } \
339 } \
340 } \
341 else \
342 for (j = 0, k = 0, k_ = 0; j < n; ++j) { \
343 kend = pp0 [j]; \
344 kend_ = pp0_[j]; \
345 while (k < kend) { \
346 if (pi0[k] > j) \
347 k = kend; \
348 else { \
349 while (k_ < kend_ && pi0_[k_] < pi0[k]) { \
350 l = iwork[pi0_[k_]]++; \
351 c##ASSIGN_IDEN(px1[pp1[ j]], px0[l]); \
352 c##ASSIGN_IDEN(px1[pp1[pi0_[k_]]], px0[l]); \
353 c##MULTIPLY (px1[pp1[ j]], -0.5); \
354 c##MULTIPLY (px1[pp1[pi0_[k_]]], 0.5); \
355 pi1[pp1[ j]++] = pi0_[k_]; \
356 pi1[pp1[pi0_[k_]]++] = j; \
357 ++k_; \
358 } \
359 l = iwork[j]++; \
360 if (pi0[k] == j) \
361 ++k_; \
362 else { \
363 c##ASSIGN_IDEN(px1[pp1[ j]], px0[k]); \
364 c##ASSIGN_IDEN(px1[pp1[pi0[k]]], px0[k]); \
365 if (k_ < kend_ && pi0_[k_] == pi0[k]) { \
366 l = iwork[pi0[k]]++; \
367 c##DECREMENT_IDEN(px1[pp1[ j]], px0[l]); \
368 c##DECREMENT_IDEN(px1[pp1[pi0[k]]], px0[l]); \
369 ++k_; \
370 } \
371 c##MULTIPLY (px1[pp1[ j]], 0.5); \
372 c##MULTIPLY (px1[pp1[pi0[k]]], -0.5); \
373 pi1[pp1[ j]++] = pi0[k]; \
374 pi1[pp1[pi0[k]]++] = j; \
375 } \
376 ++k; \
377 } \
378 } \
379 while (k_ < kend_) { \
380 if (pi0_[k_] >= j) \
381 k_ = kend_; \
382 else { \
383 l = iwork[pi0_[k_]]++; \
384 c##ASSIGN_IDEN(px1[pp1[ j]], px0[l]); \
385 c##ASSIGN_IDEN(px1[pp1[pi0_[k_]]], px0[l]); \
386 c##MULTIPLY (px1[pp1[ j]], -0.5); \
387 c##MULTIPLY (px1[pp1[pi0_[k_]]], 0.5); \
388 pi1[pp1[ j]++] = pi0_[k_]; \
389 pi1[pp1[pi0_[k_]]++] = j; \
390 ++k_; \
391 } \
392 } \
393 } \
394 } while (0)
395
396 SWITCH2(cl[0], TEMPLATE);
397
398#undef TEMPLATE
399
400 Matrix_Free(iwork, liwork);
401 UNPROTECT(3); /* x1, i1, p1 */
402
403 }
404
405 UNPROTECT(3); /* x0, i0, p0 */
406
407 } else {
408
409 SEXP i0 = PROTECT(GET_SLOT(from, Matrix_iSym)),
410 j0 = PROTECT(GET_SLOT(from, Matrix_jSym)),
411 x0 = PROTECT(GET_SLOT(from, Matrix_xSym));
412 int *pi0 = INTEGER(i0), *pj0 = INTEGER(j0);
413 R_xlen_t k, nnz0 = XLENGTH(i0), nnz1;
414
415 if (class[1] == 's') {
416
417 SET_SLOT(to, Matrix_iSym, i0);
418 SET_SLOT(to, Matrix_jSym, j0);
419
420 /* Skew-symmetric part of Hermitian matrix is imaginary part */
421 /* Skew-Hermitian part of symmetric matrix is imaginary part */
422 SEXP x1 = PROTECT(Rf_allocVector(CPLXSXP, nnz0));
423 zvimag(COMPLEX(x1), COMPLEX(x0), (size_t) nnz0);
424 SET_SLOT(to, Matrix_xSym, x1);
425 UNPROTECT(1); /* x1 */
426
427 } else {
428
429 nnz1 = nnz0;
430 for (k = 0; k < nnz0; ++k)
431 if (pi0[k] == pj0[k])
432 --nnz1;
433 nnz1 *= 2;
434
435 SEXP i1 = PROTECT(Rf_allocVector(INTSXP, nnz1)),
436 j1 = PROTECT(Rf_allocVector(INTSXP, nnz1)),
437 x1 = PROTECT(Rf_allocVector(TYPEOF(x0), nnz1));
438 int *pi1 = INTEGER(i1), *pj1 = INTEGER(j1);
439 SET_SLOT(to, Matrix_iSym, i1);
440 SET_SLOT(to, Matrix_jSym, j1);
441 SET_SLOT(to, Matrix_xSym, x1);
442
443#define TEMPLATE(c) \
444 do { \
445 c##TYPE *px0 = c##PTR(x0), *px1 = c##PTR(x1); \
446 for (k = 0; k < nnz0; ++k) { \
447 if (*pi0 != *pj0) { \
448 *pi1 = *pi0; \
449 *pj1 = *pj0; \
450 c##ASSIGN_IDEN(*px1, *px0); \
451 c##MULTIPLY(*px1, 0.5); \
452 ++pi1; ++pj1; ++px1; \
453 *pi1 = *pj0; \
454 *pj1 = *pi0; \
455 if (op_ct == 'C') \
456 c##ASSIGN_CONJ(*px1, *px0); \
457 else \
458 c##ASSIGN_IDEN(*px1, *px0); \
459 c##MULTIPLY(*px1, -0.5); \
460 ++pi1; ++pj1; ++px1; \
461 } \
462 ++pi0; ++pj0; ++px0; \
463 } \
464 } while (0)
465
466 SWITCH2(cl[0], TEMPLATE);
467
468#undef TEMPLATE
469
470 UNPROTECT(3); /* x1, j1, i1 */
471
472 }
473
474 UNPROTECT(3); /* x0, j0, i0 */
475
476 }
477
478 UNPROTECT(2); /* to, from */
479 return to;
480}
481
482SEXP R_dense_skewpart(SEXP s_from, SEXP s_trans)
483{
484 const char *class = Matrix_class(s_from, valid_dense, 6, __func__);
485
486 char ct;
487 VALID_TRANS(s_trans, ct);
488
489 return dense_skewpart(s_from, class, ct);
490}
491
492SEXP R_sparse_skewpart(SEXP s_from, SEXP s_trans)
493{
494 const char *class = Matrix_class(s_from, valid_sparse, 6, __func__);
495
496 char ct;
497 VALID_TRANS(s_trans, ct);
498
499 return sparse_skewpart(s_from, class, ct);
500}
SWITCH2
#define SWITCH2(c, template)
Definition M5.h:291
SET_DIM
#define SET_DIM(x, m, n)
Definition Mdefines.h:87
valid_dense
const char * valid_dense[]
Definition objects.c:3
Matrix_Calloc
#define Matrix_Calloc(p, n, t)
Definition Mdefines.h:45
_
#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
allocZero
SEXP allocZero(SEXPTYPE, R_xlen_t)
Definition utils.c:69
Matrix_class
const char * Matrix_class(SEXP, const char **, int, const char *)
Definition objects.c:112
Matrix_Free
#define Matrix_Free(p, n)
Definition Mdefines.h:56
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
zvimag
void zvimag(Rcomplex *x, const Rcomplex *y, size_t n)
Definition idz.c:1287
ncsptrans
void ncsptrans(int *, int *, int *, const int *, const int *, const int *, int, int, char, int *)
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
R_sparse_skewpart
SEXP R_sparse_skewpart(SEXP s_from, SEXP s_trans)
Definition skewpart.c:492
TEMPLATE
#define TEMPLATE(c)
sparse_as_kind
SEXP sparse_as_kind(SEXP, const char *, char)
Definition coerce.c:2076
dense_skewpart
SEXP dense_skewpart(SEXP from, const char *class, char op_ct)
Definition skewpart.c:11
sparse_skewpart
SEXP sparse_skewpart(SEXP from, const char *class, char op_ct)
Definition skewpart.c:169
R_dense_skewpart
SEXP R_dense_skewpart(SEXP s_from, SEXP s_trans)
Definition skewpart.c:482
dense_as_kind
SEXP dense_as_kind(SEXP, const char *, char, int)
Definition coerce.c:2000
Generated by doxygen 1.12.0