Matrix r4655
Loading...
Searching...
No Matches
coerce.c
Go to the documentation of this file.
1#include <math.h> /* trunc */
2#include "Mdefines.h"
3#include "idz.h"
4#include "coerce.h"
5
6SEXP vector_as_dense(SEXP from, const char *zzz, char ul, char di,
7 int m, int n, int byrow, SEXP dimnames)
8{
9 SEXPTYPE tf = TYPEOF(from);
10 char cl[] = "...Matrix";
11 cl[0] = (zzz[0] == '.')
12 ? typeToKind(tf)
13 : ((zzz[0] == ',') ? ((tf == CPLXSXP) ? 'z' : 'd') : zzz[0]);
14 cl[1] = zzz[1];
15 cl[2] = zzz[2];
16#ifndef MATRIX_ENABLE_IMATRIX
17 if (cl[0] == 'i')
18 cl[0] = 'd';
19#endif
20 SEXPTYPE tt = kindToType(cl[0]);
21 PROTECT(from = coerceVector(from, tt));
22
23 if (cl[1] != 'g' && m != n)
24 error(_("attempt to construct non-square %s"),
25 (cl[1] == 's') ? "symmetricMatrix" : "triangularMatrix");
26
27 Matrix_int_fast64_t mn = (Matrix_int_fast64_t) m * n;
28 if (((cl[2] != 'p') ? mn : (mn + n) / 2) > R_XLEN_T_MAX)
29 error(_("attempt to allocate vector of length exceeding %s"),
30 "R_XLEN_T_MAX");
31
32 SEXP to = PROTECT(newObject(cl));
33
34 SEXP dim = GET_SLOT(to, Matrix_DimSym);
35 int *pdim = INTEGER(dim);
36 pdim[0] = m;
37 pdim[1] = n;
38
39 if (cl[1] != 's')
40 SET_SLOT(to, Matrix_DimNamesSym, dimnames);
41 else
42 set_symmetrized_DimNames(to, dimnames, -1);
43
44 if (cl[1] != 'g' && ul != 'U') {
45 SEXP uplo = PROTECT(mkString("L"));
46 SET_SLOT(to, Matrix_uploSym, uplo);
47 UNPROTECT(1);
48 }
49
50 if (cl[1] == 't' && di != 'N') {
51 SEXP diag = PROTECT(mkString("U"));
52 SET_SLOT(to, Matrix_diagSym, diag);
53 UNPROTECT(1);
54 }
55
56 /* FIXME: add argument 'new' and conditionally avoid allocation */
57 SEXP x = PROTECT(allocVector(tt, (cl[2] != 'p') ? mn : (mn + n) / 2));
58 R_xlen_t k, r = XLENGTH(from);
59 int i, j, recycle = r < mn;
60
61#define VAD_SUBCASES(_PREFIX_, _CTYPE_, _PTR_, _NA_) \
62 do { \
63 _CTYPE_ *dest = _PTR_(x), *src = _PTR_(from); \
64 if (r == 0) { \
65 while (mn-- > 0) \
66 *(dest++) = _NA_; \
67 } else if (r == 1) { \
68 while (mn-- > 0) \
69 *(dest++) = *src; \
70 } else if (cl[2] != 'p') { \
71 if (!recycle) { \
72 if (!byrow) \
73 Matrix_memcpy(dest, src, mn, sizeof(_CTYPE_)); \
74 else \
75 _PREFIX_ ## transpose2(dest, src, n, m); \
76 } else { \
77 if (!byrow) { \
78 k = 0; \
79 while (mn-- > 0) { \
80 if (k == r) k = 0; \
81 *(dest++) = src[k++]; \
82 } \
83 } else { \
84 for (j = 0; j < n; ++j) { \
85 k = j; \
86 for (i = 0; i < m; ++i) { \
87 k %= r; \
88 *(dest++) = src[k]; \
89 k += n; \
90 } \
91 } \
92 } \
93 } \
94 } else if (ul == 'U') { \
95 if (!byrow) { \
96 k = 0; \
97 for (j = 0; j < n; ++j) { \
98 for (i = 0; i <= j; ++i) { \
99 if (recycle) k %= r; \
100 *(dest++) = src[k++]; \
101 } \
102 k += n - j - 1; \
103 } \
104 } else { \
105 for (j = 0; j < n; ++j) { \
106 k = j; \
107 for (i = 0; i <= j; ++i) { \
108 if (recycle) k %= r; \
109 *(dest++) = src[k]; \
110 k += n; \
111 } \
112 } \
113 } \
114 } else { \
115 if (!byrow) { \
116 k = 0; \
117 for (j = 0; j < n; ++j) { \
118 for (i = j; i < n; ++i) { \
119 if (recycle) k %= r; \
120 *(dest++) = src[k++]; \
121 } \
122 k += j + 1; \
123 } \
124 } else { \
125 R_xlen_t d = 0; \
126 for (j = 0; j < n; ++j) { \
127 k = j + d; \
128 for (i = 0; i <= j; ++i) { \
129 if (recycle) k %= r; \
130 *(dest++) = src[k]; \
131 k += n; \
132 } \
133 d += n; \
134 } \
135 } \
136 } \
137 } while (0)
138
139 switch (tt) {
140 case LGLSXP:
141 VAD_SUBCASES(i, int, LOGICAL, NA_LOGICAL);
142 break;
143 case INTSXP:
144 VAD_SUBCASES(i, int, INTEGER, NA_INTEGER);
145 break;
146 case REALSXP:
147 VAD_SUBCASES(d, double, REAL, NA_REAL);
148 break;
149 case CPLXSXP:
150 VAD_SUBCASES(z, Rcomplex, COMPLEX, Matrix_zna);
151 break;
152 default:
153 break;
154 }
155
156#undef VAD_SUBCASES
157
158 SET_SLOT(to, Matrix_xSym, x);
159
160 UNPROTECT(3);
161 return to;
162}
163
164SEXP R_vector_as_dense(SEXP from, SEXP zzz, SEXP uplo, SEXP diag,
165 SEXP m, SEXP n, SEXP byrow, SEXP dimnames)
166{
167 switch (TYPEOF(from)) {
168 case LGLSXP:
169 case INTSXP:
170 case REALSXP:
171 case CPLXSXP:
172 break;
173 default:
174 ERROR_INVALID_TYPE(from, __func__);
175 break;
176 }
177
178 const char *zzz_;
179 if (TYPEOF(zzz) != STRSXP || LENGTH(zzz) < 1 ||
180 (zzz = STRING_ELT(zzz, 0)) == NA_STRING ||
181 (zzz_ = CHAR(zzz))[0] == '\0' ||
182 (zzz_ )[1] == '\0' ||
183 !((zzz_[1] == 'g' && (zzz_[2] == 'e' )) ||
184 (zzz_[1] == 's' && (zzz_[2] == 'y' || zzz_[2] == 'p')) ||
185 (zzz_[1] == 't' && (zzz_[2] == 'r' || zzz_[2] == 'p'))))
186 error(_("second argument of '%s' does not specify a subclass of %s"),
187 __func__, "denseMatrix");
188
189 char ul = 'U', di = 'N';
190 if (zzz_[1] != 'g') {
191 if (TYPEOF(uplo) != STRSXP || LENGTH(uplo) < 1 ||
192 (uplo = STRING_ELT(uplo, 0)) == NA_STRING ||
193 ((ul = *CHAR(uplo)) != 'U' && ul != 'L'))
194 error(_("'%s' must be \"%s\" or \"%s\""), "uplo", "U", "L");
195 }
196 if (zzz_[1] == 't') {
197 if (TYPEOF(diag) != STRSXP || LENGTH(diag) < 1 ||
198 (diag = STRING_ELT(diag, 0)) == NA_STRING ||
199 ((di = *CHAR(diag)) != 'N' && di != 'U'))
200 error(_("'%s' must be \"%s\" or \"%s\""), "diag", "N", "U");
201 }
202
203 int m_ = -1;
204 if (m != R_NilValue) {
205 if (TYPEOF(m) == INTSXP) {
206 int tmp;
207 if (LENGTH(m) >= 1 && (tmp = INTEGER(m)[0]) != NA_INTEGER &&
208 tmp >= 0)
209 m_ = tmp;
210 } else if (TYPEOF(m) == REALSXP) {
211 double tmp;
212 if (LENGTH(m) >= 1 && !ISNAN(tmp = REAL(m)[0]) &&
213 tmp >= 0.0) {
214 if (trunc(tmp) > INT_MAX)
215 error(_("dimensions cannot exceed %s"), "2^31-1");
216 m_ = (int) tmp;
217 }
218 }
219 if (m_ < 0)
220 error(_("invalid '%s' to '%s'"), "m", __func__);
221 }
222
223 int n_ = -1;
224 if (n != R_NilValue) {
225 if (TYPEOF(n) == INTSXP) {
226 int tmp;
227 if (LENGTH(n) >= 1 && (tmp = INTEGER(n)[0]) != NA_INTEGER &&
228 tmp >= 0)
229 n_ = tmp;
230 } else if (TYPEOF(n) == REALSXP) {
231 double tmp;
232 if (LENGTH(n) >= 1 && !ISNAN(tmp = REAL(n)[0]) &&
233 tmp >= 0.0) {
234 if (trunc(tmp) > INT_MAX)
235 error(_("dimensions cannot exceed %s"), "2^31-1");
236 n_ = (int) tmp;
237 }
238 }
239 if (n_ < 0)
240 error(_("invalid '%s' to '%s'"), "n", __func__);
241 }
242
243 int byrow_;
244 if (TYPEOF(byrow) != LGLSXP || LENGTH(byrow) < 1 ||
245 (byrow_ = LOGICAL(byrow)[0]) == NA_LOGICAL)
246 error(_("'%s' must be %s or %s"), "byrow", "TRUE", "FALSE");
247
248 if (dimnames != R_NilValue)
249 if (TYPEOF(dimnames) != VECSXP || LENGTH(dimnames) != 2)
250 error(_("invalid '%s' to '%s'"), "dimnames", __func__);
251
252 R_xlen_t vlen_ = XLENGTH(from);
253 if (zzz_[1] != 'g' && (m_ < 0) != (n_ < 0)) {
254 if (m_ < 0)
255 m_ = n_;
256 else
257 n_ = m_;
258 } else if (m_ < 0 && n_ < 0) {
259 if (vlen_ > INT_MAX)
260 error(_("dimensions cannot exceed %s"), "2^31-1");
261 m_ = (int) vlen_;
262 n_ = 1;
263 } else if (m_ < 0) {
264 if (vlen_ > (Matrix_int_fast64_t) INT_MAX * n_) {
265 if (n_ == 0)
266 error(_("nonempty vector supplied for empty matrix"));
267 else
268 error(_("dimensions cannot exceed %s"), "2^31-1");
269 }
270 m_ = (n_ == 0) ? 0 : vlen_ / n_ + (vlen_ % n_ != 0);
271 } else if (n_ < 0) {
272 if (vlen_ > (Matrix_int_fast64_t) m_ * INT_MAX) {
273 if (m_ == 0)
274 error(_("nonempty vector supplied for empty matrix"));
275 else
276 error(_("dimensions cannot exceed %s"), "2^31-1");
277 }
278 n_ = (m_ == 0) ? 0 : vlen_ / m_ + (vlen_ % m_ != 0);
279 }
280
281 Matrix_int_fast64_t mlen_ = (Matrix_int_fast64_t) m_ * n_;
282 if (vlen_ <= 1)
283 /* do nothing */ ;
284 else if (mlen_ == 0)
285 warning(_("nonempty vector supplied for empty matrix"));
286 else if (vlen_ > mlen_)
287 warning(_("vector length (%lld) exceeds matrix length (%d * %d)"),
288 (long long) vlen_, m_, n_);
289 else if (mlen_ % vlen_ != 0)
290 warning(_("matrix length (%d * %d) is not a multiple of vector length (%lld)"),
291 m_, n_, (long long) vlen_);
292
293 return
294 vector_as_dense(from, zzz_, ul, di, m_, n_, byrow_, dimnames);
295}
296
297SEXP matrix_as_dense(SEXP from, const char *zzz, char ul, char di,
298 int trans, int new)
299{
300 SEXPTYPE tf = TYPEOF(from);
301 char cl[] = "...Matrix";
302 cl[0] = (zzz[0] == '.')
303 ? typeToKind(tf)
304 : ((zzz[0] == ',') ? ((tf == CPLXSXP) ? 'z' : 'd') : zzz[0]);
305 cl[1] = zzz[1];
306 cl[2] = zzz[2];
307#ifndef MATRIX_ENABLE_IMATRIX
308 if (cl[0] == 'i')
309 cl[0] = 'd';
310#endif
311 SEXPTYPE tt = kindToType(cl[0]);
312 PROTECT(from = coerceVector(from, tt));
313
314 SEXP to = PROTECT(newObject(cl));
315 int nprotect = 2;
316
317 SEXP dim = getAttrib(from, R_DimSymbol), dimnames;
318 int *pdim, isM, m, n, doDN;
319 R_xlen_t mn = XLENGTH(from);
320
321 isM = TYPEOF(dim) == INTSXP && LENGTH(dim) == 2;
322 if (isM) {
323
324 pdim = INTEGER(dim);
325 m = pdim[0];
326 n = pdim[1];
327 if (m != n || n > 0) {
328 dim = GET_SLOT(to, Matrix_DimSym);
329 pdim = INTEGER(dim);
330 pdim[0] = m;
331 pdim[1] = n;
332 }
333
334 PROTECT(dimnames = getAttrib(from, R_DimNamesSymbol));
335 ++nprotect;
336 doDN = dimnames != R_NilValue;
337
338 } else {
339
340 if (mn > INT_MAX)
341 error(_("dimensions cannot exceed %s"), "2^31-1");
342 dim = GET_SLOT(to, Matrix_DimSym);
343 pdim = INTEGER(dim);
344 if (trans) {
345 pdim[0] = m = 1;
346 pdim[1] = n = (int) mn;
347 } else {
348 pdim[0] = m = (int) mn;
349 pdim[1] = n = 1;
350 }
351
352 SEXP nms = PROTECT(getAttrib(from, R_NamesSymbol));
353 ++nprotect;
354 doDN = nms != R_NilValue;
355 if (doDN) {
356 PROTECT(dimnames = allocVector(VECSXP, 2));
357 ++nprotect;
358 SET_VECTOR_ELT(dimnames, trans ? 1 : 0, nms);
359 }
360
361 }
362
363 if (cl[1] != 'g' && m != n)
364 error(_("attempt to construct non-square %s"),
365 (cl[1] == 's') ? "symmetricMatrix" : "triangularMatrix");
366
367 if (doDN) {
368 if (cl[1] != 's')
369 SET_SLOT(to, Matrix_DimNamesSym, dimnames);
370 else
371 set_symmetrized_DimNames(to, dimnames, -1);
372 }
373
374 if (cl[1] != 'g' && ul != 'U') {
375 SEXP uplo = PROTECT(mkString("L"));
376 SET_SLOT(to, Matrix_uploSym, uplo);
377 UNPROTECT(1); /* uplo */
378 }
379
380 if (cl[1] == 't' && di != 'N') {
381 SEXP diag = PROTECT(mkString("U"));
382 SET_SLOT(to, Matrix_diagSym, diag);
383 UNPROTECT(1); /* diag */
384 }
385
386 SEXP x;
387
388 if (cl[2] != 'p') {
389
390 if (new <= 0 || (new <= 1 && ATTRIB(from) == R_NilValue) ||
391 !MAYBE_REFERENCED(from)) {
392
393 if (ATTRIB(from) != R_NilValue && new >= 1) {
394 /* 'from' has attributes and no references : */
395 SET_ATTRIB(from, R_NilValue);
396 if (OBJECT(from))
397 SET_OBJECT(from, 0);
398 }
399 x = from;
400
401 } else {
402
403 PROTECT(x = allocVector(tt, mn));
404 ++nprotect;
405 switch (tt) {
406 case LGLSXP:
407 Matrix_memcpy(LOGICAL(x), LOGICAL(from), mn, sizeof(int));
408 break;
409 case INTSXP:
410 Matrix_memcpy(INTEGER(x), INTEGER(from), mn, sizeof(int));
411 break;
412 case REALSXP:
413 Matrix_memcpy(REAL(x), REAL(from), mn, sizeof(double));
414 break;
415 case CPLXSXP:
416 Matrix_memcpy(COMPLEX(x), COMPLEX(from), mn, sizeof(Rcomplex));
417 break;
418 default:
419 break;
420 }
421
422 }
423
424 } else {
425
426 PROTECT(x = allocVector(tt, (mn - n) / 2 + n));
427 ++nprotect;
428 switch (tt) {
429 case LGLSXP:
430 ipack2(LOGICAL(x), LOGICAL(from), n, ul, di);
431 break;
432 case INTSXP:
433 ipack2(INTEGER(x), INTEGER(from), n, ul, di);
434 break;
435 case REALSXP:
436 dpack2(REAL(x), REAL(from), n, ul, di);
437 break;
438 case CPLXSXP:
439 zpack2(COMPLEX(x), COMPLEX(from), n, ul, di);
440 break;
441 default:
442 break;
443 }
444
445 }
446
447 SET_SLOT(to, Matrix_xSym, x);
448
449 UNPROTECT(nprotect);
450 return to;
451}
452
453/* as(<matrix>, ".(ge|sy|sp|tr|tp)Matrix") */
454SEXP R_matrix_as_dense(SEXP from, SEXP zzz, SEXP uplo, SEXP diag,
455 SEXP trans)
456{
457 switch (TYPEOF(from)) {
458 case LGLSXP:
459 case INTSXP:
460 case REALSXP:
461 case CPLXSXP:
462 break;
463 default:
464 ERROR_INVALID_TYPE(from, __func__);
465 break;
466 }
467
468 const char *zzz_;
469 if (TYPEOF(zzz) != STRSXP || LENGTH(zzz) < 1 ||
470 (zzz = STRING_ELT(zzz, 0)) == NA_STRING ||
471 (zzz_ = CHAR(zzz))[0] == '\0' ||
472 (zzz_ )[1] == '\0' ||
473 !((zzz_[1] == 'g' && (zzz_[2] == 'e' )) ||
474 (zzz_[1] == 's' && (zzz_[2] == 'y' || zzz_[2] == 'p')) ||
475 (zzz_[1] == 't' && (zzz_[2] == 'r' || zzz_[2] == 'p'))))
476 error(_("second argument of '%s' does not specify a subclass of %s"),
477 __func__, "denseMatrix");
478
479 char ul = 'U', di = 'N';
480 if (zzz_[1] != 'g') {
481 if (TYPEOF(uplo) != STRSXP || LENGTH(uplo) < 1 ||
482 (uplo = STRING_ELT(uplo, 0)) == NA_STRING ||
483 ((ul = *CHAR(uplo)) != 'U' && ul != 'L'))
484 error(_("'%s' must be \"%s\" or \"%s\""), "uplo", "U", "L");
485 }
486 if (zzz_[1] == 't') {
487 if (TYPEOF(diag) != STRSXP || LENGTH(diag) < 1 ||
488 (diag = STRING_ELT(diag, 0)) == NA_STRING ||
489 ((di = *CHAR(diag)) != 'N' && di != 'U'))
490 error(_("'%s' must be \"%s\" or \"%s\""), "diag", "N", "U");
491 }
492
493 int trans_;
494 if (TYPEOF(trans) != LGLSXP || LENGTH(trans) < 1 ||
495 (trans_ = LOGICAL(trans)[0]) == NA_LOGICAL)
496 error(_("'%s' must be %s or %s"), "trans", "TRUE", "FALSE");
497
498 return matrix_as_dense(from, zzz_, ul, di, trans_, 1);
499}
500
501SEXP sparse_as_dense(SEXP from, const char *class, int packed)
502{
503 packed = packed && class[1] != 'g';
504
505 char cl[] = "...Matrix";
506 cl[0] = class[0];
507 cl[1] = class[1];
508 cl[2] = (packed) ? 'p' :
509 ((class[1] == 'g') ? 'e' : ((class[1] == 's') ? 'y' : 'r'));
510 SEXP to = PROTECT(newObject(cl));
511
512 SEXP dim = PROTECT(GET_SLOT(from, Matrix_DimSym));
513 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1];
514 Matrix_int_fast64_t len = (Matrix_int_fast64_t) m * n;
515 if (packed)
516 len = (len + n) / 2;
517 if (len > R_XLEN_T_MAX)
518 error(_("attempt to allocate vector of length exceeding %s"),
519 "R_XLEN_T_MAX");
520 if (class[2] != 'C' && packed && len > R_XLEN_T_MAX)
521 error(_("coercing n-by-n %s to %s is not supported for n*n exceeding %s"),
522 "[RT]sparseMatrix", "packedMatrix", "R_XLEN_T_MAX");
523 double bytes = (double) len * kindToSize(cl[0]);
524 if (bytes > 0x1.0p+30 /* 1 GiB */)
525 warning(_("sparse->dense coercion: allocating vector of size %0.1f GiB"),
526 0x1.0p-30 * bytes);
527 if (m != n || n > 0)
528 SET_SLOT(to, Matrix_DimSym, dim);
529 UNPROTECT(1); /* dim */
530
531 SEXP dimnames = PROTECT(GET_SLOT(from, Matrix_DimNamesSym));
532 SET_SLOT(to, Matrix_DimNamesSym, dimnames);
533 UNPROTECT(1); /* dimnames */
534
535 char ul = 'U', di = 'N';
536 if (class[1] != 'g') {
537 SEXP uplo = PROTECT(GET_SLOT(from, Matrix_uploSym));
538 ul = *CHAR(STRING_ELT(uplo, 0));
539 if (ul != 'U')
540 SET_SLOT(to, Matrix_uploSym, uplo);
541 UNPROTECT(1); /* uplo */
542
543 if (class[1] == 't') {
544 SEXP diag = PROTECT(GET_SLOT(from, Matrix_diagSym));
545 di = *CHAR(STRING_ELT(diag, 0));
546 if (di != 'N')
547 SET_SLOT(to, Matrix_diagSym, diag);
548 UNPROTECT(1); /* diag */
549 }
550 }
551
552 /* It remains to fill 'x' ... */
553
554 SEXP x1 = PROTECT(allocVector(kindToType(class[0]), (R_xlen_t) len)),
555 p0, i0, j0;
556 int *pp, *pi, *pj, nprotect = 2;
557 p0 = i0 = j0 = NULL;
558 pp = pi = pj = NULL;
559 SET_SLOT(to, Matrix_xSym, x1);
560
561 if (class[2] != 'T') {
562 PROTECT(p0 = GET_SLOT(from, Matrix_pSym));
563 ++nprotect;
564 pp = INTEGER(p0) + 1;
565 }
566 if (class[2] != 'R') {
567 PROTECT(i0 = GET_SLOT(from, Matrix_iSym));
568 ++nprotect;
569 pi = INTEGER(i0);
570 }
571 if (class[2] != 'C') {
572 PROTECT(j0 = GET_SLOT(from, Matrix_jSym));
573 ++nprotect;
574 pj = INTEGER(j0);
575 }
576
577#define SAD_CASES \
578 do { \
579 switch (class[0]) { \
580 case 'l': \
581 SAD_SUBCASES(int, LOGICAL, SHOW, FIRSTOF, INCREMENT_LOGICAL, 1); \
582 break; \
583 case 'i': \
584 SAD_SUBCASES(int, INTEGER, SHOW, FIRSTOF, INCREMENT_INTEGER, 1); \
585 break; \
586 case 'd': \
587 SAD_SUBCASES(double, REAL, SHOW, FIRSTOF, INCREMENT_REAL, 1.0); \
588 break; \
589 case 'z': \
590 SAD_SUBCASES(Rcomplex, COMPLEX, SHOW, FIRSTOF, INCREMENT_COMPLEX, Matrix_zone); \
591 break; \
592 default: \
593 break; \
594 } \
595 } while (0)
596
597#define SAD_SUBCASES(_CTYPE_, _PTR_, _MASK_, _REPLACE_, _INCREMENT_, _ONE_) \
598 do { \
599 _MASK_(_CTYPE_ *px0 = _PTR_(x0)); \
600 _CTYPE_ *px1 = _PTR_(x1) ; \
601 Matrix_memset(px1, 0, len, sizeof(_CTYPE_)); \
602 if (!packed) { \
603 /* .(ge|sy|tr)Matrix */ \
604 SAD_SUBSUBCASES(SAD_LOOP_C2NP, SAD_LOOP_R2NP, SAD_LOOP_T2NP, \
605 _MASK_, _REPLACE_, _INCREMENT_); \
606 if (class[1] == 't' && di != 'N') { \
607 px1 = _PTR_(x1); \
608 R_xlen_t n1a = (R_xlen_t) n + 1; \
609 for (int d = 0; d < n; ++d, px1 += n1a) \
610 *px1 = _ONE_; \
611 } \
612 } else if (ul == 'U') { \
613 /* upper triangular .(sp|tp)Matrix */ \
614 SAD_SUBSUBCASES(SAD_LOOP_C2UP, SAD_LOOP_R2UP, SAD_LOOP_T2UP, \
615 _MASK_, _REPLACE_, _INCREMENT_); \
616 if (class[1] == 't' && di != 'N') { \
617 px1 = _PTR_(x1); \
618 for (int d = 0; d < n; px1 += (++d) + 1) \
619 *px1 = _ONE_; \
620 } \
621 } else { \
622 /* lower triangular .(sp|tp)Matrix */ \
623 SAD_SUBSUBCASES(SAD_LOOP_C2LP, SAD_LOOP_R2LP, SAD_LOOP_T2LP, \
624 _MASK_, _REPLACE_, _INCREMENT_); \
625 if (class[1] == 't' && di != 'N') { \
626 px1 = _PTR_(x1); \
627 for (int d = 0; d < n; px1 += n - (d++)) \
628 *px1 = _ONE_; \
629 } \
630 } \
631 } while (0)
632
633#define SAD_SUBSUBCASES(_LOOP_C_, _LOOP_R_, _LOOP_T_, _MASK_, _REPLACE_, _INCREMENT_) \
634 do { \
635 switch (class[2]) { \
636 case 'C': \
637 { \
638 int j, k, kend; \
639 _LOOP_C_(_MASK_, _REPLACE_, _INCREMENT_); \
640 break; \
641 } \
642 case 'R': \
643 { \
644 int i, k, kend; \
645 _LOOP_R_(_MASK_, _REPLACE_, _INCREMENT_); \
646 break; \
647 } \
648 case 'T': \
649 { \
650 R_xlen_t index, k, nnz = XLENGTH(i0); \
651 _LOOP_T_(_MASK_, _REPLACE_, _INCREMENT_); \
652 break; \
653 } \
654 default: \
655 break; \
656 } \
657 } while (0)
658
659#define SAD_LOOP_C2NP(_MASK_, _REPLACE_, _INCREMENT_) \
660 do { \
661 for (j = 0, k = 0; j < n; ++j, px1 += m) { \
662 kend = pp[j]; \
663 while (k < kend) { \
664 px1[*pi] = _REPLACE_(*px0, 1); \
665 ++k; ++pi; _MASK_(++px0); \
666 } \
667 } \
668 } while (0)
669
670#define SAD_LOOP_C2UP(_MASK_, _REPLACE_, _INCREMENT_) \
671 do { \
672 for (j = 0, k = 0; j < n; px1 += (++j)) { \
673 kend = pp[j]; \
674 while (k < kend) { \
675 px1[*pi] = _REPLACE_(*px0, 1); \
676 ++k; ++pi; _MASK_(++px0); \
677 } \
678 } \
679 } while (0)
680
681#define SAD_LOOP_C2LP(_MASK_, _REPLACE_, _INCREMENT_) \
682 do { \
683 for (j = 0, k = 0; j < n; px1 += n - (j++)) { \
684 kend = pp[j]; \
685 while (k < kend) { \
686 px1[*pi - j] = _REPLACE_(*px0, 1); \
687 ++k; ++pi; _MASK_(++px0); \
688 } \
689 } \
690 } while (0)
691
692#define SAD_LOOP_R2NP(_MASK_, _REPLACE_, _INCREMENT_) \
693 do { \
694 R_xlen_t m1 = (R_xlen_t) m; \
695 for (i = 0, k = 0; i < m; ++i, ++px1) { \
696 kend = pp[i]; \
697 while (k < kend) { \
698 px1[m1 * *pj] = _REPLACE_(*px0, 1); \
699 ++k; ++pj; _MASK_(++px0); \
700 } \
701 } \
702 } while (0)
703
704#define SAD_LOOP_R2UP(_MASK_, _REPLACE_, _INCREMENT_) \
705 do { \
706 for (i = 0, k = 0; i < n; ++i) { \
707 kend = pp[i]; \
708 while (k < kend) { \
709 px1[PACKED_AR21_UP(i, *pj)] = _REPLACE_(*px0, 1); \
710 ++k; ++pj; _MASK_(++px0); \
711 } \
712 } \
713 } while (0)
714
715#define SAD_LOOP_R2LP(_MASK_, _REPLACE_, _INCREMENT_) \
716 do { \
717 R_xlen_t n2 = (R_xlen_t) n * 2; \
718 for (i = 0, k = 0; i < n; ++i) { \
719 kend = pp[i]; \
720 while (k < kend) { \
721 px1[PACKED_AR21_LO(i, *pj, n2)] = _REPLACE_(*px0, 1); \
722 ++k; ++pj; _MASK_(++px0); \
723 } \
724 } \
725 } while (0)
726
727#define SAD_LOOP_T2NP(_MASK_, _REPLACE_, _INCREMENT_) \
728 do { \
729 R_xlen_t m1 = (R_xlen_t) m; \
730 for (k = 0; k < nnz; ++k) { \
731 index = m1 * *pj + *pi; \
732 _INCREMENT_(px1[index], (*px0)); \
733 ++pi; ++pj; _MASK_(++px0); \
734 } \
735 } while (0)
736
737#define SAD_LOOP_T2UP(_MASK_, _REPLACE_, _INCREMENT_) \
738 do { \
739 for (k = 0; k < nnz; ++k) { \
740 index = PACKED_AR21_UP(*pi, *pj); \
741 _INCREMENT_(px1[index], (*px0)); \
742 ++pi; ++pj; _MASK_(++px0); \
743 } \
744 } while (0)
745
746#define SAD_LOOP_T2LP(_MASK_, _REPLACE_, _INCREMENT_) \
747 do { \
748 R_xlen_t n2 = (R_xlen_t) n * 2; \
749 for (k = 0; k < nnz; ++k) { \
750 index = PACKED_AR21_LO(*pi, *pj, n2); \
751 _INCREMENT_(px1[index], (*px0)); \
752 ++pi; ++pj; _MASK_(++px0); \
753 } \
754 } while (0)
755
756 if (class[0] == 'n')
757 SAD_SUBCASES(int, LOGICAL, HIDE, SECONDOF, INCREMENT_PATTERN, 1);
758 else {
759 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym));
760 SAD_CASES;
761 UNPROTECT(1); /* x0 */
762 }
763
764#undef SAD_CASES
765#undef SAD_SUBCASES
766#undef SAD_SUBSUBCASES
767#undef SAD_LOOP_C2NP
768#undef SAD_LOOP_C2UP
769#undef SAD_LOOP_C2LP
770#undef SAD_LOOP_R2NP
771#undef SAD_LOOP_R2UP
772#undef SAD_LOOP_R2LP
773#undef SAD_LOOP_T2NP
774#undef SAD_LOOP_T2UP
775#undef SAD_LOOP_T2LP
776
777 UNPROTECT(nprotect);
778 return to;
779}
780
781/* as(<[CRT]sparseMatrix>, "(un)?packedMatrix") */
782SEXP R_sparse_as_dense(SEXP from, SEXP packed)
783{
784 static const char *valid[] = {
785 VALID_CSPARSE, VALID_RSPARSE, VALID_TSPARSE, "" };
786 int ivalid = R_check_class_etc(from, valid);
787 if (ivalid < 0)
788 ERROR_INVALID_CLASS(from, __func__);
789
790 int packed_;
791 if (TYPEOF(packed) != LGLSXP || LENGTH(packed) < 1 ||
792 (packed_ = LOGICAL(packed)[0]) == NA_LOGICAL)
793 error(_("'%s' must be %s or %s"), "packed", "TRUE", "FALSE");
794
795 return sparse_as_dense(from, valid[ivalid], packed_);
796}
797
798SEXP diagonal_as_dense(SEXP from, const char *class,
799 char kind, char shape, int packed, char ul)
800{
801 char cl[] = "...Matrix";
802 cl[0] = (kind == '.') ? class[0] : ((kind == ',') ? ((class[0] == 'z') ? 'z' : 'd') : kind);
803 cl[1] = shape;
804 cl[2] = (cl[1] == 'g') ? 'e' : ((packed) ? 'p' : ((cl[1] == 's') ? 'y' : 'r'));
805 SEXP to = PROTECT(newObject(cl));
806
807 SEXP dim = PROTECT(GET_SLOT(from, Matrix_DimSym));
808 int n = INTEGER(dim)[0];
809 Matrix_int_fast64_t len = (Matrix_int_fast64_t) n * n;
810 if (len > R_XLEN_T_MAX)
811 error(_("attempt to allocate vector of length exceeding %s"),
812 "R_XLEN_T_MAX");
813 double bytes = (double) len * kindToSize(cl[0]);
814 if (bytes > 0x1.0p+30 /* 1 GiB */)
815 warning(_("sparse->dense coercion: allocating vector of size %0.1f GiB"),
816 0x1.0p-30 * bytes);
817 if (n > 0)
818 SET_SLOT(to, Matrix_DimSym, dim);
819 UNPROTECT(1); /* dim */
820
821 SEXP dimnames = PROTECT(GET_SLOT(from, Matrix_DimNamesSym));
822 if (cl[1] == 's')
823 set_symmetrized_DimNames(to, dimnames, -1);
824 else
825 SET_SLOT(to, Matrix_DimNamesSym, dimnames);
826 UNPROTECT(1); /* dimnames */
827
828 if (cl[1] != 'g' && ul != 'U') {
829 SEXP uplo = PROTECT(mkString("L"));
830 SET_SLOT(to, Matrix_uploSym, uplo);
831 UNPROTECT(1); /* uplo */
832 }
833
834 SEXP diag = PROTECT(GET_SLOT(from, Matrix_diagSym));
835 char di = *CHAR(STRING_ELT(diag, 0));
836 if (cl[1] == 't' && di != 'N')
837 SET_SLOT(to, Matrix_diagSym, diag);
838 UNPROTECT(1); /* diag */
839
840 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym));
841 if (class[0] != cl[0]) {
842 if (class[0] == 'n' && cl[0] == 'l')
843 x0 = duplicate(x0);
844 else
845 x0 = coerceVector(x0, kindToType(cl[0]));
846 if (class[0] == 'n')
847 naToOne(x0);
848 UNPROTECT(1); /* x0 */
849 PROTECT(x0);
850 }
851
852 SEXP x1 = PROTECT(allocVector(TYPEOF(x0), (R_xlen_t) len));
853 SET_SLOT(to, Matrix_xSym, x1);
854
855#define DAD_SUBCASES(_PREFIX_, _CTYPE_, _PTR_) \
856 do { \
857 _CTYPE_ *px0 = _PTR_(x0), *px1 = _PTR_(x1); \
858 Matrix_memset(px1, 0, (R_xlen_t) len, sizeof(_CTYPE_)); \
859 if (di == 'N' || cl[1] != 't') { \
860 if (cl[2] != 'p') \
861 _PREFIX_ ## dcpy2(px1, px0, n, n, ul, di); \
862 else \
863 _PREFIX_ ## dcpy1(px1, px0, n, n, ul, ul, di); \
864 } \
865 } while (0)
866
867 switch (cl[0]) {
868 case 'n':
869 case 'l':
870 DAD_SUBCASES(i, int, LOGICAL);
871 break;
872 case 'i':
873 DAD_SUBCASES(i, int, INTEGER);
874 break;
875 case 'd':
876 DAD_SUBCASES(d, double, REAL);
877 break;
878 case 'z':
879 DAD_SUBCASES(z, Rcomplex, COMPLEX);
880 break;
881 default:
882 break;
883 }
884
885#undef DAD_CASES
886#undef DAD_SUBCASES
887
888 UNPROTECT(3); /* x1, x0, to */
889 return to;
890}
891
892/* as(<diagonalMatrix>, ".(ge|sy|sp|tr|tp)Matrix") */
893SEXP R_diagonal_as_dense(SEXP from,
894 SEXP kind, SEXP shape, SEXP packed, SEXP uplo)
895{
896 static const char *valid[] = { VALID_DIAGONAL, "" };
897 int ivalid = R_check_class_etc(from, valid);
898 if (ivalid < 0)
899 ERROR_INVALID_CLASS(from, __func__);
900
901 char kind_;
902 if (TYPEOF(kind) != STRSXP || LENGTH(kind) < 1 ||
903 (kind = STRING_ELT(kind, 0)) == NA_STRING ||
904 (kind_ = CHAR(kind)[0]) == '\0')
905 error(_("invalid '%s' to '%s'"), "kind", __func__);
906
907 char shape_;
908 if (TYPEOF(shape) != STRSXP || LENGTH(shape) < 1 ||
909 (shape = STRING_ELT(shape, 0)) == NA_STRING ||
910 ((shape_ = CHAR(shape)[0]) != 'g' && shape_ != 's' && shape_ != 't'))
911 error(_("invalid '%s' to '%s'"), "shape", __func__);
912
913 int packed_ = 0;
914 if (shape_ != 'g') {
915 if (TYPEOF(packed) != LGLSXP || LENGTH(packed) < 1 ||
916 (packed_ = LOGICAL(packed)[0]) == NA_LOGICAL)
917 error(_("'%s' must be %s or %s"), "packed", "TRUE", "FALSE");
918 }
919
920 char ul = 'U';
921 if (shape_ != 'g') {
922 if (TYPEOF(uplo) != STRSXP || LENGTH(uplo) < 1 ||
923 (uplo = STRING_ELT(uplo, 0)) == NA_STRING ||
924 ((ul = *CHAR(uplo)) != 'U' && ul != 'L'))
925 error(_("'%s' must be \"%s\" or \"%s\""), "uplo", "U", "L");
926 }
927
928 return diagonal_as_dense(from, valid[ivalid], kind_, shape_, packed_, ul);
929}
930
931SEXP index_as_dense(SEXP from, const char *class, char kind)
932{
933 SEXP margin = PROTECT(GET_SLOT(from, Matrix_marginSym));
934 int mg = INTEGER(margin)[0] - 1;
935 UNPROTECT(1); /* margin */
936
937 char cl[] = ".geMatrix";
938 cl[0] = (kind == '.') ? 'n' : ((kind == ',') ? 'd' : kind);
939 SEXP to = PROTECT(newObject(cl));
940
941 SEXP dim = PROTECT(GET_SLOT(from, Matrix_DimSym));
942 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1];
943 Matrix_int_fast64_t len = (Matrix_int_fast64_t) m * n;
944 if (len > R_XLEN_T_MAX)
945 error(_("attempt to allocate vector of length exceeding %s"),
946 "R_XLEN_T_MAX");
947 double bytes = (double) len * kindToSize(cl[0]);
948 if (bytes > 0x1.0p+30 /* 1 GiB */)
949 warning(_("sparse->dense coercion: allocating vector of size %0.1f GiB"),
950 0x1.0p-30 * bytes);
951 if (m != n || n > 0)
952 SET_SLOT(to, Matrix_DimSym, dim);
953 UNPROTECT(1); /* dim */
954
955 SEXP dimnames = PROTECT(GET_SLOT(from, Matrix_DimNamesSym));
956 SET_SLOT(to, Matrix_DimNamesSym, dimnames);
957 UNPROTECT(1); /* dimnames */
958
959 SEXP perm = PROTECT(GET_SLOT(from, Matrix_permSym));
960 int *pperm = INTEGER(perm);
961
962 SEXP x = PROTECT(allocVector(kindToType(cl[0]), (R_xlen_t) len));
963 SET_SLOT(to, Matrix_xSym, x);
964
965#define IAD_SUBCASES(_CTYPE_, _PTR_, _ONE_) \
966 do { \
967 _CTYPE_ *px = _PTR_(x); \
968 Matrix_memset(px, 0, (R_xlen_t) len, sizeof(_CTYPE_)); \
969 if (mg == 0) { \
970 R_xlen_t m1 = (R_xlen_t) m; \
971 for (int i = 0; i < m; ++i) \
972 px[i + m1 * (*(pperm++) - 1)] = _ONE_; \
973 } else { \
974 for (int j = 0; j < n; ++j, px += m) \
975 px[ *(pperm++) - 1 ] = _ONE_; \
976 } \
977 } while (0)
978
979 switch (cl[0]) {
980 case 'n':
981 case 'l':
982 IAD_SUBCASES(int, LOGICAL, 1);
983 break;
984 case 'i':
985 IAD_SUBCASES(int, INTEGER, 1);
986 break;
987 case 'd':
988 IAD_SUBCASES(double, REAL, 1.0);
989 break;
990 case 'z':
991 IAD_SUBCASES(Rcomplex, COMPLEX, Matrix_zone);
992 break;
993 default:
994 break;
995 }
996
997#undef IAD_SUBCASES
998
999 UNPROTECT(3); /* x, perm, to */
1000 return to;
1001}
1002
1003/* as(<indMatrix>, ".geMatrix") */
1004SEXP R_index_as_dense(SEXP from, SEXP kind)
1005{
1006 static const char *valid[] = { "indMatrix", "pMatrix" };
1007 int ivalid = R_check_class_etc(from, valid);
1008 if (ivalid < 0)
1009 ERROR_INVALID_CLASS(from, __func__);
1010
1011 char kind_;
1012 if (TYPEOF(kind) != STRSXP || LENGTH(kind) < 1 ||
1013 (kind = STRING_ELT(kind, 0)) == NA_STRING ||
1014 (kind_ = CHAR(kind)[0]) == '\0')
1015 error(_("invalid '%s' to '%s'"), "kind", __func__);
1016
1017 return index_as_dense(from, valid[ivalid], kind_);
1018}
1019
1020SEXP vector_as_sparse(SEXP from, const char *zzz, char ul, char di,
1021 int m, int n, int byrow, SEXP dimnames)
1022{
1023 SEXP length0 = GET_SLOT(from, Matrix_lengthSym);
1024 Matrix_int_fast64_t r = (Matrix_int_fast64_t)
1025 ((TYPEOF(length0) == INTSXP) ? INTEGER(length0)[0] : REAL(length0)[0]);
1026
1027 SEXP i0 = PROTECT(GET_SLOT(from, Matrix_iSym)),
1028 x0 = getAttrib(from, Matrix_xSym);
1029
1030 SEXPTYPE tf = TYPEOF(x0);
1031 char cl[] = "...Matrix";
1032 cl[0] = (zzz[0] == '.')
1033 ? ((x0 == R_NilValue) ? 'n' : typeToKind(tf))
1034 : ((zzz[0] == ',') ? ((tf == CPLXSXP) ? 'z' : 'd') : zzz[0]);
1035 cl[1] = zzz[1];
1036 cl[2] = (byrow) ? 'R' : 'C';
1037#ifndef MATRIX_ENABLE_IMATRIX
1038 if (cl[0] == 'i')
1039 cl[0] = 'd';
1040#endif
1041 SEXPTYPE tt = kindToType(cl[0]);
1042 if (x0 != R_NilValue) {
1043 PROTECT(x0);
1044 x0 = coerceVector(x0, tt);
1045 UNPROTECT(1); /* x0 */
1046 }
1047 PROTECT(x0);
1048
1049 if (cl[1] != 'g' && m != n)
1050 error(_("attempt to construct non-square %s"),
1051 (cl[1] == 's') ? "symmetricMatrix" : "triangularMatrix");
1052
1053 SEXP to = PROTECT(newObject(cl));
1054
1055 SEXP dim = GET_SLOT(to, Matrix_DimSym);
1056 int *pdim = INTEGER(dim);
1057 pdim[0] = m;
1058 pdim[1] = n;
1059
1060 if (cl[1] != 's')
1061 SET_SLOT(to, Matrix_DimNamesSym, dimnames);
1062 else
1063 set_symmetrized_DimNames(to, dimnames, -1);
1064
1065 if (cl[1] != 'g' && ul != 'U') {
1066 SEXP uplo = PROTECT(mkString("L"));
1067 SET_SLOT(to, Matrix_uploSym, uplo);
1068 UNPROTECT(1); /* uplo */
1069 }
1070
1071 if (cl[1] == 't' && di != 'N') {
1072 SEXP diag = PROTECT(mkString("U"));
1073 SET_SLOT(to, Matrix_diagSym, diag);
1074 UNPROTECT(1); /* diag */
1075 }
1076
1077 Matrix_int_fast64_t pos, mn = (Matrix_int_fast64_t) m * n, nnz1 = 0;
1078 R_xlen_t k = 0, nnz0 = XLENGTH(i0);
1079
1080#define VAS_SUBCASES(...) \
1081 do { \
1082 switch (TYPEOF(i0)) { \
1083 case INTSXP: \
1084 { \
1085 int *pi0 = INTEGER(i0); \
1086 VAS_SUBSUBCASES(__VA_ARGS__); \
1087 break; \
1088 } \
1089 case REALSXP: \
1090 { \
1091 double *pi0 = REAL(i0); \
1092 VAS_SUBSUBCASES(__VA_ARGS__); \
1093 break; \
1094 } \
1095 default: \
1096 break; \
1097 } \
1098 } while (0)
1099
1100#define VAS_SUBSUBCASES() \
1101 do { \
1102 if (nnz0 == 0) \
1103 /* do nothing */ ; \
1104 else if (cl[1] == 'g') { \
1105 if (r == 0) \
1106 nnz1 = mn; \
1107 else if (r == mn) \
1108 nnz1 = nnz0; \
1109 else if (r > mn) \
1110 while (k < nnz0 && (Matrix_int_fast64_t) pi0[k++] <= mn) \
1111 nnz1++; \
1112 else { \
1113 Matrix_int_fast64_t mn_mod_r = mn % r; \
1114 nnz1 = nnz0 * (mn / r); \
1115 while (k < nnz0 && (Matrix_int_fast64_t) pi0[k++] <= mn_mod_r) \
1116 nnz1++; \
1117 } \
1118 } \
1119 else if (cl[1] == 's' || di == 'N') { \
1120 if (r == 0) \
1121 nnz1 = (mn + n) / 2; \
1122 else if (r >= mn) { \
1123 if ((ul == 'U') == !byrow) { \
1124 while (k < nnz0 && (pos = (Matrix_int_fast64_t) pi0[k++] - 1) < mn) \
1125 if (pos % n <= pos / n) \
1126 ++nnz1; \
1127 } else { \
1128 while (k < nnz0 && (pos = (Matrix_int_fast64_t) pi0[k++] - 1) < mn) \
1129 if (pos % n >= pos / n) \
1130 ++nnz1; \
1131 } \
1132 } \
1133 else { \
1134 Matrix_int_fast64_t a = 0; \
1135 if ((ul == 'U') == !byrow) { \
1136 while (a < mn) { \
1137 k = 0; \
1138 while (k < nnz0 && (pos = a + pi0[k++] - 1) < mn) \
1139 if (pos % n <= pos / n) \
1140 ++nnz1; \
1141 a += r; \
1142 } \
1143 } else { \
1144 while (a < mn) { \
1145 k = 0; \
1146 while (k < nnz0 && (pos = a + pi0[k++] - 1) < mn) \
1147 if (pos % n >= pos / n) \
1148 ++nnz1; \
1149 a += r; \
1150 } \
1151 } \
1152 } \
1153 } \
1154 else { \
1155 if (r == 0) \
1156 nnz1 = (mn - n) / 2; \
1157 else if (r >= mn) { \
1158 if ((ul == 'U') == !byrow) { \
1159 while (k < nnz0 && (pos = (Matrix_int_fast64_t) pi0[k++] - 1) < mn) \
1160 if (pos % n < pos / n) \
1161 ++nnz1; \
1162 } else { \
1163 while (k < nnz0 && (pos = (Matrix_int_fast64_t) pi0[k++] - 1) < mn) \
1164 if (pos % n > pos / n) \
1165 ++nnz1; \
1166 } \
1167 } \
1168 else { \
1169 Matrix_int_fast64_t a = 0; \
1170 if ((ul == 'U') == !byrow) { \
1171 while (a < mn) { \
1172 k = 0; \
1173 while (k < nnz0 && (pos = a + pi0[k++] - 1) < mn) \
1174 if (pos % n < pos / n) \
1175 ++nnz1; \
1176 a += r; \
1177 } \
1178 } else { \
1179 while (a < mn) { \
1180 k = 0; \
1181 while (k < nnz0 && (pos = a + pi0[k++] - 1) < mn) \
1182 if (pos % n > pos / n) \
1183 ++nnz1; \
1184 a += r; \
1185 } \
1186 } \
1187 } \
1188 } \
1189 } while (0)
1190
1191 VAS_SUBCASES();
1192
1193#undef VAS_SUBSUBCASES
1194
1195 if (nnz1 > INT_MAX)
1196 error(_("attempt to construct %s with more than %s nonzero entries"),
1197 "sparseMatrix", "2^31-1");
1198
1199 int i_, j_, m_ = (byrow) ? n : m, n_ = (byrow) ? m : n;
1200 SEXP iSym = (byrow) ? Matrix_jSym : Matrix_iSym,
1201 p1 = PROTECT(allocVector(INTSXP, (R_xlen_t) n_ + 1)),
1202 i1 = PROTECT(allocVector(INTSXP, nnz1));
1203 SET_SLOT(to, Matrix_pSym, p1);
1204 SET_SLOT(to, iSym, i1);
1205 int *pp1 = INTEGER(p1) + 1, *pi1 = INTEGER(i1);
1206 Matrix_memset(pp1 - 1, 0, (R_xlen_t) n + 1, sizeof(int));
1207 k = 0;
1208
1209#define VAS_SUBSUBCASES(_MASK0_, _MASK1_, _REPLACE_, _CTYPE_, _PTR_, _ONE_, _NA_) \
1210 do { \
1211 _MASK0_(_CTYPE_ *px0 = _PTR_(x0)); \
1212 _MASK1_(_CTYPE_ *px1 = _PTR_(x1)); \
1213 if (nnz1 == 0) \
1214 /* do nothing */ ; \
1215 else if (cl[1] == 'g') { \
1216 if (r == 0) { \
1217 for (j_ = 0; j_ < n_; ++j_) { \
1218 pp1[j_] = m; \
1219 for (i_ = 0; i_ < m_; ++i_) { \
1220 *(pi1++) = i_; \
1221 _MASK1_(*(px1++) = _NA_); \
1222 } \
1223 } \
1224 } \
1225 else if (r >= mn) { \
1226 while (k < nnz0 && (pos = (Matrix_int_fast64_t) pi0[k] - 1) < mn) { \
1227 ++pp1[pos / m_]; \
1228 *(pi1++) = pos % m_; \
1229 _MASK1_(*(px1++) = _REPLACE_(px0[k], _ONE_)); \
1230 ++k; \
1231 } \
1232 } \
1233 else { \
1234 Matrix_int_fast64_t a = 0; \
1235 while (a < mn) { \
1236 k = 0; \
1237 while (k < nnz0 && (pos = a + pi0[k] - 1) < mn) { \
1238 ++pp1[pos / m_]; \
1239 *(pi1++) = pos % m_; \
1240 _MASK1_(*(px1++) = _REPLACE_(px0[k], _ONE_)); \
1241 ++k; \
1242 } \
1243 a += r; \
1244 } \
1245 } \
1246 } \
1247 else if (cl[1] == 's' || di == 'N') { \
1248 if (r == 0) { \
1249 if ((ul == 'U') == !byrow) { \
1250 for (j_ = 0; j_ < n_; ++j_) { \
1251 pp1[j_] = j_ + 1; \
1252 for (i_ = 0; i_ <= j_; ++i_) { \
1253 *(pi1++) = i_; \
1254 _MASK1_(*(px1++) = _NA_); \
1255 } \
1256 } \
1257 } else { \
1258 for (j_ = 0; j_ < n_; ++j_) { \
1259 pp1[j_] = n_ - j_; \
1260 for (i_ = j_; i_ < n_; ++i_) { \
1261 *(pi1++) = i_; \
1262 _MASK1_(*(px1++) = _NA_); \
1263 } \
1264 } \
1265 } \
1266 } \
1267 else if (r >= mn) { \
1268 if ((ul == 'U') == !byrow) { \
1269 while (k < nnz0 && (pos = (Matrix_int_fast64_t) pi0[k] - 1) < mn) { \
1270 if ((i_ = pos % n_) <= (j_ = pos / n_)) { \
1271 ++pp1[j_]; \
1272 *(pi1++) = i_; \
1273 _MASK1_(*(px1++) = _REPLACE_(px0[k], _ONE_)); \
1274 } \
1275 ++k; \
1276 } \
1277 } else { \
1278 while (k < nnz0 && (pos = (Matrix_int_fast64_t) pi0[k] - 1) < mn) { \
1279 if ((i_ = pos % n_) >= (j_ = pos / n_)) { \
1280 ++pp1[j_]; \
1281 *(pi1++) = i_; \
1282 _MASK1_(*(px1++) = _REPLACE_(px0[k], _ONE_)); \
1283 } \
1284 ++k; \
1285 } \
1286 } \
1287 } \
1288 else { \
1289 Matrix_int_fast64_t a = 0; \
1290 if ((ul == 'U') == !byrow) { \
1291 while (a < mn) { \
1292 k = 0; \
1293 while (k < nnz0 && (pos = a + pi0[k] - 1) < mn) { \
1294 if ((i_ = pos % n) <= (j_ = pos / n)) { \
1295 ++pp1[j_]; \
1296 *(pi1++) = i_; \
1297 _MASK1_(*(px1++) = _REPLACE_(px0[k], _ONE_)); \
1298 } \
1299 ++k; \
1300 } \
1301 a += r; \
1302 } \
1303 } else { \
1304 while (a < mn) { \
1305 k = 0; \
1306 while (k < nnz0 && (pos = a + pi0[k] - 1) < mn) { \
1307 if ((i_ = pos % n) >= (j_ = pos / n)) { \
1308 ++pp1[j_]; \
1309 *(pi1++) = i_; \
1310 _MASK1_(*(px1++) = _REPLACE_(px0[k], _ONE_)); \
1311 } \
1312 ++k; \
1313 } \
1314 a += r; \
1315 } \
1316 } \
1317 } \
1318 } \
1319 else { \
1320 if (r == 0) { \
1321 if ((ul == 'U') == !byrow) { \
1322 for (j_ = 0; j_ < n_; ++j_) { \
1323 pp1[j_] = j_; \
1324 for (i_ = 0; i_ < j_; ++i_) { \
1325 *(pi1++) = i_; \
1326 _MASK1_(*(px1++) = _NA_); \
1327 } \
1328 } \
1329 } else { \
1330 for (j_ = 0; j_ < n_; ++j_) { \
1331 pp1[j_] = n_ - j_ - 1; \
1332 for (i_ = j_ + 1; i_ < n_; ++i_) { \
1333 *(pi1++) = i_; \
1334 _MASK1_(*(px1++) = _NA_); \
1335 } \
1336 } \
1337 } \
1338 } \
1339 else if (r >= mn) { \
1340 if ((ul == 'U') == !byrow) { \
1341 while (k < nnz0 && (pos = (Matrix_int_fast64_t) pi0[k] - 1) < mn) { \
1342 if ((i_ = pos % n_) < (j_ = pos / n_)) { \
1343 ++pp1[j_]; \
1344 *(pi1++) = i_; \
1345 _MASK1_(*(px1++) = _REPLACE_(px0[k], _ONE_)); \
1346 } \
1347 ++k; \
1348 } \
1349 } else { \
1350 while (k < nnz0 && (pos = (Matrix_int_fast64_t) pi0[k] - 1) < mn) { \
1351 if ((i_ = pos % n_) > (j_ = pos / n_)) { \
1352 ++pp1[j_]; \
1353 *(pi1++) = i_; \
1354 _MASK1_(*(px1++) = _REPLACE_(px0[k], _ONE_)); \
1355 } \
1356 ++k; \
1357 } \
1358 } \
1359 } \
1360 else { \
1361 Matrix_int_fast64_t a = 0; \
1362 if ((ul == 'U') == !byrow) { \
1363 while (a < mn) { \
1364 k = 0; \
1365 while (k < nnz0 && (pos = a + pi0[k] - 1) < mn) { \
1366 if ((i_ = pos % n) < (j_ = pos / n)) { \
1367 ++pp1[j_]; \
1368 *(pi1++) = i_; \
1369 _MASK1_(*(px1++) = _REPLACE_(px0[k], _ONE_)); \
1370 } \
1371 ++k; \
1372 } \
1373 a += r; \
1374 } \
1375 } else { \
1376 while (a < mn) { \
1377 k = 0; \
1378 while (k < nnz0 && (pos = a + pi0[k] - 1) < mn) { \
1379 if ((i_ = pos % n) > (j_ = pos / n)) { \
1380 ++pp1[j_]; \
1381 *(pi1++) = i_; \
1382 _MASK1_(*(px1++) = _REPLACE_(px0[k], _ONE_)); \
1383 } \
1384 ++k; \
1385 } \
1386 a += r; \
1387 } \
1388 } \
1389 } \
1390 } \
1391 } while (0)
1392
1393 if (cl[0] == 'n')
1394 VAS_SUBCASES(HIDE, HIDE, , , , , );
1395 else {
1396 SEXP x1 = PROTECT(allocVector(kindToType(cl[0]), nnz1));
1397 switch (cl[0]) {
1398 case 'l':
1399 if (x0 == R_NilValue)
1400 VAS_SUBCASES(HIDE, SHOW, SECONDOF, int, LOGICAL, 1, NA_LOGICAL);
1401 else
1402 VAS_SUBCASES(SHOW, SHOW, FIRSTOF, int, LOGICAL, 1, NA_LOGICAL);
1403 break;
1404 case 'i':
1405 if (x0 == R_NilValue)
1406 VAS_SUBCASES(HIDE, SHOW, SECONDOF, int, INTEGER, 1, NA_INTEGER);
1407 else
1408 VAS_SUBCASES(SHOW, SHOW, FIRSTOF, int, INTEGER, 1, NA_INTEGER);
1409 break;
1410 case 'd':
1411 if (x0 == R_NilValue)
1412 VAS_SUBCASES(HIDE, SHOW, SECONDOF, double, REAL, 1.0, NA_REAL);
1413 else
1414 VAS_SUBCASES(SHOW, SHOW, FIRSTOF, double, REAL, 1.0, NA_REAL);
1415 break;
1416 case 'z':
1417 if (x0 == R_NilValue)
1418 VAS_SUBCASES(HIDE, SHOW, SECONDOF, Rcomplex, COMPLEX, Matrix_zone, Matrix_zna);
1419 else
1420 VAS_SUBCASES(SHOW, SHOW, FIRSTOF, Rcomplex, COMPLEX, Matrix_zone, Matrix_zna);
1421 break;
1422 default:
1423 break;
1424 }
1425 SET_SLOT(to, Matrix_xSym, x1);
1426 UNPROTECT(1); /* x1 */
1427 }
1428
1429#undef VAS_SUBCASES
1430#undef VAS_SUBSUBCASES
1431
1432 for (j_ = 0; j_ < n_; ++j_)
1433 pp1[j_] += pp1[j_ - 1];
1434
1435 switch (zzz[2]) {
1436 case 'C':
1437 to = sparse_as_Csparse(to, cl);
1438 break;
1439 case 'R':
1440 to = sparse_as_Rsparse(to, cl);
1441 break;
1442 case 'T':
1443 to = sparse_as_Tsparse(to, cl);
1444 break;
1445 default:
1446 break;
1447 }
1448
1449 UNPROTECT(5); /* i1, p1, to, x0, i0 */
1450 return to;
1451}
1452
1453SEXP R_vector_as_sparse(SEXP from, SEXP zzz, SEXP uplo, SEXP diag,
1454 SEXP m, SEXP n, SEXP byrow, SEXP dimnames)
1455{
1456 static const char *valid[] = { VALID_NONVIRTUAL_VECTOR, "" };
1457 int ivalid = R_check_class_etc(from, valid);
1458 if (ivalid < 0)
1459 ERROR_INVALID_CLASS(from, __func__);
1460
1461 const char *zzz_;
1462 if (TYPEOF(zzz) != STRSXP || LENGTH(zzz) < 1 ||
1463 (zzz = STRING_ELT(zzz, 0)) == NA_STRING ||
1464 (zzz_ = CHAR(zzz))[0] == '\0' ||
1465 (zzz_[1] != 'g' && zzz_[1] != 's' && zzz_[1] != 't') ||
1466 (zzz_[2] != 'C' && zzz_[2] != 'R' && zzz_[2] != 'T'))
1467 error(_("second argument of '%s' does not specify a subclass of %s"),
1468 __func__, "[CRT]sparseMatrix");
1469
1470 char ul = 'U', di = 'N';
1471 if (zzz_[1] != 'g') {
1472 if (TYPEOF(uplo) != STRSXP || LENGTH(uplo) < 1 ||
1473 (uplo = STRING_ELT(uplo, 0)) == NA_STRING ||
1474 ((ul = *CHAR(uplo)) != 'U' && ul != 'L'))
1475 error(_("'%s' must be \"%s\" or \"%s\""), "uplo", "U", "L");
1476 }
1477 if (zzz_[1] == 't') {
1478 if (TYPEOF(diag) != STRSXP || LENGTH(diag) < 1 ||
1479 (diag = STRING_ELT(diag, 0)) == NA_STRING ||
1480 ((di = *CHAR(diag)) != 'N' && di != 'U'))
1481 error(_("'%s' must be \"%s\" or \"%s\""), "diag", "N", "U");
1482 }
1483
1484 int m_ = -1;
1485 if (m != R_NilValue) {
1486 if (TYPEOF(m) == INTSXP) {
1487 int tmp;
1488 if (LENGTH(m) >= 1 && (tmp = INTEGER(m)[0]) != NA_INTEGER &&
1489 tmp >= 0)
1490 m_ = tmp;
1491 } else if (TYPEOF(m) == REALSXP) {
1492 double tmp;
1493 if (LENGTH(m) >= 1 && !ISNAN(tmp = REAL(m)[0]) &&
1494 tmp >= 0.0) {
1495 if (trunc(tmp) > INT_MAX)
1496 error(_("dimensions cannot exceed %s"), "2^31-1");
1497 m_ = (int) tmp;
1498 }
1499 }
1500 if (m_ < 0)
1501 error(_("invalid '%s' to '%s'"), "m", __func__);
1502 }
1503
1504 int n_ = -1;
1505 if (n != R_NilValue) {
1506 if (TYPEOF(n) == INTSXP) {
1507 int tmp;
1508 if (LENGTH(n) >= 1 && (tmp = INTEGER(n)[0]) != NA_INTEGER &&
1509 tmp >= 0)
1510 n_ = tmp;
1511 } else if (TYPEOF(n) == REALSXP) {
1512 double tmp;
1513 if (LENGTH(n) >= 1 && !ISNAN(tmp = REAL(n)[0]) &&
1514 tmp >= 0.0) {
1515 if (trunc(tmp) > INT_MAX)
1516 error(_("dimensions cannot exceed %s"), "2^31-1");
1517 n_ = (int) tmp;
1518 }
1519 }
1520 if (n_ < 0)
1521 error(_("invalid '%s' to '%s'"), "n", __func__);
1522 }
1523
1524 int byrow_;
1525 if (TYPEOF(byrow) != LGLSXP || LENGTH(byrow) < 1 ||
1526 (byrow_ = LOGICAL(byrow)[0]) == NA_LOGICAL)
1527 error(_("'%s' must be %s or %s"), "byrow", "TRUE", "FALSE");
1528
1529 if (dimnames != R_NilValue)
1530 if (TYPEOF(dimnames) != VECSXP || LENGTH(dimnames) != 2)
1531 error(_("invalid '%s' to '%s'"), "dimnames", __func__);
1532
1533 SEXP tmp = GET_SLOT(from, Matrix_lengthSym);
1534 Matrix_int_fast64_t vlen_ = (Matrix_int_fast64_t)
1535 ((TYPEOF(tmp) == INTSXP) ? INTEGER(tmp)[0] : REAL(tmp)[0]);
1536 if (zzz_[1] != 'g' && (m_ < 0) != (n_ < 0)) {
1537 if (m_ < 0)
1538 m_ = n_;
1539 else
1540 n_ = m_;
1541 } else if (m_ < 0 && n_ < 0) {
1542 if (vlen_ > INT_MAX)
1543 error(_("dimensions cannot exceed %s"), "2^31-1");
1544 m_ = (int) vlen_;
1545 n_ = 1;
1546 } else if (m_ < 0) {
1547 if (vlen_ > (Matrix_int_fast64_t) INT_MAX * n_) {
1548 if (n_ == 0)
1549 error(_("nonempty vector supplied for empty matrix"));
1550 else
1551 error(_("dimensions cannot exceed %s"), "2^31-1");
1552 }
1553 m_ = (n_ == 0) ? 0 : vlen_ / n_ + (vlen_ % n_ != 0);
1554 } else if (n_ < 0) {
1555 if (vlen_ > (Matrix_int_fast64_t) m_ * INT_MAX) {
1556 if (m_ == 0)
1557 error(_("nonempty vector supplied for empty matrix"));
1558 else
1559 error(_("dimensions cannot exceed %s"), "2^31-1");
1560 }
1561 n_ = (m_ == 0) ? 0 : vlen_ / m_ + (vlen_ % m_ != 0);
1562 }
1563
1564 Matrix_int_fast64_t mlen_ = (Matrix_int_fast64_t) m_ * n_;
1565 if (vlen_ <= 1)
1566 /* do nothing */ ;
1567 else if (mlen_ == 0)
1568 warning(_("nonempty vector supplied for empty matrix"));
1569 else if (vlen_ > mlen_)
1570 warning(_("vector length (%lld) exceeds matrix length (%d * %d)"),
1571 (long long) vlen_, m_, n_);
1572 else if (mlen_ % vlen_ != 0)
1573 warning(_("matrix length (%d * %d) is not a multiple of vector length (%lld)"),
1574 m_, n_, (long long) vlen_);
1575
1576 return
1577 vector_as_sparse(from, zzz_, ul, di, m_, n_, byrow_, dimnames);
1578}
1579
1580SEXP matrix_as_sparse(SEXP from, const char *zzz, char ul, char di,
1581 int trans)
1582{
1583 char cl[] = "...Matrix";
1584 cl[0] = typeToKind(TYPEOF(from));
1585 cl[1] = zzz[1];
1586 cl[2] = (zzz[1] == 'g') ? 'e' : ((zzz[1] == 's') ? 'y' : 'r');
1587#ifndef MATRIX_ENABLE_IMATRIX
1588 if (cl[0] == 'i')
1589 cl[0] = 'd';
1590#endif
1591 PROTECT_INDEX pid;
1592 PROTECT_WITH_INDEX(from, &pid);
1593 REPROTECT(from = matrix_as_dense(from, cl, ul, di, trans, 0), pid);
1594 REPROTECT(from = dense_as_sparse(from, cl, zzz[2]), pid);
1595 cl[2] = zzz[2];
1596 REPROTECT(from = sparse_as_kind(from, cl, zzz[0]), pid);
1597 UNPROTECT(1);
1598 return from;
1599}
1600
1601/* as(<matrix>, ".[gst][CRT]Matrix") */
1602SEXP R_matrix_as_sparse(SEXP from, SEXP zzz, SEXP uplo, SEXP diag,
1603 SEXP trans)
1604{
1605 switch (TYPEOF(from)) {
1606 case LGLSXP:
1607 case INTSXP:
1608 case REALSXP:
1609 case CPLXSXP:
1610 break;
1611 default:
1612 ERROR_INVALID_TYPE(from, __func__);
1613 break;
1614 }
1615
1616 const char *zzz_;
1617 if (TYPEOF(zzz) != STRSXP || LENGTH(zzz) < 1 ||
1618 (zzz = STRING_ELT(zzz, 0)) == NA_STRING ||
1619 (zzz_ = CHAR(zzz))[0] == '\0' ||
1620 (zzz_[1] != 'g' && zzz_[1] != 's' && zzz_[1] != 't') ||
1621 (zzz_[2] != 'C' && zzz_[2] != 'R' && zzz_[2] != 'T'))
1622 error(_("second argument of '%s' does not specify a subclass of %s"),
1623 __func__, "[CRT]sparseMatrix");
1624
1625 char ul = 'U', di = 'N';
1626 if (zzz_[1] != 'g') {
1627 if (TYPEOF(uplo) != STRSXP || LENGTH(uplo) < 1 ||
1628 (uplo = STRING_ELT(uplo, 0)) == NA_STRING ||
1629 ((ul = *CHAR(uplo)) != 'U' && ul != 'L'))
1630 error(_("'%s' must be \"%s\" or \"%s\""), "uplo", "U", "L");
1631 }
1632 if (zzz_[1] == 't') {
1633 if (TYPEOF(diag) != STRSXP || LENGTH(diag) < 1 ||
1634 (diag = STRING_ELT(diag, 0)) == NA_STRING ||
1635 ((di = *CHAR(diag)) != 'N' && di != 'U'))
1636 error(_("'%s' must be \"%s\" or \"%s\""), "diag", "N", "U");
1637 }
1638
1639 int trans_;
1640 if (TYPEOF(trans) != LGLSXP || LENGTH(trans) < 1 ||
1641 (trans_ = LOGICAL(trans)[0]) == NA_LOGICAL)
1642 error(_("'%s' must be %s or %s"), "trans", "TRUE", "FALSE");
1643
1644 return matrix_as_sparse(from, zzz_, ul, di, trans_);
1645}
1646
1647SEXP dense_as_sparse(SEXP from, const char *class, char repr)
1648{
1649 char cl[] = "...Matrix";
1650 cl[0] = class[0];
1651 cl[1] = class[1];
1652 cl[2] = repr;
1653 SEXP to = PROTECT(newObject(cl));
1654
1655 SEXP dim = PROTECT(GET_SLOT(from, Matrix_DimSym));
1656 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1];
1657 if (m != n || n > 0)
1658 SET_SLOT(to, Matrix_DimSym, dim);
1659 UNPROTECT(1); /* dim */
1660
1661 SEXP dimnames = PROTECT(GET_SLOT(from, Matrix_DimNamesSym));
1662 SET_SLOT(to, Matrix_DimNamesSym, dimnames);
1663 UNPROTECT(1); /* dimnames */
1664
1665 char ul = 'U', di = 'N';
1666 if (class[1] != 'g') {
1667 SEXP uplo = PROTECT(GET_SLOT(from, Matrix_uploSym));
1668 ul = *CHAR(STRING_ELT(uplo, 0));
1669 if (ul != 'U')
1670 SET_SLOT(to, Matrix_uploSym, uplo);
1671 UNPROTECT(1); /* uplo */
1672
1673 if (class[1] == 't') {
1674 SEXP diag = PROTECT(GET_SLOT(from, Matrix_diagSym));
1675 di = *CHAR(STRING_ELT(diag, 0));
1676 if (di != 'N')
1677 SET_SLOT(to, Matrix_diagSym, diag);
1678 UNPROTECT(1); /* diag */
1679 }
1680 }
1681
1682 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym)),
1683 p1, i1, j1;
1684 int i, j, *pp, *pi, *pj;
1685 R_xlen_t nnz = 0;
1686 p1 = i1 = j1 = NULL;
1687 pp = pi = pj = NULL;
1688
1689#define DAS_CASES(_MASK_) \
1690 do { \
1691 switch (class[0]) { \
1692 case 'l': \
1693 DAS_SUBCASES(int, LOGICAL, _MASK_, ISNZ_LOGICAL); \
1694 break; \
1695 case 'i': \
1696 DAS_SUBCASES(int, INTEGER, _MASK_, ISNZ_INTEGER); \
1697 break; \
1698 case 'd': \
1699 DAS_SUBCASES(double, REAL, _MASK_, ISNZ_REAL); \
1700 break; \
1701 case 'z': \
1702 DAS_SUBCASES(Rcomplex, COMPLEX, _MASK_, ISNZ_COMPLEX); \
1703 break; \
1704 default: \
1705 break; \
1706 } \
1707 } while (0)
1708
1709#define DAS_SUBCASES(_CTYPE_, _PTR_, _MASK_, _ISNZ_) \
1710 do { \
1711 _CTYPE_ *px0 = _PTR_(x0) ; \
1712 _MASK_(_CTYPE_ *px1 = _PTR_(x1)); \
1713 if (class[1] == 'g') \
1714 /* .geMatrix */ \
1715 DAS_SUBSUBCASES(DAS_LOOP_GEN2C, DAS_LOOP_GEN2R, DAS_LOOP_GEN2C, \
1716 _MASK_, _ISNZ_); \
1717 else if (class[2] != 'p' && di == 'N') \
1718 /* .syMatrix, non-unit diagonal .trMatrix */ \
1719 DAS_SUBSUBCASES(DAS_LOOP_TRN2C, DAS_LOOP_TRN2R, DAS_LOOP_TRN2C, \
1720 _MASK_, _ISNZ_); \
1721 else if (class[2] != 'p') \
1722 /* unit diagonal .trMatrix */ \
1723 DAS_SUBSUBCASES(DAS_LOOP_TRU2C, DAS_LOOP_TRU2R, DAS_LOOP_TRU2C, \
1724 _MASK_, _ISNZ_); \
1725 else if (di == 'N') \
1726 /* .spMatrix, non-unit diagonal .tpMatrix */ \
1727 DAS_SUBSUBCASES(DAS_LOOP_TPN2C, DAS_LOOP_TPN2R, DAS_LOOP_TPN2C, \
1728 _MASK_, _ISNZ_); \
1729 else \
1730 /* unit diagonal .tpMatrix */ \
1731 DAS_SUBSUBCASES(DAS_LOOP_TPU2C, DAS_LOOP_TPU2R, DAS_LOOP_TPU2C, \
1732 _MASK_, _ISNZ_); \
1733 } while (0)
1734
1735#undef DAS_SUBSUBCASES
1736#define DAS_SUBSUBCASES(_LOOP_C_, _LOOP_R_, _LOOP_T_, _MASK_, _ISNZ_) \
1737 do { \
1738 switch (cl[2]) { \
1739 case 'C': \
1740 _LOOP_C_(_ISNZ_, ++nnz, DAS_VALID2C); \
1741 break; \
1742 case 'R': \
1743 _LOOP_R_(_ISNZ_, ++nnz, DAS_VALID2R); \
1744 break; \
1745 case 'T': \
1746 _LOOP_T_(_ISNZ_, ++nnz, DAS_VALID2T); \
1747 break; \
1748 default: \
1749 break; \
1750 } \
1751 } while (0)
1752
1753#define DAS_LOOP_GEN2C(_ISNZ_, _DO_INNER_, _DO_OUTER_) \
1754 do { \
1755 for (j = 0; j < n; ++j) { \
1756 for (i = 0; i < m; ++i, ++px0) \
1757 if (_ISNZ_(*px0)) _DO_INNER_; \
1758 _DO_OUTER_; \
1759 } \
1760 } while (0)
1761
1762#define DAS_LOOP_GEN2R(_ISNZ_, _DO_INNER_, _DO_OUTER_) \
1763 do { \
1764 R_xlen_t mn1s = (R_xlen_t) m * n - 1; \
1765 for (i = 0; i < m; ++i, px0 -= mn1s) { \
1766 for (j = 0; j < n; ++j, px0 += m) \
1767 if (_ISNZ_(*px0)) _DO_INNER_; \
1768 _DO_OUTER_; \
1769 } \
1770 } while (0)
1771
1772#define DAS_LOOP_TRN2C(_ISNZ_, _DO_INNER_, _DO_OUTER_) \
1773 do { \
1774 if (ul == 'U') { \
1775 for (j = 0; j < n; px0 += n - (++j)) { \
1776 for (i = 0; i <= j; ++i, ++px0) \
1777 if (_ISNZ_(*px0)) _DO_INNER_; \
1778 _DO_OUTER_; \
1779 } \
1780 } else { \
1781 for (j = 0; j < n; px0 += (++j)) { \
1782 for (i = j; i < n; ++i, ++px0) \
1783 if (_ISNZ_(*px0)) _DO_INNER_; \
1784 _DO_OUTER_; \
1785 } \
1786 } \
1787 } while (0)
1788
1789#define DAS_LOOP_TRN2R(_ISNZ_, _DO_INNER_, _DO_OUTER_) \
1790 do { \
1791 R_xlen_t d; \
1792 if (ul == 'U') { \
1793 d = (R_xlen_t) n * n - 1; \
1794 for (i = 0; i < n; ++i, px0 -= (d -= n)) { \
1795 for (j = i; j < n; ++j, px0 += n) \
1796 if (_ISNZ_(*px0)) _DO_INNER_; \
1797 _DO_OUTER_; \
1798 } \
1799 } else { \
1800 d = -1; \
1801 for (i = 0; i < n; ++i, px0 -= (d += n)) { \
1802 for (j = 0; j <= i; ++j, px0 += n) \
1803 if (_ISNZ_(*px0)) _DO_INNER_; \
1804 _DO_OUTER_; \
1805 } \
1806 } \
1807 } while (0)
1808
1809#define DAS_LOOP_TRU2C(_ISNZ_, _DO_INNER_, _DO_OUTER_) \
1810 do { \
1811 if (ul == 'U') { \
1812 px0 += n; \
1813 for (j = 1; j < n; ++j) { \
1814 for (i = 0; i < j; ++i, ++px0) \
1815 if (_ISNZ_(*px0)) _DO_INNER_; \
1816 _DO_OUTER_; \
1817 px0 += n - j; \
1818 } \
1819 } else { \
1820 for (j = 0; j < n; ++j) { \
1821 px0 += j + 1; \
1822 for (i = j + 1; i < n; ++i, ++px0) \
1823 if (_ISNZ_(*px0)) _DO_INNER_; \
1824 _DO_OUTER_; \
1825 } \
1826 } \
1827 } while (0)
1828
1829#define DAS_LOOP_TRU2R(_ISNZ_, _DO_INNER_, _DO_OUTER_) \
1830 do { \
1831 R_xlen_t d; \
1832 if (ul == 'U') { \
1833 d = (R_xlen_t) n * (n - 1) - 1; \
1834 for (i = 0; i < n; ++i) { \
1835 for (j = i + 1; j < n; ++j) { \
1836 px0 += n; \
1837 if (_ISNZ_(*px0)) _DO_INNER_; \
1838 } \
1839 _DO_OUTER_; \
1840 px0 -= (d -= n); \
1841 } \
1842 } else { \
1843 ++px0; \
1844 d = -1; \
1845 for (i = 1; i < n; ++i) { \
1846 for (j = 0; j < i; ++j) { \
1847 if (_ISNZ_(*px0)) _DO_INNER_; \
1848 px0 += n; \
1849 } \
1850 _DO_OUTER_; \
1851 px0 -= (d += n); \
1852 } \
1853 } \
1854 } while (0)
1855
1856#define DAS_LOOP_TPN2C(_ISNZ_, _DO_INNER_, _DO_OUTER_) \
1857 do { \
1858 if (ul == 'U') { \
1859 for (j = 0; j < n; ++j) { \
1860 for (i = 0; i <= j; ++i, ++px0) \
1861 if (_ISNZ_(*px0)) _DO_INNER_; \
1862 _DO_OUTER_; \
1863 } \
1864 } else { \
1865 for (j = 0; j < n; ++j) { \
1866 for (i = j; i < n; ++i, ++px0) \
1867 if (_ISNZ_(*px0)) _DO_INNER_; \
1868 _DO_OUTER_; \
1869 } \
1870 } \
1871 } while (0)
1872
1873#define DAS_LOOP_TPN2R(_ISNZ_, _DO_INNER_, _DO_OUTER_) \
1874 do { \
1875 R_xlen_t d; \
1876 if (ul == 'U') { \
1877 d = PACKED_LENGTH(n) - 1; \
1878 for (i = 0; i < n; px0 -= (d -= (++i))) { \
1879 for (j = i; j < n; px0 += (++j)) \
1880 if (_ISNZ_(*px0)) _DO_INNER_; \
1881 _DO_OUTER_; \
1882 } \
1883 } else { \
1884 d = -1; \
1885 for (i = 0; i < n; px0 -= (d += n - (++i))) { \
1886 for (j = 0; j <= i; px0 += n - (++j)) \
1887 if (_ISNZ_(*px0)) _DO_INNER_; \
1888 _DO_OUTER_; \
1889 } \
1890 } \
1891 } while (0)
1892
1893#define DAS_LOOP_TPU2C(_ISNZ_, _DO_INNER_, _DO_OUTER_) \
1894 do { \
1895 if (ul == 'U') { \
1896 for (j = 1; j < n; ++j) { \
1897 ++px0; \
1898 for (i = 0; i < j; ++i, ++px0) \
1899 if (_ISNZ_(*px0)) _DO_INNER_; \
1900 _DO_OUTER_; \
1901 } \
1902 } else { \
1903 for (j = 0; j < n; ++j) { \
1904 ++px0; \
1905 for (i = j + 1; i < n; ++i, ++px0) \
1906 if (_ISNZ_(*px0)) _DO_INNER_; \
1907 _DO_OUTER_; \
1908 } \
1909 } \
1910 } while (0)
1911
1912#define DAS_LOOP_TPU2R(_ISNZ_, _DO_INNER_, _DO_OUTER_) \
1913 do { \
1914 R_xlen_t d; \
1915 if (ul == 'U') { \
1916 d = PACKED_LENGTH(n - 1) - 1; \
1917 for (i = 0; i < n; ++i) { \
1918 for (j = i + 1; j < n; ++j) { \
1919 px0 += j; \
1920 if (_ISNZ_(*px0)) _DO_INNER_; \
1921 } \
1922 _DO_OUTER_; \
1923 px0 -= (d -= i + 1); \
1924 } \
1925 } else { \
1926 ++px0; \
1927 d = -1; \
1928 for (i = 1; i < n; ++i) { \
1929 for (j = 0; j < i; ++j) { \
1930 if (_ISNZ_(*px0)) _DO_INNER_; \
1931 px0 += n - j - 1; \
1932 } \
1933 _DO_OUTER_; \
1934 px0 -= (d += n - i); \
1935 } \
1936 } \
1937 } while (0)
1938
1939#define DAS_VALID2T \
1940 if (nnz > INT_MAX) \
1941 error(_("attempt to construct %s with more than %s nonzero entries"), \
1942 "sparseMatrix", "2^31-1")
1943
1944#define DAS_VALID2C \
1945 do { DAS_VALID2T; else *(pp++) = (int) nnz; } while (0)
1946
1947#define DAS_VALID2R DAS_VALID2C
1948
1949 /* First we loop over the _nontrivial part_ of the denseMatrix 'from',
1950 by row ('R' case) or by column ('C' and 'T' cases), counting the
1951 nonzero entries and filling the 'p' slot of the result accordingly
1952 ('C' and 'R' cases) ...
1953 */
1954
1955 int nprotect = 2;
1956
1957 if (cl[2] != 'T') {
1958 int r = (cl[2] == 'C') ? n : m;
1959 PROTECT(p1 = allocVector(INTSXP, (R_xlen_t) r + 1));
1960 ++nprotect;
1961 SET_SLOT(to, Matrix_pSym, p1);
1962 pp = INTEGER(p1);
1963 *(pp++) = 0;
1964 if (r > 0 && di != 'N' && ul == ((cl[2] == 'C') ? 'U' : 'L'))
1965 *(pp++) = 0; /* first row or column skipped in these loops */
1966 }
1967 if (class[0] == 'n')
1968 DAS_SUBCASES(int, LOGICAL, HIDE, ISNZ_LOGICAL);
1969 else
1970 DAS_CASES(HIDE);
1971 if (cl[2] != 'R') {
1972 PROTECT(i1 = allocVector(INTSXP, nnz));
1973 ++nprotect;
1974 SET_SLOT(to, Matrix_iSym, i1);
1975 pi = INTEGER(i1);
1976 }
1977 if (cl[2] != 'C') {
1978 PROTECT(j1 = allocVector(INTSXP, nnz));
1979 ++nprotect;
1980 SET_SLOT(to, Matrix_jSym, j1);
1981 pj = INTEGER(j1);
1982 }
1983
1984#undef DAS_SUBSUBCASES
1985#define DAS_SUBSUBCASES(_LOOP_C_, _LOOP_R_, _LOOP_T_, _MASK_, _ISNZ_) \
1986 do { \
1987 switch (repr) { \
1988 case 'C': \
1989 _LOOP_C_(_ISNZ_, \
1990 do { \
1991 *(pi++) = i; \
1992 _MASK_(*(px1++) = *px0); \
1993 } while (0), ); \
1994 break; \
1995 case 'R': \
1996 _LOOP_R_(_ISNZ_, \
1997 do { \
1998 *(pj++) = j; \
1999 _MASK_(*(px1++) = *px0); \
2000 } while (0), ); \
2001 break; \
2002 case 'T': \
2003 _LOOP_T_(_ISNZ_, \
2004 do { \
2005 *(pi++) = i; \
2006 *(pj++) = j; \
2007 _MASK_(*(px1++) = *px0); \
2008 } while (0), ); \
2009 break; \
2010 default: \
2011 break; \
2012 } \
2013 } while (0)
2014
2015 /* Then we loop back over the same entries in order to fill
2016 the 'i', 'j', and 'x' slots of the result (whichever exist) ...
2017 */
2018
2019 if (class[0] == 'n')
2020 DAS_SUBCASES(int, LOGICAL, HIDE, ISNZ_LOGICAL);
2021 else {
2022 SEXP x1 = PROTECT(allocVector(TYPEOF(x0), nnz));
2023 SET_SLOT(to, Matrix_xSym, x1);
2024 DAS_CASES(SHOW);
2025 UNPROTECT(1); /* x1 */
2026 }
2027
2028#undef DAS_CASES
2029#undef DAS_SUBCASES
2030#undef DAS_SUBSUBCASES
2031#undef DAS_VALID2C
2032#undef DAS_VALID2R
2033#undef DAS_VALID2T
2034#undef DAS_LOOP_GEN2C
2035#undef DAS_LOOP_GEN2R
2036#undef DAS_LOOP_TRN2C
2037#undef DAS_LOOP_TRN2R
2038#undef DAS_LOOP_TRU2C
2039#undef DAS_LOOP_TRU2R
2040#undef DAS_LOOP_TPN2C
2041#undef DAS_LOOP_TPN2R
2042#undef DAS_LOOP_TPU2C
2043#undef DAS_LOOP_TPU2R
2044
2045 UNPROTECT(nprotect);
2046 return to;
2047}
2048
2049/* as(<denseMatrix>, "[CRT]sparseMatrix") */
2050SEXP R_dense_as_sparse(SEXP from, SEXP repr)
2051{
2052 static const char *valid[] = { VALID_DENSE, "" };
2053 int ivalid = R_check_class_etc(from, valid);
2054 if (ivalid < 0)
2055 ERROR_INVALID_CLASS(from, __func__);
2056
2057 char repr_;
2058 if (TYPEOF(repr) != STRSXP || LENGTH(repr) < 1 ||
2059 (repr = STRING_ELT(repr, 0)) == NA_STRING ||
2060 ((repr_ = CHAR(repr)[0]) != 'C' && repr_ != 'R' && repr_ != 'T'))
2061 error(_("invalid '%s' to '%s'"), "repr", __func__);
2062
2063 return dense_as_sparse(from, valid[ivalid], repr_);
2064}
2065
2066SEXP diagonal_as_sparse(SEXP from, const char *class,
2067 char kind, char shape, char repr, char ul)
2068{
2069 char cl[] = "...Matrix";
2070 cl[0] = (kind == '.') ? class[0] : ((kind == ',') ? ((class[0] == 'z') ? 'z' : 'd') : kind);
2071 cl[1] = shape;
2072 cl[2] = repr;
2073 SEXP to = PROTECT(newObject(cl));
2074
2075 SEXP dim = PROTECT(GET_SLOT(from, Matrix_DimSym));
2076 int n = INTEGER(dim)[0];
2077 if (n > 0)
2078 SET_SLOT(to, Matrix_DimSym, dim);
2079 UNPROTECT(1); /* dim */
2080
2081 SEXP dimnames = PROTECT(GET_SLOT(from, Matrix_DimNamesSym));
2082 if (cl[1] == 's')
2083 set_symmetrized_DimNames(to, dimnames, -1);
2084 else
2085 SET_SLOT(to, Matrix_DimNamesSym, dimnames);
2086 UNPROTECT(1); /* dimnames */
2087
2088 if (cl[1] != 'g' && ul != 'U') {
2089 SEXP uplo = PROTECT(mkString("L"));
2090 SET_SLOT(to, Matrix_uploSym, uplo);
2091 UNPROTECT(1); /* uplo */
2092 }
2093
2094 SEXP diag = PROTECT(GET_SLOT(from, Matrix_diagSym));
2095 char di = *CHAR(STRING_ELT(diag, 0));
2096 if (cl[1] == 't' && di != 'N')
2097 SET_SLOT(to, Matrix_diagSym, diag);
2098 UNPROTECT(1); /* diag */
2099
2100 if (cl[1] == 't' && di != 'N') {
2101 if (cl[2] != 'T') {
2102 SEXP p = PROTECT(allocVector(INTSXP, (R_xlen_t) n + 1));
2103 SET_SLOT(to, Matrix_pSym, p);
2104 Matrix_memset(INTEGER(p), 0, (R_xlen_t) n + 1, sizeof(int));
2105 UNPROTECT(1); /* p */
2106 }
2107 UNPROTECT(1); /* to */
2108 return to;
2109 }
2110
2111#define DAS_CASES(_MASK_) \
2112 do { \
2113 switch (cl[0]) { \
2114 case 'l': \
2115 DAS_LOOP(int, LOGICAL, _MASK_, ISNZ_LOGICAL, 1); \
2116 break; \
2117 case 'i': \
2118 DAS_LOOP(int, INTEGER, _MASK_, ISNZ_INTEGER, 1); \
2119 break; \
2120 case 'd': \
2121 DAS_LOOP(double, REAL, _MASK_, ISNZ_REAL, 1.0); \
2122 break; \
2123 case 'z': \
2124 DAS_LOOP(Rcomplex, COMPLEX, _MASK_, ISNZ_COMPLEX, Matrix_zone); \
2125 break; \
2126 default: \
2127 break; \
2128 } \
2129 } while (0)
2130
2131 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym));
2132 if (class[0] != cl[0]) {
2133 if (class[0] == 'n' && cl[0] == 'l')
2134 x0 = duplicate(x0);
2135 else
2136 x0 = coerceVector(x0, kindToType(cl[0]));
2137 if (class[0] == 'n')
2138 naToOne(x0);
2139 UNPROTECT(1); /* x0 */
2140 PROTECT(x0);
2141 }
2142
2143 int d, nnz;
2144 if (cl[2] != 'T') {
2145 SEXP p = PROTECT(allocVector(INTSXP, (R_xlen_t) n + 1));
2146 SET_SLOT(to, Matrix_pSym, p);
2147 int *pp = INTEGER(p);
2148 *(pp++) = 0;
2149 if (di == 'N') {
2150 nnz = 0;
2151
2152#undef DAS_LOOP
2153#define DAS_LOOP(_CTYPE_, _PTR_, _MASK_, _ISNZ_, _ONE_) \
2154 do { \
2155 _CTYPE_ *px0 = _PTR_(x0); \
2156 for (d = 0; d < n; ++d) { \
2157 if (_ISNZ_(*px0)) \
2158 ++nnz; \
2159 *(pp++) = nnz; \
2160 ++px0; \
2161 } \
2162 } while (0)
2163
2164 if (cl[0] == 'n')
2165 DAS_LOOP(int, LOGICAL, HIDE, ISNZ_LOGICAL, 1);
2166 else
2167 DAS_CASES(SHOW);
2168 } else {
2169 nnz = n;
2170 for (d = 1; d <= n; ++d)
2171 *(pp++) = d;
2172 }
2173 UNPROTECT(1); /* p */
2174 } else {
2175 if (di == 'N') {
2176 nnz = 0;
2177
2178#undef DAS_LOOP
2179#define DAS_LOOP(_CTYPE_, _PTR_, _MASK_, _ISNZ_, _ONE_) \
2180 do { \
2181 _CTYPE_ *px0 = _PTR_(x0); \
2182 for (d = 0; d < n; ++d) { \
2183 if (_ISNZ_(*px0)) \
2184 ++nnz; \
2185 ++px0; \
2186 } \
2187 } while (0)
2188
2189 if (cl[0] == 'n')
2190 DAS_LOOP(int, LOGICAL, HIDE, ISNZ_LOGICAL, 1);
2191 else
2192 DAS_CASES(SHOW);
2193 } else
2194 nnz = n;
2195 }
2196
2197 SEXP i1 = PROTECT(allocVector(INTSXP, nnz));
2198 if (cl[2] != 'T')
2199 SET_SLOT(to, (cl[2] == 'C') ? Matrix_iSym : Matrix_jSym, i1);
2200 else {
2201 SET_SLOT(to, Matrix_iSym, i1);
2202 SET_SLOT(to, Matrix_jSym, i1);
2203 }
2204 int *pi1 = INTEGER(i1);
2205
2206#undef DAS_LOOP
2207#define DAS_LOOP(_CTYPE_, _PTR_, _MASK_, _ISNZ_, _ONE_) \
2208 do { \
2209 _MASK_(_CTYPE_ *px1 = _PTR_(x1)); \
2210 if (di == 'N') { \
2211 _CTYPE_ *px0 = _PTR_(x0); \
2212 for (d = 0; d < n; ++d) { \
2213 if (_ISNZ_(*px0)) { \
2214 *(pi1++) = d; \
2215 _MASK_(*(px1++) = *px0); \
2216 } \
2217 ++px0; \
2218 } \
2219 } else { \
2220 for (d = 0; d < n; ++d) { \
2221 *(pi1++) = d; \
2222 _MASK_(*(px1++) = _ONE_); \
2223 } \
2224 } \
2225 } while (0)
2226
2227 if (cl[0] == 'n')
2228 DAS_LOOP(int, LOGICAL, HIDE, ISNZ_LOGICAL, 1);
2229 else if (di == 'N' && nnz == n) {
2230 SET_SLOT(to, Matrix_xSym, x0);
2231 DAS_CASES(HIDE);
2232 } else {
2233 SEXP x1 = PROTECT(allocVector(TYPEOF(x0), nnz));
2234 SET_SLOT(to, Matrix_xSym, x1);
2235 DAS_CASES(SHOW);
2236 UNPROTECT(1); /* x1 */
2237 }
2238
2239 UNPROTECT(3); /* i1, x0, to */
2240 return to;
2241}
2242
2243/* as(<diagonalMatrix>, ".[gst][CRT]Matrix") */
2244SEXP R_diagonal_as_sparse(SEXP from,
2245 SEXP kind, SEXP shape, SEXP repr, SEXP uplo)
2246{
2247 static const char *valid[] = { VALID_DIAGONAL, "" };
2248 int ivalid = R_check_class_etc(from, valid);
2249 if (ivalid < 0)
2250 ERROR_INVALID_CLASS(from, __func__);
2251
2252 char kind_;
2253 if (TYPEOF(kind) != STRSXP || LENGTH(kind) < 1 ||
2254 (kind = STRING_ELT(kind, 0)) == NA_STRING ||
2255 (kind_ = CHAR(kind)[0]) == '\0')
2256 error(_("invalid '%s' to '%s'"), "kind", __func__);
2257
2258 char shape_;
2259 if (TYPEOF(shape) != STRSXP || LENGTH(shape) < 1 ||
2260 (shape = STRING_ELT(shape, 0)) == NA_STRING ||
2261 ((shape_ = CHAR(shape)[0]) != 'g' && shape_ != 's' && shape_ != 't'))
2262 error(_("invalid '%s' to '%s'"), "shape", __func__);
2263
2264 char repr_;
2265 if (TYPEOF(repr) != STRSXP || LENGTH(repr) < 1 ||
2266 (repr = STRING_ELT(repr, 0)) == NA_STRING ||
2267 ((repr_ = CHAR(repr)[0]) != 'C' && repr_ != 'R' && repr_ != 'T'))
2268 error(_("invalid '%s' to '%s'"), "repr", __func__);
2269
2270 char ul = 'U';
2271 if (shape_ != 'g') {
2272 if (TYPEOF(uplo) != STRSXP || LENGTH(uplo) < 1 ||
2273 (uplo = STRING_ELT(uplo, 0)) == NA_STRING ||
2274 ((ul = *CHAR(uplo)) != 'U' && ul != 'L'))
2275 error(_("'%s' must be \"%s\" or \"%s\""), "uplo", "U", "L");
2276 }
2277
2278 return diagonal_as_sparse(from, valid[ivalid], kind_, shape_, repr_, ul);
2279}
2280
2281SEXP index_as_sparse(SEXP from, const char *class, char kind, char repr)
2282{
2283 SEXP margin = PROTECT(GET_SLOT(from, Matrix_marginSym));
2284 int mg = INTEGER(margin)[0] - 1;
2285 UNPROTECT(1); /* margin */
2286
2287 char cl[] = ".g.Matrix";
2288 cl[0] = (kind == '.') ? 'n' : ((kind == ',') ? 'd' : kind);
2289 cl[2] = (repr == '.') ? ((mg == 0) ? 'R' : 'C') : repr;
2290 SEXP to = PROTECT(newObject(cl));
2291
2292 SEXP dim = PROTECT(GET_SLOT(from, Matrix_DimSym));
2293 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1],
2294 r = (mg == 0) ? m : n, s = (mg == 0) ? n : m;
2295 if (m != n || n > 0)
2296 SET_SLOT(to, Matrix_DimSym, dim);
2297 UNPROTECT(1); /* dim */
2298
2299 SEXP dimnames = PROTECT(GET_SLOT(from, Matrix_DimNamesSym));
2300 SET_SLOT(to, Matrix_DimNamesSym, dimnames);
2301 UNPROTECT(1); /* dimnames */
2302
2303 SEXP perm = PROTECT(GET_SLOT(from, Matrix_permSym));
2304 int *pperm = INTEGER(perm);
2305
2306 if (cl[2] == 'T') {
2307 SEXP i = PROTECT(allocVector(INTSXP, r)),
2308 j = PROTECT(allocVector(INTSXP, r));
2309 int k, *pi = INTEGER(i), *pj = INTEGER(j);
2310 for (k = 0; k < r; ++k) {
2311 *(pi++) = k;
2312 *(pj++) = *(pperm++) - 1;
2313 }
2314 SET_SLOT(to, Matrix_iSym, (mg == 0) ? i : j);
2315 SET_SLOT(to, Matrix_jSym, (mg == 0) ? j : i);
2316 } else if ((cl[2] == 'C') == (mg != 0)) {
2317 SEXP p = PROTECT(allocVector(INTSXP, (R_xlen_t) r + 1)),
2318 i = PROTECT(allocVector(INTSXP, r));
2319 int k, *pp = INTEGER(p), *pi = INTEGER(i);
2320 for (k = 0; k < r; ++k) {
2321 *(pp++) = k;
2322 *(pi++) = *(pperm++) - 1;
2323 }
2324 *pp = r;
2325 SET_SLOT(to, Matrix_pSym, p);
2326 SET_SLOT(to, (mg != 0) ? Matrix_iSym : Matrix_jSym, i);
2327 } else {
2328 SEXP p = PROTECT(allocVector(INTSXP, (R_xlen_t) s + 1));
2329 int k, *pp = INTEGER(p);
2330 Matrix_memset(pp, 0, (R_xlen_t) s + 1, sizeof(int));
2331 for (k = 0; k < r; ++k)
2332 ++pp[pperm[k]];
2333 for (k = 0; k < s; ++k)
2334 pp[k + 1] += pp[k];
2335 SEXP j = PROTECT(allocVector(INTSXP, r));
2336 int *pj = INTEGER(j), *work;
2337 Matrix_Calloc(work, s, int);
2338 Matrix_memcpy(work, pp, s, sizeof(int));
2339 --work;
2340 for (k = 0; k < r; ++k)
2341 pj[work[pperm[k]]++] = k;
2342 ++work;
2343 Matrix_Free(work, s);
2344 SET_SLOT(to, Matrix_pSym, p);
2345 SET_SLOT(to, (mg != 0) ? Matrix_jSym : Matrix_iSym, j);
2346 }
2347 UNPROTECT(2);
2348
2349 if (cl[0] != 'n') {
2350 SEXP x = PROTECT(allocVector(kindToType(cl[0]), r));
2351 SET_SLOT(to, Matrix_xSym, x);
2352
2353#define IAS_SUBCASES(_CTYPE_, _PTR_, _ONE_) \
2354 do { \
2355 _CTYPE_ *px = _PTR_(x); \
2356 for (int k = 0; k < r; ++k) \
2357 *(px++) = _ONE_; \
2358 } while (0)
2359
2360 switch (cl[0]) {
2361 case 'l':
2362 IAS_SUBCASES(int, LOGICAL, 1);
2363 break;
2364 case 'i':
2365 IAS_SUBCASES(int, INTEGER, 1);
2366 break;
2367 case 'd':
2368 IAS_SUBCASES(double, REAL, 1.0);
2369 break;
2370 case 'z':
2371 IAS_SUBCASES(Rcomplex, COMPLEX, Matrix_zone);
2372 break;
2373 default:
2374 break;
2375 }
2376 UNPROTECT(1); /* x */
2377 }
2378
2379#undef IAS_SUBCASES
2380
2381 UNPROTECT(2); /* perm, to */
2382 return to;
2383}
2384
2385/* as(<indMatrix>, ".g[CRT]Matrix") */
2386SEXP R_index_as_sparse(SEXP from, SEXP kind, SEXP repr)
2387{
2388 static const char *valid[] = { "indMatrix", "pMatrix" };
2389 int ivalid = R_check_class_etc(from, valid);
2390 if (ivalid < 0)
2391 ERROR_INVALID_CLASS(from, __func__);
2392
2393 char kind_;
2394 if (TYPEOF(kind) != STRSXP || LENGTH(kind) < 1 ||
2395 (kind = STRING_ELT(kind, 0)) == NA_STRING ||
2396 (kind_ = CHAR(kind)[0]) == '\0')
2397 error(_("invalid '%s' to '%s'"), "kind", __func__);
2398
2399 char repr_;
2400 if (TYPEOF(repr) != STRSXP || LENGTH(repr) < 1 ||
2401 (repr = STRING_ELT(repr, 0)) == NA_STRING ||
2402 ((repr_ = CHAR(repr)[0]) != '.' &&
2403 repr_ != 'C' && repr_ != 'R' && repr_ != 'T'))
2404 error(_("invalid '%s' to '%s'"), "repr", __func__);
2405
2406 return index_as_sparse(from, valid[ivalid], kind_, repr_);
2407}
2408
2409SEXP dense_as_kind(SEXP from, const char *class, char kind, int new)
2410{
2411 if (kind == '.')
2412 kind = class[0];
2413 else if (kind == ',')
2414 kind = (class[0] == 'z') ? 'z' : 'd';
2415 if (kind == class[0])
2416 return from;
2417 SEXPTYPE tt = kindToType(kind);
2418
2419 char cl[] = "...Matrix";
2420 cl[0] = kind;
2421 cl[1] = class[1];
2422 cl[2] = class[2];
2423 SEXP to = PROTECT(newObject(cl));
2424
2425 SEXP dim = PROTECT(GET_SLOT(from, Matrix_DimSym));
2426 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1];
2427 if (m != n || n > 0)
2428 SET_SLOT(to, Matrix_DimSym, dim);
2429 UNPROTECT(1); /* dim */
2430
2431 SEXP dimnames = PROTECT(GET_SLOT(from, Matrix_DimNamesSym));
2432 SET_SLOT(to, Matrix_DimNamesSym, dimnames);
2433 UNPROTECT(1); /* dimnames */
2434
2435 if (class[1] != 'g') {
2436 SEXP uplo = PROTECT(GET_SLOT(from, Matrix_uploSym));
2437 char ul = *CHAR(STRING_ELT(uplo, 0));
2438 if (ul != 'U')
2439 SET_SLOT(to, Matrix_uploSym, uplo);
2440 UNPROTECT(1); /* uplo */
2441 if (class[1] == 't') {
2442 SEXP diag = PROTECT(GET_SLOT(from, Matrix_diagSym));
2443 char di = *CHAR(STRING_ELT(diag, 0));
2444 if (di != 'N')
2445 SET_SLOT(to, Matrix_diagSym, diag);
2446 UNPROTECT(1); /* diag */
2447 }
2448 }
2449
2450 PROTECT_INDEX pid;
2451 SEXP x;
2452 PROTECT_WITH_INDEX(x = GET_SLOT(from, Matrix_xSym), &pid);
2453 if (TYPEOF(x) != tt) {
2454 REPROTECT(x = coerceVector(x, tt), pid);
2455 if (class[0] == 'n')
2456 /* n->[idz] */
2457 naToOne(x);
2458 } else if (new) {
2459 REPROTECT(x = duplicate(x), pid);
2460 if (class[0] == 'n')
2461 /* n->l */
2462 naToOne(x);
2463 } else {
2464 if (class[0] == 'n') {
2465 /* n->l */
2466 R_xlen_t i, len = XLENGTH(x);
2467 int *px = LOGICAL(x);
2468 for (i = 0; i < len; ++i, ++px) {
2469 if (*px == NA_LOGICAL) {
2470 REPROTECT(x = duplicate(x), pid);
2471 naToOne(x);
2472 break;
2473 }
2474 }
2475 }
2476 }
2477 SET_SLOT(to, Matrix_xSym, x);
2478 UNPROTECT(2); /* x, to */
2479 return to;
2480}
2481
2482/* as(<denseMatrix>, "[nlidz]Matrix") */
2483SEXP R_dense_as_kind(SEXP from, SEXP kind)
2484{
2485 static const char *valid[] = { VALID_DENSE, "" };
2486 int ivalid = R_check_class_etc(from, valid);
2487 if (ivalid < 0)
2488 ERROR_INVALID_CLASS(from, __func__);
2489
2490 char kind_;
2491 if (TYPEOF(kind) != STRSXP || LENGTH(kind) < 1 ||
2492 (kind = STRING_ELT(kind, 0)) == NA_STRING ||
2493 (kind_ = CHAR(kind)[0]) == '\0')
2494 error(_("invalid '%s' to '%s'"), "kind", __func__);
2495
2496 return dense_as_kind(from, valid[ivalid], kind_, 0);
2497}
2498
2499SEXP sparse_as_kind(SEXP from, const char *class, char kind)
2500{
2501 if (kind == '.')
2502 kind = class[0];
2503 else if (kind == ',')
2504 kind = (class[0] == 'z') ? 'z' : 'd';
2505 if (kind == class[0])
2506 return from;
2507 SEXPTYPE tt = kindToType(kind);
2508
2509 if (class[2] == 'T' && (class[0] == 'n' || class[0] == 'l') &&
2510 kind != 'n' && kind != 'l') {
2511 /* defined in ./sparse.c : */
2512 SEXP Tsparse_aggregate(SEXP);
2513 from = Tsparse_aggregate(from);
2514 }
2515 PROTECT(from);
2516
2517 char cl[] = "...Matrix";
2518 cl[0] = kind;
2519 cl[1] = class[1];
2520 cl[2] = class[2];
2521 SEXP to = PROTECT(newObject(cl));
2522
2523 SEXP dim = PROTECT(GET_SLOT(from, Matrix_DimSym));
2524 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1];
2525 if (m != n || n > 0)
2526 SET_SLOT(to, Matrix_DimSym, dim);
2527 UNPROTECT(1); /* dim */
2528
2529 SEXP dimnames = PROTECT(GET_SLOT(from, Matrix_DimNamesSym));
2530 SET_SLOT(to, Matrix_DimNamesSym, dimnames);
2531 UNPROTECT(1); /* dimnames */
2532
2533 if (class[1] != 'g') {
2534 SEXP uplo = PROTECT(GET_SLOT(from, Matrix_uploSym));
2535 char ul = *CHAR(STRING_ELT(uplo, 0));
2536 if (ul != 'U')
2537 SET_SLOT(to, Matrix_uploSym, uplo);
2538 UNPROTECT(1); /* uplo */
2539 if (class[1] == 't') {
2540 SEXP diag = PROTECT(GET_SLOT(from, Matrix_diagSym));
2541 char di = *CHAR(STRING_ELT(diag, 0));
2542 if (di != 'N')
2543 SET_SLOT(to, Matrix_diagSym, diag);
2544 UNPROTECT(1); /* diag */
2545 }
2546 }
2547
2548 R_xlen_t nnz = -1;
2549 if (class[2] != 'T') {
2550 SEXP p = PROTECT(GET_SLOT(from, Matrix_pSym));
2551 SET_SLOT(to, Matrix_pSym, p);
2552 nnz = INTEGER(p)[XLENGTH(p) - 1];
2553 UNPROTECT(1); /* p */
2554 }
2555 if (class[2] != 'R') {
2556 SEXP i = PROTECT(GET_SLOT(from, Matrix_iSym));
2557 SET_SLOT(to, Matrix_iSym, i);
2558 if (nnz < 1)
2559 nnz = XLENGTH(i);
2560 UNPROTECT(1); /* i */
2561 }
2562 if (class[2] != 'C') {
2563 SEXP j = PROTECT(GET_SLOT(from, Matrix_jSym));
2564 SET_SLOT(to, Matrix_jSym, j);
2565 UNPROTECT(1); /* j */
2566 }
2567 if (class[0] == 'n') {
2568 SEXP x = PROTECT(allocVector(tt, nnz));
2569 switch (tt) {
2570 case LGLSXP:
2571 {
2572 int *px = LOGICAL(x);
2573 while (nnz--)
2574 *(px++) = 1;
2575 break;
2576 }
2577 case INTSXP:
2578 {
2579 int *px = INTEGER(x);
2580 while (nnz--)
2581 *(px++) = 1;
2582 break;
2583 }
2584 case REALSXP:
2585 {
2586 double *px = REAL(x);
2587 while (nnz--)
2588 *(px++) = 1.0;
2589 break;
2590 }
2591 case CPLXSXP:
2592 {
2593 Rcomplex *px = COMPLEX(x);
2594 while (nnz--)
2595 *(px++) = Matrix_zone;
2596 break;
2597 }
2598 default:
2599 break;
2600 }
2601 SET_SLOT(to, Matrix_xSym, x);
2602 UNPROTECT(1); /* x */
2603 } else if (kind != 'n') {
2604 PROTECT_INDEX pid;
2605 SEXP x;
2606 PROTECT_WITH_INDEX(x = GET_SLOT(from, Matrix_xSym), &pid);
2607 REPROTECT(x = coerceVector(x, tt), pid);
2608 SET_SLOT(to, Matrix_xSym, x);
2609 UNPROTECT(1); /* x */
2610 }
2611
2612 UNPROTECT(2); /* to, from */
2613 return to;
2614}
2615
2616/* as(<denseMatrix>, "[nlidz]Matrix") */
2617SEXP R_sparse_as_kind(SEXP from, SEXP kind)
2618{
2619 static const char *valid[] = {
2620 VALID_CSPARSE, VALID_RSPARSE, VALID_TSPARSE, "" };
2621 int ivalid = R_check_class_etc(from, valid);
2622 if (ivalid < 0)
2623 ERROR_INVALID_CLASS(from, __func__);
2624
2625 char kind_;
2626 if (TYPEOF(kind) != STRSXP || LENGTH(kind) < 1 ||
2627 (kind = STRING_ELT(kind, 0)) == NA_STRING ||
2628 (kind_ = CHAR(kind)[0]) == '\0')
2629 error(_("invalid '%s' to '%s'"), "kind", __func__);
2630
2631 return sparse_as_kind(from, valid[ivalid], kind_);
2632}
2633
2634SEXP diagonal_as_kind(SEXP from, const char *class, char kind)
2635{
2636 if (kind == '.')
2637 kind = class[0];
2638 else if (kind == ',')
2639 kind = (class[0] == 'z') ? 'z' : 'd';
2640 if (kind == class[0])
2641 return from;
2642 SEXPTYPE tt = kindToType(kind);
2643
2644 char cl[] = ".diMatrix";
2645 cl[0] = kind;
2646 SEXP to = PROTECT(newObject(cl));
2647
2648 SEXP dim = PROTECT(GET_SLOT(from, Matrix_DimSym));
2649 int n = INTEGER(dim)[0];
2650 if (n > 0)
2651 SET_SLOT(to, Matrix_DimSym, dim);
2652 UNPROTECT(1); /* dim */
2653
2654 SEXP dimnames = PROTECT(GET_SLOT(from, Matrix_DimNamesSym));
2655 SET_SLOT(to, Matrix_DimNamesSym, dimnames);
2656 UNPROTECT(1); /* dimnames */
2657
2658 SEXP diag = PROTECT(GET_SLOT(from, Matrix_diagSym));
2659 char di = *CHAR(STRING_ELT(diag, 0));
2660 if (di != 'N')
2661 SET_SLOT(to, Matrix_diagSym, diag);
2662 UNPROTECT(1); /* diag */
2663
2664 if (di == 'N') {
2665 PROTECT_INDEX pid;
2666 SEXP x;
2667 PROTECT_WITH_INDEX(x = GET_SLOT(from, Matrix_xSym), &pid);
2668 if (TYPEOF(x) == tt) {
2669 if (class[0] == 'n') {
2670 /* n->l */
2671 R_xlen_t i, len = XLENGTH(x);
2672 int *px = LOGICAL(x);
2673 for (i = 0; i < len; ++i, ++px) {
2674 if (*px == NA_LOGICAL) {
2675 REPROTECT(x = duplicate(x), pid);
2676 naToOne(x);
2677 break;
2678 }
2679 }
2680 }
2681 } else {
2682 REPROTECT(x = coerceVector(x, tt), pid);
2683 if (class[0] == 'n')
2684 /* n->[idz] */
2685 naToOne(x);
2686 }
2687 SET_SLOT(to, Matrix_xSym, x);
2688 UNPROTECT(1); /* x */
2689 }
2690
2691 UNPROTECT(1); /* to */
2692 return to;
2693}
2694
2695/* as(<diagonalMatrix>, "[nlidz]Matrix") */
2696SEXP R_diagonal_as_kind(SEXP from, SEXP kind)
2697{
2698 static const char *valid[] = { VALID_DIAGONAL, "" };
2699 int ivalid = R_check_class_etc(from, valid);
2700 if (ivalid < 0)
2701 ERROR_INVALID_CLASS(from, __func__);
2702
2703 char kind_;
2704 if (TYPEOF(kind) != STRSXP || LENGTH(kind) < 1 ||
2705 (kind = STRING_ELT(kind, 0)) == NA_STRING ||
2706 (kind_ = CHAR(kind)[0]) == '\0')
2707 error(_("invalid '%s' to '%s'"), "kind", __func__);
2708
2709 return diagonal_as_kind(from, valid[ivalid], kind_);
2710}
2711
2712SEXP index_as_kind(SEXP from, const char *class, char kind)
2713{
2714 return index_as_sparse(from, class, kind, '.');
2715}
2716
2717/* as(<indMatrix>, "[nlidz]Matrix") */
2718SEXP R_index_as_kind(SEXP from, SEXP kind)
2719{
2720 static const char *valid[] = { "indMatrix", "pMatrix" };
2721 int ivalid = R_check_class_etc(from, valid);
2722 if (ivalid < 0)
2723 ERROR_INVALID_CLASS(from, __func__);
2724
2725 char kind_;
2726 if (TYPEOF(kind) != STRSXP || LENGTH(kind) < 1 ||
2727 (kind = STRING_ELT(kind, 0)) == NA_STRING ||
2728 (kind_ = CHAR(kind)[0]) == '\0')
2729 error(_("invalid '%s' to '%s'"), "kind", __func__);
2730
2731 return index_as_kind(from, valid[ivalid], kind_);
2732}
2733
2734SEXP dense_as_general(SEXP from, const char *class, int new)
2735{
2736 if (class[1] == 'g')
2737 return from;
2738
2739 char cl[] = ".geMatrix";
2740 cl[0] = class[0];
2741 SEXP to = PROTECT(newObject(cl));
2742
2743 SEXP dim = PROTECT(GET_SLOT(from, Matrix_DimSym));
2744 int n = INTEGER(dim)[0];
2745 if (n > 0)
2746 SET_SLOT(to, Matrix_DimSym, dim);
2747 UNPROTECT(1); /* dim */
2748
2749 SEXP dimnames = PROTECT(GET_SLOT(from, Matrix_DimNamesSym));
2750 if (class[1] != 's')
2751 SET_SLOT(to, Matrix_DimNamesSym, dimnames);
2752 else
2753 set_symmetrized_DimNames(to, dimnames, -1);
2754 UNPROTECT(1); /* dimnames */
2755
2756 SEXP uplo = PROTECT(GET_SLOT(from, Matrix_uploSym));
2757 char ul = *CHAR(STRING_ELT(uplo, 0)), di = 'N';
2758 UNPROTECT(1); /* uplo */
2759
2760 if (class[1] == 's') {
2761 SEXP factors = PROTECT(GET_SLOT(from, Matrix_factorsSym));
2762 if (LENGTH(factors) > 0)
2763 SET_SLOT(to, Matrix_factorsSym, factors);
2764 UNPROTECT(1); /* factors */
2765 } else {
2766 SEXP diag = PROTECT(GET_SLOT(from, Matrix_diagSym));
2767 di = *CHAR(STRING_ELT(diag, 0));
2768 UNPROTECT(1); /* diag */
2769 }
2770
2771 if ((Matrix_int_fast64_t) n * n > R_XLEN_T_MAX)
2772 error(_("attempt to allocate vector of length exceeding %s"),
2773 "R_XLEN_T_MAX");
2774 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym)), x1 = x0;
2775 int nprotect = 2;
2776 if (class[2] == 'p' || new) {
2777 PROTECT(x1 = allocVector(TYPEOF(x0), (R_xlen_t) n * n));
2778 ++nprotect;
2779 }
2780 SET_SLOT(to, Matrix_xSym, x1);
2781
2782#define DAG_SUBCASES(_PREFIX_, _CTYPE_, _PTR_) \
2783 do { \
2784 _CTYPE_ *px0 = _PTR_(x0), *px1 = _PTR_(x1); \
2785 if (class[2] == 'p') \
2786 _PREFIX_ ## unpack1(px1, px0, n, ul, di); \
2787 else if (new) \
2788 Matrix_memcpy(px1, px0, (R_xlen_t) n * n, sizeof(_CTYPE_)); \
2789 if (class[1] == 's') \
2790 _PREFIX_ ## syforce2(px1, n, ul); \
2791 else \
2792 _PREFIX_ ## trforce2(px1, n, n, ul, di); \
2793 } while (0)
2794
2795 switch (class[0]) {
2796 case 'n':
2797 case 'l':
2798 DAG_SUBCASES(i, int, LOGICAL);
2799 break;
2800 case 'i':
2801 DAG_SUBCASES(i, int, INTEGER);
2802 break;
2803 case 'd':
2804 DAG_SUBCASES(d, double, REAL);
2805 break;
2806 case 'z':
2807 DAG_SUBCASES(z, Rcomplex, COMPLEX);
2808 break;
2809 default:
2810 break;
2811 }
2812
2813#undef DAG_CASES
2814#undef DAG_SUBCASES
2815
2816 UNPROTECT(nprotect);
2817 return to;
2818}
2819
2820/* as(<denseMatrix>, "generaMatrix") */
2821SEXP R_dense_as_general(SEXP from)
2822{
2823 static const char *valid[] = { VALID_DENSE, "" };
2824 int ivalid = R_check_class_etc(from, valid);
2825 if (ivalid < 0)
2826 ERROR_INVALID_CLASS(from, __func__);
2827
2828 return dense_as_general(from, valid[ivalid], 1);
2829}
2830
2831SEXP sparse_as_general(SEXP from, const char *class)
2832{
2833 if (class[1] == 'g')
2834 return from;
2835
2836 char cl[] = ".g.Matrix";
2837 cl[0] = class[0];
2838 cl[2] = class[2];
2839 SEXP to = PROTECT(newObject(cl));
2840
2841 SEXP dim = PROTECT(GET_SLOT(from, Matrix_DimSym));
2842 int n = INTEGER(dim)[0];
2843 if (n > 0)
2844 SET_SLOT(to, Matrix_DimSym, dim);
2845 UNPROTECT(1); /* dim */
2846
2847 SEXP dimnames = PROTECT(GET_SLOT(from, Matrix_DimNamesSym));
2848 if (class[1] == 's')
2849 set_symmetrized_DimNames(to, dimnames, -1);
2850 else
2851 SET_SLOT(to, Matrix_DimNamesSym, dimnames);
2852 UNPROTECT(1); /* dimnames */
2853
2854 if (class[1] == 's') {
2855 SEXP factors = PROTECT(GET_SLOT(from, Matrix_factorsSym));
2856 if (LENGTH(factors) > 0)
2857 SET_SLOT(to, Matrix_factorsSym, factors);
2858 UNPROTECT(1); /* factors */
2859 } else {
2860 SEXP diag = PROTECT(GET_SLOT(from, Matrix_diagSym));
2861 char di = *CHAR(STRING_ELT(diag, 0));
2862 UNPROTECT(1); /* diag */
2863
2864 if (di == 'N') {
2865 if (class[2] != 'T') {
2866 SEXP p = PROTECT(GET_SLOT(from, Matrix_pSym));
2867 SET_SLOT(to, Matrix_pSym, p);
2868 UNPROTECT(1); /* p */
2869 }
2870 if (class[2] != 'R') {
2871 SEXP i = PROTECT(GET_SLOT(from, Matrix_iSym));
2872 SET_SLOT(to, Matrix_iSym, i);
2873 UNPROTECT(1); /* i */
2874 }
2875 if (class[2] != 'C') {
2876 SEXP j = PROTECT(GET_SLOT(from, Matrix_jSym));
2877 SET_SLOT(to, Matrix_jSym, j);
2878 UNPROTECT(1); /* j */
2879 }
2880 if (class[0] != 'n') {
2881 SEXP x = PROTECT(GET_SLOT(from, Matrix_xSym));
2882 SET_SLOT(to, Matrix_xSym, x);
2883 UNPROTECT(1); /* x */
2884 }
2885 UNPROTECT(1); /* to */
2886 return to;
2887 }
2888 }
2889
2890 SEXP uplo = PROTECT(GET_SLOT(from, Matrix_uploSym));
2891 char ul = *CHAR(STRING_ELT(uplo, 0));
2892 UNPROTECT(1); /* uplo */
2893
2894#define ASSIGN_COMPLEX_JJ(_X_, _Y_) \
2895 do { _X_.r = _Y_.r; _X_.i = 0.0; } while (0)
2896
2897#define ASSIGN_COMPLEX_JI(_X_, _Y_) \
2898 do { _X_.r = _Y_.r; _X_.i = -_Y_.i; } while (0)
2899
2900#define SAG_CASES \
2901 do { \
2902 switch (class[0]) { \
2903 case 'l': \
2904 SAG_SUBCASES(int, LOGICAL, SHOW, 1, \
2905 ASSIGN_REAL, ASSIGN_REAL); \
2906 break; \
2907 case 'i': \
2908 SAG_SUBCASES(int, INTEGER, SHOW, 1, \
2909 ASSIGN_REAL, ASSIGN_REAL); \
2910 break; \
2911 case 'd': \
2912 SAG_SUBCASES(double, REAL, SHOW, 1.0, \
2913 ASSIGN_REAL, ASSIGN_REAL); \
2914 break; \
2915 case 'z': \
2916 SAG_SUBCASES(Rcomplex, COMPLEX, SHOW, Matrix_zone, \
2917 ASSIGN_COMPLEX_JJ, ASSIGN_COMPLEX_JI); \
2918 break; \
2919 default: \
2920 break; \
2921 } \
2922 } while (0)
2923
2924 if (class[2] != 'T') {
2925
2926 SEXP iSym = (class[2] == 'C') ? Matrix_iSym : Matrix_jSym,
2927 p0 = PROTECT(GET_SLOT(from, Matrix_pSym)),
2928 p1 = PROTECT(allocVector(INTSXP, (R_xlen_t) n + 1)),
2929 i0 = PROTECT(GET_SLOT(from, iSym));
2930 int j, k, kend,
2931 *pp0 = INTEGER(p0), *pp1 = INTEGER(p1), *pi0 = INTEGER(i0);
2932 SET_SLOT(to, Matrix_pSym, p1);
2933 pp0++; *(pp1++) = 0;
2934
2935 if (class[1] == 's') {
2936 Matrix_memset(pp1, 0, n, sizeof(int));
2937 for (j = 0, k = 0; j < n; ++j) {
2938 kend = pp0[j];
2939 while (k < kend) {
2940 if (pi0[k] != j)
2941 ++pp1[pi0[k]];
2942 ++k;
2943 }
2944 }
2945 for (j = 1; j < n; ++j)
2946 pp1[j] += pp1[j - 1];
2947 if (pp1[n - 1] > INT_MAX - pp0[n - 1])
2948 error(_("attempt to construct %s with more than %s nonzero entries"),
2949 "sparseMatrix", "2^31-1");
2950 for (j = 0; j < n; ++j)
2951 pp1[j] += pp0[j];
2952 } else {
2953 if (n > INT_MAX - pp0[n - 1])
2954 error(_("attempt to construct %s with more than %s nonzero entries"),
2955 "sparseMatrix", "2^31-1");
2956 for (j = 0; j < n; ++j)
2957 pp1[j] = pp0[j] + j + 1;
2958 }
2959
2960 SEXP i1 = PROTECT(allocVector(INTSXP, pp1[n - 1]));
2961 int *pi1 = INTEGER(i1);
2962 SET_SLOT(to, iSym, i1);
2963
2964#undef SAG_SUBCASES
2965#define SAG_SUBCASES(_CTYPE_, _PTR_, _MASK_, _ONE_, _ASSIGN_JJ_, _ASSIGN_JI_) \
2966 do { \
2967 _MASK_(_CTYPE_ *px0 = _PTR_(x0), *px1 = _PTR_(x1)); \
2968 if (class[1] == 's') { \
2969 int *pp1_; \
2970 Matrix_Calloc(pp1_, n, int); \
2971 Matrix_memcpy(pp1_, pp1 - 1, n, sizeof(int)); \
2972 for (j = 0, k = 0; j < n; ++j) { \
2973 kend = pp0[j]; \
2974 while (k < kend) { \
2975 if (pi0[k] == j) { \
2976 pi1[pp1_[j]] = pi0[k]; \
2977 _MASK_(_ASSIGN_JJ_(px1[pp1_[j]], px0[k])); \
2978 ++pp1_[j]; \
2979 } else { \
2980 pi1[pp1_[j]] = pi0[k]; \
2981 _MASK_(px1[pp1_[j]] = px0[k]); \
2982 ++pp1_[j]; \
2983 pi1[pp1_[pi0[k]]] = j; \
2984 _MASK_(_ASSIGN_JI_(px1[pp1_[pi0[k]]], px0[k])); \
2985 ++pp1_[pi0[k]]; \
2986 } \
2987 ++k; \
2988 } \
2989 } \
2990 Matrix_Free(pp1_, n); \
2991 } else if (ul == ((class[2] == 'C') ? 'U' : 'L')) { \
2992 for (j = 0, k = 0; j < n; ++j) { \
2993 kend = pp0[j]; \
2994 while (k < kend) { \
2995 *(pi1++) = *(pi0++); \
2996 _MASK_(*(px1++) = *(px0++)); \
2997 ++k; \
2998 } \
2999 *(pi1++) = j; \
3000 _MASK_(*(px1++) = _ONE_); \
3001 } \
3002 } else { \
3003 for (j = 0, k = 0; j < n; ++j) { \
3004 kend = pp0[j]; \
3005 *(pi1++) = j; \
3006 _MASK_(*(px1++) = _ONE_); \
3007 while (k < kend) { \
3008 *(pi1++) = *(pi0++); \
3009 _MASK_(*(px1++) = *(px0++)); \
3010 ++k; \
3011 } \
3012 } \
3013 } \
3014 } while (0)
3015
3016 if (class[0] == 'n')
3017 SAG_SUBCASES(int, LOGICAL, HIDE, 1, ASSIGN_REAL, ASSIGN_REAL);
3018 else {
3019 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym)),
3020 x1 = PROTECT(allocVector(TYPEOF(x0), pp1[n - 1]));
3021 SET_SLOT(to, Matrix_xSym, x1);
3022 SAG_CASES;
3023 UNPROTECT(2); /* x1, x0 */
3024 }
3025
3026 } else {
3027
3028 SEXP i0 = PROTECT(GET_SLOT(from, Matrix_iSym)),
3029 j0 = PROTECT(GET_SLOT(from, Matrix_jSym));
3030 int *pi0 = INTEGER(i0), *pj0 = INTEGER(j0);
3031 R_xlen_t nnz0 = XLENGTH(i0), nnz1;
3032
3033 if (class[1] == 's') {
3034 nnz1 = nnz0;
3035 for (R_xlen_t k = 0; k < nnz0; ++k)
3036 if (pi0[k] == pj0[k])
3037 --nnz1;
3038 } else
3039 nnz1 = n;
3040 if (nnz1 > R_XLEN_T_MAX - nnz0)
3041 error(_("attempt to allocate vector of length exceeding %s"),
3042 "R_XLEN_T_MAX");
3043 nnz1 += nnz0;
3044
3045 SEXP i1 = PROTECT(allocVector(INTSXP, nnz1)),
3046 j1 = PROTECT(allocVector(INTSXP, nnz1));
3047 int *pi1 = INTEGER(i1), *pj1 = INTEGER(j1);
3048 SET_SLOT(to, Matrix_iSym, i1);
3049 SET_SLOT(to, Matrix_jSym, j1);
3050
3051#undef SAG_SUBCASES
3052#define SAG_SUBCASES(_CTYPE_, _PTR_, _MASK_, _ONE_, _ASSIGN_JJ_, _ASSIGN_JI_) \
3053 do { \
3054 _MASK_(_CTYPE_ *px0 = _PTR_(x0), *px1 = _PTR_(x1)); \
3055 if (class[1] == 's') { \
3056 for (R_xlen_t k = 0; k < nnz0; ++k) { \
3057 if (*pi0 == *pj0) { \
3058 *(pi1++) = *pi0; \
3059 *(pj1++) = *pj0; \
3060 _MASK_(_ASSIGN_JJ_((*px1), (*px0))); \
3061 _MASK_(++px1); \
3062 } else { \
3063 *(pi1++) = *pi0; \
3064 *(pj1++) = *pj0; \
3065 _MASK_(*px1 = *px0); \
3066 _MASK_(++px1); \
3067 *(pi1++) = *pj0; \
3068 *(pj1++) = *pi0; \
3069 _MASK_(_ASSIGN_JI_((*px1), (*px0))); \
3070 _MASK_(++px1); \
3071 } \
3072 ++pi0; ++pj0; _MASK_(++px0); \
3073 } \
3074 } else { \
3075 Matrix_memcpy(pi1, pi0, nnz0, sizeof(int)); \
3076 Matrix_memcpy(pj1, pj0, nnz0, sizeof(int)); \
3077 _MASK_(Matrix_memcpy(px1, px0, nnz0, sizeof(_CTYPE_))); \
3078 pi1 += nnz0; \
3079 pj1 += nnz0; \
3080 _MASK_(px1 += nnz0); \
3081 for (int d = 0; d < n; ++d) { \
3082 *(pi1++) = *(pj1++) = d; \
3083 _MASK_(*(px1++) = _ONE_); \
3084 } \
3085 } \
3086 } while (0)
3087
3088 if (class[0] == 'n')
3089 SAG_SUBCASES(int, LOGICAL, HIDE, 1, ASSIGN_REAL, ASSIGN_REAL);
3090 else {
3091 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym)),
3092 x1 = PROTECT(allocVector(TYPEOF(x0), nnz1));
3093 SET_SLOT(to, Matrix_xSym, x1);
3094 SAG_CASES;
3095 UNPROTECT(2); /* x1, x0 */
3096 }
3097
3098 }
3099
3100#undef SAG_CASES
3101#undef SAG_SUBCASES
3102
3103 UNPROTECT(5);
3104 return to;
3105}
3106
3107/* as(<[CRT]sparseMatrix>, "generalMatrix") */
3108SEXP R_sparse_as_general(SEXP from)
3109{
3110 static const char *valid[] = {
3111 VALID_CSPARSE, VALID_RSPARSE, VALID_TSPARSE, "" };
3112 int ivalid = R_check_class_etc(from, valid);
3113 if (ivalid < 0)
3114 ERROR_INVALID_CLASS(from, __func__);
3115
3116 return sparse_as_general(from, valid[ivalid]);
3117}
3118
3119SEXP dense_as_unpacked(SEXP from, const char *class)
3120{
3121 if (class[2] != 'p')
3122 return from;
3123
3124 char cl[] = "...Matrix";
3125 cl[0] = class[0];
3126 cl[1] = class[1];
3127 switch (class[1]) {
3128 case 's':
3129 cl[2] = 'y';
3130 break;
3131 case 'p':
3132 cl[2] = 'o';
3133 break;
3134 case 'o':
3135 case 't':
3136 cl[2] = 'r';
3137 break;
3138 default:
3139 break;
3140 }
3141 SEXP to = PROTECT(newObject(cl));
3142
3143 SEXP dim = PROTECT(GET_SLOT(from, Matrix_DimSym));
3144 int n = INTEGER(dim)[0];
3145 if ((Matrix_int_fast64_t) n * n > R_XLEN_T_MAX)
3146 error(_("attempt to allocate vector of length exceeding %s"),
3147 "R_XLEN_T_MAX");
3148 if (n > 0)
3149 SET_SLOT(to, Matrix_DimSym, dim);
3150 UNPROTECT(1); /* dim */
3151
3152 SEXP dimnames = PROTECT(GET_SLOT(from, Matrix_DimNamesSym));
3153 SET_SLOT(to, Matrix_DimNamesSym, dimnames);
3154 UNPROTECT(1); /* dimnames */
3155
3156 SEXP uplo = PROTECT(GET_SLOT(from, Matrix_uploSym));
3157 char ul = *CHAR(STRING_ELT(uplo, 0));
3158 if (ul != 'U')
3159 SET_SLOT(to, Matrix_uploSym, uplo);
3160 UNPROTECT(1); /* uplo */
3161
3162 if (cl[1] != 't') {
3163 SEXP factors = PROTECT(GET_SLOT(from, Matrix_factorsSym));
3164 if (LENGTH(factors) > 0)
3165 SET_SLOT(to, Matrix_factorsSym, factors);
3166 UNPROTECT(1); /* factors */
3167
3168 if (cl[1] == 'o') {
3169 SEXP sd = PROTECT(GET_SLOT(from, Matrix_sdSym));
3170 if (LENGTH(sd) > 0)
3171 SET_SLOT(to, Matrix_sdSym, sd);
3172 UNPROTECT(1); /* sd */
3173 }
3174 } else {
3175 SEXP diag = PROTECT(GET_SLOT(from, Matrix_diagSym));
3176 char di = *CHAR(STRING_ELT(diag, 0));
3177 if (di != 'N')
3178 SET_SLOT(to, Matrix_diagSym, diag);
3179 UNPROTECT(1); /* diag */
3180 }
3181
3182 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym)),
3183 x1 = PROTECT(allocVector(TYPEOF(x0), (R_xlen_t) n * n));
3184 SET_SLOT(to, Matrix_xSym, x1);
3185
3186#define UNPACK(_PREFIX_, _CTYPE_, _PTR_) \
3187 do { \
3188 _CTYPE_ *px0 = _PTR_(x0), *px1 = _PTR_(x1); \
3189 Matrix_memset(px1, 0, (R_xlen_t) n * n, sizeof(_CTYPE_)); \
3190 _PREFIX_ ## unpack1(px1, px0, n, ul, 'N'); \
3191 } while (0)
3192
3193 switch (cl[0]) {
3194 case 'n':
3195 case 'l':
3196 UNPACK(i, int, LOGICAL);
3197 break;
3198 case 'i':
3199 UNPACK(i, int, INTEGER);
3200 break;
3201 case 'c':
3202 case 'd':
3203 UNPACK(d, double, REAL);
3204 break;
3205 case 'z':
3206 UNPACK(z, Rcomplex, COMPLEX);
3207 break;
3208 default:
3209 break;
3210 }
3211
3212#undef UNPACK
3213
3214 UNPROTECT(3); /* x1, x0, to */
3215 return to;
3216}
3217
3218/* as(<denseMatrix>, "unpackedMatrix") */
3219SEXP R_dense_as_unpacked(SEXP from)
3220{
3221 static const char *valid[] = {
3222 "dpoMatrix", "dppMatrix", "corMatrix", "copMatrix",
3223 VALID_DENSE, "" };
3224 int ivalid = R_check_class_etc(from, valid);
3225 if (ivalid < 0)
3226 ERROR_INVALID_CLASS(from, __func__);
3227
3228 return dense_as_unpacked(from, valid[ivalid]);
3229}
3230
3231SEXP dense_as_packed(SEXP from, const char *class, char ul, char di)
3232{
3233 if (class[2] == 'p')
3234 return from;
3235 int ge = class[1] == 'g';
3236
3237 char cl[] = "...Matrix";
3238 cl[0] = class[0];
3239 cl[1] = (ge) ? ((di == '\0') ? 's' : 't') : class[1];
3240 cl[2] = 'p';
3241 SEXP to = PROTECT(newObject(cl));
3242
3243 SEXP dim = PROTECT(GET_SLOT(from, Matrix_DimSym));
3244 int *pdim = INTEGER(dim), n = pdim[0];
3245 if (pdim[1] != n)
3246 error(_("attempt to pack non-square matrix"));
3247 if (n > 0)
3248 SET_SLOT(to, Matrix_DimSym, dim);
3249 UNPROTECT(1); /* dim */
3250
3251 SEXP dimnames = PROTECT(GET_SLOT(from, Matrix_DimNamesSym));
3252 SET_SLOT(to, Matrix_DimNamesSym, dimnames);
3253 UNPROTECT(1); /* dimnames */
3254
3255 if (ge) {
3256 if (ul != 'U') {
3257 SEXP uplo = PROTECT(mkString("L"));
3258 SET_SLOT(to, Matrix_uploSym, uplo);
3259 UNPROTECT(1); /* uplo */
3260 }
3261 if (cl[1] == 't' && di != 'N') {
3262 SEXP diag = PROTECT(mkString("U"));
3263 SET_SLOT(to, Matrix_diagSym, diag);
3264 UNPROTECT(1); /* diag */
3265 }
3266 } else {
3267 SEXP uplo = PROTECT(GET_SLOT(from, Matrix_uploSym));
3268 ul = *CHAR(STRING_ELT(uplo, 0));
3269 if (ul != 'U')
3270 SET_SLOT(to, Matrix_uploSym, uplo);
3271 UNPROTECT(1); /* uplo */
3272
3273 if (cl[1] == 't') {
3274 SEXP diag = PROTECT(GET_SLOT(from, Matrix_diagSym));
3275 di = *CHAR(STRING_ELT(diag, 0));
3276 if (di != 'N')
3277 SET_SLOT(to, Matrix_diagSym, diag);
3278 UNPROTECT(1); /* diag */
3279 } else {
3280 SEXP factors = PROTECT(GET_SLOT(from, Matrix_factorsSym));
3281 if (LENGTH(factors) > 0)
3282 SET_SLOT(to, Matrix_factorsSym, factors);
3283 UNPROTECT(1); /* factors */
3284
3285 if (cl[1] == 'o') {
3286 SEXP sd = PROTECT(GET_SLOT(from, Matrix_sdSym));
3287 if (LENGTH(sd) > 0)
3288 SET_SLOT(to, Matrix_sdSym, sd);
3289 UNPROTECT(1); /* sd */
3290 }
3291 }
3292 }
3293
3294 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym)),
3295 x1 = PROTECT(allocVector(TYPEOF(x0), PACKED_LENGTH(n)));
3296 SET_SLOT(to, Matrix_xSym, x1);
3297
3298#define PACK(_PREFIX_, _CTYPE_, _PTR_) \
3299 do { \
3300 _CTYPE_ *px0 = _PTR_(x0), *px1 = _PTR_(x1); \
3301 _PREFIX_ ## pack2(px1, px0, n, ul, 'N'); \
3302 } while (0)
3303
3304 switch (cl[0]) {
3305 case 'n':
3306 case 'l':
3307 PACK(i, int, LOGICAL);
3308 break;
3309 case 'i':
3310 PACK(i, int, INTEGER);
3311 break;
3312 case 'c':
3313 case 'd':
3314 PACK(d, double, REAL);
3315 break;
3316 case 'z':
3317 PACK(z, Rcomplex, COMPLEX);
3318 break;
3319 default:
3320 break;
3321 }
3322
3323#undef PACK
3324
3325 UNPROTECT(3); /* x1, x0, to */
3326 return to;
3327}
3328
3329/* as(<denseMatrix>, "packedMatrix") */
3330SEXP R_dense_as_packed(SEXP from, SEXP uplo, SEXP diag)
3331{
3332 static const char *valid[] = {
3333 "dpoMatrix", "dppMatrix", "corMatrix", "copMatrix",
3334 VALID_DENSE, "" };
3335 int ivalid = R_check_class_etc(from, valid);
3336 if (ivalid < 0)
3337 ERROR_INVALID_CLASS(from, __func__);
3338
3339 char ul = 'U', di = '\0';
3340 if (valid[ivalid][1] == 'g') {
3341 if (TYPEOF(uplo) != STRSXP || LENGTH(uplo) < 1 ||
3342 (uplo = STRING_ELT(uplo, 0)) == NA_STRING ||
3343 ((ul = *CHAR(uplo)) != 'U' && ul != 'L'))
3344 error(_("'%s' must be \"%s\" or \"%s\""), "uplo", "U", "L");
3345 if (diag != R_NilValue) {
3346 if (TYPEOF(diag) != STRSXP || LENGTH(diag) < 1 ||
3347 (diag = STRING_ELT(diag, 0)) == NA_STRING ||
3348 ((di = *CHAR(diag)) != 'N' && ul != 'U'))
3349 error(_("'%s' must be \"%s\" or \"%s\""), "diag", "N", "U");
3350 }
3351 }
3352
3353 return dense_as_packed(from, valid[ivalid], ul, di);
3354}
3355
3356void trans(SEXP p0, SEXP i0, SEXP x0, SEXP p1, SEXP i1, SEXP x1,
3357 int m, int n)
3358{
3359 int *pp0 = INTEGER(p0), *pp1 = INTEGER(p1),
3360 *pi0 = INTEGER(i0), *pi1 = INTEGER(i1),
3361 i, j, k, kend, nnz = pp0[n];
3362 Matrix_memset(pp1, 0, (R_xlen_t) m + 1, sizeof(int));
3363 ++pp0; ++pp1;
3364 for (k = 0; k < nnz; ++k)
3365 ++pp1[pi0[k]];
3366 for (i = 0; i < m; ++i)
3367 pp1[i] += pp1[i - 1];
3368
3369 int *pp1_;
3370 Matrix_Calloc(pp1_, m, int);
3371 Matrix_memcpy(pp1_, pp1 - 1, m, sizeof(int));
3372
3373#define TRANS_LOOP(_CTYPE_, _PTR_, _MASK_) \
3374 do { \
3375 _MASK_(_CTYPE_ *px0 = _PTR_(x0), *px1 = _PTR_(x1)); \
3376 for (j = 0, k = 0; j < n; ++j) { \
3377 kend = pp0[j]; \
3378 while (k < kend) { \
3379 i = pi0[k]; \
3380 pi1[pp1_[i]] = j; \
3381 _MASK_(px1[pp1_[i]] = px0[k]); \
3382 ++pp1_[i]; \
3383 ++k; \
3384 } \
3385 } \
3386 } while (0)
3387
3388 if (!(x0 && x1) || TYPEOF(x0) != TYPEOF(x1))
3389 TRANS_LOOP(int, LOGICAL, HIDE);
3390 else {
3391 switch (TYPEOF(x0)) {
3392 case LGLSXP:
3393 TRANS_LOOP(int, LOGICAL, SHOW);
3394 break;
3395 case INTSXP:
3396 TRANS_LOOP(int, INTEGER, SHOW);
3397 break;
3398 case REALSXP:
3399 TRANS_LOOP(double, REAL, SHOW);
3400 break;
3401 case CPLXSXP:
3402 TRANS_LOOP(Rcomplex, COMPLEX, SHOW);
3403 break;
3404 default:
3405 break;
3406 }
3407 }
3408
3409#undef TRANS_LOOP
3410
3411 Matrix_Free(pp1_, m);
3412 return;
3413}
3414
3415void tsort(SEXP i0, SEXP j0, SEXP x0, SEXP *p1, SEXP *i1, SEXP *x1,
3416 int m, int n)
3417{
3418 R_xlen_t nnz0 = XLENGTH(i0), nnz1 = 0;
3419 if (nnz0 > INT_MAX)
3420 error(_("unable to aggregate %s with '%s' and '%s' slots of length exceeding %s"),
3421 "TsparseMatrix", "i", "j", "2^31-1");
3422
3423 /* FIXME: test for overflow and throw error only in that case */
3424
3425 PROTECT(*p1 = allocVector(INTSXP, (R_xlen_t) n + 1));
3426 int *pi0 = INTEGER(i0), *pj0 = INTEGER(j0), *pp1 = INTEGER(*p1), *pi1,
3427 i, j, r = (m < n) ? n : m;
3428 R_xlen_t k, kstart, kend, kend_;
3429 *(pp1++) = 0;
3430
3431 int *workA, *workB, *workC, *pj_;
3432 size_t lwork = (size_t) m + r + m + nnz0;
3433 Matrix_Calloc(workA, lwork, int);
3434 workB = workA + m;
3435 workC = workB + r;
3436 pj_ = workC + m;
3437
3438#define TSORT_LOOP(_CTYPE_, _PTR_, _MASK_, _INCREMENT_) \
3439 do { \
3440 _MASK_(_CTYPE_ *px0 = _PTR_(x0), *px1, *px_); \
3441 _MASK_(Matrix_Calloc(px_, nnz0, _CTYPE_)); \
3442 \
3443 /* 1. Tabulate column indices in workA[i] */ \
3444 \
3445 for (k = 0; k < nnz0; ++k) \
3446 ++workA[pi0[k]]; \
3447 \
3448 /* 2. Compute cumulative sum in workA[i], copying to workB[i] */ \
3449 /* */ \
3450 /* workA[i]: number of column indices listed for row i, */ \
3451 /* incl. duplicates */ \
3452 \
3453 for (i = 1; i < m; ++i) \
3454 workA[i] += (workB[i] = workA[i - 1]); \
3455 \
3456 /* 3. Group column indices and data by row in pj_[k], px_[k] */ \
3457 /* */ \
3458 /* workA[i]: number of column indices listed for row <= i, */ \
3459 /* incl. duplicates */ \
3460 /* workB[i]: number of column indices listed for row < i, */ \
3461 /* incl. duplicates */ \
3462 \
3463 for (k = 0; k < nnz0; ++k) { \
3464 pj_[workB[pi0[k]]] = pj0[k] ; \
3465 _MASK_(px_[workB[pi0[k]]] = px0[k]); \
3466 ++workB[pi0[k]]; \
3467 } \
3468 \
3469 /* 4. Gather _unique_ column indices at the front of each group, */ \
3470 /* aggregating data accordingly; record in workC[i] where */ \
3471 /* the unique column indices stop and the duplicates begin */ \
3472 /* */ \
3473 /* workB[.]: unused */ \
3474 /* pj_[k]: column indices grouped by row, */ \
3475 /* incl. duplicates, unsorted */ \
3476 /* px_[k]: corresponding data */ \
3477 \
3478 k = 0; \
3479 for (j = 0; j < n; ++j) \
3480 workB[j] = -1; \
3481 for (i = 0; i < m; ++i) { \
3482 kstart = kend_ = k; \
3483 kend = workA[i]; \
3484 while (k < kend) { \
3485 if (workB[pj_[k]] < kstart) { \
3486 /* Have not yet seen this column index */ \
3487 workB[pj_[k]] = kend_; \
3488 pj_[kend_] = pj_[k] ; \
3489 _MASK_(px_[kend_] = px_[k]); \
3490 ++kend_; \
3491 } else { \
3492 /* Have already seen this column index */ \
3493 _MASK_(_INCREMENT_(px_[workB[pj_[k]]], px_[k])); \
3494 } \
3495 ++k; \
3496 } \
3497 workC[i] = kend_; \
3498 nnz1 += kend_ - kstart; \
3499 } \
3500 \
3501 /* 5. Tabulate _unique_ column indices in workB[j] */ \
3502 /* */ \
3503 /* workC[i]: pointer to first non-unique column index in row i */ \
3504 /* pi_[k]: column indices grouped by row, */ \
3505 /* with unique indices in front, */ \
3506 /* i.e., in positions workA[i - 1] <= k < workC[i] */ \
3507 /* px_[k]: corresponding data, "cumulated" appropriately */ \
3508 \
3509 k = 0; \
3510 Matrix_memset(workB, 0, n, sizeof(int)); \
3511 for (i = 0; i < m; ++i) { \
3512 kend_ = workC[i]; \
3513 while (k < kend_) { \
3514 ++workB[pj_[k]]; \
3515 ++k; \
3516 } \
3517 k = workA[i]; \
3518 } \
3519 \
3520 /* 6. Compute cumulative sum in pp1[j], copying to workB[j] */ \
3521 /* */ \
3522 /* workB[j]: number of nonzero elements in column j */ \
3523 \
3524 for (j = 0; j < n; ++j) { \
3525 pp1[j] = pp1[j - 1] + workB[j]; \
3526 workB[j] = pp1[j - 1]; \
3527 } \
3528 \
3529 PROTECT(*i1 = allocVector( INTSXP, nnz1)) ; \
3530 _MASK_(PROTECT(*x1 = allocVector(TYPEOF(x0), nnz1))); \
3531 pi1 = INTEGER(*i1) ; \
3532 _MASK_(px1 = _PTR_(*x1)); \
3533 \
3534 /* 7. Pop unique (i,j) pairs from the unsorted stacks 0 <= i < m */ \
3535 /* onto new stacks 0 <= j < n, which will be sorted */ \
3536 /* */ \
3537 /* workB[j]: number of nonzero elements in columns < j */ \
3538 /* pp1[j]: number of nonzero elements in columns <= j */ \
3539 \
3540 k = 0; \
3541 for (i = 0; i < m; ++i) { \
3542 kend_ = workC[i]; \
3543 while (k < kend_) { \
3544 pi1[workB[pj_[k]]] = i; \
3545 _MASK_(px1[workB[pj_[k]]] = px_[k]); \
3546 ++workB[pj_[k]]; \
3547 ++k; \
3548 } \
3549 k = workA[i]; \
3550 } \
3551 \
3552 _MASK_(Matrix_Free(px_, nnz0)); \
3553 _MASK_(UNPROTECT(1)); /* *px1 */ \
3554 UNPROTECT(1) ; /* *pi1 */ \
3555 } while (0)
3556
3557 if (!x0)
3558 TSORT_LOOP(int, LOGICAL, HIDE, INCREMENT_PATTERN);
3559 else {
3560 switch (TYPEOF(x0)) {
3561 case LGLSXP:
3562 TSORT_LOOP(int, LOGICAL, SHOW, INCREMENT_LOGICAL);
3563 break;
3564 case INTSXP:
3565 TSORT_LOOP(int, INTEGER, SHOW, INCREMENT_INTEGER);
3566 break;
3567 case REALSXP:
3568 TSORT_LOOP(double, REAL, SHOW, INCREMENT_REAL);
3569 break;
3570 case CPLXSXP:
3571 TSORT_LOOP(Rcomplex, COMPLEX, SHOW, INCREMENT_COMPLEX);
3572 break;
3573 default:
3574 break;
3575 }
3576 }
3577
3578#undef TSORT_LOOP
3579
3580 Matrix_Free(workA, lwork);
3581 UNPROTECT(1); /* *p1 */
3582 return;
3583}
3584
3585void taggr(SEXP i0, SEXP j0, SEXP x0, SEXP *i1, SEXP *j1, SEXP *x1,
3586 int m, int n)
3587{
3588 R_xlen_t nnz0 = XLENGTH(i0), nnz1 = 0;
3589 if (nnz0 > INT_MAX)
3590 error(_("unable to aggregate %s with '%s' and '%s' slots of length exceeding %s"),
3591 "TsparseMatrix", "i", "j", "2^31-1");
3592
3593 /* FIXME: test for overflow and throw error only in that case */
3594
3595 int *pi0 = INTEGER(i0), *pj0 = INTEGER(j0), *pi1, *pj1,
3596 i, j, r = (m < n) ? n : m;
3597 R_xlen_t k, kstart, kend, kend_;
3598
3599 int *workA, *workB, *workC, *pj_;
3600 size_t lwork = (size_t) m + r + m + nnz0;
3601 Matrix_Calloc(workA, lwork, int);
3602 workB = workA + m;
3603 workC = workB + r;
3604 pj_ = workC + m;
3605
3606#define TAGGR_LOOP(_CTYPE_, _PTR_, _MASK_, _INCREMENT_) \
3607 do { \
3608 _MASK_(_CTYPE_ *px0 = _PTR_(x0), *px1, *px_); \
3609 _MASK_(Matrix_Calloc(px_, nnz0, _CTYPE_)); \
3610 \
3611 /* 1. Tabulate column indices in workA[i] */ \
3612 \
3613 for (k = 0; k < nnz0; ++k) \
3614 ++workA[pi0[k]]; \
3615 \
3616 /* 2. Compute cumulative sum in workA[i], copying to workB[i] */ \
3617 /* */ \
3618 /* workA[i]: number of column indices listed for row i, */ \
3619 /* incl. duplicates */ \
3620 \
3621 for (i = 1; i < m; ++i) \
3622 workA[i] += (workB[i] = workA[i - 1]); \
3623 \
3624 /* 3. Group column indices and data by row in pj_[k], px_[k] */ \
3625 /* */ \
3626 /* workA[i]: number of column indices listed for row <= i, */ \
3627 /* incl. duplicates */ \
3628 /* workB[i]: number of column indices listed for row < i, */ \
3629 /* incl. duplicates */ \
3630 \
3631 for (k = 0; k < nnz0; ++k) { \
3632 pj_[workB[pi0[k]]] = pj0[k] ; \
3633 _MASK_(px_[workB[pi0[k]]] = px0[k]); \
3634 ++workB[pi0[k]]; \
3635 } \
3636 \
3637 /* 4. Gather _unique_ column indices at the front of each group, */ \
3638 /* aggregating data accordingly; record in workC[i] where */ \
3639 /* the unique column indices stop and the duplicates begin */ \
3640 /* */ \
3641 /* workB[.]: unused */ \
3642 /* pj_[k]: column indices grouped by row, */ \
3643 /* incl. duplicates, unsorted */ \
3644 /* px_[k]: corresponding data */ \
3645 \
3646 k = 0; \
3647 for (j = 0; j < n; ++j) \
3648 workB[j] = -1; \
3649 for (i = 0; i < m; ++i) { \
3650 kstart = kend_ = k; \
3651 kend = workA[i]; \
3652 while (k < kend) { \
3653 if (workB[pj_[k]] < kstart) { \
3654 /* Have not yet seen this column index */ \
3655 workB[pj_[k]] = kend_; \
3656 pj_[kend_] = pj_[k] ; \
3657 _MASK_(px_[kend_] = px_[k]); \
3658 ++kend_; \
3659 } else { \
3660 /* Have already seen this column index */ \
3661 _MASK_(_INCREMENT_(px_[workB[pj_[k]]], px_[k])); \
3662 } \
3663 ++k; \
3664 } \
3665 workC[i] = kend_; \
3666 nnz1 += kend_ - kstart; \
3667 } \
3668 if (nnz1 != nnz0) { \
3669 PROTECT(*i1 = allocVector( INTSXP, nnz1)) ; \
3670 PROTECT(*j1 = allocVector( INTSXP, nnz1)) ; \
3671 _MASK_(PROTECT(*x1 = allocVector(TYPEOF(x0), nnz1))); \
3672 pi1 = INTEGER(*i1) ; \
3673 pj1 = INTEGER(*j1) ; \
3674 _MASK_(px1 = _PTR_(*x1)); \
3675 \
3676 k = 0; \
3677 for (i = 0; i < m; ++i) { \
3678 kend_ = workC[i]; \
3679 while (k < kend_) { \
3680 *(pi1++) = i ; \
3681 *(pj1++) = pj_[k] ; \
3682 _MASK_(*(px1++) = px_[k]); \
3683 ++k; \
3684 } \
3685 k = workA[i]; \
3686 } \
3687 \
3688 _MASK_(UNPROTECT(1)); /* *px1 */ \
3689 UNPROTECT(2) ; /* *pj1, *px1 */ \
3690 } \
3691 _MASK_(Matrix_Free(px_, nnz0)); \
3692 } while (0)
3693
3694 if (!x0)
3695 TAGGR_LOOP(int, LOGICAL, HIDE, );
3696 else {
3697 switch (TYPEOF(x0)) {
3698 case LGLSXP:
3699 TAGGR_LOOP(int, LOGICAL, SHOW, INCREMENT_LOGICAL);
3700 break;
3701 case INTSXP:
3702 TAGGR_LOOP(int, INTEGER, SHOW, INCREMENT_INTEGER);
3703 break;
3704 case REALSXP:
3705 TAGGR_LOOP(double, REAL, SHOW, INCREMENT_REAL);
3706 break;
3707 case CPLXSXP:
3708 TAGGR_LOOP(Rcomplex, COMPLEX, SHOW, INCREMENT_COMPLEX);
3709 break;
3710 default:
3711 break;
3712 }
3713 }
3714
3715#undef TAGGR_LOOP
3716
3717 Matrix_Free(workA, lwork);
3718 return;
3719}
3720
3721SEXP sparse_as_Csparse(SEXP from, const char *class)
3722{
3723 if (class[2] == 'C')
3724 return from;
3725
3726 char cl[] = "..CMatrix";
3727 cl[0] = class[0];
3728 cl[1] = class[1];
3729 SEXP to = PROTECT(newObject(cl));
3730
3731 SEXP dim = PROTECT(GET_SLOT(from, Matrix_DimSym));
3732 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1];
3733 if (m != n || n > 0)
3734 SET_SLOT(to, Matrix_DimSym, dim);
3735 UNPROTECT(1); /* dim */
3736
3737 SEXP dimnames = PROTECT(GET_SLOT(from, Matrix_DimNamesSym));
3738 SET_SLOT(to, Matrix_DimNamesSym, dimnames);
3739 UNPROTECT(1); /* dimnames */
3740
3741 if (class[1] != 'g') {
3742 SEXP uplo = PROTECT(GET_SLOT(from, Matrix_uploSym));
3743 char ul = *CHAR(STRING_ELT(uplo, 0));
3744 if (ul != 'U')
3745 SET_SLOT(to, Matrix_uploSym, uplo);
3746 UNPROTECT(1); /* uplo */
3747 }
3748 if (class[1] == 't') {
3749 SEXP diag = PROTECT(GET_SLOT(from, Matrix_diagSym));
3750 char di = *CHAR(STRING_ELT(diag, 0));
3751 if (di != 'N')
3752 SET_SLOT(to, Matrix_diagSym, diag);
3753 UNPROTECT(1); /* diag */
3754 } else {
3755 SEXP factors = PROTECT(GET_SLOT(from, Matrix_factorsSym));
3756 if (LENGTH(factors) > 0)
3757 SET_SLOT(to, Matrix_factorsSym, factors);
3758 UNPROTECT(1); /* factors */
3759 }
3760
3761 if (class[2] == 'R') {
3762 SEXP p0 = PROTECT(GET_SLOT(from, Matrix_pSym)),
3763 j0 = PROTECT(GET_SLOT(from, Matrix_jSym)),
3764 p1 = PROTECT(allocVector(INTSXP, (R_xlen_t) n + 1)),
3765 i1 = PROTECT(allocVector(INTSXP, INTEGER(p0)[m]));
3766 SET_SLOT(to, Matrix_pSym, p1);
3767 SET_SLOT(to, Matrix_iSym, i1);
3768 if (class[0] == 'n')
3769 trans(p0, j0, NULL, p1, i1, NULL, n, m);
3770 else {
3771 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym)),
3772 x1 = PROTECT(allocVector(TYPEOF(x0), INTEGER(p0)[m]));
3773 SET_SLOT(to, Matrix_xSym, x1);
3774 trans(p0, j0, x0, p1, i1, x1, n, m);
3775 UNPROTECT(2); /* x1, x0 */
3776 }
3777 UNPROTECT(4); /* i1, p1, j0, p0 */
3778 } else {
3779 SEXP i0 = PROTECT(GET_SLOT(from, Matrix_iSym)),
3780 j0 = PROTECT(GET_SLOT(from, Matrix_jSym)),
3781 p1 = NULL, i1 = NULL;
3782 if (class[0] == 'n') {
3783 tsort(i0, j0, NULL, &p1, &i1, NULL, m, n);
3784 PROTECT(p1);
3785 PROTECT(i1);
3786 SET_SLOT(to, Matrix_pSym, p1);
3787 SET_SLOT(to, Matrix_iSym, i1);
3788 UNPROTECT(2); /* i1, p1 */
3789 } else {
3790 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym)),
3791 x1 = NULL;
3792 tsort(i0, j0, x0, &p1, &i1, &x1, m, n);
3793 PROTECT(p1);
3794 PROTECT(i1);
3795 PROTECT(x1);
3796 SET_SLOT(to, Matrix_pSym, p1);
3797 SET_SLOT(to, Matrix_iSym, i1);
3798 SET_SLOT(to, Matrix_xSym, x1);
3799 UNPROTECT(4); /* x1, i1, p1, x0 */
3800 }
3801 UNPROTECT(2); /* j0, i0 */
3802 }
3803
3804 UNPROTECT(1); /* to */
3805 return to;
3806}
3807
3808/* as(<[CRT]sparseMatrix>, "CsparseMatrix") */
3809SEXP R_sparse_as_Csparse(SEXP from)
3810{
3811 static const char *valid[] = {
3812 VALID_CSPARSE, VALID_RSPARSE, VALID_TSPARSE, "" };
3813 int ivalid = R_check_class_etc(from, valid);
3814 if (ivalid < 0)
3815 ERROR_INVALID_CLASS(from, __func__);
3816
3817 return sparse_as_Csparse(from, valid[ivalid]);
3818}
3819
3820SEXP sparse_as_Rsparse(SEXP from, const char *class)
3821{
3822 if (class[2] == 'R')
3823 return from;
3824
3825 char cl[] = "..RMatrix";
3826 cl[0] = class[0];
3827 cl[1] = class[1];
3828 SEXP to = PROTECT(newObject(cl));
3829
3830 SEXP dim = PROTECT(GET_SLOT(from, Matrix_DimSym));
3831 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1];
3832 if (m != n || n > 0)
3833 SET_SLOT(to, Matrix_DimSym, dim);
3834 UNPROTECT(1); /* dim */
3835
3836 SEXP dimnames = PROTECT(GET_SLOT(from, Matrix_DimNamesSym));
3837 SET_SLOT(to, Matrix_DimNamesSym, dimnames);
3838 UNPROTECT(1); /* dimnames */
3839
3840 if (class[1] != 'g') {
3841 SEXP uplo = PROTECT(GET_SLOT(from, Matrix_uploSym));
3842 char ul = *CHAR(STRING_ELT(uplo, 0));
3843 if (ul != 'U')
3844 SET_SLOT(to, Matrix_uploSym, uplo);
3845 UNPROTECT(1); /* uplo */
3846 }
3847 if (class[1] == 't') {
3848 SEXP diag = PROTECT(GET_SLOT(from, Matrix_diagSym));
3849 char di = *CHAR(STRING_ELT(diag, 0));
3850 if (di != 'N')
3851 SET_SLOT(to, Matrix_diagSym, diag);
3852 UNPROTECT(1); /* diag */
3853 } else {
3854 SEXP factors = PROTECT(GET_SLOT(from, Matrix_factorsSym));
3855 if (LENGTH(factors) > 0)
3856 SET_SLOT(to, Matrix_factorsSym, factors);
3857 UNPROTECT(1); /* factors */
3858 }
3859
3860 if (class[2] == 'C') {
3861 SEXP p0 = PROTECT(GET_SLOT(from, Matrix_pSym)),
3862 i0 = PROTECT(GET_SLOT(from, Matrix_iSym)),
3863 p1 = PROTECT(allocVector(INTSXP, (R_xlen_t) m + 1)),
3864 j1 = PROTECT(allocVector(INTSXP, INTEGER(p0)[n]));
3865 SET_SLOT(to, Matrix_pSym, p1);
3866 SET_SLOT(to, Matrix_jSym, j1);
3867 if (class[0] == 'n')
3868 trans(p0, i0, NULL, p1, j1, NULL, m, n);
3869 else {
3870 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym)),
3871 x1 = PROTECT(allocVector(TYPEOF(x0), INTEGER(p0)[n]));
3872 SET_SLOT(to, Matrix_xSym, x1);
3873 trans(p0, i0, x0, p1, j1, x1, m, n);
3874 UNPROTECT(2); /* x1, x0 */
3875 }
3876 UNPROTECT(4); /* j1, p1, i0, p0 */
3877 } else {
3878 SEXP i0 = PROTECT(GET_SLOT(from, Matrix_iSym)),
3879 j0 = PROTECT(GET_SLOT(from, Matrix_jSym)),
3880 p1 = NULL, j1 = NULL;
3881 if (class[0] == 'n') {
3882 tsort(j0, i0, NULL, &p1, &j1, NULL, n, m);
3883 PROTECT(p1);
3884 PROTECT(j1);
3885 SET_SLOT(to, Matrix_pSym, p1);
3886 SET_SLOT(to, Matrix_jSym, j1);
3887 UNPROTECT(2); /* j1, p1 */
3888 } else {
3889 SEXP x0 = PROTECT(GET_SLOT(from, Matrix_xSym)),
3890 x1 = NULL;
3891 tsort(j0, i0, x0, &p1, &j1, &x1, n, m);
3892 PROTECT(p1);
3893 PROTECT(j1);
3894 PROTECT(x1);
3895 SET_SLOT(to, Matrix_pSym, p1);
3896 SET_SLOT(to, Matrix_jSym, j1);
3897 SET_SLOT(to, Matrix_xSym, x1);
3898 UNPROTECT(4); /* x1, j1, p1, x0 */
3899 }
3900 UNPROTECT(2); /* j0, i0 */
3901 }
3902
3903 UNPROTECT(1); /* to */
3904 return to;
3905}
3906
3907/* as(<[CRT]sparseMatrix>, "RsparseMatrix") */
3908SEXP R_sparse_as_Rsparse(SEXP from)
3909{
3910 static const char *valid[] = {
3911 VALID_CSPARSE, VALID_RSPARSE, VALID_TSPARSE, "" };
3912 int ivalid = R_check_class_etc(from, valid);
3913 if (ivalid < 0)
3914 ERROR_INVALID_CLASS(from, __func__);
3915
3916 return sparse_as_Rsparse(from, valid[ivalid]);
3917}
3918
3919SEXP sparse_as_Tsparse(SEXP from, const char *class)
3920{
3921 if (class[2] == 'T')
3922 return from;
3923
3924 char cl[] = "..TMatrix";
3925 cl[0] = class[0];
3926 cl[1] = class[1];
3927 SEXP to = PROTECT(newObject(cl));
3928
3929 SEXP dim = PROTECT(GET_SLOT(from, Matrix_DimSym));
3930 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1];
3931 if (m != n || n > 0)
3932 SET_SLOT(to, Matrix_DimSym, dim);
3933 UNPROTECT(1); /* dim */
3934
3935 SEXP dimnames = PROTECT(GET_SLOT(from, Matrix_DimNamesSym));
3936 SET_SLOT(to, Matrix_DimNamesSym, dimnames);
3937 UNPROTECT(1); /* dimnames */
3938
3939 if (class[1] != 'g') {
3940 SEXP uplo = PROTECT(GET_SLOT(from, Matrix_uploSym));
3941 char ul = *CHAR(STRING_ELT(uplo, 0));
3942 if (ul != 'U')
3943 SET_SLOT(to, Matrix_uploSym, uplo);
3944 UNPROTECT(1); /* uplo */
3945 }
3946 if (class[1] == 't') {
3947 SEXP diag = PROTECT(GET_SLOT(from, Matrix_diagSym));
3948 char di = *CHAR(STRING_ELT(diag, 0));
3949 if (di != 'N')
3950 SET_SLOT(to, Matrix_diagSym, diag);
3951 UNPROTECT(1); /* diag */
3952 } else {
3953 SEXP factors = PROTECT(GET_SLOT(from, Matrix_factorsSym));
3954 if (LENGTH(factors) > 0)
3955 SET_SLOT(to, Matrix_factorsSym, factors);
3956 UNPROTECT(1); /* factors */
3957 }
3958
3959 SEXP iSym = (class[2] == 'C') ? Matrix_iSym : Matrix_jSym,
3960 jSym = (class[2] == 'C') ? Matrix_jSym : Matrix_iSym,
3961 p = PROTECT(GET_SLOT(from, Matrix_pSym)),
3962 i = PROTECT(GET_SLOT(from, iSym));
3963 int *pp = INTEGER(p), *pi, r = (class[2] == 'C') ? n : m, nnz = pp[r];
3964 if (XLENGTH(i) == nnz) {
3965 SET_SLOT(to, iSym, i);
3966 UNPROTECT(1); /* i */
3967 } else {
3968 SEXP i_ = PROTECT(allocVector(INTSXP, nnz));
3969 Matrix_memcpy(INTEGER(i_), INTEGER(i), nnz, sizeof(int));
3970 SET_SLOT(to, iSym, i_);
3971 UNPROTECT(2); /* i_, i */
3972 }
3973 PROTECT(i = allocVector(INTSXP, nnz));
3974 SET_SLOT(to, jSym, i);
3975 pi = INTEGER(i); ++pp;
3976 int j, k, kend;
3977 for (j = 0, k = 0; j < r; ++j) {
3978 kend = pp[j];
3979 while (k < kend)
3980 pi[k++] = j;
3981 }
3982 UNPROTECT(2); /* i, p */
3983
3984 if (class[0] != 'n') {
3985 SEXP x = PROTECT(GET_SLOT(from, Matrix_xSym));
3986 if (XLENGTH(x) == nnz)
3987 SET_SLOT(to, Matrix_xSym, x);
3988 else {
3989 SEXP x_ = PROTECT(allocVector(TYPEOF(x), nnz));
3990 SET_SLOT(to, Matrix_xSym, x_);
3991 switch (class[0]) {
3992 case 'l':
3993 Matrix_memcpy(LOGICAL(x_), LOGICAL(x), nnz, sizeof(int));
3994 break;
3995 case 'i':
3996 Matrix_memcpy(INTEGER(x_), INTEGER(x), nnz, sizeof(int));
3997 break;
3998 case 'd':
3999 Matrix_memcpy(REAL(x_), REAL(x), nnz, sizeof(double));
4000 break;
4001 case 'z':
4002 Matrix_memcpy(COMPLEX(x_), COMPLEX(x), nnz, sizeof(Rcomplex));
4003 break;
4004 default:
4005 break;
4006 }
4007 UNPROTECT(1); /* x_ */
4008 }
4009 UNPROTECT(1); /* x */
4010 }
4011
4012 UNPROTECT(1); /* to */
4013 return to;
4014}
4015
4016/* as(<[CRT]sparseMatrix>, "TsparseMatrix") */
4017SEXP R_sparse_as_Tsparse(SEXP from)
4018{
4019 static const char *valid[] = {
4020 VALID_CSPARSE, VALID_RSPARSE, VALID_TSPARSE, "" };
4021 int ivalid = R_check_class_etc(from, valid);
4022 if (ivalid < 0)
4023 ERROR_INVALID_CLASS(from, __func__);
4024
4025 return sparse_as_Tsparse(from, valid[ivalid]);
4026}
4027
4028/* as(<Matrix>, "vector") */
4029SEXP R_Matrix_as_vector(SEXP from)
4030{
4031 static const char *valid[] = { VALID_NONVIRTUAL_MATRIX, "" };
4032 int ivalid = R_check_class_etc(from, valid);
4033 if (ivalid < 0)
4034 ERROR_INVALID_CLASS(from, __func__);
4035 const char *cl = valid[ivalid + VALID_NONVIRTUAL_SHIFT(ivalid, 1)];
4036
4037 SEXP to = NULL;
4038 PROTECT_INDEX pid;
4039 PROTECT_WITH_INDEX(from, &pid);
4040
4041 switch (cl[2]) {
4042 case 'e':
4043 to = GET_SLOT(from, Matrix_xSym);
4044 if (cl[0] == 'n') {
4045 R_xlen_t len = XLENGTH(to);
4046 int *px = LOGICAL(to);
4047 while (len--)
4048 if (*(px++) == NA_LOGICAL) {
4049 PROTECT(to);
4050 to = duplicate(to);
4051 UNPROTECT(1);
4052 break;
4053 }
4054 }
4055 break;
4056 case 'y':
4057 case 'r':
4058 case 'p':
4059 REPROTECT(from = dense_as_general(from, cl, 1), pid);
4060 to = GET_SLOT(from, Matrix_xSym);
4061 break;
4062 case 'C':
4063 case 'R':
4064 case 'T':
4065 REPROTECT(from = sparse_as_dense(from, cl, 0), pid);
4066 REPROTECT(from = dense_as_general(from, cl, 0), pid);
4067 to = GET_SLOT(from, Matrix_xSym);
4068 break;
4069 case 'i':
4070 REPROTECT(from = diagonal_as_dense(from, cl, '.', 'g', 0, '\0'), pid);
4071 to = GET_SLOT(from, Matrix_xSym);
4072 break;
4073 case 'd':
4074 REPROTECT(from = index_as_dense(from, cl, 'n'), pid);
4075 to = GET_SLOT(from, Matrix_xSym);
4076 break;
4077 default:
4078 break;
4079 }
4080
4081 switch (cl[2]) {
4082 case 'e':
4083 case 'y':
4084 case 'r':
4085 case 'p':
4086 case 'i':
4087 if (cl[0] == 'n') {
4088 PROTECT(to);
4089 naToOne(to);
4090 UNPROTECT(1);
4091 }
4092 break;
4093 default:
4094 break;
4095 }
4096
4097 UNPROTECT(1); /* from */
4098 return to;
4099}
4100
4101/* as(<Matrix>, "matrix") */
4102SEXP R_Matrix_as_matrix(SEXP from)
4103{
4104 static const char *valid[] = { VALID_NONVIRTUAL_MATRIX, "" };
4105 int ivalid = R_check_class_etc(from, valid);
4106 if (ivalid < 0)
4107 ERROR_INVALID_CLASS(from, __func__);
4108 const char *cl = valid[ivalid + VALID_NONVIRTUAL_SHIFT(ivalid, 1)];
4109
4110 SEXP to = NULL;
4111 PROTECT_INDEX pid;
4112 PROTECT_WITH_INDEX(from, &pid);
4113
4114 switch (cl[2]) {
4115 case 'e':
4116 PROTECT(to = GET_SLOT(from, Matrix_xSym));
4117 to = duplicate(to);
4118 UNPROTECT(1);
4119 break;
4120 case 'y':
4121 case 'r':
4122 case 'p':
4123 REPROTECT(from = dense_as_general(from, cl, 1), pid);
4124 to = GET_SLOT(from, Matrix_xSym);
4125 break;
4126 case 'C':
4127 case 'R':
4128 case 'T':
4129 REPROTECT(from = sparse_as_dense(from, cl, 0), pid);
4130 REPROTECT(from = dense_as_general(from, cl, 0), pid);
4131 to = GET_SLOT(from, Matrix_xSym);
4132 break;
4133 case 'i':
4134 REPROTECT(from = diagonal_as_dense(from, cl, '.', 'g', 0, '\0'), pid);
4135 to = GET_SLOT(from, Matrix_xSym);
4136 break;
4137 case 'd':
4138 REPROTECT(from = index_as_dense(from, cl, 'n'), pid);
4139 to = GET_SLOT(from, Matrix_xSym);
4140 break;
4141 default:
4142 break;
4143 }
4144 PROTECT(to);
4145
4146 SEXP dim = PROTECT(GET_SLOT(from, Matrix_DimSym));
4147 setAttrib(to, R_DimSymbol, dim);
4148 UNPROTECT(1); /* dim */
4149
4150 SEXP dimnames = PROTECT(GET_SLOT(from, Matrix_DimNamesSym));
4151 if (!DimNames_is_trivial(dimnames))
4152 setAttrib(to, R_DimNamesSymbol, dimnames);
4153 UNPROTECT(1); /* dimnames */
4154
4155 switch (cl[2]) {
4156 case 'e':
4157 case 'y':
4158 case 'r':
4159 case 'p':
4160 case 'i':
4161 if (cl[0] == 'n')
4162 naToOne(to);
4163 break;
4164 default:
4165 break;
4166 }
4167
4168 UNPROTECT(2); /* to, from */
4169 return to;
4170}
4171
4172/* as(<Matrix>, "unpackedMatrix") */
4173SEXP R_Matrix_as_unpacked(SEXP from)
4174{
4175 static const char *valid[] = { VALID_NONVIRTUAL_MATRIX, "" };
4176 int ivalid = R_check_class_etc(from, valid);
4177 if (ivalid < 0)
4178 ERROR_INVALID_CLASS(from, __func__);
4179 const char *cl = valid[ivalid + VALID_NONVIRTUAL_SHIFT(ivalid, 1)];
4180
4181 switch (cl[2]) {
4182 case 'e':
4183 case 'y':
4184 case 'r':
4185 return from;
4186 case 'p':
4187 return dense_as_unpacked(from, valid[ivalid]);
4188 case 'C':
4189 case 'R':
4190 case 'T':
4191 return sparse_as_dense(from, cl, 0);
4192 case 'i':
4193 return diagonal_as_dense(from, cl, '.', 't', 0, 'U');
4194 case 'd':
4195 return index_as_dense(from, cl, 'n');
4196 default:
4197 return R_NilValue;
4198 }
4199}
4200
4201/* as(<Matrix>, "packedMatrix") */
4202SEXP R_Matrix_as_packed(SEXP from)
4203{
4204 static const char *valid[] = { VALID_NONVIRTUAL_MATRIX, "" };
4205 int ivalid = R_check_class_etc(from, valid);
4206 if (ivalid < 0)
4207 ERROR_INVALID_CLASS(from, __func__);
4208 const char *cl = valid[ivalid + VALID_NONVIRTUAL_SHIFT(ivalid, 1)];
4209
4210 if (cl[1] == 'g' || cl[2] == 'd')
4211 error(_("attempt to pack a %s"), "generalMatrix");
4212
4213 switch (cl[2]) {
4214 case 'y':
4215 case 'r':
4216 return dense_as_packed(from, valid[ivalid], '\0', '\0');
4217 case 'p':
4218 return from;
4219 case 'C':
4220 case 'R':
4221 case 'T':
4222 return sparse_as_dense(from, cl, 1);
4223 case 'i':
4224 return diagonal_as_dense(from, cl, '.', 't', 1, 'U');
4225 default:
4226 return R_NilValue;
4227 }
4228}
4229
4230/* as(<Matrix>, "CsparseMatrix") */
4231SEXP R_Matrix_as_Csparse(SEXP from)
4232{
4233 static const char *valid[] = { VALID_NONVIRTUAL_MATRIX, "" };
4234 int ivalid = R_check_class_etc(from, valid);
4235 if (ivalid < 0)
4236 ERROR_INVALID_CLASS(from, __func__);
4237 const char *cl = valid[ivalid + VALID_NONVIRTUAL_SHIFT(ivalid, 1)];
4238
4239 switch (cl[2]) {
4240 case 'e':
4241 case 'y':
4242 case 'r':
4243 case 'p':
4244 return dense_as_sparse(from, cl, 'C');
4245 case 'C':
4246 case 'R':
4247 case 'T':
4248 return sparse_as_Csparse(from, cl);
4249 case 'i':
4250 return diagonal_as_sparse(from, cl, '.', 't', 'C', 'U');
4251 case 'd':
4252 return index_as_sparse(from, cl, 'n', 'C');
4253 default:
4254 return R_NilValue;
4255 }
4256}
4257
4258/* as(<Matrix>, "RsparseMatrix") */
4259SEXP R_Matrix_as_Rsparse(SEXP from)
4260{
4261 static const char *valid[] = { VALID_NONVIRTUAL_MATRIX, "" };
4262 int ivalid = R_check_class_etc(from, valid);
4263 if (ivalid < 0)
4264 ERROR_INVALID_CLASS(from, __func__);
4265 const char *cl = valid[ivalid + VALID_NONVIRTUAL_SHIFT(ivalid, 1)];
4266
4267 switch (cl[2]) {
4268 case 'e':
4269 case 'y':
4270 case 'r':
4271 case 'p':
4272 return dense_as_sparse(from, cl, 'R');
4273 case 'C':
4274 case 'R':
4275 case 'T':
4276 return sparse_as_Rsparse(from, cl);
4277 case 'i':
4278 return diagonal_as_sparse(from, cl, '.', 't', 'R', 'U');
4279 case 'd':
4280 return index_as_sparse(from, cl, 'n', 'R');
4281 default:
4282 return R_NilValue;
4283 }
4284}
4285
4286/* as(<Matrix>, "TsparseMatrix") */
4287SEXP R_Matrix_as_Tsparse(SEXP from)
4288{
4289 static const char *valid[] = { VALID_NONVIRTUAL_MATRIX, "" };
4290 int ivalid = R_check_class_etc(from, valid);
4291 if (ivalid < 0)
4292 ERROR_INVALID_CLASS(from, __func__);
4293 const char *cl = valid[ivalid + VALID_NONVIRTUAL_SHIFT(ivalid, 1)];
4294
4295 switch (cl[2]) {
4296 case 'e':
4297 case 'y':
4298 case 'r':
4299 case 'p':
4300 return dense_as_sparse(from, cl, 'T');
4301 case 'C':
4302 case 'R':
4303 case 'T':
4304 return sparse_as_Tsparse(from, cl);
4305 case 'i':
4306 return diagonal_as_sparse(from, cl, '.', 't', 'T', 'U');
4307 case 'd':
4308 return index_as_sparse(from, cl, 'n', 'T');
4309 default:
4310 return R_NilValue;
4311 }
4312}
4313
4314/* as(<Matrix>, "[nldiz]Matrix") */
4315SEXP R_Matrix_as_kind(SEXP from, SEXP kind, SEXP sparse)
4316{
4317 static const char *valid[] = { VALID_NONVIRTUAL_MATRIX, "" };
4318 int ivalid = R_check_class_etc(from, valid);
4319 if (ivalid < 0)
4320 ERROR_INVALID_CLASS(from, __func__);
4321 const char *cl = valid[ivalid + VALID_NONVIRTUAL_SHIFT(ivalid, 1)];
4322
4323 char kind_;
4324 if (TYPEOF(kind) != STRSXP || LENGTH(kind) < 1 ||
4325 (kind = STRING_ELT(kind, 0)) == NA_STRING ||
4326 (kind_ = CHAR(kind)[0]) == '\0')
4327 error(_("invalid '%s' to '%s'"), "kind", __func__);
4328
4329 if (TYPEOF(sparse) != LGLSXP || LENGTH(sparse) < 1)
4330 error(_("'%s' must be %s or %s or %s"), "sparse", "TRUE", "FALSE", "NA");
4331 int sparse_ = LOGICAL(sparse)[0];
4332
4333 switch (cl[2]) {
4334 case 'e':
4335 case 'y':
4336 case 'r':
4337 case 'p':
4338 if (sparse_ == NA_LOGICAL || !sparse_)
4339 from = dense_as_kind(from, cl, kind_, 0);
4340 else {
4341 from = dense_as_sparse(from, cl, 'C');
4342 PROTECT(from);
4343 char cl_[] = "..CMatrix";
4344 cl_[0] = cl[0];
4345 cl_[1] = cl[1];
4346 from = sparse_as_kind(from, cl_, kind_);
4347 UNPROTECT(1);
4348 }
4349 return from;
4350 case 'C':
4351 case 'R':
4352 case 'T':
4353 from = sparse_as_kind(from, cl, kind_);
4354 if (sparse_ != NA_LOGICAL && !sparse_) {
4355 PROTECT(from);
4356 char cl_[] = "...Matrix";
4357 cl_[0] = (kind_ == '.') ? cl[0] : ((kind_ == ',') ? ((cl[0] == 'z') ? 'z' : 'd') : kind_);
4358 cl_[1] = cl[1];
4359 cl_[2] = cl[2];
4360 from = sparse_as_dense(from, cl_, 0);
4361 UNPROTECT(1);
4362 }
4363 return from;
4364 case 'i':
4365 if (sparse_ == NA_LOGICAL)
4366 from = diagonal_as_kind(from, cl, kind_);
4367 else if (sparse_)
4368 from = diagonal_as_sparse(from, cl, kind_, 't', 'C', 'U');
4369 else
4370 from = diagonal_as_dense(from, cl, kind_, 't', 0, 'U');
4371 return from;
4372 case 'd':
4373 if (sparse_ == NA_LOGICAL || sparse_)
4374 from = index_as_sparse(from, cl, kind_, '.');
4375 else
4376 from = index_as_dense(from, cl, kind_);
4377 return from;
4378 default:
4379 return R_NilValue;
4380 }
4381}
4382
4383/* as(as(<Matrix>, "[nlidz]Matrix"), "generalMatrix") */
4384SEXP R_Matrix_as_general(SEXP from, SEXP kind)
4385{
4386 static const char *valid[] = { VALID_NONVIRTUAL_MATRIX, "" };
4387 int ivalid = R_check_class_etc(from, valid);
4388 if (ivalid < 0)
4389 ERROR_INVALID_CLASS(from, __func__);
4390 const char *cl = valid[ivalid + VALID_NONVIRTUAL_SHIFT(ivalid, 1)];
4391
4392 char kind_;
4393 if (TYPEOF(kind) != STRSXP || LENGTH(kind) < 1 ||
4394 (kind = STRING_ELT(kind, 0)) == NA_STRING ||
4395 (kind_ = CHAR(kind)[0]) == '\0')
4396 error(_("invalid '%s' to '%s'"), "kind", __func__);
4397
4398 switch (cl[2]) {
4399 case 'e':
4400 case 'y':
4401 case 'r':
4402 case 'p':
4403 {
4404 char cl_[] = "...Matrix";
4405 cl_[0] = (kind_ == '.') ? cl[0] : ((kind_ == ',') ? ((cl[0] == 'z') ? 'z' : 'd') : kind_);
4406 cl_[1] = cl[1];
4407 cl_[2] = cl[2];
4408 from = dense_as_kind(from, cl, cl_[0], 1);
4409 PROTECT(from);
4410 from = dense_as_general(from, cl_, cl[0] == cl_[0]);
4411 UNPROTECT(1);
4412 return from;
4413 }
4414 case 'C':
4415 case 'R':
4416 case 'T':
4417 {
4418 char cl_[] = "...Matrix";
4419 cl_[0] = (kind_ == '.') ? cl[0] : ((kind_ == ',') ? ((cl[0] == 'z') ? 'z' : 'd') : kind_);
4420 cl_[1] = cl[1];
4421 cl_[2] = cl[2];
4422 from = sparse_as_kind(from, cl, cl_[0]);
4423 PROTECT(from);
4424 from = sparse_as_general(from, cl_);
4425 UNPROTECT(1);
4426 return from;
4427 }
4428 case 'i':
4429 return diagonal_as_sparse(from, cl, kind_, 'g', 'C', '\0');
4430 case 'd':
4431 /* indMatrix extends generalMatrix, but we typically do want this: */
4432 return index_as_sparse(from, cl, kind_, '.');
4433 default:
4434 return R_NilValue;
4435 }
4436}
VALID_RSPARSE
#define VALID_RSPARSE
Definition Mdefines.h:264
Matrix_int_fast64_t
long long Matrix_int_fast64_t
Definition Mdefines.h:27
VALID_TSPARSE
#define VALID_TSPARSE
Definition Mdefines.h:271
ERROR_INVALID_CLASS
#define ERROR_INVALID_CLASS(_X_, _FUNC_)
Definition Mdefines.h:208
FIRSTOF
#define FIRSTOF(x, y)
Definition Mdefines.h:99
VALID_DENSE
#define VALID_DENSE
Definition Mdefines.h:250
_
#define _(String)
Definition Mdefines.h:44
INCREMENT_REAL
#define INCREMENT_REAL(_X_, _Y_)
Definition Mdefines.h:169
typeToKind
char typeToKind(SEXPTYPE)
Definition objects.c:11
INCREMENT_COMPLEX
#define INCREMENT_COMPLEX(_X_, _Y_)
Definition Mdefines.h:173
INCREMENT_INTEGER
#define INCREMENT_INTEGER(_X_, _Y_)
Definition Mdefines.h:155
VALID_NONVIRTUAL_VECTOR
#define VALID_NONVIRTUAL_VECTOR
Definition Mdefines.h:240
ASSIGN_REAL
#define ASSIGN_REAL(_X_, _Y_)
Definition Mdefines.h:179
VALID_NONVIRTUAL_MATRIX
#define VALID_NONVIRTUAL_MATRIX
Definition Mdefines.h:220
PACKED_LENGTH
#define PACKED_LENGTH(m)
Definition Mdefines.h:198
kindToType
SEXPTYPE kindToType(char)
Definition objects.c:28
HIDE
#define HIDE(...)
Definition Mdefines.h:202
VALID_DIAGONAL
#define VALID_DIAGONAL
Definition Mdefines.h:278
ERROR_INVALID_TYPE
#define ERROR_INVALID_TYPE(_X_, _FUNC_)
Definition Mdefines.h:204
ISNZ_LOGICAL
#define ISNZ_LOGICAL(_X_)
Definition Mdefines.h:109
SHOW
#define SHOW(...)
Definition Mdefines.h:201
VALID_CSPARSE
#define VALID_CSPARSE
Definition Mdefines.h:257
Matrix_Free
#define Matrix_Free(_VAR_, _N_)
Definition Mdefines.h:77
INCREMENT_PATTERN
#define INCREMENT_PATTERN(_X_, _Y_)
Definition Mdefines.h:143
SET_SLOT
#define SET_SLOT(x, what, value)
Definition Mdefines.h:86
GET_SLOT
#define GET_SLOT(x, what)
Definition Mdefines.h:85
Matrix_Calloc
#define Matrix_Calloc(_VAR_, _N_, _CTYPE_)
Definition Mdefines.h:66
newObject
SEXP newObject(const char *)
Definition objects.c:4
kindToSize
size_t kindToSize(char)
Definition objects.c:46
SECONDOF
#define SECONDOF(x, y)
Definition Mdefines.h:100
INCREMENT_LOGICAL
#define INCREMENT_LOGICAL(_X_, _Y_)
Definition Mdefines.h:147
VALID_NONVIRTUAL_SHIFT
#define VALID_NONVIRTUAL_SHIFT(i, pToInd)
Definition Mdefines.h:247
Matrix_sdSym
SEXP Matrix_sdSym
Definition Msymbols.h:20
Matrix_permSym
SEXP Matrix_permSym
Definition Msymbols.h:18
Matrix_DimSym
SEXP Matrix_DimSym
Definition Msymbols.h:3
Matrix_factorsSym
SEXP Matrix_factorsSym
Definition Msymbols.h:12
Matrix_marginSym
SEXP Matrix_marginSym
Definition Msymbols.h:16
Matrix_xSym
SEXP Matrix_xSym
Definition Msymbols.h:22
Matrix_lengthSym
SEXP Matrix_lengthSym
Definition Msymbols.h:15
Matrix_iSym
SEXP Matrix_iSym
Definition Msymbols.h:13
Matrix_jSym
SEXP Matrix_jSym
Definition Msymbols.h:14
Matrix_DimNamesSym
SEXP Matrix_DimNamesSym
Definition Msymbols.h:2
Matrix_diagSym
SEXP Matrix_diagSym
Definition Msymbols.h:11
Matrix_uploSym
SEXP Matrix_uploSym
Definition Msymbols.h:21
Matrix_pSym
SEXP Matrix_pSym
Definition Msymbols.h:17
DimNames_is_trivial
int DimNames_is_trivial(SEXP dn)
Definition attrib.c:6
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
R_diagonal_as_dense
SEXP R_diagonal_as_dense(SEXP from, SEXP kind, SEXP shape, SEXP packed, SEXP uplo)
Definition coerce.c:893
R_diagonal_as_kind
SEXP R_diagonal_as_kind(SEXP from, SEXP kind)
Definition coerce.c:2696
dense_as_packed
SEXP dense_as_packed(SEXP from, const char *class, char ul, char di)
Definition coerce.c:3231
sparse_as_Csparse
SEXP sparse_as_Csparse(SEXP from, const char *class)
Definition coerce.c:3721
sparse_as_Rsparse
SEXP sparse_as_Rsparse(SEXP from, const char *class)
Definition coerce.c:3820
R_sparse_as_general
SEXP R_sparse_as_general(SEXP from)
Definition coerce.c:3108
R_Matrix_as_packed
SEXP R_Matrix_as_packed(SEXP from)
Definition coerce.c:4202
diagonal_as_sparse
SEXP diagonal_as_sparse(SEXP from, const char *class, char kind, char shape, char repr, char ul)
Definition coerce.c:2066
VAS_SUBCASES
#define VAS_SUBCASES(...)
R_vector_as_dense
SEXP R_vector_as_dense(SEXP from, SEXP zzz, SEXP uplo, SEXP diag, SEXP m, SEXP n, SEXP byrow, SEXP dimnames)
Definition coerce.c:164
SAG_CASES
#define SAG_CASES
taggr
void taggr(SEXP i0, SEXP j0, SEXP x0, SEXP *i1, SEXP *j1, SEXP *x1, int m, int n)
Definition coerce.c:3585
SAD_SUBCASES
#define SAD_SUBCASES(_CTYPE_, _PTR_, _MASK_, _REPLACE_, _INCREMENT_, _ONE_)
R_sparse_as_kind
SEXP R_sparse_as_kind(SEXP from, SEXP kind)
Definition coerce.c:2617
R_Matrix_as_matrix
SEXP R_Matrix_as_matrix(SEXP from)
Definition coerce.c:4102
R_dense_as_general
SEXP R_dense_as_general(SEXP from)
Definition coerce.c:2821
sparse_as_dense
SEXP sparse_as_dense(SEXP from, const char *class, int packed)
Definition coerce.c:501
IAS_SUBCASES
#define IAS_SUBCASES(_CTYPE_, _PTR_, _ONE_)
R_sparse_as_dense
SEXP R_sparse_as_dense(SEXP from, SEXP packed)
Definition coerce.c:782
dense_as_sparse
SEXP dense_as_sparse(SEXP from, const char *class, char repr)
Definition coerce.c:1647
DAS_SUBCASES
#define DAS_SUBCASES(_CTYPE_, _PTR_, _MASK_, _ISNZ_)
sparse_as_kind
SEXP sparse_as_kind(SEXP from, const char *class, char kind)
Definition coerce.c:2499
dense_as_kind
SEXP dense_as_kind(SEXP from, const char *class, char kind, int new)
Definition coerce.c:2409
R_dense_as_packed
SEXP R_dense_as_packed(SEXP from, SEXP uplo, SEXP diag)
Definition coerce.c:3330
diagonal_as_kind
SEXP diagonal_as_kind(SEXP from, const char *class, char kind)
Definition coerce.c:2634
VAD_SUBCASES
#define VAD_SUBCASES(_PREFIX_, _CTYPE_, _PTR_, _NA_)
index_as_dense
SEXP index_as_dense(SEXP from, const char *class, char kind)
Definition coerce.c:931
R_dense_as_sparse
SEXP R_dense_as_sparse(SEXP from, SEXP repr)
Definition coerce.c:2050
DAS_CASES
#define DAS_CASES(_MASK_)
R_Matrix_as_unpacked
SEXP R_Matrix_as_unpacked(SEXP from)
Definition coerce.c:4173
SAD_CASES
#define SAD_CASES
DAG_SUBCASES
#define DAG_SUBCASES(_PREFIX_, _CTYPE_, _PTR_)
R_matrix_as_sparse
SEXP R_matrix_as_sparse(SEXP from, SEXP zzz, SEXP uplo, SEXP diag, SEXP trans)
Definition coerce.c:1602
R_Matrix_as_Rsparse
SEXP R_Matrix_as_Rsparse(SEXP from)
Definition coerce.c:4259
R_vector_as_sparse
SEXP R_vector_as_sparse(SEXP from, SEXP zzz, SEXP uplo, SEXP diag, SEXP m, SEXP n, SEXP byrow, SEXP dimnames)
Definition coerce.c:1453
TRANS_LOOP
#define TRANS_LOOP(_CTYPE_, _PTR_, _MASK_)
sparse_as_Tsparse
SEXP sparse_as_Tsparse(SEXP from, const char *class)
Definition coerce.c:3919
R_sparse_as_Rsparse
SEXP R_sparse_as_Rsparse(SEXP from)
Definition coerce.c:3908
PACK
#define PACK(_PREFIX_, _CTYPE_, _PTR_)
vector_as_sparse
SEXP vector_as_sparse(SEXP from, const char *zzz, char ul, char di, int m, int n, int byrow, SEXP dimnames)
Definition coerce.c:1020
matrix_as_dense
SEXP matrix_as_dense(SEXP from, const char *zzz, char ul, char di, int trans, int new)
Definition coerce.c:297
SAG_SUBCASES
#define SAG_SUBCASES(_CTYPE_, _PTR_, _MASK_, _ONE_, _ASSIGN_JJ_, _ASSIGN_JI_)
trans
void trans(SEXP p0, SEXP i0, SEXP x0, SEXP p1, SEXP i1, SEXP x1, int m, int n)
Definition coerce.c:3356
R_sparse_as_Tsparse
SEXP R_sparse_as_Tsparse(SEXP from)
Definition coerce.c:4017
R_Matrix_as_general
SEXP R_Matrix_as_general(SEXP from, SEXP kind)
Definition coerce.c:4384
R_dense_as_kind
SEXP R_dense_as_kind(SEXP from, SEXP kind)
Definition coerce.c:2483
R_Matrix_as_vector
SEXP R_Matrix_as_vector(SEXP from)
Definition coerce.c:4029
R_index_as_kind
SEXP R_index_as_kind(SEXP from, SEXP kind)
Definition coerce.c:2718
R_dense_as_unpacked
SEXP R_dense_as_unpacked(SEXP from)
Definition coerce.c:3219
UNPACK
#define UNPACK(_PREFIX_, _CTYPE_, _PTR_)
IAD_SUBCASES
#define IAD_SUBCASES(_CTYPE_, _PTR_, _ONE_)
R_Matrix_as_Csparse
SEXP R_Matrix_as_Csparse(SEXP from)
Definition coerce.c:4231
tsort
void tsort(SEXP i0, SEXP j0, SEXP x0, SEXP *p1, SEXP *i1, SEXP *x1, int m, int n)
Definition coerce.c:3415
R_sparse_as_Csparse
SEXP R_sparse_as_Csparse(SEXP from)
Definition coerce.c:3809
dense_as_unpacked
SEXP dense_as_unpacked(SEXP from, const char *class)
Definition coerce.c:3119
matrix_as_sparse
SEXP matrix_as_sparse(SEXP from, const char *zzz, char ul, char di, int trans)
Definition coerce.c:1580
vector_as_dense
SEXP vector_as_dense(SEXP from, const char *zzz, char ul, char di, int m, int n, int byrow, SEXP dimnames)
Definition coerce.c:6
dense_as_general
SEXP dense_as_general(SEXP from, const char *class, int new)
Definition coerce.c:2734
DAD_SUBCASES
#define DAD_SUBCASES(_PREFIX_, _CTYPE_, _PTR_)
R_index_as_dense
SEXP R_index_as_dense(SEXP from, SEXP kind)
Definition coerce.c:1004
TAGGR_LOOP
#define TAGGR_LOOP(_CTYPE_, _PTR_, _MASK_, _INCREMENT_)
sparse_as_general
SEXP sparse_as_general(SEXP from, const char *class)
Definition coerce.c:2831
R_diagonal_as_sparse
SEXP R_diagonal_as_sparse(SEXP from, SEXP kind, SEXP shape, SEXP repr, SEXP uplo)
Definition coerce.c:2244
R_Matrix_as_Tsparse
SEXP R_Matrix_as_Tsparse(SEXP from)
Definition coerce.c:4287
index_as_kind
SEXP index_as_kind(SEXP from, const char *class, char kind)
Definition coerce.c:2712
TSORT_LOOP
#define TSORT_LOOP(_CTYPE_, _PTR_, _MASK_, _INCREMENT_)
diagonal_as_dense
SEXP diagonal_as_dense(SEXP from, const char *class, char kind, char shape, int packed, char ul)
Definition coerce.c:798
R_index_as_sparse
SEXP R_index_as_sparse(SEXP from, SEXP kind, SEXP repr)
Definition coerce.c:2386
R_matrix_as_dense
SEXP R_matrix_as_dense(SEXP from, SEXP zzz, SEXP uplo, SEXP diag, SEXP trans)
Definition coerce.c:454
DAS_LOOP
#define DAS_LOOP(_CTYPE_, _PTR_, _MASK_, _ISNZ_, _ONE_)
index_as_sparse
SEXP index_as_sparse(SEXP from, const char *class, char kind, char repr)
Definition coerce.c:2281
R_Matrix_as_kind
SEXP R_Matrix_as_kind(SEXP from, SEXP kind, SEXP sparse)
Definition coerce.c:4315
Matrix_zna
Rcomplex Matrix_zna
Definition init.c:26
Matrix_zone
Rcomplex Matrix_zone
Definition init.c:26
Tsparse_aggregate
SEXP Tsparse_aggregate(SEXP from)
Definition sparse.c:3672
naToOne
void naToOne(SEXP x)
Definition utils.c:185
Matrix_memset
void * Matrix_memset(void *dest, int ch, R_xlen_t length, size_t size)
Definition utils.c:8
Matrix_memcpy
void * Matrix_memcpy(void *dest, const void *src, R_xlen_t length, size_t size)
Definition utils.c:70
Generated by doxygen 1.9.8