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