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