Subversion Repositories gelsvn

Rev

Rev 116 | Rev 354 | 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.
342 jab 19 
		const unsigned int KMAX = 1e6;
2 bj 20 
 
45 jab 21 
		// The threshold below is the smallest that seems to give reliable
22 
		// solutions.
342 jab 23 
		const double EV_THRESH = 1e-6;
2 bj 24 
}
25 
 
26 
namespace CGLA
27 
{
45 jab 28 
		template <class MT>
29 
		int power_eigensolution(const MT& Ap, MT& Q, MT& L, unsigned int max_sol)
30 
		{
31 
				L = MT(0);
2 bj 32 
 
45 jab 33 
				typedef typename MT::VectorType VT;
34 
				MT A = Ap;
35 
				unsigned int n = s_min(MT::get_v_dim(), max_sol);
2 bj 36 
 
45 jab 37 
				for(unsigned int i=0;i<n;++i)
38 
				{
39 
						// Seed the eigenvector estimate
40 
						VT q(1);
342 jab 41 
						q.normalize();
42 
						double l,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;
342 jab 47 
						do
45 jab 48 
						{
342 jab 49 
							const VT z = A * q;
50 
							double z_len = length(z);
51 
 
52 
							if(z_len < EV_THRESH) return i;
53 
 
54 
							l_old = l;
55 
							l = dot(q, z)>0 ? z_len : -z_len;
56 
							q = z/z_len;
57 
 
58 
							if(++k==KMAX)
59 
								return i;
45 jab 60 
						}
342 jab 61 
						while((fabs(l-l_old) > fabs(EV_THRESH * l)) || k<2);
2 bj 62 
 
45 jab 63 
						// Update the solution by adding the eigenvector to Q and
64 
						// the eigenvalue to the diagonal of L.
65 
						Q[i] = q;
66 
						L[i][i] = l;
2 bj 67 
 
45 jab 68 
						// Update A by subtracting the subspace represented by the
69 
						// eigensolution just found. This is called the method of
70 
						// deflation.
71 
						MT B;
72 
						outer_product(q,q,B);
73 
						A = A - l * B;
74 
				}
75 
				return n;
76 
		}
2 bj 77 
 
45 jab 78 
		/* There is no reason to put this template in a header file, since
79 
			 we will only use it on matrices defined in CGLA. Instead, we
80 
			 explicitly instantiate the function for the square matrices
81 
			 of CGLA */
82 
		template int power_eigensolution<Mat2x2f>(const Mat2x2f&,
83 
																							Mat2x2f&,Mat2x2f&,unsigned int);
2 bj 84 
 
45 jab 85 
		template int power_eigensolution<Mat3x3f>(const Mat3x3f&,
86 
																							Mat3x3f&,Mat3x3f&,unsigned int);
87 
		template int power_eigensolution<Mat4x4f>(const Mat4x4f&,
88 
																							Mat4x4f&,Mat4x4f&,unsigned int);
54 jab 89 
		template int power_eigensolution<Mat2x2d>(const Mat2x2d&,
90 
																							Mat2x2d&,Mat2x2d&,unsigned int);
91 
 
92 
		template int power_eigensolution<Mat3x3d>(const Mat3x3d&,
93 
																							Mat3x3d&,Mat3x3d&,unsigned int);
94 
		template int power_eigensolution<Mat4x4d>(const Mat4x4d&,
95 
																							Mat4x4d&,Mat4x4d&,unsigned int);
2 bj 96 
}
97