home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 41 / IOPROG_41.ISO / soft / c++ / NUMCPP11.ZIP / eigensys.hpp < prev    next >
Encoding:
C/C++ Source or Header  |  1999-05-15  |  4.7 KB  |  203 lines
  1. //====================================================================
  2. // eigensys.hpp
  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 10, 1999
  15. //====================================================================
  16.  
  17. #ifndef _EIGENSYS_HPP_
  18. #define _EIGENSYS_HPP_
  19.  
  20. #include "mathx.h"
  21. #include "matrix.hpp"
  22. #include "numerics.h"
  23. #include "utility.hpp"
  24.  
  25. NUM_BEGIN
  26.  
  27. template<class T>
  28. void tridred(Matrix<T>& a, Vector<T>& d, Vector<T>& e)
  29. //--------------------------------------------------------------------
  30. // Householder reduction of a real, symmetric matrix a.  On output, a
  31. // is replaced by the orthogonal matrix effecting the transformation.
  32. // d returns the diagonal elements of the tridiagonal matrix, and e
  33. // the off-diagonal elements, with e[0] = 0.
  34. //--------------------------------------------------------------------
  35. {
  36.     int n = a.rows();
  37.     T zero(0);
  38.     T f, g, h;
  39.     int i, j, k;
  40.     T scale, tmp;
  41.     d.resize(n);
  42.     e.resize(n);
  43.  
  44.     for(i = 0; i < n; ++i){
  45.         d[i] = a[i][n - 1];
  46.         a[i][n - 1] = a[i][i];
  47.     }
  48.     for(i = n - 1; i >= 0; --i){
  49.         h = zero;
  50.         scale = zero;
  51.         if(i < 1){
  52.             e[i] = zero;
  53.         } else {
  54.             for(k = 0; k <= i; ++k){
  55.                 scale += absoluteValue(d[k]);
  56.             }
  57.             if(scale == zero){
  58.                 for(j = 0; j < i; ++j){
  59.                     d[j] = a[j][i - 1];
  60.                     a[j][i - 1] = a[j][i];
  61.                     a[j][i] = zero;
  62.                 }
  63.                 e[i] = zero;
  64.             } else {
  65.                 for(k = 0; k < i; ++k){
  66.                     tmp = (d[k] /= scale);
  67.                     h += tmp * tmp;
  68.                 }
  69.                 f = d[i - 1];
  70.                 g = (f >= zero ? -sqrt(h) : sqrt(h));
  71.                 e[i] = scale * g;
  72.                 h -= f * g;
  73.                 d[i - 1] = f - g;
  74.                 if(i != 1){
  75.                     for(j = 0; j < i; ++j){
  76.                         e[j] = zero;
  77.                     }
  78.                     for(j = 0; j < i; ++j){
  79.                         f = d[j];
  80.                         g = e[j] + a[j][j] * f;
  81.                         k = j + 1;
  82.                         if(i > k){
  83.                             for(; k < i; ++k){
  84.                                 g += a[j][k] * d[k];
  85.                                 e[k] += a[j][k] * f;
  86.                             }
  87.                         }
  88.                         e[j] = g;
  89.                     }
  90.                     f = zero;
  91.                     for(j = 0; j < i; ++j){
  92.                         e[j] /= h;
  93.                         f += e[j] * d[j];
  94.                     }
  95.                     h = f / (h + h);
  96.                     for(j = 0; j < i; ++j){
  97.                         e[j] -= h * d[j];
  98.                     }
  99.                     for(j = 0; j < i; ++j){
  100.                         f = d[j];
  101.                         g = e[j];
  102.                         for(k = j; k < i; ++k){
  103.                             a[j][k] -= (f * e[k] + g * d[k]);
  104.                         }
  105.                     }
  106.                 }    
  107.                 for(j = 0; j < i; ++j){
  108.                     f = d[j];
  109.                     d[j] = a[j][i - 1];
  110.                     a[j][i - 1] = a[j][i];
  111.                     a[j][i] = f * scale;
  112.                 }
  113.             }
  114.         }
  115.     }
  116. };
  117.  
  118. template<class T>
  119. void trideig(Vector<T>& d, Vector<T>& e)
  120. //--------------------------------------------------------------------
  121. // Algorithm to determine the eigenvalues of a real, symetric,
  122. // tridiagonal matrix, or of a real, symmetric matrix previously
  123. // reduced by tridred().  On input, d contains the diagonal elements
  124. // of the tridaigonal matrix.  On output, it returns the eigenvalues.
  125. // e inputs the subdiagonal elements of the tridiagonal matrix, with
  126. // e[0] arbitrary.  On output, e is destroyed.
  127. //--------------------------------------------------------------------
  128. {
  129.     int n = d.size();
  130.     T zero(0), one(1), two(2);
  131.         
  132.     T b, c, f, g;
  133.     int i, j, k, m;
  134.     T p, r, s;
  135.     T tmp;
  136.     
  137.     if(n != 1) {
  138.         for(i = 1; i < n; ++i){
  139.             e[i - 1] = e[i];
  140.         }
  141.         e[n - 1] = zero;
  142.         for(k = 0; k < n; ++k){
  143.             j = 0;
  144.             for(m = k + 1; m < n; ++m){
  145.                 tmp = absoluteValue(d[m - 1]) + absoluteValue(d[m]);
  146.                 if(tmp == absoluteValue(e[m - 1]) + tmp){
  147.                     break;
  148.                 }
  149.             }
  150.             p = d[k];
  151.             if(m != k + 1){
  152.                 if(j++ == NUMERICS_ITMAX){
  153.                     throw Exception("tridql", "Routine failed to converge.");
  154.                     return;
  155.                 }
  156.                 g = (d[k + 1] - p) / (e[k] * two);
  157.                 r = pythag(g, one);
  158.                 g = d[m - 1] - p + e[k] / (g + (g >= zero ? r : -r));
  159.                 s = one;
  160.                 c = one;
  161.                 p = zero;
  162.                 for(i = m - 2; i >= 0; --i){
  163.                     tmp = e[i];
  164.                     f = s * tmp;
  165.                     b = c * tmp;
  166.                     r = pythag(f, g);
  167.                     e[i + 1] = r;
  168.                     if(r == zero){
  169.                         d[i + 1] -= p;
  170.                         e[m - 1] = zero;
  171.                         break;
  172.                     }
  173.                     s = f / r;
  174.                     c = g / r;
  175.                     g = d[i + 1] - p;
  176.                     r = (d[i] - g) * s + c * two * b;
  177.                     p = s * r;
  178.                     d[i + 1] = g + p;
  179.                     g = c * r - b;
  180.                 }
  181.                 if(i <= k + 1 || r != zero){
  182.                     d[k] -= p;
  183.                     e[k] = g;
  184.                     e[m - 1] = zero;
  185.                 }
  186.             }
  187.         }
  188.     }
  189. };
  190.  
  191. NUM_END
  192.  
  193. #endif
  194.  
  195. //===================================================================
  196. // Revision History
  197. //
  198. // Version 1.0 - 08/28/1998 - New.
  199. // Version 1.1 - 04/10/1999 - Added Numerics namespace.
  200. //                            Removed calcution of eigenvectors from
  201. //                            tridred and trideig
  202. //===================================================================
  203.