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. */
10 double Max(double a
, double b
)
12 return(a
> b
? a
: b
);
17 backsolve(double **a
, double *b
, int n
, int mode
)
20 RMatrix ra
= (RMatrix
) a
;
21 RVector rb
= (RVector
) b
;
22 CMatrix ca
= (CMatrix
) a
;
23 CVector cb
= (CVector
) b
;
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
];
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
]));
42 /** Exported function **/
44 double rcondest(double **a
, int n
)
47 double est
, xp
, xm
, temp
, tempm
, xnorm
;
50 for (i
= 0; i
< n
; i
++)
51 if (a
[i
][i
] == 0.0) return(0.0);
57 /* Set est to reciprocal of L1 matrix norm of A */
59 for (j
= 1; j
< n
; j
++) {
60 for (i
= 0, temp
= fabs(a
[j
][j
]); i
< j
; i
++)
61 temp
+= fabs(a
[i
][j
]);
66 /* Solve A^Tx = e, selecting e as you proceed */
68 for (i
= 1; i
< n
; i
++) p
[i
] = a
[0][i
] * x
[0];
69 for (j
= 1; j
< n
; j
++) {
70 /* select ej and calculate x[j] */
71 xp
= ( 1.0 - p
[j
]) / a
[j
][j
];
72 xm
= (-1.0 - p
[j
]) / a
[j
][j
];
75 for (i
= j
+ 1; i
< n
; i
++) {
76 pm
[i
] = p
[i
] + a
[j
][i
] * xm
;
77 tempm
+= fabs(pm
[i
] / a
[i
][i
]);
79 temp
+= fabs(p
[i
] / a
[i
][i
]);
81 if (temp
>= tempm
) x
[j
] = xp
;
84 for (i
= j
+ 1; i
< n
; i
++) p
[i
] = pm
[i
];
88 for (j
= 0, xnorm
= 0.0; j
< n
; j
++) xnorm
+= fabs(x
[j
]);
90 backsolve(a
, x
, n
, RE
);
91 for (j
= 0, xnorm
= 0.0; j
< n
; j
++) xnorm
+= fabs(x
[j
]);
92 if (xnorm
> 0) est
= est
/ xnorm
;