Subversion Repositories gelsvn

Rev

Rev 12 | Rev 54 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

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