Matrix r4655
Loading...
Searching...
No Matches
dense.c
Go to the documentation of this file.
1#include "Mdefines.h"
2#include "idz.h"
3#include "dense.h"
4
5SEXP dense_band(SEXP from, const char *class, int a, int b)
6{
7 /* defined in ./coerce.c : */
8 SEXP dense_as_general(SEXP, const char *, int);
9
10 SEXP dim = PROTECT(GET_SLOT(from, Matrix_DimSym));
11 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1];
12 UNPROTECT(1); /* dim */
13
14 /* Need tri[ul](<0-by-0>) and tri[ul](<1-by-1>) to be triangularMatrix */
15 if (a <= 1 - m && b >= n - 1 &&
16 (class[1] == 't' || m != n || m > 1 || n > 1))
17 return from;
18
19 int ge = 0, sy = 0, tr = 0;
20 ge = m != n || !((tr = a >= 0 || b <= 0 || class[1] == 't') ||
21 (sy = a == -b && class[1] == 's'));
22
23#define BAND_CASES(_DO_) \
24 do { \
25 switch (class[0]) { \
26 case 'n': \
27 case 'l': \
28 _DO_(i, int, LOGICAL); \
29 break; \
30 case 'i': \
31 _DO_(i, int, INTEGER); \
32 break; \
33 case 'd': \
34 _DO_(d, double, REAL); \
35 break; \
36 case 'z': \
37 _DO_(z, Rcomplex, COMPLEX); \
38 break; \
39 default: \
40 break; \
41 } \
42 } while (0)
43
44#define BAND2(_PREFIX_, _CTYPE_, _PTR_) \
45 _PREFIX_ ## band2(_PTR_(x1), m, n, a, b, di)
46
47#define BAND1(_PREFIX_, _CTYPE_, _PTR_) \
48 _PREFIX_ ## band1(_PTR_(x1), n, a, b, ul1, di)
49
50#define DCPY2(_PREFIX_, _CTYPE_, _PTR_) \
51 do { \
52 _CTYPE_ *px0 = _PTR_(x0), *px1 = _PTR_(x1); \
53 Matrix_memset(px1, 0, XLENGTH(x1), sizeof(_CTYPE_)); \
54 if (a <= 0 && b >= 0) \
55 _PREFIX_ ## dcpy2(px1, px0, n, XLENGTH(x1), 'U', di); \
56 } while (0)
57
58#define DCPY1(_PREFIX_, _CTYPE_, _PTR_) \
59 do { \
60 _CTYPE_ *px0 = _PTR_(x0), *px1 = _PTR_(x1); \
61 Matrix_memset(px1, 0, XLENGTH(x1), sizeof(_CTYPE_)); \
62 if (a <= 0 && b >= 0) \
63 _PREFIX_ ## dcpy1(px1, px0, n, XLENGTH(x1), ul1, ul0, di); \
64 } while (0)
65
66 char ul0 = 'U', ul1 = 'U', di = 'N';
67 if (class[1] != 'g') {
68 if (ge) {
69 PROTECT(from = dense_as_general(from, class, 1));
70 SEXP x1 = PROTECT(GET_SLOT(from, Matrix_xSym));
71 BAND_CASES(BAND2);
72 UNPROTECT(2); /* x1, from */
73 return from;
74 }
75
76 SEXP uplo = PROTECT(GET_SLOT(from, Matrix_uploSym));
77 ul0 = *CHAR(STRING_ELT(uplo, 0));
78 UNPROTECT(1); /* uplo */
79
80 if (class[1] == 't') {
81 /* Be fast if band contains entire triangle */
82 if ((ul0 == 'U')
83 ? (a <= 0 && b >= n - 1) : (b >= 0 && a <= 1 - m))
84 return from;
85 else if (a <= 0 && b >= 0) {
86 SEXP diag = PROTECT(GET_SLOT(from, Matrix_diagSym));
87 di = *CHAR(STRING_ELT(diag, 0));
88 UNPROTECT(1); /* diag */
89 }
90 }
91 }
92
93 char cl[] = "...Matrix";
94 cl[0] = class[0];
95 cl[1] = (ge) ? 'g' : ((sy) ? 's' : 't') ;
96 cl[2] = (ge) ? 'e' : ((class[2] != 'p') ? ((sy) ? 'y' : 'r') : 'p');
97 SEXP to = PROTECT(newObject(cl));
98
99 dim = GET_SLOT(to, Matrix_DimSym);
100 pdim = INTEGER(dim);
101 pdim[0] = m;
102 pdim[1] = n;
103
104 SEXP dimnames = PROTECT(GET_SLOT(from, Matrix_DimNamesSym));
105 if (class[1] != 's' || sy)
106 SET_SLOT(to, Matrix_DimNamesSym, dimnames);
107 else
108 set_symmetrized_DimNames(to, dimnames, -1);
109 UNPROTECT(1); /* dimnames */
110
111 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym)), x1;
112
113 if (ge) {
114 PROTECT(x1 = duplicate(x0));
115 if (ATTRIB(x1) != R_NilValue) {
116 SET_ATTRIB(x1, R_NilValue);
117 if (OBJECT(x1))
118 SET_OBJECT(x1, 0);
119 }
120 SET_SLOT(to, Matrix_xSym, x1);
121 BAND_CASES(BAND2);
122 UNPROTECT(3); /* x1, x0, to */
123 return to;
124 }
125
126 /* Returning .(sy|sp|tr|tp)Matrix ... */
127
128 ul1 = (tr && class[1] != 't') ? ((a >= 0) ? 'U' : 'L') : ul0;
129 if (ul1 != 'U') {
130 SEXP uplo = PROTECT(mkString("L"));
131 SET_SLOT(to, Matrix_uploSym, uplo);
132 UNPROTECT(1); /* uplo */
133 }
134 if (di != 'N') {
135 SEXP diag = PROTECT(mkString("U"));
136 SET_SLOT(to, Matrix_diagSym, diag);
137 UNPROTECT(1); /* diag */
138 }
139
140 if (tr && class[1] == 't') {
141 if ((ul0 == 'U') ? (b <= 0) : (a >= 0)) {
142 /* Result is either a diagonal matrix or a zero matrix : */
143 PROTECT(x1 = allocVector(TYPEOF(x0), XLENGTH(x0)));
144 if (class[2] != 'p')
145 BAND_CASES(DCPY2);
146 else
147 BAND_CASES(DCPY1);
148 } else {
149 PROTECT(x1 = duplicate(x0));
150 if (class[2] != 'p')
151 BAND_CASES(BAND2);
152 else
153 BAND_CASES(BAND1);
154 }
155 } else {
156 if (sy || (tr && (class[1] == 'g' || ul0 == ul1 || n <= 1))) {
157 PROTECT(x1 = duplicate(x0));
158 if (ATTRIB(x1) != R_NilValue) {
159 SET_ATTRIB(x1, R_NilValue);
160 if (OBJECT(x1))
161 SET_OBJECT(x1, 0);
162 }
163 } else {
164 /* Band is "opposite" the stored triangle : */
165 PROTECT(from = dense_transpose(from, class));
166 x1 = GET_SLOT(from, Matrix_xSym);
167 UNPROTECT(1);
168 PROTECT(x1);
169 }
170 if (class[2] != 'p')
171 BAND_CASES(BAND2);
172 else
173 BAND_CASES(BAND1);
174 }
175 SET_SLOT(to, Matrix_xSym, x1);
176
177#undef BAND_CASES
178#undef BAND2
179#undef BAND1
180#undef DCPY2
181#undef DCPY1
182
183 UNPROTECT(3); /* x1, x0, to */
184 return to;
185}
186
187/* band(<denseMatrix>, k1, k2), tri[ul](<denseMatrix>, k) */
188/* NB: argument validation more or less copied by R_sparse_band() */
189SEXP R_dense_band(SEXP from, SEXP k1, SEXP k2)
190{
191 if (!IS_S4_OBJECT(from)) {
192 /* defined in ./coerce.c : */
193 SEXP matrix_as_dense(SEXP, const char *, char, char, int, int);
194 from = matrix_as_dense(from, ".ge", '\0', '\0', 0, 0);
195 }
196 PROTECT(from);
197 static const char *valid[] = { VALID_DENSE, "" };
198 int ivalid = R_check_class_etc(from, valid);
199 if (ivalid < 0)
200 ERROR_INVALID_CLASS(from, __func__);
201
202 SEXP dim = PROTECT(GET_SLOT(from, Matrix_DimSym));
203 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1];
204 UNPROTECT(1);
205
206 int a, b;
207 if (k1 == R_NilValue) // tril()
208 a = -m ; // was (m > 0) ? 1 - m : 0;
209 else if ((a = asInteger(k1)) == NA_INTEGER || a < -m || a > n)
210 error(_("'%s' (%d) must be an integer from %s (%d) to %s (%d)"),
211 "k1", a, "-Dim[1]", -m, "Dim[2]", n);
212 if (k2 == R_NilValue) // triu()
213 b = n; // was (n > 0) ? n - 1 : 0;
214 else if ((b = asInteger(k2)) == NA_INTEGER || b < -m || b > n)
215 error(_("'%s' (%d) must be an integer from %s (%d) to %s (%d)"),
216 "k2", b, "-Dim[1]", -m, "Dim[2]", n);
217 else if (b < a)
218 error(_("'%s' (%d) must be less than or equal to '%s' (%d)"),
219 "k1", a, "k2", b);
220
221 from = dense_band(from, valid[ivalid], a, b);
222 UNPROTECT(1);
223 return from;
224}
225
226SEXP dense_diag_get(SEXP obj, const char *class, int names)
227{
228 SEXP dim = PROTECT(GET_SLOT(obj, Matrix_DimSym));
229 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1], r = (m < n) ? m : n, j;
230 UNPROTECT(1); /* dim */
231
232 char ul = 'U', di = 'N';
233 if (class[1] != 'g') {
234 if (class[2] == 'p') {
235 SEXP uplo = PROTECT(GET_SLOT(obj, Matrix_uploSym));
236 ul = *CHAR(STRING_ELT(uplo, 0));
237 UNPROTECT(1); /* uplo */
238 }
239 if (class[1] == 't') {
240 SEXP diag = PROTECT(GET_SLOT(obj, Matrix_diagSym));
241 di = *CHAR(STRING_ELT(diag, 0));
242 UNPROTECT(1); /* diag */
243 }
244 }
245
246 SEXP x = PROTECT(GET_SLOT(obj, Matrix_xSym)),
247 res = PROTECT(allocVector(TYPEOF(x), r));
248
249#define DG_LOOP(_CTYPE_, _PTR_, _ONE_) \
250 do { \
251 _CTYPE_ *pres = _PTR_(res), *px = _PTR_(x); \
252 if (di == 'U') \
253 for (j = 0; j < r; ++j) \
254 *(pres++) = _ONE_; \
255 else if (class[2] != 'p') { \
256 R_xlen_t m1a = (R_xlen_t) m + 1; \
257 for (j = 0; j < r; ++j, px += m1a) \
258 *(pres++) = *px; \
259 } \
260 else if (ul == 'U') \
261 for (j = 0; j < n; px += (++j) + 1) \
262 *(pres++) = *px; \
263 else \
264 for (j = 0; j < n; px += n - (j++)) \
265 *(pres++) = *px; \
266 } while (0)
267
268 switch (class[0]) {
269 case 'n':
270 case 'l':
271 DG_LOOP(int, LOGICAL, 1);
272 break;
273 case 'i':
274 DG_LOOP(int, INTEGER, 1);
275 break;
276 case 'd':
277 DG_LOOP(double, REAL, 1.0);
278 break;
279 case 'z':
280 DG_LOOP(Rcomplex, COMPLEX, Matrix_zone);
281 break;
282 default:
283 break;
284 }
285
286 if (names) {
287 /* NB: The logic here must be adjusted once the validity method
288 for 'symmetricMatrix' enforces symmetric 'Dimnames'
289 */
290 SEXP dn = PROTECT(GET_SLOT(obj, Matrix_DimNamesSym)),
291 rn = VECTOR_ELT(dn, 0),
292 cn = VECTOR_ELT(dn, 1);
293 if (cn == R_NilValue) {
294 if (class[1] == 's' && rn != R_NilValue)
295 setAttrib(res, R_NamesSymbol, rn);
296 } else {
297 if (class[1] == 's')
298 setAttrib(res, R_NamesSymbol, cn);
299 else if (rn != R_NilValue &&
300 (rn == cn || equal_character_vectors(rn, cn, r)))
301 setAttrib(res, R_NamesSymbol, (r == m) ? rn : cn);
302 }
303 UNPROTECT(1); /* dn */
304 }
305
306#undef DG_LOOP
307
308 UNPROTECT(2); /* x, res */
309 return res;
310}
311
312SEXP R_dense_diag_get(SEXP obj, SEXP names)
313{
314 static const char *valid[] = { VALID_DENSE, "" };
315 int ivalid = R_check_class_etc(obj, valid);
316 if (ivalid < 0)
317 ERROR_INVALID_CLASS(obj, __func__);
318
319 int names_;
320 if (TYPEOF(names) != LGLSXP || LENGTH(names) < 1 ||
321 (names_ = LOGICAL(names)[0]) == NA_LOGICAL)
322 error(_("'%s' must be %s or %s"), "names", "TRUE", "FALSE");
323
324 return dense_diag_get(obj, valid[ivalid], names_);
325}
326
327SEXP dense_diag_set(SEXP from, const char *class, SEXP value, int new)
328{
329 SEXP to = PROTECT(newObject(class));
330 int v = LENGTH(value) != 1;
331
332 SEXP dim = PROTECT(GET_SLOT(from, Matrix_DimSym));
333 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1], r = (m < n) ? m : n, j;
334 if (m != n || n > 0)
335 SET_SLOT(to, Matrix_DimSym, dim);
336 UNPROTECT(1); /* dim */
337
338 SEXP dimnames = PROTECT(GET_SLOT(from, Matrix_DimNamesSym));
339 SET_SLOT(to, Matrix_DimNamesSym, dimnames);
340 UNPROTECT(1); /* dimnames */
341
342 char ul = 'U';
343 if (class[1] != 'g') {
344 SEXP uplo = PROTECT(GET_SLOT(from, Matrix_uploSym));
345 ul = *CHAR(STRING_ELT(uplo, 0));
346 if (ul != 'U')
347 SET_SLOT(to, Matrix_uploSym, uplo);
348 UNPROTECT(1); /* uplo */
349 }
350
351 SEXP x = PROTECT(GET_SLOT(from, Matrix_xSym));
352 if (new) {
353 x = duplicate(x);
354 UNPROTECT(1); /* x */
355 PROTECT(x);
356 }
357 SET_SLOT(to, Matrix_xSym, x);
358
359#define DS_LOOP(_CTYPE_, _PTR_) \
360 do { \
361 _CTYPE_ *px = _PTR_(x), *pvalue = _PTR_(value); \
362 if (class[2] != 'p') { \
363 R_xlen_t m1a = (R_xlen_t) m + 1; \
364 if (v) \
365 for (j = 0; j < r; ++j, px += m1a) \
366 *px = *(pvalue++); \
367 else \
368 for (j = 0; j < r; ++j, px += m1a) \
369 *px = *pvalue; \
370 } else if (ul == 'U') { \
371 if (v) \
372 for (j = 0; j < n; px += (++j) + 1) \
373 *px = *(pvalue++); \
374 else \
375 for (j = 0; j < n; px += (++j) + 1) \
376 *px = *pvalue; \
377 } else { \
378 if (v) \
379 for (j = 0; j < n; px += n - (j++)) \
380 *px = *(pvalue++); \
381 else \
382 for (j = 0; j < n; px += n - (j++)) \
383 *px = *pvalue; \
384 } \
385 } while (0)
386
387 switch (class[0]) {
388 case 'n':
389 case 'l':
390 DS_LOOP(int, LOGICAL);
391 break;
392 case 'i':
393 DS_LOOP(int, INTEGER);
394 break;
395 case 'd':
396 DS_LOOP(double, REAL);
397 break;
398 case 'z':
399 DS_LOOP(Rcomplex, COMPLEX);
400 break;
401 default:
402 break;
403 }
404
405#undef DS_LOOP
406
407 UNPROTECT(2); /* x, to */
408 return to;
409}
410
411SEXP R_dense_diag_set(SEXP from, SEXP value)
412{
413 static const char *valid[] = { VALID_DENSE, "" };
414 int ivalid = R_check_class_etc(from, valid);
415 if (ivalid < 0)
416 ERROR_INVALID_CLASS(from, __func__);
417 const char *class = valid[ivalid];
418
419 SEXPTYPE tx = kindToType(class[0]), tv = TYPEOF(value);
420
421 switch (tv) {
422 case LGLSXP:
423 case INTSXP:
424 case REALSXP:
425 case CPLXSXP:
426 break;
427 default:
428 error(_("replacement diagonal has incompatible type \"%s\""),
429 type2char(tv));
430 break;
431 }
432
433 SEXP dim = GET_SLOT(from, Matrix_DimSym);
434 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1], r = (m < n) ? m : n;
435 R_xlen_t len = XLENGTH(value);
436 if (len != 1 && len != r)
437 error(_("replacement diagonal has wrong length"));
438
439 int new = 1;
440 if (tv <= tx) {
441 PROTECT(from);
442 PROTECT(value = coerceVector(value, tx));
443 } else {
444 /* defined in ./coerce.c : */
445 SEXP dense_as_kind(SEXP, const char *, char, int);
446#ifndef MATRIX_ENABLE_IMATRIX
447 if (tv == INTSXP) {
448 PROTECT(from = dense_as_kind(from, class, 'd', 0));
449 PROTECT(value = coerceVector(value, REALSXP));
450 } else {
451#endif
452 PROTECT(from = dense_as_kind(from, class, typeToKind(tv), 0));
453 PROTECT(value);
454#ifndef MATRIX_ENABLE_IMATRIX
455 }
456#endif
457 class = valid[R_check_class_etc(from, valid)];
458 new = 0;
459 }
460
461 from = dense_diag_set(from, class, value, new);
462 UNPROTECT(2);
463 return from;
464}
465
466SEXP dense_transpose(SEXP from, const char *class)
467{
468 SEXP to = PROTECT(newObject(class));
469
470 SEXP dim = PROTECT(GET_SLOT(from, Matrix_DimSym));
471 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1], i, j;
472 if (m != n) {
473 UNPROTECT(1); /* dim */
474 PROTECT(dim = GET_SLOT(to, Matrix_DimSym));
475 pdim = INTEGER(dim);
476 pdim[0] = n;
477 pdim[1] = m;
478 } else if (n > 0)
479 SET_SLOT(to, Matrix_DimSym, dim);
480 UNPROTECT(1); /* dim */
481
482 SEXP dimnames = PROTECT(GET_SLOT(from, Matrix_DimNamesSym));
483 if (class[1] == 's' || class[1] == 'p' || class[1] == 'o')
484 SET_SLOT(to, Matrix_DimNamesSym, dimnames);
485 else
486 set_reversed_DimNames(to, dimnames);
487 UNPROTECT(1); /* dimnames */
488
489 char ul = 'U';
490 if (class[1] != 'g') {
491 SEXP uplo = PROTECT(GET_SLOT(from, Matrix_uploSym));
492 ul = *CHAR(STRING_ELT(uplo, 0));
493 UNPROTECT(1); /* uplo */
494 if (ul == 'U') {
495 PROTECT(uplo = mkString("L"));
496 SET_SLOT(to, Matrix_uploSym, uplo);
497 UNPROTECT(1); /* uplo */
498 }
499 if (class[1] == 't') {
500 SEXP diag = PROTECT(GET_SLOT(from, Matrix_diagSym));
501 char di = *CHAR(STRING_ELT(diag, 0));
502 if (di != 'N')
503 SET_SLOT(to, Matrix_diagSym, diag);
504 UNPROTECT(1); /* diag */
505 } else {
506 SEXP factors = PROTECT(GET_SLOT(from, Matrix_factorsSym));
507 if (LENGTH(factors) > 0)
508 SET_SLOT(to, Matrix_factorsSym, factors);
509 UNPROTECT(1); /* factors */
510
511 if (class[1] == 'o' && n > 0) {
512 SEXP sd = PROTECT(GET_SLOT(from, Matrix_sdSym));
513 SET_SLOT(to, Matrix_sdSym, sd);
514 UNPROTECT(1); /* sd */
515 }
516 }
517 }
518
519 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym)),
520 x1 = PROTECT(allocVector(TYPEOF(x0), XLENGTH(x0)));
521 SET_SLOT(to, Matrix_xSym, x1);
522
523#define TRANS_LOOP(_CTYPE_, _PTR_) \
524 do { \
525 _CTYPE_ *px0 = _PTR_(x0), *px1 = _PTR_(x1); \
526 if (class[2] != 'p') { \
527 R_xlen_t mn1s = XLENGTH(x0) - 1; \
528 for (j = 0; j < m; ++j, px0 -= mn1s) \
529 for (i = 0; i < n; ++i, px0 += m) \
530 *(px1++) = *px0; \
531 } else if (ul == 'U') { \
532 for (j = 0; j < n; ++j) \
533 for (i = j; i < n; ++i) \
534 *(px1++) = *(px0 + PACKED_AR21_UP(j, i)); \
535 } else { \
536 R_xlen_t n2 = (R_xlen_t) n * 2; \
537 for (j = 0; j < n; ++j) \
538 for (i = 0; i <= j; ++i) \
539 *(px1++) = *(px0 + PACKED_AR21_LO(j, i, n2)); \
540 } \
541 } while (0)
542
543 switch (class[0]) {
544 case 'n':
545 case 'l':
546 TRANS_LOOP(int, LOGICAL);
547 break;
548 case 'i':
549 TRANS_LOOP(int, INTEGER);
550 break;
551 case 'c':
552 case 'd':
553 TRANS_LOOP(double, REAL);
554 break;
555 case 'z':
556 TRANS_LOOP(Rcomplex, COMPLEX);
557 break;
558 default:
559 break;
560 }
561
562#undef TRANS_LOOP
563
564 UNPROTECT(3); /* x1, x0, to */
565 return to;
566}
567
568SEXP R_dense_transpose(SEXP from)
569{
570 static const char *valid[] = {
571 "dpoMatrix", "dppMatrix", "corMatrix", "copMatrix",
572 VALID_DENSE, "" };
573 int ivalid = R_check_class_etc(from, valid);
574 if (ivalid < 0)
575 ERROR_INVALID_CLASS(from, __func__);
576
577 return dense_transpose(from, valid[ivalid]);
578}
579
580SEXP dense_force_symmetric(SEXP from, const char *class, char ul)
581{
582 char ul0 = 'U', ul1 = 'U', di = 'N';
583 if (class[1] != 'g') {
584 SEXP uplo = PROTECT(GET_SLOT(from, Matrix_uploSym));
585 ul0 = ul1 = *CHAR(STRING_ELT(uplo, 0));
586 UNPROTECT(1); /* uplo */
587 if (class[1] == 't') {
588 SEXP diag = PROTECT(GET_SLOT(from, Matrix_diagSym));
589 di = *CHAR(STRING_ELT(diag, 0));
590 UNPROTECT(1); /* diag */
591 }
592 }
593
594 if (ul != '\0')
595 ul1 = ul;
596
597 if (class[1] == 's') {
598 /* .s[yp]Matrix */
599 if (ul0 == ul1)
600 return from;
601 SEXP to = PROTECT(dense_transpose(from, class));
602 if (class[0] == 'z') {
603 /* Need _conjugate_ transpose */
604 SEXP x1 = PROTECT(GET_SLOT(to, Matrix_xSym));
605 conjugate(x1);
606 UNPROTECT(1); /* x1 */
607 }
608 UNPROTECT(1) /* to */;
609 return to;
610 }
611
612 /* Now handling just .(ge|tr|tp)Matrix ... */
613
614 char cl[] = ".s.Matrix";
615 cl[0] = class[0];
616 cl[2] = (class[2] != 'p') ? 'y' : 'p';
617 SEXP to = PROTECT(newObject(cl));
618
619 SEXP dim = PROTECT(GET_SLOT(from, Matrix_DimSym));
620 int *pdim = INTEGER(dim), n = pdim[0];
621 if (pdim[1] != n)
622 error(_("attempt to symmetrize a non-square matrix"));
623 if (n > 0)
624 SET_SLOT(to, Matrix_DimSym, dim);
625 UNPROTECT(1); /* dim */
626
627 SEXP dimnames = PROTECT(GET_SLOT(from, Matrix_DimNamesSym));
628 set_symmetrized_DimNames(to, dimnames, -1);
629 UNPROTECT(1); /* dimnames */
630
631 if (ul1 != 'U') {
632 SEXP uplo = PROTECT(mkString("L"));
633 SET_SLOT(to, Matrix_uploSym, uplo);
634 UNPROTECT(1); /* uplo */
635 }
636
637 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym));
638
639 if (class[1] == 'g' || ul0 == ul1)
640 SET_SLOT(to, Matrix_xSym, x0);
641 else {
642 SEXP x1 = PROTECT(allocVector(TYPEOF(x0), XLENGTH(x0)));
643 SET_SLOT(to, Matrix_xSym, x1);
644
645 R_xlen_t len = XLENGTH(x1);
646
647#define DCPY(_PREFIX_, _CTYPE_, _PTR_) \
648 do { \
649 _CTYPE_ *px0 = _PTR_(x0), *px1 = _PTR_(x1); \
650 Matrix_memset(px1, 0, len, sizeof(_CTYPE_)); \
651 if (class[2] != 'p') \
652 _PREFIX_ ## dcpy2(px1, px0, n, len, '\0', di); \
653 else \
654 _PREFIX_ ## dcpy1(px1, px0, n, len, ul1, ul0, di); \
655 } while (0)
656
657 switch (class[0]) {
658 case 'n':
659 case 'l':
660 DCPY(i, int, LOGICAL);
661 break;
662 case 'i':
663 DCPY(i, int, INTEGER);
664 break;
665 case 'd':
666 DCPY(d, double, REAL);
667 break;
668 case 'z':
669 DCPY(z, Rcomplex, COMPLEX);
670 break;
671 default:
672 break;
673 }
674
675#undef DCPY
676
677 UNPROTECT(1); /* x1 */
678 }
679
680 UNPROTECT(2); /* x0, to */
681 return to;
682}
683
684SEXP R_dense_force_symmetric(SEXP from, SEXP uplo)
685{
686 static const char *valid[] = { VALID_DENSE, "" };
687 int ivalid = R_check_class_etc(from, valid);
688 if (ivalid < 0)
689 ERROR_INVALID_CLASS(from, __func__);
690
691 char ul = '\0';
692 if (uplo != R_NilValue) {
693 if (TYPEOF(uplo) != STRSXP || LENGTH(uplo) < 1 ||
694 (uplo = STRING_ELT(uplo, 0)) == NA_STRING ||
695 ((ul = *CHAR(uplo)) != 'U' && ul != 'L'))
696 error(_("invalid '%s' to '%s'"), "uplo", __func__);
697 }
698
699 return dense_force_symmetric(from, valid[ivalid], ul);
700}
701
702SEXP dense_symmpart(SEXP from, const char *class)
703{
704 if (class[0] != 'z' && class[0] != 'd') {
705 /* defined in ./coerce.c : */
706 SEXP dense_as_kind(SEXP, const char *, char, int);
707 from = dense_as_kind(from, class, 'd', 0);
708 }
709 if (class[0] != 'z' && class[1] == 's')
710 return from;
711 PROTECT(from);
712
713 char cl[] = ".s.Matrix";
714 cl[0] = (class[0] != 'z') ? 'd' : 'z';
715 cl[2] = (class[2] != 'p') ? 'y' : 'p';
716 SEXP to = PROTECT(newObject(cl));
717
718 SEXP dim = PROTECT(GET_SLOT(from, Matrix_DimSym));
719 int *pdim = INTEGER(dim), n = pdim[0];
720 if (pdim[1] != n)
721 error(_("attempt to get symmetric part of non-square matrix"));
722 if (n > 0)
723 SET_SLOT(to, Matrix_DimSym, dim);
724 UNPROTECT(1); /* dim */
725
726 SEXP dimnames = PROTECT(GET_SLOT(from, Matrix_DimNamesSym));
727 if (class[1] == 's')
728 SET_SLOT(to, Matrix_DimNamesSym, dimnames);
729 else
730 set_symmetrized_DimNames(to, dimnames, -1);
731 UNPROTECT(1); /* dimnames */
732
733 char ul = 'U', di = 'N';
734 if (class[1] != 'g') {
735 SEXP uplo = PROTECT(GET_SLOT(from, Matrix_uploSym));
736 ul = *CHAR(STRING_ELT(uplo, 0));
737 if (ul != 'U')
738 SET_SLOT(to, Matrix_uploSym, uplo);
739 UNPROTECT(1); /* uplo */
740 if (class[1] == 't') {
741 SEXP diag = PROTECT(GET_SLOT(from, Matrix_diagSym));
742 di = *CHAR(STRING_ELT(diag, 0));
743 UNPROTECT(1); /* diag */
744 }
745 }
746
747 SEXP x = PROTECT(GET_SLOT(from, Matrix_xSym));
748 if (class[0] == 'z' || class[0] == 'd') {
749 x = duplicate(x);
750 UNPROTECT(1); /* x */
751 PROTECT(x);
752 }
753 SET_SLOT(to, Matrix_xSym, x);
754
755 if (class[1] == 's') {
756 /* Symmetric part of Hermitian matrix is real part */
757 zeroIm(x);
758 UNPROTECT(3); /* x, to, from */
759 return to;
760 }
761
762 int i, j;
763
764#define SP_LOOP(_CTYPE_, _PTR_, _INCREMENT_, _SCALE1_, _ONE_) \
765 do { \
766 _CTYPE_ *px = _PTR_(x); \
767 if (class[1] == 'g') { \
768 _CTYPE_ *py = px; \
769 for (j = 0; j < n; ++j) { \
770 for (i = j + 1; i < n; ++i) { \
771 px += n; \
772 py += 1; \
773 _INCREMENT_((*px), (*py)); \
774 _SCALE1_((*px), 0.5); \
775 } \
776 px = (py += j + 2); \
777 } \
778 } else if (class[2] != 'p') { \
779 if (ul == 'U') { \
780 for (j = 0; j < n; ++j) { \
781 for (i = 0; i < j; ++i) { \
782 _SCALE1_((*px), 0.5); \
783 px += 1; \
784 } \
785 px += n - j; \
786 } \
787 } else { \
788 for (j = 0; j < n; ++j) { \
789 px += j + 1; \
790 for (i = j + 1; i < n; ++i) { \
791 _SCALE1_((*px), 0.5); \
792 px += 1; \
793 } \
794 } \
795 } \
796 if (di != 'N') { \
797 R_xlen_t n1a = (R_xlen_t) n + 1; \
798 px = _PTR_(x); \
799 for (j = 0; j < n; ++j, px += n1a) \
800 *px = _ONE_; \
801 } \
802 } else { \
803 if (ul == 'U') { \
804 for (j = 0; j < n; ++j) { \
805 for (i = 0; i < j; ++i) { \
806 _SCALE1_((*px), 0.5); \
807 px += 1; \
808 } \
809 px += 1; \
810 } \
811 if (di != 'N') { \
812 px = _PTR_(x); \
813 for (j = 0; j < n; px += (++j) + 1) \
814 *px = _ONE_; \
815 } \
816 } else { \
817 for (j = 0; j < n; ++j) { \
818 px += 1; \
819 for (i = j + 1; i < n; ++i) { \
820 _SCALE1_((*px), 0.5); \
821 px += 1; \
822 } \
823 } \
824 if (di != 'N') { \
825 px = _PTR_(x); \
826 for (j = 0; j < n; px += n - (j++)) \
827 *px = _ONE_; \
828 } \
829 } \
830 } \
831 } while (0)
832
833 if (cl[0] == 'd')
834 SP_LOOP(double, REAL, INCREMENT_REAL, SCALE1_REAL, 1.0);
835 else
836 SP_LOOP(Rcomplex, COMPLEX, INCREMENT_COMPLEX, SCALE1_COMPLEX, Matrix_zone);
837
838#undef SP_LOOP
839
840 UNPROTECT(3); /* x, to, from */
841 return to;
842}
843
844SEXP R_dense_symmpart(SEXP from)
845{
846 static const char *valid[] = { VALID_DENSE, "" };
847 int ivalid = R_check_class_etc(from, valid);
848 if (ivalid < 0)
849 ERROR_INVALID_CLASS(from, __func__);
850
851 return dense_symmpart(from, valid[ivalid]);
852}
853
854SEXP dense_skewpart(SEXP from, const char *class)
855{
856 if (class[0] != 'z' && class[0] != 'd') {
857 /* defined in ./coerce.c : */
858 SEXP dense_as_kind(SEXP, const char *, char, int);
859 from = dense_as_kind(from, class, 'd', 0);
860 }
861 PROTECT(from);
862
863 char cl[] = "...Matrix";
864 cl[0] = (class[0] != 'z') ? 'd' : 'z';
865 cl[1] = (class[1] != 's') ? 'g' : 's';
866 cl[2] = (class[1] != 's') ? 'e' :
867 ((class[0] != 'z') ? 'C' : ((class[2] != 'p') ? 'y' : 'p'));
868 SEXP to = PROTECT(newObject(cl));
869
870 SEXP dim = PROTECT(GET_SLOT(from, Matrix_DimSym));
871 int *pdim = INTEGER(dim), n = pdim[0];
872 if (pdim[1] != n)
873 error(_("attempt to get skew-symmetric part of non-square matrix"));
874 if (n > 0)
875 SET_SLOT(to, Matrix_DimSym, dim);
876 UNPROTECT(1); /* dim */
877
878 SEXP dimnames = PROTECT(GET_SLOT(from, Matrix_DimNamesSym));
879 if (class[1] == 's')
880 SET_SLOT(to, Matrix_DimNamesSym, dimnames);
881 else
882 set_symmetrized_DimNames(to, dimnames, -1);
883 UNPROTECT(1); /* dimnames */
884
885 char ul = 'U';
886 if (class[1] != 'g') {
887 SEXP uplo = PROTECT(GET_SLOT(from, Matrix_uploSym));
888 ul = *CHAR(STRING_ELT(uplo, 0));
889 if (class[1] == 's' && ul != 'U')
890 SET_SLOT(to, Matrix_uploSym, uplo);
891 UNPROTECT(1); /* uplo */
892 }
893
894 if (class[1] == 's' && class[0] != 'z') {
895 /* Skew-symmetric part of Hermitian matrix is imaginary part */
896 SEXP p = PROTECT(allocVector(INTSXP, (R_xlen_t) n + 1));
897 int *pp = INTEGER(p);
898 Matrix_memset(pp, 0, (R_xlen_t) n + 1, sizeof(int));
899 SET_SLOT(to, Matrix_pSym, p);
900 UNPROTECT(3); /* p, to, from */
901 return to;
902 }
903
904 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym)), x1 = x0;
905
906 if (class[1] == 's') {
907 /* Skew-symmetric part of Hermitian matrix is imaginary part */
908 x1 = duplicate(x1);
909 UNPROTECT(1); /* x1 */
910 PROTECT(x1);
911 SET_SLOT(to, Matrix_xSym, x1);
912 zeroRe(x1);
913 UNPROTECT(3); /* x1, to, from */
914 return to;
915 }
916
917 if (class[2] == 'p' || class[0] == 'z' || class[0] == 'd') {
918 if ((Matrix_int_fast64_t) n * n > R_XLEN_T_MAX)
919 error(_("attempt to allocate vector of length exceeding %s"),
920 "R_XLEN_T_MAX");
921 x1 = allocVector(TYPEOF(x0), (R_xlen_t) n * n);
922 }
923 PROTECT(x1);
924 SET_SLOT(to, Matrix_xSym, x1);
925
926 int i, j;
927 R_xlen_t upos = 0, lpos = 0;
928
929#define SP_LOOP(_CTYPE_, _PTR_, _INCREMENT_, _ASSIGN_, _ZERO_) \
930 do { \
931 _CTYPE_ *px0 = _PTR_(x0), *px1 = _PTR_(x1); \
932 if (class[1] == 'g') { \
933 for (j = 0; j < n; ++j) { \
934 lpos = j; \
935 for (i = 0; i < j; ++i) { \
936 _ASSIGN_(px1[upos], 0.5 * px0[upos]); \
937 _INCREMENT_(px1[upos], -0.5 * px0[lpos]); \
938 _ASSIGN_(px1[lpos], -px1[upos]); \
939 upos += 1; \
940 lpos += n; \
941 } \
942 px1[upos] = _ZERO_; \
943 upos += n - j; \
944 } \
945 } else if (class[2] != 'p') { \
946 if (ul == 'U') { \
947 for (j = 0; j < n; ++j) { \
948 lpos = j; \
949 for (i = 0; i < j; ++i) { \
950 _ASSIGN_(px1[upos], 0.5 * px0[upos]); \
951 _ASSIGN_(px1[lpos], -px1[upos]); \
952 upos += 1; \
953 lpos += n; \
954 } \
955 px1[upos] = _ZERO_; \
956 upos += n - j; \
957 } \
958 } else { \
959 for (j = 0; j < n; ++j) { \
960 upos = lpos; \
961 px1[lpos] = _ZERO_; \
962 for (i = j + 1; i < n; ++i) { \
963 upos += n; \
964 lpos += 1; \
965 _ASSIGN_(px1[lpos], 0.5 * px0[lpos]); \
966 _ASSIGN_(px1[upos], -px1[lpos]); \
967 } \
968 lpos += j + 2; \
969 } \
970 } \
971 } else { \
972 if (ul == 'U') { \
973 for (j = 0; j < n; ++j, ++px0) { \
974 lpos = j; \
975 for (i = 0; i < j; ++i, ++px0) { \
976 _ASSIGN_(px1[upos], 0.5 * (*px0)); \
977 _ASSIGN_(px1[lpos], -px1[upos]); \
978 upos += 1; \
979 lpos += n; \
980 } \
981 px1[upos] = _ZERO_; \
982 upos += n - j; \
983 } \
984 } else { \
985 for (j = 0; j < n; ++j, ++px0) { \
986 upos = lpos; \
987 px1[lpos] = _ZERO_; \
988 for (i = j + 1; i < n; ++i, ++px0) { \
989 upos += n; \
990 lpos += 1; \
991 _ASSIGN_(px1[lpos], 0.5 * (*px0)); \
992 _ASSIGN_(px1[upos], -px1[lpos]); \
993 } \
994 lpos += j + 2; \
995 } \
996 } \
997 } \
998 } while (0)
999
1000 if (cl[0] == 'd')
1001 SP_LOOP(double, REAL, INCREMENT_REAL, ASSIGN_REAL, 0.0);
1002 else
1003 SP_LOOP(Rcomplex, COMPLEX, INCREMENT_COMPLEX, ASSIGN_COMPLEX, Matrix_zzero);
1004
1005#undef SP_LOOP
1006
1007 UNPROTECT(4); /* x1, x0, to, from */
1008 return to;
1009}
1010
1011SEXP R_dense_skewpart(SEXP from)
1012{
1013 static const char *valid[] = { VALID_DENSE, "" };
1014 int ivalid = R_check_class_etc(from, valid);
1015 if (ivalid < 0)
1016 ERROR_INVALID_CLASS(from, __func__);
1017
1018 return dense_skewpart(from, valid[ivalid]);
1019}
1020
1021int dense_is_symmetric(SEXP obj, const char *class, int checkDN)
1022{
1023 if (class[1] == 's')
1024 return 1;
1025
1026 if (checkDN) {
1027 SEXP dimnames = GET_SLOT(obj, Matrix_DimNamesSym);
1028 if (!DimNames_is_symmetric(dimnames))
1029 return 0;
1030 }
1031
1032 if (class[1] == 't')
1033 return dense_is_diagonal(obj, class);
1034
1035 SEXP dim = GET_SLOT(obj, Matrix_DimSym);
1036 int *pdim = INTEGER(dim), n = pdim[0];
1037 if (pdim[1] != n)
1038 return 0;
1039 if (n <= 1)
1040 return 1;
1041
1042 SEXP x = GET_SLOT(obj, Matrix_xSym);
1043 int i, j;
1044
1045#define IS_LOOP(_CTYPE_, _PTR_, _NOTREAL_, _NOTCONJ_) \
1046 do { \
1047 _CTYPE_ *px = _PTR_(x), *py = px; \
1048 for (j = 0; j < n; px = (py += (++j) + 1)) { \
1049 if (_NOTREAL_((*px))) \
1050 return 0; \
1051 for (i = j + 1; i < n; ++i) { \
1052 px += n; \
1053 py += 1; \
1054 if (_NOTCONJ_((*px), (*py))) \
1055 return 0; \
1056 } \
1057 } \
1058 return 1; \
1059 } while (0)
1060
1061 switch (class[0]) {
1062 case 'n':
1063 IS_LOOP(int, LOGICAL, NOTREAL_PATTERN, NOTCONJ_PATTERN);
1064 break;
1065 case 'l':
1066 IS_LOOP(int, LOGICAL, NOTREAL_LOGICAL, NOTCONJ_LOGICAL);
1067 break;
1068 case 'i':
1069 IS_LOOP(int, INTEGER, NOTREAL_INTEGER, NOTCONJ_INTEGER);
1070 break;
1071 case 'd':
1072 IS_LOOP(double, REAL, NOTREAL_REAL, NOTCONJ_REAL);
1073 break;
1074 case 'z':
1075 IS_LOOP(Rcomplex, COMPLEX, NOTREAL_COMPLEX, NOTCONJ_COMPLEX);
1076 break;
1077 default:
1078 break;
1079 }
1080
1081#undef IS_LOOP
1082
1083 return 0;
1084}
1085
1086SEXP R_dense_is_symmetric(SEXP obj, SEXP checkDN)
1087{
1088 if (!IS_S4_OBJECT(obj)) {
1089 /* defined in ./coerce.c : */
1090 SEXP matrix_as_dense(SEXP, const char *, char, char, int, int);
1091 obj = matrix_as_dense(obj, ".ge", '\0', '\0', 0, 0);
1092 }
1093 PROTECT(obj);
1094 static const char *valid[] = { VALID_DENSE, "" };
1095 int ivalid = R_check_class_etc(obj, valid);
1096 if (ivalid < 0)
1097 ERROR_INVALID_CLASS(obj, __func__);
1098
1099 int checkDN_;
1100 if (TYPEOF(checkDN) != LGLSXP || LENGTH(checkDN) < 1 ||
1101 (checkDN_ = LOGICAL(checkDN)[0]) == NA_LOGICAL)
1102 error(_("'%s' must be %s or %s"), "checkDN", "TRUE", "FALSE");
1103
1104 SEXP ans = ScalarLogical(dense_is_symmetric(obj, valid[ivalid], checkDN_));
1105 UNPROTECT(1);
1106 return ans;
1107}
1108
1109int dense_is_triangular(SEXP obj, const char *class, int upper)
1110{
1111 if (class[1] == 't') {
1112 SEXP uplo = GET_SLOT(obj, Matrix_uploSym);
1113 char ul = *CHAR(STRING_ELT(uplo, 0));
1114 if (upper == NA_LOGICAL || (upper != 0) == (ul == 'U'))
1115 return (ul == 'U') ? 1 : -1;
1116 else if (dense_is_diagonal(obj, class))
1117 return (ul == 'U') ? -1 : 1;
1118 else
1119 return 0;
1120 }
1121
1122 if (class[1] == 's') {
1123 if (!dense_is_diagonal(obj, class))
1124 return 0;
1125 SEXP uplo = GET_SLOT(obj, Matrix_uploSym);
1126 char ul = *CHAR(STRING_ELT(uplo, 0));
1127 if (upper == NA_LOGICAL)
1128 return (ul == 'U') ? 1 : -1;
1129 else
1130 return (upper != 0) ? 1 : -1;
1131 }
1132
1133 SEXP dim = GET_SLOT(obj, Matrix_DimSym);
1134 int *pdim = INTEGER(dim), n = pdim[0];
1135 if (pdim[1] != n)
1136 return 0;
1137 if (n <= 1)
1138 return (upper != 0) ? 1 : -1;
1139
1140 SEXP x = GET_SLOT(obj, Matrix_xSym);
1141 int i, j;
1142
1143#define IT_LOOP(_CTYPE_, _PTR_, _ISNZ_) \
1144 do { \
1145 _CTYPE_ *px; \
1146 if (upper == NA_LOGICAL) { \
1147 px = _PTR_(x); \
1148 for (j = 0; j < n; px += (++j)) { \
1149 px += 1; \
1150 for (i = j + 1; i < n; ++i, px += 1) { \
1151 if (_ISNZ_(*px)) { \
1152 j = n; \
1153 break; \
1154 } \
1155 } \
1156 } \
1157 if (j == n) \
1158 return 1; \
1159 px = _PTR_(x); \
1160 for (j = 0; j < n; px += n - (++j)) { \
1161 for (i = 0; i < j; ++i, px += 1) { \
1162 if (_ISNZ_(*px)) { \
1163 j = n; \
1164 break; \
1165 } \
1166 } \
1167 px += 1; \
1168 } \
1169 if (j == n) \
1170 return -1; \
1171 return 0; \
1172 } else if (upper != 0) { \
1173 px = _PTR_(x); \
1174 for (j = 0; j < n; px += (++j)) { \
1175 px += 1; \
1176 for (i = j + 1; i < n; ++i, px += 1) \
1177 if (_ISNZ_(*px)) \
1178 return 0; \
1179 } \
1180 return 1; \
1181 } else { \
1182 px = _PTR_(x); \
1183 for (j = 0; j < n; px += n - (++j)) { \
1184 for (i = 0; i < j; ++i, px += 1) \
1185 if (_ISNZ_(*px)) \
1186 return 0; \
1187 px += 1; \
1188 } \
1189 return -1; \
1190 } \
1191 } while (0)
1192
1193 switch (class[0]) {
1194 case 'n':
1195 IT_LOOP(int, LOGICAL, ISNZ_PATTERN);
1196 break;
1197 case 'l':
1198 IT_LOOP(int, LOGICAL, ISNZ_LOGICAL);
1199 break;
1200 case 'i':
1201 IT_LOOP(int, INTEGER, ISNZ_INTEGER);
1202 break;
1203 case 'd':
1204 IT_LOOP(double, REAL, ISNZ_REAL);
1205 break;
1206 case 'z':
1207 IT_LOOP(Rcomplex, COMPLEX, ISNZ_COMPLEX);
1208 break;
1209 default:
1210 break;
1211 }
1212
1213#undef IT_LOOP
1214
1215 return 0;
1216}
1217
1218SEXP R_dense_is_triangular(SEXP obj, SEXP upper)
1219{
1220 if (!IS_S4_OBJECT(obj)) {
1221 /* defined in ./coerce.c : */
1222 SEXP matrix_as_dense(SEXP, const char *, char, char, int, int);
1223 obj = matrix_as_dense(obj, ".ge", '\0', '\0', 0, 0);
1224 }
1225 PROTECT(obj);
1226 static const char *valid[] = { VALID_DENSE, "" };
1227 int ivalid = R_check_class_etc(obj, valid);
1228 if (ivalid < 0)
1229 ERROR_INVALID_CLASS(obj, __func__);
1230
1231 if (TYPEOF(upper) != LGLSXP || LENGTH(upper) < 1)
1232 error(_("'%s' must be %s or %s or %s"), "upper", "TRUE", "FALSE", "NA");
1233 int upper_ = LOGICAL(upper)[0];
1234
1235 int ans_ = dense_is_triangular(obj, valid[ivalid], upper_);
1236 SEXP ans = allocVector(LGLSXP, 1);
1237 LOGICAL(ans)[0] = ans_ != 0;
1238 if (upper_ == NA_LOGICAL && ans_ != 0) {
1239 PROTECT(ans);
1240 static
1241 SEXP kindSym = NULL;
1242 SEXP kindVal = PROTECT(mkString((ans_ > 0) ? "U" : "L"));
1243 if (!kindSym) kindSym = install("kind");
1244 setAttrib(ans, kindSym, kindVal);
1245 UNPROTECT(2);
1246 }
1247 UNPROTECT(1);
1248 return ans;
1249}
1250
1251int dense_is_diagonal(SEXP obj, const char *class)
1252{
1253 SEXP dim = GET_SLOT(obj, Matrix_DimSym);
1254 int *pdim = INTEGER(dim), n = pdim[0];
1255 if (pdim[1] != n)
1256 return 0;
1257 if (n <= 1)
1258 return 1;
1259
1260 char ul = 'U';
1261 if (class[1] != 'g') {
1262 SEXP uplo = GET_SLOT(obj, Matrix_uploSym);
1263 ul = *CHAR(STRING_ELT(uplo, 0));
1264 }
1265
1266 SEXP x = GET_SLOT(obj, Matrix_xSym);
1267 int i, j;
1268
1269#define ID_LOOP(_CTYPE_, _PTR_, _ISNZ_) \
1270 do { \
1271 _CTYPE_ *px = _PTR_(x); \
1272 if (class[1] == 'g') { \
1273 for (j = 0; j < n; ++j) { \
1274 for (i = 0; i < j; ++i) { \
1275 if (_ISNZ_(*px)) \
1276 return 0; \
1277 px += 1; \
1278 } \
1279 px += 1; \
1280 for (i = j + 1; i < n; ++i) { \
1281 if (_ISNZ_(*px)) \
1282 return 0; \
1283 px += 1; \
1284 } \
1285 } \
1286 } else if (class[2] != 'p') { \
1287 if (ul == 'U') { \
1288 for (j = 0; j < n; ++j) { \
1289 for (i = 0; i < j; ++i) { \
1290 if (_ISNZ_(*px)) \
1291 return 0; \
1292 px += 1; \
1293 } \
1294 px += n - j; \
1295 } \
1296 } else { \
1297 for (j = 0; j < n; ++j) { \
1298 px += j + 1; \
1299 for (i = j + 1; i < n; ++i) { \
1300 if (_ISNZ_(*px)) \
1301 return 0; \
1302 px += 1; \
1303 } \
1304 } \
1305 } \
1306 } else { \
1307 if (ul == 'U') { \
1308 for (j = 0; j < n; ++j) { \
1309 for (i = 0; i < j; ++i) { \
1310 if (_ISNZ_(*px)) \
1311 return 0; \
1312 px += 1; \
1313 } \
1314 px += 1; \
1315 } \
1316 } else { \
1317 for (j = 0; j < n; ++j) { \
1318 px += 1; \
1319 for (i = j + 1; i < n; ++i) { \
1320 if (_ISNZ_(*px)) \
1321 return 0; \
1322 px += 1; \
1323 } \
1324 } \
1325 } \
1326 } \
1327 return 1; \
1328 } while (0)
1329
1330 switch (class[0]) {
1331 case 'n':
1332 ID_LOOP(int, LOGICAL, ISNZ_PATTERN);
1333 break;
1334 case 'l':
1335 ID_LOOP(int, LOGICAL, ISNZ_LOGICAL);
1336 break;
1337 case 'i':
1338 ID_LOOP(int, INTEGER, ISNZ_INTEGER);
1339 break;
1340 case 'd':
1341 ID_LOOP(double, REAL, ISNZ_REAL);
1342 break;
1343 case 'z':
1344 ID_LOOP(Rcomplex, COMPLEX, ISNZ_COMPLEX);
1345 break;
1346 default:
1347 break;
1348 }
1349
1350#undef ID_LOOP
1351
1352 return 0;
1353}
1354
1355SEXP R_dense_is_diagonal(SEXP obj)
1356{
1357 if (!IS_S4_OBJECT(obj)) {
1358 /* defined in ./coerce.c : */
1359 SEXP matrix_as_dense(SEXP, const char *, char, char, int, int);
1360 obj = matrix_as_dense(obj, ".ge", '\0', '\0', 0, 0);
1361 }
1362 PROTECT(obj);
1363 static const char *valid[] = { VALID_DENSE, "" };
1364 int ivalid = R_check_class_etc(obj, valid);
1365 if (ivalid < 0)
1366 ERROR_INVALID_CLASS(obj, __func__);
1367
1368 SEXP ans = ScalarLogical(dense_is_diagonal(obj, valid[ivalid]));
1369 UNPROTECT(1);
1370 return ans;
1371}
1372
1373#define CAST_PATTERN(_X_) (_X_ != 0)
1374#define CAST_LOGICAL(_X_) (_X_ != 0)
1375#define CAST_INTEGER(_X_) _X_
1376#define CAST_REAL(_X_) _X_
1377#define CAST_COMPLEX(_X_) _X_
1378
1379#define SUM_CASES \
1380do { \
1381 switch (class[0]) { \
1382 case 'n': \
1383 if (mean) \
1384 SUM_LOOP(int, LOGICAL, double, REAL, \
1385 0.0, 1.0, NA_REAL, ISNA_PATTERN, \
1386 CAST_PATTERN, INCREMENT_REAL, SCALE2_REAL); \
1387 else \
1388 SUM_LOOP(int, LOGICAL, int, INTEGER, \
1389 0, 1, NA_INTEGER, ISNA_PATTERN, \
1390 CAST_PATTERN, INCREMENT_INTEGER, SCALE2_REAL); \
1391 break; \
1392 case 'l': \
1393 if (mean) \
1394 SUM_LOOP(int, LOGICAL, double, REAL, \
1395 0.0, 1.0, NA_REAL, ISNA_LOGICAL, \
1396 CAST_LOGICAL, INCREMENT_REAL, SCALE2_REAL); \
1397 else \
1398 SUM_LOOP(int, LOGICAL, int, INTEGER, \
1399 0, 1, NA_INTEGER, ISNA_LOGICAL, \
1400 CAST_LOGICAL, INCREMENT_INTEGER, SCALE2_REAL); \
1401 break; \
1402 case 'i': \
1403 SUM_LOOP(int, INTEGER, double, REAL, \
1404 0.0, 1.0, NA_REAL, ISNA_INTEGER, \
1405 CAST_INTEGER, INCREMENT_REAL, SCALE2_REAL); \
1406 break; \
1407 case 'd': \
1408 SUM_LOOP(double, REAL, double, REAL, \
1409 0.0, 1.0, NA_REAL, ISNA_REAL, \
1410 CAST_REAL, INCREMENT_REAL, SCALE2_REAL); \
1411 break; \
1412 case 'z': \
1413 SUM_LOOP(Rcomplex, COMPLEX, Rcomplex, COMPLEX, \
1414 Matrix_zzero, Matrix_zone, Matrix_zna, ISNA_COMPLEX, \
1415 CAST_COMPLEX, INCREMENT_COMPLEX, SCALE2_COMPLEX); \
1416 break; \
1417 default: \
1418 break; \
1419 } \
1420} while (0)
1421
1422#define SUM_TYPEOF(c) (c == 'z') ? CPLXSXP : ((mean || c == 'd' || c == 'i') ? REALSXP : INTSXP)
1423
1424static
1425void dense_colsum(SEXP x, const char *class,
1426 int m, int n, char ul, char di, int narm, int mean,
1427 SEXP res)
1428{
1429 int i, j, count = -1, narm_ = narm && mean && class[0] != 'n',
1430 unpacked = class[2] != 'p';
1431
1432#define SUM_LOOP(_CTYPE0_, _PTR0_, _CTYPE1_, _PTR1_, \
1433 _ZERO_, _ONE_, _NA_, _ISNA_, \
1434 _CAST_, _INCREMENT_, _SCALE2_) \
1435 do { \
1436 _CTYPE0_ *px0 = _PTR0_( x); \
1437 _CTYPE1_ *px1 = _PTR1_(res), tmp; \
1438 if (class[1] == 'g') { \
1439 for (j = 0; j < n; ++j) { \
1440 *px1 = _ZERO_; \
1441 SUM_KERNEL(for (i = 0; i < m; ++i), _NA_, _ISNA_, \
1442 _CAST_, _INCREMENT_, _SCALE2_); \
1443 px1 += 1; \
1444 } \
1445 } else if (di == 'N') { \
1446 if (ul == 'U') { \
1447 for (j = 0; j < n; ++j) { \
1448 *px1 = _ZERO_; \
1449 SUM_KERNEL(for (i = 0; i <= j; ++i), _NA_, _ISNA_, \
1450 _CAST_, _INCREMENT_, _SCALE2_); \
1451 if (unpacked) \
1452 px0 += n - j - 1; \
1453 px1 += 1; \
1454 } \
1455 } else { \
1456 for (j = 0; j < n; ++j) { \
1457 if (unpacked) \
1458 px0 += j; \
1459 *px1 = _ZERO_; \
1460 SUM_KERNEL(for (i = j; i < n; ++i), _NA_, _ISNA_, \
1461 _CAST_, _INCREMENT_, _SCALE2_); \
1462 px1 += 1; \
1463 } \
1464 } \
1465 } else { \
1466 if (ul == 'U') { \
1467 for (j = 0; j < n; ++j) { \
1468 *px1 = _ONE_; \
1469 SUM_KERNEL(for (i = 0; i < j; ++i), _NA_, _ISNA_, \
1470 _CAST_, _INCREMENT_, _SCALE2_); \
1471 ++px0; \
1472 if (unpacked) \
1473 px0 += n - j - 1; \
1474 px1 += 1; \
1475 } \
1476 } else { \
1477 for (j = 0; j < n; ++j) { \
1478 if (unpacked) \
1479 px0 += j; \
1480 ++px0; \
1481 *px1 = _ONE_; \
1482 SUM_KERNEL(for (i = j + 1; i < n; ++i), _NA_, _ISNA_, \
1483 _CAST_, _INCREMENT_, _SCALE2_); \
1484 px1 += 1; \
1485 } \
1486 } \
1487 } \
1488 } while (0)
1489
1490#define SUM_KERNEL(_FOR_, _NA_, _ISNA_, _CAST_, _INCREMENT_, _SCALE2_) \
1491 do { \
1492 if (mean) \
1493 count = m; \
1494 _FOR_ { \
1495 if (_ISNA_(*px0)) { \
1496 if (!narm) \
1497 *px1 = _NA_; \
1498 else if (narm_) \
1499 --count; \
1500 } else { \
1501 tmp = _CAST_(*px0); \
1502 _INCREMENT_((*px1), tmp); \
1503 } \
1504 ++px0; \
1505 } \
1506 if (mean) \
1507 _SCALE2_((*px1), count); \
1508 } while (0)
1509
1510 SUM_CASES;
1511
1512#undef SUM_LOOP
1513#undef SUM_KERNEL
1514
1515 return;
1516}
1517
1518static
1519void dense_rowsum(SEXP x, const char *class,
1520 int m, int n, char ul, char di, int narm, int mean,
1521 SEXP res)
1522{
1523 int i, j, *count = NULL, narm_ = narm && mean && class[0] != 'n',
1524 unpacked = class[2] != 'p', symmetric = class[1] == 's';
1525 if (narm_) {
1526 Matrix_Calloc(count, m, int);
1527 for (i = 0; i < m; ++i)
1528 count[i] = n;
1529 }
1530
1531#define SUM_LOOP(_CTYPE0_, _PTR0_, _CTYPE1_, _PTR1_, \
1532 _ZERO_, _ONE_, _NA_, _ISNA_, \
1533 _CAST_, _INCREMENT_, _SCALE2_) \
1534 do { \
1535 _CTYPE0_ *px0 = _PTR0_( x); \
1536 _CTYPE1_ *px1 = _PTR1_(res), tmp = (di == 'N') ? _ZERO_ : _ONE_; \
1537 for (i = 0; i < m; ++i) \
1538 px1[i] = tmp; \
1539 if (class[1] == 'g') { \
1540 for (j = 0; j < n; ++j) \
1541 SUM_KERNEL(for (i = 0; i < m; ++i), _NA_, _ISNA_, \
1542 _CAST_, _INCREMENT_); \
1543 } else if (class[1] == 's' || di == 'N') { \
1544 if (ul == 'U') { \
1545 for (j = 0; j < n; ++j) { \
1546 SUM_KERNEL(for (i = 0; i <= j; ++i), _NA_, _ISNA_, \
1547 _CAST_, _INCREMENT_); \
1548 if (unpacked) \
1549 px0 += n - j - 1; \
1550 } \
1551 } else { \
1552 for (j = 0; j < n; ++j) { \
1553 if (unpacked) \
1554 px0 += j; \
1555 SUM_KERNEL(for (i = j; i < n; ++i), _NA_, _ISNA_, \
1556 _CAST_, _INCREMENT_); \
1557 } \
1558 } \
1559 } else { \
1560 if (ul == 'U') { \
1561 for (j = 0; j < n; ++j) { \
1562 SUM_KERNEL(for (i = 0; i < j; ++i), _NA_, _ISNA_, \
1563 _CAST_, _INCREMENT_); \
1564 ++px0; \
1565 if (unpacked) \
1566 px0 += n - j - 1; \
1567 } \
1568 } else { \
1569 for (j = 0; j < n; ++j) { \
1570 if (unpacked) \
1571 px0 += j; \
1572 ++px0; \
1573 SUM_KERNEL(for (i = j + 1; i < n; ++i), _NA_, _ISNA_, \
1574 _CAST_, _INCREMENT_); \
1575 } \
1576 } \
1577 } \
1578 if (mean) { \
1579 if (narm_) \
1580 for (i = 0; i < m; ++i) \
1581 _SCALE2_(px1[i], count[i]); \
1582 else \
1583 for (i = 0; i < m; ++i) \
1584 _SCALE2_(px1[i], n); \
1585 } \
1586 } while (0)
1587
1588#define SUM_KERNEL(_FOR_, _NA_, _ISNA_, _CAST_, _INCREMENT_) \
1589 do { \
1590 _FOR_ { \
1591 int again = symmetric && i != j; \
1592 if (_ISNA_(*px0)) { \
1593 if (!narm) { \
1594 px1[i] = _NA_; \
1595 if (again) \
1596 px1[j] = _NA_; \
1597 } else if (narm_) { \
1598 --count[i]; \
1599 if (again) \
1600 --count[j]; \
1601 } \
1602 } else { \
1603 tmp = _CAST_(*px0); \
1604 _INCREMENT_(px1[i], tmp); \
1605 if (again) \
1606 _INCREMENT_(px1[j], tmp); \
1607 } \
1608 ++px0; \
1609 } \
1610 } while (0)
1611
1612 SUM_CASES;
1613
1614#undef SUM_LOOP
1615#undef SUM_KERNEL
1616
1617 if (narm_)
1618 Matrix_Free(count, m);
1619 return;
1620}
1621
1622SEXP dense_marginsum(SEXP obj, const char *class, int margin,
1623 int narm, int mean)
1624{
1625 SEXP dim = GET_SLOT(obj, Matrix_DimSym);
1626 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1],
1627 r = (margin == 0) ? m : n;
1628
1629 SEXP res = PROTECT(allocVector(SUM_TYPEOF(class[0]), r)),
1630 x = PROTECT(GET_SLOT(obj, Matrix_xSym));
1631
1632 SEXP dimnames = (class[1] != 's')
1633 ? GET_SLOT(obj, Matrix_DimNamesSym)
1634 : get_symmetrized_DimNames(obj, -1),
1635 marnames = VECTOR_ELT(dimnames, margin);
1636 if (marnames != R_NilValue) {
1637 PROTECT(marnames);
1638 setAttrib(res, R_NamesSymbol, marnames);
1639 UNPROTECT(1); /* marnames */
1640 }
1641
1642 char ul = 'U', di = 'N';
1643 if (class[1] != 'g') {
1644 SEXP uplo = GET_SLOT(obj, Matrix_uploSym);
1645 ul = *CHAR(STRING_ELT(uplo, 0));
1646 if (class[1] == 't') {
1647 SEXP diag = GET_SLOT(obj, Matrix_diagSym);
1648 di = *CHAR(STRING_ELT(diag, 0));
1649 }
1650 }
1651
1652 if (margin == 0 || class[1] == 's')
1653 dense_rowsum(x, class, m, n, ul, di, narm, mean, res);
1654 else
1655 dense_colsum(x, class, m, n, ul, di, narm, mean, res);
1656
1657 UNPROTECT(2); /* x, res */
1658 return res;
1659}
1660
1661/* (row|col)(Sums|Means)(<denseMatrix>) */
1662SEXP R_dense_marginsum(SEXP obj, SEXP margin,
1663 SEXP narm, SEXP mean)
1664{
1665 static const char *valid[] = { VALID_DENSE, "" };
1666 int ivalid = R_check_class_etc(obj, valid);
1667 if (ivalid < 0)
1668 ERROR_INVALID_CLASS(obj, __func__);
1669
1670 int margin_;
1671 if (TYPEOF(margin) != INTSXP || LENGTH(margin) < 1 ||
1672 ((margin_ = INTEGER(margin)[0]) != 0 && margin_ != 1))
1673 error(_("'%s' must be %d or %d"), "margin", 0, 1);
1674
1675 int narm_;
1676 if (TYPEOF(narm) != LGLSXP || LENGTH(narm) < 1 ||
1677 (narm_ = LOGICAL(narm)[0]) == NA_LOGICAL)
1678 error(_("'%s' must be %s or %s"), "narm", "TRUE", "FALSE");
1679
1680 int mean_;
1681 if (TYPEOF(mean) != LGLSXP || LENGTH(mean) < 1 ||
1682 (mean_ = LOGICAL(mean)[0]) == NA_LOGICAL)
1683 error(_("'%s' must be %s or %s"), "mean", "TRUE", "FALSE");
1684
1685 return dense_marginsum(obj, valid[ivalid], margin_, narm_, mean_);
1686}
1687
1688#undef SUM_CASES
1689#undef SUM_TYPEOF
1690
1691#define TRY_INCREMENT(_LABEL_) \
1692 do { \
1693 if ((s >= 0) \
1694 ? ( t <= MATRIX_INT_FAST64_MAX - s) \
1695 : (-t <= s - MATRIX_INT_FAST64_MIN)) { \
1696 s += t; \
1697 t = 0; \
1698 count = 0; \
1699 } else { \
1700 over = 1; \
1701 goto _LABEL_; \
1702 } \
1703 } while (0)
1704
1705#define LONGDOUBLE_AS_DOUBLE(v) \
1706 (v > DBL_MAX) ? R_PosInf : ((v < -DBL_MAX) ? R_NegInf : (double) v);
1707
1708SEXP dense_sum(SEXP obj, const char *class, int narm)
1709{
1710 SEXP res;
1711
1712 SEXP dim = GET_SLOT(obj, Matrix_DimSym);
1713 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1];
1714
1715 char ul = 'U', di = 'N';
1716 if (class[1] != 'g') {
1717 SEXP uplo = GET_SLOT(obj, Matrix_uploSym);
1718 ul = *CHAR(STRING_ELT(uplo, 0));
1719 if (class[1] == 't') {
1720 SEXP diag = GET_SLOT(obj, Matrix_diagSym);
1721 di = *CHAR(STRING_ELT(diag, 0));
1722 }
1723 }
1724
1725 SEXP x = GET_SLOT(obj, Matrix_xSym);
1726 int i, j, unpacked = class[2] != 'p', symmetric = class[1] == 's';
1727
1728#define SUM_LOOP \
1729 do { \
1730 if (class[1] == 'g') { \
1731 for (j = 0; j < n; ++j) \
1732 SUM_KERNEL(for (i = 0; i < m; ++i)); \
1733 } else if (class[1] == 's' || di == 'N') { \
1734 if (ul == 'U') { \
1735 for (j = 0; j < n; ++j) { \
1736 SUM_KERNEL(for (i = 0; i <= j; ++i)); \
1737 if (unpacked) \
1738 px += m - j - 1; \
1739 } \
1740 } else { \
1741 for (j = 0; j < n; ++j) { \
1742 if (unpacked) \
1743 px += j; \
1744 SUM_KERNEL(for (i = j; i < m; ++i)); \
1745 } \
1746 } \
1747 } else { \
1748 if (ul == 'U') { \
1749 for (j = 0; j < n; ++j) { \
1750 SUM_KERNEL(for (i = 0; i < j; ++i)); \
1751 ++px; \
1752 if (unpacked) \
1753 px += m - j - 1; \
1754 } \
1755 } else { \
1756 for (j = 0; j < n; ++j) { \
1757 if (unpacked) \
1758 px += j; \
1759 ++px; \
1760 SUM_KERNEL(for (i = j + 1; i < m; ++i)); \
1761 } \
1762 } \
1763 } \
1764 } while (0)
1765
1766 if (class[0] == 'n') {
1767 int *px = LOGICAL(x);
1768 Matrix_int_fast64_t s = (di == 'N') ? 0LL : n;
1769
1770#define SUM_KERNEL(_FOR_) \
1771 do { \
1772 _FOR_ { \
1773 if (*px != 0) \
1774 s += (symmetric && i != j) ? 2 : 1; \
1775 ++px; \
1776 } \
1777 } while (0)
1778
1779 SUM_LOOP;
1780
1781#undef SUM_KERNEL
1782
1783 if (s <= INT_MAX) {
1784 res = allocVector(INTSXP, 1);
1785 INTEGER(res)[0] = (int) s;
1786 } else {
1787 res = allocVector(REALSXP, 1);
1788 REAL(res)[0] = (double) s;
1789 }
1790 return res;
1791 }
1792
1793 if (!narm && (class[0] == 'l' || class[0] == 'i')) {
1794 int *px = (class[0] == 'l') ? LOGICAL(x) : INTEGER(x);
1795
1796#define SUM_KERNEL(_FOR_) \
1797 do { \
1798 _FOR_ { \
1799 if (*px == NA_INTEGER) { \
1800 res = allocVector(INTSXP, 1); \
1801 INTEGER(res)[0] = NA_INTEGER; \
1802 return res; \
1803 } \
1804 ++px; \
1805 } \
1806 } while (0)
1807
1808 SUM_LOOP;
1809
1810#undef SUM_KERNEL
1811
1812 }
1813
1814 if (class[0] == 'z') {
1815 Rcomplex *px = COMPLEX(x);
1816 long double zr = (di == 'N') ? 0.0L : n, zi = 0.0L;
1817
1818#define SUM_KERNEL(_FOR_) \
1819 do { \
1820 _FOR_ { \
1821 if (!(narm && (ISNAN((*px).r) || ISNAN((*px).i)))) { \
1822 zr += (symmetric && i != j) \
1823 ? 2.0L * (*px).r : (*px).r; \
1824 zi += (symmetric && i != j) \
1825 ? 2.0L * (*px).i : (*px).i; \
1826 } \
1827 ++px; \
1828 } \
1829 } while (0)
1830
1831 SUM_LOOP;
1832
1833#undef SUM_KERNEL
1834
1835 res = allocVector(CPLXSXP, 1);
1836 COMPLEX(res)[0].r = LONGDOUBLE_AS_DOUBLE(zr);
1837 COMPLEX(res)[0].i = LONGDOUBLE_AS_DOUBLE(zi);
1838 } else if (class[0] == 'd') {
1839 double *px = REAL(x);
1840 long double zr = (di == 'N') ? 0.0L : n;
1841
1842#define SUM_KERNEL(_FOR_) \
1843 do { \
1844 _FOR_ { \
1845 if (!(narm && ISNAN(*px))) \
1846 zr += (symmetric && i != j) \
1847 ? 2.0L * *px : *px; \
1848 ++px; \
1849 } \
1850 } while (0)
1851
1852 SUM_LOOP;
1853
1854#undef SUM_KERNEL
1855
1856 res = allocVector(REALSXP, 1);
1857 REAL(res)[0] = LONGDOUBLE_AS_DOUBLE(zr);
1858 } else {
1859 int *px = (class[0] == 'i') ? INTEGER(x) : LOGICAL(x);
1860 Matrix_int_fast64_t s = (di == 'N') ? 0LL : n, t = 0LL;
1861 unsigned int count = 0;
1862 int over = 0;
1863
1864#define SUM_KERNEL(_FOR_) \
1865 do { \
1866 _FOR_ { \
1867 if (!(narm && *px == NA_INTEGER)) { \
1868 int d = (symmetric && i != j) ? 2 : 1; \
1869 if (count > UINT_MAX - d) \
1870 TRY_INCREMENT(ifover); \
1871 t += (d == 2) ? 2LL * *px : *px; \
1872 count += d; \
1873 } \
1874 ++px; \
1875 } \
1876 } while (0)
1877
1878 SUM_LOOP;
1879
1880#undef SUM_KERNEL
1881
1882 TRY_INCREMENT(ifover);
1883 ifover:
1884 if (over) {
1885 long double zr = (di == 'N') ? 0.0L : n; /* FIXME: wasteful */
1886 px = (class[0] == 'i') ? INTEGER(x) : LOGICAL(x);
1887
1888#define SUM_KERNEL(_FOR_) \
1889 do { \
1890 _FOR_ { \
1891 if (!(narm && *px == NA_INTEGER)) \
1892 zr += (symmetric && i != j) \
1893 ? 2.0L * *px : *px; \
1894 ++px; \
1895 } \
1896 } while (0)
1897
1898 SUM_LOOP;
1899
1900#undef SUM_KERNEL
1901
1902 res = allocVector(REALSXP, 1);
1903 REAL(res)[0] = LONGDOUBLE_AS_DOUBLE(zr);
1904 } else if (s > INT_MIN && s <= INT_MAX) {
1905 res = allocVector(INTSXP, 1);
1906 INTEGER(res)[0] = (int) s;
1907 } else {
1908 res = allocVector(REALSXP, 1);
1909 REAL(res)[0] = (double) s;
1910 }
1911 }
1912
1913#undef SUM_LOOP
1914
1915 return res;
1916}
1917
1918/* sum(<denseMatrix>) */
1919SEXP R_dense_sum(SEXP obj, SEXP narm)
1920{
1921 static const char *valid[] = { VALID_DENSE, "" };
1922 int ivalid = R_check_class_etc(obj, valid);
1923 if (ivalid < 0)
1924 ERROR_INVALID_CLASS(obj, __func__);
1925
1926 int narm_;
1927 if (TYPEOF(narm) != LGLSXP || LENGTH(narm) < 1 ||
1928 (narm_ = LOGICAL(narm)[0]) == NA_LOGICAL)
1929 error(_("'%s' must be %s or %s"), "narm", "TRUE", "FALSE");
1930
1931 return dense_sum(obj, valid[ivalid], narm_);
1932}
1933
1934SEXP dense_prod(SEXP obj, const char *class, int narm)
1935{
1936 SEXP res = PROTECT(allocVector((class[0] == 'z') ? CPLXSXP : REALSXP, 1));
1937
1938 SEXP dim = GET_SLOT(obj, Matrix_DimSym);
1939 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1];
1940
1941 char ul = 'U', di = 'N';
1942 if (class[1] != 'g') {
1943 SEXP uplo = GET_SLOT(obj, Matrix_uploSym);
1944 ul = *CHAR(STRING_ELT(uplo, 0));
1945 if (class[1] == 't') {
1946 SEXP diag = GET_SLOT(obj, Matrix_diagSym);
1947 di = *CHAR(STRING_ELT(diag, 0));
1948 }
1949 }
1950
1951 SEXP x = GET_SLOT(obj, Matrix_xSym);
1952 int i, j, unpacked = class[2] != 'p', symmetric = class[1] == 's';
1953 long double zr = 1.0L, zi = 0.0L;
1954
1955#define PROD_LOOP \
1956 do { \
1957 if (class[1] == 'g') { \
1958 for (j = 0; j < n; ++j) \
1959 PROD_KERNEL(for (i = 0; i < m; ++i)); \
1960 } else if (class[1] == 's') { \
1961 if (ul == 'U') { \
1962 for (j = 0; j < n; ++j) { \
1963 PROD_KERNEL(for (i = 0; i <= j; ++i)); \
1964 if (unpacked) \
1965 px += m - j - 1; \
1966 } \
1967 } else { \
1968 for (j = 0; j < n; ++j) { \
1969 if (unpacked) \
1970 px += j; \
1971 PROD_KERNEL(for (i = j; i < m; ++i)); \
1972 } \
1973 } \
1974 } else if (di == 'N') { \
1975 if (ul == 'U') { \
1976 for (j = 0; j < n; ++j) { \
1977 if (j == 1) { zr *= 0.0L; zi *= 0.0L; } \
1978 PROD_KERNEL(for (i = 0; i <= j; ++i)); \
1979 if (unpacked) \
1980 px += m - j - 1; \
1981 } \
1982 } else { \
1983 for (j = 0; j < n; ++j) { \
1984 if (j == 1) { zr *= 0.0L; zi *= 0.0L; } \
1985 if (unpacked) \
1986 px += j; \
1987 PROD_KERNEL(for (i = j; i < m; ++i)); \
1988 } \
1989 } \
1990 } else { \
1991 if (ul == 'U') { \
1992 for (j = 0; j < n; ++j) { \
1993 if (j == 1) { zr *= 0.0L; zi *= 0.0L; } \
1994 PROD_KERNEL(for (i = 0; i < j; ++i)); \
1995 ++px; \
1996 if (unpacked) \
1997 px += m - j - 1; \
1998 } \
1999 } else { \
2000 for (j = 0; j < n; ++j) { \
2001 if (j == 1) { zr *= 0.0L; zi *= 0.0L; } \
2002 if (unpacked) \
2003 px += j; \
2004 ++px; \
2005 PROD_KERNEL(for (i = j + 1; i < m; ++i)); \
2006 } \
2007 } \
2008 } \
2009 } while (0)
2010
2011 if (class[0] == 'n') {
2012 int *px = LOGICAL(x);
2013 if (class[1] == 't')
2014 REAL(res)[0] = (n > 1 || (n == 1 && *px == 0)) ? 0.0 : 1.0;
2015 else {
2016
2017#define PROD_KERNEL(_FOR_) \
2018 do { \
2019 _FOR_ { \
2020 if (*px == 0) { \
2021 REAL(res)[0] = 0.0; \
2022 UNPROTECT(1); /* res */ \
2023 return res; \
2024 } \
2025 ++px; \
2026 } \
2027 } while (0)
2028
2029 PROD_LOOP;
2030
2031#undef PROD_KERNEL
2032
2033 REAL(res)[0] = 1.0;
2034 }
2035 UNPROTECT(1); /* res */
2036 return res;
2037 }
2038
2039 if (class[0] == 'z') {
2040 Rcomplex *px = COMPLEX(x);
2041 long double zr0, zi0;
2042
2043#define PROD_KERNEL(_FOR_) \
2044 do { \
2045 _FOR_ { \
2046 if (!(narm && (ISNAN((*px).r) || ISNAN((*px).i)))) { \
2047 zr0 = zr; zi0 = zi; \
2048 zr = zr0 * (*px).r - zi0 * (*px).i; \
2049 zi = zr0 * (*px).i + zi0 * (*px).r; \
2050 if (symmetric && i != j) { \
2051 zr0 = zr; zi0 = zi; \
2052 zr = zr0 * (*px).r - zi0 * (*px).i; \
2053 zi = zr0 * (*px).i + zi0 * (*px).r; \
2054 } \
2055 } \
2056 ++px; \
2057 } \
2058 } while (0)
2059
2060 PROD_LOOP;
2061
2062#undef PROD_KERNEL
2063
2064 } else if (class[0] == 'd') {
2065 double *px = REAL(x);
2066
2067#define PROD_KERNEL(_FOR_) \
2068 do { \
2069 _FOR_ { \
2070 if (!(narm && ISNAN(*px))) \
2071 zr *= (symmetric && i != j) \
2072 ? (long double) *px * *px : *px; \
2073 ++px; \
2074 } \
2075 } while (0)
2076
2077 PROD_LOOP;
2078
2079#undef PROD_KERNEL
2080
2081 } else {
2082 int *px = (class[0] == 'l') ? LOGICAL(x) : INTEGER(x);
2083
2084#define PROD_KERNEL(_FOR_) \
2085 do { \
2086 _FOR_ { \
2087 if (*px != NA_INTEGER) \
2088 zr *= (symmetric && i != j) \
2089 ? (long double) *px * *px : *px; \
2090 else if (!narm) \
2091 zr *= NA_REAL; \
2092 ++px; \
2093 } \
2094 } while (0)
2095
2096 PROD_LOOP;
2097
2098#undef PROD_KERNEL
2099
2100 }
2101
2102#undef PROD_LOOP
2103
2104 if (class[0] == 'z') {
2105 COMPLEX(res)[0].r = LONGDOUBLE_AS_DOUBLE(zr);
2106 COMPLEX(res)[0].i = LONGDOUBLE_AS_DOUBLE(zi);
2107 } else
2108 REAL(res)[0] = LONGDOUBLE_AS_DOUBLE(zr);
2109 UNPROTECT(1); /* res */
2110 return res;
2111}
2112
2113/* prod(<denseMatrix>) */
2114SEXP R_dense_prod(SEXP obj, SEXP narm)
2115{
2116 static const char *valid[] = { VALID_DENSE, "" };
2117 int ivalid = R_check_class_etc(obj, valid);
2118 if (ivalid < 0)
2119 ERROR_INVALID_CLASS(obj, __func__);
2120
2121 int narm_;
2122 if (TYPEOF(narm) != LGLSXP || LENGTH(narm) < 1 ||
2123 (narm_ = LOGICAL(narm)[0]) == NA_LOGICAL)
2124 error(_("'%s' must be %s or %s"), "narm", "TRUE", "FALSE");
2125
2126 return dense_prod(obj, valid[ivalid], narm_);
2127}
2128
2129#undef TRY_INCREMENT
2130#undef LONGDOUBLE_AS_DOUBLE
2131
2132/* MJ: unused */
2133#if 0
2134
2151static
2152int left_cyclic(double *x, int ldx, int j, int k,
2153 double *cosines, double *sines)
2154{
2155 if (j < 0)
2156 error(_("incorrect left cyclic shift, j (%d) < 0"),
2157 j);
2158 if (j >= k)
2159 error(_("incorrect left cyclic shift, j (%d) >= k (%d)"),
2160 j, k);
2161 if (ldx < k)
2162 error(_("incorrect left cyclic shift, k (%d) > ldx (%d)"),
2163 k, ldx);
2164
2165 double *lastcol = (double *) R_alloc((size_t) k + 1, sizeof(double));
2166 int i;
2167 /* keep a copy of column j */
2168 for (i = 0; i <= j; i++)
2169 lastcol[i] = x[i + j*ldx];
2170 /* For safety, zero the rest */
2171 for (i = j+1; i <= k; i++)
2172 lastcol[i] = 0.;
2173 for (int jj = j+1, ind = 0; jj <= k; jj++, ind++) {
2174 /* columns to be shifted */
2175 int diagind = jj*(ldx+1); // ind == (jj-j) - 1
2176 double tmp = x[diagind], cc, ss;
2177 /* Calculate the Givens rotation. */
2178 /* This modified the super-diagonal element */
2179 F77_CALL(drotg)(x+diagind-1, &tmp, cosines+ind, sines+ind);
2180 cc = cosines[ind];
2181 ss = sines[ind];
2182 /* Copy column jj+1 to column jj. */
2183 for (i = 0; i < jj; i++)
2184 x[i + (jj-1)*ldx] = x[i+jj*ldx];
2185 /* Apply rotation to columns up to k */
2186 for (i = jj; i < k; i++) {
2187 tmp = cc*x[(jj-1)+i*ldx] + ss*x[jj+i*ldx];
2188 x[jj+i*ldx] = cc*x[jj+i*ldx] - ss*x[(jj-1)+i*ldx];
2189 x[(jj-1)+i*ldx] = tmp;
2190 }
2191 /* Apply rotation to lastcol */
2192 lastcol[jj] = -ss*lastcol[jj-1]; lastcol[jj-1] *= cc;
2193 }
2194 /* Copy lastcol to column k */
2195 for(i = 0; i <= k; i++)
2196 x[i+k*ldx] = lastcol[i];
2197 return 0;
2198}
2199
2200static
2201SEXP getGivens(double *x, int ldx, int jmin, int rank)
2202{
2203 int shiftlen = (rank - jmin) - 1;
2204 SEXP ans = PROTECT(allocVector(VECSXP, 4)), nms, cosines, sines;
2205 SET_VECTOR_ELT(ans, 0, ScalarInteger(jmin));
2206 SET_VECTOR_ELT(ans, 1, ScalarInteger(rank));
2207 SET_VECTOR_ELT(ans, 2, cosines = allocVector(REALSXP, shiftlen));
2208 SET_VECTOR_ELT(ans, 3, sines = allocVector(REALSXP, shiftlen));
2209 setAttrib(ans, R_NamesSymbol, nms = allocVector(STRSXP, 4));
2210 SET_STRING_ELT(nms, 0, mkChar("jmin"));
2211 SET_STRING_ELT(nms, 1, mkChar("rank"));
2212 SET_STRING_ELT(nms, 2, mkChar("cosines"));
2213 SET_STRING_ELT(nms, 3, mkChar("sines"));
2214 if (left_cyclic(x, ldx, jmin, rank - 1, REAL(cosines), REAL(sines)))
2215 error(_("unknown error in getGivens"));
2216 UNPROTECT(1);
2217 return ans;
2218}
2219
2220static
2221SEXP checkGivens(SEXP X, SEXP jmin, SEXP rank)
2222{
2223 if (!(isReal(X) && isMatrix(X)))
2224 error(_("X must be a numeric (double precision) matrix"));
2225 SEXP ans = PROTECT(allocVector(VECSXP, 2)),
2226 Xcp = PROTECT(duplicate(X));
2227 int Xdims = INTEGER(getAttrib(X, R_DimSymbol));
2228 SET_VECTOR_ELT(ans, 0, Xcp);
2229 SET_VECTOR_ELT(ans, 1, getGivens(REAL(Xcp), Xdims[0],
2230 asInteger(jmin), asInteger(rank)));
2231 UNPROTECT(2);
2232 return ans;
2233}
2234
2235SEXP lsq_dense_chol(SEXP X, SEXP y)
2236{
2237 if (!(isReal(X) && isMatrix(X)))
2238 error(_("X must be a numeric (double precision) matrix"));
2239 if (!(isReal(y) && isMatrix(y)))
2240 error(_("y must be a numeric (double precision) matrix"));
2241 int *Xdims = INTEGER(getAttrib(X, R_DimSymbol)),
2242 *ydims = INTEGER(getAttrib(y, R_DimSymbol));
2243 if (Xdims[0] != ydim[0])
2244 error(_("number of rows in y (%d) does not match "
2245 "number of rows in X (%d)"),
2246 ydims[0], Xdims[0]);
2247 int n = Xdims[0], p = Xdims[1], k = ydims[1];
2248 if (p < 1 || k < 1)
2249 return allocMatrix(REALSXP, p, k);
2250 SEXP ans = PROTECT(allocMatrix(REALSXP, p, k));
2251 double d_one = 1.0, d_zero = 0.0,
2252 *xpx = (double *) R_alloc((size_t) p * p, sizeof(double));
2253 int info;
2254 F77_CALL(dgemm)("T", "N", &p, &k, &n, &d_one, REAL(X), &n, REAL(y),
2255 &n, &d_zero, REAL(ans), &p FCONE FCONE);
2256 F77_CALL(dsyrk)("U", "T", &p, &n, &d_one, REAL(X), &n, &d_zero,
2257 xpx, &p FCONE FCONE);
2258 F77_CALL(dposv)("U", &p, &k, xpx, &p, REAL(ans), &p, &info FCONE);
2259 if (info)
2260 error(_("LAPACK dposv returned error code %d"), info);
2261 UNPROTECT(1);
2262 return ans;
2263}
2264
2265SEXP lsq_dense_qr(SEXP X, SEXP y)
2266{
2267 if (!(isReal(X) && isMatrix(X)))
2268 error(_("X must be a numeric (double precision) matrix"));
2269 if (!(isReal(y) && isMatrix(y)))
2270 error(_("y must be a numeric (double precision) matrix"));
2271 int *Xdims = INTEGER(getAttrib(X, R_DimSymbol)),
2272 *ydims = INTEGER(getAttrib(y, R_DimSymbol));
2273 if (Xdims[0] != ydim[0])
2274 error(_("number of rows in y (%d) does not match "
2275 "number of rows in X (%d)"),
2276 ydims[0], Xdims[0]);
2277 int n = Xdims[0], p = Xdims[1], k = ydims[1];
2278 if (p < 1 || k < 1)
2279 return allocMatrix(REALSXP, p, k);
2280 SEXP ans = PROTECT(duplicate(y));
2281 double *xvals = (double *) R_alloc((size_t) n * p, sizeof(double)),
2282 *work, tmp;
2283 int lwork = -1, info;
2284 Memcpy(xvals, REAL(X), (size_t) n * p);
2285 F77_CALL(dgels)("N", &n, &p, &k, xvals, &n, REAL(ans), &n,
2286 &tmp, &lwork, &info FCONE);
2287 if (info)
2288 error(_("LAPACK dgels returned error code %d"), info);
2289 lwork = (int) tmp;
2290 work = (double *) R_alloc((size_t) lwork, sizeof(double));
2291 F77_CALL(dgels)("N", &n, &p, &k, xvals, &n, REAL(ans), &n,
2292 work, &lwork, &info FCONE);
2293 if (info)
2294 error(_("LAPACK dgels returned error code %d"), info);
2295 UNPROTECT(1);
2296 return ans;
2297}
2298
2299/* Rank-Correcting/Adapting LAPACK QR Decomposition
2300 * From Doug Bates' initial import; __unused__
2301 *
2302 * Provides a qr() with 'rcond' and rank reduction while(rcond < tol),
2303 * possibly via Givens rotations but WITHOUT PIVOTING
2304 *
2305 * .Call(Matrix:::lapack_qr, A, 1e-17)
2306 * --> ~/R/MM/Pkg-ex/Matrix/qr-rank-deficient.R
2307 *
2308 * TODO: export as Matrix::qrNoPiv() or qr1() or similar
2309 */
2310SEXP lapack_qr(SEXP Xin, SEXP tl)
2311{
2312 if (!(isReal(Xin) && isMatrix(Xin)))
2313 error(_("X must be a real (numeric) matrix"));
2314 double tol = asReal(tl);
2315 if (tol < 0.0)
2316 error(_("tol, given as %g, must be >= 0"), tol);
2317 if (tol > 1.0)
2318 error(_("tol, given as %g, must be <= 1"), tol);
2319 SEXP ans = PROTECT(allocVector(VECSXP, 5)), X, qraux, pivot;
2320 int *Xdims = INTEGER(getAttrib(Xin, R_DimSymbol)),
2321 n = Xdims[0],
2322 p = Xdims[1],
2323 /* size of triangular part of decomposition : */
2324 trsz = (n < p) ? n : p,
2325 i;
2326 SET_VECTOR_ELT(ans, 0, X = duplicate(Xin));
2327 SET_VECTOR_ELT(ans, 2, qraux = allocVector(REALSXP, trsz));
2328 SET_VECTOR_ELT(ans, 3, pivot = allocVector(INTSXP, p));
2329 for (i = 0; i < p; i++)
2330 INTEGER(pivot)[i] = i + 1;
2331 SEXP nms, Givens = PROTECT(allocVector(VECSXP, trsz - 1));
2332 setAttrib(ans, R_NamesSymbol, nms = allocVector(STRSXP, 5));
2333 SET_STRING_ELT(nms, 0, mkChar("qr"));
2334 SET_STRING_ELT(nms, 1, mkChar("rank"));
2335 SET_STRING_ELT(nms, 2, mkChar("qraux"));
2336 SET_STRING_ELT(nms, 3, mkChar("pivot"));
2337 SET_STRING_ELT(nms, 4, mkChar("Givens"));
2338 int rank = trsz, nGivens = 0;
2339 double rcond = 0.0;
2340 if (n > 0 && p > 0) {
2341 double *xpt = REAL(X), *work, tmp;
2342 int *iwork, lwork, info;
2343 lwork = -1;
2344 F77_CALL(dgeqrf)(&n, &p, xpt, &n, REAL(qraux), &tmp, &lwork,
2345 &info);
2346 if (info)
2347 error(_("LAPACK dgeqrf returned error code %d"), info);
2348 lwork = (int) tmp;
2349 work = (double *) R_alloc(((size_t) lwork < (size_t) 3 * trsz)
2350 ? (size_t) 3 * trsz : (size_t) lwork,
2351 sizeof(double));
2352 F77_CALL(dgeqrf)(&n, &p, xpt, &n, REAL(qraux), work, &lwork,
2353 &info);
2354 if (info)
2355 error(_("LAPACK dgeqrf returned error code %d"), info);
2356 iwork = (int *) R_alloc((size_t) trsz, sizeof(int));
2357 F77_CALL(dtrcon)("1", "U", "N", &rank, xpt, &n, &rcond,
2358 work, iwork, &info FCONE FCONE FCONE);
2359 if (info)
2360 error(_("LAPACK dtrcon returned error code %d"), info);
2361 while (rcond < tol) { /* check diagonal elements */
2362 double el, minabs = (xpt[0] < 0.0) ? -xpt[0]: xpt[0];
2363 int jmin = 0;
2364 for (i = 1; i < rank; i++) {
2365 el = xpt[i*n]; // had i*(n+1) which looks wrong to MM
2366 if (el < 0.0)
2367 el = -el;
2368 if (el < minabs) {
2369 jmin = i;
2370 minabs = el;
2371 }
2372 }
2373 if (jmin < (rank - 1)) {
2374 SET_VECTOR_ELT(Givens, nGivens,
2375 getGivens(xpt, n, jmin, rank));
2376 nGivens++;
2377 } // otherwise jmin == (rank - 1), so just "drop that column"
2378 rank--;
2379 // new rcond := ... for reduced rank
2380 F77_CALL(dtrcon)("1", "U", "N", &rank, xpt, &n, &rcond,
2381 work, iwork, &info FCONE FCONE FCONE);
2382 if (info)
2383 error(_("LAPACK dtrcon returned error code %d"), info);
2384 }
2385 }
2386 SEXP Gcpy;
2387 SET_VECTOR_ELT(ans, 1, ScalarInteger(rank));
2388 SET_VECTOR_ELT(ans, 4, Gcpy = allocVector(VECSXP, nGivens));
2389 for (i = 0; i < nGivens; i++)
2390 SET_VECTOR_ELT(Gcpy, i, VECTOR_ELT(Givens, i));
2391 setAttrib(ans, install("useLAPACK"), ScalarLogical(1));
2392 setAttrib(ans, install("rcond"), ScalarReal(rcond));
2393 UNPROTECT(2);
2394 return ans;
2395}
2396
2397#endif /* MJ */
FCONE
#define FCONE
Definition Lapack-etc.h:18
SCALE1_REAL
#define SCALE1_REAL(_X_, _A_)
Definition Mdefines.h:184
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
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
NOTCONJ_PATTERN
#define NOTCONJ_PATTERN(_X_, _Y_)
Definition Mdefines.h:131
VALID_DENSE
#define VALID_DENSE
Definition Mdefines.h:250
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
ISNZ_REAL
#define ISNZ_REAL(_X_)
Definition Mdefines.h:111
ASSIGN_COMPLEX
#define ASSIGN_COMPLEX(_X_, _Y_)
Definition Mdefines.h:181
ISNZ_PATTERN
#define ISNZ_PATTERN(_X_)
Definition Mdefines.h:108
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
ISNZ_LOGICAL
#define ISNZ_LOGICAL(_X_)
Definition Mdefines.h:109
Matrix_Free
#define Matrix_Free(_VAR_, _N_)
Definition Mdefines.h:77
ISNZ_INTEGER
#define ISNZ_INTEGER(_X_)
Definition Mdefines.h:110
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
ISNZ_COMPLEX
#define ISNZ_COMPLEX(_X_)
Definition Mdefines.h:112
SCALE1_COMPLEX
#define SCALE1_COMPLEX(_X_, _A_)
Definition Mdefines.h:186
newObject
SEXP newObject(const char *)
Definition objects.c:4
NOTCONJ_INTEGER
#define NOTCONJ_INTEGER(_X_, _Y_)
Definition Mdefines.h:135
Matrix_sdSym
SEXP Matrix_sdSym
Definition Msymbols.h:20
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_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
dense_as_kind
SEXP dense_as_kind(SEXP from, const char *class, char kind, int new)
Definition coerce.c:2409
matrix_as_dense
SEXP matrix_as_dense(SEXP from, const char *zzz, char ul, char di, int trans, int new)
Definition coerce.c:297
dense_as_general
SEXP dense_as_general(SEXP from, const char *class, int new)
Definition coerce.c:2734
ID_LOOP
#define ID_LOOP(_CTYPE_, _PTR_, _ISNZ_)
R_dense_is_diagonal
SEXP R_dense_is_diagonal(SEXP obj)
Definition dense.c:1355
dense_force_symmetric
SEXP dense_force_symmetric(SEXP from, const char *class, char ul)
Definition dense.c:580
TRANS_LOOP
#define TRANS_LOOP(_CTYPE_, _PTR_)
R_dense_skewpart
SEXP R_dense_skewpart(SEXP from)
Definition dense.c:1011
SUM_CASES
#define SUM_CASES
Definition dense.c:1379
dense_is_symmetric
int dense_is_symmetric(SEXP obj, const char *class, int checkDN)
Definition dense.c:1021
SP_LOOP
#define SP_LOOP(_CTYPE_, _PTR_, _INCREMENT_, _SCALE1_, _ONE_)
R_dense_transpose
SEXP R_dense_transpose(SEXP from)
Definition dense.c:568
R_dense_diag_get
SEXP R_dense_diag_get(SEXP obj, SEXP names)
Definition dense.c:312
SUM_TYPEOF
#define SUM_TYPEOF(c)
Definition dense.c:1422
dense_band
SEXP dense_band(SEXP from, const char *class, int a, int b)
Definition dense.c:5
dense_diag_set
SEXP dense_diag_set(SEXP from, const char *class, SEXP value, int new)
Definition dense.c:327
IS_LOOP
#define IS_LOOP(_CTYPE_, _PTR_, _NOTREAL_, _NOTCONJ_)
R_dense_is_triangular
SEXP R_dense_is_triangular(SEXP obj, SEXP upper)
Definition dense.c:1218
DS_LOOP
#define DS_LOOP(_CTYPE_, _PTR_)
R_dense_is_symmetric
SEXP R_dense_is_symmetric(SEXP obj, SEXP checkDN)
Definition dense.c:1086
DCPY1
#define DCPY1(_PREFIX_, _CTYPE_, _PTR_)
R_dense_symmpart
SEXP R_dense_symmpart(SEXP from)
Definition dense.c:844
dense_symmpart
SEXP dense_symmpart(SEXP from, const char *class)
Definition dense.c:702
R_dense_band
SEXP R_dense_band(SEXP from, SEXP k1, SEXP k2)
Definition dense.c:189
dense_sum
SEXP dense_sum(SEXP obj, const char *class, int narm)
Definition dense.c:1708
dense_diag_get
SEXP dense_diag_get(SEXP obj, const char *class, int names)
Definition dense.c:226
BAND2
#define BAND2(_PREFIX_, _CTYPE_, _PTR_)
dense_colsum
static void dense_colsum(SEXP x, const char *class, int m, int n, char ul, char di, int narm, int mean, SEXP res)
Definition dense.c:1425
IT_LOOP
#define IT_LOOP(_CTYPE_, _PTR_, _ISNZ_)
BAND1
#define BAND1(_PREFIX_, _CTYPE_, _PTR_)
dense_is_triangular
int dense_is_triangular(SEXP obj, const char *class, int upper)
Definition dense.c:1109
dense_skewpart
SEXP dense_skewpart(SEXP from, const char *class)
Definition dense.c:854
dense_rowsum
static void dense_rowsum(SEXP x, const char *class, int m, int n, char ul, char di, int narm, int mean, SEXP res)
Definition dense.c:1519
BAND_CASES
#define BAND_CASES(_DO_)
R_dense_marginsum
SEXP R_dense_marginsum(SEXP obj, SEXP margin, SEXP narm, SEXP mean)
Definition dense.c:1662
LONGDOUBLE_AS_DOUBLE
#define LONGDOUBLE_AS_DOUBLE(v)
Definition dense.c:1705
DCPY2
#define DCPY2(_PREFIX_, _CTYPE_, _PTR_)
PROD_LOOP
#define PROD_LOOP
SUM_LOOP
#define SUM_LOOP(_CTYPE0_, _PTR0_, _CTYPE1_, _PTR1_, _ZERO_, _ONE_, _NA_, _ISNA_, _CAST_, _INCREMENT_, _SCALE2_)
dense_transpose
SEXP dense_transpose(SEXP from, const char *class)
Definition dense.c:466
R_dense_force_symmetric
SEXP R_dense_force_symmetric(SEXP from, SEXP uplo)
Definition dense.c:684
R_dense_diag_set
SEXP R_dense_diag_set(SEXP from, SEXP value)
Definition dense.c:411
dense_marginsum
SEXP dense_marginsum(SEXP obj, const char *class, int margin, int narm, int mean)
Definition dense.c:1622
R_dense_sum
SEXP R_dense_sum(SEXP obj, SEXP narm)
Definition dense.c:1919
R_dense_prod
SEXP R_dense_prod(SEXP obj, SEXP narm)
Definition dense.c:2114
TRY_INCREMENT
#define TRY_INCREMENT(_LABEL_)
Definition dense.c:1691
dense_is_diagonal
int dense_is_diagonal(SEXP obj, const char *class)
Definition dense.c:1251
dense_prod
SEXP dense_prod(SEXP obj, const char *class, int narm)
Definition dense.c:1934
DG_LOOP
#define DG_LOOP(_CTYPE_, _PTR_, _ONE_)
DCPY
#define DCPY(_PREFIX_, _CTYPE_, _PTR_)
Matrix_zone
Rcomplex Matrix_zone
Definition init.c:26
Matrix_zzero
Rcomplex Matrix_zzero
Definition init.c:26
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
zeroRe
void zeroRe(SEXP x)
Definition utils.c:163
Generated by doxygen 1.9.8