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, size_t n, int mode)
  18 {
  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;
  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, size_t n)
  45 {
  46   RVector p, pm, x;
  47   double est, xp, xm, temp, tempm, xnorm;
  48   long 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
  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);
  64   }
  65   est = 1.0 / est;
  66
  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]);
  81     }
  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];
  86     }
  87   }
  88
  89   for (j = 0, xnorm = 0.0; j < n; j++) xnorm += fabs(x[j]);
  90   est = est * xnorm;
  91
  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;
  95
  96   free_vector(p);
  97   free_vector(pm);
  98   free_vector(x);
  99
 100   return(est);
 101 }
 102