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