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