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, size_t n, IVector indx, int mode, double *d)
  12 {
  13   int i, imax, j, k;
  14   int singular = FALSE;
  15   double big, temp;
  16   Complex cdum, csum;
  17   double rdum, rsum;
  18   CMatrix cmat = (CMatrix) mat;
  19   RMatrix rmat = (RMatrix) mat;
  20   RVector vv;
  21
  22   vv = rvector(n);
  23   *d = 1.0;
  24
  25   /* set up the pivot permutation vector */
  26   for (i = 0; i < n; i++) indx[i] = i;
  27
  28   /* get scaling information for implicit pivoting */
  29   for (i = 0; i < n; i++) {
  30     big = 0.0;
  31     for (j = 0; j < n; j++) {
  32       temp = (mode == RE) ? fabs(rmat[i][j]) : modulus(cmat[i][j]);
  33       if (temp > big) big = temp;
  34     }
  35     if (big == 0.0) {
  36       vv[i] = 1.0;                  /* no scaling for zero rows */
  37       singular = TRUE;
  38     }
  39     else vv[i] = 1.0 / big;
  40   }
  41
  42   /* loop over columns for Crout's method */
  43   for (j = 0; j < n; j++) {
  44     for (i = 0; i < j; i++) {
  45       if (mode == RE) rsum = rmat[i][j];
  46       else csum = cmat[i][j];
  47
  48       for (k = 0; k < i; k++)
  49       {
  50         if (mode == RE) rsum -= rmat[i][k] * rmat[k][j];
  51         else csum = csub(csum, cmul(cmat[i][k], cmat[k][j]));
  52       }
  53
  54       if (mode == RE) rmat[i][j] = rsum;
  55       else cmat[i][j] = csum;
  56     }
  57     big = 0.0;
  58     for (i = j; i < n; i++) {
  59       if (mode == RE) rsum = rmat[i][j];
  60       else csum = cmat[i][j];
  61
  62       for (k = 0; k < j; k++)
  63       {
  64         if (mode == RE) rsum -= rmat[i][k] * rmat[k][j];
  65         else csum = csub(csum, cmul(cmat[i][k], cmat[k][j]));
  66       }
  67
  68       if (mode == RE) rmat[i][j] = rsum;
  69       else cmat[i][j] = csum;
  70
  71       temp = vv[i] * ((mode == RE) ? fabs(rsum) : modulus(csum));
  72       if (temp >= big) { big = temp; imax = i; }
  73     }
  74
  75     /* interchange rows if needed */
  76     if (j != imax) {
  77       for (k = 0; k < n; k++) {
  78         if (mode == RE) {
  79           rdum = rmat[imax][k];
  80           rmat[imax][k] = rmat[j][k];
  81           rmat[j][k] = rdum;
  82         }
  83         else {
  84           cdum = cmat[imax][k];
  85           cmat[imax][k] = cmat[j][k];
  86           cmat[j][k] = cdum;
  87         }
  88       }
  89       *d = -(*d);
  90       vv[imax] = vv[j];
  91     }
  92     indx[j] = imax;
  93
  94     /* divide by the pivot element */
  95     temp = (mode == RE) ? fabs(rmat[j][j]) : modulus(cmat[j][j]);
  96     if (temp == 0.0) singular = TRUE;
  97     else if (j < n - 1) {
  98       if (mode == RE) {
  99         rdum = 1.0 / rmat[j][j];
 100         for (i = j + 1; i < n; i++) rmat[i][j] *= rdum;
 101       }
 102       else {
 103         cdum = cdiv(cart2complex(1.0, 0.0), cmat[j][j]);
 104         for (i = j + 1; i < n; i++)
 105         {
 106           cmat[i][j] = cmul(cmat[i][j], cdum);
 107         }
 108       }
 109     }
 110   }
 111
 112   free_vector(vv);
 113   return(singular);
 114 }
 115
 116 int
 117 crlubksb(double **a, int n, int *indx, double *b, int mode)
 118 {
 119   int i, ii, ip, j, singular = FALSE;
 120   CMatrix ca = (CMatrix) a;
 121   CVector cb = (CVector) b;
 122   RMatrix ra = (RMatrix) a;
 123   RVector rb = (RVector) b;
 124   double rsum;
 125   Complex csum;
 126
 127   /* forward substitute using L part */
 128   for (i = 0, ii = -1; i < n; i++) {
 129     ip = indx[i];
 130     if (mode == RE) {
 131       rsum = rb[ip];
 132       rb[ip] = rb[i];
 133     }
 134     else {
 135       csum = cb[ip];
 136       cb[ip] = cb[i];
 137     }
 138     if (ii >= 0)
 139       for (j = ii; j <= i - 1; j++)
 140         if (mode == RE) rsum -= ra[i][j] * rb[j];
 141         else csum = csub(csum, cmul(ca[i][j], cb[j]));
 142     else {
 143       if (mode == RE) {
 144         if (rsum != 0.0) ii = i;
 145       }
 146       else if (csum.real != 0.0 || csum.imag != 0.0) ii = i;
 147     }
 148     if (mode == RE) rb[i] = rsum;
 149     else cb[i] = csum;
 150   }
 151
 152   /* back substitute using the U part */
 153   for (i = n - 1; i >= 0; i--) {
 154     if (mode == RE) {
 155       rsum = rb[i];
 156       for (j = i + 1; j < n; j++) rsum -= ra[i][j] * rb[j];
 157       if (ra[i][i] == 0.0) {
 158         singular = TRUE;
 159         break;
 160       }
 161       else rb[i] = rsum / ra[i][i];
 162     }
 163     else {
 164       csum = cb[i];
 165       for (j = i + 1; j < n; j++) csum = csub(csum, cmul(ca[i][j], cb[j]));
 166       if (modulus(ca[i][i]) == 0.0) {
 167         singular = TRUE;
 168         break;
 169       }
 170       else cb[i] = cdiv(csum, ca[i][i]);
 171     }
 172   }
 173
 174   return(singular);
 175 }
 176