Subversion Repositories gelsvn

Rev

Rev 12 | Rev 58 | 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 
{
45 jab 13 
		// During experiments 925 iterations were observed for a hard problem
14 
		// Ten times that should be safe.
15 
		const unsigned int KMAX = 10000;
2 bj 16 
 
45 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 
 
45 jab 25 
		template <class MT>
26 
		int power_eigensolution(const MT& Ap, MT& Q, MT& L, unsigned int max_sol)
27 
		{
28 
				L = MT(0);
2 bj 29 
 
45 jab 30 
				typedef typename MT::VectorType VT;
31 
				MT A = Ap;
32 
				unsigned int n = s_min(MT::get_v_dim(), max_sol);
2 bj 33 
 
45 jab 34 
				for(unsigned 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 
 
45 jab 40 
						// As long as we haven't reached the max iterations and the
41 
						// eigenvalue has not converged, do
42 
						unsigned 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 
 
45 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 
 
45 jab 53 
								// Normalize to get the new eigenvector
54 
								q = z/z_len;
12 jab 55 
 
45 jab 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 
 
45 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 
 
45 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 
 
45 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&,unsigned int);
2 bj 85 
 
45 jab 86 
		template int power_eigensolution<Mat3x3f>(const Mat3x3f&,
87 
																							Mat3x3f&,Mat3x3f&,unsigned int);
88 
		template int power_eigensolution<Mat4x4f>(const Mat4x4f&,
89 
																							Mat4x4f&,Mat4x4f&,unsigned int);
2 bj 90 
}
91