Matrix r5059
Loading...
Searching...
No Matches
symmpart.c
Go to the documentation of this file.
1/* C implementation of methods for symmpart */
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_symmpart(SEXP from, const char *class,
12 char op_ul, char op_ct)
13{
14 PROTECT(from = dense_as_kind(from, class, ',', 0));
15
16 if (class[1] != 'g')
17 op_ul = '\0';
18 if (class[0] != 'z')
19 op_ct = '\0';
20
21 char ct = '\0';
22 if (class[1] == 's' && class[0] == 'z')
23 ct = TRANS(from);
24
25 if (class[1] == 's' && op_ct == ct) {
26 UNPROTECT(1); /* from */
27 return from;
28 }
29
30 int packed = class[2] == 'p';
31
32 char cl[] = ".s.Matrix";
33 cl[0] = (class[0] == 'z') ? 'z' : 'd';
34 cl[2] = (packed) ? 'p' : 'y';
35 SEXP to = PROTECT(newObject(cl));
36
37 int *pdim = DIM(from), n = pdim[1];
38 if (pdim[0] != n)
39 Rf_error((op_ct == 'C')
40 ? _("attempt to get Hermitian part of non-square matrix")
41 : _("attempt to get symmetric part of non-square matrix"));
42 SET_DIM(to, n, n);
43 SET_DIMNAMES(to, -(class[1] != 's'), DIMNAMES(from, 0));
44
45 char ul = '\0', nu = '\0';
46 if (class[1] != 'g' && (ul = UPLO(from)) != 'U')
47 SET_UPLO(to);
48 if (class[1] == 't' && (nu = DIAG(from)) != 'N')
49 ;
50
51 if (op_ul != '\0' && op_ul != 'U')
52 SET_UPLO(to);
53 if (op_ct != '\0' && op_ct != 'C')
54 SET_TRANS(to);
55
56 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym)),
57 x1 = PROTECT(Rf_allocVector(TYPEOF(x0), XLENGTH(x0)));
58
59 if (class[1] == 's') {
60
61 /* Symmetric part of Hermitian matrix is real part */
62 /* Hermitian part of symmetric matrix is real part */
63 zvreal(COMPLEX(x1), COMPLEX(x0), (size_t) XLENGTH(x0));
64
65 } else {
66
67 int i, j;
68
69#define TEMPLATE(c) \
70 do { \
71 c##TYPE *px0 = c##PTR(x0), *pu0 = px0, *pl0 = px0; \
72 c##TYPE *px1 = c##PTR(x1), *pu1 = px1, *pl1 = px1; \
73 if (!packed) \
74 memset(px1, 0, sizeof(c##TYPE) * (size_t) XLENGTH(x1)); \
75 if (class[1] == 'g') { \
76 if (op_ul == 'U') { \
77 if (op_ct == 'C') \
78 for (j = 0; j < n; ++j) { \
79 for (i = 0; i < j; ++i) { \
80 c##ASSIGN_IDEN(*pu1, *pu0); \
81 c##INCREMENT_CONJ(*pu1, *pl0); \
82 c##MULTIPLY(*pu1, 0.5); \
83 pu0 += 1; \
84 pu1 += 1; \
85 pl0 += n; \
86 } \
87 c##ASSIGN_PROJ_REAL(*pu1, *pu0); \
88 pu0 += 1; \
89 pu1 += 1; \
90 pu0 += n - j - 1; \
91 pu1 += n - j - 1; \
92 pl0 = px0 + j + 1; \
93 } \
94 else \
95 for (j = 0; j < n; ++j) { \
96 for (i = 0; i < j; ++i) { \
97 c##ASSIGN_IDEN(*pu1, *pu0); \
98 c##INCREMENT_IDEN(*pu1, *pl0); \
99 c##MULTIPLY(*pu1, 0.5); \
100 pu0 += 1; \
101 pu1 += 1; \
102 pl0 += n; \
103 } \
104 c##ASSIGN_IDEN(*pu1, *pu0); \
105 pu0 += 1; \
106 pu1 += 1; \
107 pu0 += n - j - 1; \
108 pu1 += n - j - 1; \
109 pl0 = px0 + j + 1; \
110 } \
111 } else { \
112 if (op_ct == 'C') \
113 for (j = 0; j < n; ++j) { \
114 pl0 += j; \
115 pl1 += j; \
116 pu0 = pl0 + n; \
117 c##ASSIGN_PROJ_REAL(*pl1, *pl0); \
118 pl0 += 1; \
119 pl1 += 1; \
120 for (i = j + 1; i < n; ++i) { \
121 c##ASSIGN_IDEN(*pl1, *pl0); \
122 c##INCREMENT_CONJ(*pl1, *pu0); \
123 c##MULTIPLY(*pl1, 0.5); \
124 pl0 += 1; \
125 pl1 += 1; \
126 pu0 += n; \
127 } \
128 } \
129 else \
130 for (j = 0; j < n; ++j) { \
131 pl0 += j; \
132 pl1 += j; \
133 pu0 = pl0 + n; \
134 c##ASSIGN_IDEN(*pl1, *pl0); \
135 pl0 += 1; \
136 pl1 += 1; \
137 for (i = j + 1; i < n; ++i) { \
138 c##ASSIGN_IDEN(*pl1, *pl0); \
139 c##INCREMENT_IDEN(*pl1, *pu0); \
140 c##MULTIPLY(*pl1, 0.5); \
141 pl0 += 1; \
142 pl1 += 1; \
143 pu0 += n; \
144 } \
145 } \
146 } \
147 } else { \
148 if (ul == 'U') \
149 for (j = 0; j < n; ++j) { \
150 for (i = 0; i < j; ++i) { \
151 c##ASSIGN_IDEN(*pu1, *pu0); \
152 c##MULTIPLY(*pu1, 0.5); \
153 pu0 += 1; \
154 pu1 += 1; \
155 } \
156 if (nu != 'N') \
157 c##SET_UNIT(*pu1); \
158 else if (op_ct == 'C') \
159 c##ASSIGN_PROJ_REAL(*pu1, *pu0); \
160 else \
161 c##ASSIGN_IDEN(*pu1, *pu0); \
162 pu0 += 1; \
163 pu1 += 1; \
164 if (!packed) { \
165 pu0 += n - j - 1; \
166 pu1 += n - j - 1; \
167 } \
168 } \
169 else \
170 for (j = 0; j < n; ++j) { \
171 if (!packed) { \
172 pl0 += j; \
173 pl1 += j; \
174 } \
175 if (nu != 'N') \
176 c##SET_UNIT(*pl1); \
177 else if (op_ct == 'C') \
178 c##ASSIGN_PROJ_REAL(*pl1, *pl0); \
179 else \
180 c##ASSIGN_IDEN(*pl1, *pl0); \
181 pl0 += 1; \
182 pl1 += 1; \
183 for (i = j + 1; i < n; ++i) { \
184 c##ASSIGN_IDEN(*pl1, *pl0); \
185 c##MULTIPLY(*pl1, 0.5); \
186 pl0 += 1; \
187 pl1 += 1; \
188 } \
189 } \
190 } \
191 } while (0)
192
193 SWITCH2(class[0], TEMPLATE);
194
195#undef TEMPLATE
196
197 }
198
199 SET_SLOT(to, Matrix_xSym, x1);
200
201 UNPROTECT(4); /* x1, x0, to, from */
202 return to;
203}
204
205SEXP sparse_symmpart(SEXP from, const char *class,
206 char op_ul, char op_ct)
207{
208 PROTECT(from = sparse_as_kind(from, class, ','));
209
210 if (class[1] != 'g')
211 op_ul = '\0';
212 if (class[0] != 'z')
213 op_ct = '\0';
214
215 char ct = '\0';
216 if (class[1] == 's' && class[0] == 'z')
217 ct = TRANS(from);
218
219 if (class[1] == 's' && op_ct == ct) {
220 UNPROTECT(1); /* from */
221 return from;
222 }
223
224 char cl[] = ".s.Matrix";
225 cl[0] = (class[0] == 'z') ? 'z' : 'd';
226 cl[2] = class[2];
227 SEXP to = PROTECT(newObject(cl));
228
229 int *pdim = DIM(from), n = pdim[1];
230 if (pdim[0] != n)
231 Rf_error((op_ct == 'C')
232 ? _("attempt to get Hermitian part of non-square matrix")
233 : _("attempt to get symmetric part of non-square matrix"));
234 SET_DIM(to, n, n);
235 SET_DIMNAMES(to, -(class[1] != 's'), DIMNAMES(from, 0));
236
237 char ul = '\0', nu = '\0';
238 if (class[1] != 'g' && (ul = UPLO(from)) != 'U')
239 SET_UPLO(to);
240 if (class[1] == 't' && (nu = DIAG(from)) != 'N')
241 ;
242
243 if (op_ul != '\0' && op_ul != 'U')
244 SET_UPLO(to);
245 if (op_ct != '\0' && op_ct != 'C')
246 SET_TRANS(to);
247
248 int up = (class[0] != 'R') == (op_ul == 'U' || ul == 'U');
249
250 if (class[2] != 'T') {
251
252 SEXP iSym = (class[2] == 'C') ? Matrix_iSym : Matrix_jSym,
253 p0 = PROTECT(GET_SLOT(from, Matrix_pSym)),
254 i0 = PROTECT(GET_SLOT(from, iSym)),
255 x0 = PROTECT(GET_SLOT(from, Matrix_xSym));
256 int *pp0 = INTEGER(p0), *pi0 = INTEGER(i0),
257 j, k, kend, nnz0 = pp0[n], nnz1;
258 pp0++;
259
260 if (class[1] == 'g') {
261
262 int *iwork = NULL;
263 size_t liwork = (size_t) ((int_fast64_t) n + n + 1 + nnz0);
264 Matrix_Calloc(iwork, liwork, int);
265
266 int *pp0_ = iwork + n + 1, *pi0_ = iwork + n + n + 1;
267 ncsptrans(pp0_ - 1, pi0_, NULL, pp0 - 1, pi0, NULL, n, n, 'T', iwork);
268 memcpy(iwork, pp0 - 1, sizeof(int) * (size_t) n);
269
270 SEXP p1 = PROTECT(Rf_allocVector(INTSXP, XLENGTH(p0)));
271 int *pp1 = INTEGER(p1), k_, kend_;
272 *(pp1++) = 0;
273 SET_SLOT(to, Matrix_pSym, p1);
274
275 for (j = 0, k = 0, k_ = 0; j < n; ++j) {
276 kend = pp0 [j];
277 kend_ = pp0_[j];
278 pp1[j] = pp1[j - 1];
279 while (k < kend) {
280 if (pi0[k] > j)
281 k = kend;
282 else {
283 while (k_ < kend_ && pi0_[k_] < pi0[k]) {
284 ++pp1[j];
285 ++k_;
286 }
287 ++pp1[j];
288 if (k_ < kend_ && pi0_[k_] == pi0[k])
289 ++k_;
290 ++k;
291 }
292 }
293 while (k_ < kend_) {
294 if (pi0_[k_] > j)
295 k_ = kend_;
296 else {
297 ++pp1[j];
298 ++k_;
299 }
300 }
301 }
302 nnz1 = pp1[n - 1];
303
304 SEXP i1 = PROTECT(Rf_allocVector(INTSXP, nnz1)),
305 x1 = PROTECT(Rf_allocVector(TYPEOF(x0), nnz1));
306 int *pi1 = INTEGER(i1), l;
307 SET_SLOT(to, iSym, i1);
308 SET_SLOT(to, Matrix_xSym, x1);
309
310#define TEMPLATE(c) \
311 do { \
312 c##TYPE *px0 = c##PTR(x0), *px1 = c##PTR(x1); \
313 if (op_ct == 'C') \
314 for (j = 0, k = 0, k_ = 0; j < n; ++j) { \
315 kend = pp0 [j]; \
316 kend_ = pp0_[j]; \
317 while (k < kend) { \
318 if (up && pi0[k] > j) \
319 k = kend; \
320 else if (!up && pi0[k] < j) \
321 ++k; \
322 else { \
323 while (k_ < kend_ && pi0_[k_] < pi0[k]) { \
324 l = iwork[pi0_[k_]]++; \
325 *pi1 = pi0_[k_]; \
326 c##ASSIGN_CONJ(*px1, px0[l]); \
327 c##MULTIPLY(*px1, 0.5); \
328 ++k_; ++pi1; ++px1; \
329 } \
330 l = iwork[j]++; \
331 *pi1 = pi0[k]; \
332 if (pi0[k] == j) { \
333 c##ASSIGN_PROJ_REAL(*px1, px0[k]); \
334 ++k_; \
335 } else { \
336 c##ASSIGN_IDEN(*px1, px0[k]); \
337 if (k_ < kend_ && pi0_[k_] == pi0[k]) { \
338 l = iwork[pi0[k]]++; \
339 c##INCREMENT_CONJ(*px1, px0[l]); \
340 ++k_; \
341 } \
342 c##MULTIPLY(*px1, 0.5); \
343 } \
344 ++k; ++pi1; ++px1; \
345 } \
346 } \
347 while (k_ < kend_) { \
348 if (up && pi0_[k_] > j) \
349 k_ = kend_; \
350 else if (!up && pi0_[k_] < j) \
351 ++k_; \
352 else { \
353 l = iwork[pi0_[k_]]++; \
354 *pi1 = pi0_[k_]; \
355 c##ASSIGN_CONJ(*px1, px0[l]); \
356 c##MULTIPLY(*px1, 0.5); \
357 ++k_; ++pi1; ++px1; \
358 } \
359 } \
360 } \
361 else \
362 for (j = 0, k = 0, k_ = 0; j < n; ++j) { \
363 kend = pp0 [j]; \
364 kend_ = pp0_[j]; \
365 while (k < kend) { \
366 if (up && pi0[k] > j) \
367 k = kend; \
368 else if (!up && pi0[k] < j) \
369 ++k; \
370 else { \
371 while (k_ < kend_ && pi0_[k_] < pi0[k]) { \
372 l = iwork[pi0_[k_]]++; \
373 *pi1 = pi0_[k_]; \
374 c##ASSIGN_IDEN(*px1, px0[l]); \
375 c##MULTIPLY(*px1, 0.5); \
376 ++k_; ++pi1; ++px1; \
377 } \
378 l = iwork[j]++; \
379 *pi1 = pi0[k]; \
380 if (pi0[k] == j) { \
381 c##ASSIGN_IDEN(*px1, px0[k]); \
382 ++k_; \
383 } else { \
384 c##ASSIGN_IDEN(*px1, px0[k]); \
385 if (k_ < kend_ && pi0_[k_] == pi0[k]) { \
386 l = iwork[pi0[k]]++; \
387 c##INCREMENT_IDEN(*px1, px0[l]); \
388 ++k_; \
389 } \
390 c##MULTIPLY(*px1, 0.5); \
391 } \
392 ++k; ++pi1; ++px1; \
393 } \
394 } \
395 while (k_ < kend_) { \
396 if (up && pi0_[k_] > j) \
397 k_ = kend_; \
398 else if (!up && pi0_[k_] < j) \
399 ++k_; \
400 else { \
401 l = iwork[pi0_[k_]]++; \
402 *pi1 = pi0_[k_]; \
403 c##ASSIGN_IDEN(*px1, px0[l]); \
404 c##MULTIPLY(*px1, 0.5); \
405 ++k_; ++pi1; ++px1; \
406 } \
407 } \
408 } \
409 } while (0)
410
411 SWITCH2(cl[0], TEMPLATE);
412
413#undef TEMPLATE
414
415 Matrix_Free(iwork, liwork);
416 UNPROTECT(3); /* x1, i1, p1 */
417
418 } else if (class[1] == 's') {
419
420 SET_SLOT(to, Matrix_pSym, p0);
421 SET_SLOT(to, iSym, i0);
422
423 /* Symmetric part of Hermitian matrix is real part */
424 /* Hermitian part of symmetric matrix is real part */
425 SEXP x1 = PROTECT(Rf_allocVector(CPLXSXP, nnz0));
426 zvreal(COMPLEX(x1), COMPLEX(x0), (size_t) nnz0);
427 SET_SLOT(to, Matrix_xSym, x1);
428 UNPROTECT(1); /* x1 */
429
430 } else if (nu == 'N') {
431
432 SET_SLOT(to, Matrix_pSym, p0);
433 SET_SLOT(to, iSym, i0);
434
435 SEXP x1 = PROTECT(Rf_allocVector(TYPEOF(x0), nnz0));
436 SET_SLOT(to, Matrix_xSym, x1);
437
438#define TEMPLATE(c) \
439 do { \
440 c##TYPE *px0 = c##PTR(x0), *px1 = c##PTR(x1); \
441 for (j = 0, k = 0; j < n; ++j) { \
442 kend = pp0[j]; \
443 while (k < kend) { \
444 if (pi0[k] != j) { \
445 c##ASSIGN_IDEN(*px1, px0[k]); \
446 c##MULTIPLY(*px1, 0.5); \
447 } \
448 else if (op_ct == 'C') \
449 c##ASSIGN_PROJ_REAL(*px1, px0[k]); \
450 else \
451 c##ASSIGN_IDEN(*px1, px0[k]); \
452 ++k; ++px1; \
453 } \
454 } \
455 } while (0)
456
457 SWITCH2(cl[0], TEMPLATE);
458
459#undef TEMPLATE
460
461 UNPROTECT(1); /* x1 */
462
463 } else {
464
465 nnz1 = nnz0 + n;
466
467 SEXP p1 = PROTECT(Rf_allocVector(INTSXP, XLENGTH(p0))),
468 i1 = PROTECT(Rf_allocVector(INTSXP, nnz1)),
469 x1 = PROTECT(Rf_allocVector(TYPEOF(x0), nnz1));
470 int *pp1 = INTEGER(p1), *pi1 = INTEGER(i1);
471 *(pp1++) = 0;
472 SET_SLOT(to, Matrix_pSym, p1);
473 SET_SLOT(to, iSym, i1);
474 SET_SLOT(to, Matrix_xSym, x1);
475
476#define TEMPLATE(c) \
477 do { \
478 c##TYPE *px0 = c##PTR(x0), *px1 = c##PTR(x1); \
479 for (j = 0, k = 0; j < n; ++j) { \
480 kend = pp0[j]; \
481 if (!up) { \
482 *pi1 = j; \
483 *px1 = c##UNIT; \
484 ++pi1; ++px1; \
485 } \
486 while (k < kend) { \
487 *pi1 = pi0[k]; \
488 c##ASSIGN_IDEN(*px1, px0[k]); \
489 c##MULTIPLY(*px1, 0.5); \
490 ++k; ++pi1; ++px1; \
491 } \
492 if (up) { \
493 *pi1 = j; \
494 *px1 = c##UNIT; \
495 ++pi1; ++px1; \
496 } \
497 pp1[j] = kend + j + 1; \
498 } \
499 } while (0)
500
501 SWITCH2(cl[0], TEMPLATE);
502
503#undef TEMPLATE
504
505 UNPROTECT(3); /* x1, i1, p1 */
506
507 }
508
509 UNPROTECT(3); /* x0, i0, p0 */
510
511 } else {
512
513 SEXP i0 = PROTECT(GET_SLOT(from, Matrix_iSym)),
514 j0 = PROTECT(GET_SLOT(from, Matrix_jSym)),
515 x0 = PROTECT(GET_SLOT(from, Matrix_xSym));
516 int *pi0 = INTEGER(i0), *pj0 = INTEGER(j0);
517 R_xlen_t k, nnz0 = XLENGTH(i0), nnz1;
518
519 if (class[1] == 'g') {
520
521 SEXP i1 = PROTECT(Rf_allocVector(INTSXP, nnz0)),
522 j1 = PROTECT(Rf_allocVector(INTSXP, nnz0)),
523 x1 = PROTECT(Rf_allocVector(TYPEOF(x0), nnz0));
524 int *pi1 = INTEGER(i1), *pj1 = INTEGER(j1);
525 SET_SLOT(to, Matrix_iSym, i1);
526 SET_SLOT(to, Matrix_jSym, j1);
527 SET_SLOT(to, Matrix_xSym, x1);
528
529#define TEMPLATE(c) \
530 do { \
531 c##TYPE *px0 = c##PTR(x0), *px1 = c##PTR(x1); \
532 if (op_ct == 'C') \
533 for (k = 0; k < nnz0; ++k) { \
534 if (*pi0 != *pj0) { \
535 if ((up) ? *pi0 < *pj0 : *pi0 > *pj0) { \
536 *pi1 = *pi0; \
537 *pj1 = *pj0; \
538 c##ASSIGN_IDEN(*px1, *px0); \
539 } else { \
540 *pi1 = *pj0; \
541 *pj1 = *pi0; \
542 c##ASSIGN_CONJ(*px1, *px0); \
543 } \
544 c##MULTIPLY(*px1, 0.5); \
545 } else { \
546 *pi1 = *pi0; \
547 *pj1 = *pj0; \
548 c##ASSIGN_PROJ_REAL(*px1, *px0); \
549 } \
550 ++pi0; ++pi1; ++pj0; ++pj1; ++px0; ++px1; \
551 } \
552 else \
553 for (k = 0; k < nnz0; ++k) { \
554 if ((up) ? *pi0 <= *pj0 : *pi0 >= *pj0) { \
555 *pi1 = *pi0; \
556 *pj1 = *pj0; \
557 } else { \
558 *pi1 = *pj0; \
559 *pj1 = *pi0; \
560 } \
561 c##ASSIGN_IDEN(*px1, *px0); \
562 if (*pi0 != *pj0) \
563 c##MULTIPLY(*px1, 0.5); \
564 ++pi0; ++pi1; ++pj0; ++pj1; ++px0; ++px1; \
565 } \
566 } while (0)
567
568 SWITCH2(cl[0], TEMPLATE);
569
570#undef TEMPLATE
571
572 UNPROTECT(3); /* x1, j1, i1 */
573
574 } else if (class[1] == 's') {
575
576 SET_SLOT(to, Matrix_iSym, i0);
577 SET_SLOT(to, Matrix_jSym, j0);
578
579 /* Symmetric part of Hermitian matrix is real part */
580 /* Hermitian part of symmetric matrix is real part */
581 SEXP x1 = PROTECT(Rf_allocVector(CPLXSXP, nnz0));
582 zvreal(COMPLEX(x1), COMPLEX(x0), (size_t) nnz0);
583 SET_SLOT(to, Matrix_xSym, x1);
584 UNPROTECT(1); /* x1 */
585
586 } else if (nu == 'N') {
587
588 SET_SLOT(to, Matrix_iSym, i0);
589 SET_SLOT(to, Matrix_jSym, j0);
590
591 SEXP x1 = PROTECT(Rf_allocVector(TYPEOF(x0), nnz0));
592 SET_SLOT(to, Matrix_xSym, x1);
593
594#define TEMPLATE(c) \
595 do { \
596 c##TYPE *px0 = c##PTR(x0), *px1 = c##PTR(x1); \
597 for (k = 0; k < nnz0; ++k) { \
598 if (*pi0 != *pj0) { \
599 c##ASSIGN_IDEN(*px1, *px0); \
600 c##MULTIPLY(*px1, 0.5); \
601 } \
602 else if (op_ct == 'C') \
603 c##ASSIGN_PROJ_REAL(*px1, *px0); \
604 else \
605 c##ASSIGN_IDEN(*px1, *px0); \
606 ++pi0; ++pj0; ++px0; ++px1; \
607 } \
608 } while (0)
609
610 SWITCH2(cl[0], TEMPLATE);
611
612#undef TEMPLATE
613
614 UNPROTECT(1); /* x1 */
615
616 } else {
617
618 nnz1 = nnz0 + n;
619
620 SEXP i1 = PROTECT(Rf_allocVector(INTSXP, nnz1)),
621 j1 = PROTECT(Rf_allocVector(INTSXP, nnz1)),
622 x1 = PROTECT(Rf_allocVector(TYPEOF(x0), nnz1));
623 int *pi1 = INTEGER(i1), *pj1 = INTEGER(j1), j;
624 SET_SLOT(to, Matrix_iSym, i1);
625 SET_SLOT(to, Matrix_jSym, j1);
626 SET_SLOT(to, Matrix_xSym, x1);
627
628#define TEMPLATE(c) \
629 do { \
630 c##TYPE *px0 = c##PTR(x0), *px1 = c##PTR(x1); \
631 for (k = 0; k < nnz0; ++k) { \
632 *pi1 = *pi0; \
633 *pj1 = *pj0; \
634 c##ASSIGN_IDEN(*px1, *px0); \
635 c##MULTIPLY(*px1, 0.5); \
636 ++pi0; ++pi1; ++pj0; ++pj1; ++px0; ++px1; \
637 } \
638 for (j = 0; j < n; ++j) { \
639 *pi1 = *pj1 = j; \
640 *px1 = c##UNIT; \
641 ++pi1; ++pj1; ++px1; \
642 } \
643 } while (0)
644
645 SWITCH2(cl[0], TEMPLATE);
646
647#undef TEMPLATE
648
649 UNPROTECT(3); /* x1, j1, i1 */
650
651 }
652
653 UNPROTECT(3); /* x0, j0, i1 */
654
655 }
656
657 UNPROTECT(2); /* to, from */
658 return to;
659}
660
661SEXP R_dense_symmpart(SEXP s_from, SEXP s_uplo, SEXP s_trans)
662{
663 const char *class = Matrix_class(s_from, valid_dense, 6, __func__);
664
665 char ul, ct;
666 VALID_UPLO (s_uplo , ul);
667 VALID_TRANS(s_trans, ct);
668
669 return dense_symmpart(s_from, class, ul, ct);
670}
671
672SEXP R_sparse_symmpart(SEXP s_from, SEXP s_uplo, SEXP s_trans)
673{
674 const char *class = Matrix_class(s_from, valid_sparse, 6, __func__);
675
676 char ul, ct;
677 VALID_UPLO (s_uplo , ul);
678 VALID_TRANS(s_trans, ct);
679
680 return sparse_symmpart(s_from, class, ul, ct);
681}
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
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
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
zvreal
void zvreal(Rcomplex *x, const Rcomplex *y, size_t n)
Definition idz.c:1274
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_symmpart
SEXP R_sparse_symmpart(SEXP s_from, SEXP s_uplo, SEXP s_trans)
Definition symmpart.c:672
sparse_symmpart
SEXP sparse_symmpart(SEXP from, const char *class, char op_ul, char op_ct)
Definition symmpart.c:205
TEMPLATE
#define TEMPLATE(c)
R_dense_symmpart
SEXP R_dense_symmpart(SEXP s_from, SEXP s_uplo, SEXP s_trans)
Definition symmpart.c:661
sparse_as_kind
SEXP sparse_as_kind(SEXP, const char *, char)
Definition coerce.c:2076
dense_symmpart
SEXP dense_symmpart(SEXP from, const char *class, char op_ul, char op_ct)
Definition symmpart.c:11
dense_as_kind
SEXP dense_as_kind(SEXP, const char *, char, int)
Definition coerce.c:2000
Generated by doxygen 1.12.0