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

testBug6181.cc
Go to the documentation of this file.
1 // testBug6181.cc
2 //
3 // The bug reported an erroneous result in inverting a 7x7 matrix.
4 //
5 // 7x7 represents the break point between the low-n methods and
6 // the general-size method. This test program tests the case triggering
7 // the bug report, and also inversion for matrices (both symmetric and
8 // general) of sizes 1 - 9.
9 
10 #include <iostream>
11 #include <math.h>
12 
13 #include "CLHEP/Matrix/Matrix.h"
14 #include "CLHEP/Matrix/SymMatrix.h"
15 
16 int test_inversion (int N) {
17 
18  int i;
19  int j;
20  CLHEP::HepMatrix M(N,N,0);
21  CLHEP::HepSymMatrix S(N,0);
22  for(i=1;i<=N;++i) {
23  for(j=1;j<=N;++j) {
24  if(i<=j) {
25  S (i,j) = 10*i+j;
26  M (i,j) = 10*i+j;
27  M (j,i) = M(i,j) + .1 * (i-j)*(i-j);
28  }
29  }
30  }
31  CLHEP::HepMatrix MM(N,N,0);
32  CLHEP::HepSymMatrix SS(N,0);
33  MM = M;
34  SS = S;
35  int ierr = 0;
36  MM.invert(ierr);
37  if (ierr) {
38  std::cout<<"MM.invert failed!!!! N = " << N <<
39  " ierr = "<< ierr <<std::endl;
40  return 1+100*N;
41  }
42  ierr = 0;
43  SS.invert(ierr);
44  if (ierr) {
45  std::cout<<"SS.invert failed!!!! N = " << N <<
46  " ierr = "<< ierr <<std::endl;
47  return 2+100*N;
48  }
49  CLHEP::HepMatrix MI(N,N,0);
50  CLHEP::HepMatrix SI(N,N,0);
51  CLHEP::HepMatrix MS(N,N,0);
52  CLHEP::HepMatrix MSS(N,N,0);
53  MI = MM*M;
54  MS = S;
55  MSS = SS;
56  SI = MSS*MS;
57  for(i=1;i<=N;++i) {
58  for(j=1;j<=N;++j) {
59  if(i!=j) {
60  if (fabs(MI(i,j)) > 1.0e-12) {
61  std::cout<<"MM.invert incorrect N = " << N <<
62  " error = "<< fabs(MI(i,j)) <<std::endl;
63  return 3+100*N;
64  }
65  if (fabs(SI(i,j)) > 1.0e-12) {
66  std::cout<<"SS.invert incorrect N = " << N <<
67  " error = "<< fabs(SI(i,j)) <<std::endl;
68  return 4+100*N;
69  }
70  }
71  if(i==j) {
72  if (fabs(1-MI(i,j)) > 1.0e-12) {
73  std::cout<<"MM.invert incorrect N = " << N <<
74  " error = "<< fabs(1-MI(i,j)) <<std::endl;
75  return 3+100*N;
76  }
77  if (fabs(1-SI(i,j)) > 1.0e-12) {
78  std::cout<<"SS.invert incorrect N = " << N <<
79  " error = "<< fabs(1-SI(i,j)) <<std::endl;
80  return 4+100*N;
81  }
82  }
83  }
84  }
85  // Finally, the case that actually (for N=7) triggered bug report 6181:
86  M = S;
87  MM = M;
88  MM.invert(ierr);
89  if (ierr) {
90  std::cout<<"MM.invert for symmetric case failed!!!! N = " << N <<
91  " ierr = "<< ierr <<std::endl;
92  return 5+100*N;
93  }
94  MI = MM*M;
95  for(i=1;i<=N;++i) {
96  for(j=1;j<=N;++j) {
97  if(i!=j) {
98  if (fabs(MI(i,j)) > 1.0e-12) {
99  std::cout<<"MM.invert incorrect N = " << N <<
100  " error = "<< fabs(MI(i,j)) <<std::endl;
101  return 6+100*N;
102  }
103  }
104  if(i==j) {
105  if (fabs(1-MI(i,j)) > 1.0e-12) {
106  std::cout<<"MM.invert incorrect N = " << N <<
107  " error = "<< fabs(1-MI(i,j)) <<std::endl;
108  return 7+100*N;
109  }
110  }
111  }
112  }
113  return 0;
114 }
115 
116 
117 int main(int, char **) {
118  int ret=0;
119  for (int i=1; i<10; i++) {
120  ret = test_inversion(i);
121  if (ret) return ret;
122  }
123  return ret;
124 }
125 
CLHEP::HepMatrix::invert
virtual void invert(int &ierr)
Definition: Matrix.cc:707
main
int main(int, char **)
Definition: testBug6181.cc:117
CLHEP::HepMatrix
Definition: Matrix/CLHEP/Matrix/Matrix.h:209
N
the goal is to keep the overall false rejection probability down at the to level For each validated we discuss which of course is by necessity relative timing We take the time for a single random via one of the fastest good and at any rate the ratios will vary by around depending on the processor and memory configuration used A timing for a distribution of units would mean no time used beyond the uniform random Summary Distribution Validated Validation Rejected Past N RandGauss N
Definition: validation.doc:48
CLHEP::HepSymMatrix
Definition: Matrix/CLHEP/Matrix/SymMatrix.h:89
j
long j
Definition: JamesRandomSeeding.txt:28
i
long i
Definition: JamesRandomSeeding.txt:27
CLHEP::HepSymMatrix::invert
void invert(int &ifail)
Definition: SymMatrix.cc:845
test_inversion
int test_inversion(int N)
Definition: testBug6181.cc:16