10void scanArgs(SEXP args, SEXP exprs,
int margin,
int level,
11 int *rdim,
int *rdimnames,
char *kind,
char *repr)
14 int nS4 = 0, nDense = 0,
15 anyCsparse = 0, anyRsparse = 0, anyTsparse = 0, anyDiagonal = 0,
16 anyN = 0, anyL = 0, anyI = 0, anyD = 0, anyZ = 0,
23 rdimnames[0] = rdimnames[1] = 0;
25 for (a = args; a != R_NilValue; a = CDR(a)) {
35 if (rdim[!margin] < 0)
36 rdim[!margin] = sdim[!margin];
37 else if (sdim[!margin] != rdim[!margin]) {
39 Rf_error(
_(
"number of rows of matrices must match"));
41 Rf_error(
_(
"number of columns of matrices must match"));
43 if (sdim[margin] > INT_MAX - rdim[margin])
44 Rf_error(
_(
"dimensions cannot exceed %s"),
"2^31-1");
45 rdim[margin] += sdim[margin];
47 if (!rdimnames[0] || !rdimnames[1]) {
50 if (VECTOR_ELT(tmp, 0) != R_NilValue ||
51 VECTOR_ELT(tmp, 1) != R_NilValue)
52 rdimnames[0] = rdimnames[1] = 1;
54 for (i = 0; i < 2; ++i)
56 VECTOR_ELT(tmp, i) != R_NilValue)
106 if (
MARGIN(s) != margin) {
136 tmp = Rf_getAttrib(s, R_DimSymbol);
137 if (
TYPEOF(tmp) == INTSXP && LENGTH(tmp) == 2) {
139 if (rdim[!margin] < 0)
140 rdim[!margin] = sdim[!margin];
141 else if (rdim[!margin] != sdim[!margin]) {
143 Rf_error(
_(
"number of rows of matrices must match"));
145 Rf_error(
_(
"number of columns of matrices must match"));
147 if (sdim[margin] > INT_MAX - rdim[margin])
148 Rf_error(
_(
"dimensions cannot exceed %s"),
"2^31-1");
149 rdim[margin] += sdim[margin];
151 if (!rdimnames[0] || !rdimnames[1]) {
152 tmp = Rf_getAttrib(s, R_DimNamesSymbol);
153 if (tmp != R_NilValue)
154 for (i = 0; i < 2; ++i)
156 VECTOR_ELT(tmp, i) != R_NilValue)
163 if (rdim[!margin] < 0) {
165 R_xlen_t maxlen = -1;
166 for (a = args; a != R_NilValue; a = CDR(a)) {
172 Rf_error(
_(
"dimensions cannot exceed %s"),
"2^31-1");
173 else if (slen > maxlen)
179 rdim[!margin] = (int) maxlen;
182 for (a = args, e = exprs; a != R_NilValue; a = CDR(a), e = CDR(e)) {
184 if ((s == R_NilValue && rdim[!margin] > 0) ||
TYPEOF(s) == OBJSXP)
189 tmp = Rf_getAttrib(s, R_DimSymbol);
190 if (
TYPEOF(tmp) == INTSXP && LENGTH(tmp) == 2)
193 if (slen == 0 && rdim[!margin] > 0)
195 if (rdim[margin] == INT_MAX)
196 Rf_error(
_(
"dimensions cannot exceed %s"),
"2^31-1");
198 if (slen > rdim[!margin] || rdim[!margin] % (
int) slen) {
200 Rf_warning(
_(
"number of rows of result is not a multiple of vector length"));
202 Rf_warning(
_(
"number of columns of result is not a multiple of vector length"));
204 if (!rdimnames[!margin] && slen == rdim[!margin]) {
205 tmp = Rf_getAttrib(s, R_NamesSymbol);
206 if (tmp != R_NilValue)
207 rdimnames[!margin] = 1;
210 if (!rdimnames[margin]) {
211 if (TAG(a) != R_NilValue ||
212 level == 2 || (level == 1 &&
TYPEOF(CAR(e)) == SYMSXP))
213 rdimnames[margin] = 1;
219#ifndef MATRIX_ENABLE_IMATRIX
220 else if (anyD || anyI)
237 else if (nDense == 0) {
238 if (anyCsparse && anyRsparse)
239 *repr = (margin) ?
'C' :
'R';
246 else if (anyDiagonal)
247 *repr = (margin) ?
'C' :
'R';
255 int_fast64_t nnz = 0, len = 0, snnz = 0, slen = 0;
256 for (a = args; a != R_NilValue && nnz < INT_MAX; a = CDR(a)) {
264 slen = (int_fast64_t) sdim[0] * sdim[1];
271 snnz = (scl[1] !=
't') ? slen : sdim[0] + (slen - sdim[0]) / 2;
277 int *pp = INTEGER(p), n = sdim[(scl[2] ==
'C') ? 1 : 0];
282 int *pi = INTEGER(i), j;
284 if (
UPLO(s) ==
'U') {
285 for (j = 0; j < n; ++j)
286 if (pp[j] < pp[j + 1] && pi[pp[j + 1] - 1] == j)
289 for (j = 0; j < n; ++j)
290 if (pp[j] < pp[j + 1] && pi[pp[j]] == j)
294 }
else if (scl[1] ==
't' &&
DIAG(s) !=
'N')
305 int *pi = INTEGER(i), *pj = INTEGER(j);
306 R_xlen_t k = XLENGTH(i);
309 if (*(pi++) == *(pj++))
312 }
else if (scl[1] ==
't' &&
DIAG(s) !=
'N')
332 if (nnz > INT_MAX || nnz > len / 2)
334 else if (anyCsparse && anyRsparse)
335 *repr = (margin) ?
'C' :
'R';
343 *repr = (margin) ?
'C' :
'R';
351 int *rdim,
char kind,
char repr)
355 char scl_[] =
"...Matrix";
358 for (a = args; a != R_NilValue; a = CDR(a)) {
361 SET_TAG(a, R_NilValue);
365 PROTECT_WITH_INDEX(s, &pid);
366 if (
TYPEOF(s) == OBJSXP) {
376 scl_[0] = kind; scl_[1] = scl[1]; scl_[2] = scl[2];
384 scl_[0] = scl[0]; scl_[1] = scl[1]; scl_[2] = repr;
397 scl_[0] = kind; scl_[1] = scl[1]; scl_[2] = scl[2];
449 tmp = Rf_getAttrib(s, R_DimSymbol);
450 isM =
TYPEOF(tmp) == INTSXP && LENGTH(tmp) == 2;
452 if (rdim[!margin] > 0 && XLENGTH(s) == 0) {
459 REPROTECT(s = Rf_coerceVector(s,
kindToType(kind)), pid);
461 if (!isM && XLENGTH(s) != rdim[!margin]) {
462 static SEXP replen = NULL;
464 replen = Rf_install(
"rep_len");
465 SEXP lengthout = PROTECT(Rf_ScalarInteger(rdim[!margin])),
466 call = PROTECT(Rf_lang3(replen, s, lengthout));
467 REPROTECT(s = Rf_eval(call, R_GlobalEnv), pid);
484 int *rdim,
char kind,
char repr)
490 if (rdim[0] == 0 || rdim[1] == 0)
494 int k, m = rdim[0], n = rdim[1];
495 R_xlen_t mn = (R_xlen_t) m * n;
499 SEXP x = PROTECT(Rf_allocVector(c##TYPESXP, mn)); \
500 c##TYPE *px = c##PTR(x), *ps; \
501 for (a = args; a != R_NilValue; a = CDR(a)) { \
503 if (s == R_NilValue) \
505 if (TYPEOF(s) != OBJSXP) \
506 tmp = Rf_getAttrib(s, R_DimSymbol); \
508 s = GET_SLOT(s, Matrix_xSym); \
514 if (!tmp || (TYPEOF(tmp) == INTSXP && LENGTH(tmp) == 2)) { \
515 memcpy(px, ps, sizeof(c##TYPE) * (size_t) mn); \
517 } else if (mn >= m) { \
518 memcpy(px, ps, sizeof(c##TYPE) * (size_t) m); \
520 } else if (mn == 1) { \
522 for (k = 0; k < m; ++k) \
525 int mn_ = (int) mn; \
526 for (k = 0; k < rdim[0]; ++k) \
527 *(px++) = ps[k % mn_]; \
531 if (!tmp || (TYPEOF(tmp) == INTSXP && LENGTH(tmp) == 2)) { \
532 m = (int) (mn / n); \
533 for (k = 0; k < n; ++k) { \
534 memcpy(py, ps, sizeof(c##TYPE) * (size_t) m); \
539 } else if (mn >= n) { \
540 for (k = 0; k < n; ++k) { \
546 } else if (mn == 1) { \
548 for (k = 0; k < n; ++k) { \
554 int mn_ = (int) mn; \
555 for (k = 0; k < n; ++k) { \
563 SET_SLOT(ans, Matrix_xSym, x); \
571 }
else if ((repr ==
'C' && margin) || (repr ==
'R' && !margin)) {
573 SEXP p = PROTECT(Rf_allocVector(INTSXP, (R_xlen_t) rdim[margin] + 1));
574 int *pp = INTEGER(p);
577 if (rdim[0] == 0 || rdim[1] == 0) {
578 memset(pp, 0,
sizeof(
int) * ((
size_t) rdim[margin] + 1));
584 int *psp, j, n, nnz = 0;
586 for (a = args; a != R_NilValue; a = CDR(a)) {
592 n = (int) (XLENGTH(sp) - 1);
593 if (psp[n] > INT_MAX - nnz)
594 Rf_error(
_(
"%s cannot exceed %s"),
"p[length(p)]",
"2^31-1");
595 for (j = 0; j < n; ++j)
596 *(pp++) = nnz = nnz + (psp[j + 1] - psp[j]);
599 SEXP i = PROTECT(Rf_allocVector(INTSXP, nnz)), si,
601 int *pi = INTEGER(i), *psi;
607 SEXP x = PROTECT(Rf_allocVector(kindToType(kind), nnz)), sx; \
608 c##TYPE *px = c##PTR(x), *psx; \
610 for (a = args; a != R_NilValue; a = CDR(a)) { \
612 if (s == R_NilValue) \
614 PROTECT(sp = GET_SLOT(s, Matrix_pSym)); \
615 PROTECT(si = GET_SLOT(s, iSym)); \
617 PROTECT(sx = GET_SLOT(s, Matrix_xSym)); \
624 n = (int) (XLENGTH(sp) - 1); \
625 memcpy(pi, psi, sizeof( int) * (size_t) psp[n]); \
627 memcpy(px, psx, sizeof(c##TYPE) * (size_t) psp[n]); \
639 SET_SLOT(ans, Matrix_xSym, x); \
650 }
else if ((repr ==
'C' && !margin) || (repr ==
'R' && margin)) {
652 SEXP p = PROTECT(Rf_allocVector(INTSXP, (R_xlen_t) rdim[!margin] + 1));
653 int *pp = INTEGER(p);
655 memset(pp, 0,
sizeof(
int) * ((
size_t) rdim[!margin] + 1));
657 if (rdim[0] == 0 || rdim[1] == 0) {
663 int *psp, j, n = rdim[!margin];
665 for (a = args; a != R_NilValue; a = CDR(a)) {
670 psp = INTEGER(sp) + 1;
671 if (n > 0 && psp[n - 1] > INT_MAX - pp[n - 1])
672 Rf_error(
_(
"%s cannot exceed %s"),
"p[length(p)]",
"2^31-1");
673 for (j = 0; j < n; ++j)
679 SEXP i = PROTECT(Rf_allocVector(INTSXP, nnz)), si,
681 int *pi = INTEGER(i), *psi, *work, k, kend, pos = 0;
684 memcpy(work, pp,
sizeof(
int) * (
size_t) n);
689 SEXP x = PROTECT(Rf_allocVector(c##TYPESXP, nnz)), sx; \
690 c##TYPE *px = c##PTR(x), *psx; \
692 for (a = args; a != R_NilValue; a = CDR(a)) { \
694 if (s == R_NilValue) \
696 PROTECT(sp = GET_SLOT(s, Matrix_pSym)); \
697 PROTECT(si = GET_SLOT(s, iSym)); \
699 PROTECT(sx = GET_SLOT(s, Matrix_xSym)); \
706 for (j = 0, k = 0; j < n; ++j) { \
709 pi[work[j]] = *(psi++) + pos; \
711 px[work[j]] = *(psx++); \
721 pos += INTEGER(GET_SLOT(s, Matrix_DimSym))[margin]; \
724 SET_SLOT(ans, Matrix_xSym, x); \
736 }
else if (repr ==
'T') {
738 if (rdim[0] == 0 || rdim[1] == 0)
742 for (a = args; a != R_NilValue; a = CDR(a)) {
747 if (k > R_XLEN_T_MAX - nnz)
748 Rf_error(
_(
"attempt to allocate vector of length exceeding %s"),
754 i = PROTECT(Rf_allocVector(INTSXP, nnz)),
755 j = PROTECT(Rf_allocVector(INTSXP, nnz));
756 int *psi, *psj, *pi = INTEGER(i), *pj = INTEGER(j), pos = 0;
763 SEXP x = PROTECT(Rf_allocVector(c##TYPESXP, nnz)), sx; \
764 c##TYPE *px = c##PTR(x), *psx; \
766 for (a = args; a != R_NilValue; a = CDR(a)) { \
768 if (s == R_NilValue) \
770 PROTECT(si = GET_SLOT(s, Matrix_iSym)); \
771 PROTECT(sj = GET_SLOT(s, Matrix_jSym)); \
773 PROTECT(sx = GET_SLOT(s, Matrix_xSym)); \
783 *(pi++) = *(psi++); \
784 *(pj++) = *(psj++) + pos; \
786 *(px++) = *(psx++); \
791 *(pi++) = *(psi++) + pos; \
792 *(pj++) = *(psj++); \
794 *(px++) = *(psx++); \
802 pos += INTEGER(GET_SLOT(s, Matrix_DimSym))[margin]; \
805 SET_SLOT(ans, Matrix_xSym, x); \
818 SEXP p = PROTECT(Rf_allocVector(INTSXP, rdim[margin])), sp;
819 int *pp = INTEGER(p);
820 for (a = args; a != R_NilValue; a = CDR(a)) {
825 memcpy(pp, INTEGER(sp),
sizeof(
int) * (
size_t) LENGTH(sp));
839SEXP
bind(SEXP args, SEXP exprs,
int margin,
int level)
844 int rdim[2], rdimnames[2];
845 char kind =
'\0', repr =
'\0';
846 scanArgs(args, exprs, margin, level,
847 rdim, rdimnames, &kind, &repr);
848 if (rdim[!margin] < 0)
851 if (repr ==
'e' && (int_fast64_t) rdim[0] * rdim[1] > R_XLEN_T_MAX)
852 Rf_error(
_(
"attempt to allocate vector of length exceeding %s"),
854 char rcl[] =
"...Matrix";
855 if (kind ==
'\0' || repr ==
'\0') {
857 Rf_error(
_(
"should never happen ..."));
868 bindArgs(args, margin, ans, rdim, kind, repr);
871 INTEGER(dim)[0] = rdim[0];
872 INTEGER(dim)[1] = rdim[1];
875 if (rdimnames[0] || rdimnames[1]) {
877 marnames = R_NilValue, nms[2], nms_, a, e, s, tmp;
878 int i, r = -1, pos = 0, nprotect = 1;
880 if (rdimnames[margin]) {
881 PROTECT(marnames = Rf_allocVector(STRSXP, rdim[margin]));
883 SET_VECTOR_ELT(dimnames, margin, marnames);
885 for (a = args, e = exprs; a != R_NilValue; a = CDR(a), e = CDR(e)) {
887 if (s == R_NilValue && rdim[!margin] > 0)
889 nms[0] = nms[1] = R_NilValue;
890 if (
TYPEOF(s) == OBJSXP) {
893 r = INTEGER(tmp)[margin];
896 if ((nms_ = VECTOR_ELT(tmp, 1)) != R_NilValue ||
897 (nms_ = VECTOR_ELT(tmp, 0)) != R_NilValue)
898 nms[0] = nms[1] = nms_;
900 for (i = 0; i < 2; ++i)
901 nms[i] = VECTOR_ELT(tmp, i);
903 tmp = Rf_getAttrib(s, R_DimSymbol);
904 if (
TYPEOF(tmp) == INTSXP && LENGTH(tmp) == 2) {
905 r = INTEGER(tmp)[margin];
906 tmp = Rf_getAttrib(s, R_DimNamesSymbol);
907 if (tmp != R_NilValue)
908 for (i = 0; i < 2; ++i)
909 nms[i] = VECTOR_ELT(tmp, i);
910 }
else if (rdim[!margin] == 0 || XLENGTH(s) > 0) {
912 if (rdim[!margin] > 0 && XLENGTH(s) == rdim[!margin])
913 nms[!margin] = Rf_getAttrib(s, R_NamesSymbol);
916 if (TAG(a) != R_NilValue) {
918 nms[margin] = Rf_coerceVector(TAG(a), STRSXP);
919 else if (level == 2) {
920 PROTECT(nms_ = Rf_allocVector(EXPRSXP, 1));
921 SET_VECTOR_ELT(nms_, 0, CAR(e));
922 nms[margin] = Rf_coerceVector(nms_, STRSXP);
924 }
else if (level == 1 &&
TYPEOF(CAR(e)) == SYMSXP)
925 nms[margin] = Rf_coerceVector(CAR(e), STRSXP);
927 if (rdimnames[!margin] && nms[!margin] != R_NilValue) {
928 SET_VECTOR_ELT(dimnames, !margin, nms[!margin]);
929 rdimnames[!margin] = 0;
930 if (!rdimnames[margin])
933 if (rdimnames[ margin] && nms[ margin] != R_NilValue)
934 for (i = 0; i < r; ++i)
935 SET_STRING_ELT(marnames, pos + i,
936 STRING_ELT(nms[margin], i));
948 SEXP level, margin, exprs;
949 args = CDR(args); level = CAR(args);
950 args = CDR(args); margin = CAR(args);
951 args = CDR(args); exprs = CAR(args);
952 return bind(CDR(args), CDR(exprs), Rf_asInteger(margin), Rf_asInteger(level));
SEXP dense_as_kind(SEXP, const char *, char, int)
SEXP sparse_as_general(SEXP, const char *)
SEXP sparse_as_kind(SEXP, const char *, char)
SEXP sparse_as_Csparse(SEXP, const char *)
#define SWITCH5(c, template)
#define Matrix_Calloc(p, n, t)
const char * Matrix_class(SEXP, const char **, int, const char *)
SEXPTYPE kindToType(char)
#define Matrix_Free(p, n)
#define ERROR_INVALID_TYPE(_X_, _FUNC_)
#define GET_SLOT(x, name)
SEXP newObject(const char *)
const char * valid_matrix[]
#define SET_SLOT(x, name, value)
SEXP sparse_aggregate(SEXP from, const char *class)
static void coerceArgs(SEXP args, int margin, int *rdim, char kind, char repr)
static void bindArgs(SEXP args, int margin, SEXP ans, int *rdim, char kind, char repr)
static void scanArgs(SEXP args, SEXP exprs, int margin, int level, int *rdim, int *rdimnames, char *kind, char *repr)
static SEXP bind(SEXP args, SEXP exprs, int margin, int level)
SEXP sparse_as_Rsparse(SEXP from, const char *class)
SEXP diagonal_as_dense(SEXP from, const char *class, char kind, char shape, int packed, char ul, char ct)
SEXP sparse_as_dense(SEXP from, const char *class, int packed)
SEXP dense_as_sparse(SEXP from, const char *class, char repr)
SEXP index_as_dense(SEXP from, const char *class, char kind)
SEXP diagonal_as_sparse(SEXP from, const char *class, char kind, char shape, char repr, char ul, char ct)
SEXP matrix_as_sparse(SEXP from, const char *zzz, char ul, char ct, char nu, int mg)
SEXP sparse_as_Tsparse(SEXP from, const char *class)
SEXP dense_as_general(SEXP from, const char *class, int new)
SEXP index_as_sparse(SEXP from, const char *class, char kind, char repr)