home *** CD-ROM | disk | FTP | other *** search
/ The C Users' Group Library 1994 August / wc-cdrom-cusersgrouplibrary-1994-08.iso / vol_300 / 344_01 / mtxsolv.c < prev    next >
Text File  |  1991-05-29  |  2KB  |  81 lines
  1. /*-----------------------------------------------------------------------
  2.    Main File : mtx.exe
  3.    File Name : mtxsolv.c
  4.  
  5.    Purpose - G-J computations on an A|b matrix
  6. -------------------------------------------------------------------------*/
  7. #include "mtx.h"
  8. #include "mtxcle.h"
  9.  
  10.  
  11.   /* Indexing macro for temporary array */
  12. #define Pivot(i)            *(pivot+(i))
  13.  
  14. extern int MSize;
  15.  
  16.   /* Head pointers */
  17. Mtype *Aptr=NULL;
  18. Mtype *bptr=NULL;
  19.  
  20.  
  21. void dswap(Mtype *p1, Mtype *p2)
  22. {
  23.   Mtype temp;
  24.  
  25.   temp=*p1;
  26.   *p1=*p2;
  27.   *p2=temp;
  28. }
  29.  
  30.  
  31. void solve_matrix()
  32. {
  33.   int i,j,k;
  34.   const long int n = MSize;
  35.   Mtype *pivot,mabs,mmax,prod,sum;
  36.  
  37.   if((pivot=(Mtype *)alloc(n,sizeof(Mtype)))==NULL)
  38.     error(errHISIZE);
  39.  
  40.   for(k=0; k < n-1; ++k)
  41.   {
  42.     Pivot(k) = k;
  43.     mmax = (Mtype)absMtype(A(k,k));
  44.     for(i=k+1; i < n; ++i)
  45.     {
  46.       mabs = (Mtype)absMtype(A(i,k));
  47.       if(mabs > mmax)
  48.       {
  49.         Pivot(k) = i;
  50.         mmax = mabs;
  51.       }
  52.     }
  53.     if(mmax == 0)
  54.       error(errSINGULARITY);
  55.     if(Pivot(k) != k)
  56.     {
  57.       i = Pivot(k);
  58.       dswap(bptr+k, bptr+i);
  59.       for(j=k; j < n ; ++j)
  60.         dswap(Aptr+(k*n)+j, Aptr+(i*n)+j);
  61.     }
  62.     for(i=k+1; i < n; i++)
  63.     {
  64.       prod = A(i,k) / A(k,k);
  65.       A(i,k) = prod;
  66.       b(i) -= prod * b(k);
  67.       for(j=k+1; j < n; ++j)
  68.         A(i,j) -= prod * A(k,j);
  69.     }
  70.   }
  71.   if(A(n-1,n-1) == 0)
  72.     error(errSINGULARITY);
  73.   for(i=n-1; i >= 0; --i)
  74.   {
  75.     for(sum=0.0, j=i+1; j < n; ++j)
  76.       sum += A(i,j) * b(j);
  77.     b(i) = (b(i)-sum) / A(i,i);
  78.   }
  79. }
  80.  
  81.