home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 41 / IOPROG_41.ISO / soft / c++ / NUMCPP11.ZIP / linalg.hpp < prev    next >
Encoding:
C/C++ Source or Header  |  1999-05-17  |  8.1 KB  |  297 lines

  1. //===================================================================
  2. // linalg.h
  3. //
  4. // Version 1.1
  5. //
  6. // Written by:
  7. //   Brent Worden
  8. //   WordenWare
  9. //   email:  Brent@Worden.org
  10. //
  11. // Copyright (c) 1998-1999 WordenWare
  12. //
  13. // Created:  August 28, 1998
  14. // Revised:  April 17, 1999
  15. //===================================================================
  16.  
  17. #ifndef _LINALG_H_
  18. #define _LINALG_H_
  19.  
  20. #include "numerror.h"
  21. #include "matrix.hpp"
  22.  
  23. NUM_BEGIN
  24.  
  25. template<class T>
  26. bool choldcmp(Matrix<T>& a, Vector<T>& p)
  27. //-------------------------------------------------------------------
  28. // Given a positive-definite symmetric matrix a, this routine
  29. // constructs its Colesky decomposition, a = LL'.  On input, only the
  30. // upper triangle of a need be given; it is not modified.  The
  31. // Cholesky factor L is returned in the lower traingle of a, except
  32. // for its diagonal elements which are returned in p.  The routine
  33. // returns true if the inversion was successful, otherwise it returns
  34. // false.
  35. //-------------------------------------------------------------------
  36. {
  37.     int i, j, k;
  38.     T sum;
  39.     int n = a.rows();
  40.  
  41.     p[0] = sqrt(a[0][0]);
  42.     for(i = 1; i < n; i++){
  43.         a[i][0] /= p[0];
  44.     }
  45.     for(i = 1; i < n; ++i){
  46.         sum = a[i][i];
  47.         for(k = 0; k < i; ++k){
  48.             sum -= a[i][k] * a[i][k];
  49.         }
  50.         if(sum <= 0.0){
  51.             return false;
  52.         }
  53.         p[i] = sqrt(sum);
  54.         for(j = n - 1; j > i; --j){
  55.             sum = a[j][i];
  56.             for(k = 0; k < i; ++k){
  57.                 sum -= a[j][k] * a[i][k];
  58.             }
  59.             a[j][i] = sum / p[i];
  60.         }
  61.     }
  62.     
  63.     return true;
  64. };
  65.  
  66. template<class T>
  67. void cholsl(Matrix<T>& a, Vector<T>& p, Vector<T>& x, Vector<T>& b)
  68. //-------------------------------------------------------------------
  69. // Solves the set of n linear equations a.x = b, where a is a
  70. // positive-definite symmetric matrix.  a and p are input as the
  71. // output of the routine choldcmp.  Only the lower triangle of a is
  72. // accessed.  b is input as the right-hand side vector.  The solution
  73. // vector is returned in x[0..n-1].  a, and p are not modified and
  74. // can be left in place for successive calls with different
  75. // right-hand sides b.  b is not modified unless you identify b as x
  76. // in the calling sequence, which is allowed.
  77. //-------------------------------------------------------------------
  78. {
  79.     int i, k;
  80.     T sum;
  81.     int n = a.rows();
  82.  
  83.     for(i = 0; i < n; i++){
  84.         for(sum = b[i], k = i - 1; k >= 0; --k){
  85.             sum -= a[i][k] * x[k];
  86.         }
  87.         x[i] = sum / p[i];
  88.     }
  89.     for(i = n - 1; i >= 0; --i){
  90.         for(sum = x[i], k = i + 1; k < n; ++k){
  91.             sum -= a[k][i] * x[k];
  92.         }
  93.         x[i] = sum / p[i];
  94.     }
  95. }
  96.  
  97. template<class T>
  98. T determinant(Matrix<T>& a)
  99. //-------------------------------------------------------------------
  100. // Returns the determinant of the matrix a.
  101. //-------------------------------------------------------------------
  102. {
  103.     Matrix<T> cp(a);
  104.     T ret;
  105.     int i;
  106.     int n = cp.rows();
  107.     Vector<int> indx(n);
  108.     int sgn;
  109.  
  110.     if(ludcmp(cp, indx, sgn)){
  111.         ret = T(sgn);
  112.         for(i = 0; i < n; i++){
  113.             ret *= cp[i][i];
  114.         }
  115.     } else {
  116.         ret = 0.0;
  117.     }
  118.     return ret;
  119. };
  120.  
  121. template<class T>
  122. bool inverse(Matrix<T>& a)
  123. //-------------------------------------------------------------------
  124. // Replaces the matrix a with its inverse.  The routine returns true
  125. // if the inversion was successful, otherwise it returns false.
  126. //-------------------------------------------------------------------
  127. {
  128.     int i, j;
  129.     bool ret = false;
  130.     int n = a.rows();
  131.     Matrix<T> ai(n, n);
  132.     Vector<int> indx(n);
  133.     Vector<T> col(n);
  134.     int d;
  135.  
  136.     if(ludcmp(a, indx, d)){
  137.         for(j = 0; j < n; j++){
  138.             for(i = 0; i < n; i++){
  139.                 col[i] = T(0);
  140.             }
  141.             col[j] = T(1);
  142.             lubksb(a, indx, col);
  143.             for(i = 0; i < n; i++){
  144.                 ai[i][j] = col[i];
  145.             }
  146.         }
  147.         a = ai;
  148.         ret = true;
  149.     }
  150.     return ret;
  151. };
  152.  
  153. template<class T>
  154. void lubksb(Matrix<T>& a, Vector<int>& indx, Vector<T>& b)
  155. //-------------------------------------------------------------------
  156. // Given a matrix a and permutation vector indx returned from ludcmp,
  157. // this routines solves the set of linear equations a.x = b.  On 
  158. // input b holds the right-hand side vector.  On output, it holds the
  159. // solution vector x.  a and indx are not modified by this routine
  160. // and can be left in place for successive colls with different
  161. // right-hand sides b.
  162. //-------------------------------------------------------------------
  163. {
  164.     int i, ii = -1, ip, j;
  165.     T sum;
  166.     int n = a.rows();
  167.  
  168.     for(i = 0; i < n; i++){
  169.         ip = indx[i];
  170.         sum = b[ip];
  171.         b[ip] = b[i];
  172.         if(ii > -1){
  173.             for(j = ii; j <= i - 1; j++){
  174.                 sum -= a[i][j] * b[j];
  175.             }
  176.         } else if(sum){
  177.             ii = i;
  178.         }
  179.         b[i] = sum;
  180.     }
  181.     for(i = n - 1; i >= 0; i--){
  182.         sum = b[i];
  183.         for(j = i + 1; j < n; j++){
  184.             sum -= a[i][j] * b[j];
  185.         }
  186.         b[i] = sum / a[i][i];
  187.     }
  188. };
  189.  
  190. template<class T>
  191. bool ludcmp(Matrix<T>& a, Vector<int>& indx, int& d)
  192. //-------------------------------------------------------------------
  193. // See above comment.  d is output as 1 or -1 depending on whether
  194. // the number of row interchanges was even or odd, respectively.
  195. // This routine is used in combination with lubksb to solve linear
  196. // equations or invert a matrix.
  197. //-------------------------------------------------------------------
  198. {
  199.     int i, imax, j, k;
  200.     T big, dum, sum, temp;
  201.     int n = a.rows();
  202.     Vector<T> vv(n);
  203.     T zero(0);
  204.     T one(1);
  205.  
  206.     d = 1;
  207.     for(i = 0; i < n; i++){
  208.         big = zero;
  209.         for(j = 0; j < n; j++){
  210.             if((temp = absoluteValue(a[i][j])) > big){
  211.                 big = temp;
  212.             }
  213.         }
  214.         if(big == zero){
  215.             throw Exception("num::ludcmp", "Singular Matrix");
  216.             return false;
  217.         }
  218.         vv[i] = one / big;
  219.     }
  220.     for(j = 0; j < n; j++){
  221.         for(i = 0; i < j; i++){
  222.             sum = a[i][j];
  223.             for(k = 0; k < i; k++){
  224.                 sum -= a[i][k] * a[k][j];
  225.             }
  226.             a[i][j] = sum;
  227.         }
  228.         big = zero;
  229.         for(i = j; i < n; i++){
  230.             sum = a[i][j];
  231.             for(k = 0; k < j; k++){
  232.                 sum -= a[i][k] * a[k][j];
  233.             }
  234.             a[i][j] = sum;
  235.             if((dum = vv[i] * absoluteValue(sum)) >= big){
  236.                 big = dum;
  237.                 imax = i;
  238.             }
  239.         }
  240.         if(j != imax){
  241.             for(k = 0; k < n; k++){
  242.                 dum = a[imax][k];
  243.                 a[imax][k] = a[j][k];
  244.                 a[j][k] = dum;
  245.             }
  246.             d = -d;
  247.             vv[imax] = vv[j];
  248.         }
  249.         indx[j] = imax;
  250.         if(a[j][j] == 0.0){
  251.             a[j][j] = NUMERICS_FLOAT_MIN;
  252.         }
  253.         if(j != n - 1){
  254.             dum = one / (a[j][j]);
  255.             for(i = j + 1; i < n; i++){
  256.                 a[i][j] *= dum;
  257.             }
  258.         }
  259.     }
  260.     return true;
  261. };
  262.  
  263. template<class T>
  264. bool linsolve(Matrix<T>& a, Vector<T>& b, Vector<T>& x)
  265. //-------------------------------------------------------------------
  266. // Solves the set of linear equations a.x = b.  The routine returns
  267. // true if the solution vector is successfully found, otherwise it
  268. // returns false.
  269. //-------------------------------------------------------------------
  270. {
  271.     Matrix<T> ac(a);
  272.     Vector<T> bc(b);
  273.     bool ret = false;
  274.     int d;
  275.     Vector<int> indx(a.rows());
  276.  
  277.     if(ludcmp(ac, indx, d)){
  278.         lubksb(ac, indx, bc);
  279.         x = bc;
  280.         ret = true;
  281.     }
  282.     return ret;
  283. };
  284.  
  285. NUM_END
  286.  
  287. #endif
  288.  
  289. //===================================================================
  290. // Revision History
  291. //
  292. // Version 1.0 - 08/28/1998 - New.
  293. // Version 1.1 - 04/10/1999 - Added Numerics namespace.
  294. //               04/17/1999 - Use matrix and vector classes.
  295. //               04/17/1999 - Removed dot product functions.
  296. //===================================================================
  297.