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 
void invert(int &ifail)
Definition: SymMatrix.cc:845
int test_inversion(int N)
Definition: testBug6181.cc:16
int main(int, char **)
Definition: testBug6181.cc:117
virtual void invert(int &ierr)
Definition: Matrix.cc:707