MatrixTest.h
1/***************************************************************************
2 *
3 * Copyright (C) 2003 by Willem van Straten
4 * Licensed under the Academic Free License version 2.1
5 *
6 ***************************************************************************/
7#include <iostream>
8#include <string>
9
10#include "random.h"
11
12template <class Type1, class Type2, class Element>
13class MatrixTest {
14
15 public:
16
17 MatrixTest ()
18 { verbose = false; scale=1e3; tolerance = 1e-12; random_init(); }
19
20 virtual ~MatrixTest () {}
21
22 virtual void runtest (unsigned nloops);
23
24 bool verbose;
25
26 double scale;
27
28 double tolerance;
29
30};
31
32template <class T, class U>
33T max (const T& a, const U& b)
34{
35 if (a>b)
36 return a;
37 else
38 return b;
39}
40
41template <class Type1, class Type2, class Element>
42void test_matrix (const Type1& t1_1, const Type1& t1_2, const Type2& t2_null,
43 const Element& p, double scale, double tolerance, bool verbose)
44{
45 Type1 temp1;
46 Type1 temp2;
47
48 Element one (1.0);
49
50 //
51 // test operator =
52 //
53 if (verbose)
54 std::cerr << "\n*************** Testing operator =" << std::endl;
55
56 if (verbose)
57 std::cerr << "Type1 = Type1" << std::endl;
58 temp1 = t1_1;
59
60 if (verbose) std::cerr << "Type2 = Type1" << std::endl;
61 Type2 t2_1 = temp1;
62
63 if (verbose) std::cerr << "Type2 = Type2" << std::endl;
64 Type2 t2_2 = t2_1;
65
66 if (verbose) std::cerr << "Type1 = Type2" << std::endl;
67 temp1 = t2_2;
68
69 if (verbose) std::cerr << "Type1 == Type1" << std::endl;
70 if (t1_1 != temp1)
71 throw std::string ("test_matrix "
72 "Error casting from Type1 to Type2");
73
74 if (verbose) std::cerr << "Type1 == Type2" << std::endl;
75 if (t1_1 != t2_1)
76 throw std::string ("test_matrix "
77 "Error casting from Type1 to Type2");
78
79 //
80 // test operator + and - (implicitly tests operator+= and -=)
81 //
82 if (verbose)
83 std::cerr << "\n*************** Testing operator + and -" << std::endl;
84
85 t2_2 = t1_2;
86
87 if (verbose) std::cerr << "Type1 + Type2" << std::endl;
88 temp1 = t1_1 + t2_2;
89
90 t2_2 = temp1;
91
92 if (verbose) std::cerr << "Type2 - Type1" << std::endl;
93 t2_1 = t2_2 - t1_2;
94
95 if (verbose) std::cerr << "Type1 == Type2" << std::endl;
96
97 if (norm(t1_1 - t2_1)/norm(t2_1) > tolerance)
98 {
99 std::cerr << "t1_1 = " << t1_1 << " != "
100 << "t2_1 = " << t2_1 << std::endl;
101
102 throw std::string ("test_matrix "
103 "Error adding and subtracting Type1 and Type2");
104 }
105
106
107 //
108 // test operator * and / Element (implicitly tests operator*= and /=)
109 //
110 if (verbose)
111 std::cerr << "\n*************** Testing operator * and / Element"
112 << std::endl;
113
114 Element dfactor;
115 random_value (dfactor, scale);
116
117 t2_1 = t1_1;
118
119 if (verbose) std::cerr << "Element * Type2" << std::endl;
120 t2_2 = dfactor * t2_1;
121
122 if (verbose) std::cerr << "Type2 / Element" << std::endl;
123 t2_1 = t2_2 / dfactor;
124
125 if (verbose) std::cerr << "Type1 == Type2" << std::endl;
126 if (norm(t2_1 - t1_1)/norm(t2_1) > tolerance)
127 {
128 std::cerr << "Error Type2=" << t2_1
129 << " * / Element=" << dfactor
130 << " != Type1=" << t1_1 << std::endl;
131 throw std::string ("test_matrix Error in Type2 * / Element");
132 }
133
134 //
135 // test operator Type1 * Type1 (implicitly tests operator*=)
136 //
137 if (verbose) std::cerr << "\n*************** Testing operator * Type1" << std::endl;
138
139 if (verbose) std::cerr << "Type1::identity=" << Type1::identity() << std::endl;
140
141 if (verbose) std::cerr << "Type1 * Type1::identity" << std::endl;
142 t2_1 = t1_1 * Type1::identity();
143
144 if (verbose) std::cerr << "Type1 == Type2" << std::endl;
145 if (t1_1 != t2_1)
146 throw std::string ("test_matrix Error Type1 * Type1::identity");
147
148
149 //
150 // test complex conjugate
151 //
152 if (verbose) std::cerr << "\n*************** Testing conjugate" << std::endl;
153
154 if (verbose) std::cerr << "conj (t1_1 * t1_2)" << std::endl;
155 t2_1 = conj (t1_1 * t1_2);
156
157 if (verbose) std::cerr << "conj (t1_1) * conj (t1_2)" << std::endl;
158 t2_2 = conj (t1_1) * conj (t1_2);
159
160 if (verbose) std::cerr << "Type2 == Type2" << std::endl;
161 if (norm(t2_1 - t2_2)/norm(t2_1) > tolerance)
162 {
163 std::cerr << "1=" << t2_1 << "\n2=" << t2_2 << std::endl;
164 throw std::string ("test_matrix Error Conjugate");
165 }
166
167 //
168 // test Hermitian
169 //
170 if (verbose) std::cerr << "\n*************** Testing Hermitian" << std::endl;
171
172 t2_1 = t1_1;
173 t2_2 = t1_2;
174
175 if (verbose) std::cerr << "herm (t2_1 * t2_2)" << std::endl;
176 temp1 = herm (t2_1 * t2_2);
177
178 if (verbose) std::cerr << "herm (t2_2) * herm (t2_1)" << std::endl;
179 temp2 = herm (t2_2) * herm (t2_1);
180
181 if (verbose) std::cerr << "Type2 == Type2" << std::endl;
182
183 if ( norm(temp2 - temp1)/norm(temp2) > tolerance) {
184 std::cerr << "Error Hermitian"
185 "\n M1=" << t1_1 <<
186 "\n M2=" << t1_2 <<
187 "\n herm(M1*M2)=" << temp1 <<
188 "\n herm(M2)*herm(M1)=" << temp2 << std::endl;
189
190 throw std::string ("test_matrix Error Hermitian");
191 }
192
193 //
194 // test det(Type1) function
195 //
196 if (verbose) std::cerr << "\n*************** Testing determinant" << std::endl;
197
198 t2_1 = t1_1;
199
200 if (verbose) std::cerr << "det(Type2)" << std::endl;
201 Element det1 (det(t2_1));
202
203 if (verbose) std::cerr << "Type2 /= sqrt(det)" << std::endl;
204 t2_1 /= sqrt(det1);
205
206 if (verbose) std::cerr << "det(Type2)" << std::endl;
207 Element det2 (det(t2_1));
208
209 if (norm(det2 - one) > tolerance)
210 throw std::string ("test_matrix "
211 "Error normalizing by square root of determinant");
212
213
214 //
215 // test norm and trace
216 //
217 if (verbose)
218 std::cerr << "\n*************** Testing norm and trace" << std::endl;
219
220 t2_2 = t1_1;
221
222 if (verbose) std::cerr << "norm (Type2)" << std::endl;
223 double variance1 = norm(t2_2);
224
225 if (verbose) std::cerr << "herm (Type2)" << std::endl;
226 t2_1 = t2_2*herm(t2_2);
227
228 if (verbose) std::cerr << "trace (Type2)" << std::endl;
229 Element variance2 (trace (t2_1));
230
231 if ( norm(variance1 - variance2) > tolerance )
232 {
233 std::cerr << "Q=" << t2_2 << std::endl;
234 std::cerr << "herm(Q)=" << herm(t2_2) << std::endl;
235 std::cerr << "Q*herm(Q)=" << t2_1 << std::endl;
236 std::cerr << "norm(Q) = " << variance1 << std::endl;
237 std::cerr << "trace (Q * herm(Q)) = " << variance2 << std::endl;
238
239 throw std::string ("test_matrix "
240 "Error norm(Q) != trace (Q * herm(Q))");
241 }
242
243 //
244 // test Type1::inv method
245 //
246 if (verbose) std::cerr << "\n*************** Testing inverse" << std::endl;
247
248 if (verbose) std::cerr << "Type1::inv" << std::endl;
249 t2_2 = inv(t1_1);
250
251 if (verbose) std::cerr << "Type2 * Type2" << std::endl;
252 t2_1 = Type2(t1_1) * t2_2;
253
254 double threshold = tolerance * max (norm(t1_1), norm(t2_2));
255
256 if (verbose) std::cerr << "norm (Type2)" << std::endl;
257 variance1 = norm (Type2::identity() - t2_1);
258
259 if ( variance1 > threshold ) {
260 std::cerr << "Error multipying Type1 by its inverse"
261 "\n Q=" << t1_1 <<
262 "\n inv(Q)=" << t2_2 <<
263 "\n Q*inv(Q)=" << t2_1 << std::endl;
264
265 throw std::string ("test_matrix "
266 "Error multipying Type1 by its inverse");
267 }
268}
269
270template <class Type1, class Type2, class Element>
271void MatrixTest<Type1, Type2, Element>::runtest (unsigned nloops)
272{
273 //
274 // test constructors
275 //
276 if (verbose)
277 std::cerr << "\n*************** Testing constructors" << std::endl;
278
279 if (verbose)
280 std::cerr << "Type1 null construct" << std::endl;
281 Type1 temp;
282
283 if (verbose)
284 std::cerr << "Element array (static)" << std::endl;
285
286 Element random_elements[4];
287 for (int i=0; i<4; i++)
288 random_value (random_elements[i], scale);
289
290 if (verbose)
291 std::cerr << "Type1 construct from 4 Elements" << std::endl;
292 Type1 t1_1 (random_elements[0], random_elements[1],
293 random_elements[2], random_elements[3]);
294
295 if (verbose)
296 std::cerr << "Type1 construct from Type1" << std::endl;
297 Type1 t1_2 (t1_1);
298
299 if (verbose)
300 std::cerr << "TEST: Type1 == Type1" << std::endl;
301 if (t1_1 != t1_2)
302 throw std::string ("MatrixTest::runtest "
303 "Error in constructing Type1 from Type1");
304
305 if (verbose)
306 std::cerr << "Type1 = Type1" << std::endl;
307 temp = t1_1;
308
309 if (verbose)
310 std::cerr << "TEST: Type1 == Type1" << std::endl;
311 if (t1_1 != temp)
312 throw std::string ("MatrixTest::runtest "
313 "Error setting Type1 equal to Type1");
314
315 for (unsigned iloop=0; iloop<nloops; iloop++) {
316
317 Type1 t1_1, t1_2;
318
319 random_vector (t1_1, scale);
320 random_vector (t1_2, scale);
321
322 Type2 t2_null;
323 Element null;
324
325 try {
326 test_matrix (t1_1, t1_2, t2_null, null, scale, tolerance, verbose);
327 }
328 catch (std::string& error) {
329 std::cerr << "MatrixTest::runtest error on loop "<< iloop <<"/"<< nloops
330 << std::endl;
331 throw error;
332 }
333
334 }
335}
Generated using doxygen 1.14.0