lib/rcondest.c

   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.                       */
   6
   7 #include "linalg.h"
   8
   9 static
  10 double Max(double a, double b)
  11 {
  12   return(a > b ? a : b);
  13 }
  14
  15
  16 void
  17 backsolve(double **a, double *b, int n, int mode)
  18 {
  19   int i, j;
  20   RMatrix ra = (RMatrix) a;
  21   RVector rb = (RVector) b;
  22   CMatrix ca = (CMatrix) a;
  23   CVector cb = (CVector) b;
  24
  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];
  30     }
  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]));
  37     }
  38     break;
  39   }
  40 }
  41
  42 /** Exported function **/
  43
  44 double rcondest(double **a, int n)
  45 {
  46   RVector p, pm, x;
  47   double est, xp, xm, temp, tempm, xnorm;
  48   int i, j;
  49
  50   for (i = 0; i < n; i++)
  51     if (a[i][i] == 0.0) return(0.0);
  52
  53   p = rvector(n);
  54   pm = rvector(n);
  55   x = rvector(n);
  56
  57   /* Set est to reciprocal of L1 matrix norm of A */
  58   est = fabs(a[0][0]);
  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]);
  62     est = Max(est, temp);
  63   }
  64   est = 1.0 / est;
  65
  66   /* Solve A^Tx = e, selecting e as you proceed */
  67   x[0] = 1.0 / a[0][0];
  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];
  73     temp = fabs(xp);
  74     tempm = fabs(xm);
  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]);
  78       p[i] += a[j][i] * xp;
  79       temp += fabs(p[i] / a[i][i]);
  80     }
  81     if (temp >= tempm) x[j] = xp;
  82     else {
  83       x[j] = xm;
  84       for (i = j + 1; i < n; i++) p[i] = pm[i];
  85     }
  86   }
  87
  88   for (j = 0, xnorm = 0.0; j < n; j++) xnorm += fabs(x[j]);
  89   est = est * xnorm;
  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;
  93
  94   free_vector(p);
  95   free_vector(pm);
  96   free_vector(x);
  97
  98   return(est);
  99 }
 100