Subversion Repositories gelsvn

Rev

Rev 326 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 326 Rev 343
Line 12... Line 12...
12 
#include "CGLA/Vec2i.h"
 12 
#include "CGLA/Vec2i.h"
13 
#include "CGLA/Vec3i.h"
 13 
#include "CGLA/Vec3i.h"
14 
#include "CGLA/Vec3f.h"
 14 
#include "CGLA/Vec3f.h"
15 
#include "CGLA/Vec3Hf.h"
 15 
#include "CGLA/Vec3Hf.h"
16 
 
 16 
 
- 
 
 17 
#include "LinAlg/LapackFunc.h"
- 
 
 18 
 
17 
using namespace std;
 19 
using namespace std;
18 
using namespace CGLA;
 20 
using namespace CGLA;
- 
 
 21 
using namespace LinAlg;
- 
 
 22 
 
- 
 
 23 
 
- 
 
 24 
double percentdiff(double a, double b)
- 
 
 25 
{
- 
 
 26 
	return 100.0*fabs(a-b)/fabs(a);
- 
 
 27 
}
- 
 
 28 
 
- 
 
 29 
void power_eigensolution_test()
- 
 
 30 
{
- 
 
 31 
	srand(0);
- 
 
 32 
	double avg_trace_diff=0;
- 
 
 33 
	double avg_det_diff=0;
- 
 
 34 
	for (int k=0;k<1000;++k)
- 
 
 35 
	{
- 
 
 36 
		Mat4x4d M;
- 
 
 37 
		for(int i=0;i<4;++i)
- 
 
 38 
			for(int j=i;j<4;++j)
- 
 
 39 
				M[i][j] = M[j][i] = 1e-5*(rand()-RAND_MAX/2);
- 
 
 40
- 
 
 41 
		Mat4x4d Q,L;
- 
 
 42 
		power_eigensolution(M,Q,L);
- 
 
 43
- 
 
 44 
		double det_diff = percentdiff(determinant(M), L[0][0]*L[1][1]*L[2][2]*L[3][3]) ;
- 
 
 45 
		double trace_diff= percentdiff(trace(M), L[0][0]+L[1][1]+L[2][2]+L[3][3]);
- 
 
 46
- 
 
 47 
		if(det_diff > 1)
- 
 
 48 
		{
- 
 
 49 
			cout << "Matrix " << M << " has Q matrix " <<  Q << " and L " << L << endl;
- 
 
 50 
			CMatrix CM(M);
- 
 
 51 
			CVector CV = SVD(CM);
- 
 
 52 
			cout << "Abs. eigenvalues computed using SVD ";
- 
 
 53 
			for(int i=0;i<4;++i)
- 
 
 54 
				cout << (CV[i]) <<  " " ;
- 
 
 55 
			cout << endl;
- 
 
 56 
			cout << "Det(M) = " << determinant(M) << " trace(m) = " << trace(M) <<endl;
- 
 
 57
- 
 
 58 
			cout << "Det(M) vs product of eigenvalues, difference = "
- 
 
 59 
				<< det_diff << " %" << endl;
- 
 
 60 
			avg_det_diff += det_diff;
- 
 
 61 
			cout << "Trace of M vs sum of eigenvalues, difference = "
- 
 
 62 
				<< trace_diff << " %" << endl;
- 
 
 63 
			avg_trace_diff += trace_diff;
- 
 
 64 
			cout << endl << "-----------------------------------" << endl;
- 
 
 65 
		}
- 
 
 66 
	}
- 
 
 67 
	cout << "average trace diff = " << (avg_trace_diff / 1000) << " %" << endl;
- 
 
 68 
	cout << "average det diff = " << (avg_det_diff / 1000) << " %" << endl;
- 
 
 69 
 
- 
 
 70 
}
19 
 
 71 
 
20 
 
 72 
 
21 
/* This is a non-exhaustive test program for CGLA */
 73 
/* This is a non-exhaustive test program for CGLA */
22 
int main()
 74 
int main()
23 
{
 75 
{
Line 94... Line 146...
94 
	Vec2f v2(2,3);
 146 
	Vec2f v2(2,3);
95 
	Vec3f v3(1,2,3);
 147 
	Vec3f v3(1,2,3);
96 
	Mat2x3f m23_2;
 148 
	Mat2x3f m23_2;
97 
	outer_product(v2,v3,m23_2);
 149 
	outer_product(v2,v3,m23_2);
98 
	cout << "outer product of " << v2 << " , " << v3 << " = " << m23_2 << endl;
 150 
	cout << "outer product of " << v2 << " , " << v3 << " = " << m23_2 << endl;
99 
 
 - 
 
100 
	Mat4x4f A(Vec4f(1, 2, 3, 4),
- 
 
101 
						Vec4f(2, 9,-7,2),
- 
 
102 
						Vec4f(3,-7, 8, 1),
 - 
 
103 
						Vec4f(4,2, 31,12)
 - 
 
104 
						);
 - 
 
105 
	Mat4x4f Q,L;
 - 
 
106
151
107 
	int n = power_eigensolution(A, Q, L);
 - 
 
108 
 
 - 
 
109 
	cout << "The " << n << " eigensolutions of " << A << " are ";
 - 
 
110 
	cout << Q << L << endl;
 - 
 
111 
 
 - 
 
112 
	m44 = identity_Mat4x4f();
 152 
	m44 = identity_Mat4x4f();
113 
 
 153 
 
114 
	Mat3x3f m33(Vec3f(1,2,3), Vec3f(4,5,6), Vec3f(7,8,9));
 154 
	Mat3x3f m33(Vec3f(1,2,3), Vec3f(4,5,6), Vec3f(7,8,9));
115 
 
 155 
 
116 
	copy_matrix(m33, m44);
 156 
	copy_matrix(m33, m44);
Line 141... Line 181...
141 
								Vec4d(0.0579,0.1389,0.2722,0.4451));
 181 
								Vec4d(0.0579,0.1389,0.2722,0.4451));
142
182
143 
	cout << "det " << determinant(Mrand) << endl;
 183 
	cout << "det " << determinant(Mrand) << endl;
144 
	cout << "det of transpose " << determinant(transpose(Mrand)) << endl;
 184 
	cout << "det of transpose " << determinant(transpose(Mrand)) << endl;
145
185
- 
 
 186 
	cout << "Grand test of power eigensolution scheme " << endl;
- 
 
 187 
	power_eigensolution_test();
- 
 
 188 
 
146 
	return 0;
 189 
	return 0;
147 
}
 190 
}