CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

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 
HepMatrix inverse(int &ierr) const
Definition: excDblThrow.cc:8
Definition: excDblThrow.cc:17
int main()