Matrix r5059
Loading...
Searching...
No Matches
Summary.c
Go to the documentation of this file.
1/* C implementation of methods for sum, prod */
2
3#include "Mdefines.h"
4#include "M5.h"
5
6/* defined in ./sparse.c : */
7SEXP sparse_aggregate(SEXP, const char *);
8
9#define LONGDOUBLE_AS_DOUBLE(x) \
10(((x) > DBL_MAX) ? R_PosInf : (((x) < -DBL_MAX) ? R_NegInf : (double) (x)))
11
12#define TRY_INCREMENT(x, y, label) \
13do { \
14 if ((x >= 0) \
15 ? ( y <= INT_FAST64_MAX - x) \
16 : (-y <= x - INT_FAST64_MIN)) { \
17 x += y; \
18 y = 0; \
19 count = 0; \
20 } else \
21 goto label; \
22} while (0)
23
24SEXP dense_sum(SEXP obj, const char *class, int narm)
25{
26 int *pdim = DIM(obj), m = pdim[0], n = pdim[1];
27
28 char ul = '\0', ct = '\0', nu = '\0';
29 if (class[1] != 'g')
30 ul = UPLO(obj);
31 if (class[1] == 's' && class[0] == 'z')
32 ct = TRANS(obj);
33 if (class[1] == 't')
34 nu = DIAG(obj);
35
36 SEXP x = GET_SLOT(obj, Matrix_xSym),
37 ans = R_NilValue;
38 int i, j, packed = class[2] == 'p',
39 sy = class[1] == 's', he = sy && ct == 'C',
40 un = class[1] == 't' && nu != 'N';
41
42#define SUM \
43 do { \
44 if (class[1] == 'g') { \
45 for (j = 0; j < n; ++j) \
46 SUM_KERNEL(for (i = 0; i < m; ++i)); \
47 } else if (class[1] == 's' || nu == 'N') { \
48 if (ul == 'U') \
49 for (j = 0; j < n; ++j) { \
50 SUM_KERNEL(for (i = 0; i <= j; ++i)); \
51 if (!packed) \
52 px += n - j - 1; \
53 } \
54 else \
55 for (j = 0; j < n; ++j) { \
56 if (!packed) \
57 px += j; \
58 SUM_KERNEL(for (i = j; i < n; ++i)); \
59 } \
60 } else { \
61 if (ul == 'U') \
62 for (j = 0; j < n; ++j) { \
63 SUM_KERNEL(for (i = 0; i < j; ++i)); \
64 px += 1; \
65 if (!packed) \
66 px += n - j - 1; \
67 } \
68 else \
69 for (j = 0; j < n; ++j) { \
70 if (!packed) \
71 px += j; \
72 px += 1; \
73 SUM_KERNEL(for (i = j + 1; i < n; ++i)); \
74 } \
75 } \
76 } while (0)
77
78 switch (class[0]) {
79 case 'n':
80 {
81 int *px = LOGICAL(x);
82 int_fast64_t s = (un) ? n : 0;
83
84#define SUM_KERNEL(__for__) \
85 do { \
86 __for__ { \
87 if (*px != 0) \
88 s += (sy && i != j) ? 2 : 1; \
89 px += 1; \
90 } \
91 } while (0)
92
93 SUM;
94
95#undef SUM_KERNEL
96
97 if (s <= INT_MAX)
98 ans = Rf_ScalarInteger((int) s);
99 else
100 ans = Rf_ScalarReal((double) s);
101 break;
102 }
103 case 'l':
104 case 'i':
105 {
106 int *px = (class[0] == 'l') ? LOGICAL(x) : INTEGER(x);
107 int_fast64_t s = (un) ? n : 0, t = 0;
108 unsigned int count = 0;
109
110#define SUM_KERNEL(__for__) \
111 do { \
112 __for__ { \
113 if (*px != NA_INTEGER) { \
114 unsigned int d = (sy && i != j) ? 2 : 1; \
115 if (count > UINT_MAX - d) \
116 TRY_INCREMENT(s, t, over); \
117 t += (sy && i != j) ? 2LL * *px : *px; \
118 count += d; \
119 } \
120 else if (!narm) \
121 return Rf_ScalarInteger(NA_INTEGER); \
122 px += 1; \
123 } \
124 } while (0)
125
126 SUM;
127
128#undef SUM_KERNEL
129
130 TRY_INCREMENT(s, t, over);
131
132 if (s > INT_MIN && s <= INT_MAX)
133 ans = Rf_ScalarInteger((int) s);
134 else
135 ans = Rf_ScalarReal((double) s);
136 break;
137over:
138 px = (class[0] == 'l') ? LOGICAL(x) : INTEGER(x);
139 long double lr = (un) ? (long double) n : 0.0;
140
141#define SUM_KERNEL(__for__) \
142 do { \
143 __for__ { \
144 if (*px != NA_INTEGER) \
145 lr += (sy && i != j) ? 2.0L * *px : *px; \
146 else if (!narm) \
147 return Rf_ScalarInteger(NA_INTEGER); \
148 px += 1; \
149 } \
150 } while (0)
151
152 SUM;
153
154#undef SUM_KERNEL
155
156 ans = Rf_ScalarReal(LONGDOUBLE_AS_DOUBLE(lr));
157 break;
158 }
159 case 'd':
160 {
161 double *px = REAL(x);
162 long double lr = (un) ? (long double) n : 0.0;
163
164#define SUM_KERNEL(__for__) \
165 do { \
166 __for__ { \
167 if (!(narm && ISNAN(*px))) \
168 lr += (sy && i != j) ? 2.0L * *px : *px; \
169 px += 1; \
170 } \
171 } while (0)
172
173 SUM;
174
175#undef SUM_KERNEL
176
177 ans = Rf_ScalarReal(LONGDOUBLE_AS_DOUBLE(lr));
178 break;
179 }
180 case 'z':
181 {
182 Rcomplex *px = COMPLEX(x), tmp;
183 long double lr = (un) ? (long double) n : 0.0;
184 long double li = 0.0;
185
186#define SUM_KERNEL(__for__) \
187 do { \
188 __for__ { \
189 if (!(narm && (ISNAN((*px).r) || (!he && ISNAN((*px).i))))) { \
190 lr += (sy && i != j) ? 2.0L * (*px).r : (*px).r; \
191 if (!he) \
192 li += (sy && i != j) ? 2.0L * (*px).i : (*px).i; \
193 } \
194 px += 1; \
195 } \
196 } while (0)
197
198 SUM;
199
200#undef SUM_KERNEL
201
202 tmp.r = LONGDOUBLE_AS_DOUBLE(lr);
203 tmp.i = LONGDOUBLE_AS_DOUBLE(li);
204 ans = Rf_ScalarComplex(tmp);
205 break;
206 }
207 default:
208 break;
209 }
210
211#undef SUM
212
213 return ans;
214}
215
216SEXP sparse_sum(SEXP obj, const char *class, int narm)
217{
218 PROTECT(obj = sparse_aggregate(obj, class));
219
220 int *pdim = DIM(obj), m = pdim[0], n = pdim[1];
221
222 char ul = '\0', ct = '\0', nu = '\0';
223 if (class[1] != 'g')
224 ul = UPLO(obj);
225 if (class[1] == 's' && class[0] == 'z')
226 ct = TRANS(obj);
227 if (class[1] == 't')
228 nu = DIAG(obj);
229
230 SEXP x = PROTECT((class[0] == 'n') ? R_NilValue : GET_SLOT(obj, Matrix_xSym)),
231 ans = R_NilValue;
232 int sy = class[1] == 's', he = sy && ct == 'C',
233 un = class[1] == 't' && nu != 'N';
234
235 if (class[2] != 'T') {
236
237 if (class[2] == 'R')
238 SWAP(m, n, int, );
239
240 SEXP iSym = (class[2] == 'C') ? Matrix_iSym : Matrix_jSym,
241 p = PROTECT(GET_SLOT(obj, Matrix_pSym)),
242 i = PROTECT(GET_SLOT(obj, iSym));
243 int *pp = INTEGER(p), *pi = INTEGER(i), j, k, kend;
244 pp++;
245
246 UNPROTECT(4); /* i, p, x, obj */
247
248 switch (class[0]) {
249 case 'n':
250 {
251 int_fast64_t s = (int_fast64_t) pp[n - 1] + ((un) ? n : 0);
252 if (sy) {
253 int up = (class[2] == 'C') == (ul == 'U');
254 s *= 2;
255 for (j = 0, k = 0; j < n; ++j) {
256 kend = pp[j];
257 if (k < kend && pi[(up) ? kend - 1 : k] == j)
258 --s;
259 k = kend;
260 }
261 }
262 if (s <= INT_MAX)
263 ans = Rf_ScalarInteger((int) s);
264 else
265 ans = Rf_ScalarReal((double) s);
266 break;
267 }
268 case 'l':
269 case 'i':
270 {
271 int *px = (class[0] == 'l') ? LOGICAL(x) : INTEGER(x);
272 int_fast64_t s = (un) ? n : 0, t = 0;
273 unsigned int count = 0;
274 for (j = 0, k = 0; j < n; ++j) {
275 kend = pp[j];
276 while (k < kend) {
277 if (px[k] != NA_INTEGER) {
278 unsigned int d = (sy && pi[k] != j) ? 2 : 1;
279 if (count > UINT_MAX - d)
280 TRY_INCREMENT(s, t, overC);
281 t += (sy && pi[k] != j) ? 2LL * px[k] : px[k];
282 count += d;
283 }
284 else if (!narm)
285 return Rf_ScalarInteger(NA_INTEGER);
286 ++k;
287 }
288 }
289 TRY_INCREMENT(s, t, overC);
290 if (s > INT_MIN && s <= INT_MAX)
291 ans = Rf_ScalarInteger((int) s);
292 else
293 ans = Rf_ScalarReal((double) s);
294 break;
295overC:
296 ;
297 long double lr = (long double) s + (long double) t;
298 for (; j < n; ++j) {
299 kend = pp[j];
300 while (k < kend) {
301 if (px[k] != NA_INTEGER)
302 lr += (sy && pi[k] != j) ? 2.0L * px[k] : px[k];
303 else if (!narm)
304 return Rf_ScalarInteger(NA_INTEGER);
305 ++k;
306 }
307 }
308 ans = Rf_ScalarReal(LONGDOUBLE_AS_DOUBLE(lr));
309 break;
310 }
311 case 'd':
312 {
313 double *px = REAL(x);
314 long double lr = (un) ? (long double) n : 0.0;
315 for (j = 0, k = 0; j < n; ++j) {
316 kend = pp[j];
317 while (k < kend) {
318 if (!(narm && ISNAN(px[k])))
319 lr += (sy && pi[k] != j) ? 2.0L * px[k] : px[k];
320 ++k;
321 }
322 }
323 ans = Rf_ScalarReal(LONGDOUBLE_AS_DOUBLE(lr));
324 break;
325 }
326 case 'z':
327 {
328 Rcomplex *px = COMPLEX(x), tmp;
329 long double lr = (un) ? (long double) n : 0.0;
330 long double li = 0.0;
331 for (j = 0, k = 0; j < n; ++j) {
332 kend = pp[j];
333 while (k < kend) {
334 if (!(narm && (ISNAN(px[k].r) || (!he && ISNAN(px[k].i))))) {
335 lr += (sy && pi[k] != j) ? 2.0L * px[k].r : px[k].r;
336 if (!he)
337 li += (sy && pi[k] != j) ? 2.0L * px[k].i : px[k].i;
338 }
339 ++k;
340 }
341 }
342 tmp.r = LONGDOUBLE_AS_DOUBLE(lr);
343 tmp.i = LONGDOUBLE_AS_DOUBLE(li);
344 ans = Rf_ScalarComplex(tmp);
345 break;
346 }
347 default:
348 break;
349 }
350
351 } else {
352
353 SEXP i = PROTECT(GET_SLOT(obj, Matrix_iSym)),
354 j = PROTECT(GET_SLOT(obj, Matrix_jSym));
355 int *pi = INTEGER(i), *pj = INTEGER(j);
356 R_xlen_t k, kend = XLENGTH(i);
357
358 UNPROTECT(4); /* j, i, x, obj */
359
360 switch (class[0]) {
361 case 'n':
362 {
363 int_fast64_t s = (int_fast64_t) XLENGTH(i) + ((un) ? n : 0);
364 if (sy) {
365 s *= 2;
366 for (k = 0; k < kend; ++k)
367 if (pi[k] == pj[k])
368 --s;
369 }
370 if (s <= INT_MAX)
371 ans = Rf_ScalarInteger((int) s);
372 else
373 ans = Rf_ScalarReal((double) s);
374 break;
375 }
376 case 'l':
377 case 'i':
378 {
379 int *px = (class[0] == 'l') ? LOGICAL(x) : INTEGER(x);
380 int_fast64_t s = (un) ? n : 0, t = 0;
381 unsigned int count = 0;
382 for (k = 0; k < kend; ++k)
383 if (px[k] != NA_INTEGER) {
384 unsigned int d = (sy && pi[k] != pj[k]) ? 2 : 1;
385 if (count > UINT_MAX - d)
386 TRY_INCREMENT(s, t, overT);
387 t += (sy && pi[k] != pj[k]) ? 2LL * px[k] : px[k];
388 count += d;
389 }
390 else if (!narm)
391 return Rf_ScalarInteger(NA_INTEGER);
392 TRY_INCREMENT(s, t, overT);
393 if (s > INT_MIN && s <= INT_MAX)
394 ans = Rf_ScalarInteger((int) s);
395 else
396 ans = Rf_ScalarReal((double) s);
397 break;
398overT:
399 ;
400 long double lr = (long double) s + (long double) t;
401 for (; k < kend; ++k) {
402 if (px[k] != NA_INTEGER)
403 lr += (sy && pi[k] != pj[k]) ? 2.0L * px[k] : px[k];
404 else if (!narm)
405 return Rf_ScalarInteger(NA_INTEGER);
406 }
407 ans = Rf_ScalarReal(LONGDOUBLE_AS_DOUBLE(lr));
408 break;
409 }
410 case 'd':
411 {
412 double *px = REAL(x);
413 long double lr = (un) ? (long double) n : 0.0;
414 for (k = 0; k < kend; ++k)
415 if (!(narm && ISNAN(px[k])))
416 lr += (sy && pi[k] != pj[k]) ? 2.0L * px[k] : px[k];
417 ans = Rf_ScalarReal(LONGDOUBLE_AS_DOUBLE(lr));
418 break;
419 }
420 case 'z':
421 {
422 Rcomplex *px = COMPLEX(x), tmp;
423 long double lr = (un) ? (long double) n : 0.0;
424 long double li = 0.0;
425 for (k = 0; k < kend; ++k)
426 if (!(narm && (ISNAN(px[k].r) || (!he && ISNAN(px[k].i))))) {
427 lr += (sy && pi[k] != pj[k]) ? 2.0L * px[k].r : px[k].r;
428 if (!he)
429 li += (sy && pi[k] != pj[k]) ? 2.0L * px[k].i : px[k].i;
430 }
431 tmp.r = LONGDOUBLE_AS_DOUBLE(lr);
432 tmp.i = LONGDOUBLE_AS_DOUBLE(li);
433 ans = Rf_ScalarComplex(tmp);
434 break;
435 }
436 default:
437 break;
438 }
439
440 }
441
442 return ans;
443}
444
445SEXP R_dense_sum(SEXP s_obj, SEXP s_narm)
446{
447 const char *class = Matrix_class(s_obj, valid_dense, 6, __func__);
448
449 int narm;
450 VALID_LOGIC2(s_narm, narm);
451
452 return dense_sum(s_obj, class, narm);
453}
454
455SEXP R_sparse_sum(SEXP s_obj, SEXP s_narm)
456{
457 const char *class = Matrix_class(s_obj, valid_sparse, 6, __func__);
458
459 int narm;
460 VALID_LOGIC2(s_narm, narm);
461
462 return sparse_sum(s_obj, class, narm);
463}
464
465SEXP dense_prod(SEXP obj, const char *class, int narm)
466{
467 SEXP dim = GET_SLOT(obj, Matrix_DimSym);
468 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1];
469
470 char ul = '\0', ct = '\0', nu = '\0';
471 if (class[1] != 'g')
472 ul = UPLO(obj);
473 if (class[1] == 's' && class[0] == 'z')
474 ct = TRANS(obj);
475 if (class[1] == 't')
476 nu = DIAG(obj);
477
478 SEXP x = GET_SLOT(obj, Matrix_xSym);
479 int i, j, packed = class[2] == 'p',
480 sy = class[1] == 's', he = sy && ct == 'C',
481 un = class[1] == 't' && nu != 'N';
482 long double lr = 1.0, li = 0.0;
483
484#define PROD \
485 do { \
486 if (class[1] == 'g') { \
487 for (j = 0; j < n; ++j) \
488 PROD_KERNEL(for (i = 0; i < m; ++i)); \
489 } else if (class[1] == 's') { \
490 if (ul == 'U') \
491 for (j = 0; j < n; ++j) { \
492 PROD_KERNEL(for (i = 0; i <= j; ++i)); \
493 if (!packed) \
494 px += n - j - 1; \
495 } \
496 else \
497 for (j = 0; j < n; ++j) { \
498 if (!packed) \
499 px += j; \
500 PROD_KERNEL(for (i = j; i < n; ++i)); \
501 } \
502 } else if (nu == 'N') { \
503 if (ul == 'U') \
504 for (j = 0; j < n; ++j) { \
505 if (j == 1) { lr *= 0.0L; li *= 0.0L; } \
506 PROD_KERNEL(for (i = 0; i <= j; ++i)); \
507 if (!packed) \
508 px += n - j - 1; \
509 } \
510 else \
511 for (j = 0; j < n; ++j) { \
512 if (j == 1) { lr *= 0.0L; li *= 0.0L; } \
513 if (!packed) \
514 px += j; \
515 PROD_KERNEL(for (i = j; i < n; ++i)); \
516 } \
517 } else { \
518 if (ul == 'U') \
519 for (j = 0; j < n; ++j) { \
520 if (j == 1) { lr *= 0.0L; li *= 0.0L; } \
521 PROD_KERNEL(for (i = 0; i < j; ++i)); \
522 px += 1; \
523 if (!packed) \
524 px += n - j - 1; \
525 } \
526 else \
527 for (j = 0; j < n; ++j) { \
528 if (j == 1) { lr *= 0.0L; li *= 0.0L; } \
529 if (!packed) \
530 px += j; \
531 px += 1; \
532 PROD_KERNEL(for (i = j + 1; i < n; ++i)); \
533 } \
534 } \
535 } while (0)
536
537 switch (class[0]) {
538 case 'n':
539 {
540 int *px = LOGICAL(x);
541 if (class[1] == 't') {
542 if (n > 1 || (n == 1 && !un && *px == 0))
543 return Rf_ScalarReal(0.0);
544 break;
545 }
546
547#define PROD_KERNEL(__for__) \
548 do { \
549 __for__ { \
550 if (*px == 0) \
551 return Rf_ScalarReal(0.0); \
552 px += 1; \
553 } \
554 } while (0)
555
556 PROD;
557
558#undef PROD_KERNEL
559
560 break;
561 }
562 case 'l':
563 case 'i':
564 {
565 int *px = (class[0] == 'l') ? LOGICAL(x) : INTEGER(x);
566
567#define PROD_KERNEL(__for__) \
568 do { \
569 __for__ { \
570 if (*px != NA_INTEGER) \
571 lr *= (sy && i != j) ? (long double) *px * *px : *px; \
572 else if (!narm) \
573 lr *= NA_REAL; \
574 px += 1; \
575 } \
576 } while (0)
577
578 PROD;
579
580#undef PROD_KERNEL
581
582 break;
583 }
584 case 'd':
585 {
586 double *px = REAL(x);
587
588#define PROD_KERNEL(__for__) \
589 do { \
590 __for__ { \
591 if (!(narm && ISNAN(*px))) \
592 lr *= (sy && i != j) ? (long double) *px * *px : *px; \
593 px += 1; \
594 } \
595 } while (0)
596
597 PROD;
598
599#undef PROD_KERNEL
600
601 break;
602 }
603 case 'z':
604 {
605 Rcomplex *px = COMPLEX(x);
606 long double lr0, li0;
607
608#define PROD_KERNEL(__for__) \
609 do { \
610 __for__ { \
611 if (!(narm && (ISNAN((*px).r) || (!(he && i == j) && ISNAN((*px).i))))) { \
612 if (he) { \
613 lr0 = (*px).r; \
614 if (i != j) { \
615 li0 = (*px).i; \
616 lr0 = lr0 * lr0 + li0 * li0; \
617 } \
618 lr *= lr0; \
619 li *= lr0; \
620 } else { \
621 lr0 = lr; li0 = li; \
622 lr = lr0 * (*px).r - li0 * (*px).i; \
623 li = li0 * (*px).r + lr0 * (*px).i; \
624 if (i != j) { \
625 lr0 = lr; li0 = li; \
626 lr = lr0 * (*px).r - li0 * (*px).i; \
627 li = li0 * (*px).r + lr0 * (*px).i; \
628 } \
629 } \
630 } \
631 px += 1; \
632 } \
633 } while (0)
634
635 PROD;
636
637#undef PROD_KERNEL
638
639 break;
640 }
641 default:
642 break;
643 }
644
645#undef PROD
646
647 SEXP ans;
648 if (class[0] != 'z')
649 ans = Rf_ScalarReal(LONGDOUBLE_AS_DOUBLE(lr));
650 else {
651 Rcomplex tmp;
652 tmp.r = LONGDOUBLE_AS_DOUBLE(lr);
653 tmp.i = LONGDOUBLE_AS_DOUBLE(li);
654 ans = Rf_ScalarComplex(tmp);
655 }
656 return ans;
657}
658
659SEXP sparse_prod(SEXP obj, const char *class, int narm)
660{
661 PROTECT(obj = sparse_aggregate(obj, class));
662
663 int *pdim = DIM(obj), m = pdim[0], n = pdim[1];
664
665 char ul = '\0', ct = '\0', nu = '\0';
666 if (class[1] != 'g')
667 ul = UPLO(obj);
668 if (class[1] == 's' && class[0] == 'z')
669 ct = TRANS(obj);
670 if (class[1] == 't')
671 nu = DIAG(obj);
672
673 SEXP x = PROTECT((class[0] == 'n') ? R_NilValue : GET_SLOT(obj, Matrix_xSym));
674 int sy = class[1] == 's', he = sy && ct == 'C',
675 un = class[1] == 't' && nu != 'N';
676 long double lr = 1.0, li = 0.0;
677
678 int_fast64_t mn = (int_fast64_t) m * n,
679 q = (sy) ? (mn + n) / 2 : ((un) ? mn - n : mn);
680
681 if (class[2] != 'T') {
682
683 if (class[2] == 'R')
684 SWAP(m, n, int, );
685
686 SEXP iSym = (class[2] == 'C') ? Matrix_iSym : Matrix_jSym,
687 p = PROTECT(GET_SLOT(obj, Matrix_pSym)),
688 i = PROTECT(GET_SLOT(obj, iSym));
689 int *pp = INTEGER(p), *pi = INTEGER(i), j, k, kend;
690 pp++;
691
692 UNPROTECT(4); /* i, p, x, obj */
693
694 int seen0 = 0, i0, i1,
695 up = !sy || (class[2] == 'C') == (ul == 'U'),
696 lo = !sy || (class[2] == 'C') == (ul != 'U');
697
698 switch (class[0]) {
699 case 'n':
700 if (pp[n - 1] < q)
701 lr = 0.0;
702 break;
703 case 'l':
704 case 'i':
705 {
706 int *px = (class[0] == 'l') ? LOGICAL(x) : INTEGER(x);
707 for (j = 0, k = 0; j < n; ++j) {
708 kend = pp[j];
709 if (!seen0) {
710 i0 = (up) ? 0 : j;
711 i1 = (lo) ? m : j + 1;
712 }
713 while (k < kend) {
714 if (!seen0) {
715 if (pi[k] == i0 || (un && i0 == j && pi[k] == ++i0))
716 ++i0;
717 else {
718 lr *= 0.0;
719 seen0 = 1;
720 }
721 }
722 if (px[k] != NA_INTEGER)
723 lr *= (sy && pi[k] != j)
724 ? (long double) px[k] * px[k] : px[k];
725 else if (!narm)
726 lr *= NA_REAL;
727 ++k;
728 }
729 if (!(seen0 || i1 == i0 || (un && i0 == j && i1 == ++i0))) {
730 lr *= 0.0;
731 seen0 = 1;
732 }
733 }
734 break;
735 }
736 case 'd':
737 {
738 double *px = REAL(x);
739 for (j = 0, k = 0; j < n; ++j) {
740 kend = pp[j];
741 if (!seen0) {
742 i0 = (up) ? 0 : j;
743 i1 = (lo) ? m : j + 1;
744 }
745 while (k < kend) {
746 if (!seen0) {
747 if (pi[k] == i0 || (un && i0 == j && pi[k] == ++i0))
748 ++i0;
749 else {
750 lr *= 0.0;
751 seen0 = 1;
752 }
753 }
754 if (!(narm && ISNAN(px[k])))
755 lr *= (sy && pi[k] != j)
756 ? (long double) px[k] * px[k] : px[k];
757 ++k;
758 }
759 if (!(seen0 || i1 == i0 || (un && i0 == j && i1 == ++i0))) {
760 lr *= 0.0;
761 seen0 = 1;
762 }
763 }
764 break;
765 }
766 case 'z':
767 {
768 Rcomplex *px = COMPLEX(x);
769 long double lr0, li0;
770 for (j = 0, k = 0; j < n; ++j) {
771 kend = pp[j];
772 if (!seen0) {
773 i0 = (up) ? 0 : j;
774 i1 = (lo) ? m : j + 1;
775 }
776 while (k < kend) {
777 if (!seen0) {
778 if (pi[k] == i0 || (un && i0 == j && pi[k] == ++i0))
779 ++i0;
780 else {
781 lr *= 0.0;
782 li *= 0.0;
783 seen0 = 1;
784 }
785 }
786 if (!(narm && (ISNAN(px[k].r) || (!(he && pi[k] == j) && ISNAN(px[k].i))))) {
787 if (he) {
788 lr0 = px[k].r;
789 if (pi[k] != j) {
790 li0 = px[k].i;
791 lr0 = lr0 * lr0 + li0 * li0;
792 }
793 lr *= lr0;
794 li *= lr0;
795 } else {
796 lr0 = lr; li0 = li;
797 lr = lr0 * px[k].r - li0 * px[k].i;
798 li = li0 * px[k].r + lr0 * px[k].i;
799 if (pi[k] != j) {
800 lr0 = lr; li0 = li;
801 lr = lr0 * px[k].r - li0 * px[k].i;
802 li = li0 * px[k].r + lr0 * px[k].i;
803 }
804 }
805 }
806 ++k;
807 }
808 if (!(seen0 || i1 == i0 || (un && i0 == j && i1 == ++i0))) {
809 lr *= 0.0;
810 li *= 0.0;
811 seen0 = 1;
812 }
813 }
814 break;
815 }
816 default:
817 break;
818 }
819
820 } else {
821
822 SEXP i = PROTECT(GET_SLOT(obj, Matrix_iSym)),
823 j = PROTECT(GET_SLOT(obj, Matrix_jSym));
824 int *pi = INTEGER(i), *pj = INTEGER(j);
825 R_xlen_t k, kend = XLENGTH(i);
826 if (XLENGTH(i) < q)
827 lr = 0.0;
828
829 UNPROTECT(4); /* j, i, x, obj */
830
831 switch (class[0]) {
832 case 'l':
833 case 'i':
834 {
835 int *px = (class[0] == 'l') ? LOGICAL(x) : INTEGER(x);
836 for (k = 0; k < kend; ++k)
837 if (px[k] != NA_INTEGER)
838 lr *= (sy && pi[k] != pj[k])
839 ? (long double) px[k] * px[k] : px[k];
840 else if (!narm)
841 lr *= NA_REAL;
842 break;
843 }
844 case 'd':
845 {
846 double *px = REAL(x);
847 for (k = 0; k < kend; ++k)
848 if (!(narm && ISNAN(px[k])))
849 lr *= (sy && pi[k] != pj[k])
850 ? (long double) px[k] * px[k] : px[k];
851 break;
852 }
853 case 'z':
854 {
855 Rcomplex *px = COMPLEX(x);
856 long double lr0, li0;
857 for (k = 0; k < kend; ++k)
858 if (!(narm && (ISNAN(px[k].r) || (!(he && pi[k] == pj[k]) && ISNAN(px[k].i))))) {
859 if (he) {
860 lr0 = px[k].r;
861 if (pi[k] != pj[k]) {
862 li0 = px[k].i;
863 lr0 = lr0 * lr0 + li0 * li0;
864 }
865 lr *= lr0;
866 li *= lr0;
867 } else {
868 lr0 = lr; li0 = li;
869 lr = lr0 * px[k].r - li0 * px[k].i;
870 li = li0 * px[k].r + lr0 * px[k].i;
871 if (pi[k] != pj[k]) {
872 lr0 = lr; li0 = li;
873 lr = lr0 * px[k].r - li0 * px[k].i;
874 li = li0 * px[k].r + lr0 * px[k].i;
875 }
876 }
877 }
878 break;
879 }
880 default:
881 break;
882 }
883
884 }
885
886 SEXP ans;
887 if (class[0] != 'z')
888 ans = Rf_ScalarReal(LONGDOUBLE_AS_DOUBLE(lr));
889 else {
890 Rcomplex tmp;
891 tmp.r = LONGDOUBLE_AS_DOUBLE(lr);
892 tmp.i = LONGDOUBLE_AS_DOUBLE(li);
893 ans = Rf_ScalarComplex(tmp);
894 }
895 return ans;
896}
897
898SEXP R_dense_prod(SEXP s_obj, SEXP s_narm)
899{
900 const char *class = Matrix_class(s_obj, valid_dense, 6, __func__);
901
902 int narm;
903 VALID_LOGIC2(s_narm, narm);
904
905 return dense_prod(s_obj, class, narm);
906}
907
908SEXP R_sparse_prod(SEXP s_obj, SEXP s_narm)
909{
910 const char *class = Matrix_class(s_obj, valid_sparse, 6, __func__);
911
912 int narm;
913 VALID_LOGIC2(s_narm, narm);
914
915 return sparse_prod(s_obj, class, narm);
916}
SWAP
#define SWAP(a, b, t, op)
Definition Mdefines.h:138
valid_dense
const char * valid_dense[]
Definition objects.c:3
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
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_sparse
const char * valid_sparse[]
Definition Mdefines.h:328
VALID_LOGIC2
#define VALID_LOGIC2(s, d)
Definition Mdefines.h:216
SUM
#define SUM
sparse_sum
SEXP sparse_sum(SEXP obj, const char *class, int narm)
Definition Summary.c:216
R_sparse_prod
SEXP R_sparse_prod(SEXP s_obj, SEXP s_narm)
Definition Summary.c:908
R_dense_sum
SEXP R_dense_sum(SEXP s_obj, SEXP s_narm)
Definition Summary.c:445
LONGDOUBLE_AS_DOUBLE
#define LONGDOUBLE_AS_DOUBLE(x)
Definition Summary.c:9
dense_sum
SEXP dense_sum(SEXP obj, const char *class, int narm)
Definition Summary.c:24
PROD
#define PROD
sparse_aggregate
SEXP sparse_aggregate(SEXP, const char *)
Definition aggregate.c:5
R_sparse_sum
SEXP R_sparse_sum(SEXP s_obj, SEXP s_narm)
Definition Summary.c:455
R_dense_prod
SEXP R_dense_prod(SEXP s_obj, SEXP s_narm)
Definition Summary.c:898
TRY_INCREMENT
#define TRY_INCREMENT(x, y, label)
Definition Summary.c:12
dense_prod
SEXP dense_prod(SEXP obj, const char *class, int narm)
Definition Summary.c:465
sparse_prod
SEXP sparse_prod(SEXP obj, const char *class, int narm)
Definition Summary.c:659
Matrix_DimSym
SEXP Matrix_DimSym
Definition init.c:598
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