lib/ludecomp.c

   1 /* ludecomp - LU decomposition and backsolving routines.               */
   2 /* Taken from Numerical Recipies.                                      */
   3 /* XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney                  */
   4 /* Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz    */
   5 /* You may give out copies of this software; for conditions see the    */
   6 /* file COPYING included with this distribution.                       */
   7
   8 #include "linalg.h"
   9
  10 int
  11 crludcmp(Matrix mat, int n, IVector indx, int mode, double *d)
  12 {
  13   int i, imax, j, k, singular = FALSE;
  14   double big, temp;
  15   Complex cdum, csum;
  16   double rdum, rsum;
  17   CMatrix cmat = (CMatrix) mat;
  18   RMatrix rmat = (RMatrix) mat;
  19   RVector vv;
  20
  21   vv = rvector(n);
  22   *d = 1.0;
  23
  24   /* set up the pivot permutation vector */
  25   for (i = 0; i < n; i++) indx[i] = i;
  26
  27   /* get scaling information for implicit pivoting */
  28   for (i = 0; i < n; i++) {
  29     big = 0.0;
  30     for (j = 0; j < n; j++) {
  31       temp = (mode == RE) ? fabs(rmat[i][j]) : modulus(cmat[i][j]);
  32       if (temp > big) big = temp;
  33     }
  34     if (big == 0.0) {
  35       vv[i] = 1.0;                  /* no scaling for zero rows */
  36       singular = TRUE;
  37     }
  38     else vv[i] = 1.0 / big;
  39   }
  40
  41   /* loop over columns for Crout's method */
  42   for (j = 0; j < n; j++) {
  43     for (i = 0; i < j; i++) {
  44       if (mode == RE) rsum = rmat[i][j];
  45       else csum = cmat[i][j];
  46
  47       for (k = 0; k < i; k++)
  48         if (mode == RE) rsum -= rmat[i][k] * rmat[k][j];
  49         else csum = csub(csum, cmul(cmat[i][k], cmat[k][j]));
  50
  51       if (mode == RE) rmat[i][j] = rsum;
  52       else cmat[i][j] = csum;
  53     }
  54     big = 0.0;
  55     for (i = j; i < n; i++) {
  56       if (mode == RE) rsum = rmat[i][j];
  57       else csum = cmat[i][j];
  58
  59       for (k = 0; k < j; k++)
  60         if (mode == RE) rsum -= rmat[i][k] * rmat[k][j];
  61         else csum = csub(csum, cmul(cmat[i][k], cmat[k][j]));
  62
  63       if (mode == RE) rmat[i][j] = rsum;
  64       else cmat[i][j] = csum;
  65
  66       temp = vv[i] * ((mode == RE) ? fabs(rsum) : modulus(csum));
  67       if (temp >= big) { big = temp; imax = i; }
  68     }
  69
  70     /* interchange rows if needed */
  71     if (j != imax) {
  72       for (k = 0; k < n; k++) {
  73         if (mode == RE) {
  74           rdum = rmat[imax][k];
  75           rmat[imax][k] = rmat[j][k];
  76           rmat[j][k] = rdum;
  77         }
  78         else {
  79           cdum = cmat[imax][k];
  80           cmat[imax][k] = cmat[j][k];
  81           cmat[j][k] = cdum;
  82         }
  83       }
  84       *d = -(*d);
  85       vv[imax] = vv[j];
  86     }
  87     indx[j] = imax;
  88
  89     /* divide by the pivot element */
  90     temp = (mode == RE) ? fabs(rmat[j][j]) : modulus(cmat[j][j]);
  91     if (temp == 0.0) singular = TRUE;
  92     else if (j < n - 1) {
  93       if (mode == RE) {
  94         rdum = 1.0 / rmat[j][j];
  95         for (i = j + 1; i < n; i++) rmat[i][j] *= rdum;
  96       }
  97       else {
  98         cdum = cdiv(cart2complex(1.0, 0.0), cmat[j][j]);
  99         for (i = j + 1; i < n; i++) cmat[i][j] = cmul(cmat[i][j], cdum);
 100       }
 101     }
 102   }
 103   free_vector(vv);
 104   return(singular);
 105 }
 106
 107 int
 108 crlubksb(double **a, int n, int *indx, double *b, int mode)
 109 {
 110   int i, ii, ip, j, singular = FALSE;
 111   CMatrix ca = (CMatrix) a;
 112   CVector cb = (CVector) b;
 113   RMatrix ra = (RMatrix) a;
 114   RVector rb = (RVector) b;
 115   double rsum;
 116   Complex csum;
 117
 118   /* forward substitute using L part */
 119   for (i = 0, ii = -1; i < n; i++) {
 120     ip = indx[i];
 121     if (mode == RE) {
 122       rsum = rb[ip];
 123       rb[ip] = rb[i];
 124     }
 125     else {
 126       csum = cb[ip];
 127       cb[ip] = cb[i];
 128     }
 129     if (ii >= 0)
 130       for (j = ii; j <= i - 1; j++)
 131         if (mode == RE) rsum -= ra[i][j] * rb[j];
 132         else csum = csub(csum, cmul(ca[i][j], cb[j]));
 133     else {
 134       if (mode == RE) {
 135         if (rsum != 0.0) ii = i;
 136       }
 137       else if (csum.real != 0.0 || csum.imag != 0.0) ii = i;
 138     }
 139     if (mode == RE) rb[i] = rsum;
 140     else cb[i] = csum;
 141   }
 142
 143   /* back substitute using the U part */
 144   for (i = n - 1; i >= 0; i--) {
 145     if (mode == RE) {
 146       rsum = rb[i];
 147       for (j = i + 1; j < n; j++) rsum -= ra[i][j] * rb[j];
 148       if (ra[i][i] == 0.0) {
 149         singular = TRUE;
 150         break;
 151       }
 152       else rb[i] = rsum / ra[i][i];
 153     }
 154     else {
 155       csum = cb[i];
 156       for (j = i + 1; j < n; j++) csum = csub(csum, cmul(ca[i][j], cb[j]));
 157       if (modulus(ca[i][i]) == 0.0) {
 158         singular = TRUE;
 159         break;
 160       }
 161       else cb[i] = cdiv(csum, ca[i][i]);
 162     }
 163   }
 164
 165   return(singular);
 166 }
 167