Matrix r5059
Loading...
Searching...
No Matches
diag.c
Go to the documentation of this file.
1/* C implementation of methods for diag, diag<- */
2
3#include <math.h> /* sqrt */
4#include "cholmod-etc.h"
5#include "Mdefines.h"
6#include "M5.h"
7#include "idz.h"
8
9SEXP dense_diag_get(SEXP obj, const char *class, int names)
10{
11 int *pdim = DIM(obj), m = pdim[0], n = pdim[1],
12 r = (m < n) ? m : n, packed = class[2] == 'p';
13
14 char ul = '\0', ct = '\0', nu = '\0';
15 if (class[1] != 'g')
16 ul = UPLO(obj);
17 if (class[1] == 's' && class[0] == 'z')
18 ct = TRANS(obj);
19 if (class[1] == 't')
20 nu = DIAG(obj);
21
22 SEXP ans;
23
24 if (nu != '\0' && nu != 'N') {
25
26 PROTECT(ans = allocUnit(kindToType(class[0]), r));
27
28 } else {
29
30 PROTECT(ans = Rf_allocVector(kindToType(class[0]), r));
31
32 SEXP x = GET_SLOT(obj, Matrix_xSym);
33 size_t m_ = (size_t) m, n_ = (size_t) n, r_ = (size_t) r;
34
35#define TEMPLATE(c) \
36 do { \
37 c##TYPE *pa = c##PTR(ans), *px = c##PTR(x); \
38 if (!packed) \
39 c##NAME(copy2)(r_, pa, 1, px, m_ + 1); \
40 else if (ul == 'U') \
41 c##NAME(copy1)(n_, pa, 1, 0, 0, px, 2 , 1, 0); \
42 else \
43 c##NAME(copy1)(n_, pa, 1, 0, 0, px, n_, 1, 1); \
44 } while (0)
45
46 SWITCH4(class[0], TEMPLATE);
47
48#undef TEMPLATE
49
50 if (class[0] == 'n')
51 naToUnit(ans);
52 if (class[0] == 'z' && class[1] == 's' && ct == 'C')
53 zvreal(COMPLEX(ans), NULL, r_);
54
55 }
56
57 if (names) {
58 /* NB: Logic here must be adjusted once the validity method */
59 /* for 'symmetricMatrix' enforces symmetric 'Dimnames' */
60 SEXP dn = PROTECT(DIMNAMES(obj, 0)),
61 rn = VECTOR_ELT(dn, 0),
62 cn = VECTOR_ELT(dn, 1);
63 if (cn == R_NilValue) {
64 if (class[1] == 's' && rn != R_NilValue)
65 Rf_setAttrib(ans, R_NamesSymbol, rn);
66 } else {
67 if (class[1] == 's')
68 Rf_setAttrib(ans, R_NamesSymbol, cn);
69 else if (rn != R_NilValue &&
70 (rn == cn || equalString(rn, cn, r)))
71 Rf_setAttrib(ans, R_NamesSymbol, (r == m) ? rn : cn);
72 }
73 UNPROTECT(1); /* dn */
74 }
75
76 UNPROTECT(1); /* ans */
77 return ans;
78}
79
80SEXP sparse_diag_get(SEXP obj, const char *class, int names)
81{
82 int *pdim = DIM(obj), m = pdim[0], n = pdim[1],
83 r = (m < n) ? m : n;
84
85 char ul = '\0', ct = '\0', nu = '\0';
86 if (class[1] != 'g')
87 ul = UPLO(obj);
88 if (class[1] == 's' && class[0] == 'z')
89 ct = TRANS(obj);
90 if (class[1] == 't')
91 nu = DIAG(obj);
92
93 SEXP ans;
94
95 if (nu != '\0' && nu != 'N') {
96
97 PROTECT(ans = allocUnit(kindToType(class[0]), r));
98
99 } else if (class[2] != 'T') {
100
101 PROTECT(ans = Rf_allocVector(kindToType(class[0]), r));
102
103 SEXP iSym = (class[2] == 'C') ? Matrix_iSym : Matrix_jSym,
104 p = PROTECT(GET_SLOT(obj, Matrix_pSym)),
105 i = PROTECT(GET_SLOT(obj, iSym));
106 int *pp = INTEGER(p), *pi = INTEGER(i), j, k, kend,
107 up = (class[2] == 'C') == (ul == 'U');
108 pp++;
109
110#define TEMPLATE(c) \
111 do { \
112 c##TYPE *pa = c##PTR(ans); \
113 c##IF_NPATTERN( \
114 SEXP x = GET_SLOT(obj, Matrix_xSym); \
115 c##TYPE *px = c##PTR(x); \
116 ); \
117 if (class[1] == 'g') \
118 for (j = 0, k = 0; j < r; ++j) { \
119 pa[j] = c##ZERO; \
120 kend = pp[j]; \
121 while (k < kend) { \
122 if (pi[k] != j) \
123 ++k; \
124 else { \
125 pa[j] = c##IFELSE_NPATTERN(px[k], c##UNIT); \
126 k = kend; \
127 } \
128 } \
129 } \
130 else \
131 for (j = 0, k = 0; j < r; ++j) { \
132 kend = pp[j]; \
133 pa[j] = (k < kend && pi[(up) ? kend - 1 : k] == j) \
134 ? c##IFELSE_NPATTERN(px[(up) ? kend - 1 : k], c##UNIT) \
135 : c##ZERO; \
136 k = kend; \
137 } \
138 } while (0)
139
140 SWITCH5(class[0], TEMPLATE);
141
142#undef TEMPLATE
143
144 UNPROTECT(2); /* i, p */
145
146 if (class[0] == 'z' && class[1] == 's' && ct == 'C')
147 zvreal(COMPLEX(ans), NULL, (size_t) r);
148
149 } else {
150
151 PROTECT(ans = Rf_allocVector(kindToType(class[0]), r));
152
153 SEXP i = PROTECT(GET_SLOT(obj, Matrix_iSym)),
154 j = PROTECT(GET_SLOT(obj, Matrix_jSym));
155 int *pi = INTEGER(i), *pj = INTEGER(j);
156 R_xlen_t k, kend = XLENGTH(i);
157
158#define TEMPLATE(c) \
159 do { \
160 c##TYPE *pa = c##PTR(ans); \
161 memset(pa, 0, sizeof(c##TYPE) * (size_t) r); \
162 c##IF_NPATTERN( \
163 SEXP x = GET_SLOT(obj, Matrix_xSym); \
164 c##TYPE *px = c##PTR(x); \
165 ); \
166 for (k = 0; k < kend; ++k) \
167 if (pi[k] == pj[k]) \
168 c##INCREMENT_IDEN(pa[pi[k]], px[k]); \
169 } while (0)
170
171 SWITCH5(class[0], TEMPLATE);
172
173#undef TEMPLATE
174
175 UNPROTECT(2); /* j0, i0 */
176
177 if (class[0] == 'z' && class[1] == 's' && ct == 'C')
178 zvreal(COMPLEX(ans), NULL, (size_t) r);
179
180 }
181
182 if (names) {
183 /* NB: Logic here must be adjusted once the validity method */
184 /* for 'symmetricMatrix' enforces symmetric 'Dimnames' */
185 SEXP dn = PROTECT(DIMNAMES(obj, 0)),
186 rn = VECTOR_ELT(dn, 0),
187 cn = VECTOR_ELT(dn, 1);
188 if (cn == R_NilValue) {
189 if (class[1] == 's' && rn != R_NilValue)
190 Rf_setAttrib(ans, R_NamesSymbol, rn);
191 } else {
192 if (class[1] == 's')
193 Rf_setAttrib(ans, R_NamesSymbol, cn);
194 else if (rn != R_NilValue &&
195 (rn == cn || equalString(rn, cn, r)))
196 Rf_setAttrib(ans, R_NamesSymbol, (r == m) ? rn : cn);
197 }
198 UNPROTECT(1); /* dn */
199 }
200
201 UNPROTECT(1); /* ans */
202 return ans;
203}
204
205SEXP R_dense_diag_get(SEXP s_obj, SEXP s_names)
206{
207 const char *class = Matrix_class(s_obj, valid_dense, 6, __func__);
208
209 int names;
210 VALID_LOGIC2(s_names, names);
211
212 return dense_diag_get(s_obj, class, names);
213}
214
215SEXP R_sparse_diag_get(SEXP s_obj, SEXP s_names)
216{
217 const char *class = Matrix_class(s_obj, valid_sparse, 6, __func__);
218
219 int names;
220 VALID_LOGIC2(s_names, names);
221
222 return sparse_diag_get(s_obj, class, names);
223}
224
225SEXP dense_diag_set(SEXP from, const char *class, SEXP value, int new)
226{
227 SEXP to = PROTECT(newObject(class));
228
229 int *pdim = DIM(from), m = pdim[0], n = pdim[1],
230 r = (m < n) ? m : n, packed = class[2] == 'p';
231 SET_DIM(to, m, n);
232 SET_DIMNAMES(to, 0, DIMNAMES(from, 0));
233
234 char ul = '\0', ct = '\0';
235 if (class[1] != 'g' && (ul = UPLO(from)) != 'U')
236 SET_UPLO(to);
237 if (class[1] == 's' && class[0] == 'z' && (ct = TRANS(from)) != 'C')
238 SET_TRANS(to);
239
240 SEXP x = PROTECT(GET_SLOT(from, Matrix_xSym));
241 if (new) {
242 x = duplicateVector(x);
243 UNPROTECT(1); /* x */
244 PROTECT(x);
245 }
246
247 size_t m_ = (size_t) m, n_ = (size_t) n, r_ = (size_t) r,
248 d_ = (LENGTH(value) == r) ? 1 : 0;
249
250#define TEMPLATE(c) \
251 do { \
252 c##TYPE *px = c##PTR(x), *pv = c##PTR(value); \
253 if (!packed) \
254 c##NAME(copy2)(r_, px, m_ + 1, pv, d_); \
255 else if (ul == 'U') \
256 c##NAME(copy1)(n_, px, 2 , 1, 0, pv, d_, 0, 0); \
257 else \
258 c##NAME(copy1)(n_, px, n_, 1, 1, pv, d_, 0, 0); \
259 } while (0)
260
261 SWITCH4(class[0], TEMPLATE);
262
263#undef TEMPLATE
264
265 SET_SLOT(to, Matrix_xSym, x);
266
267 UNPROTECT(2); /* x, to */
268 return to;
269}
270
271SEXP sparse_diag_set(SEXP from, const char *class, SEXP value)
272{
273 SEXP to = PROTECT(newObject(class));
274
275 int *pdim = DIM(from), m = pdim[0], n = pdim[1],
276 r = (m < n) ? m : n;
277 SET_DIM(to, m, n);
278 SET_DIMNAMES(to, 0, DIMNAMES(from, 0));
279
280 char ul = '\0', ct = '\0', nu = '\0';
281 if (class[1] != 'g' && (ul = UPLO(from)) != 'U')
282 SET_UPLO(to);
283 if (class[1] == 's' && class[0] == 'z' && (ct = TRANS(from)) != 'C')
284 SET_TRANS(to);
285 if (class[1] == 't' && (nu = DIAG(from)) != 'N')
286 ;
287
288 int mode = 0;
289
290 if (class[2] != 'T') {
291
292 if (class[2] == 'R')
293 SWAP(m, n, int, );
294
295 SEXP iSym = (class[2] == 'C') ? Matrix_iSym : Matrix_jSym,
296 p0 = PROTECT(GET_SLOT(from, Matrix_pSym)),
297 i0 = PROTECT(GET_SLOT(from, iSym)),
298 p1 = PROTECT(Rf_allocVector(INTSXP, XLENGTH(p0)));
299 int *pp0 = INTEGER(p0), *pp1 = INTEGER(p1), *pi0 = INTEGER(i0),
300 j, k, kend, nd0 = 0, nd1 = 0,
301 up = (class[2] == 'C') == (ul == 'U');
302 pp0++; *(pp1++) = 0;
303 SET_SLOT(to, Matrix_pSym, p1);
304
305 if (class[1] == 'g') {
306 for (j = 0, k = 0; j < r; ++j) {
307 kend = pp0[j];
308 while (k < kend) {
309 if (pi0[k] < j)
310 ++k;
311 else {
312 if (pi0[k] == j)
313 ++nd0;
314 k = kend;
315 }
316 }
317 pp1[j] = kend - nd0;
318 }
319 for (j = r; j < n; ++j)
320 pp1[j] = pp0[j] - nd0;
321 }
322 else if (class[1] == 's' || nu == 'N')
323 for (j = 0, k = 0; j < n; ++j) {
324 kend = pp0[j];
325 if (k < kend && pi0[(up) ? kend - 1 : k] == j)
326 ++nd0;
327 k = kend;
328 pp1[j] = kend - nd0;
329 }
330 else
331 for (j = 0; j < n; ++j)
332 pp1[j] = pp0[j];
333
334#define TEMPLATE(c) \
335 do { \
336 c##TYPE *pv = c##PTR(value); \
337 if (LENGTH(value) == r) { \
338 mode = 1; \
339 for (j = 0; j < r; ++j) { \
340 if (c##NOT_ZERO(pv[j])) \
341 ++nd1; \
342 pp1[j] += nd1; \
343 } \
344 for (j = r; j < n; ++j) \
345 pp1[j] += nd1; \
346 } else if (c##NOT_ZERO(pv[0])) { \
347 mode = 2; \
348 for (j = 0; j < r; ++j) \
349 pp1[j] += ++nd1; \
350 for (j = r; j < n; ++j) \
351 pp1[j] += nd1; \
352 } \
353 } while (0)
354
355 SWITCH4(class[0], TEMPLATE);
356
357#undef TEMPLATE
358
359 if (nd1 - nd0 > INT_MAX - pp0[n - 1])
360 Rf_error(_("%s cannot exceed %s"), "p[length(p)]", "2^31-1");
361
362 SEXP i1 = PROTECT(Rf_allocVector(INTSXP, pp1[n - 1]));
363 int *pi1 = INTEGER(i1);
364 SET_SLOT(to, iSym, i1);
365
366#define TEMPLATE(c) \
367 do { \
368 c##IF_NPATTERN( \
369 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym)), \
370 x1 = PROTECT(Rf_allocVector(c##TYPESXP, pp1[n - 1])); \
371 c##TYPE *px0 = c##PTR(x0), *px1 = c##PTR(x1); \
372 SET_SLOT(to, Matrix_xSym, x1); \
373 UNPROTECT(2); /* x1, x0 */ \
374 ); \
375 c##TYPE *pv = c##PTR(value); \
376 for (j = 0, k = 0; j < r; ++j) { \
377 kend = pp0[j]; \
378 while (k < kend && pi0[k] < j) { \
379 *(pi1++) = pi0[k]; \
380 c##IF_NPATTERN( \
381 *(px1++) = px0[k]; \
382 ); \
383 ++k; \
384 } \
385 if (k < kend && pi0[k] == j) \
386 ++k; \
387 if (mode > 1 || (mode > 0 && c##NOT_ZERO(pv[j]))) { \
388 *(pi1++) = j; \
389 c##IF_NPATTERN( \
390 *(px1++) = pv[(mode > 1) ? 0 : j]; \
391 ); \
392 } \
393 while (k < kend) { \
394 *(pi1++) = pi0[k]; \
395 c##IF_NPATTERN( \
396 *(px1++) = px0[k]; \
397 ); \
398 ++k; \
399 } \
400 } \
401 for (j = r; j < n; ++j) { \
402 kend = pp0[j]; \
403 while (k < kend) { \
404 *(pi1++) = pi0[k]; \
405 c##IF_NPATTERN( \
406 *(px1++) = px0[k]; \
407 ); \
408 ++k; \
409 } \
410 } \
411 } while (0)
412
413 SWITCH5(class[0], TEMPLATE);
414
415#undef TEMPLATE
416
417 UNPROTECT(4); /* i1, p1, i0, p0 */
418
419 } else {
420
421 SEXP i0 = PROTECT(GET_SLOT(from, Matrix_iSym)),
422 j0 = PROTECT(GET_SLOT(from, Matrix_jSym));
423 int *pi0 = INTEGER(i0), *pj0 = INTEGER(j0), j;
424 R_xlen_t k, nd0 = 0, nd1 = 0, nnz0 = XLENGTH(i0), nnz1 = nnz0;
425
426 if (nu == '\0' || nu == 'N')
427 for (k = 0; k < nnz0; ++k)
428 if (pi0[k] == pj0[k])
429 ++nd0;
430
431#define TEMPLATE(c) \
432 do { \
433 c##TYPE *pv = c##PTR(value); \
434 if (LENGTH(value) == r) { \
435 mode = 1; \
436 for (j = 0; j < r; ++j) \
437 if (c##NOT_ZERO(pv[j])) \
438 ++nd1; \
439 } else if (c##NOT_ZERO(pv[0])) { \
440 mode = 2; \
441 nd1 = r; \
442 } \
443 } while (0)
444
445 SWITCH4(class[0], TEMPLATE);
446
447#undef TEMPLATE
448
449 if (nd1 - nd0 > R_XLEN_T_MAX - nnz0)
450 Rf_error(_("%s cannot exceed %s"), "length(i)", "R_XLEN_T_MAX");
451 nnz1 += nd1 - nd0;
452
453 SEXP i1 = PROTECT(Rf_allocVector(INTSXP, nnz1)),
454 j1 = PROTECT(Rf_allocVector(INTSXP, nnz1));
455 int *pi1 = INTEGER(i1), *pj1 = INTEGER(j1);
456 SET_SLOT(to, Matrix_iSym, i1);
457 SET_SLOT(to, Matrix_jSym, j1);
458
459#define TEMPLATE(c) \
460 do { \
461 c##IF_NPATTERN( \
462 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym)), \
463 x1 = PROTECT(Rf_allocVector(c##TYPESXP, nnz1)); \
464 c##TYPE *px0 = c##PTR(x0), *px1 = c##PTR(x1); \
465 SET_SLOT(to, Matrix_xSym, x1); \
466 UNPROTECT(2); /* x1, x0 */ \
467 ); \
468 c##TYPE *pv = c##PTR(value); \
469 for (k = 0; k < nnz0; ++k) { \
470 if (pi0[k] != pj0[k]) { \
471 *(pi1++) = pi0[k]; \
472 *(pj1++) = pj0[k]; \
473 c##IF_NPATTERN( \
474 *(px1++) = px0[k]; \
475 ); \
476 } \
477 } \
478 if (mode > 1) \
479 for (j = 0; j < r; ++j) { \
480 *(pi1++) = *(pj1++) = j; \
481 c##IF_NPATTERN( \
482 *(px1++) = pv[0]; \
483 ); \
484 } \
485 else if (mode > 0) \
486 for (j = 0; j < r; ++j) { \
487 if (c##NOT_ZERO(pv[j])) { \
488 *(pi1++) = *(pj1++) = j; \
489 c##IF_NPATTERN( \
490 *(px1++) = pv[j]; \
491 ); \
492 } \
493 } \
494 } while (0)
495
496 SWITCH5(class[0], TEMPLATE);
497
498#undef TEMPLATE
499
500 UNPROTECT(4); /* j1, i1, j0, i0 */
501
502 }
503
504 UNPROTECT(1); /* to */
505 return to;
506}
507
508SEXP R_dense_diag_set(SEXP s_from, SEXP s_value)
509{
510 const char *class = Matrix_class(s_from, valid_dense, 6, __func__);
511 SEXPTYPE tx = kindToType(class[0]), tv = TYPEOF(s_value);
512
513 switch (tv) {
514 case LGLSXP:
515 case INTSXP:
516 case REALSXP:
517 case CPLXSXP:
518 break;
519 default:
520 Rf_error(_("replacement diagonal has incompatible type \"%s\""),
521 Rf_type2char(tv));
522 break;
523 }
524
525 int *pdim = DIM(s_from), m = pdim[0], n = pdim[1],
526 r = (m < n) ? m : n;
527 if (XLENGTH(s_value) != 1 && XLENGTH(s_value) != r)
528 Rf_error(_("replacement diagonal has wrong length"));
529
530 int new = 1;
531 if (tv <= tx) {
532 /* defined in ./coerce.c : */
533 SEXP dense_as_general(SEXP, const char *, int);
534 if (class[1] == 's' && class[0] == 'z' && tv == tx &&
535 TRANS(s_from) == 'C') {
536 PROTECT(s_from = dense_as_general(s_from, class, 1));
537 class = Matrix_class(s_from, valid_dense, 6, __func__);
538 UNPROTECT(1);
539 new = 0;
540 }
541 PROTECT(s_from);
542 PROTECT(s_value = Rf_coerceVector(s_value, tx));
543 } else {
544 /* defined in ./coerce.c : */
545 SEXP dense_as_kind(SEXP, const char *, char, int);
546#ifndef MATRIX_ENABLE_IMATRIX
547 if (tv == INTSXP) {
548 PROTECT(s_from = dense_as_kind(s_from, class, 'd', 0));
549 PROTECT(s_value = Rf_coerceVector(s_value, REALSXP));
550 } else {
551#endif
552 PROTECT(s_from = dense_as_kind(s_from, class, typeToKind(tv), 0));
553 PROTECT(s_value);
554#ifndef MATRIX_ENABLE_IMATRIX
555 }
556#endif
557 class = Matrix_class(s_from, valid_dense, 6, __func__);
558 new = 0;
559 }
560
561 s_from = dense_diag_set(s_from, class, s_value, new);
562 UNPROTECT(2);
563 return s_from;
564}
565
566SEXP R_sparse_diag_set(SEXP s_from, SEXP s_value)
567{
568 const char *class = Matrix_class(s_from, valid_sparse, 6, __func__);
569 SEXPTYPE tx = kindToType(class[0]), tv = TYPEOF(s_value);
570
571 switch (tv) {
572 case LGLSXP:
573 case INTSXP:
574 case REALSXP:
575 case CPLXSXP:
576 break;
577 default:
578 Rf_error(_("replacement diagonal has incompatible type \"%s\""),
579 Rf_type2char(tv));
580 break;
581 }
582
583 int *pdim = DIM(s_from), m = pdim[0], n = pdim[1],
584 r = (m < n) ? m : n;
585 if (XLENGTH(s_value) != 1 && XLENGTH(s_value) != r)
586 Rf_error(_("replacement diagonal has wrong length"));
587
588 if (tv <= tx) {
589 /* defined in ./coerce.c : */
590 SEXP sparse_as_general(SEXP, const char *);
591 if (class[1] == 's' && class[0] == 'z' && tv == tx &&
592 TRANS(s_from) == 'C') {
593 PROTECT(s_from = sparse_as_general(s_from, class));
594 class = Matrix_class(s_from, valid_sparse, 6, __func__);
595 UNPROTECT(1);
596 }
597 PROTECT(s_from);
598 PROTECT(s_value = Rf_coerceVector(s_value, tx));
599 } else {
600 /* defined in ./coerce.c : */
601 SEXP sparse_as_kind(SEXP, const char *, char);
602#ifndef MATRIX_ENABLE_IMATRIX
603 if (tv == INTSXP) {
604 PROTECT(s_from = sparse_as_kind(s_from, class, 'd'));
605 PROTECT(s_value = Rf_coerceVector(s_value, REALSXP));
606 } else {
607#endif
608 PROTECT(s_from = sparse_as_kind(s_from, class, typeToKind(tv)));
609 PROTECT(s_value);
610#ifndef MATRIX_ENABLE_IMATRIX
611 }
612#endif
613 class = Matrix_class(s_from, valid_sparse, 6, __func__);
614 }
615
616 s_from = sparse_diag_set(s_from, class, s_value);
617 UNPROTECT(2);
618 return s_from;
619}
620
621SEXP sparse_diag_U2N(SEXP from, const char *class)
622{
623 if (class[1] != 't' || DIAG(from) == 'N')
624 return from;
625 SEXP value = PROTECT(Rf_ScalarLogical(1)),
626 to = R_sparse_diag_set(from, value);
627 UNPROTECT(1); /* value */
628 return to;
629}
630
631SEXP sparse_diag_N2U(SEXP from, const char *class)
632{
633 /* defined in ./band.c : */
634 SEXP sparse_band(SEXP, const char *, int, int);
635
636 if (class[1] != 't' || DIAG(from) != 'N')
637 return from;
638 SEXP to;
639 int n = DIM(from)[1];
640 if (n == 0) {
641 PROTECT(to = newObject(class));
642 if (UPLO(from) != 'U')
643 SET_UPLO(to);
644 }
645 else if (UPLO(from) == 'U')
646 PROTECT(to = sparse_band(from, class, 1, n));
647 else
648 PROTECT(to = sparse_band(from, class, -n, -1));
649 SET_DIAG(to);
650 UNPROTECT(1); /* to */
651 return from;
652}
653
654SEXP R_sparse_diag_U2N(SEXP s_from)
655{
656 const char *class = Matrix_class(s_from, valid_sparse, 2, __func__);
657 return sparse_diag_U2N(s_from, class);
658}
659
660SEXP R_sparse_diag_N2U(SEXP s_from)
661{
662 const char *class = Matrix_class(s_from, valid_sparse, 2, __func__);
663 return sparse_diag_N2U(s_from, class);
664}
665
666SEXP denseCholesky_diag_get(SEXP s_trf, SEXP s_root)
667{
668 SEXP x = PROTECT(GET_SLOT(s_trf, Matrix_xSym));
669 int n = DIM(s_trf)[1];
670 char ul = UPLO(s_trf);
671 SEXP ans = Rf_allocVector(TYPEOF(x), n);
672 if (n > 0) {
673 int j, root = Rf_asLogical(s_root),
674 packed = XLENGTH(x) != (int_fast64_t) n * n;
675 R_xlen_t dx = (!packed) ? (R_xlen_t) n + 1 : (ul == 'U') ? 1 : (R_xlen_t) n + 1,
676 ddx = (!packed) ? 0 : (ul == 'U') ? 1 : -1;
677 if (TYPEOF(x) == CPLXSXP) {
678 Rcomplex *pa = COMPLEX(ans), *px = COMPLEX(x);
679 for (j = 0; j < n; ++j) {
680 (*pa).r = (*px).r;
681 (*pa).i = 0.0;
682 if (!root)
683 (*pa).r *= (*pa).r;
684 pa += 1;
685 px += (dx += ddx);
686 }
687 } else {
688 double *pa = REAL(ans), *px = REAL(x);
689 for (j = 0; j < n; ++j) {
690 *pa = *px;
691 if (!root)
692 *pa *= *pa;
693 pa += 1;
694 px += (dx += ddx);
695 }
696 }
697 }
698 UNPROTECT(1);
699 return ans;
700}
701
702SEXP sparseCholesky_diag_get(SEXP s_trf, SEXP s_root)
703{
704 cholmod_factor *L = M2CHF(s_trf, 1);
705 int n = (int) L->n;
706 SEXP ans = Rf_allocVector((L->xtype == CHOLMOD_COMPLEX) ? CPLXSXP : REALSXP, n);
707 if (n > 0) {
708 int j, root = Rf_asLogical(s_root);
709 if (L->is_super) {
710 int k, nc,
711 nsuper = (int) L->nsuper,
712 *psuper = (int *) L->super,
713 *ppi = (int *) L->pi,
714 *ppx = (int *) L->px;
715 R_xlen_t nr1a;
716 if (L->xtype == CHOLMOD_COMPLEX) {
717 Rcomplex *pa = COMPLEX(ans), *px = (Rcomplex *) L->x, *py;
718 for (k = 0; k < nsuper; ++k) {
719 nc = psuper[k + 1] - psuper[k];
720 nr1a = (R_xlen_t) (ppi[k + 1] - ppi[k]) + 1;
721 py = px + ppx[k];
722 for (j = 0; j < nc; ++j) {
723 (*pa).r = (*py).r;
724 (*pa).i = 0.0;
725 if (!root)
726 (*pa).r *= (*pa).r;
727 pa += 1;
728 py += nr1a;
729 }
730 }
731 } else {
732 double *pa = REAL(ans), *px = (double *) L->x, *py;
733 for (k = 0; k < nsuper; ++k) {
734 nc = psuper[k + 1] - psuper[k];
735 nr1a = (R_xlen_t) (ppi[k + 1] - ppi[k]) + 1;
736 py = px + ppx[k];
737 for (j = 0; j < nc; ++j) {
738 *pa = *py;
739 if (!root)
740 *pa *= *pa;
741 pa += 1;
742 py += nr1a;
743 }
744 }
745 }
746 } else {
747 int *pp = (int *) L->p;
748 if (L->xtype == CHOLMOD_COMPLEX) {
749 Rcomplex *pa = COMPLEX(ans), *px = (Rcomplex *) L->x;
750 if (L->is_ll) {
751 for (j = 0; j < n; ++j) {
752 pa[j].r = px[pp[j]].r;
753 pa[j].i = 0.0;
754 if (!root)
755 pa[j].r *= pa[j].r;
756 }
757 } else {
758 for (j = 0; j < n; ++j) {
759 pa[j].r = px[pp[j]].r;
760 pa[j].i = 0.0;
761 if (root)
762 pa[j].r = (pa[j].r >= 0.0) ? sqrt(pa[j].r) : R_NaN;
763 }
764 }
765 } else {
766 double *pa = REAL(ans), *px = (double *) L->x;
767 if (L->is_ll) {
768 for (j = 0; j < n; ++j) {
769 pa[j] = px[pp[j]];
770 if (!root)
771 pa[j] *= pa[j];
772 }
773 } else {
774 for (j = 0; j < n; ++j) {
775 pa[j] = px[pp[j]];
776 if (root)
777 pa[j] = (pa[j] >= 0.0) ? sqrt(pa[j]) : R_NaN;
778 }
779 }
780 }
781 }
782 }
783 return ans;
784}
dense_as_kind
SEXP dense_as_kind(SEXP, const char *, char, int)
Definition coerce.c:2000
sparse_as_general
SEXP sparse_as_general(SEXP, const char *)
Definition coerce.c:2298
sparse_as_kind
SEXP sparse_as_kind(SEXP, const char *, char)
Definition coerce.c:2076
SWITCH5
#define SWITCH5(c, template)
Definition M5.h:327
SWITCH4
#define SWITCH4(c, template)
Definition M5.h:315
SET_DIM
#define SET_DIM(x, m, n)
Definition Mdefines.h:87
SWAP
#define SWAP(a, b, t, op)
Definition Mdefines.h:138
allocUnit
SEXP allocUnit(SEXPTYPE, R_xlen_t)
Definition utils.c:94
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
duplicateVector
SEXP duplicateVector(SEXP)
Definition utils.c:42
typeToKind
char typeToKind(SEXPTYPE)
Definition objects.c:20
Matrix_class
const char * Matrix_class(SEXP, const char **, int, const char *)
Definition objects.c:112
kindToType
SEXPTYPE kindToType(char)
Definition objects.c:38
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
SET_DIAG
#define SET_DIAG(x)
Definition Mdefines.h:113
naToUnit
void naToUnit(SEXP)
Definition utils.c:151
DIM
#define DIM(x)
Definition Mdefines.h:85
GET_SLOT
#define GET_SLOT(x, name)
Definition Mdefines.h:72
equalString
int equalString(SEXP, SEXP, R_xlen_t)
Definition utils.c:28
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
VALID_LOGIC2
#define VALID_LOGIC2(s, d)
Definition Mdefines.h:216
TYPEOF
#define TYPEOF(s)
Definition Mdefines.h:123
sparse_band
SEXP sparse_band(SEXP from, const char *class, int a, int b)
Definition band.c:95
M2CHF
cholmod_factor * M2CHF(SEXP obj, int values)
Definition cholmod-etc.c:39
dense_as_general
SEXP dense_as_general(SEXP from, const char *class, int new)
Definition coerce.c:2237
R_sparse_diag_set
SEXP R_sparse_diag_set(SEXP s_from, SEXP s_value)
Definition diag.c:566
R_dense_diag_get
SEXP R_dense_diag_get(SEXP s_obj, SEXP s_names)
Definition diag.c:205
R_dense_diag_set
SEXP R_dense_diag_set(SEXP s_from, SEXP s_value)
Definition diag.c:508
sparse_diag_N2U
SEXP sparse_diag_N2U(SEXP from, const char *class)
Definition diag.c:631
dense_diag_set
SEXP dense_diag_set(SEXP from, const char *class, SEXP value, int new)
Definition diag.c:225
sparse_diag_U2N
SEXP sparse_diag_U2N(SEXP from, const char *class)
Definition diag.c:621
TEMPLATE
#define TEMPLATE(c)
dense_diag_get
SEXP dense_diag_get(SEXP obj, const char *class, int names)
Definition diag.c:9
R_sparse_diag_U2N
SEXP R_sparse_diag_U2N(SEXP s_from)
Definition diag.c:654
sparseCholesky_diag_get
SEXP sparseCholesky_diag_get(SEXP s_trf, SEXP s_root)
Definition diag.c:702
R_sparse_diag_N2U
SEXP R_sparse_diag_N2U(SEXP s_from)
Definition diag.c:660
sparse_diag_set
SEXP sparse_diag_set(SEXP from, const char *class, SEXP value)
Definition diag.c:271
denseCholesky_diag_get
SEXP denseCholesky_diag_get(SEXP s_trf, SEXP s_root)
Definition diag.c:666
R_sparse_diag_get
SEXP R_sparse_diag_get(SEXP s_obj, SEXP s_names)
Definition diag.c:215
sparse_diag_get
SEXP sparse_diag_get(SEXP obj, const char *class, int names)
Definition diag.c:80
zvreal
void zvreal(Rcomplex *x, const Rcomplex *y, size_t n)
Definition idz.c:1274
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