Subversion Repositories gelsvn

Rev

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