Subversion Repositories gelsvn

Rev

Rev 2 | Rev 54 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
2 bj 1 
#include "eigensolution.h"
2 
#include "Mat2x2f.h"
3 
#include "Mat3x3f.h"
4 
#include "Mat4x4f.h"
5 
 
6 
#include <iostream>
7 
 
8 
using namespace std;
9 
 
10 
 
11 
namespace
12 
{
12 jab 13 
  // During experiments 925 iterations were observed for a hard problem
14 
  // Ten times that should be safe.
15 
  const int KMAX = 10000;
2 bj 16 
 
12 jab 17 
  // The threshold below is the smallest that seems to give reliable
18 
  // solutions.
19 
  const float EV_THRESH = 0.0000001;
2 bj 20 
}
21 
 
22 
namespace CGLA
23 
{
24 
 
12 jab 25 
  template <class MT>
26 
  int power_eigensolution(const MT& Ap, MT& Q, MT& L, int max_sol)
27 
  {
28 
    L = MT(0);
2 bj 29 
 
12 jab 30 
    typedef typename MT::VectorType VT;
31 
    MT A = Ap;
32 
    int n = s_min(MT::get_v_dim(), max_sol);
2 bj 33 
 
12 jab 34 
    for(int i=0;i<n;++i)
35 
      {
36 
	// Seed the eigenvector estimate
37 
	VT q(1);
38 
	float l=0,l_old;
2 bj 39 
 
12 jab 40 
	// As long as we haven't reached the max iterations and the
41 
	// eigenvalue has not converged, do
42 
	int k=0;
43 
	for(; k<2 || k<KMAX && (l-l_old > EV_THRESH * l) ; ++k)
44 
	  {
45 
	    // Multiply the eigenvector estimate onto A
46 
	    const VT z = A * q;
2 bj 47 
 
12 jab 48 
	    // Check that we did not get the null vector.
49 
	    float z_len = z.length();
50 
	    if(z_len < EV_THRESH)
51 
	      return i;
2 bj 52 
 
12 jab 53 
	    // Normalize to get the new eigenvector
54 
	    q = z/z_len;
55 
 
56 
	    // Record the old eigenvalue estimate and get a new estimate.
57 
	    l_old = l;
58 
	    l = dot(q, A * q);
59 
	  }
60 
	// If we hit the max iterations, we also don't trust the eigensolution
61 
	if(k==KMAX)
62 
	  return i;
2 bj 63 
 
12 jab 64 
	// Update the solution by adding the eigenvector to Q and
65 
	// the eigenvalue to the diagonal of L.
66 
	Q[i] = q;
67 
	L[i][i] = l;
2 bj 68 
 
12 jab 69 
	// Update A by subtracting the subspace represented by the
70 
	// eigensolution just found. This is called the method of
71 
	// deflation.
72 
	MT B;
73 
	outer_product(q,q,B);
74 
	A = A - l * B;
75 
      }
76 
    return n;
77 
  }
2 bj 78 
 
12 jab 79 
  /* There is no reason to put this template in a header file, since
80 
     we will only use it on matrices defined in CGLA. Instead, we
81 
     explicitly instantiate the function for the square matrices
82 
     of CGLA */
83 
  template int power_eigensolution<Mat2x2f>(const Mat2x2f&,
84 
					    Mat2x2f&,Mat2x2f&,int);
2 bj 85 
 
12 jab 86 
  template int power_eigensolution<Mat3x3f>(const Mat3x3f&,
87 
					    Mat3x3f&,Mat3x3f&,int);
88 
  template int power_eigensolution<Mat4x4f>(const Mat4x4f&,
89 
					    Mat4x4f&,Mat4x4f&,int);
2 bj 90 
}
91