Subversion Repositories gelsvn

/* LinAlg EigenSolutionSym test program. Andreas Bærentzen 2008. 

	This simple program generates 30 moderately sized dense symmetric real matrices and computes

	the full eigensolution (eigenvalues and eigenvectors). For each matrix and each eigenvector it

	multiplies that eigenvector back onto the matrix and checks the difference between the new vector

	and the eigenvector times the eigenvalue. It also computes the difference between the eigenvalue

	returned and the eigenvalue computed as the norm of a unit length eigenvector multiplied onto the

	matrix.

*/

//disable the stl warning of debug information conflicting with names of longer than 255 char.

#pragma warning( disable :4786 )  

#include <iostream>

#include <vector>

#include "LinAlg/Matrix.h"

#include "LinAlg/Vector.h"

#include "LinAlg/LapackFunc.h"

using namespace std;

using namespace LinAlg;

using namespace CGLA;

int main()

{

		gel_srand(0);

		double max_err_eigen=0;

		double max_err_eigen_vec=0;

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

		{

			int N = gel_rand()%50+2;

			CMatrix mat(N,N,0);

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

				for(int j=0;j<N;++j)

				{

					double x = gel_rand()/double(GEL_RAND_MAX)-0.5;

					mat.set(i,j,x);

					mat.set(j,i,x);

				}

			CVector b;

			CMatrix Q=mat;

			int info = EigenSolutionsSym(Q,b);

			cout << "Info = " << info << endl;

			for(int i=0;i<b.Length(); ++i)

			{

				CVector q;

				Q.GetRow(q, i);

				q /= q.Norm();

				CVector q2;

				q2 = mat * q;

				double eigen_check = q2.Norm();

				double err_eigen =   (fabs(b[i])-fabs(eigen_check));

				double err_eigen_vec = (q*b[i]-q2).Norm();

				cout << "l(i)-Norm(Q*e(i)) " << err_eigen << endl;

				cout << "Norm(e(i)*l(i)-Q*e(i)) " << err_eigen_vec << endl;

				max_err_eigen = max(err_eigen, max_err_eigen);

				max_err_eigen_vec = max(err_eigen_vec, max_err_eigen_vec);

			}

			cout << "-------" << endl;

		}

		cout << "Maximum eigenvalue error :" << max_err_eigen << endl;

		cout << "Maximum eigenvector error :" << max_err_eigen_vec << endl;

}