Subversion Repositories gelsvn

Rev

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

Rev 54 Rev 58
Line 18... Line 18...
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 
 
Line 37... Line 37...
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;