home *** CD-ROM | disk | FTP | other *** search
/ vis-ftp.cs.umass.edu / vis-ftp.cs.umass.edu.tar / vis-ftp.cs.umass.edu / pub / Software / ASCENDER / ascendMar8.tar / UMass / Triangulate / gaussj_check.c < prev    next >
C/C++ Source or Header  |  1995-04-13  |  4KB  |  137 lines
  1. /* @(#)gaussj.c    1.5 6/17/94 */
  2. /* Gauss Jordan Elimination */
  3.  
  4. #include "cvar.h"
  5. #include <math.h>
  6. #include "alloc.h"
  7. #include "error_message.s"
  8.  
  9. #define SWAP(a,b) {double temp = (a); (a)=(b); (b) = temp;}
  10.  
  11. FUNCTION_DEF ( int gaussj_check, (
  12.         double **a,
  13.         int n,
  14.         double *b))
  15.    {
  16.    /* Linear equation solution to  Gauss-Jordan elimination.
  17.       a[1..n][1..n] is an input matrix of n by n elements.
  18.       b[1..n] is an input matrix of size n containing the 
  19.       right-hand side vector.
  20.       On output, a is replaced by its matrix inverse, and b is replaced 
  21.       by the corresponding solution vector.
  22.     */
  23.    int *indxc, *indxr, *ipiv;    /* Integer arrays used for book-keeping on the
  24.                  * pivoting */
  25.    int i, icol, irow, j, k, l, ll;
  26.    double big, dum, pivinv;
  27.  
  28.    /* Declare some temporary arrays */
  29.    indxc = VECTOR(1, n, int);
  30.    indxr = VECTOR(1, n, int);
  31.    ipiv =  VECTOR(1, n, int);
  32.  
  33.    for (j=1; j<=n; j++) ipiv[j] = 0;
  34.    for (i=1; i<=n; i++)
  35.       {
  36.       /* This is the main loop over the columns to be reduced */
  37.       big = 0.0;
  38.       for (j=1; j<=n; j++)
  39.      {
  40.      /* This is the outer loop of the search for the pivot element */
  41.      if (ipiv[j] != 1)
  42.         {
  43.         for (k=1; k<=n; k++)
  44.            {
  45.            if (ipiv[k] == 0)
  46.           {
  47.           if (fabs(a[j][k]) >= big)
  48.              {
  49.                 big = fabs(a[j][k]);
  50.              irow = j;
  51.              icol = k;
  52.              }    
  53.           }
  54.            else if (ipiv[k] > 1) 
  55.          {
  56.            error_message("GAUSSJ : Singular Matrix - 1");
  57.            FREE (indxc, 1);
  58.            FREE (indxr, 1);
  59.            FREE (ipiv, 1);
  60.            return(FALSE);
  61.          }
  62.            }
  63.         }
  64.      }
  65.       ++(ipiv[icol]);
  66.  
  67.       /* We now have the pivot element, so we interchange rows, if 
  68.        * needed, to put the pivot element on the diagonal.  The columns 
  69.        * are not physically interchanged, only relabelled : indxr[i], 
  70.        * the colun of the i-th pivot element, is the i-th column, while 
  71.        * indxr[i] is the row in which that pivot element was originally 
  72.        * located.  If indxr[i] != indxc[i] there is an implied column 
  73.        * interchange.  With this form of book-keeping, the solutions b 
  74.        * will end up in the correct order and the inverse matrix will 
  75.        * be scrambled by columns.
  76.        */
  77.  
  78.       if (irow != icol)
  79.          {
  80.          for (l=1; l<=n; l++) SWAP(a[irow][l], a[icol][l]);
  81.          SWAP(b[irow], b[icol]);
  82.          }
  83.  
  84.       /* We are now ready to divide the pivot row by
  85.        * the pivot element, located ate irow, icol */
  86.       indxr[i] = irow;
  87.       indxc[i] = icol;
  88.       if (a[icol][icol] == 0.0) 
  89.     {
  90.       error_message("GAUSSJ : Singular Matrix - 2");
  91.       FREE (indxc, 1);
  92.       FREE (indxr, 1);
  93.       FREE (ipiv, 1);
  94.       return(FALSE);
  95.     }
  96.       pivinv = 1.0 / a[icol][icol];
  97.       a[icol][icol] = 1.0;
  98.       for (l=1; l<=n; l++) a[icol][l] *= pivinv;
  99.       b[icol] *= pivinv;
  100.       for (ll=1; ll<=n; ll++)
  101.          {
  102.          /* Next we reduce the rows */
  103.          if (ll != icol)
  104.             {
  105.             /* .. except for the pivot one of course */
  106.             dum = a[ll][icol];
  107.             a[ll][icol] = 0.0;
  108.         if (dum != 0.0)
  109.            {
  110.                for (l=1; l<=n; l++) a[ll][l] -= a[icol][l]*dum;
  111.                b[ll] -= b[icol]*dum;
  112.            }
  113.             }
  114.          }
  115.       }
  116.  
  117.    /* This is the end of the main loop over columns of the reduction.
  118.     * It only remains to unscramble the solution in view of the column
  119.     * interchanges.  We do this by interchanging pairs of columns in the 
  120.     * inverse order that the permutation was built up */
  121.  
  122.    for (l=n; l>=1; l--)
  123.       {
  124.       if (indxr[l] != indxc[l])
  125.      for (k=1; k<=n; k++)
  126.         SWAP(a[k][indxr[l]], a[k][indxc[l]]);
  127.       }
  128.  
  129.    /* Free temporaries */
  130.    FREE (indxc, 1);
  131.    FREE (indxr, 1);
  132.    FREE (ipiv, 1);
  133.  
  134.    return(TRUE);
  135.    }
  136.  
  137.