Subversion Repositories gelsvn

#include "eigensolution.h"

#include "Mat2x2f.h"

#include "Mat3x3f.h"

#include "Mat4x4f.h"

#include <iostream>

using namespace std;

namespace

{

	// During experiments 925 iterations were observed for a hard problem

	// Ten times that should be safe.

	const int KMAX = 10000;

	// The threshold below is the smallest that seems to give reliable

	// solutions.

	const float EV_THRESH = 0.0000001;

}

namespace CGLA

{

	template <class MT>

	int power_eigensolution(const MT& Ap, MT& Q, MT& L, int max_sol)

	{

		typedef typename MT::VectorType VT;

		MT A = Ap;

		int n = s_min(MT::get_v_dim(), max_sol);

		for(int i=0;i<n;++i)

			{

				// Seed the eigenvector estimate

				VT q(1);

				float l=0,l_old;

				// As long as we haven't reached the max iterations and the

				// eigenvalue has not converged, do

				int k=0;

				for(; k<2 || k<KMAX && (l-l_old > EV_THRESH * l) ; ++k)

					{

						// Multiply the eigenvector estimate onto A

						const VT z = A * q;

						// Check that we did not get the null vector.

						float z_len = z.length();

						if(z_len < EV_THRESH)

							return i;

						// Normalize to get the new eigenvector

						q = z/z_len;

						// Record the old eigenvalue estimate and get a new estimate.

						l_old = l;

						l = dot(q, A * q);

					}

				// If we hit the max iterations, we also don't trust the eigensolution

				if(k==KMAX)

					return i;

				// Update the solution by adding the eigenvector to Q and

				// the eigenvalue to the diagonal of L.

				Q[i] = q;

				L[i][i] = l;

				// Update A by subtracting the subspace represented by the 

				// eigensolution just found. This is called the method of 

				// deflation.

				MT B;

				outer_product(q,q,B);

				A = A - l * B;

			}

		return n;

	}

	/* There is no reason to put this template in a header file, since 

		 we will only use it on matrices defined in CGLA. Instead, we 

		 explicitly instantiate the function for the square matrices

		 of CGLA */

	template int power_eigensolution<Mat2x2f>(const Mat2x2f&,

																						Mat2x2f&,Mat2x2f&,int);

	template int power_eigensolution<Mat3x3f>(const Mat3x3f&,

																						Mat3x3f&,Mat3x3f&,int);

	template int power_eigensolution<Mat4x4f>(const Mat4x4f&,

																						Mat4x4f&,Mat4x4f&,int);

}