Subversion Repositories gelsvn

Rev

Rev 417 | Only display areas with differences | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 417 Rev 630
1 
//disable the stl warning of debug information conflicting with names of longer than 255 char.
 1 
//disable the stl warning of debug information conflicting with names of longer than 255 char.
2 
#pragma warning( disable :4786 )
2 
#pragma warning( disable :4786 )
3 
 
 3 
 
4 
#include <iostream>
 4 
#include <iostream>
5 
#include <vector>
 5 
#include <vector>
6 
#include "CGLA/CGLA.h"
 6 
#include "CGLA/CGLA.h"
7 
#include "LinAlg/Matrix.h"
 7 
#include "LinAlg/Matrix.h"
8 
#include "LinAlg/Vector.h"
 8 
#include "LinAlg/Vector.h"
9 
#include "LinAlg/LapackFunc.h"
 9 
#include "LinAlg/LapackFunc.h"
10 
 
 10 
 
11 
using namespace std;
 11 
using namespace std;
12 
using namespace LinAlg;
 12 
using namespace LinAlg;
13 
using namespace CGLA;
 13 
using namespace CGLA;
14 
 
 14 
 
15 
int main()
 15 
int main()
16 
{
 16 
{
17 
		// Matrix row and column dimensions
 17 
		// Matrix row and column dimensions
18 
		int N=4,M=4;
 18 
		int N=4,M=4;
19 
 
 19 
 
20 
		// Fill the matrix with random values.
 20 
		// Fill the matrix with random values.
21 
		gel_srand(0);
 21 
		gel_srand(0);
22 
		CMatrix mat(N,M,0);
 22 
		CMatrix mat(N,M,0);
23 
		for(int i=0;i<N;++i)
 23 
		for(int i=0;i<N;++i)
24 
				for(int j=0;j<M;++j)
 24 
				for(int j=0;j<M;++j)
25 
						mat.set(i,j,gel_rand()/double(GEL_RAND_MAX)-0.5);
 25 
						mat.set(i,j,gel_rand()/double(GEL_RAND_MAX)-0.5);
26 
 
 26 
 
27 
		// Create a vector of random values
 27 
		// Create a vector of random values
28 
		CVector b(N);
 28 
		CVector b(N);
29 
		for(int i=0;i<N;++i)
 29 
		for(int i=0;i<N;++i)
30 
				b.set(i, gel_rand()/double(GEL_RAND_MAX)-0.5);
 30 
				b.set(i, gel_rand()/double(GEL_RAND_MAX)-0.5);
31
31
32 
		// Some output
 32 
		// Some output
33 
		cout << "Matrix rows, cols = ";
33 
		cout << "Matrix rows, cols = ";
34 
		cout << mat.Rows() << " , " << mat.Cols() << endl;
 34 
		cout << mat.Rows() << " , " << mat.Cols() << endl;
35 
		cout << "Matrix: " << mat << endl;
 35 
		cout << "Matrix: " << mat << endl;
36 
		cout << "Vector: " << b << endl;
 36 
		cout << "Vector: " << b << endl;
37
37
38 
		// Use lapack function to do least squares solution
 38 
		// Use lapack function to do least squares solution
39 
		CVector x(M);
 39 
		CVector x(M);
40 
		x = LinearLSSolve(mat,b);
 40 
		x = LinearLSSolve(mat,b);
41 
		cout << "Least squares solution: " << x << endl;
 41 
		cout << "Least squares solution: " << x << endl;
42
42
43 
		// Of course if the system has full rank, we can use
 43 
		// Of course if the system has full rank, we can use
44 
		//CVector LinearSolve(const CMatrix& A,const CVector&b);
 44 
		//CVector LinearSolve(const CMatrix& A,const CVector&b);
45 
 
 45 
 
46 
		// Manual (and less efficient least squares solution ...
 46 
		// Manual (and less efficient least squares solution ...
47 
		CMatrix mat_T = mat.Transposed();
47 
		CMatrix mat_T = mat.Transposed();
48 
		CMatrix mat2 = mat_T * mat;
 48 
		CMatrix mat2 = mat_T * mat;
49 
		Invert(mat2);
 49 
		Invert(mat2);
50 
		mat2 = mat2 * mat_T;
 50 
		mat2 = mat2 * mat_T;
51 
 
 51 
 
52 
		CVector x2(M);
 52 
		CVector x2(M);
53 
		x2 = mat2 * b;
 53 
		x2 = mat2 * b;
54 
		cout << "Least squares solution: " << x2 << endl;
 54 
		cout << "Least squares solution: " << x2 << endl;
55 
 
 55 
 
56 
}
 56 
}
57 
 
 57