CLHEP VERSION Reference Documentation
CLHEP Home Page
CLHEP Documentation
CLHEP Bug Reports
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
B
Definition:
excDblThrow.cc:8
main
int main()
Definition:
testInversion.cc:23
D
Definition:
excDblThrow.cc:17
CLHEP::detail::n
n
Definition:
Ranlux64Engine.cc:85
CLHEP::HepMatrix
Definition:
Matrix/CLHEP/Matrix/Matrix.h:209
CLHEP
Definition:
ClhepVersion.h:13
CLHEP::HepSymMatrix
Definition:
Matrix/CLHEP/Matrix/SymMatrix.h:89
j
long j
Definition:
JamesRandomSeeding.txt:28
CLHEP::HepMatrix::inverse
HepMatrix inverse(int &ierr) const
Generated by
1.8.17