CLHEP VERSION Reference Documentation
CLHEP Home Page
CLHEP Documentation
CLHEP Bug Reports
Main Page
Namespaces
Classes
Files
File List
File Members
Matrix
test
testInversion.cc
Go to the documentation of this file.
1
// -*- C++ -*-
2
// $Id: testInversion.cc,v 1.2 2003/08/13 20:00:12 garren Exp $
3
//
4
// This file is a part of CLHEP - a Class Library for High Energy Physics.
5
//
6
// This is a collection of parts of programs that are useful
7
// for testing matrix inversion algorithms
8
// 9/97, Mario Stanke
9
10
#include <time.h>
11
#include <iostream>
12
13
#include "CLHEP/Matrix/defs.h"
14
#include "CLHEP/Matrix/Matrix.h"
15
#include "CLHEP/Matrix/SymMatrix.h"
16
#include "CLHEP/Matrix/DiagMatrix.h"
17
18
using
std::cout;
19
using
std::endl;
20
21
using namespace
CLHEP;
22
23
int
main
()
24
{
25
//int n , i, j, k, ierr1, ierr2;
26
int
n
, j, ierr1, ierr2;
27
time_t zeit1, zeit2;
28
29
// ****compare the speed of inversion algorithms
30
31
HepSymMatrix
A(5,1);
32
33
// for (j=1;j <= 100; j++)
34
// for (k=1; k<=j; k++)
35
// A(j,k)=rand()%9-5;
36
37
A(1,1)=6;
38
A(1,2)=.8545;
39
A(1,3)=-.684651;
40
A(1,4)=.372547;
41
A(1,5)=.68151;
42
//A(1,6)=.068151;
43
A(2,2)=4;
44
A(2,3)=.5151;
45
A(2,4)=.5151;
46
A(2,5)=.651651;
47
//A(2,6)=.9651651;
48
A(3,3)=5;
49
A(3,4)=.3;
50
A(3,5)=.363;
51
//A(3,6)=.7363;
52
A(4,4)=-3;
53
A(4,5)=-.23753;
54
//A(4,6)=-.23753;
55
A(5,5)=-5;
56
//A(5,6)=-2;
57
//A(6,6)=-3;
58
59
HepSymMatrix
B
(A);
60
HepSymMatrix
D
(A);
61
HepSymMatrix
C(5,0);
62
HepMatrix
M(A);
63
64
cout <<
"M inverse"
<< M.
inverse
(ierr2) << endl;
65
66
C = B.
inverse
(ierr1);
67
D.
invert
(ierr2);
68
69
cout <<
"B "
<< B << endl;
70
cout <<
"B inverse"
<< C << endl;
71
#ifndef INSTALLATION_CHECK
72
cout <<
"B * inverse"
<< B * C << endl;
73
#endif
74
cout <<
"ierr1: "
<< ierr1 << endl;
75
76
cout <<
"D * inverse"
<< D * C << endl;
77
cout <<
"ierr2: "
<< ierr2 << endl;
78
cout <<
"M inverse"
<< M.
inverse
(ierr2) << endl;
79
#ifndef INSTALLATION_CHECK
80
cout <<
"M * inverse"
<< M * M.
inverse
(ierr2)<< endl;
81
#endif
82
cout <<
"ierr2: "
<< ierr2 << endl;
83
84
#ifndef INSTALLATION_CHECK
85
n = 100000;
86
#else
87
n = 10;
88
#endif
89
zeit1 = time(NULL);
90
cout <<
"Executing "
<< n <<
" inversions ..."
<< endl;
91
for
(j=0; j<
n
; j++)
92
{
93
B.
invert
(ierr1);
94
if
(ierr1)
95
cout <<
"B: error in invert"
<< endl;
96
}
97
zeit2 = time(NULL);
98
cout <<
"B: duration of inversion: "
<< zeit2-zeit1 <<
" secs"
<< endl;
99
100
zeit1 = time(NULL);
101
cout <<
"Executing "
<< n <<
" inversions ..."
<< endl;
102
for
(j=0; j<
n
; j++)
103
{
104
D.
invert
(ierr2);
105
if
(ierr2)
106
cout <<
"D: error in invert"
<< endl;
107
}
108
zeit2 = time(NULL);
109
cout << D << endl;
110
cout <<
"D: duration of inversion: "
<< zeit2-zeit1 <<
" secs"
<< endl;
111
112
113
/***** check correctness and compare results of inversion algorithms
114
115
double dist1, dist2;
116
HepSymMatrix A(5,1), B(5), C(5), I(5,1);
117
HepSymMatrix M(5);
118
n = 200000;
119
120
for (j=1; j <= n ; j++)
121
{
122
A(1,1)=(((rand()%7-3)))+ (rand()%3-1)*1e-10;
123
A(1,2)=(((rand()%7-3)))+ (rand()%3-1)*1e-10;
124
A(1,3)=(((rand()%7-3)))+ (rand()%3-1)*1e-10;
125
A(1,4)=(((rand()%7-3)))+ (rand()%3-1)*1e-10;
126
A(1,5)=(((rand()%7-3)))+ (rand()%3-1)*1e-10;
127
//A(1,6)=(((rand()%7-3)))+ (rand()%3-1)*1e-10;
128
A(2,2)=(((rand()%7-3)))+ (rand()%3-1)*1e-10;
129
A(2,3)=(((rand()%7-3)))+ (rand()%3-1)*1e-10;
130
A(2,4)=(((rand()%7-3)))+ (rand()%3-1)*1e-10;
131
A(2,5)=(((rand()%7-3)))+ (rand()%3-1)*1e-10;
132
//A(2,6)=(((rand()%7-3)))+ (rand()%3-1)*1e-10;
133
A(3,3)=(((rand()%7-3)))+ (rand()%3-1)*1e-10;
134
A(3,4)=(((rand()%7-3)))+ (rand()%3-1)*1e-10;
135
A(3,5)=(((rand()%7-3)))+ (rand()%3-1)*1e-10;
136
//A(3,6)=(((rand()%7-3)))+ (rand()%3-1)*1e-10;
137
A(4,4)=(((rand()%7-3)))+ (rand()%3-1)*1e-10;
138
A(4,5)=(((rand()%7-3)))+ (rand()%3-1)*1e-10;
139
//A(4,6)=(((rand()%7-3)))+ (rand()%3-1)*1e-10;
140
A(5,5)=(((rand()%7-3)))+ (rand()%3-1)*1e-10;
141
//A(5,6)=(((rand()%7-3)))+ (rand()%3-1)*1e-10;
142
//A(6,6)=(((rand()%7-3)))+ (rand()%3-1)*1e-10;
143
144
M=B=C=A;
145
146
B.invert(ierr2);
147
M.old_invert(ierr1);
148
149
dist2 = norm_infinity(B*C-I);
150
dist1 = norm_infinity(M*C-I);
151
152
153
if (ierr1 != ierr2)
154
{
155
cout << C << endl;
156
cout << "B " << B << endl;
157
cout << "B*C" << B*C << endl;
158
cout << "M*C" << M*C << endl;
159
cout << "M " << M << endl;
160
cout << "determinant " << C.determinant() << endl;
161
cout << "dist2 " << dist2 << endl;
162
cout << "dist1 " << dist1 << endl;
163
164
cout << "ierr2 " << ierr2 <<endl;
165
cout << "ierr1 " << ierr1 <<endl;
166
167
cout << "j " << j << endl;
168
}
169
170
if (ierr2==0 && dist2 > 1e-4)
171
{
172
cout << "bunch failed to invert but did not indicate" << endl;
173
cout << C << endl;
174
cout << "B " << B << endl;
175
cout << "B*C" << B*C << endl;
176
cout << "M*C" << M*C << endl;
177
cout << "determinant " << C.determinant() << endl;
178
cout << "dist2 " << dist2 << endl;
179
cout << "dist1 " << dist1 << endl;
180
181
cout << "ierr2 " << ierr2 <<endl;
182
cout << "ierr1 " << ierr1 <<endl;
183
184
cout << "j " << j << endl;
185
}
186
if (ierr2==1)
187
{
188
// bunch thinks it is singular
189
if (norm_infinity(M*C-I) < 1e-6)
190
{
191
cout << "bunch said it is singular but old could invert"
192
<< endl;
193
cout << C << endl;
194
cout << "B*C" << B*C << endl;
195
cout << "M*C" << M*C << endl;
196
cout << "determinant " << C.determinant() << endl;
197
198
cout << "dist2" << dist2 << endl;
199
cout << "ierr2 " << ierr2 <<endl;
200
cout << "ierr1 " << ierr1 <<endl;
201
202
cout << "j " << j << endl;
203
}
204
}
205
206
if (ierr1==0 && dist1 > 1e-4 && ierr2==0)
207
{
208
cout << "old failed to invert but did not indicate" << endl;
209
210
cout << C << endl;
211
cout << "B*C" << B*C << endl;
212
cout << "M*C" << M*C << endl;
213
cout << "determinant " << C.determinant() << endl;
214
215
cout << "ierr2 " << ierr2 <<endl;
216
cout << "ierr1 " << ierr1 <<endl;
217
218
cout << "dist1 " << dist1 << endl;
219
cout << "j " << j << endl;
220
return 0;
221
222
}
223
if (ierr1==1)
224
{
225
// old thinks it is singular
226
if (norm_infinity(B*C-I) < 1e-6)
227
{
228
cout << "old said it is singular but bunch could invert"
229
<< endl;
230
231
cout << C << endl;
232
cout << "B*C" << B*C << endl;
233
cout << "M*C" << M*C << endl;
234
cout << "determinant " << C.determinant() << endl;
235
236
cout << "dist1" << dist1 << endl;
237
cout << "dist2" << dist2 << endl;
238
cout << "ierr2 " << ierr2 <<endl;
239
cout << "ierr1 " << ierr1 <<endl;
240
241
cout << "j " << j << endl;
242
return 0;
243
244
}
245
}
246
247
}
248
*/
249
250
/*** one tough symmetric matrix from real physical data
251
252
sm(1,1)=5347.51;
253
sm(1,2)=-142756;
254
sm(1,3)= -1.86624e+06;
255
sm(1,4)=0.0164743;
256
sm(1,5)=0.0915348;
257
sm(1,6)=0.421851;
258
sm(2,2)=3.81277;
259
sm(2,3)=4.98668e+07;
260
sm(2,4)=-0.0697533;
261
sm(2,5)=12.8555;
262
sm(2,6)=-2.01124;
263
sm(3,3)=6.52498e+08;
264
sm(3,4)=3.87491;
265
sm(3,5)=365.965;
266
sm(3,6)=93.3686;
267
sm(4,4)=7.77672e-05;
268
sm(4,5)=0.0032134;
269
sm(4,6)=0.00194407;
270
sm(5,5)=0.132845;
271
sm(5,6)=0.0803294;
272
sm(6,6)=0.0485992;
273
274
*/
275
276
/**** a group of near singular (nonsingular) matrices
277
int n=5;
278
HepSymMatrix sm(n); // nxn Hilbert Matrix
279
for(i=1;i<=n;i++)
280
for(k=i;k<=n;k++)
281
sm(i,k)=1./(i+k-1);
282
*/
283
return
0;
284
}
// end of main
285
286
287
Generated on Mon May 6 2013 04:04:11 for CLHEP by
1.8.1.2