12#define TOLBASED_ISNZ_REAL(_X_) \
13 (ISNA_REAL(_X_) || fabs(_X_) > tol)
15#define TOLBASED_ISNZ_COMPLEX(_X_) \
16 (ISNA_COMPLEX(_X_) || hypot((_X_).r, (_X_).i) > tol)
18#define DROP0_CASES(_DO_) \
22 _DO_(int, LOGICAL, ISNZ_LOGICAL); \
25 _DO_(int, INTEGER, ISNZ_INTEGER); \
29 _DO_(double, REAL, TOLBASED_ISNZ_REAL); \
31 _DO_(double, REAL, ISNZ_REAL); \
35 _DO_(Rcomplex, COMPLEX, TOLBASED_ISNZ_COMPLEX); \
37 _DO_(Rcomplex, COMPLEX, ISNZ_COMPLEX); \
44 if (
class[2] !=
'T') {
47 int *pp0 = INTEGER(p0), k, n = (int) (XLENGTH(p0) - 1),
48 nnz0 = pp0[n], nnz1 = 0;
51#define DROP0_LOOP1(_CTYPE_, _PTR_, _ISNZ_) \
53 _CTYPE_ *px0 = _PTR_(x0); \
54 for (k = 0; k < nnz0; ++k) { \
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),
81#define DROP0_LOOP2(_CTYPE_, _PTR_, _ISNZ_) \
83 _CTYPE_ *px0 = _PTR_(x0), *px1 = _PTR_(x1); \
84 for (j = 0, k = 0; j < n; ++j) { \
85 pp1[j] = pp1[j - 1]; \
103 R_xlen_t k, nnz0 = XLENGTH(x0), nnz1 = 0;
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);
124#define DROP0_LOOP2(_CTYPE_, _PTR_, _ISNZ_) \
126 _CTYPE_ *px0 = _PTR_(x0), *px1 = _PTR_(x1); \
127 for (k = 0; k < nnz0; ++k) { \
128 if (_ISNZ_(*px0)) { \
133 ++pi0; ++pj0; ++px0; \
145 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1];
154 if (
class[1] !=
'g') {
156 char ul = *CHAR(STRING_ELT(uplo, 0));
161 if (
class[1] ==
't') {
163 char di = *CHAR(STRING_ELT(diag, 0));
169 if (LENGTH(factors) > 0)
174#undef TOLBASED_ISNZ_REAL
175#undef TOLBASED_ISNZ_COMPLEX
187 static const char *
valid[] = {
189 int ivalid = R_check_class_etc(from,
valid);
194 if (TYPEOF(tol) != REALSXP || LENGTH(tol) < 1 ||
195 ISNAN(tol_ = REAL(tol)[0]))
196 error(
_(
"'%s' is not a number"),
"tol");
207 char di = *CHAR(STRING_ELT(diag, 0));
212 SEXP val = PROTECT(ScalarLogical(1));
224 static const char *
valid[] = {
226 int ivalid = R_check_class_etc(from,
valid);
239 char di = *CHAR(STRING_ELT(diag, 0));
245 int n = INTEGER(dim)[0];
249 PROTECT(from = duplicate(from));
252 char ul = *CHAR(STRING_ELT(uplo, 0));
255 PROTECT(from =
sparse_band(from,
class, 1, n - 1));
257 PROTECT(from =
sparse_band(from,
class, 1 - n, -1));
260 PROTECT(diag = mkString(
"U"));
272 static const char *
valid[] = {
274 int ivalid = R_check_class_etc(from,
valid);
284 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1];
288 if (a <= 1 - m && b >= n - 1 &&
289 (
class[1] ==
't' || m != n || m > 1 || n > 1))
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'));
296 char ul0 =
'U', ul1 =
'U', di =
'N';
297 if (
class[1] !=
'g') {
299 ul0 = *CHAR(STRING_ELT(uplo, 0));
302 if (
class[1] ==
't') {
304 if ((ul0 ==
'U') ? (a <= 0 && b >= n - 1) : (b >= 0 && a <= 1 - m))
306 else if (a <= 0 && b >= 0) {
308 di = *CHAR(STRING_ELT(diag, 0));
316 if (
class[2] ==
'R') {
319 r = a; a = -b; b = -r;
320 ul0 = (ul0 ==
'U') ?
'L' :
'U';
325 char cl[] =
"...Matrix";
327 cl[1] = (ge) ?
'g' : ((tr) ?
't' :
's');
328 cl[2] = (
class[2] ==
'R') ?
'C' :
class[2];
337 if (
class[1] !=
's' || sy)
344 ul1 = (tr &&
class[1] !=
't') ? ((a >= 0) ?
'U' :
'L') : ul0;
346 SEXP uplo = PROTECT(mkString(
"L"));
351 SEXP diag = PROTECT(mkString(
"U"));
361 switch (class[0]) { \
363 BAND_SUBCASES(int, LOGICAL, SHOW); \
366 BAND_SUBCASES(int, INTEGER, SHOW); \
369 BAND_SUBCASES(double, REAL, SHOW); \
372 BAND_SUBCASES(Rcomplex, COMPLEX, SHOW); \
379 if (
class[2] !=
'T') {
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;
388 if (
class[1] ==
's' && !sy) {
390 for (j = 0, k = 0; j < n; ++j) {
393 if ((d = j - pi0[k]) >= a && d <= b)
395 if (d != 0 && -d >= a && -d <= b)
400 for (j = 0; j < n; ++j) {
405 for (j = 0, k = 0; j < n; ++j) {
408 if ((d = j - pi0[k]) >= a && d <= b)
416 if (nnz1 == nnz0 && (
class[1] !=
's' || sy)) {
419 if (
class[0] !=
'n') {
430 SEXP i1 = PROTECT(allocVector(INTSXP, nnz1));
431 int *pi1 = INTEGER(i1);
435#define BAND_SUBCASES(_CTYPE_, _PTR_, _MASK_) \
437 _MASK_(_CTYPE_ *px0 = _PTR_(x0), *px1 = _PTR_(x1)); \
438 if (class[1] == 's' && !sy) { \
440 Matrix_Calloc(pp1_, n, int); \
441 Matrix_memcpy(pp1_, pp1 - 1, n, sizeof(int)); \
442 for (j = 0, k = 0; j < n; ++j) { \
445 if ((d = j - pi0[k]) >= a && d <= b) { \
446 pi1[pp1_[j]] = pi0[k]; \
447 _MASK_(px1[pp1_[j]] = px0[k]); \
450 if (d != 0 && -d >= a && -d <= b) { \
451 pi1[pp1_[pi0[k]]] = j; \
452 _MASK_(px1[pp1_[pi0[k]]] = px0[k]); \
458 Matrix_Free(pp1_, n); \
460 for (j = 0, k = 0; j < n; ++j) { \
463 if ((d = j - pi0[k]) >= a && d <= b) { \
465 _MASK_(*(px1++) = px0[k]); \
477 x1 = PROTECT(allocVector(TYPEOF(x0), nnz1));
489 int *pi0 = INTEGER(i0), *pj0 = INTEGER(j0), d;
490 R_xlen_t k, nnz0 = XLENGTH(i0), nnz1 = 0;
492 if (
class[1] ==
's' && !sy) {
493 for (k = 0; k < nnz0; ++k) {
494 if ((d = pj0[k] - pi0[k]) >= a && d <= b)
496 if (d != 0 && -d >= a && -d <= b)
500 for (k = 0; k < nnz0; ++k) {
501 if ((d = pj0[k] - pi0[k]) >= a && d <= b)
506 if (nnz1 == nnz0 && (
class[1] !=
's' || sy)) {
510 if (
class[0] !=
'n') {
519 SEXP i1 = PROTECT(allocVector(INTSXP, nnz1)),
520 j1 = PROTECT(allocVector(INTSXP, nnz1));
521 int *pi1 = INTEGER(i1), *pj1 = INTEGER(j1);
526#define BAND_SUBCASES(_CTYPE_, _PTR_, _MASK_) \
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) { \
534 _MASK_(*(px1++) = px0[k]); \
536 if (d != 0 && -d >= a && -d <= b) { \
539 _MASK_(*(px1++) = px0[k]); \
543 for (k = 0; k < nnz0; ++k) { \
544 if ((d = pj0[k] - pi0[k]) >= a && d <= b) { \
547 _MASK_(*(px1++) = px0[k]); \
557 x1 = PROTECT(allocVector(TYPEOF(x0), nnz1));
576 static const char *
valid[] = {
578 int ivalid = R_check_class_etc(from,
valid);
583 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1];
587 if (k1 == R_NilValue)
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)
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);
598 error(
_(
"'%s' (%d) must be less than or equal to '%s' (%d)"),
607 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1], r = (m < n) ? m : n;
610 char ul =
'U', di =
'N';
611 if (
class[1] !=
'g') {
613 ul = *CHAR(STRING_ELT(uplo, 0));
616 if (
class[1] ==
't') {
618 di = *CHAR(STRING_ELT(diag, 0));
623 SEXP res = PROTECT(allocVector(
kindToType(
class[0]), r));
627 switch (class[0]) { \
630 DG_LOOP(int, LOGICAL, SHOW, FIRSTOF, INCREMENT_LOGICAL, 0, 1); \
633 DG_LOOP(int, INTEGER, SHOW, FIRSTOF, INCREMENT_INTEGER, 0, 1); \
636 DG_LOOP(double, REAL, SHOW, FIRSTOF, INCREMENT_REAL, 0.0, 1.0); \
639 DG_LOOP(Rcomplex, COMPLEX, SHOW, FIRSTOF, INCREMENT_COMPLEX, Matrix_zzero, Matrix_zone); \
651#define DG_LOOP(_CTYPE_, _PTR_, _MASK_, _REPLACE_, _INCREMENT_, _ZERO_, _ONE_) \
653 _CTYPE_ *pres = _PTR_(res); \
654 for (j = 0; j < r; ++j) \
660 }
else if (
class[2] !=
'T') {
665 int j, k, kend, *pp0 = INTEGER(p0), *pi0 = INTEGER(i0);
669#define DG_LOOP(_CTYPE_, _PTR_, _MASK_, _REPLACE_, _INCREMENT_, _ZERO_, _ONE_) \
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) { \
681 pres[j] = _REPLACE_(px0[k], 1); \
686 } else if ((class[2] == 'C') == (ul != 'U')) { \
687 for (j = 0, k = 0; j < r; ++j) { \
689 pres[j] = (k < kend && pi0[k] == j) \
690 ? _REPLACE_(px0[k], 1) : _ZERO_; \
694 for (j = 0, k = 0; j < r; ++j) { \
696 pres[j] = (k < kend && pi0[kend - 1] == j) \
697 ? _REPLACE_(px0[kend - 1], 1) : _ZERO_; \
703 if (
class[0] ==
'n') {
717 int *pi0 = INTEGER(i0), *pj0 = INTEGER(j0);
718 R_xlen_t k, nnz0 = XLENGTH(i0);
721#define DG_LOOP(_CTYPE_, _PTR_, _MASK_, _REPLACE_, _INCREMENT_, _ZERO_, _ONE_) \
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) { \
728 _INCREMENT_(pres[*pi0], (*px0)); \
729 ++pi0; ++pj0; _MASK_(++px0); \
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);
760 setAttrib(res, R_NamesSymbol, cn);
761 else if (rn != R_NilValue &&
763 setAttrib(res, R_NamesSymbol, (r == m) ? rn : cn);
775 static const char *
valid[] = {
777 int ivalid = R_check_class_etc(obj,
valid);
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");
792 int v = LENGTH(value) != 1, full = 0;
795 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1], r = (m < n) ? m : n;
804 char ul =
'U', di =
'N';
805 if (
class[1] !=
'g') {
807 ul = *CHAR(STRING_ELT(uplo, 0));
812 if (
class[1] ==
't') {
814 di = *CHAR(STRING_ELT(diag, 0));
821 switch (class[0]) { \
824 DS_LOOP(int, LOGICAL, SHOW, ISNZ_LOGICAL); \
827 DS_LOOP(int, INTEGER, SHOW, ISNZ_INTEGER); \
830 DS_LOOP(double, REAL, SHOW, ISNZ_REAL); \
833 DS_LOOP(Rcomplex, COMPLEX, SHOW, ISNZ_COMPLEX); \
840 if (
class[2] !=
'T') {
842 int n_ = (
class[2] ==
'C') ? n : m;
846 p1 = PROTECT(allocVector(INTSXP, (R_xlen_t) n_ + 1));
848 int *pp0 = INTEGER(p0), *pp1 = INTEGER(p1), *pi0 = INTEGER(i0),
849 j, k, kend, nd0 = 0, nd1 = 0;
852 if (
class[1] ==
'g') {
853 for (j = 0, k = 0; j < r; ++j) {
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)
871 }
else if ((
class[2] ==
'C') == (ul ==
'U')) {
872 for (j = 0, k = 0; j < n_; ++j) {
874 if (k < kend && pi0[kend - 1] == j)
880 for (j = 0, k = 0; j < n_; ++j) {
882 if (k < kend && pi0[k] == j)
890#define DS_LOOP(_CTYPE_, _PTR_, _MASK_, _ISNZ_) \
892 _CTYPE_ *pvalue = _PTR_(value); \
894 for (j = 0; j < r; ++j) { \
895 if (_ISNZ_(pvalue[j])) \
899 for (j = r; j < n_; ++j) \
901 } else if (_ISNZ_(pvalue[0])) { \
903 for (j = 0; j < r; ++j) \
905 for (j = r; j < n_; ++j) \
912 if (nd1 - nd0 > INT_MAX - pp0[n_ - 1])
913 error(
_(
"%s cannot exceed %s"),
"p[length(p)]",
"2^31-1");
915 SEXP i1 = PROTECT(allocVector(INTSXP, pp1[n_ - 1]));
917 int *pi1 = INTEGER(i1);
920#define DS_LOOP(_CTYPE_, _PTR_, _MASK_, _ISNZ_) \
922 _MASK_(_CTYPE_ *px0 = _PTR_(x0), *px1 = _PTR_(x1)); \
923 _CTYPE_ *pvalue = _PTR_(value); \
924 for (j = 0, k = 0; j < r; ++j) { \
926 while (k < kend && pi0[k] < j) { \
927 *(pi1++) = pi0[k] ; \
928 _MASK_(*(px1++) = px0[k]); \
931 if (k < kend && pi0[k] == j) \
933 if ((v) ? _ISNZ_(pvalue[j]) : full) { \
935 _MASK_(*(px1++) = (v) ? pvalue[j] : pvalue[0]); \
938 *(pi1++) = pi0[k] ; \
939 _MASK_(*(px1++) = px0[k]); \
943 for (j = r; j < n_; ++j) { \
946 *(pi1++) = pi0[k] ; \
947 _MASK_(*(px1++) = px0[k]); \
957 x1 = PROTECT(allocVector(TYPEOF(x0), pp1[n_ - 1]));
969 int *pi0 = INTEGER(i0), *pj0 = INTEGER(j0), j, nd0 = 0, nd1 = 0;
970 R_xlen_t k, nnz0 = XLENGTH(i0), nnz1 = nnz0;
973 for (k = 0; k < nnz0; ++k)
974 if (pi0[k] == pj0[k])
978#define DS_LOOP(_CTYPE_, _PTR_, _MASK_, _ISNZ_) \
980 _CTYPE_ *pvalue = _PTR_(value); \
982 for (j = 0; j < r; ++j) \
983 if (_ISNZ_(pvalue[j])) \
985 } else if (_ISNZ_(pvalue[0])) { \
993 if (nd1 - nd0 > R_XLEN_T_MAX - nnz0)
994 error(
_(
"%s cannot exceed %s"),
"length(i)",
"R_XLEN_T_MAX");
997 SEXP i1 = PROTECT(allocVector(INTSXP, nnz1)),
998 j1 = PROTECT(allocVector(INTSXP, nnz1));
1001 int *pi1 = INTEGER(i1), *pj1 = INTEGER(j1);
1004#define DS_LOOP(_CTYPE_, _PTR_, _MASK_, _ISNZ_) \
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]); \
1016 for (j = 0; j < r; ++j) { \
1017 if (_ISNZ_(pvalue[j])) { \
1018 *(pi1++) = *(pj1++) = j; \
1019 _MASK_(*(px1++) = pvalue[j]); \
1022 } else if (full) { \
1023 for (j = 0; j < r; ++j) { \
1024 *(pi1++) = *(pj1++) = j; \
1025 _MASK_(*(px1++) = pvalue[0]); \
1030 if (
class[0] ==
'n')
1034 x1 = PROTECT(allocVector(TYPEOF(x0), nnz1));
1054 static const char *
valid[] = {
1056 int ivalid = R_check_class_etc(from,
valid);
1059 const char *
class =
valid[ivalid];
1061 SEXPTYPE tx =
kindToType(
class[0]), tv = TYPEOF(value);
1070 error(
_(
"replacement diagonal has incompatible type \"%s\""),
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"));
1083 PROTECT(value = coerceVector(value, tx));
1087#ifndef MATRIX_ENABLE_IMATRIX
1090 PROTECT(value = coerceVector(value, REALSXP));
1095#ifndef MATRIX_ENABLE_IMATRIX
1098 class =
valid[R_check_class_etc(from,
valid)];
1109 if (
class[2] ==
'T' || !lazy)
1112 char cl[] =
"...Matrix";
1115 cl[2] = (
class[2] ==
'C') ?
'R' :
'C';
1120 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1];
1124 pdim = INTEGER(dim);
1132 if (
class[1] ==
's')
1138 if (
class[1] !=
'g') {
1140 char ul = *CHAR(STRING_ELT(uplo, 0));
1143 PROTECT(uplo = mkString(
"L"));
1147 if (
class[1] ==
't') {
1149 char di = *CHAR(STRING_ELT(diag, 0));
1155 if (LENGTH(factors) > 0)
1163 if (
class[2] ==
'T') {
1170 if (
class[0] !=
'n') {
1183 i0 = PROTECT(
GET_SLOT(from, iSym));
1191 if (
class[0] !=
'n') {
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_]));
1207 void trans(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP,
int,
int);
1209 if (
class[0] ==
'n')
1210 trans(p0, i0, NULL, p1, i1, NULL, m_, n_);
1213 x1 = PROTECT(allocVector(TYPEOF(x0), INTEGER(p0)[n_]));
1215 trans(p0, i0, x0, p1, i1, x1, m_, n_);
1225 static const char *
valid[] = {
1227 int ivalid = R_check_class_etc(from,
valid);
1232 if (TYPEOF(lazy) != LGLSXP || LENGTH(lazy) < 1 ||
1233 (lazy_ = LOGICAL(lazy)[0]) == NA_LOGICAL)
1234 error(
_(
"invalid '%s' to '%s'"),
"lazy", __func__);
1241 char ul0 =
'U', ul1 =
'U';
1242 if (
class[1] !=
'g') {
1244 ul0 = ul1 = *CHAR(STRING_ELT(uplo, 0));
1250 if (
class[1] ==
's') {
1255 if (
class[0] ==
'z') {
1267 char cl[] =
".s.Matrix";
1273 int *pdim = INTEGER(dim), n = pdim[0];
1275 error(
_(
"attempt to symmetrize a non-square matrix"));
1285 SEXP uplo = PROTECT(mkString(
"L"));
1291 if (
class[1] ==
't') {
1293 di = *CHAR(STRING_ELT(diag, 0));
1301 switch (class[0]) { \
1303 FS_SUBCASES(int, LOGICAL, SHOW, 1); \
1306 FS_SUBCASES(int, INTEGER, SHOW, 1); \
1309 FS_SUBCASES(double, REAL, SHOW, 1.0); \
1312 FS_SUBCASES(Rcomplex, COMPLEX, SHOW, Matrix_zone); \
1319 if (
class[1] ==
't' && di ==
'N' && ul0 == ul1) {
1322 if (
class[2] !=
'T') {
1327 if (
class[2] !=
'R') {
1332 if (
class[2] !=
'C') {
1337 if (
class[0] !=
'n') {
1345 }
else if (
class[2] !=
'T') {
1351 p1 = PROTECT(allocVector(INTSXP, (R_xlen_t) n + 1)),
1352 i0 = PROTECT(
GET_SLOT(from, iSym));
1359 pp0++; *(pp1++) = 0;
1364 if (
class[1] ==
't') {
1369 for (j = 0; j < n; ++j)
1373 for (j = 0; j < n; ++j)
1374 pp1[j] = ++nnz1 + pp0[j];
1377 }
else if (ul0 == ((
class[2] ==
'C') ?
'U' :
'L')) {
1380 for (j = 0; j < n; ++j) {
1381 if (pp0[j - 1] < pp0[j] && pi0[pp0[j] - 1] == j)
1388 for (j = 0; j < n; ++j) {
1389 if (pp0[j - 1] < pp0[j] && pi0[pp0[j - 1]] == j)
1394 }
else if (ul1 == ((
class[2] ==
'C') ?
'U' :
'L')) {
1396 for (j = 0, k = 0; j < n; ++j) {
1407 for (j = 0, k = 0; j < n; ++j) {
1420 SEXP i1 = PROTECT(allocVector(INTSXP, nnz1));
1421 int *pi1 = INTEGER(i1);
1425#define FS_SUBCASES(_CTYPE_, _PTR_, _MASK_, _ONE_) \
1427 _MASK_(_CTYPE_ *px0 = _PTR_(x0), *px1 = _PTR_(x1)); \
1428 if (class[1] == 't') { \
1433 for (j = 0; j < n; ++j) { \
1435 _MASK_(*(px1++) = _ONE_); \
1437 } else if (ul0 == ((class[2] == 'C') ? 'U' : 'L')) { \
1440 for (j = 0, k = 0; j < n; ++j) { \
1442 while (k < kend) { \
1443 *(pi1++) = pi0[k]; \
1444 _MASK_(*(px1++) = px0[k]); \
1448 _MASK_(*(px1++) = _ONE_); \
1453 for (j = 0, k = 0; j < n; ++j) { \
1455 _MASK_(*(px1++) = _ONE_); \
1457 while (k < kend) { \
1458 *(pi1++) = pi0[k]; \
1459 _MASK_(*(px1++) = px0[k]); \
1464 } else if (ul0 == ((class[2] == 'C') ? 'U' : 'L')) { \
1467 for (j = 0; j < n; ++j) { \
1468 if (pp0[j - 1] < pp0[j] && pi0[pp0[j] - 1] == j) { \
1470 _MASK_(*(px1++) = px0[pp0[j] - 1]); \
1476 for (j = 0; j < n; ++j) { \
1477 if (pp0[j - 1] < pp0[j] && pi0[pp0[j - 1]] == j) { \
1479 _MASK_(*(px1++) = px0[pp0[j - 1]]); \
1483 } else if (ul1 == ((class[2] == 'C') ? 'U' : 'L')) { \
1485 for (j = 0, k = 0; j < n; ++j) { \
1487 while (k < kend) { \
1488 if (pi0[k] <= j) { \
1489 *(pi1++) = pi0[k]; \
1490 _MASK_(*(px1++) = px0[k]); \
1497 for (j = 0, k = 0; j < n; ++j) { \
1499 while (k < kend) { \
1500 if (pi0[k] >= j) { \
1501 *(pi1++) = pi0[k]; \
1502 _MASK_(*(px1++) = px0[k]); \
1510 if (
class[0] ==
'n')
1514 x1 = PROTECT(allocVector(TYPEOF(x0), nnz1));
1526 int *pi0 = INTEGER(i0),
1528 R_xlen_t j, k, nnz0 = XLENGTH(i0), nnz1 = 0;
1532 if (
class[1] ==
't' && di !=
'N')
1533 nnz1 = (ul0 == ul1) ? n + nnz0 : n;
1536 for (k = 0; k < nnz0; ++k)
1537 if (pi0[k] <= pj0[k])
1540 for (k = 0; k < nnz0; ++k)
1541 if (pi0[k] >= pj0[k])
1548 SEXP i1 = PROTECT(allocVector(INTSXP, nnz1)),
1549 j1 = PROTECT(allocVector(INTSXP, nnz1));
1550 int *pi1 = INTEGER(i1),
1556#define FS_SUBCASES(_CTYPE_, _PTR_, _MASK_, _ONE_) \
1558 _MASK_(_CTYPE_ *px0 = _PTR_(x0), *px1 = _PTR_(x1)); \
1559 if (class[1] == 't' && di != 'N') { \
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_))); \
1566 _MASK_(px1 += nnz0); \
1568 for (j = 0; j < n; ++j) { \
1569 *(pi1++) = *(pj1++) = j; \
1570 _MASK_(*(px1++) = _ONE_); \
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]); \
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]); \
1593 if (
class[0] ==
'n')
1597 x1 = PROTECT(allocVector(TYPEOF(x0), nnz1));
1615 static const char *
valid[] = {
1617 int ivalid = R_check_class_etc(from,
valid);
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__);
1634 if (
class[0] !=
'z' &&
class[0] !=
'd') {
1639 if (
class[0] !=
'z' &&
class[1] ==
's')
1643 PROTECT_WITH_INDEX(from, &pid);
1645 char cl[] =
".s.Matrix";
1646 cl[0] = (
class[0] !=
'z') ?
'd' :
'z';
1651 int *pdim = INTEGER(dim), n = pdim[0];
1653 error(
_(
"attempt to get symmetric part of non-square matrix"));
1659 if (
class[1] ==
's')
1665 char ul =
'U', di =
'N';
1666 if (
class[1] !=
'g') {
1668 ul = *CHAR(STRING_ELT(uplo, 0));
1672 if (
class[1] ==
't') {
1674 di = *CHAR(STRING_ELT(diag, 0));
1677 }
else if (
class[2] ==
'R') {
1678 SEXP uplo = PROTECT(mkString(
"L"));
1684 if (
class[2] !=
'T') {
1688 i0 = PROTECT(
GET_SLOT(from, iSym)),
1690 int j, k, kend, *pp0 = INTEGER(p0) + 1, *pi0 = INTEGER(i0),
1693 if (
class[1] ==
'g') {
1698 SEXP p1 = PROTECT(allocVector(INTSXP, (R_xlen_t) n + 1)),
1700 i0_ = PROTECT(
GET_SLOT(from, iSym)),
1702 int k_, kend_, *pp0_ = INTEGER(p0_) + 1, *pi0_ = INTEGER(i0_),
1706 for (j = 0, k = 0, k_ = 0; j < n; ++j) {
1709 pp1[j] = pp1[j - 1];
1714 while (k_ < kend_ && pi0_[k_] < pi0[k]) {
1719 if (k_ < kend_ && pi0_[k_] == pi0[k])
1724 while (k_ < kend_) {
1734 SEXP i1 = PROTECT(allocVector(INTSXP, pp1[n - 1])),
1735 x1 = PROTECT(allocVector(
kindToType(
cl[0]), pp1[n - 1]));
1736 int *pi1 = INTEGER(i1);
1739#define SP_LOOP(_CTYPE_, _PTR_, _ASSIGN_, _INCREMENT_) \
1741 _CTYPE_ *px0 = _PTR_(x0), *px1 = _PTR_(x1), \
1742 *px0_ = _PTR_(x0_); \
1743 for (j = 0, k = 0, k_ = 0; j < n; ++j) { \
1746 while (k < kend) { \
1750 while (k_ < kend_ && pi0_[k_] < pi0[k]) { \
1752 _ASSIGN_((*px1), 0.5 * px0_[k_]); \
1753 ++k_; ++pi1; ++px1; \
1756 _ASSIGN_((*px1), 0.5 * px0[k]); \
1757 if (k_ < kend_ && pi0_[k_] == pi0[k]) { \
1758 _INCREMENT_((*px1), 0.5 * px0_[k_]); \
1761 ++k; ++pi1; ++px1; \
1764 while (k_ < kend_) { \
1769 _ASSIGN_((*px1), 0.5 * px0_[k_]); \
1770 ++k_; ++pi1; ++px1; \
1786 }
else if (
class[1] ==
't') {
1788 int leading = (
class[2] ==
'C') == (ul !=
'U');
1792 SEXP x1 = PROTECT(allocVector(
kindToType(
cl[0]), nnz));
1795#define SP_LOOP(_CTYPE_, _PTR_, _ASSIGN_) \
1797 _CTYPE_ *px0 = _PTR_(x0), *px1 = _PTR_(x1); \
1799 for (j = 0, k = 0; j < n; ++j) { \
1805 _ASSIGN_((*px1), 0.5 * (*px0)); \
1806 ++k, ++px0; ++px1; \
1807 while (k < kend) { \
1808 _ASSIGN_((*px1), 0.5 * (*px0)); \
1809 ++k; ++px0; ++px1; \
1814 for (j = 0, k = 0; j < n; ++j) { \
1817 while (k < kend - 1) { \
1818 _ASSIGN_((*px1), 0.5 * (*px0)); \
1819 ++k; ++px0; ++px1; \
1824 _ASSIGN_((*px1), 0.5 * (*px0)); \
1825 ++k; ++px0; ++px1; \
1844 SEXP p1 = PROTECT(allocVector(INTSXP, (R_xlen_t) n + 1)),
1845 i1 = PROTECT(allocVector(INTSXP, nnz)),
1847 int *pp1 = INTEGER(p1), *pi1 = INTEGER(i1);
1851#define SP_LOOP(_CTYPE_, _PTR_, _ASSIGN_, _ONE_) \
1853 _CTYPE_ *px0 = _PTR_(x0), *px1 = _PTR_(x1); \
1855 for (j = 0, k = 0; j < n; ++j) { \
1857 pp1[j] = pp1[j - 1] + kend - k + 1; \
1859 _ASSIGN_((*px1), _ONE_); \
1861 while (k < kend) { \
1863 _ASSIGN_((*px1), 0.5 * (*px0)); \
1864 ++k; ++pi0; ++pi1; ++px0; ++px1; \
1868 for (j = 0, k = 0; j < n; ++j) { \
1870 pp1[j] = pp1[j - 1] + kend - k + 1; \
1871 while (k < kend) { \
1873 _ASSIGN_((*px1), 0.5 * (*px0)); \
1874 ++k; ++pi0; ++pi1; ++px0; ++px1; \
1877 _ASSIGN_((*px1), _ONE_); \
1904 SEXP x1 = PROTECT(duplicate(x0));
1919 int *pi0 = INTEGER(i0), *pj0 = INTEGER(j0);
1920 R_xlen_t k, nnz = XLENGTH(i0);
1922 if (
class[1] ==
'g') {
1924 SEXP i1 = PROTECT(allocVector(INTSXP, nnz)),
1925 j1 = PROTECT(allocVector(INTSXP, nnz)),
1927 int *pi1 = INTEGER(i1), *pj1 = INTEGER(j1);
1930#define SP_LOOP(_CTYPE_, _PTR_, _ASSIGN_) \
1932 _CTYPE_ *px0 = _PTR_(x0), *px1 = _PTR_(x1); \
1933 for (k = 0; k < nnz; ++k) { \
1934 if (*pi0 == *pj0) { \
1938 } else if (*pi0 < *pj0) { \
1941 _ASSIGN_((*px1), 0.5 * (*px0)); \
1945 _ASSIGN_((*px1), 0.5 * (*px0)); \
1947 ++pi0; ++pi1; ++pj0; ++pj1; ++px0; ++px1; \
1961 }
else if (
class[1] ==
't') {
1965 SEXP x1 = PROTECT(allocVector(
kindToType(
cl[0]), nnz));
1968#define SP_LOOP(_CTYPE_, _PTR_, _ASSIGN_) \
1970 _CTYPE_ *px0 = _PTR_(x0), *px1 = _PTR_(x1); \
1971 for (k = 0; k < nnz; ++k) { \
1975 _ASSIGN_((*px1), 0.5 * (*px0)); \
1992 SEXP i1 = PROTECT(allocVector(INTSXP, nnz + n)),
1993 j1 = PROTECT(allocVector(INTSXP, nnz + n)),
1995 int j, *pi1 = INTEGER(i1), *pj1 = INTEGER(j1);
1998#define SP_LOOP(_CTYPE_, _PTR_, _ASSIGN_, _ONE_) \
2000 _CTYPE_ *px0 = _PTR_(x0), *px1 = _PTR_(x1); \
2001 for (k = 0; k < nnz; ++k) { \
2004 _ASSIGN_((*px1), 0.5 * (*px0)); \
2005 ++pi0; ++pi1; ++pj0; ++pj1; ++px0; ++px1; \
2007 for (j = 0; j < n; ++j) { \
2009 _ASSIGN_((*px1), _ONE_); \
2010 ++pi1; ++pj1; ++px1; \
2035 SEXP x1 = PROTECT(duplicate(x0));
2054 static const char *
valid[] = {
2056 int ivalid = R_check_class_etc(from,
valid);
2065 if (
class[0] !=
'z' &&
class[0] !=
'd') {
2072 PROTECT_WITH_INDEX(from, &pid);
2074 char cl[] =
"...Matrix";
2075 cl[0] = (
class[0] !=
'z') ?
'd' :
'z';
2076 cl[1] = (
class[1] !=
's') ?
'g' :
's';
2081 int *pdim = INTEGER(dim), n = pdim[0];
2083 error(
_(
"attempt to get skew-symmetric part of non-square matrix"));
2089 if (
class[1] ==
's')
2095 if (
class[1] ==
's') {
2098 char ul = *CHAR(STRING_ELT(uplo, 0));
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);
2113 if (
class[2] !=
'T') {
2118 if (
class[2] !=
'R') {
2123 if (
class[2] !=
'C') {
2129 x1 = PROTECT(duplicate(x0));
2135 }
else if (
class[2] !=
'T') {
2139 i0 = PROTECT(
GET_SLOT(from, iSym)),
2141 p1 = PROTECT(allocVector(INTSXP, (R_xlen_t) n + 1));
2142 int j, k, kend, k_, kend_, *pp0 = INTEGER(p0), *pi0 = INTEGER(i0),
2144 pp0++; *(pp1++) = 0;
2150 i0_ = PROTECT(
GET_SLOT(from, iSym)),
2152 int *pp0_ = INTEGER(p0_), *pi0_ = INTEGER(i0_), *pp1_;
2156 for (j = 0, k = 0, k_ = 0; j < n; ++j) {
2163 while (k_ < kend_ && pi0_[k_] < pi0[k]) {
2170 if (k_ < kend_ && pi0_[k_] == pi0[k])
2175 while (k_ < kend_) {
2186 for (j = 0; j < n; ++j) {
2187 pp1[j] = pp1[j - 1] + pp1_[j];
2188 pp1_[j] = pp1[j - 1];
2191 SEXP i1 = PROTECT(allocVector(INTSXP, pp1[n - 1])),
2192 x1 = PROTECT(allocVector(
kindToType(
cl[0]), pp1[n - 1]));
2193 int *pi1 = INTEGER(i1);
2196#define SP_LOOP(_CTYPE_, _PTR_, _ASSIGN_, _INCREMENT_) \
2198 _CTYPE_ *px0 = _PTR_(x0), *px1 = _PTR_(x1), \
2199 *px0_ = _PTR_(x0_); \
2200 for (j = 0, k = 0, k_ = 0; j < n; ++j) { \
2203 while (k < kend) { \
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]]); \
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_]); \
2222 pi1[pp1_[pi0[k]]] = j; \
2223 _ASSIGN_(px1[pp1_[pi0[k]]], -px1[pp1_[j]]); \
2229 while (k_ < kend_) { \
2230 if (pi0_[k_] >= j) \
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]]); \
2261 int *pi0 = INTEGER(i0), *pj0 = INTEGER(j0);
2262 R_xlen_t k, nnz0 = XLENGTH(i0), nnz1 = nnz0;
2264 for (k = 0; k < nnz0; ++k)
2265 if (pi0[k] == pj0[k])
2269 SEXP i1 = PROTECT(allocVector(INTSXP, nnz1)),
2270 j1 = PROTECT(allocVector(INTSXP, nnz1)),
2272 int *pi1 = INTEGER(i1), *pj1 = INTEGER(j1);
2275#define SP_LOOP(_CTYPE_, _PTR_, _ASSIGN_) \
2277 _CTYPE_ *px0 = _PTR_(x0), *px1 = _PTR_(x1); \
2278 for (k = 0; k < nnz0; ++k) { \
2279 if (*pi0 != *pj0) { \
2282 _ASSIGN_((*px1), 0.5 * (*px0)); \
2283 ++pi1; ++pj1; ++px1; \
2286 _ASSIGN_((*px1), -0.5 * (*px0)); \
2287 ++pi1; ++pj1; ++px1; \
2289 ++pi0; ++pj0; ++px0; \
2314 static const char *
valid[] = {
2316 int ivalid = R_check_class_etc(from,
valid);
2325 if (
class[1] ==
's')
2334 if (
class[1] ==
't')
2338 int *pdim = INTEGER(dim), n = pdim[0];
2344 if (
class[2] ==
'T') {
2354 int i, j, k, kend, *pp_, *pp0 = INTEGER(p0) + 1, *pi0 = INTEGER(i0);
2360#define IS_LOOP(_CTYPE_, _PTR_, _MASK_, _NOTREAL_, _NOTCONJ_) \
2362 _MASK_(_CTYPE_ *px0 = _PTR_(x0)); \
2363 for (j = 0, k = 0; j < n; ++j) { \
2365 while (k < kend) { \
2369 if (_NOTREAL_(px0[k])) \
2375 if (pp_[i] == pp0[i] || pi0[pp_[i]] != j || \
2376 _NOTCONJ_(px0[k], px0[pp_[i]])) \
2386#undef NOTCONJ_PATTERN
2387#define NOTCONJ_PATTERN(_X_, _Y_) 0
2394 if (
class[0] ==
'n')
2419 for (j = 0; j < n; ++j)
2420 if (pp_[j] != pp0[j])
2440 static const char *
valid[] = {
2442 int ivalid = R_check_class_etc(obj,
valid);
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");
2456 if (
class[1] ==
't') {
2458 char ul = *CHAR(STRING_ELT(uplo, 0));
2459 if (upper == NA_LOGICAL || (upper != 0) == (ul ==
'U'))
2460 return (ul ==
'U') ? 1 : -1;
2462 return (ul ==
'U') ? -1 : 1;
2467 if (
class[1] ==
's') {
2471 char ul = *CHAR(STRING_ELT(uplo, 0));
2472 if (upper == NA_LOGICAL)
2473 return (ul ==
'U') ? 1 : -1;
2475 return (upper != 0) ? 1 : -1;
2479 int *pdim = INTEGER(dim), n = pdim[0];
2483 return (upper != 0) ? 1 : -1;
2485 if (
class[2] !=
'T') {
2491 int j, k, kend, *pp0 = INTEGER(p0) + 1, *pi0 = INTEGER(i0);
2493 if (upper == NA_LOGICAL) {
2495 for (j = 0, k = 0; j < n; ++j) {
2497 if (k < kend && pi0[kend - 1] > j)
2502 return (
class[2] ==
'C') ? 1 : -1;
2504 for (j = 0, k = 0; j < n; ++j) {
2506 if (k < kend && pi0[k] < j)
2511 return (
class[2] ==
'C') ? -1 : 1;
2513 }
else if ((
class[2] ==
'C') == (upper != 0)) {
2515 for (j = 0, k = 0; j < n; ++j) {
2517 if (k < kend && pi0[kend - 1] > j)
2521 return (
class[2] ==
'C') ? 1 : -1;
2524 for (j = 0, k = 0; j < n; ++j) {
2526 if (k < kend && pi0[k] < j)
2530 return (
class[2] ==
'C') ? -1 : 1;
2538 int *pi0 = INTEGER(i0), *pj0 = INTEGER(j0);
2539 R_xlen_t k, nnz0 = XLENGTH(i0);
2541 if (upper == NA_LOGICAL) {
2542 for (k = 0; k < nnz0; ++k)
2543 if (pi0[k] > pj0[k])
2547 for (k = 0; k < nnz0; ++k)
2548 if (pi0[k] < pj0[k])
2553 }
else if (upper != 0) {
2554 for (k = 0; k < nnz0; ++k)
2555 if (pi0[k] > pj0[k])
2559 for (k = 0; k < nnz0; ++k)
2560 if (pi0[k] < pj0[k])
2573 static const char *
valid[] = {
2575 int ivalid = R_check_class_etc(obj,
valid);
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];
2584 SEXP ans = allocVector(LGLSXP, 1);
2585 LOGICAL(ans)[0] = ans_ != 0;
2586 if (upper_ == NA_LOGICAL && ans_ != 0) {
2589 SEXP kindSym = NULL;
2590 SEXP kindVal = PROTECT(mkString((ans_ > 0) ?
"U" :
"L"));
2591 if (!kindSym) kindSym = install(
"kind");
2592 setAttrib(ans, kindSym, kindVal);
2601 int *pdim = INTEGER(dim), n = pdim[0];
2607 if (
class[2] !=
'T') {
2613 int j, k, kend, *pp0 = INTEGER(p0) + 1, *pi0 = INTEGER(i0);
2615 for (j = 0, k = 0; j < n; ++j) {
2617 if (kend - k > 1 || (kend - k == 1 && pi0[k] != j))
2628 int *pi0 = INTEGER(i0), *pj0 = INTEGER(j0);
2629 R_xlen_t k, nnz0 = XLENGTH(i0);
2631 for (k = 0; k < nnz0; ++k)
2632 if (*(pi0++) != *(pj0++))
2644 static const char *
valid[] = {
2646 int ivalid = R_check_class_etc(obj,
valid);
2653#define MAP(_I_) work[_I_]
2654#define NOMAP(_I_) _I_
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_
2662#define SUM_CASES(_MAP_) \
2664 if (class[0] == 'n') { \
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); \
2670 SUM_LOOP(int, LOGICAL, int, INTEGER, HIDE, \
2671 0, 1, NA_INTEGER, ISNA_PATTERN, \
2672 _MAP_, CAST_PATTERN, INCREMENT_INTEGER, SCALE2_REAL); \
2674 SEXP x0 = PROTECT(GET_SLOT(obj, Matrix_xSym)); \
2675 switch (class[0]) { \
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); \
2682 SUM_LOOP(int, LOGICAL, int, INTEGER, SHOW, \
2683 0, 1, NA_INTEGER, ISNA_LOGICAL, \
2684 _MAP_, CAST_LOGICAL, INCREMENT_INTEGER, SCALE2_REAL); \
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); \
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); \
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); \
2708#define SUM_TYPEOF(c) (c == 'z') ? CPLXSXP : ((mean || c == 'd' || c == 'i') ? REALSXP : INTSXP)
2712 int m,
int n,
char di,
int narm,
int mean,
2715 int narm_ = narm && mean &&
class[0] !=
'n';
2718 int *pp0 = INTEGER(p0) + 1, j, k, kend, nnz1 = n, count = -1;
2720 if (IS_S4_OBJECT(res)) {
2724 for (j = 0; j < n; ++j)
2725 if (pp0[j - 1] < pp0[j])
2730 j1 = PROTECT(allocVector(INTSXP, nnz1)),
2731 x1 = PROTECT(allocVector(
SUM_TYPEOF(
class[0]), nnz1));
2735 int *pj1 = INTEGER(j1);
2737 for (j = 0; j < n; ++j)
2740 for (j = 0; j < n; ++j)
2741 if (pp0[j - 1] < pp0[j])
2744#define SUM_LOOP(_CTYPE0_, _PTR0_, _CTYPE1_, _PTR1_, _MASK_, \
2745 _ZERO_, _ONE_, _NA_, _ISNA_, \
2746 _MAP_, _CAST_, _INCREMENT_, _SCALE2_) \
2748 _MASK_(_CTYPE0_ *px0 = _PTR0_(x0)); \
2749 _CTYPE1_ *px1 = _PTR1_(x1) , tmp; \
2750 for (j = 0, k = 0; j < n; ++j) { \
2752 if (k < kend || nnz1 == n) { \
2753 *px1 = (di != 'N') ? _ONE_ : _ZERO_; \
2756 while (k < kend) { \
2757 if (_ISNA_(*px0)) { \
2763 tmp = _CAST_(*px0); \
2764 _INCREMENT_((*px1), tmp); \
2770 _SCALE2_((*px1), count); \
2794 int m,
int n,
char di,
int narm,
int mean,
2795 SEXP res, SEXP iSym)
2797 int narm_ = narm && mean &&
class[0] !=
'n';
2801 int *pp0 = INTEGER(p0) + 1, *pi0 = INTEGER(i0), i, j, k, kend,
2802 nnz0 = pp0[n - 1], nnz1 = m;
2804 if (IS_S4_OBJECT(res)) {
2809 for (i = 0; i < m; ++i)
2812 if (
class[1] !=
's') {
2813 for (k = 0; k < nnz0; ++k)
2816 for (j = 0, k = 0; j < n; ++j) {
2827 for (i = 0; i < m; ++i)
2828 work[i] = (work[i]) ? nnz1++ : -1;
2832 i1 = PROTECT(allocVector(INTSXP, nnz1)),
2833 x1 = PROTECT(allocVector(
SUM_TYPEOF(
class[0]), nnz1));
2836 int *pi1 = INTEGER(i1);
2838 for (i = 0; i < nnz1; ++i)
2841#define SUM_LOOP(_CTYPE0_, _PTR0_, _CTYPE1_, _PTR1_, _MASK_, \
2842 _ZERO_, _ONE_, _NA_, _ISNA_, \
2843 _MAP_, _CAST_, _INCREMENT_, _SCALE2_) \
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) \
2850 if (class[1] != 's') { \
2851 for (k = 0; k < nnz0; ++k) { \
2852 if (_ISNA_(px0[k])) { \
2854 px1[_MAP_(pi0[k])] = _NA_; \
2856 --pi1[_MAP_(pi0[k])]; \
2858 tmp = _CAST_(px0[k]); \
2859 _INCREMENT_(px1[_MAP_(pi0[k])], tmp); \
2864 for (j = 0, k = 0; j < n; ++j) { \
2866 while (k < kend) { \
2867 off = pi0[k] != j; \
2868 if (_ISNA_(px0[k])) { \
2870 px1[_MAP_(pi0[k])] = _NA_; \
2872 px1[_MAP_( j)] = _NA_; \
2873 } else if (narm_) { \
2874 --pi1[_MAP_(pi0[k])]; \
2879 tmp = _CAST_(px0[k]); \
2880 _INCREMENT_(px1[_MAP_(pi0[k])], tmp); \
2882 _INCREMENT_(px1[_MAP_( j)], tmp); \
2890 for (i = 0; i < nnz1; ++i) \
2891 _SCALE2_(px1[i], pi1[i]); \
2893 for (i = 0; i < nnz1; ++i) \
2894 _SCALE2_(px1[i], n); \
2899 for (i = 0; i < m; ++i)
2911 for (i = 0; i < m; ++i)
2928 int m,
int n,
char di,
int narm,
int mean,
2929 SEXP res, SEXP iSym, SEXP jSym)
2931 int narm_ = narm && mean &&
class[0] !=
'n';
2936 SEXP i0 = PROTECT(
GET_SLOT(obj, iSym)),
2938 int *pi0 = INTEGER(i0), *pj0 = INTEGER(j0), j, nnz1 = n;
2939 R_xlen_t k, nnz0 = XLENGTH(i0);
2941 if (IS_S4_OBJECT(res)) {
2946 for (j = 0; j < n; ++j)
2949 if (
class[1] !=
's') {
2950 for (k = 0; k < nnz0; ++k)
2953 for (k = 0; k < nnz0; ++k) {
2955 if (pi0[k] != pj0[k])
2960 for (j = 0; j < n; ++j)
2961 work[j] = (work[j]) ? nnz1++ : -1;
2965 j1 = PROTECT(allocVector(INTSXP, nnz1)),
2966 x1 = PROTECT(allocVector(
SUM_TYPEOF(
class[0]), nnz1));
2969 int *pj1 = INTEGER(j1);
2971 for (j = 0; j < nnz1; ++j)
2974#define SUM_LOOP(_CTYPE0_, _PTR0_, _CTYPE1_, _PTR1_, _MASK_, \
2975 _ZERO_, _ONE_, _NA_, _ISNA_, \
2976 _MAP_, _CAST_, _INCREMENT_, _SCALE2_) \
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) \
2983 if (class[1] != 's') { \
2984 for (k = 0; k < nnz0; ++k) { \
2985 if (_ISNA_(px0[k])) { \
2987 px1[_MAP_(pj0[k])] = _NA_; \
2989 --pj1[_MAP_(pj0[k])]; \
2991 tmp = _CAST_(px0[k]); \
2992 _INCREMENT_(px1[_MAP_(pj0[k])], tmp); \
2997 for (k = 0; k < nnz0; ++k) { \
2998 off = pi0[k] != pj0[k]; \
2999 if (_ISNA_(px0[k])) { \
3001 px1[_MAP_(pj0[k])] = _NA_; \
3003 px1[_MAP_(pi0[k])] = _NA_; \
3004 } else if (narm_) { \
3005 --pj1[_MAP_(pj0[k])]; \
3007 --pj1[_MAP_(pi0[k])]; \
3010 tmp = _CAST_(px0[k]); \
3011 _INCREMENT_(px1[_MAP_(pj0[k])], tmp); \
3013 _INCREMENT_(px1[_MAP_(pi0[k])], tmp); \
3019 for (j = 0; j < nnz1; ++j) \
3020 _SCALE2_(px1[j], pj1[j]); \
3022 for (j = 0; j < nnz1; ++j) \
3023 _SCALE2_(px1[j], m); \
3028 for (j = 0; j < n; ++j)
3053 int narm,
int mean,
int sparse)
3056 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1],
3057 r = (margin == 0) ? m : n;
3062 char cl[] =
".sparseVector";
3063 cl[0] = (type == CPLXSXP) ?
'z' : ((type == REALSXP) ?
'd' :
'i');
3066 SEXP length = PROTECT(ScalarInteger(r));
3070 PROTECT(res = allocVector(type, r));
3072 SEXP dimnames = (
class[1] !=
's')
3075 marnames = VECTOR_ELT(dimnames, margin);
3076 if (marnames != R_NilValue) {
3078 setAttrib(res, R_NamesSymbol, marnames);
3084 if (
class[1] ==
't') {
3086 di = *CHAR(STRING_ELT(diag, 0));
3090 if (
class[2] ==
'C') {
3093 }
else if (
class[2] ==
'R') {
3094 if (
class[1] !=
's')
3104 if (
class[2] ==
'C') {
3105 if (
class[1] !=
's')
3110 }
else if (
class[2] ==
'R') {
3125 SEXP narm, SEXP mean, SEXP sparse)
3127 static const char *
valid[] = {
3129 int ivalid = R_check_class_etc(obj,
valid);
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);
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");
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");
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");
3154 margin_, narm_, mean_, sparse_);
3160#define TRY_INCREMENT(_LABEL_) \
3163 ? ( t <= MATRIX_INT_FAST64_MAX - s) \
3164 : (-t <= s - MATRIX_INT_FAST64_MIN)) { \
3174#define LONGDOUBLE_AS_DOUBLE(v) \
3175 (v > DBL_MAX) ? R_PosInf : ((v < -DBL_MAX) ? R_NegInf : (double) v);
3179 if (
class[2] ==
'T')
3185 if (!narm && (
class[0] ==
'l' ||
class[0] ==
'i')) {
3187 int *px = (
class[0] ==
'l') ? LOGICAL(x) : INTEGER(x);
3188 R_xlen_t nx = XLENGTH(x);
3190 if (*px == NA_INTEGER) {
3191 res = allocVector(INTSXP, 1);
3192 INTEGER(res)[0] = NA_INTEGER;
3201 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1];
3204 if (
class[1] ==
't') {
3206 di = *CHAR(STRING_ELT(diag, 0));
3209 int symmetric =
class[1] ==
's';
3211 if (
class[2] !=
'T') {
3216 int *pp = INTEGER(p) + 1, *pi = INTEGER(i), j_, k = 0, kend,
3217 n_ = (
class[2] ==
'C') ? n : m;
3219 if (
class[0] ==
'n') {
3225 char ul = *CHAR(STRING_ELT(uplo, 0));
3228 for (j_ = 0; j_ < n_; ++j_) {
3230 if (k < kend && pi[(ul ==
'U') ? kend - 1 : k] == j_)
3235 if (nnz <= INT_MAX) {
3236 res = allocVector(INTSXP, 1);
3237 INTEGER(res)[0] = (int) nnz;
3239 res = allocVector(REALSXP, 1);
3240 REAL(res)[0] = (double) nnz;
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_) {
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;
3264 res = allocVector(CPLXSXP, 1);
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_) {
3273 if (!(narm && ISNAN(px[k])))
3274 zr += (symmetric && pi[k] != j_)
3275 ? 2.0L * px[k] : px[k];
3279 res = allocVector(REALSXP, 1);
3282 int *px = (
class[0] ==
'l') ? LOGICAL(x) : INTEGER(x);
3284 unsigned int count = 0;
3286 for (j_ = 0; j_ < n_; ++j_) {
3289 if (!narm || px[k] != NA_INTEGER) {
3290 int d = (symmetric && pi[k] != j_) ? 2 : 1;
3291 if (count > UINT_MAX - d)
3293 t += (d == 2) ? 2LL * px[k] : px[k];
3302 long double zr = (
long double) s + (
long double) t;
3303 for (; j_ < n_; ++j_) {
3306 if (!narm || px[k] != NA_INTEGER)
3307 zr += (symmetric && pi[k] != j_)
3308 ? 2.0L * px[k] : px[k];
3312 res = allocVector(REALSXP, 1);
3314 }
else if (s > INT_MIN && s <= INT_MAX) {
3315 res = allocVector(INTSXP, 1);
3316 INTEGER(res)[0] = (int) s;
3318 res = allocVector(REALSXP, 1);
3319 REAL(res)[0] = (double) s;
3327 int *pi = INTEGER(i), *pj = INTEGER(j);
3328 R_xlen_t k, kend = XLENGTH(i);
3330 if (
class[0] ==
'n') {
3336 for (k = 0; k < kend; ++k)
3340 if (nnz <= INT_MAX) {
3341 res = allocVector(INTSXP, 1);
3342 INTEGER(res)[0] = (int) nnz;
3344 res = allocVector(REALSXP, 1);
3345 REAL(res)[0] = (double) nnz;
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;
3364 res = allocVector(CPLXSXP, 1);
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);
3377 int *px = (
class[0] ==
'i') ? INTEGER(x) : LOGICAL(x);
3379 unsigned int count = 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)
3386 t += (d == 2) ? 2LL * px[k] : px[k];
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);
3400 }
else if (s > INT_MIN && s <= INT_MAX) {
3401 res = allocVector(INTSXP, 1);
3402 INTEGER(res)[0] = (int) s;
3404 res = allocVector(REALSXP, 1);
3405 REAL(res)[0] = (double) s;
3418 static const char *
valid[] = {
3420 int ivalid = R_check_class_etc(obj,
valid);
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");
3434 if (
class[2] ==
'T')
3438 SEXP res = PROTECT(allocVector((
class[0] ==
'z') ? CPLXSXP : REALSXP, 1));
3441 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1];
3443 char ul =
'U', di =
'N';
3444 if (
class[1] !=
'g') {
3446 ul = *CHAR(STRING_ELT(uplo, 0));
3447 if (
class[1] ==
't') {
3449 di = *CHAR(STRING_ELT(diag, 0));
3453 int symmetric = (
class[1] !=
's')
3454 ? 0 : (((
class[2] ==
'C') == (ul ==
'U')) ? 1 : -1);
3455 long double zr = 1.0L, zi = 0.0L;
3458 nnz, nnzmax = (symmetric) ? (mn + n) / 2 : mn;
3460 if (
class[2] !=
'T') {
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,
3472 if (
class[0] ==
'n') {
3473 REAL(res)[0] = (nnz == nnzmax) ? 1.0 : 0.0;
3481 if (
class[0] ==
'z') {
3482 Rcomplex *px = COMPLEX(x);
3483 long double zr0, zi0;
3484 for (j_ = 0; j_ < n_; ++j_) {
3486 if (seen0 || kend - k == m_) {
3488 if (!(narm && (ISNAN(px[k].r) || ISNAN(px[k].i)))) {
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_) {
3494 zr = zr0 * px[k].r - zi0 * px[k].i;
3495 zi = zr0 * px[k].i + zi0 * px[k].r;
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_)) {
3507 if (!(narm && (ISNAN(px[k].r) || ISNAN(px[k].i)))) {
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_) {
3513 zr = zr0 * px[k].r - zi0 * px[k].i;
3514 zi = zr0 * px[k].i + zi0 * px[k].r;
3518 }
else if (di ==
'N' || i_ != j_) {
3526 }
else if (
class[0] ==
'd') {
3527 double *px = REAL(x);
3528 for (j_ = 0; j_ < n_; ++j_) {
3530 if (seen0 || kend - k == m_) {
3532 if (!(narm && ISNAN(px[k])))
3533 zr *= (symmetric && pi[k] != j_)
3534 ? (
long double) px[k] * px[k] : px[k];
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_)) {
3544 if (!(narm && ISNAN(px[k])))
3545 zr *= (symmetric && pi[k] != j_)
3546 ? (
long double) px[k] * px[k] : px[k];
3548 }
else if (di ==
'N' || i_ != j_) {
3556 int *px = (
class[0] ==
'l') ? LOGICAL(x) : INTEGER(x);
3557 for (j_ = 0; j_ < n_; ++j_) {
3559 if (seen0 || kend - k == m_) {
3561 if (px[k] != NA_INTEGER)
3562 zr *= (symmetric && pi[k] != j_)
3563 ? (
long double) px[k] * px[k] : px[k];
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_)) {
3575 if (px[k] != NA_INTEGER)
3576 zr *= (symmetric && pi[k] != j_)
3577 ? (
long double) px[k] * px[k] : px[k];
3581 }
else if (di ==
'N' || i_ != j_) {
3594 int *pi = INTEGER(i), *pj = INTEGER(j);
3595 R_xlen_t k, kend = XLENGTH(i);
3600 if (
class[0] ==
'n') {
3601 REAL(res)[0] = (nnz == nnzmax) ? 1.0 : 0.0;
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)))) {
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]) {
3621 zr = zr0 * px[k].r - zi0 * px[k].i;
3622 zi = zr0 * px[k].i + zi0 * px[k].r;
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];
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];
3643 if (
class[0] ==
'z') {
3655 static const char *
valid[] = {
3657 int ivalid = R_check_class_etc(obj,
valid);
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");
3670#undef LONGDOUBLE_AS_DOUBLE
3675 int ivalid = R_check_class_etc(from,
valid);
3678 const char *
cl =
valid[ivalid];
3681 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1];
3687 i1 = NULL, j1 = NULL;
3690 void taggr(SEXP, SEXP, SEXP, SEXP *, SEXP *, SEXP *,
int,
int);
3693 taggr(j0, i0, NULL, &j1, &i1, NULL, n, m);
3707 taggr(j0, i0, x0, &j1, &i1, &x1, n, m);
3724 if (m != n || n > 0) {
3726 pdim = INTEGER(dim);
3738 char ul = *CHAR(STRING_ELT(uplo, 0));
3745 char di = *CHAR(STRING_ELT(diag, 0));
3751 if (LENGTH(factors) > 0)
long long Matrix_int_fast64_t
#define NOTCONJ_LOGICAL(_X_, _Y_)
#define NOTREAL_LOGICAL(_X_)
#define ERROR_INVALID_CLASS(_X_, _FUNC_)
#define NOTREAL_INTEGER(_X_)
#define NOTREAL_PATTERN(_X_)
#define NOTREAL_COMPLEX(_X_)
#define NOTREAL_REAL(_X_)
#define INCREMENT_REAL(_X_, _Y_)
char typeToKind(SEXPTYPE)
#define ASSIGN_COMPLEX(_X_, _Y_)
#define INCREMENT_COMPLEX(_X_, _Y_)
#define ASSIGN_REAL(_X_, _Y_)
SEXPTYPE kindToType(char)
#define NOTCONJ_REAL(_X_, _Y_)
#define NOTCONJ_COMPLEX(_X_, _Y_)
#define ISNZ_LOGICAL(_X_)
#define Matrix_Free(_VAR_, _N_)
#define INCREMENT_PATTERN(_X_, _Y_)
#define SET_SLOT(x, what, value)
#define GET_SLOT(x, what)
#define Matrix_Calloc(_VAR_, _N_, _CTYPE_)
SEXP newObject(const char *)
#define NOTCONJ_INTEGER(_X_, _Y_)
SEXP get_symmetrized_DimNames(SEXP obj, int J)
int DimNames_is_symmetric(SEXP dn)
void set_reversed_DimNames(SEXP obj, SEXP dn)
void set_symmetrized_DimNames(SEXP obj, SEXP dn, int J)
static const char * valid[]
SEXP sparse_as_Csparse(SEXP from, const char *class)
void taggr(SEXP i0, SEXP j0, SEXP x0, SEXP *i1, SEXP *j1, SEXP *x1, int m, int n)
SEXP sparse_as_kind(SEXP from, const char *class, char kind)
void trans(SEXP p0, SEXP i0, SEXP x0, SEXP p1, SEXP i1, SEXP x1, int m, int n)
#define FS_SUBCASES(_CTYPE_, _PTR_, _MASK_, _ONE_)
#define DROP0_LOOP2(_CTYPE_, _PTR_, _ISNZ_)
SEXP R_sparse_band(SEXP from, SEXP k1, SEXP k2)
SEXP sparse_sum(SEXP obj, const char *class, int narm)
#define NOTCONJ_PATTERN(_X_, _Y_)
SEXP R_sparse_is_diagonal(SEXP obj)
static void Csparse_colsum(SEXP obj, const char *class, int m, int n, char di, int narm, int mean, SEXP res)
SEXP sparse_diag_N2U(SEXP from, const char *class)
SEXP R_sparse_skewpart(SEXP from)
SEXP sparse_force_symmetric(SEXP from, const char *class, char ul)
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)
#define DROP0_LOOP1(_CTYPE_, _PTR_, _ISNZ_)
SEXP sparse_diag_U2N(SEXP from, const char *class)
SEXP R_sparse_diag_U2N(SEXP from)
SEXP R_sparse_diag_N2U(SEXP from)
static void Csparse_rowsum(SEXP obj, const char *class, int m, int n, char di, int narm, int mean, SEXP res, SEXP iSym)
#define DS_LOOP(_CTYPE_, _PTR_, _MASK_, _ISNZ_)
#define DROP0_CASES(_DO_)
SEXP R_sparse_transpose(SEXP from, SEXP lazy)
#define SP_LOOP(_CTYPE_, _PTR_, _ASSIGN_, _INCREMENT_)
SEXP R_sparse_drop0(SEXP from, SEXP tol)
SEXP R_sparse_sum(SEXP obj, SEXP narm)
int sparse_is_triangular(SEXP obj, const char *class, int upper)
SEXP R_sparse_force_symmetric(SEXP from, SEXP uplo)
#define LONGDOUBLE_AS_DOUBLE(v)
SEXP sparse_transpose(SEXP from, const char *class, int lazy)
SEXP R_sparse_prod(SEXP obj, SEXP narm)
SEXP sparse_symmpart(SEXP from, const char *class)
SEXP R_sparse_is_triangular(SEXP obj, SEXP upper)
SEXP R_sparse_symmpart(SEXP from)
SEXP R_sparse_diag_set(SEXP from, SEXP value)
SEXP R_sparse_diag_get(SEXP obj, SEXP names)
SEXP sparse_band(SEXP from, const char *class, int a, int b)
#define DG_LOOP(_CTYPE_, _PTR_, _MASK_, _REPLACE_, _INCREMENT_, _ZERO_, _ONE_)
SEXP sparse_diag_set(SEXP from, const char *class, SEXP value)
#define BAND_SUBCASES(_CTYPE_, _PTR_, _MASK_)
#define TRY_INCREMENT(_LABEL_)
SEXP sparse_marginsum(SEXP obj, const char *class, int margin, int narm, int mean, int sparse)
#define IS_LOOP(_CTYPE_, _PTR_, _MASK_, _NOTREAL_, _NOTCONJ_)
SEXP Tsparse_aggregate(SEXP from)
int sparse_is_symmetric(SEXP obj, const char *class, int checkDN)
int sparse_is_diagonal(SEXP obj, const char *class)
SEXP sparse_skewpart(SEXP from, const char *class)
SEXP sparse_drop0(SEXP from, const char *class, double tol)
SEXP sparse_diag_get(SEXP obj, const char *class, int names)
SEXP R_sparse_is_symmetric(SEXP obj, SEXP checkDN)
SEXP R_sparse_marginsum(SEXP obj, SEXP margin, SEXP narm, SEXP mean, SEXP sparse)
SEXP sparse_prod(SEXP obj, const char *class, int narm)
void * Matrix_memset(void *dest, int ch, R_xlen_t length, size_t size)
int equal_character_vectors(SEXP s1, SEXP s2, int n)
void * Matrix_memcpy(void *dest, const void *src, R_xlen_t length, size_t size)