Matrix r4655
Loading...
Searching...
No Matches
sparse.c
Go to the documentation of this file.
1#include <math.h> /* fabs, hypot */
2#include "Mdefines.h"
3#include "sparse.h"
4
5SEXP sparse_drop0(SEXP from, const char *class, double tol)
6{
7 if (class[0] == 'n')
8 return from;
9
10 SEXP to, x0 = PROTECT(GET_SLOT(from, Matrix_xSym));
11
12#define TOLBASED_ISNZ_REAL(_X_) \
13 (ISNA_REAL(_X_) || fabs(_X_) > tol)
14
15#define TOLBASED_ISNZ_COMPLEX(_X_) \
16 (ISNA_COMPLEX(_X_) || hypot((_X_).r, (_X_).i) > tol)
17
18#define DROP0_CASES(_DO_) \
19 do { \
20 switch (class[0]) { \
21 case 'l': \
22 _DO_(int, LOGICAL, ISNZ_LOGICAL); \
23 break; \
24 case 'i': \
25 _DO_(int, INTEGER, ISNZ_INTEGER); \
26 break; \
27 case 'd': \
28 if (tol > 0.0) \
29 _DO_(double, REAL, TOLBASED_ISNZ_REAL); \
30 else \
31 _DO_(double, REAL, ISNZ_REAL); \
32 break; \
33 case 'z': \
34 if (tol > 0.0) \
35 _DO_(Rcomplex, COMPLEX, TOLBASED_ISNZ_COMPLEX); \
36 else \
37 _DO_(Rcomplex, COMPLEX, ISNZ_COMPLEX); \
38 break; \
39 default: \
40 break; \
41 } \
42 } while (0)
43
44 if (class[2] != 'T') {
45
46 SEXP p0 = PROTECT(GET_SLOT(from, Matrix_pSym));
47 int *pp0 = INTEGER(p0), k, n = (int) (XLENGTH(p0) - 1),
48 nnz0 = pp0[n], nnz1 = 0;
49
50#undef DROP0_LOOP1
51#define DROP0_LOOP1(_CTYPE_, _PTR_, _ISNZ_) \
52 do { \
53 _CTYPE_ *px0 = _PTR_(x0); \
54 for (k = 0; k < nnz0; ++k) { \
55 if (_ISNZ_(*px0)) \
56 ++nnz1; \
57 ++px0; \
58 } \
59 } while (0)
60
61 DROP0_CASES(DROP0_LOOP1);
62 if (nnz1 == nnz0) {
63 UNPROTECT(2); /* p0, x0 */
64 return from;
65 }
66 PROTECT(to = newObject(class));
67
68 SEXP iSym = (class[2] == 'C') ? Matrix_iSym : Matrix_jSym,
69 i0 = PROTECT(GET_SLOT(from, iSym)),
70 p1 = PROTECT(allocVector(INTSXP, (R_xlen_t) n + 1)),
71 i1 = PROTECT(allocVector(INTSXP, nnz1)),
72 x1 = PROTECT(allocVector(TYPEOF(x0), nnz1));
73 int *pi0 = INTEGER(i0), *pp1 = INTEGER(p1), *pi1 = INTEGER(i1),
74 j, kend;
75 pp0++; *(pp1++) = 0;
76 SET_SLOT(to, Matrix_pSym, p1);
77 SET_SLOT(to, iSym, i1);
78 SET_SLOT(to, Matrix_xSym, x1);
79
80#undef DROP0_LOOP2
81#define DROP0_LOOP2(_CTYPE_, _PTR_, _ISNZ_) \
82 do { \
83 _CTYPE_ *px0 = _PTR_(x0), *px1 = _PTR_(x1); \
84 for (j = 0, k = 0; j < n; ++j) { \
85 pp1[j] = pp1[j - 1]; \
86 kend = pp0[j]; \
87 while (k < kend) { \
88 if (_ISNZ_(*px0)) { \
89 ++pp1[j]; \
90 *(pi1++) = *pi0; \
91 *(px1++) = *px0; \
92 } \
93 ++k; ++pi0; ++px0; \
94 } \
95 } \
96 } while (0)
97
98 DROP0_CASES(DROP0_LOOP2);
99 UNPROTECT(7); /* x1, i1, p1, i0, to, p0, x0 */
100
101 } else {
102
103 R_xlen_t k, nnz0 = XLENGTH(x0), nnz1 = 0;
104
105 DROP0_CASES(DROP0_LOOP1);
106 if (nnz1 == nnz0) {
107 UNPROTECT(1); /* x0 */
108 return from;
109 }
110 PROTECT(to = newObject(class));
111
112 SEXP i0 = PROTECT(GET_SLOT(from, Matrix_iSym)),
113 j0 = PROTECT(GET_SLOT(from, Matrix_jSym)),
114 i1 = PROTECT(allocVector(INTSXP, nnz1)),
115 j1 = PROTECT(allocVector(INTSXP, nnz1)),
116 x1 = PROTECT(allocVector(TYPEOF(x0), nnz1));
117 int *pi0 = INTEGER(i0), *pj0 = INTEGER(j0),
118 *pi1 = INTEGER(i1), *pj1 = INTEGER(j1);
119 SET_SLOT(to, Matrix_iSym, i1);
120 SET_SLOT(to, Matrix_jSym, j1);
121 SET_SLOT(to, Matrix_xSym, x1);
122
123#undef DROP0_LOOP2
124#define DROP0_LOOP2(_CTYPE_, _PTR_, _ISNZ_) \
125 do { \
126 _CTYPE_ *px0 = _PTR_(x0), *px1 = _PTR_(x1); \
127 for (k = 0; k < nnz0; ++k) { \
128 if (_ISNZ_(*px0)) { \
129 *(pi1++) = *pi0; \
130 *(pj1++) = *pj0; \
131 *(px1++) = *px0; \
132 } \
133 ++pi0; ++pj0; ++px0; \
134 } \
135 } while (0)
136
137 DROP0_CASES(DROP0_LOOP2);
138 UNPROTECT(7); /* x1, j1, i1, j0, i0, to, x0 */
139
140 }
141
142 PROTECT(to);
143
144 SEXP dim = PROTECT(GET_SLOT(from, Matrix_DimSym));
145 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1];
146 if (m != n || n > 0)
147 SET_SLOT(to, Matrix_DimSym, dim);
148 UNPROTECT(1); /* dim */
149
150 SEXP dimnames = PROTECT(GET_SLOT(from, Matrix_DimNamesSym));
151 SET_SLOT(to, Matrix_DimNamesSym, dimnames);
152 UNPROTECT(1); /* dimnames */
153
154 if (class[1] != 'g') {
155 SEXP uplo = PROTECT(GET_SLOT(from, Matrix_uploSym));
156 char ul = *CHAR(STRING_ELT(uplo, 0));
157 if (ul != 'U')
158 SET_SLOT(to, Matrix_uploSym, uplo);
159 UNPROTECT(1); /* uplo */
160 }
161 if (class[1] == 't') {
162 SEXP diag = PROTECT(GET_SLOT(from, Matrix_diagSym));
163 char di = *CHAR(STRING_ELT(diag, 0));
164 if (di != 'N')
165 SET_SLOT(to, Matrix_diagSym, diag);
166 UNPROTECT(1); /* diag */
167 } else {
168 SEXP factors = PROTECT(GET_SLOT(from, Matrix_factorsSym));
169 if (LENGTH(factors) > 0)
170 SET_SLOT(to, Matrix_factorsSym, factors);
171 UNPROTECT(1); /* factors */
172 }
173
174#undef TOLBASED_ISNZ_REAL
175#undef TOLBASED_ISNZ_COMPLEX
176#undef DROP0_CASES
177#undef DROP0_LOOP1
178#undef DROP0_LOOP2
179
180 UNPROTECT(1); /* to */
181 return to;
182}
183
184/* drop0(<[CRT]sparseMatrix>, tol) */
185SEXP R_sparse_drop0(SEXP from, SEXP tol)
186{
187 static const char *valid[] = {
188 VALID_CSPARSE, VALID_RSPARSE, VALID_TSPARSE, "" };
189 int ivalid = R_check_class_etc(from, valid);
190 if (ivalid < 0)
191 ERROR_INVALID_CLASS(from, __func__);
192
193 double tol_;
194 if (TYPEOF(tol) != REALSXP || LENGTH(tol) < 1 ||
195 ISNAN(tol_ = REAL(tol)[0]))
196 error(_("'%s' is not a number"), "tol");
197
198 return sparse_drop0(from, valid[ivalid], tol_);
199}
200
201SEXP sparse_diag_U2N(SEXP from, const char *class)
202{
203 if (class[1] != 't')
204 return from;
205
206 SEXP diag = PROTECT(GET_SLOT(from, Matrix_diagSym));
207 char di = *CHAR(STRING_ELT(diag, 0));
208 UNPROTECT(1); /* diag */
209 if (di == 'N')
210 return from;
211
212 SEXP val = PROTECT(ScalarLogical(1));
213 from = R_sparse_diag_set(from, val);
214 UNPROTECT(1); /* val */
215
216 return from;
217}
218
219/* diagU2N(<[CRT]sparseMatrix>), parallel to R-level ..diagU2N(),
220 though that is more general, working for _all_ Matrix
221*/
222SEXP R_sparse_diag_U2N(SEXP from)
223{
224 static const char *valid[] = {
225 VALID_CSPARSE, VALID_RSPARSE, VALID_TSPARSE, "" };
226 int ivalid = R_check_class_etc(from, valid);
227 if (ivalid < 0)
228 ERROR_INVALID_CLASS(from, __func__);
229
230 return sparse_diag_U2N(from, valid[ivalid]);
231}
232
233SEXP sparse_diag_N2U(SEXP from, const char *class)
234{
235 if (class[1] != 't')
236 return from;
237
238 SEXP diag = PROTECT(GET_SLOT(from, Matrix_diagSym));
239 char di = *CHAR(STRING_ELT(diag, 0));
240 UNPROTECT(1); /* diag */
241 if (di != 'N')
242 return from;
243
244 SEXP dim = PROTECT(GET_SLOT(from, Matrix_DimSym));
245 int n = INTEGER(dim)[0];
246 UNPROTECT(1); /* dim */
247
248 if (n == 0)
249 PROTECT(from = duplicate(from));
250 else {
251 SEXP uplo = PROTECT(GET_SLOT(from, Matrix_uploSym));
252 char ul = *CHAR(STRING_ELT(uplo, 0));
253 UNPROTECT(1); /* uplo */
254 if (ul == 'U')
255 PROTECT(from = sparse_band(from, class, 1, n - 1));
256 else
257 PROTECT(from = sparse_band(from, class, 1 - n, -1));
258 }
259
260 PROTECT(diag = mkString("U"));
261 SET_SLOT(from, Matrix_diagSym, diag);
262 UNPROTECT(2); /* diag, from */
263
264 return from;
265}
266
267/* diagN2U(<[CRT]sparseMatrix>), parallel to R-level ..diagN2U(),
268 though that is more general, working for _all_ Matrix
269*/
270SEXP R_sparse_diag_N2U(SEXP from)
271{
272 static const char *valid[] = {
273 VALID_CSPARSE, VALID_RSPARSE, VALID_TSPARSE, "" };
274 int ivalid = R_check_class_etc(from, valid);
275 if (ivalid < 0)
276 ERROR_INVALID_CLASS(from, __func__);
277
278 return sparse_diag_N2U(from, valid[ivalid]);
279}
280
281SEXP sparse_band(SEXP from, const char *class, int a, int b)
282{
283 SEXP dim = PROTECT(GET_SLOT(from, Matrix_DimSym));
284 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1];
285 UNPROTECT(1); /* dim */
286
287 /* Need tri[ul](<0-by-0>) and tri[ul](<1-by-1>) to be triangularMatrix */
288 if (a <= 1 - m && b >= n - 1 &&
289 (class[1] == 't' || m != n || m > 1 || n > 1))
290 return from;
291
292 int ge = 0, sy = 0, tr = 0;
293 ge = m != n || !((tr = a >= 0 || b <= 0 || class[1] == 't') ||
294 (sy = a == -b && class[1] == 's'));
295
296 char ul0 = 'U', ul1 = 'U', di = 'N';
297 if (class[1] != 'g') {
298 SEXP uplo = PROTECT(GET_SLOT(from, Matrix_uploSym));
299 ul0 = *CHAR(STRING_ELT(uplo, 0));
300 UNPROTECT(1); /* uplo */
301
302 if (class[1] == 't') {
303 /* Be fast if band contains entire triangle */
304 if ((ul0 == 'U') ? (a <= 0 && b >= n - 1) : (b >= 0 && a <= 1 - m))
305 return from;
306 else if (a <= 0 && b >= 0) {
307 SEXP diag = PROTECT(GET_SLOT(from, Matrix_diagSym));
308 di = *CHAR(STRING_ELT(diag, 0));
309 UNPROTECT(1); /* diag */
310 }
311 }
312 }
313
314 /* band(<R>, a, b) is equivalent to t(band(t(<R>), -b, -a)) ! */
315
316 if (class[2] == 'R') {
317 int r;
318 r = m; m = n; n = r;
319 r = a; a = -b; b = -r;
320 ul0 = (ul0 == 'U') ? 'L' : 'U';
321 from = sparse_transpose(from, class, 1);
322 }
323 PROTECT(from);
324
325 char cl[] = "...Matrix";
326 cl[0] = class[0];
327 cl[1] = (ge) ? 'g' : ((tr) ? 't' : 's');
328 cl[2] = (class[2] == 'R') ? 'C' : class[2];
329 SEXP to = PROTECT(newObject(cl));
330
331 dim = GET_SLOT(to, Matrix_DimSym);
332 pdim = INTEGER(dim);
333 pdim[0] = m;
334 pdim[1] = n;
335
336 SEXP dimnames = PROTECT(GET_SLOT(from, Matrix_DimNamesSym));
337 if (class[1] != 's' || sy)
338 SET_SLOT(to, Matrix_DimNamesSym, dimnames);
339 else
340 set_symmetrized_DimNames(to, dimnames, -1);
341 UNPROTECT(1); /* dimnames */
342
343 if (!ge) {
344 ul1 = (tr && class[1] != 't') ? ((a >= 0) ? 'U' : 'L') : ul0;
345 if (ul1 != 'U') {
346 SEXP uplo = PROTECT(mkString("L"));
347 SET_SLOT(to, Matrix_uploSym, uplo);
348 UNPROTECT(1); /* uplo */
349 }
350 if (di != 'N') {
351 SEXP diag = PROTECT(mkString("U"));
352 SET_SLOT(to, Matrix_diagSym, diag);
353 UNPROTECT(1); /* diag */
354 }
355 }
356
357 /* It remains to set some subset of 'p', 'i', 'j', 'x' ... */
358
359#define BAND_CASES \
360 do { \
361 switch (class[0]) { \
362 case 'l': \
363 BAND_SUBCASES(int, LOGICAL, SHOW); \
364 break; \
365 case 'i': \
366 BAND_SUBCASES(int, INTEGER, SHOW); \
367 break; \
368 case 'd': \
369 BAND_SUBCASES(double, REAL, SHOW); \
370 break; \
371 case 'z': \
372 BAND_SUBCASES(Rcomplex, COMPLEX, SHOW); \
373 break; \
374 default: \
375 break; \
376 } \
377 } while (0)
378
379 if (class[2] != 'T') {
380 SEXP p0 = PROTECT(GET_SLOT(from, Matrix_pSym)),
381 i0 = PROTECT(GET_SLOT(from, Matrix_iSym)),
382 p1 = PROTECT(allocVector(INTSXP, (R_xlen_t) n + 1));
383 int *pp0 = INTEGER(p0), *pi0 = INTEGER(i0), *pp1 = INTEGER(p1),
384 d, j, k, kend, nnz0 = pp0[n], nnz1 = 0;
385 pp0++; *(pp1++) = 0;
386 SET_SLOT(to, Matrix_pSym, p1);
387
388 if (class[1] == 's' && !sy) {
389 Matrix_memset(pp1, 0, n, sizeof(int));
390 for (j = 0, k = 0; j < n; ++j) {
391 kend = pp0[j];
392 while (k < kend) {
393 if ((d = j - pi0[k]) >= a && d <= b)
394 ++pp1[j];
395 if (d != 0 && -d >= a && -d <= b)
396 ++pp1[pi0[k]];
397 ++k;
398 }
399 }
400 for (j = 0; j < n; ++j) {
401 nnz1 += pp1[j];
402 pp1[j] = nnz1;
403 }
404 } else {
405 for (j = 0, k = 0; j < n; ++j) {
406 kend = pp0[j];
407 while (k < kend) {
408 if ((d = j - pi0[k]) >= a && d <= b)
409 ++nnz1;
410 ++k;
411 }
412 pp1[j] = nnz1;
413 }
414 }
415
416 if (nnz1 == nnz0 && (class[1] != 's' || sy)) {
417 /* No need to allocate in this case: band has all nonzero elements */
418 SET_SLOT(to, Matrix_iSym, i0);
419 if (class[0] != 'n') {
420 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym));
421 SET_SLOT(to, Matrix_xSym, x0);
422 UNPROTECT(1); /* x0 */
423 }
424 if (class[2] == 'R')
425 to = sparse_transpose(to, cl, 1);
426 UNPROTECT(5); /* p1, i0, p0, to, from */
427 return to;
428 }
429
430 SEXP i1 = PROTECT(allocVector(INTSXP, nnz1));
431 int *pi1 = INTEGER(i1);
432 SET_SLOT(to, Matrix_iSym, i1);
433
434#undef BAND_SUBCASES
435#define BAND_SUBCASES(_CTYPE_, _PTR_, _MASK_) \
436 do { \
437 _MASK_(_CTYPE_ *px0 = _PTR_(x0), *px1 = _PTR_(x1)); \
438 if (class[1] == 's' && !sy) { \
439 int *pp1_; \
440 Matrix_Calloc(pp1_, n, int); \
441 Matrix_memcpy(pp1_, pp1 - 1, n, sizeof(int)); \
442 for (j = 0, k = 0; j < n; ++j) { \
443 kend = pp0[j]; \
444 while (k < kend) { \
445 if ((d = j - pi0[k]) >= a && d <= b) { \
446 pi1[pp1_[j]] = pi0[k]; \
447 _MASK_(px1[pp1_[j]] = px0[k]); \
448 ++pp1_[j]; \
449 } \
450 if (d != 0 && -d >= a && -d <= b) { \
451 pi1[pp1_[pi0[k]]] = j; \
452 _MASK_(px1[pp1_[pi0[k]]] = px0[k]); \
453 ++pp1_[pi0[k]]; \
454 } \
455 ++k; \
456 } \
457 } \
458 Matrix_Free(pp1_, n); \
459 } else { \
460 for (j = 0, k = 0; j < n; ++j) { \
461 kend = pp0[j]; \
462 while (k < kend) { \
463 if ((d = j - pi0[k]) >= a && d <= b) { \
464 *(pi1++) = pi0[k]; \
465 _MASK_(*(px1++) = px0[k]); \
466 } \
467 ++k; \
468 } \
469 } \
470 } \
471 } while (0)
472
473 if (class[0] == 'n')
474 BAND_SUBCASES(int, LOGICAL, HIDE);
475 else {
476 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym)),
477 x1 = PROTECT(allocVector(TYPEOF(x0), nnz1));
478 SET_SLOT(to, Matrix_xSym, x1);
479 BAND_CASES;
480 UNPROTECT(2); /* x1, x0 */
481 }
482 if (class[2] == 'R')
483 to = sparse_transpose(to, cl, 1);
484
485 } else {
486
487 SEXP i0 = PROTECT(GET_SLOT(from, Matrix_iSym)),
488 j0 = PROTECT(GET_SLOT(from, Matrix_jSym));
489 int *pi0 = INTEGER(i0), *pj0 = INTEGER(j0), d;
490 R_xlen_t k, nnz0 = XLENGTH(i0), nnz1 = 0;
491
492 if (class[1] == 's' && !sy) {
493 for (k = 0; k < nnz0; ++k) {
494 if ((d = pj0[k] - pi0[k]) >= a && d <= b)
495 ++nnz1;
496 if (d != 0 && -d >= a && -d <= b)
497 ++nnz1;
498 }
499 } else {
500 for (k = 0; k < nnz0; ++k) {
501 if ((d = pj0[k] - pi0[k]) >= a && d <= b)
502 ++nnz1;
503 }
504 }
505
506 if (nnz1 == nnz0 && (class[1] != 's' || sy)) {
507 /* No need to allocate in this case: band has all nonzero elements */
508 SET_SLOT(to, Matrix_iSym, i0);
509 SET_SLOT(to, Matrix_jSym, j0);
510 if (class[0] != 'n') {
511 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym));
512 SET_SLOT(to, Matrix_xSym, x0);
513 UNPROTECT(1); /* x0 */
514 }
515 UNPROTECT(4); /* j0, i0, to, from */
516 return to;
517 }
518
519 SEXP i1 = PROTECT(allocVector(INTSXP, nnz1)),
520 j1 = PROTECT(allocVector(INTSXP, nnz1));
521 int *pi1 = INTEGER(i1), *pj1 = INTEGER(j1);
522 SET_SLOT(to, Matrix_iSym, i1);
523 SET_SLOT(to, Matrix_jSym, j1);
524
525#undef BAND_SUBCASES
526#define BAND_SUBCASES(_CTYPE_, _PTR_, _MASK_) \
527 do { \
528 _MASK_(_CTYPE_ *px0 = _PTR_(x0), *px1 = _PTR_(x1)); \
529 if (class[1] == 's' && !sy) { \
530 for (k = 0; k < nnz0; ++k) { \
531 if ((d = pj0[k] - pi0[k]) >= a && d <= b) { \
532 *(pi1++) = pi0[k]; \
533 *(pj1++) = pj0[k]; \
534 _MASK_(*(px1++) = px0[k]); \
535 } \
536 if (d != 0 && -d >= a && -d <= b) { \
537 *(pi1++) = pj0[k]; \
538 *(pj1++) = pi0[k]; \
539 _MASK_(*(px1++) = px0[k]); \
540 } \
541 } \
542 } else { \
543 for (k = 0; k < nnz0; ++k) { \
544 if ((d = pj0[k] - pi0[k]) >= a && d <= b) { \
545 *(pi1++) = pi0[k]; \
546 *(pj1++) = pj0[k]; \
547 _MASK_(*(px1++) = px0[k]); \
548 } \
549 } \
550 } \
551 } while (0)
552
553 if (class[0] == 'n')
554 BAND_SUBCASES(int, LOGICAL, HIDE);
555 else {
556 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym)),
557 x1 = PROTECT(allocVector(TYPEOF(x0), nnz1));
558 SET_SLOT(to, Matrix_xSym, x1);
559 BAND_CASES;
560 UNPROTECT(2); /* x1, x0 */
561 }
562
563 }
564
565#undef BAND_CASES
566#undef BAND_SUBCASES
567
568 UNPROTECT(6);
569 return to;
570}
571
572/* band(<[CRT]sparseMatrix>, k1, k2), tri[ul](<[CRT]sparseMatrix>, k) */
573/* NB: argument validation more or less copied from R_dense_band() */
574SEXP R_sparse_band(SEXP from, SEXP k1, SEXP k2)
575{
576 static const char *valid[] = {
577 VALID_CSPARSE, VALID_RSPARSE, VALID_TSPARSE, "" };
578 int ivalid = R_check_class_etc(from, valid);
579 if (ivalid < 0)
580 ERROR_INVALID_CLASS(from, __func__);
581
582 SEXP dim = PROTECT(GET_SLOT(from, Matrix_DimSym));
583 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1];
584 UNPROTECT(1);
585
586 int a, b;
587 if (k1 == R_NilValue) // tril()
588 a = -m ; // was (m > 0) ? 1 - m : 0;
589 else if ((a = asInteger(k1)) == NA_INTEGER || a < -m || a > n)
590 error(_("'%s' (%d) must be an integer from %s (%d) to %s (%d)"),
591 "k1", a, "-Dim[1]", -m, "Dim[2]", n);
592 if (k2 == R_NilValue) // triu()
593 b = n; // was (n > 0) ? n - 1 : 0;
594 else if ((b = asInteger(k2)) == NA_INTEGER || b < -m || b > n)
595 error(_("'%s' (%d) must be an integer from %s (%d) to %s (%d)"),
596 "k2", b, "-Dim[1]", -m, "Dim[2]", n);
597 else if (b < a)
598 error(_("'%s' (%d) must be less than or equal to '%s' (%d)"),
599 "k1", a, "k2", b);
600
601 return sparse_band(from, valid[ivalid], a, b);
602}
603
604SEXP sparse_diag_get(SEXP obj, const char *class, int names)
605{
606 SEXP dim = PROTECT(GET_SLOT(obj, Matrix_DimSym));
607 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1], r = (m < n) ? m : n;
608 UNPROTECT(1); /* dim */
609
610 char ul = 'U', di = 'N';
611 if (class[1] != 'g') {
612 SEXP uplo = PROTECT(GET_SLOT(obj, Matrix_uploSym));
613 ul = *CHAR(STRING_ELT(uplo, 0));
614 UNPROTECT(1); /* uplo */
615
616 if (class[1] == 't') {
617 SEXP diag = PROTECT(GET_SLOT(obj, Matrix_diagSym));
618 di = *CHAR(STRING_ELT(diag, 0));
619 UNPROTECT(1); /* diag */
620 }
621 }
622
623 SEXP res = PROTECT(allocVector(kindToType(class[0]), r));
624
625#define DG_CASES \
626 do { \
627 switch (class[0]) { \
628 case 'n': \
629 case 'l': \
630 DG_LOOP(int, LOGICAL, SHOW, FIRSTOF, INCREMENT_LOGICAL, 0, 1); \
631 break; \
632 case 'i': \
633 DG_LOOP(int, INTEGER, SHOW, FIRSTOF, INCREMENT_INTEGER, 0, 1); \
634 break; \
635 case 'd': \
636 DG_LOOP(double, REAL, SHOW, FIRSTOF, INCREMENT_REAL, 0.0, 1.0); \
637 break; \
638 case 'z': \
639 DG_LOOP(Rcomplex, COMPLEX, SHOW, FIRSTOF, INCREMENT_COMPLEX, Matrix_zzero, Matrix_zone); \
640 break; \
641 default: \
642 break; \
643 } \
644 } while (0)
645
646 if (di != 'N') {
647
648 int j;
649
650#undef DG_LOOP
651#define DG_LOOP(_CTYPE_, _PTR_, _MASK_, _REPLACE_, _INCREMENT_, _ZERO_, _ONE_) \
652 do { \
653 _CTYPE_ *pres = _PTR_(res); \
654 for (j = 0; j < r; ++j) \
655 *(pres++) = _ONE_; \
656 } while (0)
657
658 DG_CASES;
659
660 } else if (class[2] != 'T') {
661
662 SEXP iSym = (class[2] == 'C') ? Matrix_iSym : Matrix_jSym,
663 p0 = PROTECT(GET_SLOT(obj, Matrix_pSym)),
664 i0 = PROTECT(GET_SLOT(obj, iSym));
665 int j, k, kend, *pp0 = INTEGER(p0), *pi0 = INTEGER(i0);
666 pp0++;
667
668#undef DG_LOOP
669#define DG_LOOP(_CTYPE_, _PTR_, _MASK_, _REPLACE_, _INCREMENT_, _ZERO_, _ONE_) \
670 do { \
671 _MASK_(_CTYPE_ *px0 = _PTR_(x0)); \
672 _CTYPE_ *pres = _PTR_(res); \
673 if (class[1] == 'g') { \
674 for (j = 0, k = 0; j < r; ++j) { \
675 pres[j] = _ZERO_; \
676 kend = pp0[j]; \
677 while (k < kend) { \
678 if (pi0[k] != j) \
679 ++k; \
680 else { \
681 pres[j] = _REPLACE_(px0[k], 1); \
682 k = kend; \
683 } \
684 } \
685 } \
686 } else if ((class[2] == 'C') == (ul != 'U')) { \
687 for (j = 0, k = 0; j < r; ++j) { \
688 kend = pp0[j]; \
689 pres[j] = (k < kend && pi0[k] == j) \
690 ? _REPLACE_(px0[k], 1) : _ZERO_; \
691 k = kend; \
692 } \
693 } else { \
694 for (j = 0, k = 0; j < r; ++j) { \
695 kend = pp0[j]; \
696 pres[j] = (k < kend && pi0[kend - 1] == j) \
697 ? _REPLACE_(px0[kend - 1], 1) : _ZERO_; \
698 k = kend; \
699 } \
700 } \
701 } while (0)
702
703 if (class[0] == 'n') {
704 DG_LOOP(int, LOGICAL, HIDE, SECONDOF, , 0, 1);
705 } else {
706 SEXP x0 = PROTECT(GET_SLOT(obj, Matrix_xSym));
707 DG_CASES;
708 UNPROTECT(1); /* x0 */
709 }
710
711 UNPROTECT(2); /* i0, p0 */
712
713 } else {
714
715 SEXP i0 = PROTECT(GET_SLOT(obj, Matrix_iSym)),
716 j0 = PROTECT(GET_SLOT(obj, Matrix_jSym));
717 int *pi0 = INTEGER(i0), *pj0 = INTEGER(j0);
718 R_xlen_t k, nnz0 = XLENGTH(i0);
719
720#undef DG_LOOP
721#define DG_LOOP(_CTYPE_, _PTR_, _MASK_, _REPLACE_, _INCREMENT_, _ZERO_, _ONE_) \
722 do { \
723 _MASK_(_CTYPE_ *px0 = _PTR_(x0)); \
724 _CTYPE_ *pres = _PTR_(res); \
725 Matrix_memset(pres, 0, r, sizeof(_CTYPE_)); \
726 for (k = 0; k < nnz0; ++k) { \
727 if (*pi0 == *pj0) \
728 _INCREMENT_(pres[*pi0], (*px0)); \
729 ++pi0; ++pj0; _MASK_(++px0); \
730 } \
731 } while (0)
732
733 if (class[0] == 'n')
734 DG_LOOP(int, LOGICAL, HIDE, SECONDOF, INCREMENT_PATTERN, 0, 1);
735 else {
736 SEXP x0 = PROTECT(GET_SLOT(obj, Matrix_xSym));
737 DG_CASES;
738 UNPROTECT(1); /* x0 */
739 }
740
741 UNPROTECT(2); /* j0, i0 */
742
743 }
744
745#undef DG_CASES
746#undef DG_LOOP
747
748 if (names) {
749 /* NB: The logic here must be adjusted once the validity method
750 for 'symmetricMatrix' enforces symmetric 'Dimnames'
751 */
752 SEXP dn = PROTECT(GET_SLOT(obj, Matrix_DimNamesSym)),
753 rn = VECTOR_ELT(dn, 0),
754 cn = VECTOR_ELT(dn, 1);
755 if (cn == R_NilValue) {
756 if (class[1] == 's' && rn != R_NilValue)
757 setAttrib(res, R_NamesSymbol, rn);
758 } else {
759 if (class[1] == 's')
760 setAttrib(res, R_NamesSymbol, cn);
761 else if (rn != R_NilValue &&
762 (rn == cn || equal_character_vectors(rn, cn, r)))
763 setAttrib(res, R_NamesSymbol, (r == m) ? rn : cn);
764 }
765 UNPROTECT(1); /* dn */
766 }
767
768 UNPROTECT(1); /* res */
769 return res;
770}
771
772/* diag(<[CRT]sparseMatrix>, names=) */
773SEXP R_sparse_diag_get(SEXP obj, SEXP names)
774{
775 static const char *valid[] = {
776 VALID_CSPARSE, VALID_RSPARSE, VALID_TSPARSE, "" };
777 int ivalid = R_check_class_etc(obj, valid);
778 if (ivalid < 0)
779 ERROR_INVALID_CLASS(obj, __func__);
780
781 int names_;
782 if (TYPEOF(names) != LGLSXP || LENGTH(names) < 1 ||
783 (names_ = LOGICAL(names)[0]) == NA_LOGICAL)
784 error(_("'%s' must be %s or %s"), "names", "TRUE", "FALSE");
785
786 return sparse_diag_get(obj, valid[ivalid], names_);
787}
788
789SEXP sparse_diag_set(SEXP from, const char *class, SEXP value)
790{
791 SEXP to = PROTECT(newObject(class));
792 int v = LENGTH(value) != 1, full = 0;
793
794 SEXP dim = PROTECT(GET_SLOT(from, Matrix_DimSym));
795 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1], r = (m < n) ? m : n;
796 if (m != n || n > 0)
797 SET_SLOT(to, Matrix_DimSym, dim);
798 UNPROTECT(1); /* dim */
799
800 SEXP dimnames = PROTECT(GET_SLOT(from, Matrix_DimNamesSym));
801 SET_SLOT(to, Matrix_DimNamesSym, dimnames);
802 UNPROTECT(1); /* dimnames */
803
804 char ul = 'U', di = 'N';
805 if (class[1] != 'g') {
806 SEXP uplo = PROTECT(GET_SLOT(from, Matrix_uploSym));
807 ul = *CHAR(STRING_ELT(uplo, 0));
808 if (ul != 'U')
809 SET_SLOT(to, Matrix_uploSym, uplo);
810 UNPROTECT(1); /* uplo */
811
812 if (class[1] == 't') {
813 SEXP diag = PROTECT(GET_SLOT(from, Matrix_diagSym));
814 di = *CHAR(STRING_ELT(diag, 0));
815 UNPROTECT(1); /* diag */
816 }
817 }
818
819#define DS_CASES \
820 do { \
821 switch (class[0]) { \
822 case 'n': \
823 case 'l': \
824 DS_LOOP(int, LOGICAL, SHOW, ISNZ_LOGICAL); \
825 break; \
826 case 'i': \
827 DS_LOOP(int, INTEGER, SHOW, ISNZ_INTEGER); \
828 break; \
829 case 'd': \
830 DS_LOOP(double, REAL, SHOW, ISNZ_REAL); \
831 break; \
832 case 'z': \
833 DS_LOOP(Rcomplex, COMPLEX, SHOW, ISNZ_COMPLEX); \
834 break; \
835 default: \
836 break; \
837 } \
838 } while (0)
839
840 if (class[2] != 'T') {
841
842 int n_ = (class[2] == 'C') ? n : m;
843 SEXP iSym = (class[2] == 'C') ? Matrix_iSym : Matrix_jSym,
844 p0 = PROTECT(GET_SLOT(from, Matrix_pSym)),
845 i0 = PROTECT(GET_SLOT(from, iSym)),
846 p1 = PROTECT(allocVector(INTSXP, (R_xlen_t) n_ + 1));
847 SET_SLOT(to, Matrix_pSym, p1);
848 int *pp0 = INTEGER(p0), *pp1 = INTEGER(p1), *pi0 = INTEGER(i0),
849 j, k, kend, nd0 = 0, nd1 = 0;
850 pp0++; *(pp1++) = 0;
851
852 if (class[1] == 'g') {
853 for (j = 0, k = 0; j < r; ++j) {
854 kend = pp0[j];
855 while (k < kend) {
856 if (pi0[k] < j)
857 ++k;
858 else {
859 if (pi0[k] == j)
860 ++nd0;
861 k = kend;
862 }
863 }
864 pp1[j] = kend - nd0;
865 }
866 for (j = r; j < n_; ++j)
867 pp1[j] = pp0[j] - nd0;
868 } else if (di != 'N') {
869 for (j = 0; j < n_; ++j)
870 pp1[j] = pp0[j];
871 } else if ((class[2] == 'C') == (ul == 'U')) {
872 for (j = 0, k = 0; j < n_; ++j) {
873 kend = pp0[j];
874 if (k < kend && pi0[kend - 1] == j)
875 ++nd0;
876 k = kend;
877 pp1[j] = kend - nd0;
878 }
879 } else {
880 for (j = 0, k = 0; j < n_; ++j) {
881 kend = pp0[j];
882 if (k < kend && pi0[k] == j)
883 ++nd0;
884 k = kend;
885 pp1[j] = kend - nd0;
886 }
887 }
888
889#undef DS_LOOP
890#define DS_LOOP(_CTYPE_, _PTR_, _MASK_, _ISNZ_) \
891 do { \
892 _CTYPE_ *pvalue = _PTR_(value); \
893 if (v) { \
894 for (j = 0; j < r; ++j) { \
895 if (_ISNZ_(pvalue[j])) \
896 ++nd1; \
897 pp1[j] += nd1; \
898 } \
899 for (j = r; j < n_; ++j) \
900 pp1[j] += nd1; \
901 } else if (_ISNZ_(pvalue[0])) { \
902 full = 1; \
903 for (j = 0; j < r; ++j) \
904 pp1[j] += ++nd1; \
905 for (j = r; j < n_; ++j) \
906 pp1[j] += nd1; \
907 } \
908 } while (0)
909
910 DS_CASES;
911
912 if (nd1 - nd0 > INT_MAX - pp0[n_ - 1])
913 error(_("%s cannot exceed %s"), "p[length(p)]", "2^31-1");
914
915 SEXP i1 = PROTECT(allocVector(INTSXP, pp1[n_ - 1]));
916 SET_SLOT(to, iSym, i1);
917 int *pi1 = INTEGER(i1);
918
919#undef DS_LOOP
920#define DS_LOOP(_CTYPE_, _PTR_, _MASK_, _ISNZ_) \
921 do { \
922 _MASK_(_CTYPE_ *px0 = _PTR_(x0), *px1 = _PTR_(x1)); \
923 _CTYPE_ *pvalue = _PTR_(value); \
924 for (j = 0, k = 0; j < r; ++j) { \
925 kend = pp0[j]; \
926 while (k < kend && pi0[k] < j) { \
927 *(pi1++) = pi0[k] ; \
928 _MASK_(*(px1++) = px0[k]); \
929 ++k; \
930 } \
931 if (k < kend && pi0[k] == j) \
932 ++k; \
933 if ((v) ? _ISNZ_(pvalue[j]) : full) { \
934 *(pi1++) = j ; \
935 _MASK_(*(px1++) = (v) ? pvalue[j] : pvalue[0]); \
936 } \
937 while (k < kend) { \
938 *(pi1++) = pi0[k] ; \
939 _MASK_(*(px1++) = px0[k]); \
940 ++k; \
941 } \
942 } \
943 for (j = r; j < n_; ++j) { \
944 kend = pp0[j]; \
945 while (k < kend) { \
946 *(pi1++) = pi0[k] ; \
947 _MASK_(*(px1++) = px0[k]); \
948 ++k; \
949 } \
950 } \
951 } while (0)
952
953 if (class[0] == 'n')
954 DS_LOOP(int, LOGICAL, HIDE, ISNZ_LOGICAL);
955 else {
956 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym)),
957 x1 = PROTECT(allocVector(TYPEOF(x0), pp1[n_ - 1]));
958 SET_SLOT(to, Matrix_xSym, x1);
959 DS_CASES;
960 UNPROTECT(2); /* x1, x0 */
961 }
962
963 UNPROTECT(4); /* i1, p1, i0, p0 */
964
965 } else {
966
967 SEXP i0 = PROTECT(GET_SLOT(from, Matrix_iSym)),
968 j0 = PROTECT(GET_SLOT(from, Matrix_jSym));
969 int *pi0 = INTEGER(i0), *pj0 = INTEGER(j0), j, nd0 = 0, nd1 = 0;
970 R_xlen_t k, nnz0 = XLENGTH(i0), nnz1 = nnz0;
971
972 if (di == 'N')
973 for (k = 0; k < nnz0; ++k)
974 if (pi0[k] == pj0[k])
975 ++nd0;
976
977#undef DS_LOOP
978#define DS_LOOP(_CTYPE_, _PTR_, _MASK_, _ISNZ_) \
979 do { \
980 _CTYPE_ *pvalue = _PTR_(value); \
981 if (v) { \
982 for (j = 0; j < r; ++j) \
983 if (_ISNZ_(pvalue[j])) \
984 ++nd1; \
985 } else if (_ISNZ_(pvalue[0])) { \
986 full = 1; \
987 nd1 = r; \
988 } \
989 } while (0)
990
991 DS_CASES;
992
993 if (nd1 - nd0 > R_XLEN_T_MAX - nnz0)
994 error(_("%s cannot exceed %s"), "length(i)", "R_XLEN_T_MAX");
995 nnz1 += nd1 - nd0;
996
997 SEXP i1 = PROTECT(allocVector(INTSXP, nnz1)),
998 j1 = PROTECT(allocVector(INTSXP, nnz1));
999 SET_SLOT(to, Matrix_iSym, i1);
1000 SET_SLOT(to, Matrix_jSym, j1);
1001 int *pi1 = INTEGER(i1), *pj1 = INTEGER(j1);
1002
1003#undef DS_LOOP
1004#define DS_LOOP(_CTYPE_, _PTR_, _MASK_, _ISNZ_) \
1005 do { \
1006 _MASK_(_CTYPE_ *px0 = _PTR_(x0), *px1 = _PTR_(x1)); \
1007 _CTYPE_ *pvalue = _PTR_(value); \
1008 for (k = 0; k < nnz0; ++k) { \
1009 if (pi0[k] != pj0[k]) { \
1010 *(pi1++) = pi0[k]; \
1011 *(pj1++) = pj0[k]; \
1012 _MASK_(*(px1++) = px0[k]); \
1013 } \
1014 } \
1015 if (v) { \
1016 for (j = 0; j < r; ++j) { \
1017 if (_ISNZ_(pvalue[j])) { \
1018 *(pi1++) = *(pj1++) = j; \
1019 _MASK_(*(px1++) = pvalue[j]); \
1020 } \
1021 } \
1022 } else if (full) { \
1023 for (j = 0; j < r; ++j) { \
1024 *(pi1++) = *(pj1++) = j; \
1025 _MASK_(*(px1++) = pvalue[0]); \
1026 } \
1027 } \
1028 } while (0)
1029
1030 if (class[0] == 'n')
1031 DS_LOOP(int, LOGICAL, HIDE, ISNZ_LOGICAL);
1032 else {
1033 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym)),
1034 x1 = PROTECT(allocVector(TYPEOF(x0), nnz1));
1035 SET_SLOT(to, Matrix_xSym, x1);
1036 DS_CASES;
1037 UNPROTECT(2); /* x1, x0 */
1038 }
1039
1040 UNPROTECT(4); /* j1, i1, j0, i0 */
1041
1042 }
1043
1044#undef DS_CASES
1045#undef DS_LOOP
1046
1047 UNPROTECT(1); /* to */
1048 return to;
1049}
1050
1051/* diag(<[CRT]sparseMatrix>) <- value */
1052SEXP R_sparse_diag_set(SEXP from, SEXP value)
1053{
1054 static const char *valid[] = {
1055 VALID_CSPARSE, VALID_RSPARSE, VALID_TSPARSE, "" };
1056 int ivalid = R_check_class_etc(from, valid);
1057 if (ivalid < 0)
1058 ERROR_INVALID_CLASS(from, __func__);
1059 const char *class = valid[ivalid];
1060
1061 SEXPTYPE tx = kindToType(class[0]), tv = TYPEOF(value);
1062
1063 switch (tv) {
1064 case LGLSXP:
1065 case INTSXP:
1066 case REALSXP:
1067 case CPLXSXP:
1068 break;
1069 default:
1070 error(_("replacement diagonal has incompatible type \"%s\""),
1071 type2char(tv));
1072 break;
1073 }
1074
1075 SEXP dim = GET_SLOT(from, Matrix_DimSym);
1076 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1], r = (m < n) ? m : n;
1077 R_xlen_t len = XLENGTH(value);
1078 if (len != 1 && len != r)
1079 error(_("replacement diagonal has wrong length"));
1080
1081 if (tv <= tx) {
1082 PROTECT(from);
1083 PROTECT(value = coerceVector(value, tx));
1084 } else {
1085 /* defined in ./coerce.c : */
1086 SEXP sparse_as_kind(SEXP, const char *, char);
1087#ifndef MATRIX_ENABLE_IMATRIX
1088 if (tv == INTSXP) {
1089 PROTECT(from = sparse_as_kind(from, class, 'd'));
1090 PROTECT(value = coerceVector(value, REALSXP));
1091 } else {
1092#endif
1093 PROTECT(from = sparse_as_kind(from, class, typeToKind(tv)));
1094 PROTECT(value);
1095#ifndef MATRIX_ENABLE_IMATRIX
1096 }
1097#endif
1098 class = valid[R_check_class_etc(from, valid)];
1099 }
1100
1101 from = sparse_diag_set(from, class, value);
1102 UNPROTECT(2);
1103 return from;
1104}
1105
1106SEXP sparse_transpose(SEXP from, const char *class, int lazy)
1107{
1108 SEXP to;
1109 if (class[2] == 'T' || !lazy)
1110 PROTECT(to = newObject(class));
1111 else {
1112 char cl[] = "...Matrix";
1113 cl[0] = class[0];
1114 cl[1] = class[1];
1115 cl[2] = (class[2] == 'C') ? 'R' : 'C';
1116 PROTECT(to = newObject(cl));
1117 }
1118
1119 SEXP dim = PROTECT(GET_SLOT(from, Matrix_DimSym));
1120 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1];
1121 if (m != n) {
1122 UNPROTECT(1); /* dim */
1123 PROTECT(dim = GET_SLOT(to, Matrix_DimSym));
1124 pdim = INTEGER(dim);
1125 pdim[0] = n;
1126 pdim[1] = m;
1127 } else if (n > 0)
1128 SET_SLOT(to, Matrix_DimSym, dim);
1129 UNPROTECT(1); /* dim */
1130
1131 SEXP dimnames = PROTECT(GET_SLOT(from, Matrix_DimNamesSym));
1132 if (class[1] == 's')
1133 SET_SLOT(to, Matrix_DimNamesSym, dimnames);
1134 else
1135 set_reversed_DimNames(to, dimnames);
1136 UNPROTECT(1); /* dimnames */
1137
1138 if (class[1] != 'g') {
1139 SEXP uplo = PROTECT(GET_SLOT(from, Matrix_uploSym));
1140 char ul = *CHAR(STRING_ELT(uplo, 0));
1141 UNPROTECT(1); /* uplo */
1142 if (ul == 'U') {
1143 PROTECT(uplo = mkString("L"));
1144 SET_SLOT(to, Matrix_uploSym, uplo);
1145 UNPROTECT(1); /* uplo */
1146 }
1147 if (class[1] == 't') {
1148 SEXP diag = PROTECT(GET_SLOT(from, Matrix_diagSym));
1149 char di = *CHAR(STRING_ELT(diag, 0));
1150 if (di != 'N')
1151 SET_SLOT(to, Matrix_diagSym, diag);
1152 UNPROTECT(1); /* diag */
1153 } else {
1154 SEXP factors = PROTECT(GET_SLOT(from, Matrix_factorsSym));
1155 if (LENGTH(factors) > 0)
1156 SET_SLOT(to, Matrix_factorsSym, factors);
1157 UNPROTECT(1); /* factors */
1158 }
1159 }
1160
1161 /* It remains to set some subset of 'p', 'i', 'j', and 'x' ... */
1162
1163 if (class[2] == 'T') {
1164 /* No need to allocate in this case: need only reverse 'i' and 'j' */
1165 SEXP i0 = PROTECT(GET_SLOT(from, Matrix_iSym)),
1166 j0 = PROTECT(GET_SLOT(from, Matrix_jSym));
1167 SET_SLOT(to, Matrix_iSym, j0);
1168 SET_SLOT(to, Matrix_jSym, i0);
1169 UNPROTECT(2); /* j0, i0 */
1170 if (class[0] != 'n') {
1171 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym));
1172 SET_SLOT(to, Matrix_xSym, x0);
1173 UNPROTECT(1); /* x */
1174 }
1175 UNPROTECT(1); /* to */
1176 return to;
1177 }
1178
1179 /* Now dealing only with [CR]sparseMatrix ... */
1180
1181 SEXP iSym = (class[2] == 'C') ? Matrix_iSym : Matrix_jSym,
1182 p0 = PROTECT(GET_SLOT(from, Matrix_pSym)),
1183 i0 = PROTECT(GET_SLOT(from, iSym));
1184
1185 if (lazy) {
1186 /* No need to allocate in this case: need only reverse 'i' and 'j' */
1187 SEXP jSym = (class[2] == 'C') ? Matrix_jSym : Matrix_iSym;
1188 SET_SLOT(to, Matrix_pSym, p0);
1189 SET_SLOT(to, jSym, i0);
1190 UNPROTECT(2); /* i0, p0 */
1191 if (class[0] != 'n') {
1192 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym));
1193 SET_SLOT(to, Matrix_xSym, x0);
1194 UNPROTECT(1); /* x */
1195 }
1196 UNPROTECT(1); /* to */
1197 return to;
1198 }
1199
1200 int m_ = (class[2] == 'C') ? m : n, n_ = (class[2] == 'C') ? n : m;
1201 SEXP p1 = PROTECT(allocVector(INTSXP, (R_xlen_t) m_ + 1)),
1202 i1 = PROTECT(allocVector(INTSXP, INTEGER(p0)[n_]));
1203 SET_SLOT(to, Matrix_pSym, p1);
1204 SET_SLOT(to, iSym, i1);
1205
1206 /* defined in ./coerce.c : */
1207 void trans(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, int, int);
1208
1209 if (class[0] == 'n')
1210 trans(p0, i0, NULL, p1, i1, NULL, m_, n_);
1211 else {
1212 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym)),
1213 x1 = PROTECT(allocVector(TYPEOF(x0), INTEGER(p0)[n_]));
1214 SET_SLOT(to, Matrix_xSym, x1);
1215 trans(p0, i0, x0, p1, i1, x1, m_, n_);
1216 UNPROTECT(2); /* x1, x0 */
1217 }
1218 UNPROTECT(5); /* i1, p1, i0, p0, to */
1219 return to;
1220}
1221
1222/* t(<[CRT]sparseMatrix>) */
1223SEXP R_sparse_transpose(SEXP from, SEXP lazy)
1224{
1225 static const char *valid[] = {
1226 VALID_CSPARSE, VALID_RSPARSE, VALID_TSPARSE, "" };
1227 int ivalid = R_check_class_etc(from, valid);
1228 if (ivalid < 0)
1229 ERROR_INVALID_CLASS(from, __func__);
1230
1231 int lazy_;
1232 if (TYPEOF(lazy) != LGLSXP || LENGTH(lazy) < 1 ||
1233 (lazy_ = LOGICAL(lazy)[0]) == NA_LOGICAL)
1234 error(_("invalid '%s' to '%s'"), "lazy", __func__);
1235
1236 return sparse_transpose(from, valid[ivalid], lazy_);
1237}
1238
1239SEXP sparse_force_symmetric(SEXP from, const char *class, char ul)
1240{
1241 char ul0 = 'U', ul1 = 'U';
1242 if (class[1] != 'g') {
1243 SEXP uplo = PROTECT(GET_SLOT(from, Matrix_uploSym));
1244 ul0 = ul1 = *CHAR(STRING_ELT(uplo, 0));
1245 UNPROTECT(1); /* uplo */
1246 }
1247 if (ul != '\0')
1248 ul1 = ul;
1249
1250 if (class[1] == 's') {
1251 /* .s[CRT]Matrix */
1252 if (ul0 == ul1)
1253 return from;
1254 SEXP to = PROTECT(sparse_transpose(from, class, 0));
1255 if (class[0] == 'z') {
1256 /* Need _conjugate_ transpose */
1257 SEXP x1 = PROTECT(GET_SLOT(to, Matrix_xSym));
1258 conjugate(x1);
1259 UNPROTECT(1); /* x1 */
1260 }
1261 UNPROTECT(1) /* to */;
1262 return to;
1263 }
1264
1265 /* Now handling just .[gt][CRT]Matrix ... */
1266
1267 char cl[] = ".s.Matrix";
1268 cl[0] = class[0];
1269 cl[2] = class[2];
1270 SEXP to = PROTECT(newObject(cl));
1271
1272 SEXP dim = PROTECT(GET_SLOT(from, Matrix_DimSym));
1273 int *pdim = INTEGER(dim), n = pdim[0];
1274 if (pdim[1] != n)
1275 error(_("attempt to symmetrize a non-square matrix"));
1276 if (n > 0)
1277 SET_SLOT(to, Matrix_DimSym, dim);
1278 UNPROTECT(1); /* dim */
1279
1280 SEXP dimnames = PROTECT(GET_SLOT(from, Matrix_DimNamesSym));
1281 set_symmetrized_DimNames(to, dimnames, -1);
1282 UNPROTECT(1); /* dimnames */
1283
1284 if (ul1 != 'U') {
1285 SEXP uplo = PROTECT(mkString("L"));
1286 SET_SLOT(to, Matrix_uploSym, uplo);
1287 UNPROTECT(1); /* uplo */
1288 }
1289
1290 char di = 'N';
1291 if (class[1] == 't') {
1292 SEXP diag = PROTECT(GET_SLOT(from, Matrix_diagSym));
1293 di = *CHAR(STRING_ELT(diag, 0));
1294 UNPROTECT(1); /* diag */
1295 }
1296
1297 /* It remains to set some subset of 'p', 'i', 'j', and 'x' ... */
1298
1299#define FS_CASES \
1300 do { \
1301 switch (class[0]) { \
1302 case 'l': \
1303 FS_SUBCASES(int, LOGICAL, SHOW, 1); \
1304 break; \
1305 case 'i': \
1306 FS_SUBCASES(int, INTEGER, SHOW, 1); \
1307 break; \
1308 case 'd': \
1309 FS_SUBCASES(double, REAL, SHOW, 1.0); \
1310 break; \
1311 case 'z': \
1312 FS_SUBCASES(Rcomplex, COMPLEX, SHOW, Matrix_zone); \
1313 break; \
1314 default: \
1315 break; \
1316 } \
1317 } while (0)
1318
1319 if (class[1] == 't' && di == 'N' && ul0 == ul1) {
1320
1321 /* No need to allocate in this case: we have the triangle we want */
1322 if (class[2] != 'T') {
1323 SEXP p0 = PROTECT(GET_SLOT(from, Matrix_pSym));
1324 SET_SLOT(to, Matrix_pSym, p0);
1325 UNPROTECT(1); /* p0 */
1326 }
1327 if (class[2] != 'R') {
1328 SEXP i0 = PROTECT(GET_SLOT(from, Matrix_iSym));
1329 SET_SLOT(to, Matrix_iSym, i0);
1330 UNPROTECT(1); /* i0 */
1331 }
1332 if (class[2] != 'C') {
1333 SEXP j0 = PROTECT(GET_SLOT(from, Matrix_jSym));
1334 SET_SLOT(to, Matrix_jSym, j0);
1335 UNPROTECT(1); /* j0 */
1336 }
1337 if (class[0] != 'n') {
1338 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym));
1339 SET_SLOT(to, Matrix_xSym, x0);
1340 UNPROTECT(1); /* x0 */
1341 }
1342 UNPROTECT(1); /* to */
1343 return to;
1344
1345 } else if (class[2] != 'T') {
1346
1347 /* Symmetrizing square .[gt][CR]Matrix ... */
1348
1349 SEXP iSym = (class[2] == 'C') ? Matrix_iSym : Matrix_jSym,
1350 p0 = PROTECT(GET_SLOT(from, Matrix_pSym)),
1351 p1 = PROTECT(allocVector(INTSXP, (R_xlen_t) n + 1)),
1352 i0 = PROTECT(GET_SLOT(from, iSym));
1353 int j, k, kend,
1354 *pp0 = INTEGER(p0),
1355 *pp1 = INTEGER(p1),
1356 *pi0 = INTEGER(i0),
1357 nnz0 = pp0[n],
1358 nnz1 = 0;
1359 pp0++; *(pp1++) = 0;
1360 SET_SLOT(to, Matrix_pSym, p1);
1361
1362 /* Counting number of nonzero elements in triangle, by "column" ... */
1363
1364 if (class[1] == 't') {
1365 if (di != 'N') {
1366 /* Have triangular matrix with unit diagonal */
1367 if (ul0 != ul1) {
1368 /* Returning identity matrix */
1369 for (j = 0; j < n; ++j)
1370 pp1[j] = ++nnz1;
1371 } else {
1372 /* Returning symmetric matrix with unit diagonal */
1373 for (j = 0; j < n; ++j)
1374 pp1[j] = ++nnz1 + pp0[j];
1375 nnz1 += nnz0;
1376 }
1377 } else if (ul0 == ((class[2] == 'C') ? 'U' : 'L')) {
1378 /* Have triangular matrix with non-unit "trailing" diagonal
1379 and returning diagonal part */
1380 for (j = 0; j < n; ++j) {
1381 if (pp0[j - 1] < pp0[j] && pi0[pp0[j] - 1] == j)
1382 ++nnz1;
1383 pp1[j] = nnz1;
1384 }
1385 } else {
1386 /* Have triangular matrix with non-unit "leading" diagonal
1387 and returning diagonal part */
1388 for (j = 0; j < n; ++j) {
1389 if (pp0[j - 1] < pp0[j] && pi0[pp0[j - 1]] == j)
1390 ++nnz1;
1391 pp1[j] = nnz1;
1392 }
1393 }
1394 } else if (ul1 == ((class[2] == 'C') ? 'U' : 'L')) {
1395 /* Have general matrix and returning upper triangle */
1396 for (j = 0, k = 0; j < n; ++j) {
1397 kend = pp0[j];
1398 while (k < kend) {
1399 if (pi0[k] <= j)
1400 ++nnz1;
1401 ++k;
1402 }
1403 pp1[j] = nnz1;
1404 }
1405 } else {
1406 /* Have general matrix and returning lower triangle */
1407 for (j = 0, k = 0; j < n; ++j) {
1408 kend = pp0[j];
1409 while (k < kend) {
1410 if (pi0[k] >= j)
1411 ++nnz1;
1412 ++k;
1413 }
1414 pp1[j] = nnz1;
1415 }
1416 }
1417
1418 /* Now allocating and filling out slots ... */
1419
1420 SEXP i1 = PROTECT(allocVector(INTSXP, nnz1));
1421 int *pi1 = INTEGER(i1);
1422 SET_SLOT(to, iSym, i1);
1423
1424#undef FS_SUBCASES
1425#define FS_SUBCASES(_CTYPE_, _PTR_, _MASK_, _ONE_) \
1426 do { \
1427 _MASK_(_CTYPE_ *px0 = _PTR_(x0), *px1 = _PTR_(x1)); \
1428 if (class[1] == 't') { \
1429 if (di != 'N') { \
1430 /* Have triangular matrix with unit diagonal */ \
1431 if (ul0 != ul1) { \
1432 /* Returning identity matrix */ \
1433 for (j = 0; j < n; ++j) { \
1434 *(pi1++) = j; \
1435 _MASK_(*(px1++) = _ONE_); \
1436 } \
1437 } else if (ul0 == ((class[2] == 'C') ? 'U' : 'L')) { \
1438 /* Returning symmetric matrix */ \
1439 /* with unit "trailing" diagonal */ \
1440 for (j = 0, k = 0; j < n; ++j) { \
1441 kend = pp0[j]; \
1442 while (k < kend) { \
1443 *(pi1++) = pi0[k]; \
1444 _MASK_(*(px1++) = px0[k]); \
1445 ++k; \
1446 } \
1447 *(pi1++) = j; \
1448 _MASK_(*(px1++) = _ONE_); \
1449 } \
1450 } else { \
1451 /* Returning symmetric matrix */ \
1452 /* with unit "leading" diagonal */ \
1453 for (j = 0, k = 0; j < n; ++j) { \
1454 *(pi1++) = j; \
1455 _MASK_(*(px1++) = _ONE_); \
1456 kend = pp0[j]; \
1457 while (k < kend) { \
1458 *(pi1++) = pi0[k]; \
1459 _MASK_(*(px1++) = px0[k]); \
1460 ++k; \
1461 } \
1462 } \
1463 } \
1464 } else if (ul0 == ((class[2] == 'C') ? 'U' : 'L')) { \
1465 /* Have triangular matrix with non-unit "trailing" */ \
1466 /* diagonal and returning diagonal part */ \
1467 for (j = 0; j < n; ++j) { \
1468 if (pp0[j - 1] < pp0[j] && pi0[pp0[j] - 1] == j) { \
1469 *(pi1++) = j; \
1470 _MASK_(*(px1++) = px0[pp0[j] - 1]); \
1471 } \
1472 } \
1473 } else { \
1474 /* Have triangular matrix with non-unit "leading" */ \
1475 /* diagonal and returning diagonal part */ \
1476 for (j = 0; j < n; ++j) { \
1477 if (pp0[j - 1] < pp0[j] && pi0[pp0[j - 1]] == j) { \
1478 *(pi1++) = j; \
1479 _MASK_(*(px1++) = px0[pp0[j - 1]]); \
1480 } \
1481 } \
1482 } \
1483 } else if (ul1 == ((class[2] == 'C') ? 'U' : 'L')) { \
1484 /* Have general matrix and returning upper triangle */ \
1485 for (j = 0, k = 0; j < n; ++j) { \
1486 kend = pp0[j]; \
1487 while (k < kend) { \
1488 if (pi0[k] <= j) { \
1489 *(pi1++) = pi0[k]; \
1490 _MASK_(*(px1++) = px0[k]); \
1491 } \
1492 ++k; \
1493 } \
1494 } \
1495 } else { \
1496 /* Have general matrix and returning lower triangle */ \
1497 for (j = 0, k = 0; j < n; ++j) { \
1498 kend = pp0[j]; \
1499 while (k < kend) { \
1500 if (pi0[k] >= j) { \
1501 *(pi1++) = pi0[k]; \
1502 _MASK_(*(px1++) = px0[k]); \
1503 } \
1504 ++k; \
1505 } \
1506 } \
1507 } \
1508 } while (0)
1509
1510 if (class[0] == 'n')
1511 FS_SUBCASES(int, LOGICAL, HIDE, 1);
1512 else {
1513 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym)),
1514 x1 = PROTECT(allocVector(TYPEOF(x0), nnz1));
1515 SET_SLOT(to, Matrix_xSym, x1);
1516 FS_CASES;
1517 UNPROTECT(2); /* x1, x0 */
1518 }
1519
1520 } else {
1521
1522 /* Symmetrizing square .[gt]TMatrix ... */
1523
1524 SEXP i0 = PROTECT(GET_SLOT(from, Matrix_iSym)),
1525 j0 = PROTECT(GET_SLOT(from, Matrix_jSym));
1526 int *pi0 = INTEGER(i0),
1527 *pj0 = INTEGER(j0);
1528 R_xlen_t j, k, nnz0 = XLENGTH(i0), nnz1 = 0;
1529
1530 /* Counting number of nonzero elements in triangle ... */
1531
1532 if (class[1] == 't' && di != 'N')
1533 nnz1 = (ul0 == ul1) ? n + nnz0 : n;
1534 else {
1535 if (ul1 == 'U') {
1536 for (k = 0; k < nnz0; ++k)
1537 if (pi0[k] <= pj0[k])
1538 ++nnz1;
1539 } else {
1540 for (k = 0; k < nnz0; ++k)
1541 if (pi0[k] >= pj0[k])
1542 ++nnz1;
1543 }
1544 }
1545
1546 /* Now allocating and filling out slots ... */
1547
1548 SEXP i1 = PROTECT(allocVector(INTSXP, nnz1)),
1549 j1 = PROTECT(allocVector(INTSXP, nnz1));
1550 int *pi1 = INTEGER(i1),
1551 *pj1 = INTEGER(j1);
1552 SET_SLOT(to, Matrix_iSym, i1);
1553 SET_SLOT(to, Matrix_jSym, j1);
1554
1555#undef FS_SUBCASES
1556#define FS_SUBCASES(_CTYPE_, _PTR_, _MASK_, _ONE_) \
1557 do { \
1558 _MASK_(_CTYPE_ *px0 = _PTR_(x0), *px1 = _PTR_(x1)); \
1559 if (class[1] == 't' && di != 'N') { \
1560 if (ul0 == ul1) { \
1561 Matrix_memcpy(pi1, pi0, nnz0, sizeof(int)); \
1562 Matrix_memcpy(pj1, pj0, nnz0, sizeof(int)); \
1563 _MASK_(Matrix_memcpy(px1, px0, nnz0, sizeof(_CTYPE_))); \
1564 pi1 += nnz0; \
1565 pj1 += nnz0; \
1566 _MASK_(px1 += nnz0); \
1567 } \
1568 for (j = 0; j < n; ++j) { \
1569 *(pi1++) = *(pj1++) = j; \
1570 _MASK_(*(px1++) = _ONE_); \
1571 } \
1572 } else { \
1573 if (ul1 == 'U') { \
1574 for (k = 0; k < nnz0; ++k) { \
1575 if (pi0[k] <= pj0[k]) { \
1576 *(pi1++) = pi0[k]; \
1577 *(pj1++) = pj0[k]; \
1578 _MASK_(*(px1++) = px0[k]); \
1579 } \
1580 } \
1581 } else { \
1582 for (k = 0; k < nnz0; ++k) { \
1583 if (pi0[k] <= pj0[k]) { \
1584 *(pi1++) = pi0[k]; \
1585 *(pj1++) = pj0[k]; \
1586 _MASK_(*(px1++) = px0[k]); \
1587 } \
1588 } \
1589 } \
1590 } \
1591 } while (0)
1592
1593 if (class[0] == 'n')
1594 FS_SUBCASES(int, LOGICAL, HIDE, 1);
1595 else {
1596 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym)),
1597 x1 = PROTECT(allocVector(TYPEOF(x0), nnz1));
1598 SET_SLOT(to, Matrix_xSym, x1);
1599 FS_CASES;
1600 UNPROTECT(2); /* x1, x0 */
1601 }
1602
1603 }
1604
1605#undef FS_CASES
1606#undef FS_SUBCASES
1607
1608 UNPROTECT(5);
1609 return to;
1610}
1611
1612/* forceSymmetric(<[CRT]sparseMatrix>, uplo) */
1613SEXP R_sparse_force_symmetric(SEXP from, SEXP uplo)
1614{
1615 static const char *valid[] = {
1616 VALID_CSPARSE, VALID_RSPARSE, VALID_TSPARSE, "" };
1617 int ivalid = R_check_class_etc(from, valid);
1618 if (ivalid < 0)
1619 ERROR_INVALID_CLASS(from, __func__);
1620
1621 char ul = '\0';
1622 if (uplo != R_NilValue) {
1623 if (TYPEOF(uplo) != STRSXP || LENGTH(uplo) < 1 ||
1624 (uplo = STRING_ELT(uplo, 0)) == NA_STRING ||
1625 ((ul = *CHAR(uplo)) != 'U' && ul != 'L'))
1626 error(_("invalid '%s' to '%s'"), "uplo", __func__);
1627 }
1628
1629 return sparse_force_symmetric(from, valid[ivalid], ul);
1630}
1631
1632SEXP sparse_symmpart(SEXP from, const char *class)
1633{
1634 if (class[0] != 'z' && class[0] != 'd') {
1635 /* defined in ./coerce.c : */
1636 SEXP sparse_as_kind(SEXP, const char *, char);
1637 from = sparse_as_kind(from, class, 'd');
1638 }
1639 if (class[0] != 'z' && class[1] == 's')
1640 return from;
1641
1642 PROTECT_INDEX pid;
1643 PROTECT_WITH_INDEX(from, &pid);
1644
1645 char cl[] = ".s.Matrix";
1646 cl[0] = (class[0] != 'z') ? 'd' : 'z';
1647 cl[2] = class[2];
1648 SEXP to = PROTECT(newObject(cl));
1649
1650 SEXP dim = PROTECT(GET_SLOT(from, Matrix_DimSym));
1651 int *pdim = INTEGER(dim), n = pdim[0];
1652 if (pdim[1] != n)
1653 error(_("attempt to get symmetric part of non-square matrix"));
1654 if (n > 0)
1655 SET_SLOT(to, Matrix_DimSym, dim);
1656 UNPROTECT(1); /* dim */
1657
1658 SEXP dimnames = PROTECT(GET_SLOT(from, Matrix_DimNamesSym));
1659 if (class[1] == 's')
1660 SET_SLOT(to, Matrix_DimNamesSym, dimnames);
1661 else
1662 set_symmetrized_DimNames(to, dimnames, -1);
1663 UNPROTECT(1); /* dimnames */
1664
1665 char ul = 'U', di = 'N';
1666 if (class[1] != 'g') {
1667 SEXP uplo = PROTECT(GET_SLOT(from, Matrix_uploSym));
1668 ul = *CHAR(STRING_ELT(uplo, 0));
1669 if (ul != 'U')
1670 SET_SLOT(to, Matrix_uploSym, uplo);
1671 UNPROTECT(1); /* uplo */
1672 if (class[1] == 't') {
1673 SEXP diag = PROTECT(GET_SLOT(from, Matrix_diagSym));
1674 di = *CHAR(STRING_ELT(diag, 0));
1675 UNPROTECT(1); /* diag */
1676 }
1677 } else if (class[2] == 'R') {
1678 SEXP uplo = PROTECT(mkString("L"));
1679 ul = 'L';
1680 SET_SLOT(to, Matrix_uploSym, uplo);
1681 UNPROTECT(1); /* uplo */
1682 }
1683
1684 if (class[2] != 'T') {
1685
1686 SEXP iSym = (class[2] == 'C') ? Matrix_iSym : Matrix_jSym,
1687 p0 = PROTECT(GET_SLOT(from, Matrix_pSym)),
1688 i0 = PROTECT(GET_SLOT(from, iSym)),
1689 x0 = PROTECT(GET_SLOT(from, Matrix_xSym));
1690 int j, k, kend, *pp0 = INTEGER(p0) + 1, *pi0 = INTEGER(i0),
1691 nnz = pp0[n - 1];
1692
1693 if (class[1] == 'g') {
1694
1695 cl[1] = 'g';
1696 REPROTECT(from = sparse_transpose(from, cl, 0), pid);
1697
1698 SEXP p1 = PROTECT(allocVector(INTSXP, (R_xlen_t) n + 1)),
1699 p0_ = PROTECT(GET_SLOT(from, Matrix_pSym)),
1700 i0_ = PROTECT(GET_SLOT(from, iSym)),
1701 x0_ = PROTECT(GET_SLOT(from, Matrix_xSym));
1702 int k_, kend_, *pp0_ = INTEGER(p0_) + 1, *pi0_ = INTEGER(i0_),
1703 *pp1 = INTEGER(p1);
1704 *(pp1++) = 0;
1705
1706 for (j = 0, k = 0, k_ = 0; j < n; ++j) {
1707 kend = pp0[j];
1708 kend_ = pp0_[j];
1709 pp1[j] = pp1[j - 1];
1710 while (k < kend) {
1711 if (pi0[k] > j)
1712 k = kend;
1713 else {
1714 while (k_ < kend_ && pi0_[k_] < pi0[k]) {
1715 ++pp1[j];
1716 ++k_;
1717 }
1718 ++pp1[j];
1719 if (k_ < kend_ && pi0_[k_] == pi0[k])
1720 ++k_;
1721 ++k;
1722 }
1723 }
1724 while (k_ < kend_) {
1725 if (pi0_[k_] > j)
1726 k_ = kend_;
1727 else {
1728 ++pp1[j];
1729 ++k_;
1730 }
1731 }
1732 }
1733
1734 SEXP i1 = PROTECT(allocVector(INTSXP, pp1[n - 1])),
1735 x1 = PROTECT(allocVector(kindToType(cl[0]), pp1[n - 1]));
1736 int *pi1 = INTEGER(i1);
1737
1738#undef SP_LOOP
1739#define SP_LOOP(_CTYPE_, _PTR_, _ASSIGN_, _INCREMENT_) \
1740 do { \
1741 _CTYPE_ *px0 = _PTR_(x0), *px1 = _PTR_(x1), \
1742 *px0_ = _PTR_(x0_); \
1743 for (j = 0, k = 0, k_ = 0; j < n; ++j) { \
1744 kend = pp0[j]; \
1745 kend_ = pp0_[j]; \
1746 while (k < kend) { \
1747 if (pi0[k] > j) \
1748 k = kend; \
1749 else { \
1750 while (k_ < kend_ && pi0_[k_] < pi0[k]) { \
1751 *pi1 = pi0_[k_]; \
1752 _ASSIGN_((*px1), 0.5 * px0_[k_]); \
1753 ++k_; ++pi1; ++px1; \
1754 } \
1755 *pi1 = pi0[k]; \
1756 _ASSIGN_((*px1), 0.5 * px0[k]); \
1757 if (k_ < kend_ && pi0_[k_] == pi0[k]) { \
1758 _INCREMENT_((*px1), 0.5 * px0_[k_]); \
1759 ++k_; \
1760 } \
1761 ++k; ++pi1; ++px1; \
1762 } \
1763 } \
1764 while (k_ < kend_) { \
1765 if (pi0_[k_] > j) \
1766 k_ = kend_; \
1767 else { \
1768 *pi1 = pi0_[k_]; \
1769 _ASSIGN_((*px1), 0.5 * px0_[k_]); \
1770 ++k_; ++pi1; ++px1; \
1771 } \
1772 } \
1773 } \
1774 } while (0)
1775
1776 if (cl[0] == 'd')
1777 SP_LOOP(double, REAL, ASSIGN_REAL, INCREMENT_REAL);
1778 else
1779 SP_LOOP(Rcomplex, COMPLEX, ASSIGN_COMPLEX, INCREMENT_COMPLEX);
1780
1781 SET_SLOT(to, Matrix_pSym, p1);
1782 SET_SLOT(to, iSym, i1);
1783 SET_SLOT(to, Matrix_xSym, x1);
1784 UNPROTECT(6); /* x1, i1, p1, x0_, i0_, p0_ */
1785
1786 } else if (class[1] == 't') {
1787
1788 int leading = (class[2] == 'C') == (ul != 'U');
1789
1790 if (di == 'N') {
1791
1792 SEXP x1 = PROTECT(allocVector(kindToType(cl[0]), nnz));
1793
1794#undef SP_LOOP
1795#define SP_LOOP(_CTYPE_, _PTR_, _ASSIGN_) \
1796 do { \
1797 _CTYPE_ *px0 = _PTR_(x0), *px1 = _PTR_(x1); \
1798 if (leading) { \
1799 for (j = 0, k = 0; j < n; ++j) { \
1800 kend = pp0[j]; \
1801 if (k < kend) { \
1802 if (pi0[k] == j) \
1803 *px1 = *px0; \
1804 else \
1805 _ASSIGN_((*px1), 0.5 * (*px0)); \
1806 ++k, ++px0; ++px1; \
1807 while (k < kend) { \
1808 _ASSIGN_((*px1), 0.5 * (*px0)); \
1809 ++k; ++px0; ++px1; \
1810 } \
1811 } \
1812 } \
1813 } else { \
1814 for (j = 0, k = 0; j < n; ++j) { \
1815 kend = pp0[j]; \
1816 if (k < kend) { \
1817 while (k < kend - 1) { \
1818 _ASSIGN_((*px1), 0.5 * (*px0)); \
1819 ++k; ++px0; ++px1; \
1820 } \
1821 if (pi0[k] == j) \
1822 *px1 = *px0; \
1823 else \
1824 _ASSIGN_((*px1), 0.5 * (*px0)); \
1825 ++k; ++px0; ++px1; \
1826 } \
1827 } \
1828 } \
1829 } while (0)
1830
1831 if (cl[0] == 'd')
1832 SP_LOOP(double, REAL, ASSIGN_REAL);
1833 else
1834 SP_LOOP(Rcomplex, COMPLEX, ASSIGN_COMPLEX);
1835
1836 SET_SLOT(to, Matrix_pSym, p0);
1837 SET_SLOT(to, iSym, i0);
1838 SET_SLOT(to, Matrix_xSym, x1);
1839 UNPROTECT(1); /* x1 */
1840
1841 } else {
1842
1843 nnz += n;
1844 SEXP p1 = PROTECT(allocVector(INTSXP, (R_xlen_t) n + 1)),
1845 i1 = PROTECT(allocVector(INTSXP, nnz)),
1846 x1 = PROTECT(allocVector(kindToType(cl[0]), nnz));
1847 int *pp1 = INTEGER(p1), *pi1 = INTEGER(i1);
1848 *(pp1++) = 0;
1849
1850#undef SP_LOOP
1851#define SP_LOOP(_CTYPE_, _PTR_, _ASSIGN_, _ONE_) \
1852 do { \
1853 _CTYPE_ *px0 = _PTR_(x0), *px1 = _PTR_(x1); \
1854 if (leading) { \
1855 for (j = 0, k = 0; j < n; ++j) { \
1856 kend = pp0[j]; \
1857 pp1[j] = pp1[j - 1] + kend - k + 1; \
1858 *pi1 = j; \
1859 _ASSIGN_((*px1), _ONE_); \
1860 ++pi1, ++px1; \
1861 while (k < kend) { \
1862 *pi1 = *pi0; \
1863 _ASSIGN_((*px1), 0.5 * (*px0)); \
1864 ++k; ++pi0; ++pi1; ++px0; ++px1; \
1865 } \
1866 } \
1867 } else { \
1868 for (j = 0, k = 0; j < n; ++j) { \
1869 kend = pp0[j]; \
1870 pp1[j] = pp1[j - 1] + kend - k + 1; \
1871 while (k < kend) { \
1872 *pi1 = *pi0; \
1873 _ASSIGN_((*px1), 0.5 * (*px0)); \
1874 ++k; ++pi0; ++pi1; ++px0; ++px1; \
1875 } \
1876 *pi1 = j; \
1877 _ASSIGN_((*px1), _ONE_); \
1878 ++pi1; ++px1; \
1879 } \
1880 } \
1881 } while (0)
1882
1883 if (cl[0] == 'd')
1884 SP_LOOP(double, REAL, ASSIGN_REAL, 1.0);
1885 else
1886 SP_LOOP(Rcomplex, COMPLEX, ASSIGN_COMPLEX, Matrix_zone);
1887
1888 SET_SLOT(to, Matrix_pSym, p1);
1889 SET_SLOT(to, iSym, i1);
1890 SET_SLOT(to, Matrix_xSym, x1);
1891 UNPROTECT(3); /* p1, i1, x1 */
1892
1893 }
1894
1895 } else {
1896
1897 SET_SLOT(to, Matrix_pSym, p0);
1898 SET_SLOT(to, iSym, i0);
1899
1900 if (cl[0] == 'd')
1901 SET_SLOT(to, Matrix_xSym, x0);
1902 else {
1903 /* Symmetric part of Hermitian matrix is real part */
1904 SEXP x1 = PROTECT(duplicate(x0));
1905 zeroIm(x1);
1906 SET_SLOT(to, Matrix_xSym, x1);
1907 UNPROTECT(1); /* x1 */
1908 }
1909
1910 }
1911
1912 UNPROTECT(3); /* x0, i0, p0 */
1913
1914 } else {
1915
1916 SEXP i0 = PROTECT(GET_SLOT(from, Matrix_iSym)),
1917 j0 = PROTECT(GET_SLOT(from, Matrix_jSym)),
1918 x0 = PROTECT(GET_SLOT(from, Matrix_xSym));
1919 int *pi0 = INTEGER(i0), *pj0 = INTEGER(j0);
1920 R_xlen_t k, nnz = XLENGTH(i0);
1921
1922 if (class[1] == 'g') {
1923
1924 SEXP i1 = PROTECT(allocVector(INTSXP, nnz)),
1925 j1 = PROTECT(allocVector(INTSXP, nnz)),
1926 x1 = PROTECT(allocVector(kindToType(cl[0]), nnz));
1927 int *pi1 = INTEGER(i1), *pj1 = INTEGER(j1);
1928
1929#undef SP_LOOP
1930#define SP_LOOP(_CTYPE_, _PTR_, _ASSIGN_) \
1931 do { \
1932 _CTYPE_ *px0 = _PTR_(x0), *px1 = _PTR_(x1); \
1933 for (k = 0; k < nnz; ++k) { \
1934 if (*pi0 == *pj0) { \
1935 *pi1 = *pi0; \
1936 *pj1 = *pj0; \
1937 *px1 = *px0; \
1938 } else if (*pi0 < *pj0) { \
1939 *pi1 = *pi0; \
1940 *pj1 = *pj0; \
1941 _ASSIGN_((*px1), 0.5 * (*px0)); \
1942 } else { \
1943 *pi1 = *pj0; \
1944 *pj1 = *pi0; \
1945 _ASSIGN_((*px1), 0.5 * (*px0)); \
1946 } \
1947 ++pi0; ++pi1; ++pj0; ++pj1; ++px0; ++px1; \
1948 } \
1949 } while (0)
1950
1951 if (cl[0] == 'd')
1952 SP_LOOP(double, REAL, ASSIGN_REAL);
1953 else
1954 SP_LOOP(Rcomplex, COMPLEX, ASSIGN_COMPLEX);
1955
1956 SET_SLOT(to, Matrix_iSym, i1);
1957 SET_SLOT(to, Matrix_jSym, j1);
1958 SET_SLOT(to, Matrix_xSym, x1);
1959 UNPROTECT(3); /* x1, j1, i1 */
1960
1961 } else if (class[1] == 't') {
1962
1963 if (di == 'N') {
1964
1965 SEXP x1 = PROTECT(allocVector(kindToType(cl[0]), nnz));
1966
1967#undef SP_LOOP
1968#define SP_LOOP(_CTYPE_, _PTR_, _ASSIGN_) \
1969 do { \
1970 _CTYPE_ *px0 = _PTR_(x0), *px1 = _PTR_(x1); \
1971 for (k = 0; k < nnz; ++k) { \
1972 if (*pi0 == *pj0) \
1973 *px1 = *px0; \
1974 else \
1975 _ASSIGN_((*px1), 0.5 * (*px0)); \
1976 ++px0; ++px1; \
1977 } \
1978 } while (0)
1979
1980 if (cl[0] == 'd')
1981 SP_LOOP(double, REAL, ASSIGN_REAL);
1982 else
1983 SP_LOOP(Rcomplex, COMPLEX, ASSIGN_COMPLEX);
1984
1985 SET_SLOT(to, Matrix_iSym, i0);
1986 SET_SLOT(to, Matrix_jSym, j0);
1987 SET_SLOT(to, Matrix_xSym, x1);
1988 UNPROTECT(1); /* x1 */
1989
1990 } else {
1991
1992 SEXP i1 = PROTECT(allocVector(INTSXP, nnz + n)),
1993 j1 = PROTECT(allocVector(INTSXP, nnz + n)),
1994 x1 = PROTECT(allocVector(kindToType(cl[0]), nnz + n));
1995 int j, *pi1 = INTEGER(i1), *pj1 = INTEGER(j1);
1996
1997#undef SP_LOOP
1998#define SP_LOOP(_CTYPE_, _PTR_, _ASSIGN_, _ONE_) \
1999 do { \
2000 _CTYPE_ *px0 = _PTR_(x0), *px1 = _PTR_(x1); \
2001 for (k = 0; k < nnz; ++k) { \
2002 *pi1 = *pi0; \
2003 *pj1 = *pj0; \
2004 _ASSIGN_((*px1), 0.5 * (*px0)); \
2005 ++pi0; ++pi1; ++pj0; ++pj1; ++px0; ++px1; \
2006 } \
2007 for (j = 0; j < n; ++j) { \
2008 *pi1 = *pj1 = j; \
2009 _ASSIGN_((*px1), _ONE_); \
2010 ++pi1; ++pj1; ++px1; \
2011 } \
2012 } while (0)
2013
2014 if (cl[0] == 'd')
2015 SP_LOOP(double, REAL, ASSIGN_REAL, 1.0);
2016 else
2017 SP_LOOP(Rcomplex, COMPLEX, ASSIGN_COMPLEX, Matrix_zone);
2018
2019 SET_SLOT(to, Matrix_iSym, i1);
2020 SET_SLOT(to, Matrix_jSym, j1);
2021 SET_SLOT(to, Matrix_xSym, x1);
2022 UNPROTECT(3); /* x1, j1, i1 */
2023
2024 }
2025
2026 } else {
2027
2028 SET_SLOT(to, Matrix_iSym, i0);
2029 SET_SLOT(to, Matrix_jSym, j0);
2030
2031 if (cl[0] == 'd')
2032 SET_SLOT(to, Matrix_xSym, x0);
2033 else {
2034 /* Symmetric part of Hermitian matrix is real part */
2035 SEXP x1 = PROTECT(duplicate(x0));
2036 zeroIm(x1);
2037 SET_SLOT(to, Matrix_xSym, x1);
2038 UNPROTECT(1); /* x1 */
2039 }
2040
2041 }
2042
2043 UNPROTECT(3); /* x0, j0, i1 */
2044
2045 }
2046
2047 UNPROTECT(2); /* to, from */
2048 return to;
2049}
2050
2051/* symmpart(<[CRT]sparseMatrix>) */
2052SEXP R_sparse_symmpart(SEXP from)
2053{
2054 static const char *valid[] = {
2055 VALID_CSPARSE, VALID_RSPARSE, VALID_TSPARSE, "" };
2056 int ivalid = R_check_class_etc(from, valid);
2057 if (ivalid < 0)
2058 ERROR_INVALID_CLASS(from, __func__);
2059
2060 return sparse_symmpart(from, valid[ivalid]);
2061}
2062
2063SEXP sparse_skewpart(SEXP from, const char *class)
2064{
2065 if (class[0] != 'z' && class[0] != 'd') {
2066 /* defined in ./coerce.c : */
2067 SEXP sparse_as_kind(SEXP, const char *, char);
2068 from = sparse_as_kind(from, class, 'd');
2069 }
2070
2071 PROTECT_INDEX pid;
2072 PROTECT_WITH_INDEX(from, &pid);
2073
2074 char cl[] = "...Matrix";
2075 cl[0] = (class[0] != 'z') ? 'd' : 'z';
2076 cl[1] = (class[1] != 's') ? 'g' : 's';
2077 cl[2] = class[2];
2078 SEXP to = PROTECT(newObject(cl));
2079
2080 SEXP dim = PROTECT(GET_SLOT(from, Matrix_DimSym));
2081 int *pdim = INTEGER(dim), n = pdim[0];
2082 if (pdim[1] != n)
2083 error(_("attempt to get skew-symmetric part of non-square matrix"));
2084 if (n > 0)
2085 SET_SLOT(to, Matrix_DimSym, dim);
2086 UNPROTECT(1); /* dim */
2087
2088 SEXP dimnames = PROTECT(GET_SLOT(from, Matrix_DimNamesSym));
2089 if (class[1] == 's')
2090 SET_SLOT(to, Matrix_DimNamesSym, dimnames);
2091 else
2092 set_symmetrized_DimNames(to, dimnames, -1);
2093 UNPROTECT(1); /* dimnames */
2094
2095 if (class[1] == 's') {
2096
2097 SEXP uplo = PROTECT(GET_SLOT(from, Matrix_uploSym));
2098 char ul = *CHAR(STRING_ELT(uplo, 0));
2099 if (ul != 'U')
2100 SET_SLOT(to, Matrix_uploSym, uplo);
2101 UNPROTECT(1); /* uplo */
2102
2103 /* Skew-symmetric part of Hermitian matrix is imaginary part */
2104 if (class[0] != 'z') {
2105 if (class[2] != 'T') {
2106 SEXP p1 = PROTECT(allocVector(INTSXP, (R_xlen_t) n + 1));
2107 int *pp1 = INTEGER(p1);
2108 Matrix_memset(pp1, 0, (R_xlen_t) n + 1, sizeof(int));
2109 SET_SLOT(to, Matrix_pSym, p1);
2110 UNPROTECT(1); /* p1 */
2111 }
2112 } else {
2113 if (class[2] != 'T') {
2114 SEXP p0 = PROTECT(GET_SLOT(from, Matrix_pSym));
2115 SET_SLOT(to, Matrix_pSym, p0);
2116 UNPROTECT(1); /* p0 */
2117 }
2118 if (class[2] != 'R') {
2119 SEXP i0 = PROTECT(GET_SLOT(from, Matrix_iSym));
2120 SET_SLOT(to, Matrix_iSym, i0);
2121 UNPROTECT(1); /* i0 */
2122 }
2123 if (class[2] != 'C') {
2124 SEXP j0 = PROTECT(GET_SLOT(from, Matrix_jSym));
2125 SET_SLOT(to, Matrix_jSym, j0);
2126 UNPROTECT(1); /* j0 */
2127 }
2128 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym)),
2129 x1 = PROTECT(duplicate(x0));
2130 zeroRe(x1);
2131 SET_SLOT(to, Matrix_xSym, x1);
2132 UNPROTECT(2); /* x1, x0 */
2133 }
2134
2135 } else if (class[2] != 'T') {
2136
2137 SEXP iSym = (class[2] == 'C') ? Matrix_iSym : Matrix_jSym,
2138 p0 = PROTECT(GET_SLOT(from, Matrix_pSym)),
2139 i0 = PROTECT(GET_SLOT(from, iSym)),
2140 x0 = PROTECT(GET_SLOT(from, Matrix_xSym)),
2141 p1 = PROTECT(allocVector(INTSXP, (R_xlen_t) n + 1));
2142 int j, k, kend, k_, kend_, *pp0 = INTEGER(p0), *pi0 = INTEGER(i0),
2143 *pp1 = INTEGER(p1);
2144 pp0++; *(pp1++) = 0;
2145
2146 REPROTECT(from = sparse_transpose(from, cl, 0), pid);
2147
2148 SEXP
2149 p0_ = PROTECT(GET_SLOT(from, Matrix_pSym)),
2150 i0_ = PROTECT(GET_SLOT(from, iSym)),
2151 x0_ = PROTECT(GET_SLOT(from, Matrix_xSym));
2152 int *pp0_ = INTEGER(p0_), *pi0_ = INTEGER(i0_), *pp1_;
2153 pp0_++;
2154 Matrix_Calloc(pp1_, n, int);
2155
2156 for (j = 0, k = 0, k_ = 0; j < n; ++j) {
2157 kend = pp0[j];
2158 kend_ = pp0_[j];
2159 while (k < kend) {
2160 if (pi0[k] >= j)
2161 k = kend;
2162 else {
2163 while (k_ < kend_ && pi0_[k_] < pi0[k]) {
2164 ++pp1_[j];
2165 ++pp1_[pi0_[k_]];
2166 ++k_;
2167 }
2168 ++pp1_[j];
2169 ++pp1_[pi0[k]];
2170 if (k_ < kend_ && pi0_[k_] == pi0[k])
2171 ++k_;
2172 ++k;
2173 }
2174 }
2175 while (k_ < kend_) {
2176 if (pi0_[k_] >= j)
2177 k_ = kend_;
2178 else {
2179 ++pp1_[j];
2180 ++pp1_[pi0_[k_]];
2181 ++k_;
2182 }
2183 }
2184 }
2185
2186 for (j = 0; j < n; ++j) {
2187 pp1[j] = pp1[j - 1] + pp1_[j];
2188 pp1_[j] = pp1[j - 1];
2189 }
2190
2191 SEXP i1 = PROTECT(allocVector(INTSXP, pp1[n - 1])),
2192 x1 = PROTECT(allocVector(kindToType(cl[0]), pp1[n - 1]));
2193 int *pi1 = INTEGER(i1);
2194
2195#undef SP_LOOP
2196#define SP_LOOP(_CTYPE_, _PTR_, _ASSIGN_, _INCREMENT_) \
2197 do { \
2198 _CTYPE_ *px0 = _PTR_(x0), *px1 = _PTR_(x1), \
2199 *px0_ = _PTR_(x0_); \
2200 for (j = 0, k = 0, k_ = 0; j < n; ++j) { \
2201 kend = pp0[j]; \
2202 kend_ = pp0_[j]; \
2203 while (k < kend) { \
2204 if (pi0[k] >= j) \
2205 k = kend; \
2206 else { \
2207 while (k_ < kend_ && pi0_[k_] < pi0[k]) { \
2208 pi1[pp1_[j]] = pi0_[k_]; \
2209 _ASSIGN_(px1[pp1_[j]], -0.5 * px0_[k_]); \
2210 pi1[pp1_[pi0_[k_]]] = j; \
2211 _ASSIGN_(px1[pp1_[pi0_[k_]]], -px1[pp1_[j]]); \
2212 ++pp1_[j]; \
2213 ++pp1_[pi0_[k_]]; \
2214 ++k_; \
2215 } \
2216 pi1[pp1_[j]] = pi0[k]; \
2217 _ASSIGN_(px1[pp1_[j]], 0.5 * px0[k]); \
2218 if (k_ < kend_ && pi0_[k_] == pi0[k]) { \
2219 _INCREMENT_(px1[pp1_[j]], -0.5 * px0_[k_]); \
2220 ++k_; \
2221 } \
2222 pi1[pp1_[pi0[k]]] = j; \
2223 _ASSIGN_(px1[pp1_[pi0[k]]], -px1[pp1_[j]]); \
2224 ++pp1_[j]; \
2225 ++pp1_[pi0[k]]; \
2226 ++k; \
2227 } \
2228 } \
2229 while (k_ < kend_) { \
2230 if (pi0_[k_] >= j) \
2231 k_ = kend_; \
2232 else { \
2233 pi1[pp1_[j]] = pi0_[k_]; \
2234 _ASSIGN_(px1[pp1_[j]], -0.5 * px0_[k_]); \
2235 pi1[pp1_[pi0_[k_]]] = j; \
2236 _ASSIGN_(px1[pp1_[pi0_[k_]]], -px1[pp1_[j]]); \
2237 ++pp1_[j]; \
2238 ++pp1_[pi0_[k_]]; \
2239 ++k_; \
2240 } \
2241 } \
2242 } \
2243 } while (0)
2244
2245 if (cl[0] == 'd')
2246 SP_LOOP(double, REAL, ASSIGN_REAL, INCREMENT_REAL);
2247 else
2248 SP_LOOP(Rcomplex, COMPLEX, ASSIGN_COMPLEX, INCREMENT_COMPLEX);
2249
2250 Matrix_Free(pp1_, n);
2251 SET_SLOT(to, Matrix_pSym, p1);
2252 SET_SLOT(to, iSym, i1);
2253 SET_SLOT(to, Matrix_xSym, x1);
2254 UNPROTECT(9); /* x1, i1, p1, x0, i0, p0, x0_, i0_, p0_ */
2255
2256 } else {
2257
2258 SEXP i0 = PROTECT(GET_SLOT(from, Matrix_iSym)),
2259 j0 = PROTECT(GET_SLOT(from, Matrix_jSym)),
2260 x0 = PROTECT(GET_SLOT(from, Matrix_xSym));
2261 int *pi0 = INTEGER(i0), *pj0 = INTEGER(j0);
2262 R_xlen_t k, nnz0 = XLENGTH(i0), nnz1 = nnz0;
2263
2264 for (k = 0; k < nnz0; ++k)
2265 if (pi0[k] == pj0[k])
2266 --nnz1;
2267 nnz1 *= 2;
2268
2269 SEXP i1 = PROTECT(allocVector(INTSXP, nnz1)),
2270 j1 = PROTECT(allocVector(INTSXP, nnz1)),
2271 x1 = PROTECT(allocVector(kindToType(cl[0]), nnz1));
2272 int *pi1 = INTEGER(i1), *pj1 = INTEGER(j1);
2273
2274#undef SP_LOOP
2275#define SP_LOOP(_CTYPE_, _PTR_, _ASSIGN_) \
2276 do { \
2277 _CTYPE_ *px0 = _PTR_(x0), *px1 = _PTR_(x1); \
2278 for (k = 0; k < nnz0; ++k) { \
2279 if (*pi0 != *pj0) { \
2280 *pi1 = *pi0; \
2281 *pj1 = *pj0; \
2282 _ASSIGN_((*px1), 0.5 * (*px0)); \
2283 ++pi1; ++pj1; ++px1; \
2284 *pi1 = *pj0; \
2285 *pj1 = *pi0; \
2286 _ASSIGN_((*px1), -0.5 * (*px0)); \
2287 ++pi1; ++pj1; ++px1; \
2288 } \
2289 ++pi0; ++pj0; ++px0; \
2290 } \
2291 } while (0)
2292
2293 if (cl[0] == 'd')
2294 SP_LOOP(double, REAL, ASSIGN_REAL);
2295 else
2296 SP_LOOP(Rcomplex, COMPLEX, ASSIGN_COMPLEX);
2297
2298 SET_SLOT(to, Matrix_iSym, i1);
2299 SET_SLOT(to, Matrix_jSym, j1);
2300 SET_SLOT(to, Matrix_xSym, x1);
2301 UNPROTECT(6); /* x1, j1, i1, x0, j0, i0 */
2302
2303 }
2304
2305#undef SP_LOOP
2306
2307 UNPROTECT(2); /* to, from */
2308 return to;
2309}
2310
2311/* skewpart(<[CRT]sparseMatrix>) */
2312SEXP R_sparse_skewpart(SEXP from)
2313{
2314 static const char *valid[] = {
2315 VALID_CSPARSE, VALID_RSPARSE, VALID_TSPARSE, "" };
2316 int ivalid = R_check_class_etc(from, valid);
2317 if (ivalid < 0)
2318 ERROR_INVALID_CLASS(from, __func__);
2319
2320 return sparse_skewpart(from, valid[ivalid]);
2321}
2322
2323int sparse_is_symmetric(SEXP obj, const char *class, int checkDN)
2324{
2325 if (class[1] == 's')
2326 return 1;
2327
2328 if (checkDN) {
2329 SEXP dimnames = GET_SLOT(obj, Matrix_DimNamesSym);
2330 if (!DimNames_is_symmetric(dimnames))
2331 return 0;
2332 }
2333
2334 if (class[1] == 't')
2335 return sparse_is_diagonal(obj, class);
2336
2337 SEXP dim = GET_SLOT(obj, Matrix_DimSym);
2338 int *pdim = INTEGER(dim), n = pdim[0];
2339 if (pdim[1] != n)
2340 return 0;
2341 if (n <= 1)
2342 return 1;
2343
2344 if (class[2] == 'T') {
2345 /* defined in ./coerce.c : */
2346 SEXP sparse_as_Csparse(SEXP, const char *);
2347 obj = sparse_as_Csparse(obj, class);
2348 }
2349 PROTECT(obj);
2350
2351 SEXP iSym = (class[2] != 'R') ? Matrix_iSym : Matrix_jSym,
2352 p0 = PROTECT(GET_SLOT(obj, Matrix_pSym)),
2353 i0 = PROTECT(GET_SLOT(obj, iSym));
2354 int i, j, k, kend, *pp_, *pp0 = INTEGER(p0) + 1, *pi0 = INTEGER(i0);
2355 Matrix_Calloc(pp_, n, int);
2356 Matrix_memcpy(pp_, pp0 - 1, n, sizeof(int));
2357
2358 int ans = 0;
2359
2360#define IS_LOOP(_CTYPE_, _PTR_, _MASK_, _NOTREAL_, _NOTCONJ_) \
2361 do { \
2362 _MASK_(_CTYPE_ *px0 = _PTR_(x0)); \
2363 for (j = 0, k = 0; j < n; ++j) { \
2364 kend = pp0[j]; \
2365 while (k < kend) { \
2366 i = pi0[k]; \
2367 if (i >= j) { \
2368 if (i == j) { \
2369 if (_NOTREAL_(px0[k])) \
2370 goto finish; \
2371 ++pp_[j]; \
2372 } \
2373 k = kend; \
2374 } else { \
2375 if (pp_[i] == pp0[i] || pi0[pp_[i]] != j || \
2376 _NOTCONJ_(px0[k], px0[pp_[i]])) \
2377 goto finish; \
2378 ++pp_[i]; \
2379 ++pp_[j]; \
2380 ++k; \
2381 } \
2382 } \
2383 } \
2384 } while (0)
2385
2386#undef NOTCONJ_PATTERN
2387#define NOTCONJ_PATTERN(_X_, _Y_) 0
2388
2389 /* For all X[i,j], i >= j, we require:
2390 o that X[j,i] exists
2391 o that X[j,i] == Conj(X[i,j])
2392 o that Im(X[j,i]) == 0 if i == j
2393 */
2394 if (class[0] == 'n')
2395 IS_LOOP(int, LOGICAL, HIDE, NOTREAL_PATTERN, NOTCONJ_PATTERN);
2396 else {
2397 SEXP x0 = GET_SLOT(obj, Matrix_xSym);
2398 switch (class[0]) {
2399 case 'l':
2400 IS_LOOP(int, LOGICAL, SHOW, NOTREAL_LOGICAL, NOTCONJ_LOGICAL);
2401 break;
2402 case 'i':
2403 IS_LOOP(int, INTEGER, SHOW, NOTREAL_INTEGER, NOTCONJ_INTEGER);
2404 break;
2405 case 'd':
2406 IS_LOOP(double, REAL, SHOW, NOTREAL_REAL, NOTCONJ_REAL);
2407 break;
2408 case 'z':
2409 IS_LOOP(Rcomplex, COMPLEX, SHOW, NOTREAL_COMPLEX, NOTCONJ_COMPLEX);
2410 break;
2411 default:
2412 break;
2413 }
2414 }
2415
2416 /* We further require that the upper and lower triangles
2417 have the same number of entries ...
2418 */
2419 for (j = 0; j < n; ++j)
2420 if (pp_[j] != pp0[j])
2421 goto finish;
2422
2423 ans = 1;
2424
2425finish:
2426 Matrix_Free(pp_, n);
2427 UNPROTECT(3); /* i0, p0, obj */
2428
2429#undef IS_LOOP
2430
2431 return ans;
2432}
2433
2434/* isSymmetric(<[CRT]sparseMatrix>, checkDN, tol = 0)
2435 NB: requires symmetric nonzero pattern
2436 TODO: support 'tol', 'scale' arguments and bypass all.equal ??
2437*/
2438SEXP R_sparse_is_symmetric(SEXP obj, SEXP checkDN)
2439{
2440 static const char *valid[] = {
2441 VALID_CSPARSE, VALID_RSPARSE, VALID_TSPARSE, "" };
2442 int ivalid = R_check_class_etc(obj, valid);
2443 if (ivalid < 0)
2444 ERROR_INVALID_CLASS(obj, __func__);
2445
2446 int checkDN_;
2447 if (TYPEOF(checkDN) != LGLSXP || LENGTH(checkDN) < 1 ||
2448 (checkDN_ = LOGICAL(checkDN)[0]) == NA_LOGICAL)
2449 error(_("'%s' must be %s or %s"), "checkDN", "TRUE", "FALSE");
2450
2451 return ScalarLogical(sparse_is_symmetric(obj, valid[ivalid], checkDN_));
2452}
2453
2454int sparse_is_triangular(SEXP obj, const char *class, int upper)
2455{
2456 if (class[1] == 't') {
2457 SEXP uplo = GET_SLOT(obj, Matrix_uploSym);
2458 char ul = *CHAR(STRING_ELT(uplo, 0));
2459 if (upper == NA_LOGICAL || (upper != 0) == (ul == 'U'))
2460 return (ul == 'U') ? 1 : -1;
2461 else if (sparse_is_diagonal(obj, class))
2462 return (ul == 'U') ? -1 : 1;
2463 else
2464 return 0;
2465 }
2466
2467 if (class[1] == 's') {
2468 if (!sparse_is_diagonal(obj, class))
2469 return 0;
2470 SEXP uplo = GET_SLOT(obj, Matrix_uploSym);
2471 char ul = *CHAR(STRING_ELT(uplo, 0));
2472 if (upper == NA_LOGICAL)
2473 return (ul == 'U') ? 1 : -1;
2474 else
2475 return (upper != 0) ? 1 : -1;
2476 }
2477
2478 SEXP dim = GET_SLOT(obj, Matrix_DimSym);
2479 int *pdim = INTEGER(dim), n = pdim[0];
2480 if (pdim[1] != n)
2481 return 0;
2482 if (n <= 1)
2483 return (upper != 0) ? 1 : -1;
2484
2485 if (class[2] != 'T') {
2486
2487 SEXP iSym = (class[2] == 'C') ? Matrix_iSym : Matrix_jSym,
2488 p0 = PROTECT(GET_SLOT(obj, Matrix_pSym)),
2489 i0 = PROTECT(GET_SLOT(obj, iSym));
2490 UNPROTECT(2); /* i0, p0 */
2491 int j, k, kend, *pp0 = INTEGER(p0) + 1, *pi0 = INTEGER(i0);
2492
2493 if (upper == NA_LOGICAL) {
2494 /* Examine last entry in each "column" */
2495 for (j = 0, k = 0; j < n; ++j) {
2496 kend = pp0[j];
2497 if (k < kend && pi0[kend - 1] > j)
2498 break;
2499 k = kend;
2500 }
2501 if (j == n)
2502 return (class[2] == 'C') ? 1 : -1;
2503 /* Examine first entry in each "column" */
2504 for (j = 0, k = 0; j < n; ++j) {
2505 kend = pp0[j];
2506 if (k < kend && pi0[k] < j)
2507 break;
2508 k = kend;
2509 }
2510 if (j == n)
2511 return (class[2] == 'C') ? -1 : 1;
2512 return 0;
2513 } else if ((class[2] == 'C') == (upper != 0)) {
2514 /* Examine last entry in each "column" */
2515 for (j = 0, k = 0; j < n; ++j) {
2516 kend = pp0[j];
2517 if (k < kend && pi0[kend - 1] > j)
2518 return 0;
2519 k = kend;
2520 }
2521 return (class[2] == 'C') ? 1 : -1;
2522 } else {
2523 /* Examine first entry in each "column" */
2524 for (j = 0, k = 0; j < n; ++j) {
2525 kend = pp0[j];
2526 if (k < kend && pi0[k] < j)
2527 return 0;
2528 k = kend;
2529 }
2530 return (class[2] == 'C') ? -1 : 1;
2531 }
2532
2533 } else {
2534
2535 SEXP i0 = PROTECT(GET_SLOT(obj, Matrix_iSym)),
2536 j0 = PROTECT(GET_SLOT(obj, Matrix_jSym));
2537 UNPROTECT(2); /* i0, j0 */
2538 int *pi0 = INTEGER(i0), *pj0 = INTEGER(j0);
2539 R_xlen_t k, nnz0 = XLENGTH(i0);
2540
2541 if (upper == NA_LOGICAL) {
2542 for (k = 0; k < nnz0; ++k)
2543 if (pi0[k] > pj0[k])
2544 break;
2545 if (k == nnz0)
2546 return 1;
2547 for (k = 0; k < nnz0; ++k)
2548 if (pi0[k] < pj0[k])
2549 break;
2550 if (k == nnz0)
2551 return -1;
2552 return 0;
2553 } else if (upper != 0) {
2554 for (k = 0; k < nnz0; ++k)
2555 if (pi0[k] > pj0[k])
2556 return 0;
2557 return 1;
2558 } else {
2559 for (k = 0; k < nnz0; ++k)
2560 if (pi0[k] < pj0[k])
2561 return 0;
2562 return -1;
2563 }
2564
2565 }
2566}
2567
2568/* isTriangular(<[CRT]sparseMatrix>, upper)
2569 NB: requires triangular nonzero pattern
2570*/
2571SEXP R_sparse_is_triangular(SEXP obj, SEXP upper)
2572{
2573 static const char *valid[] = {
2574 VALID_CSPARSE, VALID_RSPARSE, VALID_TSPARSE, "" };
2575 int ivalid = R_check_class_etc(obj, valid);
2576 if (ivalid < 0)
2577 ERROR_INVALID_CLASS(obj, __func__);
2578
2579 if (TYPEOF(upper) != LGLSXP || LENGTH(upper) < 1)
2580 error(_("'%s' must be %s or %s or %s"), "upper", "TRUE", "FALSE", "NA");
2581 int upper_ = LOGICAL(upper)[0];
2582
2583 int ans_ = sparse_is_triangular(obj, valid[ivalid], upper_);
2584 SEXP ans = allocVector(LGLSXP, 1);
2585 LOGICAL(ans)[0] = ans_ != 0;
2586 if (upper_ == NA_LOGICAL && ans_ != 0) {
2587 PROTECT(ans);
2588 static
2589 SEXP kindSym = NULL;
2590 SEXP kindVal = PROTECT(mkString((ans_ > 0) ? "U" : "L"));
2591 if (!kindSym) kindSym = install("kind");
2592 setAttrib(ans, kindSym, kindVal);
2593 UNPROTECT(2);
2594 }
2595 return ans;
2596}
2597
2598int sparse_is_diagonal(SEXP obj, const char *class)
2599{
2600 SEXP dim = GET_SLOT(obj, Matrix_DimSym);
2601 int *pdim = INTEGER(dim), n = pdim[0];
2602 if (pdim[1] != n)
2603 return 0;
2604 if (n <= 1)
2605 return 1;
2606
2607 if (class[2] != 'T') {
2608
2609 SEXP iSym = (class[2] == 'C') ? Matrix_iSym : Matrix_jSym,
2610 p0 = PROTECT(GET_SLOT(obj, Matrix_pSym)),
2611 i0 = PROTECT(GET_SLOT(obj, iSym));
2612 UNPROTECT(2); /* i0, p0 */
2613 int j, k, kend, *pp0 = INTEGER(p0) + 1, *pi0 = INTEGER(i0);
2614
2615 for (j = 0, k = 0; j < n; ++j) {
2616 kend = pp0[j];
2617 if (kend - k > 1 || (kend - k == 1 && pi0[k] != j))
2618 return 0;
2619 k = kend;
2620 }
2621 return 1;
2622
2623 } else {
2624
2625 SEXP i0 = PROTECT(GET_SLOT(obj, Matrix_iSym)),
2626 j0 = PROTECT(GET_SLOT(obj, Matrix_jSym));
2627 UNPROTECT(2); /* i0, j0 */
2628 int *pi0 = INTEGER(i0), *pj0 = INTEGER(j0);
2629 R_xlen_t k, nnz0 = XLENGTH(i0);
2630
2631 for (k = 0; k < nnz0; ++k)
2632 if (*(pi0++) != *(pj0++))
2633 return 0;
2634 return 1;
2635
2636 }
2637}
2638
2639/* isDiagonal(<[CRT]sparseMatrix>)
2640 NB: requires diagonal nonzero pattern
2641*/
2642SEXP R_sparse_is_diagonal(SEXP obj)
2643{
2644 static const char *valid[] = {
2645 VALID_CSPARSE, VALID_RSPARSE, VALID_TSPARSE, "" };
2646 int ivalid = R_check_class_etc(obj, valid);
2647 if (ivalid < 0)
2648 ERROR_INVALID_CLASS(obj, __func__);
2649
2650 return ScalarLogical(sparse_is_diagonal(obj, valid[ivalid]));
2651}
2652
2653#define MAP(_I_) work[_I_]
2654#define NOMAP(_I_) _I_
2655
2656#define CAST_PATTERN(_X_) 1
2657#define CAST_LOGICAL(_X_) (_X_ != 0)
2658#define CAST_INTEGER(_X_) _X_
2659#define CAST_REAL(_X_) _X_
2660#define CAST_COMPLEX(_X_) _X_
2661
2662#define SUM_CASES(_MAP_) \
2663do { \
2664 if (class[0] == 'n') { \
2665 if (mean) \
2666 SUM_LOOP(int, LOGICAL, double, REAL, HIDE, \
2667 0.0, 1.0, NA_REAL, ISNA_PATTERN, \
2668 _MAP_, CAST_PATTERN, INCREMENT_REAL, SCALE2_REAL); \
2669 else \
2670 SUM_LOOP(int, LOGICAL, int, INTEGER, HIDE, \
2671 0, 1, NA_INTEGER, ISNA_PATTERN, \
2672 _MAP_, CAST_PATTERN, INCREMENT_INTEGER, SCALE2_REAL); \
2673 } else { \
2674 SEXP x0 = PROTECT(GET_SLOT(obj, Matrix_xSym)); \
2675 switch (class[0]) { \
2676 case 'l': \
2677 if (mean) \
2678 SUM_LOOP(int, LOGICAL, double, REAL, SHOW, \
2679 0.0, 1.0, NA_REAL, ISNA_LOGICAL, \
2680 _MAP_, CAST_LOGICAL, INCREMENT_REAL, SCALE2_REAL); \
2681 else \
2682 SUM_LOOP(int, LOGICAL, int, INTEGER, SHOW, \
2683 0, 1, NA_INTEGER, ISNA_LOGICAL, \
2684 _MAP_, CAST_LOGICAL, INCREMENT_INTEGER, SCALE2_REAL); \
2685 break; \
2686 case 'i': \
2687 SUM_LOOP(int, INTEGER, double, REAL, SHOW, \
2688 0.0, 1.0, NA_REAL, ISNA_INTEGER, \
2689 _MAP_, CAST_INTEGER, INCREMENT_REAL, SCALE2_REAL); \
2690 break; \
2691 case 'd': \
2692 SUM_LOOP(double, REAL, double, REAL, SHOW, \
2693 0.0, 1.0, NA_REAL, ISNA_REAL, \
2694 _MAP_, CAST_REAL, INCREMENT_REAL, SCALE2_REAL); \
2695 break; \
2696 case 'z': \
2697 SUM_LOOP(Rcomplex, COMPLEX, Rcomplex, COMPLEX, SHOW, \
2698 Matrix_zzero, Matrix_zone, Matrix_zna, ISNA_COMPLEX, \
2699 _MAP_, CAST_COMPLEX, INCREMENT_COMPLEX, SCALE2_COMPLEX); \
2700 break; \
2701 default: \
2702 break; \
2703 } \
2704 UNPROTECT(1); /* x0 */ \
2705 } \
2706} while (0)
2707
2708#define SUM_TYPEOF(c) (c == 'z') ? CPLXSXP : ((mean || c == 'd' || c == 'i') ? REALSXP : INTSXP)
2709
2710static
2711void Csparse_colsum(SEXP obj, const char *class,
2712 int m, int n, char di, int narm, int mean,
2713 SEXP res)
2714{
2715 int narm_ = narm && mean && class[0] != 'n';
2716
2717 SEXP p0 = PROTECT(GET_SLOT(obj, Matrix_pSym));
2718 int *pp0 = INTEGER(p0) + 1, j, k, kend, nnz1 = n, count = -1;
2719
2720 if (IS_S4_OBJECT(res)) {
2721
2722 if (di == 'N') {
2723 nnz1 = 0;
2724 for (j = 0; j < n; ++j)
2725 if (pp0[j - 1] < pp0[j])
2726 ++nnz1;
2727 }
2728
2729 SEXP
2730 j1 = PROTECT(allocVector(INTSXP, nnz1)),
2731 x1 = PROTECT(allocVector(SUM_TYPEOF(class[0]), nnz1));
2732 SET_SLOT(res, Matrix_iSym, j1);
2733 SET_SLOT(res, Matrix_xSym, x1);
2734
2735 int *pj1 = INTEGER(j1);
2736 if (di != 'N')
2737 for (j = 0; j < n; ++j)
2738 *(pj1++) = j + 1;
2739 else
2740 for (j = 0; j < n; ++j)
2741 if (pp0[j - 1] < pp0[j])
2742 *(pj1++) = j + 1;
2743
2744#define SUM_LOOP(_CTYPE0_, _PTR0_, _CTYPE1_, _PTR1_, _MASK_, \
2745 _ZERO_, _ONE_, _NA_, _ISNA_, \
2746 _MAP_, _CAST_, _INCREMENT_, _SCALE2_) \
2747 do { \
2748 _MASK_(_CTYPE0_ *px0 = _PTR0_(x0)); \
2749 _CTYPE1_ *px1 = _PTR1_(x1) , tmp; \
2750 for (j = 0, k = 0; j < n; ++j) { \
2751 kend = pp0[j]; \
2752 if (k < kend || nnz1 == n) { \
2753 *px1 = (di != 'N') ? _ONE_ : _ZERO_; \
2754 if (mean) \
2755 count = m; \
2756 while (k < kend) { \
2757 if (_ISNA_(*px0)) { \
2758 if (!narm) \
2759 *px1 = _NA_; \
2760 else if (narm_) \
2761 --count; \
2762 } else { \
2763 tmp = _CAST_(*px0); \
2764 _INCREMENT_((*px1), tmp); \
2765 } \
2766 _MASK_(++px0); \
2767 ++k; \
2768 } \
2769 if (mean) \
2770 _SCALE2_((*px1), count); \
2771 ++px1; \
2772 } \
2773 } \
2774 } while (0)
2775
2776 SUM_CASES(MAP);
2777 UNPROTECT(2); /* x1, j1 */
2778
2779 } else {
2780
2781 SEXP x1 = res;
2782 SUM_CASES(NOMAP);
2783
2784 }
2785
2786#undef SUM_LOOP
2787
2788 UNPROTECT(1); /* p0 */
2789 return;
2790}
2791
2792static
2793void Csparse_rowsum(SEXP obj, const char *class,
2794 int m, int n, char di, int narm, int mean,
2795 SEXP res, SEXP iSym)
2796{
2797 int narm_ = narm && mean && class[0] != 'n';
2798
2799 SEXP p0 = PROTECT(GET_SLOT(obj, Matrix_pSym)),
2800 i0 = PROTECT(GET_SLOT(obj, iSym));
2801 int *pp0 = INTEGER(p0) + 1, *pi0 = INTEGER(i0), i, j, k, kend,
2802 nnz0 = pp0[n - 1], nnz1 = m;
2803
2804 if (IS_S4_OBJECT(res)) {
2805
2806 int *work;
2807 Matrix_Calloc(work, m, int);
2808 if (di != 'N')
2809 for (i = 0; i < m; ++i)
2810 work[i] = i;
2811 else {
2812 if (class[1] != 's') {
2813 for (k = 0; k < nnz0; ++k)
2814 ++work[pi0[k]];
2815 } else {
2816 for (j = 0, k = 0; j < n; ++j) {
2817 kend = pp0[j];
2818 while (k < kend) {
2819 ++work[pi0[k]];
2820 if (pi0[k] != j)
2821 ++work[ j];
2822 ++k;
2823 }
2824 }
2825 }
2826 nnz1 = 0;
2827 for (i = 0; i < m; ++i)
2828 work[i] = (work[i]) ? nnz1++ : -1;
2829 }
2830
2831 SEXP
2832 i1 = PROTECT(allocVector(INTSXP, nnz1)),
2833 x1 = PROTECT(allocVector(SUM_TYPEOF(class[0]), nnz1));
2834 SET_SLOT(res, Matrix_iSym, i1);
2835 SET_SLOT(res, Matrix_xSym, x1);
2836 int *pi1 = INTEGER(i1);
2837 if (narm_)
2838 for (i = 0; i < nnz1; ++i)
2839 pi1[i] = n;
2840
2841#define SUM_LOOP(_CTYPE0_, _PTR0_, _CTYPE1_, _PTR1_, _MASK_, \
2842 _ZERO_, _ONE_, _NA_, _ISNA_, \
2843 _MAP_, _CAST_, _INCREMENT_, _SCALE2_) \
2844 do { \
2845 _MASK_(_CTYPE0_ *px0 = _PTR0_(x0)); \
2846 _CTYPE1_ *px1 = _PTR1_(x1) ; \
2847 _CTYPE1_ tmp = (di != 'N') ? _ONE_ : _ZERO_; \
2848 for (i = 0; i < nnz1; ++i) \
2849 px1[i] = tmp; \
2850 if (class[1] != 's') { \
2851 for (k = 0; k < nnz0; ++k) { \
2852 if (_ISNA_(px0[k])) { \
2853 if (!narm) \
2854 px1[_MAP_(pi0[k])] = _NA_; \
2855 else if (narm_) \
2856 --pi1[_MAP_(pi0[k])]; \
2857 } else { \
2858 tmp = _CAST_(px0[k]); \
2859 _INCREMENT_(px1[_MAP_(pi0[k])], tmp); \
2860 } \
2861 } \
2862 } else { \
2863 int off; \
2864 for (j = 0, k = 0; j < n; ++j) { \
2865 kend = pp0[j]; \
2866 while (k < kend) { \
2867 off = pi0[k] != j; \
2868 if (_ISNA_(px0[k])) { \
2869 if (!narm) { \
2870 px1[_MAP_(pi0[k])] = _NA_; \
2871 if (off) \
2872 px1[_MAP_( j)] = _NA_; \
2873 } else if (narm_) { \
2874 --pi1[_MAP_(pi0[k])]; \
2875 if (off) \
2876 --pi1[_MAP_( j)]; \
2877 } \
2878 } else { \
2879 tmp = _CAST_(px0[k]); \
2880 _INCREMENT_(px1[_MAP_(pi0[k])], tmp); \
2881 if (off) \
2882 _INCREMENT_(px1[_MAP_( j)], tmp); \
2883 } \
2884 ++k; \
2885 } \
2886 } \
2887 } \
2888 if (mean) { \
2889 if (narm_) \
2890 for (i = 0; i < nnz1; ++i) \
2891 _SCALE2_(px1[i], pi1[i]); \
2892 else \
2893 for (i = 0; i < nnz1; ++i) \
2894 _SCALE2_(px1[i], n); \
2895 } \
2896 } while (0)
2897
2898 SUM_CASES(MAP);
2899 for (i = 0; i < m; ++i)
2900 if (work[i] >= 0)
2901 *(pi1++) = i + 1;
2902 Matrix_Free(work, m);
2903 UNPROTECT(2); /* x1, i1 */
2904
2905 } else {
2906
2907 SEXP x1 = res;
2908 int *pi1 = NULL;
2909 if (narm_) {
2910 Matrix_Calloc(pi1, m, int);
2911 for (i = 0; i < m; ++i)
2912 pi1[i] = n;
2913 }
2914 SUM_CASES(NOMAP);
2915 if (narm_)
2916 Matrix_Free(pi1, m);
2917
2918 }
2919
2920#undef SUM_LOOP
2921
2922 UNPROTECT(2); /* i0, p0 */
2923 return;
2924}
2925
2926static
2927void Tsparse_colsum(SEXP obj, const char *class,
2928 int m, int n, char di, int narm, int mean,
2929 SEXP res, SEXP iSym, SEXP jSym)
2930{
2931 int narm_ = narm && mean && class[0] != 'n';
2932 if (narm_)
2933 obj = Tsparse_aggregate(obj);
2934 PROTECT(obj);
2935
2936 SEXP i0 = PROTECT(GET_SLOT(obj, iSym)),
2937 j0 = PROTECT(GET_SLOT(obj, jSym));
2938 int *pi0 = INTEGER(i0), *pj0 = INTEGER(j0), j, nnz1 = n;
2939 R_xlen_t k, nnz0 = XLENGTH(i0);
2940
2941 if (IS_S4_OBJECT(res)) {
2942
2943 int *work;
2944 Matrix_Calloc(work, n, int);
2945 if (di != 'N')
2946 for (j = 0; j < n; ++j)
2947 work[j] = j;
2948 else {
2949 if (class[1] != 's') {
2950 for (k = 0; k < nnz0; ++k)
2951 ++work[pj0[k]];
2952 } else {
2953 for (k = 0; k < nnz0; ++k) {
2954 ++work[pj0[k]];
2955 if (pi0[k] != pj0[k])
2956 ++work[pi0[k]];
2957 }
2958 }
2959 nnz1 = 0;
2960 for (j = 0; j < n; ++j)
2961 work[j] = (work[j]) ? nnz1++ : -1;
2962 }
2963
2964 SEXP
2965 j1 = PROTECT(allocVector(INTSXP, nnz1)),
2966 x1 = PROTECT(allocVector(SUM_TYPEOF(class[0]), nnz1));
2967 SET_SLOT(res, Matrix_iSym, j1);
2968 SET_SLOT(res, Matrix_xSym, x1);
2969 int *pj1 = INTEGER(j1);
2970 if (narm_)
2971 for (j = 0; j < nnz1; ++j)
2972 pj1[j] = m;
2973
2974#define SUM_LOOP(_CTYPE0_, _PTR0_, _CTYPE1_, _PTR1_, _MASK_, \
2975 _ZERO_, _ONE_, _NA_, _ISNA_, \
2976 _MAP_, _CAST_, _INCREMENT_, _SCALE2_) \
2977 do { \
2978 _MASK_(_CTYPE0_ *px0 = _PTR0_(x0)); \
2979 _CTYPE1_ *px1 = _PTR1_(x1) ; \
2980 _CTYPE1_ tmp = (di != 'N') ? _ONE_ : _ZERO_; \
2981 for (j = 0; j < nnz1; ++j) \
2982 px1[j] = tmp; \
2983 if (class[1] != 's') { \
2984 for (k = 0; k < nnz0; ++k) { \
2985 if (_ISNA_(px0[k])) { \
2986 if (!narm) \
2987 px1[_MAP_(pj0[k])] = _NA_; \
2988 else if (narm_) \
2989 --pj1[_MAP_(pj0[k])]; \
2990 } else { \
2991 tmp = _CAST_(px0[k]); \
2992 _INCREMENT_(px1[_MAP_(pj0[k])], tmp); \
2993 } \
2994 } \
2995 } else { \
2996 int off; \
2997 for (k = 0; k < nnz0; ++k) { \
2998 off = pi0[k] != pj0[k]; \
2999 if (_ISNA_(px0[k])) { \
3000 if (!narm) { \
3001 px1[_MAP_(pj0[k])] = _NA_; \
3002 if (off) \
3003 px1[_MAP_(pi0[k])] = _NA_; \
3004 } else if (narm_) { \
3005 --pj1[_MAP_(pj0[k])]; \
3006 if (off) \
3007 --pj1[_MAP_(pi0[k])]; \
3008 } \
3009 } else { \
3010 tmp = _CAST_(px0[k]); \
3011 _INCREMENT_(px1[_MAP_(pj0[k])], tmp); \
3012 if (off) \
3013 _INCREMENT_(px1[_MAP_(pi0[k])], tmp); \
3014 } \
3015 } \
3016 } \
3017 if (mean) { \
3018 if (narm_) \
3019 for (j = 0; j < nnz1; ++j) \
3020 _SCALE2_(px1[j], pj1[j]); \
3021 else \
3022 for (j = 0; j < nnz1; ++j) \
3023 _SCALE2_(px1[j], m); \
3024 } \
3025 } while (0)
3026
3027 SUM_CASES(MAP);
3028 for (j = 0; j < n; ++j)
3029 if (work[j] >= 0)
3030 *(pj1++) = j + 1;
3031 Matrix_Free(work, n);
3032 UNPROTECT(2); /* x1, j1 */
3033
3034 } else {
3035
3036 SEXP x1 = res;
3037 int *pj1 = NULL;
3038 if (narm_)
3039 Matrix_Calloc(pj1, n, int);
3040 SUM_CASES(NOMAP);
3041 if (narm_)
3042 Matrix_Free(pj1, n);
3043
3044 }
3045
3046#undef SUM_LOOP
3047
3048 UNPROTECT(3); /* j0, i0, obj */
3049 return;
3050}
3051
3052SEXP sparse_marginsum(SEXP obj, const char *class, int margin,
3053 int narm, int mean, int sparse)
3054{
3055 SEXP dim = GET_SLOT(obj, Matrix_DimSym);
3056 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1],
3057 r = (margin == 0) ? m : n;
3058
3059 SEXP res;
3060 SEXPTYPE type = SUM_TYPEOF(class[0]);
3061 if (sparse) {
3062 char cl[] = ".sparseVector";
3063 cl[0] = (type == CPLXSXP) ? 'z' : ((type == REALSXP) ? 'd' : 'i');
3064 PROTECT(res = newObject(cl));
3065
3066 SEXP length = PROTECT(ScalarInteger(r));
3067 SET_SLOT(res, Matrix_lengthSym, length);
3068 UNPROTECT(1); /* length */
3069 } else {
3070 PROTECT(res = allocVector(type, r));
3071
3072 SEXP dimnames = (class[1] != 's')
3073 ? GET_SLOT(obj, Matrix_DimNamesSym)
3074 : get_symmetrized_DimNames(obj, -1),
3075 marnames = VECTOR_ELT(dimnames, margin);
3076 if (marnames != R_NilValue) {
3077 PROTECT(marnames);
3078 setAttrib(res, R_NamesSymbol, marnames);
3079 UNPROTECT(1); /* marnames */
3080 }
3081 }
3082
3083 char di = 'N';
3084 if (class[1] == 't') {
3085 SEXP diag = GET_SLOT(obj, Matrix_diagSym);
3086 di = *CHAR(STRING_ELT(diag, 0));
3087 }
3088
3089 if (margin == 0) {
3090 if (class[2] == 'C') {
3091 Csparse_rowsum(obj, class, m, n, di, narm, mean, res,
3092 Matrix_iSym);
3093 } else if (class[2] == 'R') {
3094 if (class[1] != 's')
3095 Csparse_colsum(obj, class, n, m, di, narm, mean, res);
3096 else
3097 Csparse_rowsum(obj, class, n, m, di, narm, mean, res,
3098 Matrix_jSym);
3099 } else {
3100 Tsparse_colsum(obj, class, n, m, di, narm, mean, res,
3101 Matrix_jSym, Matrix_iSym);
3102 }
3103 } else {
3104 if (class[2] == 'C') {
3105 if (class[1] != 's')
3106 Csparse_colsum(obj, class, m, n, di, narm, mean, res);
3107 else
3108 Csparse_rowsum(obj, class, m, n, di, narm, mean, res,
3109 Matrix_iSym);
3110 } else if (class[2] == 'R') {
3111 Csparse_rowsum(obj, class, n, m, di, narm, mean, res,
3112 Matrix_jSym);
3113 } else {
3114 Tsparse_colsum(obj, class, m, n, di, narm, mean, res,
3115 Matrix_iSym, Matrix_jSym);
3116 }
3117 }
3118
3119 UNPROTECT(1); /* res */
3120 return res;
3121}
3122
3123/* (row|col)(Sums|Means)(<[CRT]sparseMatrix>) */
3124SEXP R_sparse_marginsum(SEXP obj, SEXP margin,
3125 SEXP narm, SEXP mean, SEXP sparse)
3126{
3127 static const char *valid[] = {
3128 VALID_CSPARSE, VALID_RSPARSE, VALID_TSPARSE, "" };
3129 int ivalid = R_check_class_etc(obj, valid);
3130 if (ivalid < 0)
3131 ERROR_INVALID_CLASS(obj, __func__);
3132
3133 int margin_;
3134 if (TYPEOF(margin) != INTSXP || LENGTH(margin) < 1 ||
3135 ((margin_ = INTEGER(margin)[0]) != 0 && margin_ != 1))
3136 error(_("'%s' must be %d or %d"), "margin", 0, 1);
3137
3138 int narm_;
3139 if (TYPEOF(narm) != LGLSXP || LENGTH(narm) < 1 ||
3140 (narm_ = LOGICAL(narm)[0]) == NA_LOGICAL)
3141 error(_("'%s' must be %s or %s"), "narm", "TRUE", "FALSE");
3142
3143 int mean_;
3144 if (TYPEOF(mean) != LGLSXP || LENGTH(mean) < 1 ||
3145 (mean_ = LOGICAL(mean)[0]) == NA_LOGICAL)
3146 error(_("'%s' must be %s or %s"), "mean", "TRUE", "FALSE");
3147
3148 int sparse_;
3149 if (TYPEOF(sparse) != LGLSXP || LENGTH(sparse) < 1 ||
3150 (sparse_ = LOGICAL(sparse)[0]) == NA_LOGICAL)
3151 error(_("'%s' must be %s or %s"), "sparse", "TRUE", "FALSE");
3152
3153 return sparse_marginsum(obj, valid[ivalid],
3154 margin_, narm_, mean_, sparse_);
3155}
3156
3157#undef SUM_CASES
3158#undef SUM_TYPEOF
3159
3160#define TRY_INCREMENT(_LABEL_) \
3161 do { \
3162 if ((s >= 0) \
3163 ? ( t <= MATRIX_INT_FAST64_MAX - s) \
3164 : (-t <= s - MATRIX_INT_FAST64_MIN)) { \
3165 s += t; \
3166 t = 0; \
3167 count = 0; \
3168 } else { \
3169 over = 1; \
3170 goto _LABEL_; \
3171 } \
3172 } while (0)
3173
3174#define LONGDOUBLE_AS_DOUBLE(v) \
3175 (v > DBL_MAX) ? R_PosInf : ((v < -DBL_MAX) ? R_NegInf : (double) v);
3176
3177SEXP sparse_sum(SEXP obj, const char *class, int narm)
3178{
3179 if (class[2] == 'T')
3180 obj = Tsparse_aggregate(obj);
3181 PROTECT(obj);
3182
3183 SEXP res;
3184
3185 if (!narm && (class[0] == 'l' || class[0] == 'i')) {
3186 SEXP x = GET_SLOT(obj, Matrix_xSym);
3187 int *px = (class[0] == 'l') ? LOGICAL(x) : INTEGER(x);
3188 R_xlen_t nx = XLENGTH(x);
3189 while (nx--) {
3190 if (*px == NA_INTEGER) {
3191 res = allocVector(INTSXP, 1);
3192 INTEGER(res)[0] = NA_INTEGER;
3193 UNPROTECT(1); /* obj */
3194 return res;
3195 }
3196 ++px;
3197 }
3198 }
3199
3200 SEXP dim = GET_SLOT(obj, Matrix_DimSym);
3201 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1];
3202
3203 char di = 'N';
3204 if (class[1] == 't') {
3205 SEXP diag = GET_SLOT(obj, Matrix_diagSym);
3206 di = *CHAR(STRING_ELT(diag, 0));
3207 }
3208
3209 int symmetric = class[1] == 's';
3210
3211 if (class[2] != 'T') {
3212
3213 SEXP iSym = (class[2] == 'C') ? Matrix_iSym : Matrix_jSym,
3214 p = PROTECT(GET_SLOT(obj, Matrix_pSym)),
3215 i = PROTECT(GET_SLOT(obj, iSym));
3216 int *pp = INTEGER(p) + 1, *pi = INTEGER(i), j_, k = 0, kend,
3217 n_ = (class[2] == 'C') ? n : m;
3218
3219 if (class[0] == 'n') {
3220 Matrix_int_fast64_t nnz = pp[n_ - 1];
3221 if (di != 'N')
3222 nnz += n;
3223 if (symmetric) {
3224 SEXP uplo = GET_SLOT(obj, Matrix_uploSym);
3225 char ul = *CHAR(STRING_ELT(uplo, 0));
3226
3227 nnz *= 2;
3228 for (j_ = 0; j_ < n_; ++j_) {
3229 kend = pp[j_];
3230 if (k < kend && pi[(ul == 'U') ? kend - 1 : k] == j_)
3231 --nnz;
3232 k = kend;
3233 }
3234 }
3235 if (nnz <= INT_MAX) {
3236 res = allocVector(INTSXP, 1);
3237 INTEGER(res)[0] = (int) nnz;
3238 } else {
3239 res = allocVector(REALSXP, 1);
3240 REAL(res)[0] = (double) nnz;
3241 }
3242 UNPROTECT(3); /* i, p, obj */
3243 return res;
3244 }
3245
3246 SEXP x = GET_SLOT(obj, Matrix_xSym);
3247 UNPROTECT(2); /* i, p */
3248
3249 if (class[0] == 'z') {
3250 Rcomplex *px = COMPLEX(x);
3251 long double zr = (di == 'N') ? 0.0L : n, zi = 0.0L;
3252 for (j_ = 0; j_ < n_; ++j_) {
3253 kend = pp[j_];
3254 while (k < kend) {
3255 if (!(narm && (ISNAN(px[k].r) || ISNAN(px[k].i)))) {
3256 zr += (symmetric && pi[k] != j_)
3257 ? 2.0L * px[k].r : px[k].r;
3258 zi += (symmetric && pi[k] != j_)
3259 ? 2.0L * px[k].i : px[k].i;
3260 }
3261 ++k;
3262 }
3263 }
3264 res = allocVector(CPLXSXP, 1);
3265 COMPLEX(res)[0].r = LONGDOUBLE_AS_DOUBLE(zr);
3266 COMPLEX(res)[0].i = LONGDOUBLE_AS_DOUBLE(zi);
3267 } else if (class[0] == 'd') {
3268 double *px = REAL(x);
3269 long double zr = (di == 'N') ? 0.0L : n;
3270 for (j_ = 0; j_ < n_; ++j_) {
3271 kend = pp[j_];
3272 while (k < kend) {
3273 if (!(narm && ISNAN(px[k])))
3274 zr += (symmetric && pi[k] != j_)
3275 ? 2.0L * px[k] : px[k];
3276 ++k;
3277 }
3278 }
3279 res = allocVector(REALSXP, 1);
3280 REAL(res)[0] = LONGDOUBLE_AS_DOUBLE(zr);
3281 } else {
3282 int *px = (class[0] == 'l') ? LOGICAL(x) : INTEGER(x);
3283 Matrix_int_fast64_t s = (di == 'N') ? 0LL : n, t = 0LL;
3284 unsigned int count = 0;
3285 int over = 0;
3286 for (j_ = 0; j_ < n_; ++j_) {
3287 kend = pp[j_];
3288 while (k < kend) {
3289 if (!narm || px[k] != NA_INTEGER) {
3290 int d = (symmetric && pi[k] != j_) ? 2 : 1;
3291 if (count > UINT_MAX - d)
3292 TRY_INCREMENT(ifoverC);
3293 t += (d == 2) ? 2LL * px[k] : px[k];
3294 count += d;
3295 }
3296 ++k;
3297 }
3298 }
3299 TRY_INCREMENT(ifoverC);
3300 ifoverC:
3301 if (over) {
3302 long double zr = (long double) s + (long double) t;
3303 for (; j_ < n_; ++j_) {
3304 kend = pp[j_];
3305 while (k < kend) {
3306 if (!narm || px[k] != NA_INTEGER)
3307 zr += (symmetric && pi[k] != j_)
3308 ? 2.0L * px[k] : px[k];
3309 ++k;
3310 }
3311 }
3312 res = allocVector(REALSXP, 1);
3313 REAL(res)[0] = LONGDOUBLE_AS_DOUBLE(zr);
3314 } else if (s > INT_MIN && s <= INT_MAX) {
3315 res = allocVector(INTSXP, 1);
3316 INTEGER(res)[0] = (int) s;
3317 } else {
3318 res = allocVector(REALSXP, 1);
3319 REAL(res)[0] = (double) s;
3320 }
3321 }
3322
3323 } else {
3324
3325 SEXP i = PROTECT(GET_SLOT(obj, Matrix_iSym)),
3326 j = PROTECT(GET_SLOT(obj, Matrix_jSym));
3327 int *pi = INTEGER(i), *pj = INTEGER(j);
3328 R_xlen_t k, kend = XLENGTH(i);
3329
3330 if (class[0] == 'n') {
3331 Matrix_int_fast64_t nnz = (Matrix_int_fast64_t) kend;
3332 if (di != 'N')
3333 nnz += n;
3334 if (symmetric) {
3335 nnz *= 2;
3336 for (k = 0; k < kend; ++k)
3337 if (pi[k] == pj[k])
3338 --nnz;
3339 }
3340 if (nnz <= INT_MAX) {
3341 res = allocVector(INTSXP, 1);
3342 INTEGER(res)[0] = (int) nnz;
3343 } else {
3344 res = allocVector(REALSXP, 1);
3345 REAL(res)[0] = (double) nnz;
3346 }
3347 UNPROTECT(3); /* j, i, obj */
3348 return res;
3349 }
3350
3351 SEXP x = GET_SLOT(obj, Matrix_xSym);
3352 UNPROTECT(2); /* j, i */
3353
3354 if (class[0] == 'z') {
3355 Rcomplex *px = COMPLEX(x);
3356 long double zr = (di == 'N') ? 0.0L : n, zi = 0.0L;
3357 for (k = 0; k < kend; ++k)
3358 if (!(narm && (ISNAN(px[k].r) || ISNAN(px[k].i)))) {
3359 zr += (symmetric && pi[k] != pj[k])
3360 ? 2.0L * px[k].r : px[k].r;
3361 zi += (symmetric && pi[k] != pj[k])
3362 ? 2.0L * px[k].i : px[k].i;
3363 }
3364 res = allocVector(CPLXSXP, 1);
3365 COMPLEX(res)[0].r = LONGDOUBLE_AS_DOUBLE(zr);
3366 COMPLEX(res)[0].i = LONGDOUBLE_AS_DOUBLE(zi);
3367 } else if (class[0] == 'd') {
3368 double *px = REAL(x);
3369 long double zr = (di == 'N') ? 0.0L : n;
3370 for (k = 0; k < kend; ++k)
3371 if (!(narm && ISNAN(px[k])))
3372 zr += (symmetric && pi[k] != pj[k])
3373 ? 2.0L * px[k] : px[k];
3374 res = allocVector(REALSXP, 1);
3375 REAL(res)[0] = LONGDOUBLE_AS_DOUBLE(zr);
3376 } else {
3377 int *px = (class[0] == 'i') ? INTEGER(x) : LOGICAL(x);
3378 Matrix_int_fast64_t s = (di == 'N') ? 0LL : n, t = 0LL;
3379 unsigned int count = 0;
3380 int over = 0;
3381 for (k = 0; k < kend; ++k) {
3382 if (!narm || px[k] != NA_INTEGER) {
3383 int d = (symmetric && pi[k] != pj[k]) ? 2 : 1;
3384 if (count > UINT_MAX - d)
3385 TRY_INCREMENT(ifoverT);
3386 t += (d == 2) ? 2LL * px[k] : px[k];
3387 count += d;
3388 }
3389 }
3390 TRY_INCREMENT(ifoverT);
3391 ifoverT:
3392 if (over) {
3393 long double zr = (long double) s + (long double) t;
3394 for (; k < kend; ++k)
3395 if (!(narm && px[k] == NA_INTEGER))
3396 zr += (symmetric && pi[k] != pj[k])
3397 ? 2.0L * px[k] : px[k];
3398 res = allocVector(REALSXP, 1);
3399 REAL(res)[0] = LONGDOUBLE_AS_DOUBLE(zr);
3400 } else if (s > INT_MIN && s <= INT_MAX) {
3401 res = allocVector(INTSXP, 1);
3402 INTEGER(res)[0] = (int) s;
3403 } else {
3404 res = allocVector(REALSXP, 1);
3405 REAL(res)[0] = (double) s;
3406 }
3407 }
3408
3409 }
3410
3411 UNPROTECT(1); /* obj */
3412 return res;
3413}
3414
3415/* sum(<[CRT]sparseMatrix>) */
3416SEXP R_sparse_sum(SEXP obj, SEXP narm)
3417{
3418 static const char *valid[] = {
3419 VALID_CSPARSE, VALID_RSPARSE, VALID_TSPARSE, "" };
3420 int ivalid = R_check_class_etc(obj, valid);
3421 if (ivalid < 0)
3422 ERROR_INVALID_CLASS(obj, __func__);
3423
3424 int narm_;
3425 if (TYPEOF(narm) != LGLSXP || LENGTH(narm) < 1 ||
3426 (narm_ = LOGICAL(narm)[0]) == NA_LOGICAL)
3427 error(_("'%s' must be %s or %s"), "narm", "TRUE", "FALSE");
3428
3429 return sparse_sum(obj, valid[ivalid], narm_);
3430}
3431
3432SEXP sparse_prod(SEXP obj, const char *class, int narm)
3433{
3434 if (class[2] == 'T')
3435 obj = Tsparse_aggregate(obj);
3436 PROTECT(obj);
3437
3438 SEXP res = PROTECT(allocVector((class[0] == 'z') ? CPLXSXP : REALSXP, 1));
3439
3440 SEXP dim = GET_SLOT(obj, Matrix_DimSym);
3441 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1];
3442
3443 char ul = 'U', di = 'N';
3444 if (class[1] != 'g') {
3445 SEXP uplo = GET_SLOT(obj, Matrix_uploSym);
3446 ul = *CHAR(STRING_ELT(uplo, 0));
3447 if (class[1] == 't') {
3448 SEXP diag = GET_SLOT(obj, Matrix_diagSym);
3449 di = *CHAR(STRING_ELT(diag, 0));
3450 }
3451 }
3452
3453 int symmetric = (class[1] != 's')
3454 ? 0 : (((class[2] == 'C') == (ul == 'U')) ? 1 : -1);
3455 long double zr = 1.0L, zi = 0.0L;
3456
3457 Matrix_int_fast64_t mn = (Matrix_int_fast64_t) m * n,
3458 nnz, nnzmax = (symmetric) ? (mn + n) / 2 : mn;
3459
3460 if (class[2] != 'T') {
3461
3462 SEXP iSym = (class[2] == 'C') ? Matrix_iSym : Matrix_jSym,
3463 p = PROTECT(GET_SLOT(obj, Matrix_pSym)),
3464 i = PROTECT(GET_SLOT(obj, iSym));
3465 int *pp = INTEGER(p) + 1, *pi = INTEGER(i), i_, j_, k = 0, kend,
3466 m_ = (class[2] == 'C') ? m : n, n_ = (class[2] == 'C') ? n : m,
3467 seen0 = 0;
3468
3469 nnz = pp[n_ - 1];
3470 if (di != 'N')
3471 nnz += n;
3472 if (class[0] == 'n') {
3473 REAL(res)[0] = (nnz == nnzmax) ? 1.0 : 0.0;
3474 UNPROTECT(4); /* i, p, res, obj */
3475 return res;
3476 }
3477
3478 SEXP x = GET_SLOT(obj, Matrix_xSym);
3479 UNPROTECT(2); /* i, p */
3480
3481 if (class[0] == 'z') {
3482 Rcomplex *px = COMPLEX(x);
3483 long double zr0, zi0;
3484 for (j_ = 0; j_ < n_; ++j_) {
3485 kend = pp[j_];
3486 if (seen0 || kend - k == m_) {
3487 while (k < kend) {
3488 if (!(narm && (ISNAN(px[k].r) || ISNAN(px[k].i)))) {
3489 zr0 = zr; zi0 = zi;
3490 zr = zr0 * px[k].r - zi0 * px[k].i;
3491 zi = zr0 * px[k].i + zi0 * px[k].r;
3492 if (symmetric && pi[k] != j_) {
3493 zr0 = zr; zi0 = zi;
3494 zr = zr0 * px[k].r - zi0 * px[k].i;
3495 zi = zr0 * px[k].i + zi0 * px[k].r;
3496 }
3497 }
3498 ++k;
3499 }
3500 } else {
3501 int i0 = (symmetric >= 0) ? 0 : j_,
3502 i1 = (symmetric <= 0) ? m_ : j_ + 1;
3503 for (i_ = i0; i_ < i1; ++i_) {
3504 if (seen0 || (k < kend && pi[k] == i_)) {
3505 if (k >= kend)
3506 break;
3507 if (!(narm && (ISNAN(px[k].r) || ISNAN(px[k].i)))) {
3508 zr0 = zr; zi0 = zi;
3509 zr = zr0 * px[k].r - zi0 * px[k].i;
3510 zi = zr0 * px[k].i + zi0 * px[k].r;
3511 if (symmetric && pi[k] != j_) {
3512 zr0 = zr; zi0 = zi;
3513 zr = zr0 * px[k].r - zi0 * px[k].i;
3514 zi = zr0 * px[k].i + zi0 * px[k].r;
3515 }
3516 }
3517 ++k;
3518 } else if (di == 'N' || i_ != j_) {
3519 zr *= 0.0L;
3520 zi *= 0.0L;
3521 seen0 = 1;
3522 }
3523 }
3524 }
3525 }
3526 } else if (class[0] == 'd') {
3527 double *px = REAL(x);
3528 for (j_ = 0; j_ < n_; ++j_) {
3529 kend = pp[j_];
3530 if (seen0 || kend - k == m_) {
3531 while (k < kend) {
3532 if (!(narm && ISNAN(px[k])))
3533 zr *= (symmetric && pi[k] != j_)
3534 ? (long double) px[k] * px[k] : px[k];
3535 ++k;
3536 }
3537 } else {
3538 int i0 = (symmetric >= 0) ? 0 : j_,
3539 i1 = (symmetric <= 0) ? m_ : j_ + 1;
3540 for (i_ = i0; i_ < i1; ++i_) {
3541 if (seen0 || (k < kend && pi[k] == i_)) {
3542 if (k >= kend)
3543 break;
3544 if (!(narm && ISNAN(px[k])))
3545 zr *= (symmetric && pi[k] != j_)
3546 ? (long double) px[k] * px[k] : px[k];
3547 ++k;
3548 } else if (di == 'N' || i_ != j_) {
3549 zr *= 0.0L;
3550 seen0 = 1;
3551 }
3552 }
3553 }
3554 }
3555 } else {
3556 int *px = (class[0] == 'l') ? LOGICAL(x) : INTEGER(x);
3557 for (j_ = 0; j_ < n_; ++j_) {
3558 kend = pp[j_];
3559 if (seen0 || kend - k == m_) {
3560 while (k < kend) {
3561 if (px[k] != NA_INTEGER)
3562 zr *= (symmetric && pi[k] != j_)
3563 ? (long double) px[k] * px[k] : px[k];
3564 else if (!narm)
3565 zr *= NA_REAL;
3566 ++k;
3567 }
3568 } else {
3569 int i0 = (symmetric >= 0) ? 0 : j_,
3570 i1 = (symmetric <= 0) ? m_ : j_ + 1;
3571 for (i_ = i0; i_ < i1; ++i_) {
3572 if (seen0 || (k < kend && pi[k] == i_)) {
3573 if (k >= kend)
3574 break;
3575 if (px[k] != NA_INTEGER)
3576 zr *= (symmetric && pi[k] != j_)
3577 ? (long double) px[k] * px[k] : px[k];
3578 else if (!narm)
3579 zr *= NA_REAL;
3580 ++k;
3581 } else if (di == 'N' || i_ != j_) {
3582 zr *= 0.0L;
3583 seen0 = 1;
3584 }
3585 }
3586 }
3587 }
3588 }
3589
3590 } else {
3591
3592 SEXP i = PROTECT(GET_SLOT(obj, Matrix_iSym)),
3593 j = PROTECT(GET_SLOT(obj, Matrix_jSym));
3594 int *pi = INTEGER(i), *pj = INTEGER(j);
3595 R_xlen_t k, kend = XLENGTH(i);
3596
3597 nnz = (Matrix_int_fast64_t) kend;
3598 if (di != 'N')
3599 nnz += n;
3600 if (class[0] == 'n') {
3601 REAL(res)[0] = (nnz == nnzmax) ? 1.0 : 0.0;
3602 UNPROTECT(4); /* j, i, res, obj */
3603 return res;
3604 }
3605 if (nnz < nnzmax)
3606 zr = 0.0;
3607
3608 SEXP x = GET_SLOT(obj, Matrix_xSym);
3609 UNPROTECT(2); /* j, i */
3610
3611 if (class[0] == 'z') {
3612 Rcomplex *px = COMPLEX(x);
3613 long double zr0, zi0;
3614 for (k = 0; k < kend; ++k)
3615 if (!(narm && (ISNAN(px[k].r) || ISNAN(px[k].i)))) {
3616 zr0 = zr; zi0 = zi;
3617 zr = zr0 * px[k].r - zi0 * px[k].i;
3618 zi = zr0 * px[k].i + zi0 * px[k].r;
3619 if (symmetric && pi[k] != pj[k]) {
3620 zr0 = zr; zi0 = zi;
3621 zr = zr0 * px[k].r - zi0 * px[k].i;
3622 zi = zr0 * px[k].i + zi0 * px[k].r;
3623 }
3624 }
3625 } else if (class[0] == 'd') {
3626 double *px = REAL(x);
3627 for (k = 0; k < kend; ++k)
3628 if (!(narm && ISNAN(px[k])))
3629 zr *= (symmetric && pi[k] != pj[k])
3630 ? (long double) px[k] * px[k] : px[k];
3631 } else {
3632 int *px = (class[0] == 'l') ? LOGICAL(x) : INTEGER(x);
3633 for (k = 0; k < kend; ++k)
3634 if (px[k] != NA_INTEGER)
3635 zr *= (symmetric && pi[k] != pj[k])
3636 ? (long double) px[k] * px[k] : px[k];
3637 else if (!narm)
3638 zr *= NA_REAL;
3639 }
3640
3641 }
3642
3643 if (class[0] == 'z') {
3644 COMPLEX(res)[0].r = LONGDOUBLE_AS_DOUBLE(zr);
3645 COMPLEX(res)[0].i = LONGDOUBLE_AS_DOUBLE(zi);
3646 } else
3647 REAL(res)[0] = LONGDOUBLE_AS_DOUBLE(zr);
3648 UNPROTECT(2); /* res, obj */
3649 return res;
3650}
3651
3652/* prod(<[CRT]sparseMatrix>) */
3653SEXP R_sparse_prod(SEXP obj, SEXP narm)
3654{
3655 static const char *valid[] = {
3656 VALID_CSPARSE, VALID_RSPARSE, VALID_TSPARSE, "" };
3657 int ivalid = R_check_class_etc(obj, valid);
3658 if (ivalid < 0)
3659 ERROR_INVALID_CLASS(obj, __func__);
3660
3661 int narm_;
3662 if (TYPEOF(narm) != LGLSXP || LENGTH(narm) < 1 ||
3663 (narm_ = LOGICAL(narm)[0]) == NA_LOGICAL)
3664 error(_("'%s' must be %s or %s"), "narm", "TRUE", "FALSE");
3665
3666 return sparse_prod(obj, valid[ivalid], narm_);
3667}
3668
3669#undef TRY_INCREMENT
3670#undef LONGDOUBLE_AS_DOUBLE
3671
3672SEXP Tsparse_aggregate(SEXP from)
3673{
3674 static const char *valid[] = { VALID_TSPARSE, "" };
3675 int ivalid = R_check_class_etc(from, valid);
3676 if (ivalid < 0)
3677 ERROR_INVALID_CLASS(from, __func__);
3678 const char *cl = valid[ivalid];
3679
3680 SEXP dim = PROTECT(GET_SLOT(from, Matrix_DimSym));
3681 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1];
3682 UNPROTECT(1); /* dim */
3683
3684 SEXP to,
3685 i0 = PROTECT(GET_SLOT(from, Matrix_iSym)),
3686 j0 = PROTECT(GET_SLOT(from, Matrix_jSym)),
3687 i1 = NULL, j1 = NULL;
3688
3689 /* defined in ./coerce.c : */
3690 void taggr(SEXP, SEXP, SEXP, SEXP *, SEXP *, SEXP *, int, int);
3691
3692 if (cl[0] == 'n') {
3693 taggr(j0, i0, NULL, &j1, &i1, NULL, n, m);
3694 if (!i1) {
3695 UNPROTECT(2); /* j0, i0 */
3696 return from;
3697 }
3698 PROTECT(i1);
3699 PROTECT(j1);
3700 PROTECT(to = newObject(cl));
3701 SET_SLOT(to, Matrix_iSym, i1);
3702 SET_SLOT(to, Matrix_jSym, j1);
3703 UNPROTECT(5); /* to, j1, i1, j0, i0 */
3704 } else {
3705 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym)),
3706 x1 = NULL;
3707 taggr(j0, i0, x0, &j1, &i1, &x1, n, m);
3708 if (!i1) {
3709 UNPROTECT(3); /* x0, j0, i0 */
3710 return from;
3711 }
3712 PROTECT(i1);
3713 PROTECT(j1);
3714 PROTECT(x1);
3715 PROTECT(to = newObject(cl));
3716 SET_SLOT(to, Matrix_iSym, i1);
3717 SET_SLOT(to, Matrix_jSym, j1);
3718 SET_SLOT(to, Matrix_xSym, x1);
3719 UNPROTECT(7); /* to, x1, j1, i1, x0, j0, i0 */
3720 }
3721
3722 PROTECT(to);
3723
3724 if (m != n || n > 0) {
3725 PROTECT(dim = GET_SLOT(to, Matrix_DimSym));
3726 pdim = INTEGER(dim);
3727 pdim[0] = m;
3728 pdim[1] = n;
3729 UNPROTECT(1); /* dim */
3730 }
3731
3732 SEXP dimnames = PROTECT(GET_SLOT(from, Matrix_DimNamesSym));
3733 SET_SLOT(to, Matrix_DimNamesSym, dimnames);
3734 UNPROTECT(1); /* dimnames */
3735
3736 if (cl[1] != 'g') {
3737 SEXP uplo = PROTECT(GET_SLOT(from, Matrix_uploSym));
3738 char ul = *CHAR(STRING_ELT(uplo, 0));
3739 if (ul != 'U')
3740 SET_SLOT(to, Matrix_uploSym, uplo);
3741 UNPROTECT(1); /* uplo */
3742 }
3743 if (cl[1] == 't') {
3744 SEXP diag = PROTECT(GET_SLOT(from, Matrix_diagSym));
3745 char di = *CHAR(STRING_ELT(diag, 0));
3746 if (di != 'N')
3747 SET_SLOT(to, Matrix_diagSym, diag);
3748 UNPROTECT(1); /* diag */
3749 } else {
3750 SEXP factors = PROTECT(GET_SLOT(from, Matrix_factorsSym));
3751 if (LENGTH(factors) > 0)
3752 SET_SLOT(to, Matrix_factorsSym, factors);
3753 UNPROTECT(1); /* factors */
3754 }
3755
3756 UNPROTECT(1); /* to */
3757 return to;
3758}
VALID_RSPARSE
#define VALID_RSPARSE
Definition Mdefines.h:264
Matrix_int_fast64_t
long long Matrix_int_fast64_t
Definition Mdefines.h:27
NOTCONJ_LOGICAL
#define NOTCONJ_LOGICAL(_X_, _Y_)
Definition Mdefines.h:133
VALID_TSPARSE
#define VALID_TSPARSE
Definition Mdefines.h:271
NOTREAL_LOGICAL
#define NOTREAL_LOGICAL(_X_)
Definition Mdefines.h:126
ERROR_INVALID_CLASS
#define ERROR_INVALID_CLASS(_X_, _FUNC_)
Definition Mdefines.h:208
NOTREAL_INTEGER
#define NOTREAL_INTEGER(_X_)
Definition Mdefines.h:127
NOTREAL_PATTERN
#define NOTREAL_PATTERN(_X_)
Definition Mdefines.h:125
NOTREAL_COMPLEX
#define NOTREAL_COMPLEX(_X_)
Definition Mdefines.h:129
_
#define _(String)
Definition Mdefines.h:44
NOTREAL_REAL
#define NOTREAL_REAL(_X_)
Definition Mdefines.h:128
INCREMENT_REAL
#define INCREMENT_REAL(_X_, _Y_)
Definition Mdefines.h:169
typeToKind
char typeToKind(SEXPTYPE)
Definition objects.c:11
ASSIGN_COMPLEX
#define ASSIGN_COMPLEX(_X_, _Y_)
Definition Mdefines.h:181
INCREMENT_COMPLEX
#define INCREMENT_COMPLEX(_X_, _Y_)
Definition Mdefines.h:173
ASSIGN_REAL
#define ASSIGN_REAL(_X_, _Y_)
Definition Mdefines.h:179
kindToType
SEXPTYPE kindToType(char)
Definition objects.c:28
NOTCONJ_REAL
#define NOTCONJ_REAL(_X_, _Y_)
Definition Mdefines.h:137
NOTCONJ_COMPLEX
#define NOTCONJ_COMPLEX(_X_, _Y_)
Definition Mdefines.h:139
HIDE
#define HIDE(...)
Definition Mdefines.h:202
ISNZ_LOGICAL
#define ISNZ_LOGICAL(_X_)
Definition Mdefines.h:109
SHOW
#define SHOW(...)
Definition Mdefines.h:201
VALID_CSPARSE
#define VALID_CSPARSE
Definition Mdefines.h:257
Matrix_Free
#define Matrix_Free(_VAR_, _N_)
Definition Mdefines.h:77
INCREMENT_PATTERN
#define INCREMENT_PATTERN(_X_, _Y_)
Definition Mdefines.h:143
SET_SLOT
#define SET_SLOT(x, what, value)
Definition Mdefines.h:86
GET_SLOT
#define GET_SLOT(x, what)
Definition Mdefines.h:85
Matrix_Calloc
#define Matrix_Calloc(_VAR_, _N_, _CTYPE_)
Definition Mdefines.h:66
newObject
SEXP newObject(const char *)
Definition objects.c:4
SECONDOF
#define SECONDOF(x, y)
Definition Mdefines.h:100
NOTCONJ_INTEGER
#define NOTCONJ_INTEGER(_X_, _Y_)
Definition Mdefines.h:135
Matrix_DimSym
SEXP Matrix_DimSym
Definition Msymbols.h:3
Matrix_factorsSym
SEXP Matrix_factorsSym
Definition Msymbols.h:12
Matrix_xSym
SEXP Matrix_xSym
Definition Msymbols.h:22
Matrix_lengthSym
SEXP Matrix_lengthSym
Definition Msymbols.h:15
Matrix_iSym
SEXP Matrix_iSym
Definition Msymbols.h:13
Matrix_jSym
SEXP Matrix_jSym
Definition Msymbols.h:14
Matrix_DimNamesSym
SEXP Matrix_DimNamesSym
Definition Msymbols.h:2
Matrix_diagSym
SEXP Matrix_diagSym
Definition Msymbols.h:11
Matrix_uploSym
SEXP Matrix_uploSym
Definition Msymbols.h:21
Matrix_pSym
SEXP Matrix_pSym
Definition Msymbols.h:17
get_symmetrized_DimNames
SEXP get_symmetrized_DimNames(SEXP obj, int J)
Definition attrib.c:108
DimNames_is_symmetric
int DimNames_is_symmetric(SEXP dn)
Definition attrib.c:14
set_reversed_DimNames
void set_reversed_DimNames(SEXP obj, SEXP dn)
Definition attrib.c:142
set_symmetrized_DimNames
void set_symmetrized_DimNames(SEXP obj, SEXP dn, int J)
Definition attrib.c:132
valid
static const char * valid[]
Definition bind.c:5
cl
cholmod_common cl
Definition cholmod-etc.c:6
sparse_as_Csparse
SEXP sparse_as_Csparse(SEXP from, const char *class)
Definition coerce.c:3721
taggr
void taggr(SEXP i0, SEXP j0, SEXP x0, SEXP *i1, SEXP *j1, SEXP *x1, int m, int n)
Definition coerce.c:3585
sparse_as_kind
SEXP sparse_as_kind(SEXP from, const char *class, char kind)
Definition coerce.c:2499
trans
void trans(SEXP p0, SEXP i0, SEXP x0, SEXP p1, SEXP i1, SEXP x1, int m, int n)
Definition coerce.c:3356
Matrix_zone
Rcomplex Matrix_zone
Definition init.c:26
FS_SUBCASES
#define FS_SUBCASES(_CTYPE_, _PTR_, _MASK_, _ONE_)
DROP0_LOOP2
#define DROP0_LOOP2(_CTYPE_, _PTR_, _ISNZ_)
FS_CASES
#define FS_CASES
R_sparse_band
SEXP R_sparse_band(SEXP from, SEXP k1, SEXP k2)
Definition sparse.c:574
sparse_sum
SEXP sparse_sum(SEXP obj, const char *class, int narm)
Definition sparse.c:3177
NOTCONJ_PATTERN
#define NOTCONJ_PATTERN(_X_, _Y_)
R_sparse_is_diagonal
SEXP R_sparse_is_diagonal(SEXP obj)
Definition sparse.c:2642
Csparse_colsum
static void Csparse_colsum(SEXP obj, const char *class, int m, int n, char di, int narm, int mean, SEXP res)
Definition sparse.c:2711
sparse_diag_N2U
SEXP sparse_diag_N2U(SEXP from, const char *class)
Definition sparse.c:233
DS_CASES
#define DS_CASES
R_sparse_skewpart
SEXP R_sparse_skewpart(SEXP from)
Definition sparse.c:2312
sparse_force_symmetric
SEXP sparse_force_symmetric(SEXP from, const char *class, char ul)
Definition sparse.c:1239
SUM_TYPEOF
#define SUM_TYPEOF(c)
Definition sparse.c:2708
Tsparse_colsum
static void Tsparse_colsum(SEXP obj, const char *class, int m, int n, char di, int narm, int mean, SEXP res, SEXP iSym, SEXP jSym)
Definition sparse.c:2927
DROP0_LOOP1
#define DROP0_LOOP1(_CTYPE_, _PTR_, _ISNZ_)
sparse_diag_U2N
SEXP sparse_diag_U2N(SEXP from, const char *class)
Definition sparse.c:201
R_sparse_diag_U2N
SEXP R_sparse_diag_U2N(SEXP from)
Definition sparse.c:222
R_sparse_diag_N2U
SEXP R_sparse_diag_N2U(SEXP from)
Definition sparse.c:270
Csparse_rowsum
static void Csparse_rowsum(SEXP obj, const char *class, int m, int n, char di, int narm, int mean, SEXP res, SEXP iSym)
Definition sparse.c:2793
DS_LOOP
#define DS_LOOP(_CTYPE_, _PTR_, _MASK_, _ISNZ_)
DROP0_CASES
#define DROP0_CASES(_DO_)
R_sparse_transpose
SEXP R_sparse_transpose(SEXP from, SEXP lazy)
Definition sparse.c:1223
SP_LOOP
#define SP_LOOP(_CTYPE_, _PTR_, _ASSIGN_, _INCREMENT_)
R_sparse_drop0
SEXP R_sparse_drop0(SEXP from, SEXP tol)
Definition sparse.c:185
R_sparse_sum
SEXP R_sparse_sum(SEXP obj, SEXP narm)
Definition sparse.c:3416
sparse_is_triangular
int sparse_is_triangular(SEXP obj, const char *class, int upper)
Definition sparse.c:2454
MAP
#define MAP(_I_)
Definition sparse.c:2653
R_sparse_force_symmetric
SEXP R_sparse_force_symmetric(SEXP from, SEXP uplo)
Definition sparse.c:1613
LONGDOUBLE_AS_DOUBLE
#define LONGDOUBLE_AS_DOUBLE(v)
Definition sparse.c:3174
sparse_transpose
SEXP sparse_transpose(SEXP from, const char *class, int lazy)
Definition sparse.c:1106
R_sparse_prod
SEXP R_sparse_prod(SEXP obj, SEXP narm)
Definition sparse.c:3653
NOMAP
#define NOMAP(_I_)
Definition sparse.c:2654
sparse_symmpart
SEXP sparse_symmpart(SEXP from, const char *class)
Definition sparse.c:1632
R_sparse_is_triangular
SEXP R_sparse_is_triangular(SEXP obj, SEXP upper)
Definition sparse.c:2571
R_sparse_symmpart
SEXP R_sparse_symmpart(SEXP from)
Definition sparse.c:2052
R_sparse_diag_set
SEXP R_sparse_diag_set(SEXP from, SEXP value)
Definition sparse.c:1052
R_sparse_diag_get
SEXP R_sparse_diag_get(SEXP obj, SEXP names)
Definition sparse.c:773
sparse_band
SEXP sparse_band(SEXP from, const char *class, int a, int b)
Definition sparse.c:281
BAND_CASES
#define BAND_CASES
DG_LOOP
#define DG_LOOP(_CTYPE_, _PTR_, _MASK_, _REPLACE_, _INCREMENT_, _ZERO_, _ONE_)
sparse_diag_set
SEXP sparse_diag_set(SEXP from, const char *class, SEXP value)
Definition sparse.c:789
BAND_SUBCASES
#define BAND_SUBCASES(_CTYPE_, _PTR_, _MASK_)
TRY_INCREMENT
#define TRY_INCREMENT(_LABEL_)
Definition sparse.c:3160
sparse_marginsum
SEXP sparse_marginsum(SEXP obj, const char *class, int margin, int narm, int mean, int sparse)
Definition sparse.c:3052
DG_CASES
#define DG_CASES
IS_LOOP
#define IS_LOOP(_CTYPE_, _PTR_, _MASK_, _NOTREAL_, _NOTCONJ_)
Tsparse_aggregate
SEXP Tsparse_aggregate(SEXP from)
Definition sparse.c:3672
sparse_is_symmetric
int sparse_is_symmetric(SEXP obj, const char *class, int checkDN)
Definition sparse.c:2323
sparse_is_diagonal
int sparse_is_diagonal(SEXP obj, const char *class)
Definition sparse.c:2598
sparse_skewpart
SEXP sparse_skewpart(SEXP from, const char *class)
Definition sparse.c:2063
sparse_drop0
SEXP sparse_drop0(SEXP from, const char *class, double tol)
Definition sparse.c:5
sparse_diag_get
SEXP sparse_diag_get(SEXP obj, const char *class, int names)
Definition sparse.c:604
R_sparse_is_symmetric
SEXP R_sparse_is_symmetric(SEXP obj, SEXP checkDN)
Definition sparse.c:2438
R_sparse_marginsum
SEXP R_sparse_marginsum(SEXP obj, SEXP margin, SEXP narm, SEXP mean, SEXP sparse)
Definition sparse.c:3124
SUM_CASES
#define SUM_CASES(_MAP_)
Definition sparse.c:2662
sparse_prod
SEXP sparse_prod(SEXP obj, const char *class, int narm)
Definition sparse.c:3432
Matrix_memset
void * Matrix_memset(void *dest, int ch, R_xlen_t length, size_t size)
Definition utils.c:8
zeroIm
void zeroIm(SEXP x)
Definition utils.c:174
conjugate
void conjugate(SEXP x)
Definition utils.c:152
equal_character_vectors
int equal_character_vectors(SEXP s1, SEXP s2, int n)
Definition utils.c:143
Matrix_memcpy
void * Matrix_memcpy(void *dest, const void *src, R_xlen_t length, size_t size)
Definition utils.c:70
zeroRe
void zeroRe(SEXP x)
Definition utils.c:163
Generated by doxygen 1.9.8