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

testBug7328.cc
Go to the documentation of this file.
1 // testBug7328.cc
2 //
3 // The bug reported a memory leak in inverting symmetric matrices (of size
4 // greater than 6x6).
5 //
6 // This test verifies that the leak is no longer present, and checks for
7 // correctness. There is a break in method at N=25, so both sides are examined.
8 //
9 // A similar (though unreported) leak was present in the Determinant method;
10 // since this was also fixed, this test tests for correctness of determinant as
11 // well.
12 //
13 
14 #include <iostream>
15 //#include <iomanip>
16 #include <cmath>
17 #include <stdlib.h>
18 
19 #include "CLHEP/Matrix/Matrix.h"
20 #include "CLHEP/Matrix/SymMatrix.h"
21 
22 int test_inversion (int N) {
23 
24  int i,j;
25  CLHEP::HepSymMatrix S(N,0);
26  for(i=1;i<=N;++i) {
27  for(j=1;j<=N;++j) {
28  if(i<=j) {
29  S (i,j) = (10.0*i+j)/10;
30  }
31  }
32  }
33  CLHEP::HepSymMatrix SS(N,0);
34  SS = S;
35  int ierr = 0;
36  SS.invert(ierr);
37  if (ierr) {
38  std::cout<<"SS.invert failed!!!! N = " << N <<
39  " ierr = "<< ierr <<std::endl;
40  return 2+100*N;
41  }
42  CLHEP::HepMatrix SI(N,N,0);
43  CLHEP::HepMatrix MS(N,N,0);
44  CLHEP::HepMatrix MSS(N,N,0);
45  MS = S;
46  MSS = SS;
47  SI = MSS*MS;
48  for(i=1;i<=N;++i) {
49  for(j=1;j<=N;++j) {
50  if(i!=j) {
51  if (fabs(SI(i,j)) > 1.0e-6) {
52  std::cout<<"SS.invert incorrect N = " << N <<
53  " error = "<< fabs(SI(i,j)) <<std::endl;
54  return 3+100*N;
55  }
56  }
57  if(i==j) {
58  if (fabs(1-SI(i,j)) > 1.0e-6) {
59  std::cout<<"SS.invert incorrect N = " << N <<
60  " error = "<< fabs(1-SI(i,j)) <<std::endl;
61  return 4+100*N;
62  }
63  }
64  }
65  }
66 #define DET_ALSO
67 #ifdef DET_ALSO
68  double detS = S.determinant();
69 // std::cout<<"Determinant N = " << N <<
70 // " = " << detS <<std::endl;
71  double detSS = SS.determinant();
72 // std::cout<<"Determinant Inverse N = " << N <<
73 // " = " << detSS <<std::endl;
74  if (fabs((detS-1.0/detSS)/detS) > 1.0e-6) {
75  std::cout<<"Determinant incorrect N = " << N <<
76  " error = " << fabs((detS-1.0/detSS)/detS) <<std::endl;
77  return 5+100*N;
78  }
79 #endif
80  return 0;
81 }
82 
83 void heapAddresses ( double * &hNew,
84  double * &hMalloc,
85  double * &hNew10000,
86  double * &hMalloc80000 ) {
87  hNew = new double;
88  hMalloc = (double*) malloc(sizeof(double));
89  hNew10000 = new double[10000];
90  hMalloc80000 = (double*) malloc(10000*sizeof(double));
91 // std::cout << std::hex << hNew << " " << hMalloc<< " "
92 // << hNew10000 << " " << hMalloc80000 << std::endl;
93  free (hMalloc80000);
94  delete[] hNew10000;
95  free (hMalloc);
96  delete hNew;
97 }
98 
99 int checkHeap ( double * &hNew,
100  double * &hMalloc,
101  double * &hNew10000,
102  double * &hMalloc80000,
103  double * &xhNew,
104  double * &xhMalloc,
105  double * &xhNew10000,
106  double * &xhMalloc80000 ) {
107  int ret = 0;
108  if (hNew != xhNew) {
109  std::cout<< "Leak:\n"
110  << "xhNew - hNew = " << xhNew - hNew << "\n";
111  ret |= 1;
112  }
113  if (hMalloc != xhMalloc) {
114  std::cout<< "Leak:\n"
115  << "xhMalloc - hMalloc = " << xhMalloc - hMalloc << "\n";
116  ret |= 2;
117  }
118  if (hNew10000 != xhNew10000) {
119  std::cout<< "Leak:\n"
120  << "xhNew10000 - hNew10000 = " << xhNew10000 - hNew10000 << "\n";
121  ret |= 4;
122  }
123  if (hMalloc80000 != xhMalloc80000) {
124  std::cout<< "Leak:\n"
125  << "xhMalloc80000 - hMalloc80000 = " << xhMalloc80000 -hMalloc80000
126  << "\n";
127  ret |= 8;
128  }
129  return ret;
130 }
131 
132 int main(int, char **) {
133  int ret=0;
134  int rhp;
135  int i,j;
136  for ( i = 1; i <= 50; i++) {
137  ret = test_inversion(i);
138  if (ret) return ret;
139  }
140  double *hNew, *hMalloc, *hNew10000, *hMalloc80000;
141  double *xhNew, *xhMalloc, *xhNew10000, *xhMalloc80000;
142 
143  for (int count=0; count < 2; ++count) {
144 
145  int n1 = 400;
146  int n2 = 25;
147  heapAddresses ( hNew, hMalloc, hNew10000, hMalloc80000 );
148  for (i=0; i<n1; i++) {
149  for (j=1; j <= n2; j++) {
150  ret = test_inversion(j);
151  if (ret) return ret;
152  }
153  }
154  heapAddresses ( xhNew, xhMalloc, xhNew10000, xhMalloc80000 );
155  rhp = 0;
156  if(count != 0) rhp = checkHeap ( hNew, hMalloc, hNew10000, hMalloc80000,
157  xhNew, xhMalloc, xhNew10000, xhMalloc80000 );
158  if (rhp) std::cout << "Above Leak is after " << n1*n2 << " test inversions\n";
159  ret |= rhp;
160 
161  heapAddresses ( hNew, hMalloc, hNew10000, hMalloc80000 );
162  for (i=0; i<2; i++) {
163  for (j=1; j < 20; j++) {
164  rhp = test_inversion(25+2*j);
165  if (rhp) return rhp;
166  }
167  }
168  heapAddresses ( xhNew, xhMalloc, xhNew10000, xhMalloc80000 );
169  rhp = 0;
170  if(count != 0) rhp = checkHeap ( hNew, hMalloc, hNew10000, hMalloc80000,
171  xhNew, xhMalloc, xhNew10000, xhMalloc80000 );
172  if (rhp) std::cout << "Leak after big inversions\n";
173  ret |= rhp;
174  }
175  return ret;
176 }
177 
double
#define double(obj)
Definition: excDblThrow.cc:32
heapAddresses
void heapAddresses(double *&hNew, double *&hMalloc, double *&hNew10000, double *&hMalloc80000)
Definition: testBug7328.cc:83
main
int main(int, char **)
Definition: testBug7328.cc:132
checkHeap
int checkHeap(double *&hNew, double *&hMalloc, double *&hNew10000, double *&hMalloc80000, double *&xhNew, double *&xhMalloc, double *&xhNew10000, double *&xhMalloc80000)
Definition: testBug7328.cc:99
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
CLHEP::HepSymMatrix::determinant
double determinant() const
Definition: SymMatrix.cc:943
i
long i
Definition: JamesRandomSeeding.txt:27
test_inversion
int test_inversion(int N)
Definition: testBug7328.cc:22
count
user code seldom needs to call this function directly ZMerrno count() Return the(integer) number of ZMthrow 'n exceptions ever recorded via ZMerrno.write()
CLHEP::HepSymMatrix::invert
void invert(int &ifail)
Definition: SymMatrix.cc:845