SCM

SCM Repository

[matrix] Annotation of /pkg/src/dgeMatrix.c
ViewVC logotype

Annotation of /pkg/src/dgeMatrix.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 655 - (view) (download) (as text)

1 : bates 478 #include "dgeMatrix.h"
2 : bates 10
3 : bates 478 SEXP dgeMatrix_validate(SEXP obj)
4 : bates 10 {
5 :     SEXP x = GET_SLOT(obj, Matrix_xSym),
6 : bates 332 Dim = GET_SLOT(obj, Matrix_DimSym),
7 : bates 476 fact = GET_SLOT(obj, Matrix_factorSym),
8 : maechler 534 rc = GET_SLOT(obj, Matrix_rcondSym);
9 : bates 10 int m, n;
10 :    
11 :     if (length(Dim) != 2)
12 : bates 582 return mkString(_("Dim slot must have length 2"));
13 : bates 10 m = INTEGER(Dim)[0]; n = INTEGER(Dim)[1];
14 :     if (m < 0 || n < 0)
15 : bates 582 return mkString(_("Negative value(s) in Dim"));
16 : bates 10 if (length(x) != m * n)
17 : maechler 604 return mkString(_("length of x slot != prod(Dim)"));
18 :     if (!isReal(x))
19 :     return mkString(_("x slot must be numeric \"double\""));
20 : bates 332 if (length(fact) > 0 && getAttrib(fact, R_NamesSymbol) == R_NilValue)
21 : bates 582 return mkString(_("factors slot must be named list"));
22 : bates 332 if (length(rc) > 0 && getAttrib(rc, R_NamesSymbol) == R_NilValue)
23 : bates 582 return mkString(_("rcond slot must be named numeric vector"));
24 : bates 10 return ScalarLogical(1);
25 :     }
26 :    
27 :     static
28 :     double get_norm(SEXP obj, char *typstr)
29 : maechler 534 {
30 : bates 10 char typnm[] = {'\0', '\0'};
31 :     int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym));
32 :     double *work = (double *) NULL;
33 :    
34 :     typnm[0] = norm_type(typstr);
35 :     if (*typnm == 'I') {
36 :     work = (double *) R_alloc(dims[0], sizeof(double));
37 :     }
38 :     return F77_CALL(dlange)(typstr, dims, dims+1,
39 :     REAL(GET_SLOT(obj, Matrix_xSym)),
40 :     dims, work);
41 :     }
42 :    
43 : bates 478 SEXP dgeMatrix_norm(SEXP obj, SEXP type)
44 : bates 10 {
45 :     return ScalarReal(get_norm(obj, CHAR(asChar(type))));
46 :     }
47 :    
48 :     static
49 :     double set_rcond(SEXP obj, char *typstr)
50 :     {
51 :     char typnm[] = {'\0', '\0'};
52 : bates 332 SEXP rcv = GET_SLOT(obj, Matrix_rcondSym);
53 : bates 10 double rcond = get_double_by_name(rcv, typnm);
54 :    
55 :     typnm[0] = rcond_type(typstr);
56 :     if (R_IsNA(rcond)) {
57 : bates 478 SEXP LU = dgeMatrix_LU(obj);
58 : bates 10 int *dims = INTEGER(GET_SLOT(LU, Matrix_DimSym)), info;
59 :     double anorm = get_norm(obj, typstr);
60 :    
61 :     if (dims[0] != dims[1] || dims[0] < 1)
62 : bates 582 error(_("rcond requires a square, non-empty matrix"));
63 : bates 10 F77_CALL(dgecon)(typnm,
64 :     dims, REAL(GET_SLOT(LU, Matrix_xSym)),
65 :     dims, &anorm, &rcond,
66 :     (double *) R_alloc(4*dims[0], sizeof(double)),
67 :     (int *) R_alloc(dims[0], sizeof(int)), &info);
68 : bates 332 SET_SLOT(obj, Matrix_rcondSym,
69 : bates 10 set_double_by_name(rcv, rcond, typnm));
70 :     }
71 :     return rcond;
72 :     }
73 :    
74 : bates 478 SEXP dgeMatrix_rcond(SEXP obj, SEXP type)
75 : bates 10 {
76 :     return ScalarReal(set_rcond(obj, CHAR(asChar(type))));
77 :     }
78 :    
79 : bates 478 SEXP dgeMatrix_crossprod(SEXP x)
80 : bates 10 {
81 : bates 478 SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dpoMatrix")));
82 : bates 10 int *Dims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
83 : bates 342 *vDims;
84 : bates 343 int i, n = Dims[1];
85 :     int nsqr = n * n;
86 :     double one = 1.0, *xvals, zero = 0.0;
87 : bates 10
88 : bates 476 SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));
89 : bates 332 SET_SLOT(val, Matrix_rcondSym, allocVector(REALSXP, 0));
90 : maechler 534 SET_SLOT(val, Matrix_uploSym, mkString("U"));
91 : bates 342 SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));
92 :     vDims = INTEGER(GET_SLOT(val, Matrix_DimSym));
93 : bates 332 vDims[0] = vDims[1] = n;
94 : bates 345 SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, nsqr));
95 : bates 343 xvals = REAL(GET_SLOT(val, Matrix_xSym));
96 :     for (i = 0; i < nsqr; i++) xvals[i] = 0.; /* keep valgrind happy */
97 : bates 345 F77_CALL(dsyrk)("U", "T", vDims, Dims,
98 :     &one, REAL(GET_SLOT(x, Matrix_xSym)), Dims,
99 :     &zero, xvals, vDims);
100 : bates 10 UNPROTECT(1);
101 :     return val;
102 :     }
103 :    
104 : bates 478 SEXP dgeMatrix_dgeMatrix_crossprod(SEXP x, SEXP y)
105 : bates 10 {
106 : bates 478 SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
107 : bates 10 int *xDims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
108 :     *yDims = INTEGER(GET_SLOT(y, Matrix_DimSym)),
109 : bates 332 *vDims;
110 : bates 10 int m = xDims[1], n = yDims[1];
111 :     double one = 1.0, zero = 0.0;
112 :    
113 : bates 332 SET_SLOT(val, Matrix_rcondSym, allocVector(REALSXP, 0));
114 : bates 476 SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));
115 : bates 332 SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));
116 :     vDims = INTEGER(GET_SLOT(val, Matrix_DimSym));
117 : bates 10 if ((*xDims) > 0 && (*yDims) > 0 && n > 0 && m > 0) {
118 :     if (*xDims != *yDims)
119 : bates 582 error(_("Dimensions of x and y are not compatible for crossprod"));
120 : bates 10 vDims[0] = m; vDims[1] = n;
121 :     SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, m * n));
122 :     F77_CALL(dgemm)("T", "N", xDims + 1, yDims + 1, xDims, &one,
123 :     REAL(GET_SLOT(x, Matrix_xSym)), xDims,
124 :     REAL(GET_SLOT(y, Matrix_xSym)), yDims,
125 :     &zero, REAL(GET_SLOT(val, Matrix_xSym)),
126 :     xDims + 1);
127 :     }
128 :     UNPROTECT(1);
129 :     return val;
130 :     }
131 :    
132 : bates 478 SEXP dgeMatrix_matrix_crossprod(SEXP x, SEXP y)
133 : bates 10 {
134 : bates 478 SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
135 : bates 10 int *xDims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
136 :     *yDims = INTEGER(getAttrib(y, R_DimSymbol)),
137 : bates 332 *vDims;
138 : bates 10 int m = xDims[1], n = yDims[1];
139 :     double one = 1.0, zero = 0.0;
140 :    
141 :     if (!(isMatrix(y) && isReal(y)))
142 : bates 582 error(_("Argument y must be a numeric (real) matrix"));
143 : bates 332 SET_SLOT(val, Matrix_rcondSym, allocVector(REALSXP, 0));
144 : bates 476 SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));
145 : bates 332 SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));
146 :     vDims = INTEGER(GET_SLOT(val, Matrix_DimSym));
147 : bates 10 if ((*xDims) > 0 && (*yDims) > 0 && n > 0 && m > 0) {
148 :     if (*xDims != *yDims)
149 : bates 582 error(_("Dimensions of x and y are not compatible for crossprod"));
150 : bates 10 vDims[0] = m; vDims[1] = n;
151 :     SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, m * n));
152 :     F77_CALL(dgemm)("T", "N", xDims + 1, yDims + 1, xDims, &one,
153 :     REAL(GET_SLOT(x, Matrix_xSym)), xDims,
154 :     REAL(y), yDims,
155 :     &zero, REAL(GET_SLOT(val, Matrix_xSym)),
156 :     xDims + 1);
157 :     }
158 :     UNPROTECT(1);
159 :     return val;
160 :     }
161 :    
162 : bates 478 SEXP dgeMatrix_getDiag(SEXP x)
163 : bates 10 {
164 :     int *dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
165 :     int i, m = dims[0], nret = (m < dims[1]) ? m : dims[1];
166 :     SEXP ret = PROTECT(allocVector(REALSXP, nret)),
167 :     xv = GET_SLOT(x, Matrix_xSym);
168 :    
169 :     for (i = 0; i < nret; i++) {
170 :     REAL(ret)[i] = REAL(xv)[i * (m + 1)];
171 :     }
172 :     UNPROTECT(1);
173 :     return ret;
174 :     }
175 :    
176 : bates 478 SEXP dgeMatrix_LU(SEXP x)
177 : bates 10 {
178 : bates 476 SEXP val = get_factors(x, "LU");
179 : bates 10 int *dims, npiv, info;
180 : maechler 534
181 : bates 10 if (val != R_NilValue) return val;
182 :     dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
183 :     if (dims[0] < 1 || dims[1] < 1)
184 : bates 582 error(_("Cannot factor a matrix with zero extents"));
185 : bates 10 npiv = (dims[0] <dims[1]) ? dims[0] : dims[1];
186 :     val = PROTECT(NEW_OBJECT(MAKE_CLASS("LU")));
187 :     SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(x, Matrix_xSym)));
188 :     SET_SLOT(val, Matrix_DimSym, duplicate(GET_SLOT(x, Matrix_DimSym)));
189 :     F77_CALL(dgetrf)(dims, dims + 1, REAL(GET_SLOT(val, Matrix_xSym)),
190 : bates 652 dims,
191 :     INTEGER(ALLOC_SLOT(val, Matrix_permSym, INTSXP, npiv)),
192 : bates 10 &info);
193 : bates 652 if (info)
194 :     error(_("Lapack routine %s returned error code %d"), "dgetrf", info);
195 : bates 10 UNPROTECT(1);
196 : bates 476 return set_factors(x, val, "LU");
197 : bates 10 }
198 :    
199 : bates 478 SEXP dgeMatrix_determinant(SEXP x, SEXP logarithm)
200 : bates 10 {
201 :     int lg = asLogical(logarithm);
202 : bates 478 SEXP lu = dgeMatrix_LU(x);
203 : bates 10 int *dims = INTEGER(GET_SLOT(lu, Matrix_DimSym)),
204 : bates 652 *jpvt = INTEGER(GET_SLOT(lu, Matrix_permSym)),
205 : bates 10 i, n = dims[0], sign = 1;
206 :     double *luvals = REAL(GET_SLOT(lu, Matrix_xSym)), modulus;
207 :    
208 :     if (n != dims[1])
209 : bates 582 error(_("Determinant requires a square matrix"));
210 : bates 10 for (i = 0; i < n; i++) if (jpvt[i] != (i + 1)) sign = -sign;
211 :     if (lg) {
212 :     modulus = 0.0;
213 :     for (i = 0; i < n; i++) {
214 :     double dii = luvals[i*(n + 1)]; /* ith diagonal element */
215 :     modulus += log(dii < 0 ? -dii : dii);
216 :     if (dii < 0) sign = -sign;
217 :     }
218 :     } else {
219 :     modulus = 1.0;
220 :     for (i = 0; i < n; i++)
221 :     modulus *= luvals[i*(n + 1)];
222 :     if (modulus < 0) {
223 :     modulus = -modulus;
224 :     sign = -sign;
225 :     }
226 :     }
227 :     return as_det_obj(modulus, lg, sign);
228 :     }
229 :    
230 : bates 478 SEXP dgeMatrix_solve(SEXP a)
231 : bates 10 {
232 : bates 478 SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix"))),
233 :     lu = dgeMatrix_LU(a);
234 : bates 10 int *dims = INTEGER(GET_SLOT(lu, Matrix_DimSym)),
235 : bates 652 *pivot = INTEGER(GET_SLOT(lu, Matrix_permSym));
236 : bates 10 double *x, tmp;
237 :     int info, lwork = -1;
238 :    
239 :    
240 : bates 582 if (dims[0] != dims[1]) error(_("Solve requires a square matrix"));
241 : bates 10 SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(lu, Matrix_xSym)));
242 :     x = REAL(GET_SLOT(val, Matrix_xSym));
243 :     SET_SLOT(val, Matrix_DimSym, duplicate(GET_SLOT(lu, Matrix_DimSym)));
244 :     F77_CALL(dgetri)(dims, x, dims, pivot, &tmp, &lwork, &info);
245 :     lwork = (int) tmp;
246 :     F77_CALL(dgetri)(dims, x, dims, pivot,
247 :     (double *) R_alloc((size_t) lwork, sizeof(double)),
248 :     &lwork, &info);
249 :     UNPROTECT(1);
250 :     return val;
251 :     }
252 :    
253 : bates 652 SEXP dgeMatrix_matrix_solve(SEXP a, SEXP b, SEXP classed)
254 : bates 10 {
255 : bates 652 int cl = asLogical(classed);
256 :     SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix"))),
257 :     lu = dgeMatrix_LU(a);
258 :     int *adims = INTEGER(GET_SLOT(lu, Matrix_DimSym)),
259 :     *bdims = INTEGER(cl ? GET_SLOT(b, Matrix_DimSym) :
260 :     getAttrib(b, R_DimSymbol));
261 :     int info, n = bdims[0], nrhs = bdims[1];
262 :     int sz = n * nrhs;
263 :    
264 :     if (*adims != *bdims || bdims[1] < 1 || *adims < 1 || *adims != adims[1])
265 :     error(_("Dimensions of system to be solved are inconsistent"));
266 :     Memcpy(INTEGER(ALLOC_SLOT(val, Matrix_DimSym, INTSXP, 2)), bdims, 2);
267 :     F77_CALL(dgetrs)("N", &n, &nrhs, REAL(GET_SLOT(val, Matrix_xSym)), &n,
268 : maechler 655 INTEGER(GET_SLOT(lu, Matrix_permSym)),
269 : bates 652 Memcpy(REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, sz)),
270 :     REAL(cl ? GET_SLOT(b, Matrix_xSym):b), sz),
271 :     &n, &info);
272 :     UNPROTECT(1);
273 :     return val;
274 :     }
275 :    
276 :     SEXP dgeMatrix_matrix_mm(SEXP a, SEXP b, SEXP classed, SEXP right)
277 :     {
278 : maechler 655 int cl = asLogical(classed);
279 : bates 652 SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
280 : bates 10 int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
281 : bates 652 *bdims = INTEGER(cl ? GET_SLOT(b, Matrix_DimSym) :
282 :     getAttrib(b, R_DimSymbol)),
283 :     *cdims = INTEGER(ALLOC_SLOT(val, Matrix_DimSym, INTSXP, 2));
284 : bates 10 double one = 1., zero = 0.;
285 : maechler 534
286 : bates 652 if (asLogical(right)) {
287 :     int m = bdims[0], n = adims[1], k = bdims[1];
288 :     if (adims[0] != k)
289 :     error(_("Matrices are not conformable for multiplication"));
290 :     if (m < 1 || n < 1 || k < 1)
291 :     error(_("Matrices with zero extents cannot be multiplied"));
292 :     cdims[0] = m; cdims[1] = n;
293 :     F77_CALL(dgemm) ("N", "N", &m, &n, &k, &one,
294 :     REAL(cl ? GET_SLOT(b, Matrix_xSym) : b), &m,
295 :     REAL(GET_SLOT(a, Matrix_xSym)), &k, &zero,
296 :     REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, m * n)),
297 :     &m);
298 :     } else {
299 :     int m = adims[0], n = bdims[1], k = adims[1];
300 :    
301 :     if (bdims[0] != k)
302 :     error(_("Matrices are not conformable for multiplication"));
303 :     if (m < 1 || n < 1 || k < 1)
304 :     error(_("Matrices with zero extents cannot be multiplied"));
305 :     cdims[0] = m; cdims[1] = n;
306 :     F77_CALL(dgemm)
307 :     ("N", "N", &m, &n, &k, &one, REAL(GET_SLOT(a, Matrix_xSym)),
308 :     &m, REAL(cl ? GET_SLOT(b, Matrix_xSym) : b), &k, &zero,
309 :     REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, m * n)), &m);
310 :     }
311 : bates 10 UNPROTECT(1);
312 :     return val;
313 :     }
314 : bates 301
315 : bates 652
316 : bates 478 SEXP dgeMatrix_svd(SEXP x, SEXP nnu, SEXP nnv)
317 : bates 301 {
318 : bates 377 int /* nu = asInteger(nnu),
319 :     nv = asInteger(nnv), */
320 : bates 301 *dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
321 :     double *xx = REAL(GET_SLOT(x, Matrix_xSym));
322 :     SEXP val = PROTECT(allocVector(VECSXP, 3));
323 :    
324 :     if (dims[0] && dims[1]) {
325 :     int m = dims[0], n = dims[1], mm = (m < n)?m:n,
326 :     lwork = -1, info;
327 :     int *iwork = Calloc(8 * mm, int);
328 :     double tmp, *work;
329 : bates 306 /* int bdspac = 3*m*m + 4*m, */
330 :     /* wrkbl, maxwrk, minwrk, itmp, */
331 :     /* ione = 1, iminus1 = -1; */
332 :     /* int i1, i2, i3; */
333 : bates 301
334 :     SET_VECTOR_ELT(val, 0, allocVector(REALSXP, mm));
335 :     SET_VECTOR_ELT(val, 1, allocMatrix(REALSXP, m, mm));
336 :     SET_VECTOR_ELT(val, 2, allocMatrix(REALSXP, mm, n));
337 :     F77_CALL(dgesdd)("S", &m, &n, xx, &m,
338 :     REAL(VECTOR_ELT(val, 0)),
339 :     REAL(VECTOR_ELT(val, 1)), &m,
340 :     REAL(VECTOR_ELT(val, 2)), &mm,
341 :     &tmp, &lwork, iwork, &info);
342 :     lwork = (int) tmp;
343 : bates 306 /* F77_CALL(foo)(&i1, &i2, &i3); */
344 :     /* wrkbl = 3*m+(m+n)*i1; */
345 :     /* if (wrkbl < (itmp = 3*m + m*i2)) wrkbl = itmp; */
346 :     /* if (wrkbl < (itmp = 3*m + m*i3)) wrkbl = itmp; */
347 :     /* itmp = bdspac+3*m; */
348 :     /* maxwrk = (wrkbl > itmp) ? wrkbl : itmp; */
349 :     /* minwrk = 3*m + ((bdspac > n) ? bdspac : n); */
350 : bates 301 work = Calloc(lwork, double);
351 :     F77_CALL(dgesdd)("S", &m, &n, xx, &m,
352 :     REAL(VECTOR_ELT(val, 0)),
353 :     REAL(VECTOR_ELT(val, 1)), &m,
354 :     REAL(VECTOR_ELT(val, 2)), &mm,
355 :     work, &lwork, iwork, &info);
356 :     Free(iwork); Free(work);
357 :     }
358 :     UNPROTECT(1);
359 :     return val;
360 :     }
361 : bates 461
362 :     static double padec [] = /* Constants for matrix exponential calculation. */
363 :     {
364 :     5.0000000000000000e-1,
365 :     1.1666666666666667e-1,
366 :     1.6666666666666667e-2,
367 :     1.6025641025641026e-3,
368 :     1.0683760683760684e-4,
369 :     4.8562548562548563e-6,
370 :     1.3875013875013875e-7,
371 :     1.9270852604185938e-9,
372 :     };
373 :    
374 : maechler 534 /**
375 : bates 461 * Matrix exponential - based on the code for Octave's expm function.
376 : maechler 534 *
377 : bates 461 * @param x real square matrix to exponentiate
378 : maechler 534 *
379 : bates 461 * @return matrix exponential of x
380 :     */
381 : bates 478 SEXP dgeMatrix_exp(SEXP x)
382 : bates 461 {
383 :     SEXP val = PROTECT(duplicate(x));
384 :     int *Dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
385 :     int i, ilo, ilos, ihi, ihis, j, nc = Dims[1], sqpow;
386 :     int ncp1 = Dims[1] + 1, ncsqr = nc * nc;
387 :     int *pivot = Calloc(nc, int);
388 :     int *iperm = Calloc(nc, int);
389 :     double *dpp = Calloc(ncsqr, double), /* denominator power Pade' */
390 :     *npp = Calloc(ncsqr, double), /* numerator power Pade' */
391 :     *perm = Calloc(nc, double),
392 :     *scale = Calloc(nc, double),
393 :     *v = REAL(GET_SLOT(val, Matrix_xSym)),
394 :     *work = Calloc(ncsqr, double), inf_norm, m1_j, /* (-1)^j */
395 :     one = 1., trshift, zero = 0.;
396 : maechler 534
397 : bates 461 if (nc < 1 || Dims[0] != nc)
398 : bates 582 error(_("Matrix exponential requires square, non-null matrix"));
399 : bates 461
400 :     /* FIXME: Add special treatment for nc == 1 */
401 : maechler 534
402 : bates 461 /* Preconditioning 1. Shift diagonal by average diagonal if positive. */
403 :     trshift = 0; /* determine average diagonal element */
404 :     for (i = 0; i < nc; i++) trshift += v[i * ncp1];
405 :     trshift /= nc;
406 :     if (trshift > 0.) { /* shift diagonal by -trshift */
407 :     for (i = 0; i < nc; i++) v[i * ncp1] -= trshift;
408 :     }
409 :    
410 :     /* Preconditioning 2. Balancing with dgebal. */
411 :     F77_CALL(dgebal)("P", &nc, v, &nc, &ilo, &ihi, perm, &j);
412 : bates 582 if (j) error(_("dgeMatrix_exp: LAPACK routine dgebal returned %d"), j);
413 : bates 461 F77_CALL(dgebal)("S", &nc, v, &nc, &ilos, &ihis, scale, &j);
414 : bates 582 if (j) error(_("dgeMatrix_exp: LAPACK routine dgebal returned %d"), j);
415 : bates 461
416 :     /* Preconditioning 3. Scaling according to infinity norm */
417 :     inf_norm = F77_CALL(dlange)("I", &nc, &nc, v, &nc, work);
418 :     sqpow = (inf_norm > 0) ? (int) (1 + log(inf_norm)/log(2.)) : 0;
419 :     if (sqpow < 0) sqpow = 0;
420 :     if (sqpow > 0) {
421 :     double scale_factor = 1.0;
422 :     for (i = 0; i < sqpow; i++) scale_factor *= 2.;
423 :     for (i = 0; i < ncsqr; i++) v[i] /= scale_factor;
424 :     }
425 :    
426 :     /* Pade' approximation. Powers v^8, v^7, ..., v^1 */
427 :     AZERO(npp, ncsqr);
428 :     AZERO(dpp, ncsqr);
429 :     m1_j = -1;
430 :     for (j = 7; j >=0; j--) {
431 :     double mult = padec[j];
432 :     /* npp = m * npp + padec[j] *m */
433 :     F77_CALL(dgemm)("N", "N", &nc, &nc, &nc, &one, v, &nc, npp, &nc,
434 :     &zero, work, &nc);
435 :     for (i = 0; i < ncsqr; i++) npp[i] = work[i] + mult * v[i];
436 :     /* dpp = m * dpp * (m1_j * padec[j]) * m */
437 :     mult *= m1_j;
438 :     F77_CALL(dgemm)("N", "N", &nc, &nc, &nc, &one, v, &nc, dpp, &nc,
439 :     &zero, work, &nc);
440 :     for (i = 0; i < ncsqr; i++) dpp[i] = work[i] + mult * v[i];
441 :     m1_j *= -1;
442 :     }
443 :     /* Zero power */
444 :     for (i = 0; i < ncsqr; i++) dpp[i] *= -1.;
445 :     for (j = 0; j < nc; j++) {
446 :     npp[j * ncp1] += 1.;
447 :     dpp[j * ncp1] += 1.;
448 :     }
449 : maechler 534
450 : bates 461 /* Pade' approximation is solve(dpp, npp) */
451 :     F77_CALL(dgetrf)(&nc, &nc, dpp, &nc, pivot, &j);
452 : bates 582 if (j) error(_("dgeMatrix_exp: dgetrf returned error code %d"), j);
453 : bates 461 F77_CALL(dgetrs)("N", &nc, &nc, dpp, &nc, pivot, npp, &nc, &j);
454 : bates 582 if (j) error(_("dgeMatrix_exp: dgetrs returned error code %d"), j);
455 : bates 461 Memcpy(v, npp, ncsqr);
456 :    
457 :     /* Now undo all of the preconditioning */
458 :     /* Preconditioning 3: square the result for every power of 2 */
459 :     while (sqpow--) {
460 :     F77_CALL(dgemm)("N", "N", &nc, &nc, &nc, &one, v, &nc, v, &nc,
461 :     &zero, work, &nc);
462 :     Memcpy(v, work, ncsqr);
463 :     }
464 :     /* Preconditioning 2: apply inverse scaling */
465 :     for (j = 0; j < nc; j++)
466 :     for (i = 0; i < nc; i++)
467 :     v[i + j * nc] *= scale[i]/scale[j];
468 :     /* Construct balancing permutation vector */
469 :     for (i = 0; i < nc; i++) iperm[i] = i; /* identity permutation */
470 :     /* Leading permutations applied in forward order */
471 :     for (i = 0; i < (ilo - 1); i++) {
472 :     int swapidx = (int) (perm[i]) - 1;
473 :     int tmp = iperm[i];
474 :     iperm[i] = iperm[swapidx];
475 :     iperm[swapidx] = tmp;
476 :     }
477 :     /* Trailing permutations applied in reverse order */
478 :     for (i = nc - 1; i >= ihi; i--) {
479 :     int swapidx = (int) (perm[i]) - 1;
480 :     int tmp = iperm[i];
481 :     iperm[i] = iperm[swapidx];
482 :     iperm[swapidx] = tmp;
483 :     }
484 :     /* Construct inverse balancing permutation vector */
485 :     Memcpy(pivot, iperm, nc);
486 :     for (i = 0; i < nc; i++) iperm[pivot[i]] = i;
487 :     /* Apply inverse permutation */
488 :     Memcpy(work, v, ncsqr);
489 :     for (j = 0; j < nc; j++)
490 :     for (i = 0; i < nc; i++)
491 :     v[i + j * nc] = work[iperm[i] + iperm[j] * nc];
492 : maechler 534
493 : bates 461 /* Preconditioning 1: Trace normalization */
494 :     if (trshift > 0.) {
495 :     double mult = exp(trshift);
496 :     for (i = 0; i < ncsqr; i++) v[i] *= mult;
497 :     }
498 : maechler 534
499 : bates 461 /* Clean up */
500 :     Free(dpp); Free(npp); Free(perm); Free(iperm); Free(pivot); Free(scale); Free(work);
501 :     UNPROTECT(1);
502 :     return val;
503 :     }
504 : bates 495
505 :     SEXP dgeMatrix_Schur(SEXP x, SEXP vectors)
506 :     {
507 :     int *dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
508 :     int vecs = asLogical(vectors), info, izero = 0, lwork = -1, n = dims[0];
509 :     double *work, tmp;
510 :     char *nms[] = {"WR", "WI", "T", "Z", ""};
511 :     SEXP val = PROTECT(Matrix_make_named(VECSXP, nms));
512 : maechler 534
513 : bates 495 if (n != dims[1] || n < 1)
514 : bates 582 error(_("dgeMatrix_Schur: argument x must be a non-null square matrix"));
515 : bates 495 SET_VECTOR_ELT(val, 0, allocVector(REALSXP, n));
516 :     SET_VECTOR_ELT(val, 1, allocVector(REALSXP, n));
517 :     SET_VECTOR_ELT(val, 2, allocMatrix(REALSXP, n, n));
518 :     Memcpy(REAL(VECTOR_ELT(val, 2)), REAL(GET_SLOT(x, Matrix_xSym)), n * n);
519 :     SET_VECTOR_ELT(val, 3, allocMatrix(REALSXP, vecs ? n : 0, vecs ? n : 0));
520 :     F77_CALL(dgees)(vecs ? "V" : "N", "N", NULL, dims, (double *) NULL, dims, &izero,
521 :     (double *) NULL, (double *) NULL, (double *) NULL, dims,
522 :     &tmp, &lwork, (int *) NULL, &info);
523 : bates 582 if (info) error(_("dgeMatrix_Schur: first call to dgees failed"));
524 : bates 495 lwork = (int) tmp;
525 :     work = Calloc(lwork, double);
526 :     F77_CALL(dgees)(vecs ? "V" : "N", "N", NULL, dims, REAL(VECTOR_ELT(val, 2)), dims,
527 :     &izero, REAL(VECTOR_ELT(val, 0)), REAL(VECTOR_ELT(val, 1)),
528 : maechler 534 REAL(VECTOR_ELT(val, 3)), dims, work, &lwork,
529 : bates 495 (int *) NULL, &info);
530 : bates 582 if (info) error(_("dgeMatrix_Schur: dgees returned code %d"), info);
531 : bates 495 Free(work);
532 :     UNPROTECT(1);
533 :     return val;
534 :     }

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge