Matrix r5059
Loading...
Searching...
No Matches
matmult.c
Go to the documentation of this file.
1/* C implementation of methods for %&%, %*%, crossprod, tcrossprod */
2
3#include "Lapack-etc.h"
4#include "cholmod-etc.h"
5#include "Mdefines.h"
6#include "M5.h"
7#include "idz.h"
8#include "coerce.h"
9
10SEXP dense_transpose(SEXP, const char *, char);
11SEXP sparse_transpose(SEXP, const char *, char, int);
12SEXP sparse_dropzero(SEXP, const char *, double);
13SEXP sparse_diag_U2N(SEXP, const char *);
14
15static const char *valid_matmult[] = {
16VALID_DIAGONAL, VALID_SPARSE, VALID_DENSE, "" };
17/* ptr: 0, 5, 56, 87 */
18
19static
20void matmultDim(SEXP x, SEXP y, char *xtrans, char *ytrans, char *ztrans,
21 int *m, int *n, int *v)
22{
23 int xt = *xtrans == 'C' || *xtrans == 'T'; if (!xt) *xtrans = 'N';
24 int yt = *ytrans == 'C' || *ytrans == 'T'; if (!yt) *ytrans = 'N';
25 int zt = *ztrans == 'C' || *ztrans == 'T'; if (!zt) *ztrans = 'N';
26 if (y == R_NilValue) {
27 if (xt == yt)
28 Rf_error(_("should never happen ..."));
29 SEXP
30 xdim = (TYPEOF(x) == OBJSXP)
31 ? GET_SLOT(x, Matrix_DimSym) : Rf_getAttrib(x, R_DimSymbol);
32 if (TYPEOF(xdim) == INTSXP && LENGTH(xdim) == 2) {
33 *v = 0;
34 *m = *n = INTEGER(xdim)[(xt) ? 1 : 0];
35 } else if (XLENGTH(x) <= INT_MAX) {
36 *v = 1;
37 *m = *n = (xt) ? 1 : LENGTH(x);
38 } else
39 Rf_error(_("dimensions cannot exceed %s"), "2^31-1");
40 } else {
41 /* MJ: So that I don't lose my mind ... : */
42 if (zt) {
43 SWAP(x, y, SEXP, );
44 SWAP(xt, yt, int, !);
45 }
46 SEXP
47 xdim = (TYPEOF(x) == OBJSXP)
48 ? GET_SLOT(x, Matrix_DimSym) : Rf_getAttrib(x, R_DimSymbol),
49 ydim = (TYPEOF(y) == OBJSXP)
50 ? GET_SLOT(y, Matrix_DimSym) : Rf_getAttrib(y, R_DimSymbol);
51 int xm, xn, ym, yn, x2, y2;
52 xm = xn = ym = yn = -1;
53 x2 = TYPEOF(xdim) == INTSXP && LENGTH(xdim) == 2;
54 y2 = TYPEOF(ydim) == INTSXP && LENGTH(ydim) == 2;
55 if (x2) {
56 int *pxdim = INTEGER(xdim);
57 xm = pxdim[0];
58 xn = pxdim[1];
59 } else if (XLENGTH(x) > INT_MAX)
60 Rf_error(_("dimensions cannot exceed %s"), "2^31-1");
61 if (y2) {
62 int *pydim = INTEGER(ydim);
63 ym = pydim[0];
64 yn = pydim[1];
65 } else if (XLENGTH(y) > INT_MAX)
66 Rf_error(_("dimensions cannot exceed %s"), "2^31-1");
67 /* MJ: R's do_matprod behaves quite asymmetrically ... what a pain */
68 if (x2 && y2)
69 *v = 0;
70 else if (y2) {
71 *v = (zt) ? 2 : 1;
72 int k = (yt) ? yn : ym, xl = LENGTH(x);
73 if (k == xl || (k == 1 && !(xt))) {
74 xm = (int) xl;
75 xn = 1;
76 xt = (k == xl) ? 1 : 0;
77 }
78 } else if (x2) {
79 *v = (zt) ? 1 : 2;
80 int k = (xt) ? xm : xn, yl = LENGTH(y);
81 if (yt) {
82 if (xm == 1 || xn == 1) {
83 ym = (int) yl;
84 yn = 1;
85 yt = (((xt) ? xn : xm) == 1) ? 0 : 1;
86 }
87 } else {
88 if (k == yl || k == 1) {
89 ym = (int) yl;
90 yn = 1;
91 yt = (k == yl) ? 0 : 1;
92 }
93 }
94 } else {
95 *v = 3;
96 int xl = LENGTH(x), yl = LENGTH(y);
97 if (xt) {
98 xm = xl;
99 xn = 1;
100 ym = yl;
101 yn = 1;
102 yt = xl == 1;
103 } else if (yt) {
104 xm = xl;
105 xn = 1;
106 ym = yl;
107 yn = 1;
108 /* xt = 0; */
109 } else {
110 xm = 1;
111 xn = xl;
112 ym = (xl == 1) ? 1 : yl;
113 yn = (xl == 1) ? yl : 1;
114 }
115 }
116 if (((xt) ? xm : xn) != ((yt) ? yn : ym))
117 Rf_error(_("non-conformable arguments"));
118 *m = (xt) ? xn : xm;
119 *n = (yt) ? ym : yn;
120 if (zt) {
121 SWAP(*m, *n, int, );
122 SWAP(xt, yt, int, !);
123 }
124 }
125 if (*v % 2) *xtrans = (xt) ? 'T' : 'N';
126 if (*v > 1) *ytrans = (yt) ? 'T' : 'N';
127 return;
128}
129
130static
131void matmultDN(SEXP dest, SEXP asrc, int ai, SEXP bsrc, int bi) {
132 SEXP s;
133 if ((s = VECTOR_ELT(asrc, ai)) != R_NilValue)
134 SET_VECTOR_ELT(dest, 0, s);
135 if ((s = VECTOR_ELT(bsrc, bi)) != R_NilValue)
136 SET_VECTOR_ELT(dest, 1, s);
137 PROTECT(asrc = Rf_getAttrib(asrc, R_NamesSymbol));
138 PROTECT(bsrc = Rf_getAttrib(bsrc, R_NamesSymbol));
139 if (asrc != R_NilValue || bsrc != R_NilValue) {
140 SEXP destnms = PROTECT(Rf_allocVector(STRSXP, 2));
141 if (asrc != R_NilValue && CHAR(s = STRING_ELT(asrc, ai))[0] != '\0')
142 SET_STRING_ELT(destnms, 0, s);
143 if (bsrc != R_NilValue && CHAR(s = STRING_ELT(bsrc, bi))[0] != '\0')
144 SET_STRING_ELT(destnms, 1, s);
145 Rf_setAttrib(dest, R_NamesSymbol, destnms);
146 UNPROTECT(1);
147 }
148 UNPROTECT(2);
149 return;
150}
151
152#define CONJ2(_X_, _M_, _N_) \
153do { \
154 size_t m = (size_t) _M_, n = (size_t) _N_, xlen = m * n; \
155 Rcomplex *x = (Rcomplex *) _X_; \
156 Rcomplex *y = (Rcomplex *) R_alloc(xlen, sizeof(Rcomplex)); \
157 zvconj(x, y, xlen); \
158 _X_ = y; \
159} while (0)
160
161#define CONJ1(_X_, _N_) \
162do { \
163 size_t n = (size_t) _N_, xlen = PACKED_LENGTH(n); \
164 Rcomplex *x = (Rcomplex *) _X_; \
165 Rcomplex *y = (Rcomplex *) R_alloc(xlen, sizeof(Rcomplex)); \
166 zvconj(x, y, xlen); \
167 _X_ = y; \
168} while (0)
169
170/* op(<,ge>) * op(<,ge>) */
171static
172SEXP geMatrix_matmult(SEXP a, SEXP b, char atrans, char btrans)
173{
174 int *padim = DIM(a), am = padim[0], an = padim[1],
175 rm = (atrans != 'N') ? an : am, rk = (atrans != 'N') ? am : an;
176
177 if (b == R_NilValue) {
178
179 if ((int_fast64_t) rm * rm > R_XLEN_T_MAX)
180 Rf_error(_("attempt to allocate vector of length exceeding %s"),
181 "R_XLEN_T_MAX");
182
183 SEXP ax = PROTECT(GET_SLOT(a, Matrix_xSym));
184
185 char rct = (TYPEOF(ax) != CPLXSXP || ((atrans != 'N') ? atrans : btrans) == 'C') ? 'C' : 'T';
186
187 char rcl[] = "...Matrix";
188 rcl[0] = (TYPEOF(ax) == CPLXSXP) ? 'z' : 'd';
189 rcl[1] = (rct == 'C') ? 'p' : 's';
190 rcl[2] = (rct == 'C') ? 'o' : 'y';
191
192 SEXP r = PROTECT(newObject(rcl));
193
194 SET_DIM(r, rm, rm);
195
196 SEXP adimnames = PROTECT(DIMNAMES(a, 0)),
197 rdimnames = PROTECT(DIMNAMES(r, 0));
198 symDN(rdimnames, adimnames, (atrans != 'N') ? 1 : 0);
199 UNPROTECT(2); /* rdimnames, adimnames */
200
201 if (rct != 'C')
202 SET_TRANS(r);
203
204 if (rm > 0) {
205 SEXP rx = PROTECT(Rf_allocVector(TYPEOF(ax), (R_xlen_t) rm * rm));
206 if (TYPEOF(ax) == CPLXSXP) {
207 Rcomplex *prx = COMPLEX(rx);
208 memset(prx, 0, sizeof(Rcomplex) * (size_t) rm * (size_t) rm);
209 if (rk > 0) {
210 Rcomplex *pax = COMPLEX(ax),
211 zero = Matrix_zzero, one = Matrix_zunit;
212 if (rct == 'C')
213 F77_CALL(zherk)("U", &atrans, &rm, &rk,
214 & one.r, pax, &am,
215 &zero.r, prx, &rm FCONE FCONE);
216 else
217 F77_CALL(zsyrk)("U", &atrans, &rm, &rk,
218 & one , pax, &am,
219 &zero , prx, &rm FCONE FCONE);
220 }
221 } else {
222 double *prx = REAL(rx);
223 memset(prx, 0, sizeof(double) * (size_t) rm * (size_t) rm);
224 if (rk > 0) {
225 double *pax = REAL(ax),
226 zero = 0.0, one = 1.0;
227 F77_CALL(dsyrk)("U", &atrans, &rm, &rk,
228 & one , pax, &am,
229 &zero , prx, &rm FCONE FCONE);
230 }
231 }
232 SET_SLOT(r, Matrix_xSym, rx);
233 UNPROTECT(1); /* rx */
234 }
235
236 UNPROTECT(2); /* r, ax */
237 return r;
238
239 } else {
240
241 int *pbdim = DIM(b), bm = pbdim[0], bn = pbdim[1],
242 rn = (btrans != 'N') ? bm : bn;
243
244 if (rk != ((btrans != 'N') ? bn : bm))
245 Rf_error(_("non-conformable arguments"));
246 if ((int_fast64_t) rm * rn > R_XLEN_T_MAX)
247 Rf_error(_("attempt to allocate vector of length exceeding %s"),
248 "R_XLEN_T_MAX");
249
250 SEXP ax = PROTECT(GET_SLOT(a, Matrix_xSym));
251
252 char rcl[] = ".geMatrix";
253 rcl[0] = (TYPEOF(ax) == CPLXSXP) ? 'z' : 'd';
254 SEXP r = PROTECT(newObject(rcl));
255
256 SET_DIM(r, rm, rn);
257
258 SEXP adimnames = PROTECT(DIMNAMES(a, 0)),
259 bdimnames = PROTECT(DIMNAMES(b, 0)),
260 rdimnames = PROTECT(DIMNAMES(r, 0));
261 matmultDN(rdimnames,
262 adimnames, (atrans != 'N') ? 1 : 0,
263 bdimnames, (btrans != 'N') ? 0 : 1);
264 UNPROTECT(3); /* rdimnames, bdimnames, adimnames */
265
266 if (rm > 0 && rn > 0) {
267 SEXP rx = PROTECT(Rf_allocVector(TYPEOF(ax), (R_xlen_t) rm * rn));
268 if (TYPEOF(ax) == CPLXSXP) {
269 Rcomplex *prx = COMPLEX(rx);
270 if (rk == 0)
271 memset(prx, 0, sizeof(Rcomplex) * (size_t) rm * (size_t) rn);
272 else {
273 SEXP bx = PROTECT(GET_SLOT(b, Matrix_xSym));
274 Rcomplex *pax = COMPLEX(ax), *pbx = COMPLEX(bx),
275 zero = Matrix_zzero, one = Matrix_zunit;
276 F77_CALL(zgemm)(&atrans, &btrans, &rm, &rn, &rk,
277 & one, pax, &am, pbx, &bm,
278 &zero, prx, &rm FCONE FCONE);
279 UNPROTECT(1); /* bx */
280 }
281 } else {
282 double *prx = REAL(rx);
283 if (rk == 0)
284 memset(prx, 0, sizeof(double) * (size_t) rm * (size_t) rn);
285 else {
286 SEXP bx = PROTECT(GET_SLOT(b, Matrix_xSym));
287 double *pax = REAL(ax), *pbx = REAL(bx),
288 zero = 0.0, one = 1.0;
289 F77_CALL(dgemm)(&atrans, &btrans, &rm, &rn, &rk,
290 & one, pax, &am, pbx, &bm,
291 &zero, prx, &rm FCONE FCONE);
292 UNPROTECT(1); /* bx */
293 }
294 }
295 SET_SLOT(r, Matrix_xSym, rx);
296 UNPROTECT(1); /* rx */
297 }
298
299 UNPROTECT(2); /* r, ax */
300 return r;
301
302 }
303}
304
305/* op(<,sy>) * op(<,ge>) or op(<,ge>) * op(<,sy>) */
306static
307SEXP syMatrix_matmult(SEXP a, SEXP b, char atrans, char btrans, char aside)
308{
309 int *padim = DIM(a),
310 rk = padim[0];
311
312 int *pbdim = DIM(b), bm = pbdim[0], bn = pbdim[1],
313 rm = (btrans != 'N') ? bn : bm, rn = (btrans != 'N') ? bm : bn;
314
315 if (rk != (((aside == 'L') == (btrans != 'N')) ? bn : bm))
316 Rf_error(_("non-conformable arguments"));
317 if ((int_fast64_t) rm * rn > R_XLEN_T_MAX)
318 Rf_error(_("attempt to allocate vector of length exceeding %s"),
319 "R_XLEN_T_MAX");
320
321 SEXP ax = PROTECT(GET_SLOT(a, Matrix_xSym));
322
323 char rcl[] = ".geMatrix";
324 rcl[0] = (TYPEOF(ax) == CPLXSXP) ? 'z' : 'd';
325 SEXP r = PROTECT(newObject(rcl));
326
327 SET_DIM(r, rm, rn);
328
329 SEXP adimnames = PROTECT(DIMNAMES(a, -1)),
330 bdimnames = PROTECT(DIMNAMES(b, 0)),
331 rdimnames = PROTECT(DIMNAMES(r, 0));
332 if (aside == 'L')
333 matmultDN(rdimnames, adimnames, 0, bdimnames, btrans == 'N');
334 else
335 matmultDN(rdimnames, bdimnames, btrans != 'N', adimnames, 1);
336 UNPROTECT(3); /* rdimnames, bdimnames, adimnames */
337
338 /* use *mm */
339 /* R := A B */
340 /* R := A' B */
341 /* R := B A */
342 /* R := B A' */
343 /* use *mv and access B by row */
344 /* R := A B. */
345 /* R := A' B. */
346 /* use *mv and access R by row */
347 /* R := B. A = (A. B). */
348 /* R := B. A' = (A'. B). */
349
350 if (rm > 0 && rn > 0) {
351 SEXP bx = PROTECT(GET_SLOT(b, Matrix_xSym)),
352 rx = PROTECT(Rf_allocVector(TYPEOF(ax), (R_xlen_t) rm * rn));
353 char aul = UPLO(a);
354 int i,
355 d = (aside == 'L') ? rn : rm,
356 binc = (aside == 'L') ? bm : 1,
357 bincp = (aside == 'L') ? 1 : bm,
358 rinc = (aside == 'L') ? 1 : rm,
359 rincp = (aside == 'L') ? rm : 1;
360
361 if (TYPEOF(ax) == CPLXSXP) {
362 Rcomplex *pax = COMPLEX(ax), *pbx = COMPLEX(bx), *prx = COMPLEX(rx),
363 zero = Matrix_zzero, one = Matrix_zunit;
364 char act = TRANS(a);
365 if (btrans == 'N') {
366 if (atrans != 'N' && atrans != act)
367 CONJ2(pax, rk, rk);
368 if (act == 'C')
369 F77_CALL(zhemm)(&aside, &aul, &rm, &rn,
370 & one, pax, &rk, pbx, &bm,
371 &zero, prx, &rm FCONE FCONE);
372 else
373 F77_CALL(zsymm)(&aside, &aul, &rm, &rn,
374 & one, pax, &rk, pbx, &bm,
375 &zero, prx, &rm FCONE FCONE);
376 } else {
377 if (aside == 'L') {
378 if (atrans != 'N' && atrans != act)
379 CONJ2(pax, rk, rk);
380 if (btrans == 'C')
381 CONJ2(pbx, bm, bn);
382 } else {
383 if (((atrans != 'N') ? atrans : act) != btrans)
384 CONJ2(pax, rk, rk);
385 }
386 for (i = 0; i < d; ++i) {
387 if (act == 'C')
388 F77_CALL(zhemv)( &aul, &rk,
389 & one, pax, &rk, pbx, &binc,
390 &zero, prx, &rinc FCONE);
391 else
392 F77_CALL(zsymv)( &aul, &rk,
393 & one, pax, &rk, pbx, &binc,
394 &zero, prx, &rinc FCONE);
395 pbx += bincp;
396 prx += rincp;
397 }
398 if (aside != 'L' && btrans == 'C')
399 zvconj(pax, NULL, (size_t) rk * (size_t) rk);
400 }
401 } else {
402 double *pax = REAL(ax), *pbx = REAL(bx), *prx = REAL(rx),
403 zero = 0.0, one = 1.0;
404 if (btrans == 'N') {
405 F77_CALL(dsymm)(&aside, &aul, &rm, &rn,
406 & one, pax, &rk, pbx, &bm,
407 &zero, prx, &rm FCONE FCONE);
408 } else {
409 for (i = 0; i < d; ++i) {
410 F77_CALL(dsymv)( &aul, &rk,
411 &one, pax, &rk, pbx, &binc,
412 &zero, prx, &rinc FCONE);
413 pbx += bincp;
414 prx += rincp;
415 }
416 }
417 }
418 SET_SLOT(r, Matrix_xSym, rx);
419 UNPROTECT(2); /* rx, bx */
420 }
421
422 UNPROTECT(2); /* r, ax */
423 return r;
424}
425
426/* op(<,sp>) * op(<,ge>) or op(<,ge>) * op(<,sp>) */
427static
428SEXP spMatrix_matmult(SEXP a, SEXP b, char atrans, char btrans, char aside)
429{
430 int *padim = DIM(a),
431 rk = padim[0];
432
433 int *pbdim = DIM(b), bm = pbdim[0], bn = pbdim[1],
434 rm = (btrans != 'N') ? bn : bm, rn = (btrans != 'N') ? bm : bn;
435
436 if (rk != (((aside == 'L') == (btrans != 'N')) ? bn : bm))
437 Rf_error(_("non-conformable arguments"));
438 if ((int_fast64_t) rm * rn > R_XLEN_T_MAX)
439 Rf_error(_("attempt to allocate vector of length exceeding %s"),
440 "R_XLEN_T_MAX");
441
442 SEXP ax = PROTECT(GET_SLOT(a, Matrix_xSym));
443
444 char rcl[] = ".geMatrix";
445 rcl[0] = (TYPEOF(ax) == CPLXSXP) ? 'z' : 'd';
446 SEXP r = PROTECT(newObject(rcl));
447
448 SET_DIM(r, rm, rn);
449
450 SEXP adimnames = PROTECT(DIMNAMES(a, -1)),
451 bdimnames = PROTECT(DIMNAMES(b, 0)),
452 rdimnames = PROTECT(DIMNAMES(r, 0));
453 if (aside == 'L')
454 matmultDN(rdimnames, adimnames, 0, bdimnames, btrans == 'N');
455 else
456 matmultDN(rdimnames, bdimnames, btrans != 'N', adimnames, 1);
457 UNPROTECT(3); /* rdimnames, bdimnames, adimnames */
458
459 /* use *mv */
460 /* R := A B */
461 /* R := A' B */
462 /* use *mv and access B by row */
463 /* R := A B. */
464 /* R := A' B. */
465 /* use *mv and access R by row */
466 /* R := B. A = (A. B ). */
467 /* R := B. A' = (A'. B ). */
468 /* use *mv and access B and R by row */
469 /* R := B A = (A* B*)* */
470 /* R := B A' = (A B')' */
471
472 if (rm > 0 && rn > 0) {
473 SEXP bx = PROTECT(GET_SLOT(b, Matrix_xSym)),
474 rx = PROTECT(Rf_allocVector(REALSXP, (R_xlen_t) rm * rn));
475 char aul = UPLO(a);
476 int i,
477 d = ( aside == 'L' ) ? rn : rm,
478 binc = ((aside == 'L') == (btrans != 'N')) ? bm : 1,
479 bincp = ((aside == 'L') == (btrans != 'N')) ? 1 : bm,
480 rinc = ( aside == 'L' ) ? 1 : rm,
481 rincp = ( aside == 'L' ) ? rm : 1;
482
483 if (TYPEOF(ax) == CPLXSXP) {
484 Rcomplex *pax = COMPLEX(ax), *pbx = COMPLEX(bx), *prx = COMPLEX(rx),
485 zero = Matrix_zzero, one = Matrix_zunit;
486 char act = TRANS(a);
487 if (aside == 'L') {
488 if (atrans != 'N' && atrans != act)
489 CONJ1(pax, rk);
490 if (btrans == 'C')
491 CONJ2(pbx, bm, bn);
492 } else {
493 if (btrans != 'N' && ((atrans != 'N') ? atrans : act) != btrans)
494 CONJ1(pax, rk);
495 if (btrans == 'N' && ((atrans != 'N') ? atrans : act) == 'C')
496 CONJ2(pbx, bm, bn);
497 }
498 for (i = 0; i < d; ++i) {
499 if (act == 'C')
500 F77_CALL(zhpmv)(&aul, &rk,
501 & one, pax, pbx, &binc,
502 &zero, prx, &rinc FCONE);
503 else
504 F77_CALL(zspmv)(&aul, &rk,
505 & one, pax, pbx, &binc,
506 &zero, prx, &rinc FCONE);
507 pbx += bincp;
508 prx += rincp;
509 }
510 if (aside != 'L' &&
511 ((btrans != 'N') ? btrans : ((atrans != 'N') ? atrans : act)) == 'C')
512 zvconj(prx, NULL, (size_t) rm * (size_t) rn);
513 } else {
514 double *pax = REAL(ax), *pbx = REAL(bx), *prx = REAL(rx),
515 zero = 0.0, one = 1.0;
516 for (i = 0; i < d; ++i) {
517 F77_CALL(dspmv)(&aul, &rk,
518 & one, pax, pbx, &binc,
519 &zero, prx, &rinc FCONE);
520 pbx += bincp;
521 prx += rincp;
522 }
523 }
524 SET_SLOT(r, Matrix_xSym, rx);
525 UNPROTECT(2); /* rx, bx */
526 }
527
528 UNPROTECT(2); /* r, ax */
529 return r;
530}
531
532/* op(<,tr>) * op(<,ge>) or op(<,ge>) * op(<,tr>) */
533static
534SEXP trMatrix_matmult(SEXP a, SEXP b, char atrans, char btrans, char aside,
535 int triangular)
536{
537 int *padim = DIM(a),
538 rk = padim[0];
539
540 int *pbdim = DIM(b), bm = pbdim[0], bn = pbdim[1],
541 rm = (btrans != 'N') ? bn : bm, rn = (btrans != 'N') ? bm : bn;
542
543 if (rk != (((aside == 'L') == (btrans != 'N')) ? bn : bm))
544 Rf_error(_("non-conformable arguments"));
545 if ((int_fast64_t) rm * rn > R_XLEN_T_MAX)
546 Rf_error(_("attempt to allocate vector of length exceeding %s"),
547 "R_XLEN_T_MAX");
548
549 SEXP ax = PROTECT(GET_SLOT(a, Matrix_xSym));
550
551 char rcl[] = "...Matrix";
552 rcl[0] = (TYPEOF(ax) == CPLXSXP) ? 'z' : 'd';
553 rcl[1] = (triangular != 0) ? 't' : 'g';
554 rcl[2] = (triangular != 0) ? 'r' : 'e';
555 SEXP r = PROTECT(newObject(rcl));
556
557 SET_DIM(r, rm, rn);
558
559 SEXP adimnames = PROTECT(DIMNAMES(a, 0)),
560 bdimnames = PROTECT(DIMNAMES(b, 0)),
561 rdimnames = PROTECT(DIMNAMES(r, 0));
562 if (aside == 'L')
563 matmultDN(rdimnames, adimnames, atrans != 'N', bdimnames, btrans == 'N');
564 else
565 matmultDN(rdimnames, bdimnames, btrans != 'N', adimnames, atrans == 'N');
566 UNPROTECT(3); /* rdimnames, bdimnames, adimnames */
567
568 char aul = UPLO(a);
569 if (triangular < 0)
570 SET_UPLO(r);
571
572 char anu = DIAG(a);
573 if (triangular < -1 || triangular > 1)
574 SET_DIAG(r);
575
576 /* use *mm after copying B into R */
577 /* R := A B */
578 /* R := A' B */
579 /* R := B A */
580 /* R := B A' */
581 /* use *mm after transposing B into R */
582 /* R := A B. */
583 /* R := A' B. */
584 /* R := B. A */
585 /* R := B. A' */
586
587 if (rm > 0 && rn > 0) {
588 SEXP bx = PROTECT(GET_SLOT(b, Matrix_xSym)),
589 rx = PROTECT(Rf_allocVector(TYPEOF(ax), (R_xlen_t) rm * rn));
590 if (TYPEOF(ax) == CPLXSXP) {
591 Rcomplex *pax = COMPLEX(ax), *pbx = COMPLEX(bx), *prx = COMPLEX(rx),
592 one = Matrix_zunit;
593 ztrans2(prx, pbx, (size_t) bm, (size_t) bn, btrans);
594 F77_CALL(ztrmm)(&aside, &aul, &atrans, &anu, &rm, &rn,
595 &one, pax, &rk, prx, &rm FCONE FCONE FCONE FCONE);
596 } else {
597 double *pax = REAL(ax), *pbx = REAL(bx), *prx = REAL(rx),
598 one = 1.0;
599 dtrans2(prx, pbx, (size_t) bm, (size_t) bn, btrans);
600 F77_CALL(dtrmm)(&aside, &aul, &atrans, &anu, &rm, &rn,
601 &one, pax, &rk, prx, &rm FCONE FCONE FCONE FCONE);
602 }
603 SET_SLOT(r, Matrix_xSym, rx);
604 UNPROTECT(2); /* rx, bx */
605 }
606
607 UNPROTECT(2); /* r, ax */
608 return r;
609}
610
611/* op(<,tp>) * op(<,ge>) or op(<,ge>) * op(<,tp>) */
612static
613SEXP tpMatrix_matmult(SEXP a, SEXP b, char atrans, char btrans, char aside,
614 int triangular)
615{
616 int *padim = DIM(a),
617 rk = padim[0];
618
619 int *pbdim = DIM(b), bm = pbdim[0], bn = pbdim[1],
620 rm = (btrans != 'N') ? bn : bm, rn = (btrans != 'N') ? bm : bn;
621
622 if (rk != (((aside == 'L') == (btrans != 'N')) ? bn : bm))
623 Rf_error(_("non-conformable arguments"));
624 if ((int_fast64_t) rm * rn > R_XLEN_T_MAX)
625 Rf_error(_("attempt to allocate vector of length exceeding %s"),
626 "R_XLEN_T_MAX");
627
628 SEXP ax = PROTECT(GET_SLOT(a, Matrix_xSym));
629
630 char rcl[] = "...Matrix";
631 rcl[0] = (TYPEOF(ax) == CPLXSXP) ? 'z' : 'd';
632 rcl[1] = (triangular != 0) ? 't' : 'g';
633 rcl[2] = (triangular != 0) ? 'r' : 'e';
634 SEXP r = PROTECT(newObject(rcl));
635
636 SET_DIM(r, rm, rn);
637
638 SEXP adimnames = PROTECT(DIMNAMES(a, 0)),
639 bdimnames = PROTECT(DIMNAMES(b, 0)),
640 rdimnames = PROTECT(DIMNAMES(r, 0));
641 if (aside == 'L')
642 matmultDN(rdimnames, adimnames, atrans != 'N', bdimnames, btrans == 'N');
643 else
644 matmultDN(rdimnames, bdimnames, btrans != 'N', adimnames, atrans == 'N');
645 UNPROTECT(3); /* rdimnames, bdimnames, adimnames */
646
647 char aul = UPLO(a);
648 if (triangular < 0)
649 SET_UPLO(r);
650
651 char anu = DIAG(a);
652 if (triangular < -1 || triangular > 1)
653 SET_DIAG(r);
654
655 /* use *mv after copying B into R */
656 /* R := A B */
657 /* R := A' B */
658 /* use *mv after transposing B into R */
659 /* R := A B. */
660 /* R := A' B. */
661 /* use *mv after transposing B into R and access R by row */
662 /* R := B. A = (A. B ). */
663 /* R := B. A' = (A'. B ). */
664 /* use *mv after copying B into R and access R by row */
665 /* R := B A = (A* B*)* */
666 /* R := B A' = (A B')' */
667
668 if (rm > 0 && rn > 0) {
669 SEXP bx = PROTECT(GET_SLOT(b, Matrix_xSym)),
670 rx = PROTECT(Rf_allocVector(REALSXP, (R_xlen_t) rm * rn));
671 int i, rinc = (aside == 'L') ? 1 : rm, rincp = (aside == 'L') ? rm : 1;
672
673 char
674 atransp = (aside == 'L')
675 ? atrans : ((atrans != 'N') ? 'N' : ((btrans != 'N') ? btrans : 'T')),
676 btransp = (aside == 'L')
677 ? btrans : ((btrans != 'N') ? 'T' : 'N');
678
679 if (TYPEOF(ax) == CPLXSXP) {
680 Rcomplex *pax = COMPLEX(ax), *pbx = COMPLEX(bx), *prx = COMPLEX(rx);
681 if (aside != 'L' && atrans != 'N' && btrans != 'N' && atrans != btrans)
682 CONJ1(pax, rk);
683 ztrans2(prx, pbx, (size_t) bm, (size_t) bn, btransp);
684 if (aside != 'L' && atrans == 'C' && btrans == 'N')
685 zvconj(prx, NULL, (size_t) rm * (size_t) rn);
686 for (i = 0; i < rn; ++i) {
687 F77_CALL(ztpmv)(&aul, &atransp, &anu, &rk,
688 pax, prx, &rinc FCONE FCONE FCONE);
689 prx += rincp;
690 }
691 if (aside != 'L' &&
692 ((btrans != 'N') ? btrans : ((atrans != 'N') ? atrans : 'T')) == 'C')
693 zvconj(prx, NULL, (size_t) rm * (size_t) rn);
694 } else {
695 double *pax = REAL(ax), *pbx = REAL(bx), *prx = REAL(rx);
696 dtrans2(prx, pbx, (size_t) bm, (size_t) bn, btransp);
697 for (i = 0; i < rn; ++i) {
698 F77_CALL(dtpmv)(&aul, &atransp, &anu, &rk,
699 pax, prx, &rinc FCONE FCONE FCONE);
700 prx += rincp;
701 }
702 }
703 SET_SLOT(r, Matrix_xSym, rx);
704 UNPROTECT(2); /* rx, bx */
705 }
706
707 UNPROTECT(2); /* r, ax */
708 return r;
709}
710
711SEXP R_dense_matmult(SEXP s_x, SEXP s_y,
712 SEXP s_xtrans, SEXP s_ytrans)
713{
714 char
715 xtrans = CHAR(STRING_ELT(s_xtrans, 0))[0],
716 ytrans = CHAR(STRING_ELT(s_ytrans, 0))[0],
717 ztrans = 'N';
718 int m, n, v;
719 matmultDim(s_x, s_y, &xtrans, &ytrans, &ztrans, &m, &n, &v);
720
721 PROTECT_INDEX xpid, ypid;
722 PROTECT_WITH_INDEX(s_x, &xpid);
723 PROTECT_WITH_INDEX(s_y, &ypid);
724
725#define DO_S3(_A_, _TRANS_, _PID_, _ISV_) \
726 do { \
727 if (TYPEOF(_A_) != OBJSXP) { \
728 REPROTECT(_A_ = matrix_as_dense(_A_, ",ge", '\0', '\0', '\0', (_TRANS_ != 'N') ? 0 : 1, 0), _PID_); \
729 if (_ISV_) { \
730 /* Vector: discard names and don't transpose again */ \
731 SET_VECTOR_ELT(DIMNAMES(_A_, 0), \
732 (_TRANS_ != 'N') ? 1 : 0, R_NilValue); \
733 _TRANS_ = 'N'; \
734 } \
735 } \
736 } while (0)
737
738 DO_S3(s_x, xtrans, xpid, v % 2);
739 if (s_y != R_NilValue)
740 DO_S3(s_y, ytrans, ypid, v > 1);
741
742#undef DO_S3
743
744 const char **valid = valid_matmult + 56;
745 const char *xclass = NULL, *yclass = NULL;
746 xclass = Matrix_class(s_x, valid, 6, __func__);
747 if (s_y != R_NilValue)
748 yclass = Matrix_class(s_y, valid, 6, __func__);
749
750 char kind = (xclass[0] == 'z' || (s_y != R_NilValue && yclass[0] == 'z'))
751 ? 'z' : 'd';
752
753#define DO_AS(_A_, _CLASS_, _TRANS_, _PID_) \
754 do { \
755 if (_CLASS_[0] != kind) { \
756 REPROTECT(_A_ = dense_as_kind(_A_, _CLASS_, kind, 0), _PID_); \
757 _CLASS_ = Matrix_class(_A_, valid, 6, __func__); \
758 } \
759 } while (0)
760
761 DO_AS(s_x, xclass, xtrans, xpid);
762 if (s_y != R_NilValue)
763 DO_AS(s_y, yclass, ytrans, ypid);
764
765#undef DO_AS
766
767 if (s_y == R_NilValue) {
768 REPROTECT(s_x = dense_as_general(s_x, xclass, 1), xpid);
769 s_x = geMatrix_matmult(s_x, s_y, xtrans, ytrans);
770 } else if (xclass[1] == 'g' && yclass[1] == 'g') {
771 s_x = geMatrix_matmult(s_x, s_y, xtrans, ytrans);
772 } else if (xclass[1] == 'g' || yclass[1] == 'g') {
773 s_x = (xclass[1] == 'g')
774 ? ((yclass[1] == 's')
775 ? ((yclass[2] != 'p')
776 ? syMatrix_matmult(s_y, s_x, ytrans, xtrans, 'R')
777 : spMatrix_matmult(s_y, s_x, ytrans, xtrans, 'R'))
778 : ((yclass[2] != 'p')
779 ? trMatrix_matmult(s_y, s_x, ytrans, xtrans, 'R', 0)
780 : tpMatrix_matmult(s_y, s_x, ytrans, xtrans, 'R', 0)))
781 : ((xclass[1] == 's')
782 ? ((xclass[2] != 'p')
783 ? syMatrix_matmult(s_x, s_y, xtrans, ytrans, 'L')
784 : spMatrix_matmult(s_x, s_y, xtrans, ytrans, 'L'))
785 : ((xclass[2] != 'p')
786 ? trMatrix_matmult(s_x, s_y, xtrans, ytrans, 'L', 0)
787 : tpMatrix_matmult(s_x, s_y, xtrans, ytrans, 'L', 0)));
788 } else if (xclass[1] == 's' && yclass[1] == 's') {
789 if (xclass[2] == 'p' && yclass[2] == 'p') {
790 REPROTECT(s_y = dense_as_general(s_y, yclass, 1), ypid);
791 s_x = spMatrix_matmult(s_x, s_y, xtrans, ytrans, 'L');
792 } else if (xclass[2] == 'p') {
793 REPROTECT(s_x = dense_as_general(s_x, xclass, 1), xpid);
794 s_x = syMatrix_matmult(s_y, s_x, ytrans, xtrans, 'R');
795 } else {
796 REPROTECT(s_y = dense_as_general(s_y, yclass, 1), ypid);
797 s_x = syMatrix_matmult(s_x, s_y, xtrans, ytrans, 'L');
798 }
799 } else if (xclass[1] == 's' || yclass[1] == 's') {
800 if (xclass[1] == 's') {
801 REPROTECT(s_x = dense_as_general(s_x, xclass, 1), xpid);
802 s_x = (yclass[2] != 'p')
803 ? trMatrix_matmult(s_y, s_x, ytrans, xtrans, 'R', 0)
804 : tpMatrix_matmult(s_y, s_x, ytrans, xtrans, 'R', 0);
805 } else {
806 REPROTECT(s_y = dense_as_general(s_y, yclass, 1), ypid);
807 s_x = (xclass[2] != 'p')
808 ? trMatrix_matmult(s_x, s_y, xtrans, ytrans, 'L', 0)
809 : tpMatrix_matmult(s_x, s_y, xtrans, ytrans, 'L', 0);
810 }
811 } else {
812 int triangular = 0;
813
814#define DO_TR \
815 do { \
816 char xul = UPLO(s_x), xnu = DIAG(s_x), \
817 yul = UPLO(s_y), ynu = DIAG(s_y); \
818 if (xtrans != 'N') \
819 xul = (xul == 'U') ? 'L' : 'U'; \
820 if (ytrans != 'N') \
821 yul = (yul == 'U') ? 'L' : 'U'; \
822 triangular = (xul != yul) ? 0 : ((xnu != ynu || xnu == 'N') ? 1 : 2); \
823 if (xul != 'U') \
824 triangular = -triangular; \
825 } while (0)
826
827 DO_TR;
828
829 if (xclass[2] == 'p' && yclass[2] == 'p') {
830 REPROTECT(s_y = dense_as_general(s_y, yclass, 1), ypid);
831 s_x = tpMatrix_matmult(s_x, s_y, xtrans, ytrans, 'L', triangular);
832 } else if (xclass[2] == 'p') {
833 REPROTECT(s_x = dense_as_general(s_x, xclass, 1), xpid);
834 s_x = trMatrix_matmult(s_y, s_x, ytrans, xtrans, 'R', triangular);
835 } else {
836 REPROTECT(s_y = dense_as_general(s_y, yclass, 1), ypid);
837 s_x = trMatrix_matmult(s_x, s_y, xtrans, ytrans, 'L', triangular);
838 }
839 }
840
841 UNPROTECT(2); /* s_y, s_x */
842 return s_x;
843}
844
845/* boolean: op(op(<.gC>) & op(<.gC>)) */
846/* numeric: op(op(<,gC>) * op(<,gC>)) */
847static
848SEXP gCgCMatrix_matmult(SEXP x, SEXP y, char xtrans, char ytrans, char ztrans,
849 int triangular, int boolean)
850{
851
852#define ASMODE(_TRANS_) ((boolean) ? 0 : (((_TRANS_) != 'C') ? 1 : 2))
853
854 cholmod_sparse *X = M2CHS(x, !boolean), *Y = NULL, *Z = NULL;
855
856 SEXP z;
857 char zclass[] = "..CMatrix";
858
859 if (y == R_NilValue) {
860
861 zclass[0] = (boolean) ? 'n' : ((X->xtype != CHOLMOD_COMPLEX) ? 'd' : 'z');
862 zclass[1] = (boolean) ? 's' : ((X->xtype != CHOLMOD_COMPLEX) ? 'p' : ((((xtrans != 'N') ? xtrans : ytrans) == 'C') ? 'p' : 's'));
863
864 if (xtrans != 'N')
865 X = cholmod_transpose(X, ASMODE(xtrans), &c);
866 Z = cholmod_aat(X, NULL, 0, ASMODE((xtrans != 'N') ? xtrans : ytrans), &c);
867 if (xtrans != 'N')
868 cholmod_free_sparse(&X, &c);
869 if (!Z->sorted)
870 cholmod_sort(Z, &c);
871 X = cholmod_copy(Z, (ztrans != 'N') ? -1 : 1, !boolean, &c);
872 cholmod_free_sparse(&Z, &c);
873 Z = X;
874 PROTECT(z = CHS2M(Z, !boolean, zclass[1]));
875 cholmod_free_sparse(&Z, &c);
876
877 SEXP xdimnames = PROTECT(DIMNAMES(x, 0)),
878 zdimnames = PROTECT(DIMNAMES(z, 0));
879 symDN(zdimnames, xdimnames, (xtrans != 'N') ? 1 : 0);
880 UNPROTECT(2); /* zdimnames, xdimnames */
881
882 if (ztrans != 'N')
883 SET_UPLO(z);
884 if (zclass[1] == 's' && zclass[0] == 'z')
885 SET_TRANS(z);
886
887 } else {
888
889 Y = M2CHS(y, !boolean);
890
891 if (((xtrans != 'N') ? X->nrow : X->ncol) !=
892 ((ytrans != 'N') ? Y->ncol : Y->nrow))
893 Rf_error(_("non-conformable arguments"));
894
895 zclass[0] = (boolean) ? 'n' : ((X->xtype != CHOLMOD_COMPLEX && Y->xtype != CHOLMOD_COMPLEX) ? 'd' : 'z');
896 zclass[1] = (triangular != 0) ? 't' : 'g';
897
898 if (xtrans != 'N')
899 X = cholmod_transpose(X, ASMODE(xtrans), &c);
900 if (ytrans != 'N')
901 Y = cholmod_transpose(Y, ASMODE(ytrans), &c);
902 Z = cholmod_ssmult(X, Y, 0, !boolean, 1, &c);
903 if (xtrans != 'N')
904 cholmod_free_sparse(&X, &c);
905 if (ytrans != 'N')
906 cholmod_free_sparse(&Y, &c);
907 if (triangular < -1)
908 cholmod_band_inplace(-((int) Z->nrow), -1, !boolean, Z, &c);
909 else if (triangular > 1)
910 cholmod_band_inplace(1, (int) Z->ncol, !boolean, Z, &c);
911 PROTECT(z = CHS2M(Z, !boolean, zclass[1]));
912 cholmod_free_sparse(&Z, &c);
913
914 SEXP xdimnames = PROTECT(DIMNAMES(x, 0)),
915 ydimnames = PROTECT(DIMNAMES(y, 0)),
916 zdimnames = PROTECT(DIMNAMES(z, 0));
917 matmultDN(zdimnames,
918 xdimnames, (xtrans != 'N') ? 1 : 0,
919 ydimnames, (ytrans != 'N') ? 0 : 1);
920 UNPROTECT(3); /* zdimnames, ydimnames, xdimnames */
921
922 if (triangular < 0)
923 SET_UPLO(z);
924 if (triangular < -1 || triangular > 1)
925 SET_DIAG(z);
926
927 }
928
929 if (ztrans != 'N')
930 z = sparse_transpose(z, zclass, ztrans, 1);
931
932 UNPROTECT(1); /* z */
933 return z;
934}
935
936/* op(op(<,[gs]C>) * op(<,ge>)) */
937/* NB: currently, caller must coerce complex symmetric zsC to zgC */
938static
939SEXP gCgeMatrix_matmult(SEXP x, SEXP y, int xtrans, char ytrans, char ztrans,
940 int triangular, int symmetric)
941{
942 cholmod_sparse *X = M2CHS(x, 1);
943 X->stype = symmetric;
944
945 cholmod_dense *Y = M2CHD(y, ytrans);
946
947 if (((xtrans != 'N') ? X->nrow : X->ncol) != Y->nrow)
948 Rf_error(_("non-conformable arguments"));
949 size_t m = (xtrans != 'N') ? X->ncol : X->nrow, n = Y->ncol;
950 if (m * n > R_XLEN_T_MAX)
951 Rf_error(_("attempt to allocate vector of length exceeding %s"),
952 "R_XLEN_T_MAX");
953
954 if (X->xtype == CHOLMOD_COMPLEX && xtrans == 'T')
955 CONJ2(X->x, X->nrow, X->ncol);
956
957 char zclass[] = "...Matrix";
958 zclass[0] = (X->xtype != CHOLMOD_COMPLEX && Y->xtype != CHOLMOD_COMPLEX) ? 'd' : 'z';
959 zclass[1] = (triangular != 0) ? 't' : 'g';
960 zclass[2] = (triangular != 0) ? 'r' : 'e';
961 SEXP z = PROTECT(newObject(zclass));
962
963 SET_DIM(z,
964 (int) ((ztrans != 'N') ? n : m),
965 (int) ((ztrans != 'N') ? m : n));
966
967 SEXP xdimnames = PROTECT(DIMNAMES(x, -(symmetric != 0))),
968 ydimnames = PROTECT(DIMNAMES(y, 0)),
969 zdimnames = PROTECT(DIMNAMES(z, 0));
970 if (ztrans != 'N')
971 matmultDN(zdimnames,
972 ydimnames, (ytrans != 'N') ? 0 : 1,
973 xdimnames, (xtrans != 'N') ? 1 : 0);
974 else
975 matmultDN(zdimnames,
976 xdimnames, (xtrans != 'N') ? 1 : 0,
977 ydimnames, (ytrans != 'N') ? 0 : 1);
978 UNPROTECT(3); /* zdimnames, ydimnames, xdimnames */
979
980 if (triangular != 0 && (ztrans != 'N') == (triangular > 0))
981 SET_UPLO(z);
982 if (triangular < -1 || triangular > 1)
983 SET_DIAG(z);
984
985 double alpha[2] = { 1.0, 0.0 }, beta[2] = { 0.0, 0.0 };
986 cholmod_dense *Z = (cholmod_dense *) R_alloc(1, sizeof(cholmod_dense));
987 memset(Z, 0, sizeof(cholmod_dense));
988 Z->nrow = m;
989 Z->ncol = n;
990 Z->d = Z->nrow;
991 Z->nzmax = Z->nrow * Z->ncol;
992 Z->xtype = X->xtype;
993 Z->dtype = X->dtype;
994
995 SEXP zx;
996 if (zclass[0] == 'z') {
997 PROTECT(zx = Rf_allocVector(CPLXSXP, (R_xlen_t) (m * n)));
998 Z->x = (ztrans != 'N')
999 ? (Rcomplex *) R_alloc(m * n, sizeof(Rcomplex))
1000 : COMPLEX(zx);
1001 } else {
1002 PROTECT(zx = Rf_allocVector(REALSXP, (R_xlen_t) (m * n)));
1003 Z->x = (ztrans != 'N')
1004 ? (double *) R_alloc(m * n, sizeof(double))
1005 : REAL(zx);
1006 }
1007 cholmod_sdmult(X, xtrans != 'N', alpha, beta, Y, Z, &c);
1008 if (ztrans != 'N') {
1009 if (zclass[0] == 'z')
1010 ztrans2(COMPLEX(zx), (Rcomplex *) Z->x, m, n, ztrans);
1011 else
1012 dtrans2( REAL(zx), ( double *) Z->x, m, n, ztrans);
1013 }
1014 SET_SLOT(z, Matrix_xSym, zx);
1015 UNPROTECT(1); /* zx */
1016
1017 UNPROTECT(1); /* z */
1018 return z;
1019}
1020
1021SEXP R_sparse_matmult(SEXP s_x, SEXP s_y,
1022 SEXP s_xtrans, SEXP s_ytrans, SEXP s_ztrans,
1023 SEXP s_boolean)
1024{
1025 if (TYPEOF(s_boolean) != LGLSXP || LENGTH(s_boolean) < 1)
1026 Rf_error(_("invalid '%s' to '%s'"), "boolean", __func__);
1027 int boolean = LOGICAL(s_boolean)[0];
1028
1029 char
1030 xtrans = CHAR(STRING_ELT(s_xtrans, 0))[0],
1031 ytrans = CHAR(STRING_ELT(s_ytrans, 0))[0],
1032 ztrans = CHAR(STRING_ELT(s_ztrans, 0))[0];
1033 int m, n, v;
1034 matmultDim(s_x, s_y, &xtrans, &ytrans, &ztrans, &m, &n, &v);
1035
1036 PROTECT_INDEX xpid, ypid;
1037 PROTECT_WITH_INDEX(s_x, &xpid);
1038 PROTECT_WITH_INDEX(s_y, &ypid);
1039
1040#define DO_S3(_A_, _TRANS_, _PID_, _ISV_) \
1041 do { \
1042 if (TYPEOF(_A_) != OBJSXP) { \
1043 if (boolean == NA_LOGICAL || !boolean) \
1044 REPROTECT(_A_ = matrix_as_dense (_A_, ",ge", '\0', '\0', '\0', (_TRANS_ != 'N') ? 0 : 1, 0), _PID_); \
1045 else if (_TRANS_ != 'N') \
1046 REPROTECT(_A_ = matrix_as_sparse(_A_, "ngR", '\0', '\0', '\0', 0), _PID_); \
1047 else \
1048 REPROTECT(_A_ = matrix_as_sparse(_A_, "ngC", '\0', '\0', '\0', 0), _PID_); \
1049 if (_ISV_) { \
1050 /* Discard names and don't transpose again */ \
1051 SET_VECTOR_ELT(DIMNAMES(_A_, 0), \
1052 (_TRANS_ != 'N') ? 1 : 0, R_NilValue); \
1053 _TRANS_ = 'N'; \
1054 } \
1055 } \
1056 } while (0)
1057
1058 DO_S3(s_x, xtrans, xpid, v % 2);
1059 if (s_y != R_NilValue)
1060 DO_S3(s_y, ytrans, ypid, v > 1);
1061
1062#undef DO_S3
1063
1064 const char **valid = valid_matmult + 5;
1065 const char *xclass = NULL, *yclass = NULL;
1066 xclass = Matrix_class(s_x, valid, 6, __func__);
1067 if (s_y != R_NilValue)
1068 yclass = Matrix_class(s_y, valid, 6, __func__);
1069
1070 if (boolean == NA_LOGICAL)
1071 boolean = xclass[0] == 'n' && (s_y == R_NilValue || yclass[0] == 'n');
1072 char kind = (boolean) ? 'n' :
1073 ((xclass[0] != 'z' && (s_y == R_NilValue || yclass[0] != 'z')) ? 'd' : 'z');
1074
1075#define DO_AS(_A_, _CLASS_, _TRANS_, _PID_) \
1076 do { \
1077 if (_CLASS_[2] != 'C' && _TRANS_ != 'N') { \
1078 if (_CLASS_[2] != 'R' && _CLASS_[2] != 'T') { \
1079 REPROTECT(_A_ = dense_as_sparse(_A_, _CLASS_, 'R'), _PID_); \
1080 _CLASS_ = Matrix_class(_A_, valid, 6, __func__); \
1081 } \
1082 REPROTECT(_A_ = sparse_transpose(_A_, _CLASS_, _TRANS_, 1), _PID_); \
1083 _CLASS_ = Matrix_class(_A_, valid, 6, __func__); \
1084 _TRANS_ = 'N'; \
1085 } \
1086 if (_CLASS_[2] != 'C') { \
1087 if (_CLASS_[2] != 'R' && _CLASS_[2] != 'T') \
1088 REPROTECT(_A_ = dense_as_sparse(_A_, _CLASS_, 'C'), _PID_); \
1089 else \
1090 REPROTECT(_A_ = sparse_as_Csparse(_A_, _CLASS_), _PID_); \
1091 _CLASS_ = Matrix_class(_A_, valid, 6, __func__); \
1092 } \
1093 if (_TRANS_ != 'N' && _CLASS_[1] == 's' && \
1094 (_CLASS_[0] != 'z' || TRANS(_A_) == _TRANS_)) \
1095 _TRANS_ = 'N'; \
1096 if (_CLASS_[0] != kind) { \
1097 if (boolean) \
1098 REPROTECT(_A_ = sparse_dropzero(_A_, _CLASS_, 0.0), _PID_); \
1099 else { \
1100 REPROTECT(_A_ = sparse_as_kind(_A_, _CLASS_, kind), _PID_); \
1101 _CLASS_ = Matrix_class(_A_, valid, 6, __func__); \
1102 } \
1103 } \
1104 } while (0)
1105
1106 DO_AS(s_x, xclass, xtrans, xpid);
1107
1108 if (s_y == R_NilValue) {
1109 REPROTECT(s_x = sparse_as_general(s_x, xclass), xpid);
1110 s_x = gCgCMatrix_matmult(
1111 s_x, s_y, xtrans, ytrans, ztrans, 0, boolean);
1112 UNPROTECT(2); /* s_y, s_x */
1113 return s_x;
1114 }
1115
1116 int triangular = 0;
1117 if (xclass[1] == 't' && yclass[1] == 't')
1118 DO_TR;
1119
1120 if (!boolean && yclass[2] != 'C' && yclass[2] != 'R' && yclass[2] != 'T') {
1121 if (xclass[1] == 's' && xclass[0] == 'z' && TRANS(s_x) != 'C') {
1122 REPROTECT(s_x = sparse_as_general(s_x, xclass), xpid);
1123 xclass = Matrix_class(s_x, valid, 6, __func__);
1124 }
1125 int symmetric = xclass[1] == 's';
1126 if (symmetric && UPLO(s_x) != 'U')
1127 symmetric = -symmetric;
1128 if (yclass[0] != kind) {
1129 REPROTECT(s_y = dense_as_kind(s_y, yclass, kind, 0), ypid);
1130 yclass = Matrix_class(s_y, valid, 6, __func__);
1131 }
1132 REPROTECT(s_x = sparse_diag_U2N(s_x, xclass), xpid);
1133 REPROTECT(s_y = dense_as_general(s_y, yclass, 1), ypid);
1134 s_x = gCgeMatrix_matmult(
1135 s_x, s_y, xtrans, ytrans, ztrans, triangular, symmetric);
1136 UNPROTECT(2); /* s_y, s_x */
1137 return s_x;
1138 }
1139
1140 DO_AS(s_y, yclass, ytrans, ypid);
1141
1142#undef DO_AS
1143
1144 REPROTECT(s_x = sparse_as_general(s_x, xclass), xpid);
1145 REPROTECT(s_y = sparse_as_general(s_y, yclass), ypid);
1146 s_x = gCgCMatrix_matmult(
1147 s_x, s_y, xtrans, ytrans, ztrans, triangular, boolean);
1148 UNPROTECT(2); /* s_y, s_x */
1149 return s_x;
1150}
1151
1152#define zSCALE(x, y) \
1153 do { \
1154 Rcomplex tmp = (x); \
1155 (x).r = tmp.r * (y).r - tmp.i * (y).i; \
1156 (x).i = tmp.r * (y).i + tmp.i * (y).r; \
1157 } while (0)
1158#define dSCALE(x, y) \
1159 (x) = (x) * (y)
1160#define lSCALE(x, y) \
1161 (x) = (x) && (y)
1162
1163#define SCALE3(t, index) \
1164 do { \
1165 switch (t) { \
1166 case CPLXSXP: SCALE(z, index); break; \
1167 case REALSXP: SCALE(d, index); break; \
1168 case LGLSXP: SCALE(l, index); break; \
1169 default: break; \
1170 } \
1171 } while (0)
1172
1173static
1174void dense_colscale(SEXP obj, SEXP d, int m, int n, char ul, char nu)
1175{
1176 SEXP x = GET_SLOT(obj, Matrix_xSym);
1177 int i, j, packed = XLENGTH(x) != (int_fast64_t) m * n;
1178
1179#define SCALE(c, index) \
1180 do { \
1181 c##TYPE *px = c##PTR(x), *pd = c##PTR(d); \
1182 if (ul == '\0') { \
1183 for (j = 0; j < n; ++j) { \
1184 for (i = 0; i < m; ++i) { \
1185 c##SCALE(*px, pd[index]); \
1186 ++px; \
1187 } \
1188 } \
1189 } else if (ul == 'U') { \
1190 for (j = 0; j < n; ++j) { \
1191 for (i = 0; i <= j; ++i) { \
1192 c##SCALE(*px, pd[index]); \
1193 ++px; \
1194 } \
1195 if (!packed) \
1196 px += n - j - 1; \
1197 } \
1198 } else { \
1199 for (j = 0; j < n; ++j) { \
1200 if (!packed) \
1201 px += j; \
1202 for (i = j; i < n; ++i) { \
1203 c##SCALE(*px, pd[index]); \
1204 ++px; \
1205 } \
1206 } \
1207 } \
1208 if (nu != '\0' && nu != 'N') { \
1209 size_t n_ = (size_t) n; \
1210 if (!packed) \
1211 c##NAME(copy2)(n_, c##PTR(x), n_ + 1, pd, 1); \
1212 else if (ul == 'U') \
1213 c##NAME(copy1)(n_, c##PTR(x), 2 , 1, 0, pd, 1, 0, 0); \
1214 else \
1215 c##NAME(copy1)(n_, c##PTR(x), n_, 1, 1, pd, 1, 0, 0); \
1216 } \
1217 } while (0)
1218
1219 SCALE3(TYPEOF(d), j);
1220 return;
1221}
1222
1223static
1224void dense_rowscale(SEXP obj, SEXP d, int m, int n, char ul, char nu)
1225{
1226 SEXP x = GET_SLOT(obj, Matrix_xSym);
1227 int i, j, packed = XLENGTH(x) != (int_fast64_t) m * n;
1228 SCALE3(TYPEOF(d), i);
1229
1230#undef SCALE
1231
1232 return;
1233}
1234
1235/* boolean: <lgC> & <ldi> or <ldi> & <lgR> */
1236/* numeric: <,gC> * <,di> or <,di> * <,gR> */
1237static
1238void Csparse_colscale(SEXP obj, SEXP d)
1239{
1240 SEXP x = PROTECT(GET_SLOT(obj, Matrix_xSym)),
1241 p = PROTECT(GET_SLOT(obj, Matrix_pSym));
1242 int *pp = INTEGER(p) + 1, n = (int) (XLENGTH(p) - 1), j, k, kend;
1243 UNPROTECT(2); /* p, x */
1244
1245#define SCALE(c, index) \
1246 do { \
1247 c##TYPE *px = c##PTR(x), *pd = c##PTR(d); \
1248 for (j = 0, k = 0; j < n; ++j) { \
1249 kend = pp[j]; \
1250 while (k < kend) { \
1251 c##SCALE(*px, *pd); \
1252 ++px; \
1253 ++k; \
1254 } \
1255 ++pd; \
1256 } \
1257 } while (0)
1258
1259 SCALE3(TYPEOF(d), );
1260
1261#undef SCALE
1262
1263 return;
1264}
1265
1266/* boolean: <ldi> & <lgC> or <lgR> & <ldi> */
1267/* numeric: <,di> * <,gC> or <,gR> * <,di> */
1268static
1269void Csparse_rowscale(SEXP obj, SEXP d, SEXP iSym)
1270{
1271 SEXP x = PROTECT(GET_SLOT(obj, Matrix_xSym)),
1272 p = PROTECT(GET_SLOT(obj, Matrix_pSym)),
1273 i = PROTECT(GET_SLOT(obj, iSym));
1274 int *pi = INTEGER(i), k, kend = INTEGER(p)[XLENGTH(p) - 1];
1275 UNPROTECT(3); /* i, p, x */
1276
1277#define SCALE(c, index) \
1278 do { \
1279 c##TYPE *px = c##PTR(x), *pd = c##PTR(d); \
1280 for (k = 0; k < kend; ++k) { \
1281 c##SCALE(*px, pd[*pi]); \
1282 ++px; \
1283 ++pi; \
1284 } \
1285 } while (0)
1286
1287 SCALE3(TYPEOF(d), );
1288 return;
1289}
1290
1291/* boolean: <ldi> & <lgT> or <lgT> & <ldi> */
1292/* numeric: <,di> * <,gT> or <,gT> * <,di> */
1293static
1294void Tsparse_rowscale(SEXP obj, SEXP d, SEXP iSym)
1295{
1296 SEXP x = PROTECT(GET_SLOT(obj, Matrix_xSym)),
1297 i = PROTECT(GET_SLOT(obj, iSym));
1298 int *pi = INTEGER(i);
1299 R_xlen_t k, kend = XLENGTH(i);
1300 UNPROTECT(2); /* i, x */
1301 SCALE3(TYPEOF(d), );
1302
1303#undef SCALE
1304
1305 return;
1306}
1307
1308#undef SCALE3
1309
1310SEXP R_diagonal_matmult(SEXP s_x, SEXP s_y,
1311 SEXP s_xtrans, SEXP s_ytrans,
1312 SEXP s_boolean)
1313{
1314 SEXP _s_x = s_x, _s_y = s_y; /* for later pointer comparison */
1315
1316 if (TYPEOF(s_boolean) != LGLSXP || LENGTH(s_boolean) < 1)
1317 Rf_error(_("invalid '%s' to '%s'"), "boolean", __func__);
1318 int boolean = LOGICAL(s_boolean)[0];
1319
1320 char
1321 xtrans = CHAR(STRING_ELT(s_xtrans, 0))[0],
1322 ytrans = CHAR(STRING_ELT(s_ytrans, 0))[0],
1323 ztrans = 'N';
1324 int m, n, v;
1325 matmultDim(s_x, s_y, &xtrans, &ytrans, &ztrans, &m, &n, &v);
1326
1327 PROTECT_INDEX xpid, ypid;
1328 PROTECT_WITH_INDEX(s_x, &xpid);
1329 PROTECT_WITH_INDEX(s_y, &ypid);
1330
1331#define DO_S3(_A_, _TRANS_, _PID_, _ISV_) \
1332 do { \
1333 if (TYPEOF(_A_) != OBJSXP) { \
1334 if (boolean == NA_LOGICAL || !boolean) \
1335 REPROTECT(_A_ = matrix_as_dense(_A_, ",ge", '\0', '\0', '\0', (_TRANS_ != 'N') ? 0 : 1, 2), _PID_); \
1336 else \
1337 REPROTECT(_A_ = matrix_as_dense(_A_, "nge", '\0', '\0', '\0', (_TRANS_ != 'N') ? 0 : 1, 2), _PID_); \
1338 if (_ISV_) { \
1339 /* Vector: discard names and don't transpose again */ \
1340 SET_VECTOR_ELT(DIMNAMES(_A_, 0), \
1341 (_TRANS_ != 'N') ? 1 : 0, R_NilValue); \
1342 _TRANS_ = 'N'; \
1343 } \
1344 } \
1345 } while (0)
1346
1347 DO_S3(s_x, xtrans, xpid, v % 2);
1348 if (s_y != R_NilValue)
1349 DO_S3(s_y, ytrans, ypid, v > 1);
1350
1351#undef DO_S3
1352
1353 const char **valid = valid_matmult;
1354 const char *xclass = NULL, *yclass = NULL;
1355 xclass = Matrix_class(s_x, valid, 6, __func__);
1356 if (s_y != R_NilValue)
1357 yclass = Matrix_class(s_y, valid, 6, __func__);
1358
1359 if (boolean == NA_LOGICAL)
1360 boolean = xclass[0] == 'n' && yclass[0] == 'n';
1361 char kind = (boolean) ? 'n' :
1362 ((xclass[0] != 'z' && yclass[0] != 'z') ? 'd' : 'z');
1363 char ks = (boolean) ? 'l' : kind, kd = kind;
1364
1365 int mg = -1, id = -1;
1366 if (xclass[2] == 'i') {
1367 mg = 0;
1368 id = DIAG(s_x) != 'N';
1369 } else if (yclass[2] == 'i') {
1370 mg = 1;
1371 id = DIAG(s_y) != 'N';
1372 } else
1373 Rf_error(_("should never happen ..."));
1374
1375#define DO_AS(_A_, _CLASS_, _TRANS_, _PID_) \
1376 do { \
1377 if (_TRANS_ != 'N' && _CLASS_[1] == 's' && \
1378 (_CLASS_[0] != 'z' || TRANS(_A_) == _TRANS_)) \
1379 _TRANS_ = 'N'; \
1380 switch (_CLASS_[2]) { \
1381 case 'i': \
1382 if (!id && _CLASS_[0] != ks) { \
1383 REPROTECT(_A_ = diagonal_as_kind(_A_, _CLASS_, ks), _PID_); \
1384 _CLASS_ = Matrix_class(_A_, valid, 6, __func__); \
1385 } \
1386 break; \
1387 case 'C': \
1388 case 'R': \
1389 case 'T': \
1390 if (_CLASS_[0] != ks) { \
1391 REPROTECT(_A_ = sparse_as_kind(_A_, _CLASS_, ks), _PID_); \
1392 _CLASS_ = Matrix_class(_A_, valid, 6, __func__); \
1393 } \
1394 if (!id && _CLASS_[1] == 's') { \
1395 REPROTECT(_A_ = sparse_as_general(_A_, _CLASS_), _PID_); \
1396 _CLASS_ = Matrix_class(_A_, valid, 6, __func__); \
1397 } \
1398 if (!id && _CLASS_[1] == 't') \
1399 REPROTECT(_A_ = sparse_diag_U2N(_A_, _CLASS_), _PID_); \
1400 if (_TRANS_ != 'N') { \
1401 REPROTECT(_A_ = sparse_transpose(_A_, _CLASS_, _TRANS_, 0), _PID_); \
1402 _TRANS_ = 'N'; \
1403 } \
1404 break; \
1405 default: \
1406 if (_CLASS_[0] != kd) { \
1407 REPROTECT(_A_ = dense_as_kind(_A_, _CLASS_, kd, 1), _PID_); \
1408 _CLASS_ = Matrix_class(_A_, valid, 6, __func__); \
1409 } \
1410 if (!id && _CLASS_[1] == 's') { \
1411 REPROTECT(_A_ = dense_as_general(_A_, _CLASS_, _A_ == _ ## _A_), _PID_); \
1412 _CLASS_ = Matrix_class(_A_, valid, 6, __func__); \
1413 } \
1414 if (_TRANS_ != 'N') { \
1415 REPROTECT(_A_ = dense_transpose(_A_, _CLASS_, _TRANS_), _PID_); \
1416 _TRANS_ = 'N'; \
1417 } \
1418 break; \
1419 } \
1420 } while (0)
1421
1422 DO_AS(s_x, xclass, xtrans, xpid);
1423 DO_AS(s_y, yclass, ytrans, ypid);
1424
1425 PROTECT_INDEX zpid;
1426 const char *zclass = (mg == 0) ? yclass : xclass;
1427 SEXP z = newObject(zclass);
1428 PROTECT_WITH_INDEX(z, &zpid);
1429
1430 SET_DIM(z, m, n);
1431
1432 SEXP xdimnames = PROTECT(DIMNAMES(s_x, 0)),
1433 ydimnames = PROTECT(DIMNAMES(s_y, 0)),
1434 zdimnames = PROTECT(DIMNAMES(z, 0));
1435 matmultDN(zdimnames,
1436 xdimnames, (xtrans != 'N') ? 1 : 0,
1437 ydimnames, (ytrans != 'N') ? 0 : 1);
1438 UNPROTECT(3); /* zdimnames, ydimnames, xdimnames */
1439
1440 char ul = '\0', nu = '\0';
1441 if (zclass[1] != 'g' && (ul = UPLO((mg == 0) ? s_y : s_x)) != 'U')
1442 SET_UPLO(z);
1443 if (zclass[1] == 't' && (nu = DIAG((mg == 0) ? s_y : s_x)) != 'N' && id)
1444 SET_DIAG(z);
1445
1446 if (zclass[2] == 'C' || zclass[2] == 'R' || zclass[2] == 'T') {
1447 if (zclass[2] != 'T') {
1448 SEXP p = PROTECT(GET_SLOT((mg == 0) ? s_y : s_x, Matrix_pSym));
1449 SET_SLOT(z, Matrix_pSym, p);
1450 UNPROTECT(1); /* p */
1451 }
1452 if (zclass[2] != 'R') {
1453 SEXP i = PROTECT(GET_SLOT((mg == 0) ? s_y : s_x, Matrix_iSym));
1454 SET_SLOT(z, Matrix_iSym, i);
1455 UNPROTECT(1); /* i */
1456 }
1457 if (zclass[2] != 'C') {
1458 SEXP j = PROTECT(GET_SLOT((mg == 0) ? s_y : s_x, Matrix_jSym));
1459 SET_SLOT(z, Matrix_jSym, j);
1460 UNPROTECT(1); /* j */
1461 }
1462 }
1463
1464 SEXP x0 = PROTECT(GET_SLOT((mg == 0) ? s_y : s_x, Matrix_xSym));
1465 if (id || ((mg == 0) ? s_y != _s_y : s_x != _s_x))
1466 SET_SLOT(z, Matrix_xSym, x0);
1467 else {
1468 SEXP x1 = PROTECT(Rf_allocVector(TYPEOF(x0), XLENGTH(x0)));
1469 switch (kind) {
1470 case 'z':
1471 memcpy(COMPLEX(x1), COMPLEX(x0), sizeof(Rcomplex) * (size_t) XLENGTH(x0));
1472 break;
1473 case 'd':
1474 memcpy( REAL(x1), REAL(x0), sizeof( double) * (size_t) XLENGTH(x0));
1475 break;
1476 default:
1477 memcpy(LOGICAL(x1), LOGICAL(x0), sizeof( int) * (size_t) XLENGTH(x0));
1478 break;
1479 }
1480 SET_SLOT(z, Matrix_xSym, x1);
1481 UNPROTECT(1); /* x1 */
1482 }
1483 UNPROTECT(1); /* x0 */
1484
1485 if (!id) {
1486 SEXP d = PROTECT(GET_SLOT((mg == 0) ? s_x : s_y, Matrix_xSym));
1487 switch (zclass[2]) {
1488 case 'C':
1489 if (mg == 0)
1490 Csparse_rowscale(z, d, Matrix_iSym);
1491 else
1492 Csparse_colscale(z, d);
1493 break;
1494 case 'R':
1495 if (mg == 0)
1496 Csparse_colscale(z, d);
1497 else
1498 Csparse_rowscale(z, d, Matrix_jSym);
1499 break;
1500 case 'T':
1501 if (mg == 0)
1502 Tsparse_rowscale(z, d, Matrix_iSym);
1503 else
1504 Tsparse_rowscale(z, d, Matrix_jSym);
1505 break;
1506 default:
1507 if (mg == 0)
1508 dense_rowscale(z, d, m, n, ul, nu);
1509 else
1510 dense_colscale(z, d, m, n, ul, nu);
1511 break;
1512 }
1513 UNPROTECT(1); /* d */
1514 }
1515
1516 if (boolean && (zclass[2] == 'C' || zclass[2] == 'R' || zclass[2] == 'T')) {
1517 REPROTECT(z = sparse_dropzero(z, zclass, 0.0), zpid);
1518 REPROTECT(z = sparse_as_kind(z, zclass, 'n'), zpid);
1519 }
1520
1521 UNPROTECT(3); /* z, s_y, s_x */
1522 return z;
1523}
dense_as_kind
SEXP dense_as_kind(SEXP, const char *, char, int)
Definition coerce.c:2000
sparse_as_general
SEXP sparse_as_general(SEXP, const char *)
Definition coerce.c:2298
sparse_as_kind
SEXP sparse_as_kind(SEXP, const char *, char)
Definition coerce.c:2076
FCONE
#define FCONE
Definition Lapack-etc.h:13
SET_DIM
#define SET_DIM(x, m, n)
Definition Mdefines.h:87
SWAP
#define SWAP(a, b, t, op)
Definition Mdefines.h:138
VALID_DENSE
#define VALID_DENSE
Definition Mdefines.h:230
VALID_SPARSE
#define VALID_SPARSE
Definition Mdefines.h:257
_
#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
Matrix_class
const char * Matrix_class(SEXP, const char **, int, const char *)
Definition objects.c:112
DIMNAMES
#define DIMNAMES(x, mode)
Definition Mdefines.h:96
SET_TRANS
#define SET_TRANS(x)
Definition Mdefines.h:108
VALID_DIAGONAL
#define VALID_DIAGONAL
Definition Mdefines.h:260
TRANS
#define TRANS(x)
Definition Mdefines.h:106
SET_DIAG
#define SET_DIAG(x)
Definition Mdefines.h:113
DIM
#define DIM(x)
Definition Mdefines.h:85
GET_SLOT
#define GET_SLOT(x, name)
Definition Mdefines.h:72
newObject
SEXP newObject(const char *)
Definition objects.c:13
SET_SLOT
#define SET_SLOT(x, name, value)
Definition Mdefines.h:73
TYPEOF
#define TYPEOF(s)
Definition Mdefines.h:123
symDN
void symDN(SEXP dest, SEXP src, int J)
Definition attrib.c:69
M2CHS
cholmod_sparse * M2CHS(SEXP obj, int values)
Definition cholmod-etc.c:132
c
cholmod_common c
Definition cholmod-etc.c:5
CHS2M
SEXP CHS2M(cholmod_sparse *A, int values, char shape)
Definition cholmod-etc.c:328
M2CHD
cholmod_dense * M2CHD(SEXP obj, char trans)
Definition cholmod-etc.c:176
dense_as_general
SEXP dense_as_general(SEXP from, const char *class, int new)
Definition coerce.c:2237
zvconj
void zvconj(Rcomplex *x, const Rcomplex *y, size_t n)
Definition idz.c:1300
ztrans2
void ztrans2(Rcomplex *, const Rcomplex *, size_t, size_t, char)
dtrans2
void dtrans2(double *, const double *, size_t, size_t, char)
Matrix_DimSym
SEXP Matrix_DimSym
Definition init.c:598
Matrix_zunit
Rcomplex Matrix_zunit
Definition init.c:642
Matrix_xSym
SEXP Matrix_xSym
Definition init.c:635
Matrix_iSym
SEXP Matrix_iSym
Definition init.c:607
Matrix_jSym
SEXP Matrix_jSym
Definition init.c:610
Matrix_zzero
Rcomplex Matrix_zzero
Definition init.c:641
Matrix_pSym
SEXP Matrix_pSym
Definition init.c:622
sparse_transpose
SEXP sparse_transpose(SEXP, const char *, char, int)
Definition t.c:52
Tsparse_rowscale
static void Tsparse_rowscale(SEXP obj, SEXP d, SEXP iSym)
Definition matmult.c:1294
R_dense_matmult
SEXP R_dense_matmult(SEXP s_x, SEXP s_y, SEXP s_xtrans, SEXP s_ytrans)
Definition matmult.c:711
matmultDN
static void matmultDN(SEXP dest, SEXP asrc, int ai, SEXP bsrc, int bi)
Definition matmult.c:131
gCgeMatrix_matmult
static SEXP gCgeMatrix_matmult(SEXP x, SEXP y, int xtrans, char ytrans, char ztrans, int triangular, int symmetric)
Definition matmult.c:939
DO_AS
#define DO_AS(_A_, _CLASS_, _TRANS_, _PID_)
dense_colscale
static void dense_colscale(SEXP obj, SEXP d, int m, int n, char ul, char nu)
Definition matmult.c:1174
Csparse_rowscale
static void Csparse_rowscale(SEXP obj, SEXP d, SEXP iSym)
Definition matmult.c:1269
sparse_dropzero
SEXP sparse_dropzero(SEXP, const char *, double)
Definition dropzero.c:4
ASMODE
#define ASMODE(_TRANS_)
tpMatrix_matmult
static SEXP tpMatrix_matmult(SEXP a, SEXP b, char atrans, char btrans, char aside, int triangular)
Definition matmult.c:613
syMatrix_matmult
static SEXP syMatrix_matmult(SEXP a, SEXP b, char atrans, char btrans, char aside)
Definition matmult.c:307
dense_transpose
SEXP dense_transpose(SEXP, const char *, char)
Definition t.c:7
Csparse_colscale
static void Csparse_colscale(SEXP obj, SEXP d)
Definition matmult.c:1238
dense_rowscale
static void dense_rowscale(SEXP obj, SEXP d, int m, int n, char ul, char nu)
Definition matmult.c:1224
R_sparse_matmult
SEXP R_sparse_matmult(SEXP s_x, SEXP s_y, SEXP s_xtrans, SEXP s_ytrans, SEXP s_ztrans, SEXP s_boolean)
Definition matmult.c:1021
DO_S3
#define DO_S3(_A_, _TRANS_, _PID_, _ISV_)
matmultDim
static void matmultDim(SEXP x, SEXP y, char *xtrans, char *ytrans, char *ztrans, int *m, int *n, int *v)
Definition matmult.c:20
CONJ2
#define CONJ2(_X_, _M_, _N_)
Definition matmult.c:152
valid_matmult
static const char * valid_matmult[]
Definition matmult.c:15
trMatrix_matmult
static SEXP trMatrix_matmult(SEXP a, SEXP b, char atrans, char btrans, char aside, int triangular)
Definition matmult.c:534
CONJ1
#define CONJ1(_X_, _N_)
Definition matmult.c:161
SCALE3
#define SCALE3(t, index)
Definition matmult.c:1163
DO_TR
#define DO_TR
gCgCMatrix_matmult
static SEXP gCgCMatrix_matmult(SEXP x, SEXP y, char xtrans, char ytrans, char ztrans, int triangular, int boolean)
Definition matmult.c:848
sparse_diag_U2N
SEXP sparse_diag_U2N(SEXP, const char *)
Definition diag.c:621
R_diagonal_matmult
SEXP R_diagonal_matmult(SEXP s_x, SEXP s_y, SEXP s_xtrans, SEXP s_ytrans, SEXP s_boolean)
Definition matmult.c:1310
geMatrix_matmult
static SEXP geMatrix_matmult(SEXP a, SEXP b, char atrans, char btrans)
Definition matmult.c:172
spMatrix_matmult
static SEXP spMatrix_matmult(SEXP a, SEXP b, char atrans, char btrans, char aside)
Definition matmult.c:428
Generated by doxygen 1.12.0