home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / matrix04 / matsolve.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-05-24  |  3.8 KB  |  183 lines
  1. /*
  2. *-----------------------------------------------------------------------------
  3. *    file:    matsolve.c
  4. *    desc:    solve linear equations
  5. *    by:    ko shu pui, patrick
  6. *    date:    24 nov 91 v0.1
  7. *    revi:    14 may 92 v0.2
  8. *    ref:
  9. *       [1] Mary L.Boas, "Mathematical Methods in the Physical Sciene,"
  10. *    John Wiley & Sons, 2nd Ed., 1983. Chap 3.
  11. *
  12. *    [2] Kendall E.Atkinson, "An Introduction to Numberical Analysis,"
  13. *    John Wiley & Sons, 1978.
  14. *
  15. *-----------------------------------------------------------------------------
  16. */
  17. #include <stdio.h>
  18. #include <math.h>
  19.  
  20. #ifdef    __TURBOC__
  21. #include <alloc.h>
  22. #else
  23. #include <malloc.h>
  24. #endif
  25.  
  26. #include "matrix.h"
  27.  
  28. /*
  29. *-----------------------------------------------------------------------------
  30. *    funct:    mat_lu
  31. *    desct:    in-place LU decomposition with partial pivoting
  32. *    given:    !! A = square matrix (n x n) !ATTENTION! see commen
  33. *        P = permutation vector (n x 1)
  34. *    retrn:    number of permutation performed
  35. *        -1 means suspected singular matrix
  36. *    comen:    A will be overwritten to be a LU-composite matrix
  37. *
  38. *    note:    the LU decomposed may NOT be equal to the LU of
  39. *        the orignal matrix a. But equal to the LU of the
  40. *        rows interchanged matrix.
  41. *-----------------------------------------------------------------------------
  42. */
  43. int mat_lu( A, P )
  44. MATRIX A;
  45. MATRIX P;
  46. {
  47.     int    i, j, k, n;
  48.     int    maxi, tmp;
  49.     double    c, c1;
  50.     int    p;
  51.  
  52.     n = MatCol(A);
  53.  
  54.     for (p=0,i=0; i<n; i++)
  55.         {
  56.         P[i][0] = i;
  57.         }
  58.  
  59.     for (k=0; k<n; k++)
  60.     {
  61.     /*
  62.     * --- partial pivoting ---
  63.     */
  64.     for (i=k, maxi=k, c=0.0; i<n; i++)
  65.         {
  66.         c1 = fabs( A[P[i][0]][k] );
  67.         if (c1 > c)
  68.             {
  69.             c = c1;
  70.             maxi = i;
  71.             }
  72.         }
  73.  
  74.     /*
  75.     *    row exchange, update permutation vector
  76.     */
  77.     if (k != maxi)
  78.         {
  79.         p++;
  80.         tmp = P[k][0];
  81.         P[k][0] = P[maxi][0];
  82.         P[maxi][0] = tmp;
  83.         }
  84.  
  85.     /*
  86.     *    suspected singular matrix
  87.     */
  88.     if ( A[P[k][0]][k] == 0.0 )
  89.         return (-1);
  90.  
  91.     for (i=k+1; i<n; i++)
  92.         {
  93.         /*
  94.         * --- calculate m(i,j) ---
  95.         */
  96.         A[P[i][0]][k] = A[P[i][0]][k] / A[P[k][0]][k];
  97.  
  98.         /*
  99.         * --- elimination ---
  100.         */
  101.         for (j=k+1; j<n; j++)
  102.             {
  103.             A[P[i][0]][j] -= A[P[i][0]][k] * A[P[k][0]][j];
  104.             }
  105.         }
  106.     }
  107.  
  108.     return (p);
  109. }
  110.  
  111. /*
  112. *-----------------------------------------------------------------------------
  113. *    funct:    mat_backsubs1
  114. *    desct:    back substitution
  115. *    given:    A = square matrix A (LU composite)
  116. *        !! B = column matrix B (attention!, see comen)
  117. *        !! X = place to put the result of X
  118. *        P = Permutation vector (after calling mat_lu)
  119. *        xcol = column of x to put the result
  120. *    retrn:    column matrix X (of AX = B)
  121. *    comen:    B will be overwritten
  122. *-----------------------------------------------------------------------------
  123. */
  124. MATRIX mat_backsubs1( A, B, X, P, xcol )
  125. MATRIX A, B, X, P;
  126. int xcol;
  127. {
  128.     int    i, j, k, n;
  129.     double    sum;
  130.  
  131.     n = MatCol(A);
  132.  
  133.     for (k=0; k<n; k++)
  134.         {
  135.         for (i=k+1; i<n; i++)
  136.             B[P[i][0]][0] -= A[P[i][0]][k] * B[P[k][0]][0];
  137.         }
  138.  
  139.     X[n-1][xcol] = B[P[n-1][0]][0] / A[P[n-1][0]][n-1];
  140.     for (k=n-2; k>=0; k--)
  141.         {
  142.         sum = 0.0;
  143.         for (j=k+1; j<n; j++)
  144.             {
  145.             sum += A[P[k][0]][j] * X[j][xcol];
  146.             }
  147.         X[k][xcol] = (B[P[k][0]][0] - sum) / A[P[k][0]][k];
  148.         }
  149.  
  150.     return (X);
  151. }
  152.  
  153. /*
  154. *-----------------------------------------------------------------------------
  155. *    funct:    mat_lsolve
  156. *    desct:    solve linear equations
  157. *    given:    a = square matrix A
  158. *        b = column matrix B
  159. *    retrn:    column matrix X (of AX = B)
  160. *-----------------------------------------------------------------------------
  161. */
  162. MATRIX mat_lsolve( a, b )
  163. MATRIX a, b;
  164. {
  165.     MATRIX    A, B, X, P;
  166.     int    i, n;
  167.     double    temp;
  168.  
  169.     n = MatCol(a);
  170.     A = mat_copy(a);
  171.     B = mat_copy(b);
  172.     X = mat_creat(n, 1, ZERO_MATRIX);
  173.     P = mat_creat(n, 1, UNDEFINED);
  174.  
  175.     mat_lu( A, P );
  176.     mat_backsubs1( A, B, X, P, 0 );
  177.  
  178.     mat_free(A);
  179.     mat_free(B);
  180.     mat_free(P);
  181.     return (X);
  182. }
  183.