Subversion Repositories gelsvn

Rev

Rev 326 | Go to most recent revision | Only display areas with differences | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 326 Rev 343
1 
#include <iostream>
 1 
#include <iostream>
2 
#include <algorithm>
 2 
#include <algorithm>
3 
 
 3 
 
4 
#include "CGLA/Mat4x4d.h"
 4 
#include "CGLA/Mat4x4d.h"
5 
#include "CGLA/Mat4x4f.h"
 5 
#include "CGLA/Mat4x4f.h"
6 
#include "CGLA/Mat2x2f.h"
 6 
#include "CGLA/Mat2x2f.h"
7 
#include "CGLA/Mat3x3f.h"
 7 
#include "CGLA/Mat3x3f.h"
8 
#include "CGLA/Mat2x3f.h"
 8 
#include "CGLA/Mat2x3f.h"
9 
#include "CGLA/eigensolution.h"
 9 
#include "CGLA/eigensolution.h"
10 
 
 10 
 
11 
#include "CGLA/Vec2f.h"
 11 
#include "CGLA/Vec2f.h"
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 
{
24 
 
 76 
 
25 
	Mat2x3f m23(Vec3f(1,1,1),Vec3f(2,2,2));
 77 
	Mat2x3f m23(Vec3f(1,1,1),Vec3f(2,2,2));
26 
	Mat3x2f m32;
 78 
	Mat3x2f m32;
27
79
28 
	// Try transpose of non-sq matrix
 80 
	// Try transpose of non-sq matrix
29 
	transpose(m23, m32);
 81 
	transpose(m23, m32);
30 
	cout << "2 by 3 " << m23 << endl;
 82 
	cout << "2 by 3 " << m23 << endl;
31 
	cout << "Transpose " << m32 << endl;
 83 
	cout << "Transpose " << m32 << endl;
32
84
33
85
34 
	m32[0] = Vec2f(1,2);
 86 
	m32[0] = Vec2f(1,2);
35 
	m32[1] = Vec2f(10,20);
 87 
	m32[1] = Vec2f(10,20);
36 
	m32[2] = Vec2f(100,200);
 88 
	m32[2] = Vec2f(100,200);
37 
 
 89 
 
38 
	Mat2x2f m22;
 90 
	Mat2x2f m22;
39
91
40 
	// Multiply 2x3 and 3x2 matrices
 92 
	// Multiply 2x3 and 3x2 matrices
41 
	mul(m23, m32, m22);
 93 
	mul(m23, m32, m22);
42 
	cout << "Multiplication of non. sq. mat " << m22 << endl;
 94 
	cout << "Multiplication of non. sq. mat " << m22 << endl;
43 
 
 95 
 
44 
	// Check matrix mul using * operator for sq matrices.
 96 
	// Check matrix mul using * operator for sq matrices.
45 
	Mat2x2f m22_2(1,2,3,4);
 97 
	Mat2x2f m22_2(1,2,3,4);
46 
	cout << "multiplication of sq. mats " << m22_2 * m22 << endl;
 98 
	cout << "multiplication of sq. mats " << m22_2 * m22 << endl;
47 
	cout << "Trace " << trace(m22) << endl;
 99 
	cout << "Trace " << trace(m22) << endl;
48
100
49 
	m22_2 += m22 + m22;
 101 
	m22_2 += m22 + m22;
50 
	m22_2 -= 2* m22;
 102 
	m22_2 -= 2* m22;
51
103
52 
	Mat4x4f m44(0.0f);
 104 
	Mat4x4f m44(0.0f);
53 
	m44[0][0] = 1;
 105 
	m44[0][0] = 1;
54 
	m44[1][2] = 2;
 106 
	m44[1][2] = 2;
55 
	m44[2][1] = 3;
 107 
	m44[2][1] = 3;
56 
	m44[3][3] = 4;
 108 
	m44[3][3] = 4;
57 
 
 109 
 
58 
	Mat4x4f m44_2 = transpose(m44);
 110 
	Mat4x4f m44_2 = transpose(m44);
59 
	cout << m44 << m44_2 << endl;
111 
	cout << m44 << m44_2 << endl;
60 
 
 112 
 
61 
	// Compile should fail if line below is uncommented
 113 
	// Compile should fail if line below is uncommented
62 
	// mul(m23, m44, m22);
 114 
	// mul(m23, m44, m22);
63 
	// Compile should fail if line below is uncommented
 115 
	// Compile should fail if line below is uncommented
64 
	// transpose(m23,m44);
 116 
	// transpose(m23,m44);
65 
 
 117 
 
66 
	cout << "Determinant of 4x4 " << m44 << determinant(m44) << endl;
 118 
	cout << "Determinant of 4x4 " << m44 << determinant(m44) << endl;
67 
	cout << "Determinant of 2x2 " << m22 << determinant(m22) << endl;
 119 
	cout << "Determinant of 2x2 " << m22 << determinant(m22) << endl;
68 
 
 120 
 
69 
	Mat4x4f mna(Vec4f(10,120,10,40),
 121 
	Mat4x4f mna(Vec4f(10,120,10,40),
70 
							Vec4f(43,10,31254,10),
 122 
							Vec4f(43,10,31254,10),
71 
							Vec4f(43,11,54,10),
 123 
							Vec4f(43,11,54,10),
72 
							Vec4f(0,0,0,1));
 124 
							Vec4f(0,0,0,1));
73 
 
 125 
 
74 
	cout << fixed << endl;
 126 
	cout << fixed << endl;
75 
	try
 127 
	try
76 
		{
 128 
		{
77 
			cout << "Invert: " << invert_affine(mna) << endl;
 129 
			cout << "Invert: " << invert_affine(mna) << endl;
78 
			cout << "Test: " << invert_affine(mna)*mna << endl;
 130 
			cout << "Test: " << invert_affine(mna)*mna << endl;
79 
		}
 131 
		}
80 
	catch(const Mat4x4fException& me)
 132 
	catch(const Mat4x4fException& me)
81 
		{
 133 
		{
82 
			me.print(cout);
 134 
			me.print(cout);
83 
		}
 135 
		}
84 
	try
 136 
	try
85 
		{
137 
		{
86 
			cout << "Invert " << invert(mna) << endl;
 138 
			cout << "Invert " << invert(mna) << endl;
87 
			cout << "test " << invert(mna)*mna << endl;
 139 
			cout << "test " << invert(mna)*mna << endl;
88 
		}
 140 
		}
89 
	catch(const Mat4x4fException& me)
 141 
	catch(const Mat4x4fException& me)
90 
		{
 142 
		{
91 
			me.print(cout);
 143 
			me.print(cout);
92 
		}
 144 
		}
93 
 
 145 
 
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);
117 
 
 157 
 
