CLHEP 2.0.4.7 Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // testBug7328.cc 00002 // 00003 // The bug reported a memory leak in inverting symmetric matrices (of size 00004 // greater than 6x6). 00005 // 00006 // This test verifies that the leak is no longer present, and checks for 00007 // correctness. There is a break in method at N=25, so both sides are examined. 00008 // 00009 // A similar (though unreported) leak was present in the Determinant method; 00010 // since this was also fixed, this test tests for correctness of determinant as 00011 // well. 00012 // 00013 00014 #include <iostream> 00015 //#include <iomanip> 00016 #include <cmath> 00017 #include <stdlib.h> 00018 00019 #include "CLHEP/Matrix/Matrix.h" 00020 #include "CLHEP/Matrix/SymMatrix.h" 00021 00022 int test_inversion (int N) { 00023 00024 int i,j; 00025 CLHEP::HepSymMatrix S(N,0); 00026 for(i=1;i<=N;++i) { 00027 for(j=1;j<=N;++j) { 00028 if(i<=j) { 00029 S (i,j) = (10.0*i+j)/10; 00030 } 00031 } 00032 } 00033 CLHEP::HepSymMatrix SS(N,0); 00034 SS = S; 00035 int ierr = 0; 00036 SS.invert(ierr); 00037 if (ierr) { 00038 std::cout<<"SS.invert failed!!!! N = " << N << 00039 " ierr = "<< ierr <<std::endl; 00040 return 2+100*N; 00041 } 00042 CLHEP::HepMatrix SI(N,N,0); 00043 CLHEP::HepMatrix MS(N,N,0); 00044 CLHEP::HepMatrix MSS(N,N,0); 00045 MS = S; 00046 MSS = SS; 00047 SI = MSS*MS; 00048 for(i=1;i<=N;++i) { 00049 for(j=1;j<=N;++j) { 00050 if(i!=j) { 00051 if (fabs(SI(i,j)) > 1.0e-6) { 00052 std::cout<<"SS.invert incorrect N = " << N << 00053 " error = "<< fabs(SI(i,j)) <<std::endl; 00054 return 3+100*N; 00055 } 00056 } 00057 if(i==j) { 00058 if (fabs(1-SI(i,j)) > 1.0e-6) { 00059 std::cout<<"SS.invert incorrect N = " << N << 00060 " error = "<< fabs(1-SI(i,j)) <<std::endl; 00061 return 4+100*N; 00062 } 00063 } 00064 } 00065 } 00066 #define DET_ALSO 00067 #ifdef DET_ALSO 00068 double detS = S.determinant(); 00069 // std::cout<<"Determinant N = " << N << 00070 // " = " << detS <<std::endl; 00071 double detSS = SS.determinant(); 00072 // std::cout<<"Determinant Inverse N = " << N << 00073 // " = " << detSS <<std::endl; 00074 if (fabs((detS-1.0/detSS)/detS) > 1.0e-6) { 00075 std::cout<<"Determinant incorrect N = " << N << 00076 " error = " << fabs((detS-1.0/detSS)/detS) <<std::endl; 00077 return 5+100*N; 00078 } 00079 #endif 00080 return 0; 00081 } 00082 00083 void heapAddresses ( double * &hNew, 00084 double * &hMalloc, 00085 double * &hNew10000, 00086 double * &hMalloc80000 ) { 00087 hNew = new double; 00088 hMalloc = (double*) malloc(sizeof(double)); 00089 hNew10000 = new double[10000]; 00090 hMalloc80000 = (double*) malloc(10000*sizeof(double)); 00091 // std::cout << std::hex << hNew << " " << hMalloc<< " " 00092 // << hNew10000 << " " << hMalloc80000 << std::endl; 00093 free (hMalloc80000); 00094 delete[] hNew10000; 00095 free (hMalloc); 00096 delete hNew; 00097 } 00098 00099 int checkHeap ( double * &hNew, 00100 double * &hMalloc, 00101 double * &hNew10000, 00102 double * &hMalloc80000, 00103 double * &xhNew, 00104 double * &xhMalloc, 00105 double * &xhNew10000, 00106 double * &xhMalloc80000 ) { 00107 int ret = 0; 00108 if (hNew != xhNew) { 00109 std::cout<< "Leak:\n" 00110 << "xhNew - hNew = " << xhNew - hNew << "\n"; 00111 ret |= 1; 00112 } 00113 if (hMalloc != xhMalloc) { 00114 std::cout<< "Leak:\n" 00115 << "xhMalloc - hMalloc = " << xhMalloc - hMalloc << "\n"; 00116 ret |= 2; 00117 } 00118 if (hNew10000 != xhNew10000) { 00119 std::cout<< "Leak:\n" 00120 << "xhNew10000 - hNew10000 = " << xhNew10000 - hNew10000 << "\n"; 00121 ret |= 4; 00122 } 00123 if (hMalloc80000 != xhMalloc80000) { 00124 std::cout<< "Leak:\n" 00125 << "xhMalloc80000 - hMalloc80000 = " << xhMalloc80000 -hMalloc80000 00126 << "\n"; 00127 ret |= 8; 00128 } 00129 return ret; 00130 } 00131 00132 int main(int, char **) { 00133 int ret=0; 00134 int rhp; 00135 int i,j; 00136 for ( i = 1; i <= 50; i++) { 00137 ret = test_inversion(i); 00138 if (ret) return ret; 00139 } 00140 double *hNew, *hMalloc, *hNew10000, *hMalloc80000; 00141 double *xhNew, *xhMalloc, *xhNew10000, *xhMalloc80000; 00142 00143 for (int count=0; count < 2; ++count) { 00144 00145 int n1 = 400; 00146 int n2 = 25; 00147 heapAddresses ( hNew, hMalloc, hNew10000, hMalloc80000 ); 00148 for (i=0; i<n1; i++) { 00149 for (j=1; j <= n2; j++) { 00150 ret = test_inversion(j); 00151 if (ret) return ret; 00152 } 00153 } 00154 heapAddresses ( xhNew, xhMalloc, xhNew10000, xhMalloc80000 ); 00155 rhp = 0; 00156 if(count != 0) rhp = checkHeap ( hNew, hMalloc, hNew10000, hMalloc80000, 00157 xhNew, xhMalloc, xhNew10000, xhMalloc80000 ); 00158 if (rhp) std::cout << "Above Leak is after " << n1*n2 << " test inversions\n"; 00159 ret |= rhp; 00160 00161 heapAddresses ( hNew, hMalloc, hNew10000, hMalloc80000 ); 00162 for (i=0; i<2; i++) { 00163 for (j=1; j < 20; j++) { 00164 rhp = test_inversion(25+2*j); 00165 if (rhp) return rhp; 00166 } 00167 } 00168 heapAddresses ( xhNew, xhMalloc, xhNew10000, xhMalloc80000 ); 00169 rhp = 0; 00170 if(count != 0) rhp = checkHeap ( hNew, hMalloc, hNew10000, hMalloc80000, 00171 xhNew, xhMalloc, xhNew10000, xhMalloc80000 ); 00172 if (rhp) std::cout << "Leak after big inversions\n"; 00173 ret |= rhp; 00174 } 00175 return ret; 00176 } 00177