home *** CD-ROM | disk | FTP | other *** search
/ The C Users' Group Library 1994 August / wc-cdrom-cusersgrouplibrary-1994-08.iso / listings / v_08_12 / 8n12104a < prev    next >
Text File  |  1990-05-27  |  4KB  |  134 lines
  1.  
  2.  
  3. #include "a:\include\stdio.h"
  4. #include "a:\include\io.h"
  5. #include "a:\include\fcntl.h"
  6. #include "a:\include\sys\stat.h"
  7. #include "a:\include\math.h"
  8. #include "a:\include\float.h"
  9. #include "a:\include\malloc.h"
  10.  
  11. #define SWAP(a, b) {float temp=(a); (a)=(b); (b)=temp;}
  12.  
  13.  
  14.  
  15. void gaussj(a, n, b, m)
  16. float **a, **b;
  17. int n,m;
  18.  
  19. /* Linear equation solution by Guass-Jordan elimination, equation (2.1.1)
  20.  above. a[1..n][1..n] is an input matrix of n by n elements. b[1..n][1..m]
  21.  is an input matrix of size n x m containing the m right hand side vectors.
  22.  On output, a is replace by its matrix inverse, and b is replaced by the
  23.  corresponding set of solution vectors. */
  24.  
  25.  {
  26.     int *indxc, *indxr, *ipiv;
  27.     /* The integer arrays ipiv[1..n], indxr[1..n] and indxc[1..n] are
  28.     used for bookkeeping on the pivoting */
  29.     int i, icol, irow, j, k, l, ll, *ivector();
  30.     float big, dum, pivinv;
  31.     void nrerror(), free_ivector();
  32.  
  33.     indxc=ivector(1,n);
  34.     indxr=ivector(1,n);
  35.     ipiv=ivector(1,n);
  36.     for (j=1;j<n;j++) ipiv[j]=0;
  37.     for (i=1;i<=n;i++)  { /* This is the main loop over the columns
  38.                  to be reduced. */
  39.        big=0.0;
  40.        for (j=1;j<=n;j++)  /* This is the outer loop of the search
  41.                   for a pivot element. */
  42.       if (ipiv[j] != 1)
  43.          for (k=1;k<=n;k++) {
  44.         if (ipiv[k] == 0) {
  45.            if (fabs(a[j][k]) >= big) {
  46.               big=fabs(a[j][k]);
  47.               irow=j;
  48.               icol=k;
  49.            }
  50.         } else if (ipiv[k] > 1) nrerror("GUASSJ: Singular Matrix-1");
  51.          }
  52.        ++(ipiv[icol]);
  53.        /* We now have the pivot element, so we interchange rows, if needed,
  54.        to put the pivot element on the diagonal.  The columns are not
  55.        physically interchanged, only relabeled:indx[i], the column of
  56.        the ith pivot element, is the ith column that is reduced, while
  57.        indxr[i] is the row in which that pivot element was originally
  58.        located.  If indxr[i] != indxc[i] there is an implied column
  59.        interchange.  With this form of bookkeeping, the solution b's
  60.        will end up in the correct order, and the inverse matrix will be
  61.        scrambled by columns.
  62.        if (irow != icol) {
  63.       for (l=1;l<=n;l++) SWAP(a[irow][l],a[icol][l])
  64.       for (l=1;l<=m;l++) SWAP(b[irow][l],b[icol][l])
  65.        }
  66.        indxr[i]=irow;
  67.        indxc[i]=icol;  /* We are now ready to divide the pivot row by
  68.               the pivot element, located at irow and icol */
  69.        if (a[icol][icol] == 0.0) nrerror("GUASSJ:Singular Matrix-2");
  70.        pivinv=1.0/a[icol][icol];
  71.        a[icol][icol]=1.0;
  72.        for (l=1;l<=n;l++) a[icol][l] *= pivinv;
  73.        for (l=1;l<=m;l++) b[icol][l] *= pivinv;
  74.        for (ll=1;ll<=n;ll++) /* Next we reduce the rows, except for the
  75.                 pivot one, of course */
  76.       if (ll != icol) {
  77.          dum=a[ll][icol];
  78.          a[ll][icol]=0.0;
  79.          for (l=1;l<=n;l++) a[ll][l] -= a[icol][l]*dum;
  80.          for (l=1;l<=m;l++) b[ll][l] -= b[icol][l]*dum;
  81.       }
  82.     } /* This is the end of the main loop over columns of the reduction.
  83.      It only remains to unscramble the solution in view of the
  84.      column interchanges.  We do this by interchanging pairs of
  85.      columns in the reverse order that the permutation was
  86.      built up. */
  87.     for (l=n;l>=1;l--) {
  88.        if (indxr[l] != indxc[l])
  89.       for (k=1;k<=n;k++)
  90.          SWAP(a[k][indxr[l]],a[k][indxc[l]]);
  91.     } /* And we are done. */
  92.     free_ivector(ipiv,1,n);
  93.     free_ivector(indxr,1,n);
  94.     free_ivector(indxc,1,n);
  95.  }
  96.  
  97.  
  98.  
  99.  
  100.  void nrerror(error_text)
  101.  char error_text[];
  102.  /* Numerical Recipes standard error handler */
  103.  {
  104.     void exit();
  105.  
  106.     fprintf(stderr,"Numerical Recipes run-time error...\n");
  107.     fprintf(stderr,"%s\n",error_text);
  108.     fprintf(stderr,"...now exiting to system...\n");
  109.     exit(1);
  110.  }
  111.  
  112.  
  113.  
  114.  
  115.  int *ivector(nl,nh)
  116.  int nl,nh;
  117.  /* Allocates an int vector with range [nl..nh] */
  118.  {
  119.     int *v;
  120.  
  121.     v=(int *)malloc((unsigned) (nh-nl+1)*sizeof(int));
  122.     if (!v) nrerror("allocation failure in ivector()");
  123.     return v-nl;
  124.  }
  125.  
  126.  
  127.  
  128.  
  129.  void free_ivector(v,nl,nh)
  130.  int *v,nl,nh;
  131.  /* Frees an int vector allocated by ivector() */
  132.  {
  133.     free((char*) (v+nl));
  134.  }