home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 7 / 07.iso / c / c129 / 1.ddi / GJ.C < prev    next >
Encoding:
C/C++ Source or Header  |  1989-09-15  |  2.6 KB  |  97 lines
  1. # include <math.h>
  2.  
  3. void gjswap(float *s1,float *s2)
  4. {
  5.    float temp;
  6.  
  7.    temp = (*s1);
  8.    (*s1) = (*s2);
  9.    (*s2) = temp;
  10. }
  11.  
  12. void GaussJordan( float *coefary,float *constary,
  13.                   int numcol,float *solcoef, float *invary,
  14.                   float *det)
  15. {
  16.    int pivlst[50][2];
  17.    char pivchk[50];
  18.    int i;
  19.    int j;
  20.    int k;
  21.    int l;
  22.    int lerow;
  23.    int lecol;
  24.    int l1;
  25.    float piv;
  26.    float t;
  27.    float leval;
  28.  
  29.    (*det) = 1;
  30.    for ( i = 0; i <= numcol-1; ++i ) {
  31.       pivchk[i] = 0;
  32.       for ( j = 0; j <= numcol-1; ++j ) {
  33.          invary[i*numcol+j] = coefary[i*numcol+j];
  34.                }
  35.    }
  36.    for ( i = 0; i <= numcol-1; ++i ) {
  37.       leval = 0.0;
  38.       for ( j = 0; j <= numcol-1; ++j ) {
  39.          if ( ! (pivchk[j]) ) {
  40.             for ( k = 0; k <= numcol-1; ++k ) {
  41.                if ( ! (pivchk[k]) ) {
  42.                   if ( fabs(invary[j*numcol+k]) > leval ) {
  43.                      lerow = j;
  44.                      lecol = k;
  45.                      leval = fabs(invary[j*numcol+k]);
  46.                   }
  47.                }
  48.             }
  49.          }
  50.       }
  51.       pivchk[lecol] = 1;
  52.       pivlst[i][0] = lerow;
  53.       pivlst[i][1] = lecol;
  54.       if ( lerow != lecol ) {
  55.          (*det) = -(*det);
  56.          for ( l = 0; l <= numcol-1; ++l ) {
  57.             gjswap(&invary[lerow*numcol+l],&invary[lecol*numcol+l]);
  58.          }
  59.          gjswap(&constary[lerow],&constary[lecol]);
  60.       }
  61.       piv = invary[lecol*numcol+lecol];
  62.       (*det) = (*det) * piv;
  63.       if ( (*det) > 1.0e+20 ) {
  64.          (*det) = 1.0e+30;
  65.       }
  66.       invary[lecol*numcol+lecol] = 1.0;
  67.       for ( l = 0; l <= numcol-1; ++l ) {
  68.          invary[lecol*numcol+l] = invary[lecol*numcol+l] / piv;
  69.       }
  70.       constary[lecol] = constary[lecol] / piv;
  71.       for ( l1 = 0; l1 <= numcol-1; ++l1 ) {
  72.          if ( l1 != lecol ) {
  73.             t = invary[l1*numcol+lecol];
  74.             invary[l1*numcol+lecol] = 0.0;
  75.             for ( l = 0; l <= numcol-1; ++l ) {
  76.                invary[l1*numcol+l] = invary[l1*numcol+l] -
  77.                                      invary[lecol*numcol+l] * t;
  78.             }
  79.             constary[l1] = constary[l1] - constary[lecol] * t;
  80.          }
  81.       }
  82.    }
  83.   for ( i = 0; i <= numcol-1; ++i ) {
  84.       l = numcol - i - 1;
  85.       if ( pivlst[l][0] != pivlst[l][1] ) {
  86.          lerow = pivlst[l][0];
  87.          lecol = pivlst[l][1];
  88.          for ( k = 0; k <= numcol-1; ++k ) {
  89.             gjswap(&invary[k*numcol+lerow],&invary[k*numcol+lecol]);
  90.          }
  91.       }
  92.    }
  93.    for ( i = 0; i <= numcol-1; ++i ) {
  94.       solcoef[i] = constary[i];
  95.    }
  96. }
  97.