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