Subversion Repositories gelsvn

Rev

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

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