Matrix r5059
Loading...
Searching...
No Matches
expm.c
Go to the documentation of this file.
1/* C implementation of methods for expm */
2
3#include <complex.h> /* _Complex, cexp */
4#include "Lapack-etc.h"
5#include "Mdefines.h"
6#include "idz.h"
7
8#define SHOW(x, header) \
9do { \
10 Rprintf("%s ...\n", (header)); \
11 for (k = 0; k < nn; ++k) \
12 Rprintf(" % .8e\n", (x)[k]); \
13} while (0)
14
15/* defined in ./Schur.c : */
16SEXP dense_schur(SEXP, const char *, int, int);
17
18/* defined in ./coerce.c : */
19SEXP dense_as_kind(SEXP, const char *, char, int);
20SEXP dense_as_general(SEXP, const char *, int);
21SEXP dense_as_packed(SEXP, const char *, char, char, char);
22
23static double thetam[] = { 1.5e-2, 2.5e-1, 9.5e-1, 2.1e+0, 5.4e+0 };
24static double padecm[][14] = {
25 { 120.0, 60.0, 12.0, 1.0 },
26 { 30240.0, 15120.0, 3360.0, 420.0, 30.0, 1.0 },
27 { 17297280.0, 8648640.0, 1995840.0, 277200.0, 25200.0,
28 1512.0, 56.0, 1.0 },
29 { 17643225600.0, 8821612800.0, 2075673600.0, 302702400.0,
30 30270240.0, 2162160.0, 110880.0, 3960.0, 90.0, 1.0 },
31 { 64764752532480000.0, 32382376266240000.0,
32 7771770303897600.0, 1187353796428800.0, 129060195264000.0,
33 10559470521600.0, 670442572800.0, 33522128640.0,
34 1323241920.0, 40840800.0, 960960.0, 16380.0, 182.0, 1.0 } };
35
36SEXP dense_expm(SEXP obj, const char *class)
37{
38 int *pdim = DIM(obj), n = pdim[1];
39 if (pdim[0] != n)
40 Rf_error(_("matrix is not square"));
41 int cs = class[1] == 's' && class[0] == 'z' && TRANS(obj) != 'C';
42 char cl[] = "...Matrix";
43 cl[0] = (class[0] == 'z') ? 'z' : 'd';
44 cl[1] = (class[1] == 'g' || cs) ? 'g' : (class[1] == 's') ? 'p' : 't';
45 cl[2] = (class[1] == 'g' || cs) ? 'e' : (class[1] == 's') ? 'o' : 'r';
46
47 SEXP ans = PROTECT(newObject(cl));
48 SET_DIM(ans, n, n);
49 SET_DIMNAMES(ans, -(class[1] == 's'), DIMNAMES(obj, 0));
50 if (cl[1] != 'g' && UPLO(obj) != 'U')
51 SET_UPLO(ans);
52 if (n > 0) {
53
54 PROTECT_INDEX pid;
55 PROTECT_WITH_INDEX(obj, &pid);
56 if (class[0] != 'z' && class[0] != 'd') {
57 REPROTECT(obj = dense_as_kind(obj, class, ',', 1), pid);
58 class = Matrix_class(obj, valid_dense, 6, __func__);
59 }
60
61 if (class[1] != 's' || cs) {
62
63 /* Implementation of Algorithm 2.3 from Higham (2005). */
64
65 REPROTECT(obj = dense_as_general(obj, class, 1), pid);
66
67 SEXP x0 = PROTECT(GET_SLOT(obj, Matrix_xSym)),
68 x1 = PROTECT(Rf_allocVector(TYPEOF(x0), XLENGTH(x0)));
69 size_t n1 = (size_t) n, nn = n1 * n1, dk = n1 + 1, k;
70 int *pivot = (int *) R_alloc(n1, sizeof(int)),
71 i, j, l, s, ilo, ihi, info;
72 double *scale = (double *) R_alloc(n1, sizeof(double)),
73 norm, *b;
74
75 if (TYPEOF(x0) == CPLXSXP) {
76
77 Rcomplex *a = COMPLEX(x0), *x = COMPLEX(x1), *work, *tmp,
78 zero = Matrix_zzero, unit = Matrix_zunit, trace = zero;
79 memcpy(x, a, sizeof(Rcomplex) * nn);
80 /* 1. Shift */
81 k = 0;
82 for (j = 0; j < n; ++j) {
83 trace.r += x[k].r;
84 trace.i += x[k].i;
85 k += dk;
86 }
87 trace.r /= (double) n;
88 trace.i /= (double) n;
89 k = 0;
90 for (j = 0; j < n; ++j) {
91 x[k].r -= trace.r;
92 x[k].i -= trace.i;
93 k += dk;
94 }
95 /* 2. Balance */
96 F77_CALL(zgebal)("B", &n, x, &n, &ilo, &ihi, scale, &info FCONE);
97 ERROR_LAPACK_1(zgebal, info);
98 ilo -= 1; ihi -= 1;
99 /* 3. Scale */
100 norm = F77_CALL(zlange)("O", &n, &n, x, &n, (double *) 0 FCONE);
101 if (!R_FINITE(norm))
102 Rf_error("matrix one-norm is not finite");
103 s = (norm <= thetam[4]) ? 0 : (int) ceil(log2(norm / thetam[4]));
104 if (s > 0) {
105 double e = ldexp(1.0, -s);
106 for (k = 0; k < nn; ++k) {
107 x[k].r *= e;
108 x[k].i *= e;
109 }
110 }
111 /* 4. Pade approximation */
112 for (l = 0; l < 4; ++l)
113 if (norm <= thetam[l])
114 break;
115 b = &padecm[l][0]; /* Pade coefficients */
116 work = (Rcomplex *) R_alloc(nn * 2, sizeof(Rcomplex));
117 Rcomplex *u = work, *v = u + nn;
118 memset(u, 0, sizeof(Rcomplex) * nn * 2);
119 k = 0;
120 for (j = 0; j < n; ++j) {
121 u[k].r = b[1];
122 v[k].r = b[0];
123 k += dk;
124 }
125 if (l < 4) {
126 /* Pade approximant of degree 3, 5, 7, or 9 */
127 b += 2;
128 work = (Rcomplex *) R_alloc((l > 0) ? nn * 3 : nn, sizeof(Rcomplex));
129 Rcomplex *a2 = work, *a2m = work + nn, *a2m2 = a2m + nn;
130 F77_CALL(zgemm)("N", "N", &n, &n, &n, &unit, x, &n,
131 x, &n, &zero, a2, &n FCONE FCONE);
132 for (k = 0; k < nn; ++k) {
133 u[k].r += b[1] * a2[k].r;
134 u[k].i += b[1] * a2[k].i;
135 v[k].r += b[0] * a2[k].r;
136 v[k].i += b[0] * a2[k].i;
137 }
138 if (l > 0) {
139 b += 2;
140 memcpy(a2m, a2, sizeof(double) * nn);
141 for (i = 0; i < l; ++i) {
142 F77_CALL(zgemm)("N", "N", &n, &n, &n, &unit, a2, &n,
143 a2m, &n, &zero, a2m2, &n FCONE FCONE);
144 for (k = 0; k < nn; ++k) {
145 u[k].r += b[1] * a2m2[k].r;
146 u[k].i += b[1] * a2m2[k].i;
147 v[k].r += b[0] * a2m2[k].r;
148 v[k].i += b[0] * a2m2[k].i;
149 }
150 b += 2;
151 tmp = a2m2; a2m2 = a2m; a2m = tmp;
152 }
153 }
154 F77_CALL(zgemm)("N", "N", &n, &n, &n, &unit, x, &n,
155 u, &n, &zero, a2, &n FCONE FCONE);
156 memcpy(u, a2, sizeof(Rcomplex) * nn);
157 } else {
158 /* Pade approximant of degree 13 */
159 work = (Rcomplex *) R_alloc(nn * 5, sizeof(Rcomplex));
160 Rcomplex *a2 = work, *a4 = a2 + nn, *a6 = a4 + nn,
161 *upart = a6 + nn, *vpart = upart + nn;
162 F77_CALL(zgemm)("N", "N", &n, &n, &n, &unit, x , &n,
163 x , &n, &zero, a2, &n FCONE FCONE);
164 F77_CALL(zgemm)("N", "N", &n, &n, &n, &unit, a2, &n,
165 a2, &n, &zero, a4, &n FCONE FCONE);
166 F77_CALL(zgemm)("N", "N", &n, &n, &n, &unit, a2, &n,
167 a4, &n, &zero, a6, &n FCONE FCONE);
168 for (k = 0; k < nn; ++k) {
169 u[k].r += b[3] * a2[k].r + b[5] * a4[k].r + b[7] * a6[k].r;
170 u[k].i += b[3] * a2[k].i + b[5] * a4[k].i + b[7] * a6[k].i;
171 v[k].r += b[2] * a2[k].r + b[4] * a4[k].r + b[6] * a6[k].r;
172 v[k].i += b[2] * a2[k].i + b[4] * a4[k].i + b[6] * a6[k].i;
173 upart[k].r = b[9] * a2[k].r + b[11] * a4[k].r + b[13] * a6[k].r;
174 upart[k].i = b[9] * a2[k].i + b[11] * a4[k].i + b[13] * a6[k].i;
175 vpart[k].r = b[8] * a2[k].r + b[10] * a4[k].r + b[12] * a6[k].r;
176 vpart[k].i = b[8] * a2[k].i + b[10] * a4[k].i + b[12] * a6[k].i;
177 }
178 F77_CALL(zgemm)("N", "N", &n, &n, &n, &unit, a6, &n,
179 upart, &n, &unit, u, &n FCONE FCONE);
180 F77_CALL(zgemm)("N", "N", &n, &n, &n, &unit, a6, &n,
181 vpart, &n, &unit, v, &n FCONE FCONE);
182 F77_CALL(zgemm)("N", "N", &n, &n, &n, &unit, x, &n,
183 u, &n, &zero, a2, &n FCONE FCONE);
184 memcpy(u, a2, sizeof(Rcomplex) * nn);
185 }
186 for (k = 0; k < nn; ++k) {
187 x[k].r = u[k].r + v[k].r;
188 x[k].i = u[k].i + v[k].i;
189 u[k].r = -u[k].r + v[k].r;
190 u[k].i = -u[k].i + v[k].i;
191 }
192 F77_CALL(zgetrf)(&n, &n, u, &n, pivot, &info);
193 ERROR_LAPACK_2(zgetrf, info, 2, U);
194 F77_CALL(zgetrs)("N", &n, &n, u, &n, pivot, x, &n, &info FCONE);
195 ERROR_LAPACK_1(zgetrs, info);
196 /* 5. Square */
197 if (s > 0) {
198 u = x;
199 for (i = 0; i < s; ++i) {
200 F77_CALL(zgemm)("N", "N", &n, &n, &n, &unit, u, &n,
201 u, &n, &zero, v, &n FCONE FCONE);
202 tmp = v; v = u; u = tmp;
203 }
204 if (s % 2)
205 memcpy(x, u, sizeof(Rcomplex) * nn);
206 }
207 /* 6. Inverse balance */
208 tmp = x;
209 for (j = 0; j < n; ++j) {
210 for (i = ilo; i <= ihi; ++i) {
211 tmp[i].r *= scale[i];
212 tmp[i].i *= scale[i];
213 }
214 tmp += n;
215 }
216 tmp = x + n1 * (size_t) ilo;
217 for (j = ilo; j <= ihi; ++j) {
218 for (i = 0; i < n; ++i) {
219 tmp[i].r /= scale[j];
220 tmp[i].i /= scale[j];
221 }
222 tmp += n;
223 }
224 for (j = ilo - 1; j >= 0; --j) {
225 i = (int) scale[j] - 1;
226 if (i != j) {
227 zswap2(n1, x + n1 * (size_t) i, 1, x + n1 * (size_t) j, 1);
228 zswap2(n1, x + i, n1, x + j, n1);
229 }
230 }
231 for (j = ihi + 1; j < n; ++j) {
232 i = (int) scale[j] - 1;
233 if (i != j) {
234 zswap2(n1, x + n1 * (size_t) i, 1, x + n1 * (size_t) j, 1);
235 zswap2(n1, x + i, n1, x + j, n1);
236 }
237 }
238 /* 7. Inverse shift */
239#define asC(x) ((double _Complex *) &(x))[0]
240#define asR(x) (( Rcomplex *) &(x))[0]
241 Rcomplex x___;
242 double _Complex z = cexp(asC(trace));
243 trace = asR(z);
244 for (k = 0; k < nn; ++k) {
245 x___ = x[k];
246 x[k].r = x___.r * trace.r - x___.i * trace.i;
247 x[k].i = x___.r * trace.i + x___.i * trace.r;
248 }
249
250 } else {
251
252 double *a = REAL(x0), *x = REAL(x1), *work, *tmp,
253 zero = 0.0, unit = 1.0, trace = zero;
254 memcpy(x, a, sizeof(double) * nn);
255 /* 1. Shift */
256 k = 0;
257 for (j = 0; j < n; ++j) {
258 trace += x[k];
259 k += dk;
260 }
261 trace /= (double) n;
262 k = 0;
263 for (j = 0; j < n; ++j) {
264 x[k] -= trace;
265 k += dk;
266 }
267 /* 2. Balance */
268 F77_CALL(dgebal)("B", &n, x, &n, &ilo, &ihi, scale, &info FCONE);
269 ERROR_LAPACK_1(dgebal, info);
270 ilo -= 1; ihi -= 1;
271 /* 3. Scale */
272 norm = F77_CALL(dlange)("O", &n, &n, x, &n, (double *) 0 FCONE);
273 if (!R_FINITE(norm))
274 Rf_error("matrix one-norm is not finite");
275 s = (norm <= thetam[4]) ? 0 : (int) ceil(log2(norm / thetam[4]));
276 if (s > 0) {
277 double e = ldexp(1.0, -s);
278 for (k = 0; k < nn; ++k)
279 x[k] *= e;
280 }
281 /* 4. Pade approximation */
282 for (l = 0; l < 4; ++l)
283 if (norm <= thetam[l])
284 break;
285 b = &padecm[l][0]; /* Pade coefficients */
286 work = (double *) R_alloc(nn * 2, sizeof(double));
287 double *u = work, *v = u + nn;
288 memset(u, 0, sizeof(double) * nn * 2);
289 k = 0;
290 for (j = 0; j < n; ++j) {
291 u[k] = b[1];
292 v[k] = b[0];
293 k += dk;
294 }
295 if (l < 4) {
296 /* Pade approximant of degree 3, 5, 7, or 9 */
297 b += 2;
298 work = (double *) R_alloc((l > 0) ? nn * 3 : nn, sizeof(double));
299 double *a2 = work, *a2m = work + nn, *a2m2 = a2m + nn;
300 F77_CALL(dgemm)("N", "N", &n, &n, &n, &unit, x, &n,
301 x, &n, &zero, a2, &n FCONE FCONE);
302 for (k = 0; k < nn; ++k) {
303 u[k] += b[1] * a2[k];
304 v[k] += b[0] * a2[k];
305 }
306 if (l > 0) {
307 b += 2;
308 memcpy(a2m, a2, sizeof(double) * nn);
309 for (i = 0; i < l; ++i) {
310 F77_CALL(dgemm)("N", "N", &n, &n, &n, &unit, a2, &n,
311 a2m, &n, &zero, a2m2, &n FCONE FCONE);
312 for (k = 0; k < nn; ++k) {
313 u[k] += b[1] * a2m2[k];
314 v[k] += b[0] * a2m2[k];
315 }
316 b += 2;
317 tmp = a2m2; a2m2 = a2m; a2m = tmp;
318 }
319 }
320 F77_CALL(dgemm)("N", "N", &n, &n, &n, &unit, x, &n,
321 u, &n, &zero, a2, &n FCONE FCONE);
322 memcpy(u, a2, sizeof(double) * nn);
323 } else {
324 /* Pade approximant of degree 13 */
325 work = (double *) R_alloc(nn * 5, sizeof(double));
326 double *a2 = work, *a4 = a2 + nn, *a6 = a4 + nn,
327 *upart = a6 + nn, *vpart = upart + nn;
328 F77_CALL(dgemm)("N", "N", &n, &n, &n, &unit, x , &n,
329 x , &n, &zero, a2, &n FCONE FCONE);
330 F77_CALL(dgemm)("N", "N", &n, &n, &n, &unit, a2, &n,
331 a2, &n, &zero, a4, &n FCONE FCONE);
332 F77_CALL(dgemm)("N", "N", &n, &n, &n, &unit, a2, &n,
333 a4, &n, &zero, a6, &n FCONE FCONE);
334 for (k = 0; k < nn; ++k) {
335 u[k] += b[3] * a2[k] + b[5] * a4[k] + b[7] * a6[k];
336 v[k] += b[2] * a2[k] + b[4] * a4[k] + b[6] * a6[k];
337 upart[k] = b[9] * a2[k] + b[11] * a4[k] + b[13] * a6[k];
338 vpart[k] = b[8] * a2[k] + b[10] * a4[k] + b[12] * a6[k];
339 }
340 F77_CALL(dgemm)("N", "N", &n, &n, &n, &unit, a6, &n,
341 upart, &n, &unit, u, &n FCONE FCONE);
342 F77_CALL(dgemm)("N", "N", &n, &n, &n, &unit, a6, &n,
343 vpart, &n, &unit, v, &n FCONE FCONE);
344 F77_CALL(dgemm)("N", "N", &n, &n, &n, &unit, x, &n,
345 u, &n, &zero, a2, &n FCONE FCONE);
346 memcpy(u, a2, sizeof(double) * nn);
347 }
348 for (k = 0; k < nn; ++k) {
349 x[k] = u[k] + v[k];
350 u[k] = -u[k] + v[k];
351 }
352 F77_CALL(dgetrf)(&n, &n, u, &n, pivot, &info);
353 ERROR_LAPACK_2(dgetrf, info, 2, U);
354 F77_CALL(dgetrs)("N", &n, &n, u, &n, pivot, x, &n, &info FCONE);
355 ERROR_LAPACK_1(dgetrs, info);
356 /* 5. Square */
357 if (s > 0) {
358 u = x;
359 for (i = 0; i < s; ++i) {
360 F77_CALL(dgemm)("N", "N", &n, &n, &n, &unit, u, &n,
361 u, &n, &zero, v, &n FCONE FCONE);
362 tmp = v; v = u; u = tmp;
363 }
364 if (s % 2)
365 memcpy(x, u, sizeof(double) * nn);
366 }
367 /* 6. Inverse balance */
368 tmp = x;
369 for (j = 0; j < n; ++j) {
370 for (i = ilo; i <= ihi; ++i)
371 tmp[i] *= scale[i];
372 tmp += n;
373 }
374 tmp = x + n1 * (size_t) ilo;
375 for (j = ilo; j <= ihi; ++j) {
376 for (i = 0; i < n; ++i)
377 tmp[i] /= scale[j];
378 tmp += n;
379 }
380 for (j = ilo - 1; j >= 0; --j) {
381 i = (int) scale[j] - 1;
382 if (i != j) {
383 dswap2(n1, x + n1 * (size_t) i, 1, x + n1 * (size_t) j, 1);
384 dswap2(n1, x + i, n1, x + j, n1);
385 }
386 }
387 for (j = ihi + 1; j < n; ++j) {
388 i = (int) scale[j] - 1;
389 if (i != j) {
390 dswap2(n1, x + n1 * (size_t) i, 1, x + n1 * (size_t) j, 1);
391 dswap2(n1, x + i, n1, x + j, n1);
392 }
393 }
394 /* 7. Inverse shift */
395 trace = exp(trace);
396 for (k = 0; k < nn; ++k)
397 x[k] *= trace;
398
399 }
400
401 SET_SLOT(ans, Matrix_xSym, x1);
402 UNPROTECT(2); /* x1, x0 */
403
404 } else {
405
406 /* Use the Schur factorization as the argument is Hermitian. */
407
408 REPROTECT(obj = dense_schur(obj, class, 2, 1), pid);
409 SEXP vectors = PROTECT(GET_SLOT(obj, Matrix_vectorsSym)),
410 values = PROTECT(GET_SLOT(obj, Matrix_valuesSym)),
411 x1 = PROTECT(Rf_allocVector(TYPEOF(vectors), XLENGTH(vectors)));
412 double *w = REAL(values);
413 size_t n1 = (size_t) n, nn = n1 * n1;
414 int i, j;
415 for (j = 0; j < n; ++j)
416 w[j] = exp(w[j]);
417 if (TYPEOF(vectors) == CPLXSXP) {
418 Rcomplex *qw = (Rcomplex *) R_alloc(nn, sizeof(Rcomplex)),
419 *q = COMPLEX(vectors), *x = COMPLEX(x1), *tmp = qw,
420 zero = Matrix_zzero, unit = Matrix_zunit;
421 for (j = 0; j < n; ++j) {
422 for (i = 0; i < n; ++i) {
423 qw[i].r = q[i].r * w[j];
424 qw[i].i = q[i].i * w[j];
425 }
426 qw += n;
427 q += n;
428 }
429 qw = tmp; q = COMPLEX(vectors);
430 F77_CALL(zgemm)("N", "C", &n, &n, &n, &unit, qw, &n,
431 q, &n, &zero, x, &n FCONE FCONE);
432 } else {
433 double *qw = (double *) R_alloc(nn, sizeof(double)),
434 *q = REAL(vectors), *x = REAL(x1), *tmp = qw,
435 zero = 0.0, unit = 1.0;
436 for (j = 0; j < n; ++j) {
437 for (i = 0; i < n; ++i)
438 qw[i] = q[i] * w[j];
439 qw += n;
440 q += n;
441 }
442 qw = tmp; q = REAL(vectors);
443 F77_CALL(dgemm)("N", "C", &n, &n, &n, &unit, qw, &n,
444 q, &n, &zero, x, &n FCONE FCONE);
445 }
446
447 SET_SLOT(ans, Matrix_xSym, x1);
448 UNPROTECT(3); /* x1, values, vectors */
449
450 }
451
452 UNPROTECT(1); /* obj */
453
454 }
455
456 if (cl[1] != 'g' && class[2] == 'p') {
457 if (class[1] == 's')
458 ans = dense_as_packed(ans, cl, '\0', 'C', '\0');
459 else
460 ans = dense_as_packed(ans, cl, '\0', '\0', 'N');
461 }
462 UNPROTECT(1); /* ans */
463 return ans;
464}
465
466SEXP R_dense_expm(SEXP s_obj)
467{
468 const char *class = Matrix_class(s_obj, valid_dense, 6, __func__);
469 return dense_expm(s_obj, class);
470}
FCONE
#define FCONE
Definition Lapack-etc.h:13
ERROR_LAPACK_1
#define ERROR_LAPACK_1(_ROUTINE_, _INFO_)
Definition Lapack-etc.h:18
ERROR_LAPACK_2
#define ERROR_LAPACK_2(_ROUTINE_, _INFO_, _WARN_, _LETTER_)
Definition Lapack-etc.h:25
SET_DIM
#define SET_DIM(x, m, n)
Definition Mdefines.h:87
valid_dense
const char * valid_dense[]
Definition objects.c:3
_
#define _(String)
Definition Mdefines.h:66
SET_UPLO
#define SET_UPLO(x)
Definition Mdefines.h:103
UPLO
#define UPLO(x)
Definition Mdefines.h:101
SET_DIMNAMES
#define SET_DIMNAMES(x, mode, value)
Definition Mdefines.h:98
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
TRANS
#define TRANS(x)
Definition Mdefines.h:106
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
cl
cholmod_common cl
Definition cholmod-etc.c:6
dense_expm
SEXP dense_expm(SEXP obj, const char *class)
Definition expm.c:36
asR
#define asR(x)
thetam
static double thetam[]
Definition expm.c:23
dense_as_general
SEXP dense_as_general(SEXP, const char *, int)
Definition coerce.c:2237
dense_as_packed
SEXP dense_as_packed(SEXP, const char *, char, char, char)
Definition coerce.c:2654
R_dense_expm
SEXP R_dense_expm(SEXP s_obj)
Definition expm.c:466
dense_schur
SEXP dense_schur(SEXP, const char *, int, int)
Definition Schur.c:8
asC
#define asC(x)
padecm
static double padecm[][14]
Definition expm.c:24
dense_as_kind
SEXP dense_as_kind(SEXP, const char *, char, int)
Definition coerce.c:2000
dswap2
void dswap2(size_t, double *, size_t, double *, size_t)
zswap2
void zswap2(size_t, Rcomplex *, size_t, Rcomplex *, size_t)
Matrix_valuesSym
SEXP Matrix_valuesSym
Definition init.c:633
Matrix_zunit
Rcomplex Matrix_zunit
Definition init.c:642
Matrix_xSym
SEXP Matrix_xSym
Definition init.c:635
Matrix_zzero
Rcomplex Matrix_zzero
Definition init.c:641
Matrix_vectorsSym
SEXP Matrix_vectorsSym
Definition init.c:634
Generated by doxygen 1.12.0