Subversion Repositories gelsvn

Rev

Rev 54 | Rev 116 | Go to most recent revision | Only display areas with differences | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

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