118 
	cout << "The matrix " << m33 << " is copied to " << m44 << endl;
 158 
	cout << "The matrix " << m33 << " is copied to " << m44 << endl;
119 
 
 159 
 
120 
	cout << "Determinant of singular 4x4 matrix" << endl;
 160 
	cout << "Determinant of singular 4x4 matrix" << endl;
121 
 
 161 
 
122 
 	Mat4x4d m_odd(Vec4d(0, 373, 139129, 1),
 162 
 	Mat4x4d m_odd(Vec4d(0, 373, 139129, 1),
123 
 								Vec4d(-373, -125, 154754, 1),
 163 
 								Vec4d(-373, -125, 154754, 1),
124 
 								Vec4d(-50.6168, 21.2469, 3013.5, 100),
 164 
 								Vec4d(-50.6168, 21.2469, 3013.5, 100),
125 
 								Vec4d(-50.6168, 21.2469, 3013.5, 100));
 165 
 								Vec4d(-50.6168, 21.2469, 3013.5, 100));
126 
 
 166 
 
127 
	cout << m_odd << endl;
 167 
	cout << m_odd << endl;
128 
	cout << "Determinant: " << determinant(m_odd) << endl;
 168 
	cout << "Determinant: " << determinant(m_odd) << endl;
129 
 
 169 
 
130 
	Mat4x4d m_odd_t = transpose(m_odd);
 170 
	Mat4x4d m_odd_t = transpose(m_odd);
131 
 
 171 
 
132 
	cout << m_odd_t << endl;
 172 
	cout << m_odd_t << endl;
133 
	cout << "Det of transpose:" << determinant(m_odd_t) << endl;
 173 
	cout << "Det of transpose:" << determinant(m_odd_t) << endl;
134 
 
 174 
 
135 
	cout << "Determinant of random matrix (MatLab says 0.0491)" << endl;
 175 
	cout << "Determinant of random matrix (MatLab says 0.0491)" << endl;
136
176
137 
 
 177 
 
138 
	Mat4x4d Mrand(Vec4d(0.9169,0.3529,0.2028,0.1988),
 178 
	Mat4x4d Mrand(Vec4d(0.9169,0.3529,0.2028,0.1988),
139 
								Vec4d(0.4103,0.8132,0.1987,0.0153),
 179 
								Vec4d(0.4103,0.8132,0.1987,0.0153),
140 
								Vec4d(0.8936,0.0099,0.6038,0.7468),
 180 
								Vec4d(0.8936,0.0099,0.6038,0.7468),
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 
}
148 
 
 191