Subversion Repositories gelsvn

Rev

Rev 116 | Rev 351 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 116 Rev 342
Line 14... Line 14...
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 = 1e6;
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 double EV_THRESH = 0.00000001;
 23 
		const double EV_THRESH = 1e-6;
24 
}
 24 
}
25 
 
 25 
 
26 
namespace CGLA
 26 
namespace CGLA
27 
{
 27 
{
28 
 
 - 
 
29 
		template <class MT>
 28 
		template <class MT>
30 
		int power_eigensolution(const MT& Ap, MT& Q, MT& L, unsigned int max_sol)
 29 
		int power_eigensolution(const MT& Ap, MT& Q, MT& L, unsigned int max_sol)
31 
		{
 30 
		{
32 
				L = MT(0);
 31 
				L = MT(0);
33 
 
 32 
 
Line 37... Line 36...
37 
 
 36 
 
38 
				for(unsigned int i=0;i<n;++i)
 37 
				for(unsigned int i=0;i<n;++i)
39 
				{
 38 
				{
40 
						// Seed the eigenvector estimate
 39 
						// Seed the eigenvector estimate
41 
						VT q(1);
 40 
						VT q(1);
- 
 
 41 
						q.normalize();
42 
						double l=0,l_old;
 42 
						double l,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 
						do
48 
						{
 48 
						{
49 
								// Multiply the eigenvector estimate onto A
 - 
 
50 
								const VT z = A * q;
 49 
							const VT z = A * q;
51 
 
 - 
 
52 
								// Check that we did not get the null vector.
 - 
 
53 
								double z_len = z.length();
 50 
							double z_len = length(z);
54 
								if(z_len < EV_THRESH)
 - 
 
55 
										return i;
 51
56 
 
 - 
 
57 
								// Normalize to get the new eigenvector
 52 
							if(z_len < EV_THRESH) return i;
58 
								q = z/z_len;
 53
59 
 
 - 
 
60 
								// Record the old eigenvalue estimate and get a new estimate.
 - 
 
61 
								l_old = l;
 54 
							l_old = l;
62 
								l = dot(q, A * q);
 55 
							l = dot(q, z)>0 ? z_len : -z_len;
- 
 
 56 
							q = z/z_len;
63 
						}
 57
64 
						// If we hit the max iterations, we also don't trust the eigensolution
 - 
 
65 
						if(k==KMAX)
 58 
							if(++k==KMAX)
66 
								return i;
 59 
								return i;
- 
 
 60 
						}
- 
 
 61 
						while((fabs(l-l_old) > fabs(EV_THRESH * l)) || k<2);
67
62
68 
						// Update the solution by adding the eigenvector to Q and
 63 
						// Update the solution by adding the eigenvector to Q and
69 
						// the eigenvalue to the diagonal of L.
 64 
						// the eigenvalue to the diagonal of L.
70 
						Q[i] = q;
 65 
						Q[i] = q;
71 
						L[i][i] = l;
 66 
						L[i][i] = l;