Subversion Repositories gelsvn

Rev

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

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