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 
void invert(int &ifail)
Definition: SymMatrix.cc:845
double determinant() const
Definition: SymMatrix.cc:943
#define double(obj)
Definition: excDblThrow.cc:32
int main(int, char **)
Definition: testBug7328.cc:132
int test_inversion(int N)
Definition: testBug7328.cc:22
void heapAddresses(double *&hNew, double *&hMalloc, double *&hNew10000, double *&hMalloc80000)
Definition: testBug7328.cc:83
int checkHeap(double *&hNew, double *&hMalloc, double *&hNew10000, double *&hMalloc80000, double *&xhNew, double *&xhMalloc, double *&xhNew10000, double *&xhMalloc80000)
Definition: testBug7328.cc:99