Matrix r5059
Loading...
Searching...
No Matches
subscript.c
Go to the documentation of this file.
1/* C implementation of methods for [, [[ */
2
3#include "Mdefines.h"
4#include "M5.h"
5#include "idz.h"
6
7static
8int anyNA(int *p, int n)
9{
10 if (p && n > 0)
11 while (n--)
12 if (*(p++) == NA_INTEGER)
13 return 1;
14 return 0;
15}
16
17#define SUB1(c) \
18do { \
19 if (TYPEOF(s) == INTSXP) \
20 SUB1__(i, c); \
21 else \
22 SUB1__(d, c); \
23} while (0)
24
25#define iOOB(s, t, mn) (s == NA_INTEGER || (t = s - 1) >= mn)
26#define dOOB(s, t, mn) (ISNAN(s) || s >= 0x1.0p+63 || (t = (int_fast64_t) s - 1) >= mn)
27
28#define lCAST(x) (x)
29#define iCAST(x) (x)
30#define dCAST(x) (x)
31#define zCAST(x) (x)
32
33#undef nCAST
34#define nCAST(x) (x != 0)
35
36static
37SEXP dense_subscript_1ary(SEXP obj, const char *class, SEXP s)
38{
39 R_xlen_t slen = XLENGTH(s);
40 SEXP ans = Rf_allocVector(kindToType(class[0]), slen);
41 if (slen == 0)
42 return ans;
43 PROTECT(ans);
44
45 int *pdim = DIM(obj), m = pdim[0], n = pdim[1];
46 int_fast64_t mn = (int_fast64_t) m * n;
47
48 int packed = class[2] == 'p';
49 char ul = '\0', ct = '\0', nu = '\0';
50 if (class[1] != 'g')
51 ul = UPLO(obj);
52 if (class[1] == 's' && class[0] == 'z')
53 ct = TRANS(obj);
54 if (class[1] == 't')
55 nu = DIAG(obj);
56
57 SEXP x = GET_SLOT(obj, Matrix_xSym);
58 R_xlen_t l;
59 int_fast64_t b, i, j, k;
60
61 int ge = class[1] == 'g', sy = class[1] == 's', he = sy && ct == 'C',
62 un = class[1] == 't' && nu != 'N', up = ul == 'U';
63
64#define SUB1__(d, c) \
65 do { \
66 d##TYPE *ps = d##PTR(s); \
67 c##TYPE *pa = c##PTR(ans); \
68 c##TYPE *px = c##PTR(x); \
69 for (l = 0; l < slen; ++l) { \
70 if (d##OOB(ps[l], k, mn)) { \
71 pa[l] = c##NA; \
72 continue; \
73 } \
74 if (ge) { \
75 c##ASSIGN_IDEN(pa[l], c##CAST(px[k])); \
76 continue; \
77 } \
78 i = k % m; \
79 j = k / m; \
80 b = (up) ? j - i : i - j; \
81 if (!sy) { \
82 if (b < 0) { \
83 pa[l] = c##ZERO; \
84 continue; \
85 } \
86 if (un && b == 0) { \
87 pa[l] = c##UNIT; \
88 continue; \
89 } \
90 } \
91 if (!packed) { \
92 if (b >= 0) \
93 /* do nothing */; \
94 else \
95 k = i * m + j; \
96 } else if (up) { \
97 if (b >= 0) \
98 k = DENSE_INDEX_U(i, j, m); \
99 else \
100 k = DENSE_INDEX_U(j, i, m); \
101 } else { \
102 if (b >= 0) \
103 k = DENSE_INDEX_L(i, j, m); \
104 else \
105 k = DENSE_INDEX_L(j, i, m); \
106 } \
107 if (!he || b > 0) \
108 c##ASSIGN_IDEN (pa[l], c##CAST(px[k])); \
109 else if (b == 0) \
110 c##ASSIGN_PROJ_REAL(pa[l], c##CAST(px[k])); \
111 else \
112 c##ASSIGN_CONJ (pa[l], c##CAST(px[k])); \
113 } \
114 } while (0)
115
116 SWITCH5(class[0], SUB1);
117
118#undef SUB1__
119
120 UNPROTECT(1); /* ans */
121 return ans;
122}
123
124#undef nCAST
125#define nCAST(x) (1)
126
127static
128SEXP sparse_subscript_1ary(SEXP obj, const char *class, SEXP s, SEXP o)
129{
130 R_xlen_t slen = XLENGTH(s);
131 SEXP ans = Rf_allocVector(kindToType(class[0]), slen);
132 if (slen == 0)
133 return ans;
134 PROTECT(ans);
135
136 if (class[2] == 'T') {
137 /* defined in ./coerce.c : */
138 SEXP sparse_as_Csparse(SEXP, const char *);
139 obj = sparse_as_Csparse(obj, class);
140 }
141 PROTECT(obj);
142
143 int *pdim = DIM(obj), m = pdim[0], n = pdim[1];
144 int_fast64_t mn = (int_fast64_t) m * n;
145
146 char ul = '\0', ct = '\0', nu = '\0';
147 if (class[1] != 'g')
148 ul = UPLO(obj);
149 if (class[1] == 's' && class[0] == 'z')
150 ct = TRANS(obj);
151 if (class[1] == 't')
152 nu = DIAG(obj);
153
154 SEXP iSym = (class[2] != 'R') ? Matrix_iSym : Matrix_jSym,
155 p_ = PROTECT(GET_SLOT(obj, Matrix_pSym)),
156 i_ = PROTECT(GET_SLOT(obj, iSym));
157 int *pp = INTEGER(p_), *pi = INTEGER(i_), j_, k_, kend_;
158 pp++;
159
160 int *po_i = (TYPEOF(o) == INTSXP) ? INTEGER(o) : NULL;
161 double *po_d = (TYPEOF(o) == REALSXP) ? REAL(o) : NULL;
162
163 R_xlen_t l = 0, l_;
164 int b = 0, i, j;
165 int_fast64_t k;
166
167 int byrow = class[2] == 'R',
168 sy = class[1] == 's', he = sy && ct == 'C',
169 un = class[1] == 't' && nu != 'N', up = (!byrow) == (ul == 'U');
170
171#define MAP(l) \
172 (po_i) ? (R_xlen_t) po_i[l] - 1 : (po_d) ? (R_xlen_t) po_d[l] - 1 : l
173
174#define ADV(d) \
175 do { \
176 if (l >= slen) \
177 j = -1; \
178 else { \
179 l_ = MAP(l); \
180 if (d##OOB(ps[l_], k, mn)) \
181 j = -1; \
182 else { \
183 if (byrow) { \
184 i = (int) (k / m); \
185 j = (int) (k % m); \
186 } \
187 else { \
188 i = (int) (k % m); \
189 j = (int) (k / m); \
190 } \
191 if (sy && (b = (up) ? j - i : i - j) < 0) \
192 SWAP(i, j, int, ); \
193 ++l; \
194 } \
195 } \
196 } while (0)
197
198#define SUB1__(d, c) \
199 do { \
200 d##TYPE *ps = d##PTR(s); \
201 c##TYPE *pa = c##PTR(ans); \
202 c##IF_NPATTERN( \
203 SEXP x = GET_SLOT(obj, Matrix_xSym); \
204 c##TYPE *px = c##PTR(x); \
205 ); \
206 ADV(d); \
207 while (j >= 0) { \
208 j_ = j; \
209 k_ = pp[j_ - 1]; \
210 kend_ = pp[j_]; \
211 while (k_ < kend_ && j == j_) { \
212 if (pi[k_] < i) \
213 ++k_; \
214 else { \
215 if (pi[k_] > i) \
216 pa[l_] = (un && i == j) ? c##UNIT : c##ZERO; \
217 else if (!he || b > 0) \
218 c##ASSIGN_IDEN (pa[l_], c##CAST(px[k_])); \
219 else if (b == 0) \
220 c##ASSIGN_PROJ_REAL(pa[l_], c##CAST(px[k_])); \
221 else \
222 c##ASSIGN_CONJ (pa[l_], c##CAST(px[k_])); \
223 ADV(d); \
224 } \
225 } \
226 while (j == j_) { \
227 pa[l_] = (un && i == j) ? c##UNIT : c##ZERO; \
228 ADV(d); \
229 } \
230 } \
231 while (l < slen) { \
232 l_ = MAP(l); \
233 pa[l_] = c##NA; \
234 ++l; \
235 } \
236 } while (0)
237
238 SWITCH5(class[0], SUB1);
239
240#undef MAP
241#undef ADV
242#undef SUB1__
243
244 UNPROTECT(4); /* i, p, obj, ans */
245 return ans;
246}
247
248#undef nCAST
249#define nCAST(x) (x != 0)
250
251static
252SEXP diagonal_subscript_1ary(SEXP obj, const char *class, SEXP s)
253{
254 R_xlen_t slen = XLENGTH(s);
255 SEXP ans = Rf_allocVector(kindToType(class[0]), slen);
256 if (slen == 0)
257 return ans;
258 PROTECT(ans);
259
260 int n = DIM(obj)[1];
261 int_fast64_t nn = (int_fast64_t) n * n, n1a = (int_fast64_t) n + 1;
262
263 int un = DIAG(obj) != 'N';
264
265 SEXP x = GET_SLOT(obj, Matrix_xSym);
266 R_xlen_t l;
267 int_fast64_t k;
268
269
270#define SUB1__(d, c) \
271 do { \
272 d##TYPE *ps = d##PTR(s); \
273 c##TYPE *pa = c##PTR(ans); \
274 c##TYPE *px = c##PTR(x); \
275 for (l = 0; l < slen; ++l) { \
276 if (d##OOB(ps[l], k, nn)) \
277 pa[l] = c##NA; \
278 else if (k % n1a) \
279 pa[l] = c##ZERO; \
280 else if (un) \
281 pa[l] = c##UNIT; \
282 else \
283 pa[l] = c##CAST(px[k / n1a]); \
284 } \
285 } while (0)
286
287 SWITCH5(class[0], SUB1);
288
289#undef SUB1__
290
291 UNPROTECT(1); /* ans */
292 return ans;
293}
294
295static
296SEXP index_subscript_1ary(SEXP obj, const char *class, SEXP s)
297{
298 R_xlen_t slen = XLENGTH(s);
299 SEXP ans = Rf_allocVector(LGLSXP, slen);
300 if (slen == 0)
301 return ans;
302 PROTECT(ans);
303 int *pa = LOGICAL(ans);
304
305 int *pdim = DIM(obj), m = pdim[0], n = pdim[1];
306 int_fast64_t mn = (int_fast64_t) m * n;
307
308 SEXP perm = PROTECT(GET_SLOT(obj, Matrix_permSym));
309 int *pperm = INTEGER(perm), mg = MARGIN(obj);
310
311 R_xlen_t l;
312 int_fast64_t k;
313
314#define SUB1__(d, c) \
315 do { \
316 d##TYPE *ps = d##PTR(s); \
317 for (l = 0; l < slen; ++l) { \
318 if (d##OOB(ps[l], k, mn)) \
319 pa[l] = NA_LOGICAL; \
320 else if (mg == 0) \
321 pa[l] = k / m == pperm[k % m] - 1; \
322 else \
323 pa[l] = k % m == pperm[k / m] - 1; \
324 } \
325 } while (0)
326
327 SUB1();
328
329#undef SUB1__
330
331 UNPROTECT(2); /* perm, ans */
332 return ans;
333}
334
335/* x[s] where 's' is a vector of type "integer" or "double" */
336/* with elements greater than or equal to 1 (or NA) */
337/* NB: for [CRT], the caller must provide a permutation 'o' of */
338/* the processed subscript op(s) sorting the latter by row */
339/* then by column [R] or by column then by row [CT]. */
340/* For general and triangular matrices, 'op' is an identity. */
341/* For symmetric matrices: */
342/* */
343/* / i+j*m if (i,j) is " on-side" */
344/* op(i+j*m) = */
345/* \ j+i*m if (i,j) is "off-side" */
346/* */
347/* o=NULL <==> o=seq_along(s) */
348SEXP R_subscript_1ary(SEXP s_obj, SEXP s_s, SEXP s_o)
349{
350 const char *class = Matrix_class(s_obj, valid_matrix, 7, __func__);
351 validObject(s_obj, class);
352 switch (class[2]) {
353 case 'e':
354 case 'y':
355 case 'r':
356 case 'p':
357 return dense_subscript_1ary(s_obj, class, s_s);
358 case 'C':
359 case 'R':
360 case 'T':
361 return sparse_subscript_1ary(s_obj, class, s_s, s_o);
362 case 'i':
363 return diagonal_subscript_1ary(s_obj, class, s_s);
364 case 'd':
365 return index_subscript_1ary(s_obj, class, s_s);
366 default:
367 return R_NilValue;
368 }
369}
370
371#undef SUB1
372
373#undef nCAST
374#define nCAST(x) (x != 0)
375
376static
377SEXP dense_subscript_1ary_2col(SEXP obj, const char *class, SEXP s)
378{
379 int slen = (int) (XLENGTH(s) / 2);
380 SEXP ans = Rf_allocVector(kindToType(class[0]), slen);
381 if (slen == 0)
382 return ans;
383 PROTECT(ans);
384 int m = DIM(obj)[0], *ps0 = INTEGER(s), *ps1 = ps0 + slen;
385
386 int packed = class[2] == 'p';
387 char ul = '\0', ct = '\0', nu = '\0';
388 if (class[1] != 'g')
389 ul = UPLO(obj);
390 if (class[1] == 's' && class[0] == 'z')
391 ct = TRANS(obj);
392 if (class[1] == 't')
393 nu = DIAG(obj);
394
395 SEXP x = GET_SLOT(obj, Matrix_xSym);
396 int l;
397 int_fast64_t b, i, j, k;
398
399 int ge = class[1] == 'g', sy = class[1] == 's', he = sy && ct == 'C',
400 un = class[1] == 't' && nu != 'N', up = ul == 'U';
401
402#define SUB1(c) \
403 do { \
404 c##TYPE *pa = c##PTR(ans); \
405 c##TYPE *px = c##PTR(x); \
406 for (l = 0; l < slen; ++l) { \
407 if (ps0[l] == NA_INTEGER || ps1[l] == NA_INTEGER) { \
408 pa[l] = c##NA; \
409 continue; \
410 } \
411 i = (int_fast64_t) ps0[l] - 1; \
412 j = (int_fast64_t) ps1[l] - 1; \
413 if (ge) { \
414 c##ASSIGN_IDEN(pa[l], c##CAST(px[DENSE_INDEX_N(i, j, m)])); \
415 continue; \
416 } \
417 b = (up) ? j - i : i - j; \
418 if (!sy) { \
419 if (b < 0) { \
420 pa[l] = c##ZERO; \
421 continue; \
422 } \
423 if (un && b == 0) { \
424 pa[l] = c##UNIT; \
425 continue; \
426 } \
427 } \
428 if (!packed) { \
429 if (b >= 0) \
430 k = j * m + i; \
431 else \
432 k = i * m + j; \
433 } else if (up) { \
434 if (b >= 0) \
435 k = DENSE_INDEX_U(i, j, m); \
436 else \
437 k = DENSE_INDEX_U(j, i, m); \
438 } else { \
439 if (b >= 0) \
440 k = DENSE_INDEX_L(i, j, m); \
441 else \
442 k = DENSE_INDEX_L(j, i, m); \
443 } \
444 if (!he || b > 0) \
445 c##ASSIGN_IDEN (pa[l], c##CAST(px[k])); \
446 else if (b == 0) \
447 c##ASSIGN_PROJ_REAL(pa[l], c##CAST(px[k])); \
448 else \
449 c##ASSIGN_CONJ (pa[l], c##CAST(px[k])); \
450 } \
451 } while (0)
452
453 SWITCH5(class[0], SUB1);
454
455#undef SUB1
456
457 UNPROTECT(1); /* ans */
458 return ans;
459}
460
461#undef nCAST
462#define nCAST(x) (1)
463
464static
465SEXP sparse_subscript_1ary_2col(SEXP obj, const char *class, SEXP s, SEXP o)
466{
467 int slen = (int) (XLENGTH(s) / 2);
468 SEXP ans = Rf_allocVector(kindToType(class[0]), slen);
469 if (slen == 0)
470 return ans;
471 PROTECT(ans);
472 int *ps0 = INTEGER(s), *ps1 = ps0 + slen;
473
474 if (class[2] == 'T') {
475 /* defined in ./coerce.c : */
476 SEXP sparse_as_Csparse(SEXP, const char *);
477 obj = sparse_as_Csparse(obj, class);
478 }
479 PROTECT(obj);
480
481 char ul = '\0', ct = '\0', nu = '\0';
482 if (class[1] != 'g')
483 ul = UPLO(obj);
484 if (class[1] == 's' && class[0] == 'z')
485 ct = TRANS(obj);
486 if (class[1] == 't')
487 nu = DIAG(obj);
488
489 SEXP iSym = (class[2] != 'R') ? Matrix_iSym : Matrix_jSym,
490 p_ = PROTECT(GET_SLOT(obj, Matrix_pSym)),
491 i_ = PROTECT(GET_SLOT(obj, iSym));
492 int *pp = INTEGER(p_), *pi = INTEGER(i_), j_, k_, kend_;
493 pp++;
494
495 int *po_i = (TYPEOF(o) == INTSXP) ? INTEGER(o) : NULL;
496 double *po_d = (TYPEOF(o) == REALSXP) ? REAL(o) : NULL;
497
498 int l = 0, l_;
499 int b = 0, i, j;
500
501 int byrow = class[2] == 'R',
502 sy = class[1] == 's', he = sy && ct == 'C',
503 un = class[1] == 't' && nu != 'N', up = (!byrow) == (ul == 'U');
504
505#define MAP(l) \
506 (po_i) ? (int) po_i[l] - 1 : (po_d) ? (int) po_d[l] - 1 : l
507
508#define ADV \
509 do { \
510 if (l >= slen) \
511 j = -1; \
512 else { \
513 l_ = MAP(l); \
514 if (ps0[l_] == NA_INTEGER || ps1[l_] == NA_INTEGER) \
515 j = -1; \
516 else { \
517 if (byrow) { \
518 i = ps1[l_] - 1; \
519 j = ps0[l_] - 1; \
520 } \
521 else { \
522 i = ps0[l_] - 1; \
523 j = ps1[l_] - 1; \
524 } \
525 if (sy && (b = (up) ? j - i : i - j) < 0) \
526 SWAP(i, j, int, ); \
527 ++l; \
528 } \
529 } \
530 } while (0)
531
532#define SUB1(c) \
533 do { \
534 c##TYPE *pa = c##PTR(ans); \
535 c##IF_NPATTERN( \
536 SEXP x = GET_SLOT(obj, Matrix_xSym); \
537 c##TYPE *px = c##PTR(x); \
538 ); \
539 ADV; \
540 while (j >= 0) { \
541 j_ = j; \
542 k_ = pp[j_ - 1]; \
543 kend_ = pp[j_]; \
544 while (k_ < kend_ && j == j_) { \
545 if (pi[k_] < i) \
546 ++k_; \
547 else { \
548 if (pi[k_] > i) \
549 pa[l_] = (un && i == j) ? c##UNIT : c##ZERO; \
550 else if (!he || b > 0) \
551 c##ASSIGN_IDEN (pa[l_], c##CAST(px[k_])); \
552 else if (b == 0) \
553 c##ASSIGN_PROJ_REAL(pa[l_], c##CAST(px[k_])); \
554 else \
555 c##ASSIGN_CONJ (pa[l_], c##CAST(px[k_])); \
556 ADV; \
557 } \
558 } \
559 while (j == j_) { \
560 pa[l_] = (un && i == j) ? c##UNIT : c##ZERO; \
561 ADV; \
562 } \
563 } \
564 while (l < slen) { \
565 l_ = MAP(l); \
566 pa[l_] = c##NA; \
567 ++l; \
568 } \
569 } while (0)
570
571 SWITCH5(class[0], SUB1);
572
573#undef MAP
574#undef ADV
575#undef SUB1
576
577 UNPROTECT(4); /* i, p, obj, ans */
578 return ans;
579}
580
581#undef nCAST
582#define nCAST(x) (x != 0)
583
584static
585SEXP diagonal_subscript_1ary_2col(SEXP obj, const char *class, SEXP s)
586{
587 int slen = (int) (XLENGTH(s) / 2);
588 SEXP ans = Rf_allocVector(kindToType(class[0]), slen);
589 if (slen == 0)
590 return ans;
591 PROTECT(ans);
592 int *ps0 = INTEGER(s), *ps1 = ps0 + slen;
593
594 int un = DIAG(obj) != 'N';
595
596 SEXP x = GET_SLOT(obj, Matrix_xSym);
597 int l;
598
599#define SUB1(c) \
600 do { \
601 c##TYPE *pa = c##PTR(ans); \
602 c##TYPE *px = c##PTR(x); \
603 for (l = 0; l < slen; ++l) { \
604 if (ps0[l] == NA_INTEGER || ps1[l] == NA_INTEGER) \
605 pa[l] = c##NA; \
606 else if (ps0[l] != ps1[l]) \
607 pa[l] = c##ZERO; \
608 else if (un) \
609 pa[l] = c##UNIT; \
610 else \
611 pa[l] = c##CAST(px[ps0[l] - 1]); \
612 } \
613 } while (0)
614
615 SWITCH5(class[0], SUB1);
616
617#undef SUB1
618
619 UNPROTECT(1); /* ans */
620 return ans;
621}
622
623static
624SEXP index_subscript_1ary_2col(SEXP obj, const char *class, SEXP s)
625{
626 int slen = (int) (XLENGTH(s) / 2);
627 SEXP ans = Rf_allocVector(LGLSXP, slen);
628 if (slen == 0)
629 return ans;
630 PROTECT(ans);
631 int *ps0 = INTEGER(s), *ps1 = ps0 + slen, *pa = LOGICAL(ans);
632
633 SEXP perm = PROTECT(GET_SLOT(obj, Matrix_permSym));
634 int *pperm = INTEGER(perm), mg = MARGIN(obj);
635
636 int l;
637
638 for (l = 0; l < slen; ++l) {
639 if (ps0[l] == NA_INTEGER || ps1[l] == NA_INTEGER)
640 pa[l] = NA_LOGICAL;
641 else if (mg == 0)
642 pa[l] = ps1[l] == pperm[ps0[l]] - 1;
643 else
644 pa[l] = ps0[l] == pperm[ps1[l]] - 1;
645 }
646
647 UNPROTECT(2); /* perm, ans */
648 return ans;
649}
650
651/* x[s] where 's' is a 2-column matrix of type "integer" */
652/* with s[, 1L] in 1:m (or NA) and s[, 2L] in 1:n (or NA) */
653/* NB: for [CRT], the caller must provide a permutation 'o' of */
654/* the processed subscript op(s) sorting the latter by row */
655/* then by column [R] or by column then by row [CT]. */
656/* For general and triangular matrices, 'op' is an identity. */
657/* For symmetric matrices: */
658/* */
659/* / (i,j) if (i,j) is " on-side" */
660/* op((i,j)) = */
661/* \ (j,i) if (i,j) is "off-side" */
662/* */
663/* o=NULL <==> o=seq_len(nrow(s)) */
664SEXP R_subscript_1ary_2col(SEXP s_obj, SEXP s_s, SEXP s_o)
665{
666 const char *class = Matrix_class(s_obj, valid_matrix, 7, __func__);
667 validObject(s_obj, class);
668 switch (class[2]) {
669 case 'e':
670 case 'y':
671 case 'r':
672 case 'p':
673 return dense_subscript_1ary_2col(s_obj, class, s_s);
674 case 'C':
675 case 'R':
676 case 'T':
677 return sparse_subscript_1ary_2col(s_obj, class, s_s, s_o);
678 case 'i':
679 return diagonal_subscript_1ary_2col(s_obj, class, s_s);
680 case 'd':
681 return index_subscript_1ary_2col(s_obj, class, s_s);
682 default:
683 return R_NilValue;
684 }
685}
686
687static
688int stay_sy(int *pi, int ni, int *pj, int nj, int n, char uplo, int checkNA)
689{
690 if (!pi || !pj || ni != nj ||
691 memcmp(pi, pj, sizeof(int) * (size_t) ni) != 0)
692 return 0;
693 int ki, r = (uplo == 'U') ? 1 : -1;
694 if (checkNA && anyNA(pi, ni))
695 return r;
696 if (ni >= 2) {
697 /* triangular iff monotone */
698 if (pi[0] == pi[1])
699 return r;
700 else if (pi[0] < pi[1]) {
701 for (ki = 2; ki < ni; ++ki)
702 if (pi[ki - 1] >= pi[ki])
703 return r;
704 /* up->up, lo->lo */
705 } else {
706 for (ki = 2; ki < ni; ++ki)
707 if (pi[ki - 1] <= pi[ki])
708 return r;
709 /* up->lo, lo->up */
710 r = -r;
711 }
712 }
713 return 2 * r;
714}
715
716static
717int stay_tr(int *pi, int ni, int *pj, int nj, int n, char uplo, int checkNA)
718{
719 if (!pi || !pj || ni != nj)
720 return 0;
721 int ki, kj, iden = memcmp(pi, pj, sizeof(int) * (size_t) ni) == 0,
722 r = (uplo == 'U') ? 1 : -1;
723 if (checkNA && (anyNA(pi, ni) || (!iden && anyNA(pj, ni))))
724 return 0;
725 if (iden) {
726 /* triangular iff monotone; unit diagonal is preserved */
727 if (ni >= 2) {
728 if (pi[0] == pi[1])
729 return 0;
730 else if (pi[0] < pi[1]) {
731 for (ki = 2; ki < ni; ++ki)
732 if (pi[ki - 1] >= pi[ki])
733 return 0;
734 /* up->up, lo->lo */
735 } else {
736 for (ki = 2; ki < ni; ++ki)
737 if (pi[ki - 1] <= pi[ki])
738 return 0;
739 /* up->lo, lo->up */
740 r = -r;
741 }
742 }
743 return r * 2;
744 } else {
745 /* brute force ... */
746 int j;
747 if (uplo == 'U') {
748 for (kj = 0; kj < nj; ++kj)
749 for (ki = kj + 1, j = pj[kj]; ki < ni; ++ki)
750 if (pi[ki] <= j)
751 goto LO;
752 /* up->up */
753 return r;
754 LO:
755 for (kj = 0; kj < nj; ++kj)
756 for (ki = 0, j = pj[kj]; ki < kj; ++ki)
757 if (pi[ki] <= j)
758 return 0;
759 /* up->lo */
760 return -r;
761 } else {
762 for (kj = 0; kj < nj; ++kj)
763 for (ki = 0, j = pj[kj]; ki < kj; ++ki)
764 if (pi[ki] >= j)
765 goto UP;
766 /* lo->lo */
767 return r;
768 UP:
769 for (kj = 0; kj < nj; ++kj)
770 for (ki = kj + 1, j = pj[kj]; ki < ni; ++ki)
771 if (pi[ki] >= j)
772 return 0;
773 /* lo->up */
774 return -r;
775 }
776 }
777}
778
779static
780int stay_di(int *pi, int ni, int *pj, int nj, int n, int checkNA)
781{
782 if (!pi || !pj || ni != nj)
783 return 0;
784 int ki, kj, iden = memcmp(pi, pj, sizeof(int) * (size_t) ni) == 0;
785 if (checkNA && (anyNA(pi, ni) || (!iden && anyNA(pj, ni))))
786 return 0;
787 if (iden) {
788 /* diagonal iff no duplicates; unit diagonal is preserved */
789 char *iwork;
790 Matrix_Calloc(iwork, n, char);
791 --iwork;
792 for (ki = 0; ki < ni; ++ki) {
793 if (iwork[pi[ki]])
794 return 0;
795 iwork[pi[ki]] = 1;
796 }
797 ++iwork;
798 Matrix_Free(iwork, n);
799 return 2;
800 } else {
801 /* brute force ... */
802 int j;
803 for (kj = 0; kj < nj; ++kj) {
804 j = pj[kj];
805 for (ki = 0; ki < kj; ++ki)
806 if (pi[ki] == j)
807 return 0;
808 for (ki = kj + 1; ki < ni; ++ki)
809 if (pi[ki] == j)
810 return 0;
811 }
812 return 1;
813 }
814}
815
816static
817SEXP dense_subscript_2ary(SEXP obj, const char *class, SEXP si, SEXP sj)
818{
819 if (si == R_NilValue && sj == R_NilValue)
820 return obj;
821
822 int *pdim = DIM(obj), m = pdim[0], n = pdim[1],
823 ki, mi = si == R_NilValue, ni = (mi) ? m : LENGTH(si),
824 kj, mj = sj == R_NilValue, nj = (mj) ? n : LENGTH(sj),
825 *pi = (mi) ? NULL : INTEGER(si),
826 *pj = (mj) ? NULL : INTEGER(sj);
827
828 int packed = class[2] == 'p';
829 char ul0 = '\0', ct0 = '\0', nu0 = '\0';
830 if (class[1] != 'g')
831 ul0 = UPLO(obj);
832 if (class[1] == 's' && class[0] == 'z')
833 ct0 = TRANS(obj);
834 if (class[1] == 't')
835 nu0 = DIAG(obj);
836
837 int stay = (class[1] == 'g') ? 0 : (class[1] == 's')
838 ? stay_sy(pi, ni, pj, nj, n, ul0, 1)
839 : stay_tr(pi, ni, pj, nj, n, ul0, 1);
840 int_fast64_t ninj = (int_fast64_t) ni * nj,
841 xlen = (!packed || stay == 0) ? ninj : ni + (ninj - ni) / 2;
842 if (xlen > R_XLEN_T_MAX)
843 Rf_error(_("attempt to allocate vector of length exceeding %s"),
844 "R_XLEN_T_MAX");
845
846 char ul1 = '\0', ct1 = '\0', nu1 = '\0';
847 if (stay != 0) {
848 if (class[1] != 'g')
849 ul1 = (stay > 0) ? 'U' : 'L';
850 if (class[1] == 's' && class[0] == 'z')
851 ct1 = ct0;
852 if (class[1] == 't')
853 nu1 = (ABS(stay) == 1) ? 'N' : nu0;
854 }
855
856 int he = class[1] == 's' && ct0 == 'C' && ct1 != 'C',
857 un = class[1] == 't' && nu0 != 'N';
858
859 char cl[] = "...Matrix";
860 cl[0] = class[0];
861 cl[1] = (stay == 0) ? 'g' : class[1];
862 cl[2] = (stay == 0) ? 'e' : class[2];
863 SEXP ans = PROTECT(newObject(cl));
864
865 SET_DIM(ans, ni, nj);
866 if (cl[1] != 'g' && ul1 != 'U')
867 SET_UPLO(ans);
868 if (cl[1] == 's' && ct1 != 'C' && cl[0] == 'z')
869 SET_TRANS(ans);
870 if (cl[1] == 't' && nu1 != 'N')
871 SET_DIAG(ans);
872
873 SEXP x0 = PROTECT(GET_SLOT(obj, Matrix_xSym)),
874 x1 = PROTECT(Rf_allocVector(TYPEOF(x0), (R_xlen_t) xlen));
875 SET_SLOT(ans, Matrix_xSym, x1);
876
877 int_fast64_t i, j;
878
879#define SUB2(c) \
880 do { \
881 c##TYPE *px0 = c##PTR(x0), *px1 = c##PTR(x1); \
882 if (!packed && cl[1] != 'g') \
883 memset(px1, 0, sizeof(c##TYPE) * (size_t) XLENGTH(x1)); \
884 if (class[1] == 'g') { \
885 SUB2__(c, ASSIGN_GE, N, N, for (ki = 0; ki < ni; ++ki), , ); \
886 } else if (class[1] == 's' && !packed) { \
887 if (ul1 == '\0') { \
888 if (ul0 == 'U') \
889 SUB2__(c, ASSIGN_SY, U, N, for (ki = 0; ki < ni; ++ki), , ); \
890 else \
891 SUB2__(c, ASSIGN_SY, L, N, for (ki = 0; ki < ni; ++ki), , ); \
892 } else if (ul1 == 'U') { \
893 if (ul0 == 'U') \
894 SUB2__(c, ASSIGN_SY, U, N, for (ki = 0; ki <= kj; ++ki), \
895 , px1 += ni - kj - 1); \
896 else \
897 SUB2__(c, ASSIGN_SY, L, N, for (ki = 0; ki <= kj; ++ki), \
898 , px1 += ni - kj - 1); \
899 } else { \
900 if (ul0 == 'U') \
901 SUB2__(c, ASSIGN_SY, U, N, for (ki = kj; ki < ni; ++ki), \
902 px1 += kj, ); \
903 else \
904 SUB2__(c, ASSIGN_SY, L, N, for (ki = kj; ki < ni; ++ki), \
905 px1 += kj, ); \
906 } \
907 } else if (class[1] == 's') { \
908 if (ul1 == '\0') { \
909 if (ul0 == 'U') \
910 SUB2__(c, ASSIGN_SY, U, U, for (ki = 0; ki < ni; ++ki), , ); \
911 else \
912 SUB2__(c, ASSIGN_SY, L, L, for (ki = 0; ki < ni; ++ki), , ); \
913 } else if (ul1 == 'U') { \
914 if (ul0 == 'U') \
915 SUB2__(c, ASSIGN_SY, U, U, for (ki = 0; ki <= kj; ++ki), , ); \
916 else \
917 SUB2__(c, ASSIGN_SY, L, L, for (ki = 0; ki <= kj; ++ki), , ); \
918 } else { \
919 if (ul0 == 'U') \
920 SUB2__(c, ASSIGN_SY, U, U, for (ki = kj; ki < ni; ++ki), , ); \
921 else \
922 SUB2__(c, ASSIGN_SY, L, L, for (ki = kj; ki < ni; ++ki), , ); \
923 } \
924 } else if (!packed) { \
925 if (ul1 == '\0') { \
926 if (ul0 == 'U') \
927 SUB2__(c, ASSIGN_TR, U, N, for (ki = 0; ki < ni; ++ki), , ); \
928 else \
929 SUB2__(c, ASSIGN_TR, L, N, for (ki = 0; ki < ni; ++ki), , ); \
930 } else if (ul1 == 'U') { \
931 if (ul0 == 'U') \
932 SUB2__(c, ASSIGN_TR, U, N, for (ki = 0; ki <= kj; ++ki), \
933 , px1 += ni - kj - 1); \
934 else \
935 SUB2__(c, ASSIGN_TR, L, N, for (ki = 0; ki <= kj; ++ki), \
936 , px1 += ni - kj - 1); \
937 } else { \
938 if (ul0 == 'U') \
939 SUB2__(c, ASSIGN_TR, U, N, for (ki = kj; ki < ni; ++ki), \
940 px1 += kj, ); \
941 else \
942 SUB2__(c, ASSIGN_TR, L, N, for (ki = kj; ki < ni; ++ki), \
943 px1 += kj, ); \
944 } \
945 } else { \
946 if (ul1 == '\0') { \
947 if (ul0 == 'U') \
948 SUB2__(c, ASSIGN_TR, U, U, for (ki = 0; ki < ni; ++ki), , ); \
949 else \
950 SUB2__(c, ASSIGN_TR, L, L, for (ki = 0; ki < ni; ++ki), , ); \
951 } else if (ul1 == 'U') { \
952 if (ul0 == 'U') \
953 SUB2__(c, ASSIGN_TR, U, U, for (ki = 0; ki <= kj; ++ki), , ); \
954 else \
955 SUB2__(c, ASSIGN_TR, L, L, for (ki = 0; ki <= kj; ++ki), , ); \
956 } else { \
957 if (ul0 == 'U') \
958 SUB2__(c, ASSIGN_TR, U, U, for (ki = kj; ki < ni; ++ki), , ); \
959 else \
960 SUB2__(c, ASSIGN_TR, L, L, for (ki = kj; ki < ni; ++ki), , ); \
961 } \
962 } \
963 } while (0)
964
965#define SUB2__(c, assign, s, t, __for__, jump0, jump1) \
966 do { \
967 for (kj = 0; kj < nj; ++kj) { \
968 if (mj) \
969 j = kj; \
970 else if (pj[kj] != NA_INTEGER) \
971 j = pj[kj] - 1; \
972 else { \
973 jump0; \
974 __for__ { \
975 *(px1++) = c##NA; \
976 } \
977 jump1; \
978 continue; \
979 } \
980 jump0; \
981 __for__ { \
982 if (mi) \
983 i = ki; \
984 else if (pi[ki] != NA_INTEGER) \
985 i = pi[ki] - 1; \
986 else { \
987 *(px1++) = c##NA; \
988 continue; \
989 } \
990 assign(c, px0, px1, i, j, m, s, t); \
991 px1++; \
992 } \
993 jump1; \
994 } \
995 } while (0)
996
997#define ASSIGN_GE(c, x, y, i, j, m, s, t) \
998 do { \
999 c##ASSIGN_IDEN(*y, x[DENSE_INDEX_N(i, j, m)]); \
1000 } while (0)
1001
1002#define ASSIGN_SY(c, x, y, i, j, m, s, t) \
1003 do { \
1004 if (he && i == j) \
1005 c##ASSIGN_PROJ_REAL(*y, x[DENSE_INDEX_##t(i, j, m)]); \
1006 else if (CMP_##s(i, j)) \
1007 c##ASSIGN_IDEN (*y, x[DENSE_INDEX_##t(i, j, m)]); \
1008 else if (he) \
1009 c##ASSIGN_CONJ (*y, x[DENSE_INDEX_##t(j, i, m)]); \
1010 else \
1011 c##ASSIGN_IDEN (*y, x[DENSE_INDEX_##t(j, i, m)]); \
1012 } while (0)
1013
1014#define ASSIGN_TR(c, x, y, i, j, m, s, t) \
1015 do { \
1016 if (un && i == j) \
1017 c##ASSIGN_IDEN(*y, c##UNIT); \
1018 else if (CMP_##s(i, j)) \
1019 c##ASSIGN_IDEN(*y, x[DENSE_INDEX_##t(i, j, m)]); \
1020 else \
1021 c##ASSIGN_IDEN(*y, c##ZERO); \
1022 } while (0)
1023
1024#define CMP_U(i, j) i <= j
1025#define CMP_L(i, j) i >= j
1026
1027 SWITCH4(class[0], SUB2);
1028
1029#undef SUB2__
1030#undef SUB2
1031
1032 UNPROTECT(3); /* x1, x0, ans */
1033 return ans;
1034}
1035
1036static
1037SEXP sparse_subscript_2ary(SEXP obj, const char *class, SEXP si, SEXP sj)
1038{
1039 if (si == R_NilValue && sj == R_NilValue)
1040 return obj;
1041
1042 int mg = (class[2] == 'R') ? 0 : 1;
1043 if (mg == 0)
1044 SWAP(si, sj, SEXP, );
1045
1046 int *pdim = DIM(obj), m = pdim[mg == 0], n = pdim[mg != 0],
1047 ki, mi = si == R_NilValue, ni = (mi) ? m : LENGTH(si),
1048 kj, mj = sj == R_NilValue, nj = (mj) ? n : LENGTH(sj),
1049 *pi = (mi) ? NULL : INTEGER(si),
1050 *pj = (mj) ? NULL : INTEGER(sj);
1051 if (anyNA(pi, ni) || anyNA(pj, nj))
1052 Rf_error(_("NA subscripts in %s not supported for '%s' inheriting from %s"),
1053 "x[i, j]", "x", "sparseMatrix");
1054
1055 char ul0 = '\0', ct0 = '\0', nu0 = '\0';
1056 if (class[1] != 'g') {
1057 ul0 = UPLO(obj);
1058 if (mg == 0)
1059 ul0 = (ul0 == 'U') ? 'L' : 'U';
1060 }
1061 if (class[1] == 's' && class[0] == 'z')
1062 ct0 = TRANS(obj);
1063 if (class[1] == 't')
1064 nu0 = DIAG(obj);
1065
1066 int stay = (class[1] == 'g') ? 0 : (class[1] == 's')
1067 ? stay_sy(pi, ni, pj, nj, n, ul0, 0)
1068 : stay_tr(pi, ni, pj, nj, n, ul0, 0);
1069
1070 char ul1 = '\0', ct1 = '\0', nu1 = '\0';
1071 if (stay != 0) {
1072 if (class[1] != 'g')
1073 ul1 = (stay > 0) ? 'U' : 'L';
1074 if (class[1] == 's' && class[0] == 'z')
1075 ct1 = ct0;
1076 if (class[1] == 't')
1077 nu1 = (ABS(stay) == 1) ? 'N' : nu0;
1078 }
1079
1080 char class__[] = "...Matrix";
1081 class__[0] = class[0];
1082 class__[1] = class[1];
1083 class__[2] = class[2];
1084 PROTECT_INDEX pid_obj;
1085 PROTECT_WITH_INDEX(obj, &pid_obj);
1086 if (class[2] == 'T') {
1087 /* defined in ./coerce.c : */
1088 SEXP sparse_as_Csparse(SEXP, const char *);
1089 REPROTECT(obj = sparse_as_Csparse(obj, class__), pid_obj);
1090 class__[2] = 'C';
1091 }
1092 if (class[1] != 'g' && ABS(stay) <= 1) {
1093 /* defined in ./coerce.c : */
1094 SEXP sparse_as_general(SEXP, const char *);
1095 REPROTECT(obj = sparse_as_general(obj, class__), pid_obj);
1096 class__[1] = 'g';
1097 }
1098
1099 char cl[] = "...Matrix";
1100 cl[0] = class [0];
1101 cl[1] = class [1];
1102 cl[2] = class__[2];
1103 if (ABS(stay) <= ((class[1] != 's') ? 0 : 1))
1104 cl[1] = 'g';
1105 PROTECT_INDEX pid_ans;
1106 SEXP ans;
1107 PROTECT_WITH_INDEX(ans = newObject(cl), &pid_ans);
1108
1109 SET_DIM(ans, (mg == 0) ? nj : ni, (mg == 0) ? ni : nj);
1110 if (cl[1] != 'g' && (mg == (ul1 != 'U')))
1111 SET_UPLO(ans);
1112 if (cl[1] == 's' && ct1 != 'C' && cl[0] == 'z')
1113 SET_TRANS(ans);
1114 if (cl[1] == 't' && nu1 != 'N')
1115 SET_DIAG(ans);
1116
1117 SEXP iSym = (class[2] != 'R') ? Matrix_iSym : Matrix_jSym,
1118 p0 = PROTECT(GET_SLOT(obj, Matrix_pSym)),
1119 i0 = PROTECT(GET_SLOT(obj, iSym)),
1120 p1 = PROTECT(Rf_allocVector(INTSXP, (R_xlen_t) nj + 1));
1121 int *pp0 = INTEGER(p0), *pi0 = INTEGER(i0), *pp1 = INTEGER(p1),
1122 i, j, k, kend;
1123 pp0++; *(pp1++) = 0;
1124 SET_SLOT(ans, Matrix_pSym, p1);
1125
1126 if (mi) {
1127
1128 int_fast64_t nnz = 0;
1129 for (kj = 0; kj < nj; ++kj) {
1130 j = pj[kj] - 1;
1131 nnz += pp0[j] - pp0[j - 1];
1132 pp1[kj] = (int) nnz;
1133 }
1134 if (nnz > INT_MAX)
1135 Rf_error(_("%s too dense for %s; would have more than %s nonzero entries"),
1136 "x[i, j]", "[CR]sparseMatrix", "2^31-1");
1137
1138 SEXP i1 = PROTECT(Rf_allocVector(INTSXP, (int) nnz));
1139 int *pi1 = INTEGER(i1);
1140 SET_SLOT(ans, iSym, i1);
1141
1142#define SUB2(c) \
1143 do { \
1144 c##IF_NPATTERN( \
1145 SEXP x0 = PROTECT(GET_SLOT(obj, Matrix_xSym)), \
1146 x1 = PROTECT(Rf_allocVector(c##TYPESXP, (int) nnz)); \
1147 c##TYPE *px0 = c##PTR(x0), *px1 = c##PTR(x1); \
1148 SET_SLOT(ans, Matrix_xSym, x1); \
1149 UNPROTECT(2); /* x1, x0 */ \
1150 ); \
1151 for (kj = 0; kj < nj; ++kj) { \
1152 j = pj[kj] - 1; \
1153 k = pp0[j - 1]; \
1154 kend = pp0[j]; \
1155 while (k < kend) { \
1156 *(pi1++) = pi0[k]; \
1157 c##IF_NPATTERN( \
1158 *(px1++) = px0[k]; \
1159 ); \
1160 ++k; \
1161 } \
1162 } \
1163 } while (0)
1164
1165 SWITCH5(class[0], SUB2);
1166
1167#undef SUB2
1168
1169 UNPROTECT(1); /* i1 */
1170
1171 } else {
1172
1173 int *iwork;
1174 size_t liwork = (size_t) ((int_fast64_t) m + m + ni);
1175 Matrix_Calloc(iwork, liwork, int);
1176
1177 int *iworkA = iwork, *iworkB = iworkA + m, *iworkC = iworkB + m,
1178 prev = m, sort = 0;
1179
1180 /* workA[ i]: size of set { ki : pi[ki] - 1 == i } */
1181 /* workB[ i]: minimum of set { ki : pi[ki] - 1 == i } */
1182 /* workC[ki]: minimum of set { ki' : ki' > ki, pi[ki'] == pi[ki] } */
1183
1184 for (ki = ni - 1; ki >= 0; --ki) {
1185 i = pi[ki] - 1;
1186 ++iworkA[i];
1187 iworkC[ki] = iworkB[i];
1188 iworkB[i] = ki;
1189 if (i > prev)
1190 sort = 1;
1191 prev = i;
1192 }
1193
1194 int_fast64_t nnz = 0;
1195 for (kj = 0; kj < nj; ++kj) {
1196 j = (mj) ? kj : pj[kj] - 1;
1197 k = pp0[j - 1];
1198 kend = pp0[j];
1199 while (k < kend) {
1200 nnz += iworkA[pi0[k]];
1201 ++k;
1202 }
1203 pp1[kj] = (int) nnz;
1204 }
1205 if (nnz > INT_MAX)
1206 Rf_error(_("%s too dense for %s; would have more than %s nonzero entries"),
1207 "x[i, j]", "[CR]sparseMatrix", "2^31-1");
1208
1209 SEXP i1 = PROTECT(Rf_allocVector(INTSXP, (int) nnz));
1210 int *pi1 = INTEGER(i1), d;
1211 SET_SLOT(ans, iSym, i1);
1212
1213#define SUB2(c) \
1214 do { \
1215 c##IF_NPATTERN( \
1216 SEXP x0 = PROTECT(GET_SLOT(obj, Matrix_xSym)), \
1217 x1 = PROTECT(Rf_allocVector(c##TYPESXP, (int) nnz)); \
1218 c##TYPE *px0 = c##PTR(x0), *px1 = c##PTR(x1); \
1219 SET_SLOT(ans, Matrix_xSym, x1); \
1220 UNPROTECT(2); /* x1, x0 */ \
1221 ); \
1222 for (kj = 0; kj < nj; ++kj) { \
1223 j = (mj) ? kj : pj[kj] - 1; \
1224 k = pp0[j - 1]; \
1225 kend = pp0[j]; \
1226 while (k < kend) { \
1227 i = pi0[k]; \
1228 d = iworkA[i]; \
1229 ki = iworkB[i]; \
1230 while (d--) { \
1231 *(pi1++) = ki; \
1232 c##IF_NPATTERN( \
1233 *(px1++) = px0[k]; \
1234 ); \
1235 ki = iworkC[ki]; \
1236 } \
1237 ++k; \
1238 } \
1239 } \
1240 } while (0)
1241
1242 SWITCH5(class[0], SUB2);
1243
1244#undef SUB2
1245
1246 Matrix_Free(iwork, liwork);
1247
1248 if (sort) {
1249 liwork = (size_t) ((int_fast64_t) ni + 1 + ((ni < nj) ? nj : ni) + nnz);
1250 Matrix_Calloc(iwork, liwork, int);
1251
1252#define TEMPLATE(c) \
1253 do { \
1254 c##TYPE *px1 = NULL, *work = NULL; \
1255 c##IF_NPATTERN( \
1256 SEXP x1 = PROTECT(GET_SLOT(ans, Matrix_xSym)); \
1257 px1 = c##PTR(x1); \
1258 Matrix_Calloc(work, nnz, c##TYPE); \
1259 ); \
1260 c##cspsort(INTEGER(p1), INTEGER(i1), px1, ni, nj, iwork, work); \
1261 c##IF_NPATTERN( \
1262 Matrix_Free(work, nnz); \
1263 UNPROTECT(1); /* x1 */ \
1264 ); \
1265 } while (0)
1266
1267 SWITCH5(class[0], TEMPLATE);
1268
1269#undef TEMPLATE
1270
1271 Matrix_Free(iwork, liwork);
1272 }
1273
1274 UNPROTECT(1); /* i1 */
1275
1276 }
1277
1278 if (class[1] == 's' && ABS(stay) == 1) {
1279 /* defined in ./forceSymmetric.c : */
1280 SEXP sparse_force_symmetric(SEXP, const char *, char, char);
1281 REPROTECT(ans = sparse_force_symmetric(ans, cl, (mg == (ul1 == 'U')) ? 'U' : 'L', ct1), pid_ans);
1282 cl[1] = 's';
1283 }
1284 if (class[2] == 'T') {
1285 /* defined in ./coerce.c : */
1286 SEXP sparse_as_Tsparse(SEXP, const char *);
1287 REPROTECT(ans = sparse_as_Tsparse(ans, cl), pid_ans);
1288 cl[2] = 'T';
1289 }
1290
1291 UNPROTECT(5); /* p1, i0, p0, ans, obj */
1292 return ans;
1293}
1294
1295static
1296SEXP diagonal_subscript_2ary(SEXP obj, const char *class, SEXP si, SEXP sj)
1297{
1298 if (si == R_NilValue && sj == R_NilValue)
1299 return obj;
1300
1301 int n = DIM(obj)[1],
1302 ki, mi = si == R_NilValue, ni = (mi) ? n : LENGTH(si),
1303 kj, mj = sj == R_NilValue, nj = (mj) ? n : LENGTH(sj),
1304 *pi = (mi) ? NULL : INTEGER(si),
1305 *pj = (mj) ? NULL : INTEGER(sj);
1306 if (anyNA(pi, ni) || anyNA(pj, nj))
1307 Rf_error(_("NA subscripts in %s not supported for '%s' inheriting from %s"),
1308 "x[i, j]", "x", "sparseMatrix");
1309 int stay = stay_di(pi, ni, pj, nj, n, 0);
1310
1311 char nu0 = DIAG(obj),
1312 nu1 = (stay <= 0) ? '\0' : (stay <= 1) ? 'N' : nu0;
1313
1314 int un = nu0 != 'N';
1315
1316 char cl[] = "...Matrix";
1317 cl[0] = class[0];
1318 cl[1] = (stay <= 0) ? 'g' : 'd';
1319 cl[2] = (stay <= 0) ? 'C' : 'i';
1320 SEXP ans = PROTECT(newObject(cl));
1321
1322 SET_DIM(ans, ni, nj);
1323
1324 if (nu1 != '\0' && nu1 != 'N') {
1325
1326 SET_DIAG(ans);
1327
1328 } else if (nu1 != '\0') {
1329
1330 SEXP x0 = PROTECT(GET_SLOT(obj, Matrix_xSym)),
1331 x1 = PROTECT(Rf_allocVector(TYPEOF(x0), ni));
1332 SET_SLOT(ans, Matrix_xSym, x1);
1333
1334 int j;
1335
1336#define SUB2(c) \
1337 do { \
1338 c##TYPE *px0 = c##PTR(x0) - 1, *px1 = c##PTR(x1); \
1339 for (ki = 0; ki < ni; ++ki) \
1340 if (*(pi++) != (j = *(pj++))) \
1341 *(px1++) = c##ZERO; \
1342 else \
1343 *(px1++) = (un) ? c##UNIT : px0[j]; \
1344 } while (0)
1345
1346 SWITCH4(class[0], SUB2);
1347
1348#undef SUB2
1349
1350 UNPROTECT(2); /* x1, x0 */
1351
1352 } else {
1353
1354 SEXP x0 = PROTECT(GET_SLOT(obj, Matrix_xSym));
1355
1356 SEXP p1 = PROTECT(Rf_allocVector(INTSXP, (R_xlen_t) nj + 1));
1357 int *pp1 = INTEGER(p1), j;
1358 *(pp1++) = 0;
1359 SET_SLOT(ans, Matrix_pSym, p1);
1360
1361#define TEMPLATE(c) \
1362 do { \
1363 c##TYPE *px0 = c##PTR(x0); \
1364 for (kj = 0; kj < nj; ++kj) { \
1365 j = (mj) ? kj : pj[kj] - 1; \
1366 pp1[kj] = 0; \
1367 if (un || c##NOT_ZERO(px0[j])) { \
1368 if (mi) { \
1369 for (ki = 0; ki < ni; ++ki) \
1370 if (ki == j) \
1371 ++pp1[kj]; \
1372 } else { \
1373 for (ki = 0; ki < ni; ++ki) \
1374 if (pi[ki] - 1 == j) \
1375 ++pp1[kj]; \
1376 } \
1377 if (pp1[kj] > INT_MAX - pp1[kj - 1]) \
1378 break; \
1379 } \
1380 pp1[kj] += pp1[kj - 1]; \
1381 } \
1382 } while (0)
1383
1384 kj = -1;
1385 SWITCH4(class[0], TEMPLATE);
1386
1387#undef TEMPLATE
1388
1389 if (kj < nj)
1390 Rf_error(_("%s too dense for %s; would have more than %s nonzero entries"),
1391 "x[i, j]", "[CR]sparseMatrix", "2^31-1");
1392
1393 SEXP i1 = PROTECT(Rf_allocVector(INTSXP, pp1[nj - 1]));
1394 int *pi1 = INTEGER(i1);
1395 SET_SLOT(ans, Matrix_iSym, i1);
1396
1397#define SUB2(c) \
1398 do { \
1399 c##TYPE *px0 = c##PTR(x0); \
1400 c##IF_NPATTERN( \
1401 SEXP x1 = PROTECT(Rf_allocVector(c##TYPESXP, pp1[nj - 1])); \
1402 c##TYPE *px1 = c##PTR(x1); \
1403 SET_SLOT(ans, Matrix_xSym, x1); \
1404 UNPROTECT(1); /* x1 */ \
1405 ); \
1406 for (kj = 0; kj < nj; ++kj) { \
1407 j = (mj) ? kj : pj[kj] - 1; \
1408 if (un || c##NOT_ZERO(px0[j])) { \
1409 if (mi) { \
1410 for (ki = 0; ki < ni; ++ki) \
1411 if (ki == j) { \
1412 *(pi1++) = ki; \
1413 c##IF_NPATTERN( \
1414 *(px1++) = (un) ? c##UNIT : px0[j]; \
1415 ); \
1416 } \
1417 } else { \
1418 for (ki = 0; ki < ni; ++ki) \
1419 if (pi[ki] - 1 == j) { \
1420 *(pi1++) = ki; \
1421 c##IF_NPATTERN( \
1422 *(px1++) = (un) ? c##UNIT : px0[j]; \
1423 ); \
1424 } \
1425 } \
1426 } \
1427 } \
1428 } while (0)
1429
1430 SWITCH4(class[0], SUB2);
1431
1432#undef SUB2
1433
1434 UNPROTECT(3); /* i1, p1, x0 */
1435
1436 }
1437
1438 UNPROTECT(1); /* ans */
1439 return ans;
1440}
1441
1442static
1443SEXP index_subscript_2ary(SEXP obj, const char *class, SEXP si, SEXP sj)
1444{
1445 if (si == R_NilValue && sj == R_NilValue)
1446 return obj;
1447
1448 /* defined in ./perm.c : */
1449 int isPerm(const int *, int, int);
1450 void invertPerm(const int *, int *, int, int, int);
1451
1452 int mg = MARGIN(obj);
1453 if (mg != 0)
1454 SWAP(si, sj, SEXP, );
1455
1456 int *pdim = DIM(obj), m = pdim[mg != 0], n = pdim[mg == 0],
1457 ki, mi = si == R_NilValue, ni = (mi) ? m : LENGTH(si),
1458 kj, mj = sj == R_NilValue, nj = (mj) ? n : LENGTH(sj),
1459 *pi = (mi) ? NULL : INTEGER(si),
1460 *pj = (mj) ? NULL : INTEGER(sj);
1461 if (anyNA(pi, ni) || anyNA(pj, nj))
1462 Rf_error(_("NA subscripts in %s not supported for '%s' inheriting from %s"),
1463 "x[i, j]", "x", "sparseMatrix");
1464 int stay = class[0] == 'p';
1465
1466 PROTECT_INDEX pid_perm;
1467 SEXP perm = GET_SLOT(obj, Matrix_permSym);
1468 int *pperm = INTEGER(perm);
1469 PROTECT_WITH_INDEX(perm, &pid_perm);
1470
1471 PROTECT_INDEX pid_ans;
1472 SEXP ans = obj;
1473 PROTECT_WITH_INDEX(ans, &pid_ans);
1474
1475 if (!mi) {
1476 stay = stay && ni == m && isPerm(pi, m, 1);
1477 REPROTECT(ans = newObject((stay) ? "pMatrix" : "indMatrix"), pid_ans);
1478
1479 SET_DIM(ans, (mg == 0) ? ni : n, (mg == 0) ? n : ni);
1480 SET_MARGIN(ans, mg);
1481
1482 int *tmp = pperm;
1483 REPROTECT(perm = Rf_allocVector(INTSXP, ni), pid_perm);
1484 pperm = INTEGER(perm);
1485
1486 --tmp; /* 1-indexed */
1487 for (ki = 0; ki < ni; ++ki)
1488 pperm[ki] = tmp[pi[ki]];
1489 ++tmp; /* 0-indexed */
1490
1491 SET_SLOT(ans, Matrix_permSym, perm);
1492 m = ni;
1493 }
1494
1495 if (!mj) {
1496 stay = stay && nj == n && isPerm(pj, n, 1);
1497 REPROTECT(ans = newObject((stay) ? "pMatrix" : (mg == 0) ? "ngCMatrix" : "ngRMatrix"), pid_ans);
1498
1499 SET_DIM(ans, (mg == 0) ? m : nj, (mg == 0) ? nj : m);
1500
1501 int *iwork;
1502 size_t liwork = (size_t) ((stay) ? nj : (int_fast64_t) n + 1 + n + m);
1503 Matrix_Calloc(iwork, liwork, int);
1504
1505 if (stay) {
1506
1507 invertPerm(pj, iwork, nj, 1, 1);
1508
1509 int *tmp = pperm;
1510 REPROTECT(perm = Rf_allocVector(INTSXP, nj), pid_perm);
1511 pperm = INTEGER(perm);
1512
1513 --iwork; /* 1-indexed */
1514 for (kj = 0; kj < nj; ++kj)
1515 pperm[kj] = iwork[tmp[kj]];
1516 ++iwork; /* 0-indexed */
1517
1518 SET_SLOT(ans, Matrix_permSym, perm);
1519
1520 } else {
1521
1522 int *pp0 = iwork, *pp_ = pp0 + n + 1, *pi0 = pp_ + n;
1523
1524 SEXP p1 = PROTECT(Rf_allocVector(INTSXP, (R_xlen_t) nj + 1));
1525 int *pp1 = INTEGER(p1), i, j, k, kend, k_;
1526 *(pp1++) = 0;
1527 SET_SLOT(ans, Matrix_pSym, p1);
1528
1529 /* 1. Compute old column counts */
1530 for (i = 0; i < m; ++i)
1531 ++pp0[pperm[i]];
1532
1533 /* 2. Compute new column pointers */
1534 int_fast64_t nnz = 0;
1535 for (kj = 0; kj < nj; ++kj) {
1536 nnz += pp0[pj[kj]];
1537 pp1[kj] = (int) nnz;
1538 }
1539 if (nnz > INT_MAX)
1540 Rf_error(_("%s too dense for %s; would have more than %s nonzero entries"),
1541 "x[i, j]", "[CR]sparseMatrix", "2^31-1");
1542
1543 SEXP i1 = PROTECT(Rf_allocVector(INTSXP, pp1[nj - 1]));
1544 int *pi1 = INTEGER(i1);
1545 SET_SLOT(ans, (mg == 0) ? Matrix_iSym : Matrix_jSym, i1);
1546
1547 /* 3. Compute old column pointers */
1548 for (j = 0; j < n; ++j)
1549 pp0[j + 1] += (pp_[j] = pp0[j]);
1550
1551 /* 4. Sort old row indices by column */
1552 --pp0;
1553 for (i = 0; i < m; ++i)
1554 pi0[pp0[pperm[i]]++] = i;
1555
1556 /* 5. Copy old row indices */
1557 --pp_;
1558 for (kj = 0, k = 0; kj < nj; ++kj) {
1559 k_ = pp_[pj[kj]];
1560 kend = pp1[kj];
1561 while (k < kend)
1562 pi1[k++] = pi0[k_++];
1563 }
1564
1565 UNPROTECT(2); /* i1, p1 */
1566
1567 }
1568
1569 Matrix_Free(iwork, liwork);
1570 n = nj;
1571 }
1572
1573 UNPROTECT(2); /* perm, ans */
1574 return ans;
1575}
1576
1577/* x[i, j, drop = FALSE] where 'i' and 'j' are vectors of type "integer" */
1578/* of length at most 2^31-1 with 'i' in 1:m (or NA) and 'j' in 1:n (or NA) */
1579/* ... _not_ handling 'Dimnames' here */
1580SEXP R_subscript_2ary(SEXP s_obj, SEXP s_si, SEXP s_sj)
1581{
1582 const char *class = Matrix_class(s_obj, valid_matrix, 6, __func__);
1583 validObject(s_obj, class);
1584 switch (class[2]) {
1585 case 'e':
1586 case 'y':
1587 case 'r':
1588 case 'p':
1589 return dense_subscript_2ary(s_obj, class, s_si, s_sj);
1590 case 'C':
1591 case 'R':
1592 case 'T':
1593 return sparse_subscript_2ary(s_obj, class, s_si, s_sj);
1594 case 'i':
1595 return diagonal_subscript_2ary(s_obj, class, s_si, s_sj);
1596 case 'd':
1597 case 'a':
1598 return index_subscript_2ary(s_obj, class, s_si, s_sj);
1599 default:
1600 return R_NilValue;
1601 }
1602}
sparse_as_general
SEXP sparse_as_general(SEXP, const char *)
Definition coerce.c:2298
sparse_force_symmetric
SEXP sparse_force_symmetric(SEXP, const char *, char, char)
Definition forceSymmetric.c:83
sparse_as_Csparse
SEXP sparse_as_Csparse(SEXP, const char *)
Definition coerce.c:2731
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
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
Matrix_class
const char * Matrix_class(SEXP, const char **, int, const char *)
Definition objects.c:112
kindToType
SEXPTYPE kindToType(char)
Definition objects.c:38
Matrix_Free
#define Matrix_Free(p, n)
Definition Mdefines.h:56
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
ABS
#define ABS(i)
Definition Mdefines.h:135
DIM
#define DIM(x)
Definition Mdefines.h:85
GET_SLOT
#define GET_SLOT(x, name)
Definition Mdefines.h:72
SET_MARGIN
#define SET_MARGIN(x, j)
Definition Mdefines.h:118
validObject
void validObject(SEXP, const char *)
Definition validity.c:2004
newObject
SEXP newObject(const char *)
Definition objects.c:13
valid_matrix
const char * valid_matrix[]
Definition Mdefines.h:333
MARGIN
#define MARGIN(x)
Definition Mdefines.h:116
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_as_Tsparse
SEXP sparse_as_Tsparse(SEXP from, const char *class)
Definition coerce.c:2963
Matrix_permSym
SEXP Matrix_permSym
Definition init.c:623
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
isPerm
int isPerm(const int *p, int n, int off)
Definition perm.c:3
invertPerm
void invertPerm(const int *p, int *ip, int n, int off, int ioff)
Definition perm.c:47
R_subscript_1ary
SEXP R_subscript_1ary(SEXP s_obj, SEXP s_s, SEXP s_o)
Definition subscript.c:348
stay_di
static int stay_di(int *pi, int ni, int *pj, int nj, int n, int checkNA)
Definition subscript.c:780
SUB2
#define SUB2(c)
stay_tr
static int stay_tr(int *pi, int ni, int *pj, int nj, int n, char uplo, int checkNA)
Definition subscript.c:717
dense_subscript_1ary_2col
static SEXP dense_subscript_1ary_2col(SEXP obj, const char *class, SEXP s)
Definition subscript.c:377
index_subscript_2ary
static SEXP index_subscript_2ary(SEXP obj, const char *class, SEXP si, SEXP sj)
Definition subscript.c:1443
dense_subscript_2ary
static SEXP dense_subscript_2ary(SEXP obj, const char *class, SEXP si, SEXP sj)
Definition subscript.c:817
stay_sy
static int stay_sy(int *pi, int ni, int *pj, int nj, int n, char uplo, int checkNA)
Definition subscript.c:688
diagonal_subscript_1ary
static SEXP diagonal_subscript_1ary(SEXP obj, const char *class, SEXP s)
Definition subscript.c:252
sparse_subscript_1ary_2col
static SEXP sparse_subscript_1ary_2col(SEXP obj, const char *class, SEXP s, SEXP o)
Definition subscript.c:465
TEMPLATE
#define TEMPLATE(c)
R_subscript_2ary
SEXP R_subscript_2ary(SEXP s_obj, SEXP s_si, SEXP s_sj)
Definition subscript.c:1580
dense_subscript_1ary
static SEXP dense_subscript_1ary(SEXP obj, const char *class, SEXP s)
Definition subscript.c:37
index_subscript_1ary_2col
static SEXP index_subscript_1ary_2col(SEXP obj, const char *class, SEXP s)
Definition subscript.c:624
SUB1
#define SUB1(c)
Definition subscript.c:17
sparse_subscript_2ary
static SEXP sparse_subscript_2ary(SEXP obj, const char *class, SEXP si, SEXP sj)
Definition subscript.c:1037
diagonal_subscript_2ary
static SEXP diagonal_subscript_2ary(SEXP obj, const char *class, SEXP si, SEXP sj)
Definition subscript.c:1296
sparse_subscript_1ary
static SEXP sparse_subscript_1ary(SEXP obj, const char *class, SEXP s, SEXP o)
Definition subscript.c:128
R_subscript_1ary_2col
SEXP R_subscript_1ary_2col(SEXP s_obj, SEXP s_s, SEXP s_o)
Definition subscript.c:664
index_subscript_1ary
static SEXP index_subscript_1ary(SEXP obj, const char *class, SEXP s)
Definition subscript.c:296
anyNA
static int anyNA(int *p, int n)
Definition subscript.c:8
diagonal_subscript_1ary_2col
static SEXP diagonal_subscript_1ary_2col(SEXP obj, const char *class, SEXP s)
Definition subscript.c:585
Generated by doxygen 1.12.0