Subversion Repositories gelsvn

Rev

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