SCM

SCM Repository

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

Annotation of /pkg/src/geMatrix.c

Parent Directory Parent Directory | Revision Log Revision Log


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

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