home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS - Coast to Coast / simteldosarchivecoasttocoast2.iso / c / lpc05b.zip / MATSOLVE.C < prev    next >
C/C++ Source or Header  |  1992-05-22  |  4KB  |  175 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 (attention! see commen)
  33. *        ri = place to put row index (size of n)
  34. *    retrn:    A = successful, NULL = fail
  35. *    comen:    A will be overwritten to be a LU-composite matrix
  36. *    note:    the LU decomposed may NOT be equal to the LU of
  37. *        the orignal matrix a. But equal to the LU of the
  38. *        rows interchanged matrix.
  39. *-----------------------------------------------------------------------------
  40. */
  41. MATRIX mat_lu( A, ri )
  42. MATRIX A;
  43. int *ri;
  44. {
  45.     int    i, j, k, n;
  46.     int    maxi, tmp;
  47.     double    c, c1;
  48.  
  49.     n = MatCol(A);
  50.  
  51.     for (i=0; i<n; i++)
  52.         {
  53.         ri[i] = i;
  54.         }
  55.  
  56.     for (k=0; k<n; k++)
  57.     {
  58.     /*
  59.     * --- partial pivoting ---
  60.     */
  61.     for (i=k, maxi=k, c=0.0; i<n; i++)
  62.         {
  63.         c1 = fabs( A[ri[i]][k] );
  64.         if (c1 > c)
  65.             {
  66.             c = c1;
  67.             maxi = i;
  68.             }
  69.         }
  70.     tmp = ri[k];
  71.     ri[k] = ri[maxi];
  72.     ri[maxi] = tmp;
  73.  
  74.     /*
  75.     *    suspected singular matrix
  76.     */
  77.     if ( A[ri[k]][k] == 0.0 )
  78.         return ((MATRIX)mat_error(MAT_SINGULAR));
  79.  
  80.     for (i=k+1; i<n; i++)
  81.         {
  82.         /*
  83.         * --- calculate m(i,j) ---
  84.         */
  85.         A[ri[i]][k] = A[ri[i]][k] / A[ri[k]][k];
  86.  
  87.         /*
  88.         * --- elimination ---
  89.         */
  90.         for (j=k+1; j<n; j++)
  91.             {
  92.             A[ri[i]][j] -= A[ri[i]][k] * A[ri[k]][j];
  93.             }
  94.         }
  95.     }
  96.  
  97. }
  98.  
  99. /*
  100. *-----------------------------------------------------------------------------
  101. *    funct:    mat_backsubs1
  102. *    desct:    back substitution
  103. *    given:    A = square matrix A (LU composite)
  104. *        !! B = column matrix B (attention!, see comen)
  105. *        !! X = place to put the result of X
  106. *        xcol = column of x to put the result
  107. *        ri = row index (in case after partial pivoting)
  108. *    retrn:    column matrix X (of AX = B)
  109. *    comen:    B will be overwritten
  110. *-----------------------------------------------------------------------------
  111. */
  112. MATRIX mat_backsubs1( A, B, X, xcol, ri )
  113. MATRIX A, B, X;
  114. int xcol;
  115. int *ri;
  116. {
  117.     int    i, j, k, n;
  118.     double    sum;
  119.  
  120.     n = MatCol(A);
  121.  
  122.     for (k=0; k<n; k++)
  123.         {
  124.         for (i=k+1; i<n; i++)
  125.             B[ri[i]][0] -= A[ri[i]][k] * B[ri[k]][0];
  126.         }
  127.  
  128.     X[n-1][xcol] = B[ri[n-1]][0] / A[ri[n-1]][n-1];
  129.     for (k=n-2; k>=0; k--)
  130.         {
  131.         sum = 0.0;
  132.         for (j=k+1; j<n; j++)
  133.             {
  134.             sum += A[ri[k]][j] * X[j][xcol];
  135.             }
  136.         X[k][xcol] = (B[ri[k]][0] - sum) / A[ri[k]][k];
  137.         }
  138.  
  139.     return (X);
  140. }
  141.  
  142. /*
  143. *-----------------------------------------------------------------------------
  144. *    funct:    mat_lsolve
  145. *    desct:    solve linear equations
  146. *    given:    a = square matrix A
  147. *        b = column matrix B
  148. *    retrn:    column matrix X (of AX = B)
  149. *-----------------------------------------------------------------------------
  150. */
  151. MATRIX mat_lsolve( a, b )
  152. MATRIX a, b;
  153. {
  154.     MATRIX    A, B, X;
  155.     int    *ri, i, n;
  156.     double    temp;
  157.  
  158.     n = MatCol(a);
  159.     A = mat_copy(a);
  160.     B = mat_copy(b);
  161.     X = mat_creat(n, 1, ZERO_MATRIX);
  162.  
  163.     if ((ri = (int *)malloc(sizeof(int) * n)) == NULL)
  164.         return (NULL);
  165.  
  166.     mat_lu( A, ri );
  167.     mat_backsubs1( A, B, X, 0, ri );
  168.  
  169.     free(ri);
  170.     mat_free(A);
  171.     mat_free(B);
  172.  
  173.     return (X);
  174. }
  175.