blob 03c5e1775dbef4176119c312050fae6bd177a9db
1 /* rcondest - Estimates reciprocal of condition number. */
2 /* XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney */
3 /* Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz */
4 /* You may give out copies of this software; for conditions see the */
5 /* file COPYING included with this distribution. */
7 #include "linalg.h"
9 static
10 double Max(double a, double b)
12 return(a > b ? a : b);
16 void
17 backsolve(double **a, double *b, size_t n, int mode)
19 long i, j; /* Must be signed, since counting down in loops... */
20 RMatrix ra = (RMatrix) a;
21 RVector rb = (RVector) b;
22 CMatrix ca = (CMatrix) a;
23 CVector cb = (CVector) b;
25 switch (mode) {
26 case RE:
27 for (i = n - 1; i >= 0; i--) {
28 if (ra[i][i] != 0.0) rb[i] = rb[i] / ra[i][i];
29 for (j = i + 1; j < n; j++) rb[i] -= ra[i][j] * rb[j];
31 break;
32 case CX:
33 for (i = n - 1; i >= 0; i--) {
34 if (modulus(ca[i][i]) != 0.0) cb[i] = cdiv(cb[i], ca[i][i]);
35 for (j = i + 1; j < n; j++)
36 cb[i] = csub(cb[i], cmul(ca[i][j], cb[j]));
38 break;
42 /** Exported function **/
44 double rcondest(double **a, size_t n)
46 RVector p, pm, x;
47 double est, xp, xm, temp, tempm, xnorm;
48 long i, j;
50 for (i = 0; i < n; i++)
51 if (a[i][i] == 0.0) return(0.0);
53 p = rvector(n);
54 pm = rvector(n);
55 x = rvector(n);
57 /* Set est to reciprocal of L1 matrix norm of A */
58 est = fabs(a[0][0]);
60 for (j = 1; j < n; j++) {
61 for (i = 0, temp = fabs(a[j][j]); i < j; i++)
62 temp += fabs(a[i][j]);
63 est = Max(est, temp);
65 est = 1.0 / est;
67 /* Solve A^Tx = e, selecting e as you proceed */
68 x[0] = 1.0 / a[0][0];
69 for (i = 1; i < n; i++) p[i] = a[0][i] * x[0];
70 for (j = 1; j < n; j++) {
71 /* select ej and calculate x[j] */
72 xp = ( 1.0 - p[j]) / a[j][j];
73 xm = (-1.0 - p[j]) / a[j][j];
74 temp = fabs(xp);
75 tempm = fabs(xm);
76 for (i = j + 1; i < n; i++) {
77 pm[i] = p[i] + a[j][i] * xm;
78 tempm += fabs(pm[i] / a[i][i]);
79 p[i] += a[j][i] * xp;
80 temp += fabs(p[i] / a[i][i]);
82 if (temp >= tempm) x[j] = xp;
83 else {
84 x[j] = xm;
85 for (i = j + 1; i < n; i++) p[i] = pm[i];
89 for (j = 0, xnorm = 0.0; j < n; j++) xnorm += fabs(x[j]);
90 est = est * xnorm;
92 backsolve(a, x, n, RE);
93 for (j = 0, xnorm = 0.0; j < n; j++) xnorm += fabs(x[j]);
94 if (xnorm > 0) est = est / xnorm;
96 free_vector(p);
97 free_vector(pm);
98 free_vector(x);
100 return(est);