Subversion Repositories gelsvn

//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;

int main()

{

		// Matrix row and column dimensions

		int N=4,M=4;

		// Fill the matrix with random values.

		srand(time(0));

		CMatrix mat(N,M,0);

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

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

						mat.set(i,j,rand()/double(RAND_MAX)-0.5);

		// Create a vector of random values

		CVector b(N);

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

				b.set(i, rand()/double(RAND_MAX)-0.5);

		// Some output

		cout << "Matrix rows, cols = "; 

		cout << mat.Rows() << " , " << mat.Cols() << endl;

		cout << "Matrix: " << mat << endl;

		cout << "Vector: " << b << endl;

		// Use lapack function to do least squares solution

		CVector x(M);

		x = LinearLSSolve(mat,b);

		cout << "Least squares solution: " << x << endl;

		// Of course if the system has full rank, we can use

		//CVector LinearSolve(const CMatrix& A,const CVector&b);

		// Manual (and less efficient least squares solution ...

		CMatrix mat_T = mat.Transposed(); 

		CMatrix mat2 = mat_T * mat;

		Invert(mat2);

		mat2 = mat2 * mat_T;

		CVector x2(M);

		x2 = mat2 * b;

		cout << "Least squares solution: " << x2 << endl;

}