home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / remes / simq.c < prev   
Encoding:
Text File  |  1992-11-17  |  2.9 KB  |  178 lines
  1. /*                            simq.c
  2.  *
  3.  *    Solution of simultaneous linear equations AX = B
  4.  *    by Gaussian elimination with partial pivoting
  5.  *
  6.  *
  7.  *
  8.  * SYNOPSIS:
  9.  *
  10.  * double A[n*n], B[n], X[n];
  11.  * int n, flag;
  12.  * int IPS[];
  13.  * int simq();
  14.  *
  15.  * ercode = simq( A, B, X, n, flag, IPS );
  16.  *
  17.  *
  18.  *
  19.  * DESCRIPTION:
  20.  *
  21.  * B, X, IPS are vectors of length n.
  22.  * A is an n x n matrix (i.e., a vector of length n*n),
  23.  * stored row-wise: that is, A(i,j) = A[ij],
  24.  * where ij = i*n + j, which is the transpose of the normal
  25.  * column-wise storage.
  26.  *
  27.  * The contents of matrix A are destroyed.
  28.  *
  29.  * Set flag=0 to solve.
  30.  * Set flag=-1 to do a new back substitution for different B vector
  31.  * using the same A matrix previously reduced when flag=0.
  32.  *
  33.  * The routine returns nonzero on error; messages are printed.
  34.  *
  35.  *
  36.  * ACCURACY:
  37.  *
  38.  * Depends on the conditioning (range of eigenvalues) of matrix A.
  39.  *
  40.  *
  41.  * REFERENCE:
  42.  *
  43.  * Computer Solution of Linear Algebraic Systems,
  44.  * by George E. Forsythe and Cleve B. Moler; Prentice-Hall, 1967.
  45.  *
  46.  */
  47.  
  48. /*                            simq    2 */
  49.  
  50. int simq( A, B, X, n, flag, IPS )
  51. double A[], B[], X[];
  52. int n, flag;
  53. int IPS[];
  54. {
  55. int i, j, ij, ip, ipj, ipk, ipn;
  56. int idxpiv, iback;
  57. int k, kp, kp1, kpj, kpk, kpn;
  58. int nip, nkp, nm1;
  59. double em, q, rownrm, big, size, pivot, sum;
  60. double fabs();
  61.  
  62. if( flag < 0 )
  63.     goto solve;
  64.  
  65. /*    Initialize IPS and X    */
  66.  
  67. ij=0;
  68. for( i=0; i<n; i++ )
  69.     {
  70.     IPS[i] = i;
  71.     rownrm = 0.0;
  72.     for( j=0; j<n; j++ )
  73.         {
  74.         q = fabs( A[ij] );
  75.         if( rownrm < q )
  76.             rownrm = q;
  77.         ++ij;
  78.         }
  79.     if( rownrm == 0.0 )
  80.         {
  81.         puts("SIMQ ROWNRM=0");
  82.         return(1);
  83.         }
  84.     X[i] = 1.0/rownrm;
  85.     }
  86.  
  87. /*                            simq    3 */
  88. /*    Gaussian elimination with partial pivoting     */
  89.  
  90. nm1 = n-1;
  91. for( k=0; k<nm1; k++ )
  92.     {
  93.     big= 0.0;
  94.     for( i=k; i<n; i++ )
  95.         {
  96.         ip = IPS[i];
  97.         ipk = n*ip + k;
  98.         size = fabs( A[ipk] ) * X[ip];
  99.         if( size > big )
  100.             {
  101.             big = size;
  102.             idxpiv = i;
  103.             }
  104.         }
  105.  
  106.     if( big == 0.0 )
  107.         {
  108.         puts( "SIMQ BIG=0" );
  109.         return(2);
  110.         }
  111.     if( idxpiv != k )
  112.         {
  113.         j = IPS[k];
  114.         IPS[k] = IPS[idxpiv];
  115.         IPS[idxpiv] = j;
  116.         }
  117.     kp = IPS[k];
  118.     kpk = n*kp + k;
  119.     pivot = A[kpk];
  120.     kp1 = k+1;
  121.     for( i=kp1; i<n; i++ )
  122.         {
  123.         ip = IPS[i];
  124.         ipk = n*ip + k;
  125.         em = -A[ipk]/pivot;
  126.         A[ipk] = -em;
  127.         nip = n*ip;
  128.         nkp = n*kp;
  129.         for( j=kp1; j<n; j++ )
  130.             {
  131.             ipj = nip + j;
  132.             A[ipj] = A[ipj] + em * A[nkp + j];
  133.             }
  134.         }
  135.     }
  136. kpn = n * IPS[n-1] + n - 1;    /* last element of IPS[n] th row */
  137. if( A[kpn] == 0.0 )
  138.     {
  139.     puts( "SIMQ A[kpn]=0");
  140.     return(3);
  141.     }
  142.  
  143. /*                            simq 4 */
  144. /*    back substitution    */
  145.  
  146. solve:
  147. ip = IPS[0];
  148. X[0] = B[ip];
  149. for( i=1; i<n; i++ )
  150.     {
  151.     ip = IPS[i];
  152.     ipj = n * ip;
  153.     sum = 0.0;
  154.     for( j=0; j<i; j++ )
  155.         {
  156.         sum += A[ipj] * X[j];
  157.         ++ipj;
  158.         }
  159.     X[i] = B[ip] - sum;
  160.     }
  161.  
  162. ipn = n * IPS[n-1] + n - 1;
  163. X[n-1] = X[n-1]/A[ipn];
  164.  
  165. for( iback=1; iback<n; iback++ )
  166.     {
  167. /* i goes (n-1),...,1    */
  168.     i = nm1 - iback;
  169.     ip = IPS[i];
  170.     nip = n*ip;
  171.     sum = 0.0;
  172.     for( j=i+1; j<n; j++ )
  173.         sum += A[nip+j] * X[j];
  174.     X[i] = (X[i] - sum)/A[nip+i];
  175.     }
  176. return(0);
  177. }
  178.