Commit 8ad527fc authored by Mark Fishcler's avatar Mark Fishcler
Browse files

Fix for memory leak described in bug 7328; test program for that as well.

parent 1d5ecc29
2005-03-15 Mark Fischler <mf@fnal.gov>
* SymMatrix.cc, testBug7328.cc
Repair of bug 7328, a memory leak encountered when inverting symmetric
matrices above the size of 6x6.
2005-02-25 Mark Fischler <mf@fnal.gov>
* Matrix.cc, testBug6181.cc
Repair of bug 6181, a serious error in inverting matrices above the
size of 6x6.
2005-02-14 Mark Fischler <mf@fnal.gov>
* SymMatrixInvert.cc
......
// -*- C++ -*-
// $Id: SymMatrix.cc,v 1.3.2.5 2005/02/01 20:21:24 garren Exp $
// $Id: SymMatrix.cc,v 1.3.2.6 2005/03/15 20:39:06 fischler Exp $
// ---------------------------------------------------------------------------
//
// This file is a part of the CLHEP - a Class Library for High Energy Physics.
......@@ -952,14 +952,10 @@ void HepSymMatrix::invert(int &ifail) {
}
double HepSymMatrix::determinant() const {
static int max_array = 20;
static int *ir = new int [max_array+1];
if (nrow > max_array) {
delete [] ir;
max_array = nrow;
ir = new int [max_array+1];
}
static const int max_array = 20;
static std::vector<int> ir_vec (max_array);
ir_vec.resize(nrow);
int * ir = &ir_vec[0];
double det;
HepMatrix mt(*this);
......@@ -986,18 +982,32 @@ void HepSymMatrix::invertBunchKaufman(int &ifail) {
int i, j, k, s;
int pivrow;
int max_array = 25;
static int * piv = new int[max_array];
// used to store details of exchanges
// Establish the two working-space arrays needed: x and piv are
// used as pointers to arrays of doubles and ints respectively, each
// of length nrow. We do not want to reallocate each time through
// unless the size needs to grow. We do not want to leak memory, even
// by having a new without a delete that is only done once.
static const int max_array = 25;
#ifdef DISABLE_ALLOC
std::vector<double >* xvec = new std::vector<double>(max_array);
static std::vector<double> xvec (max_array);
static std::vector<int> pivv (max_array);
typedef std::vector<int>::iterator pivIter;
#else
std::vector<double,Alloc<double,25> >* xvec = new std::vector<double,Alloc<double,25> >(max_array);
#endif
static mIter x = xvec->begin();
// helper storage, needs to have at least size nrow,
// which will be less than or equal 6 most of the time
static std::vector<double,Alloc<double,25> > xvec (max_array);
static std::vector<int, Alloc<int, 25> > pivv (max_array);
typedef std::vector<int,Alloc<int,25> >::iterator pivIter;
#endif
xvec.resize(nrow);
pivv.resize(nrow);
// Note - resize does nothing if the size is already larger than nrow.
// Note - the data elements in a vector are guaranteed to be contiguous,
// so x[i] and piv[i] are optimally fast.
mIter x = xvec.begin();
// x[i] is used as helper storage, needs to have at least size nrow.
pivIter piv = pivv.begin();
// piv[i] is used to store details of exchanges
double temp1, temp2;
HepMatrix::mIter ip, mjj, iq;
......@@ -1008,16 +1018,6 @@ void HepSymMatrix::invertBunchKaufman(int &ifail) {
// it is set to zero.
// this constant could be set to zero but then the algorithm
// doesn't neccessarily detect that a matrix is singular
if (nrow > max_array) // need more space than expected
{
max_array = nrow;
delete [] piv;
xvec->resize(nrow);
x = xvec->begin();
piv = new int[nrow];
}
for (i = 0; i < nrow; i++)
piv[i] = i+1;
......
......@@ -21,7 +21,7 @@ CXXFLAGS = @CXXFLAGS@ -DINSTALLATION_CHECK
# Identify executables needed during testing:
check_PROGRAMS = \
testMatrix testInversion \
testBug6176 testBug6181
testBug6176 testBug6181 testBug7328
check_SCRIPTS = \
testMatrix.sh testInversion.sh \
......@@ -30,7 +30,7 @@ check_SCRIPTS = \
# Identify test(s) to run when 'make check' is requested:
TESTS = \
testMatrix.sh testInversion.sh \
testBug6176.sh testBug6181
testBug6176.sh testBug6181 testBug7328
# Identify the test(s) for which failure is the intended outcome:
XFAIL_TESTS =
......@@ -40,6 +40,7 @@ testMatrix_SOURCES = testMatrix.cc
testInversion_SOURCES = testInversion.cc
testBug6176_SOURCES = testBug6176.cc
testBug6181_SOURCES = testBug6181.cc
testBug7328_SOURCES = testBug7328.cc
# Identify input data file(s) and prototype output file(s):
EXTRA_DIST = \
......
// testBug7328.cc
//
// The bug reported a memory leak in inverting symmetric matrices (of size
// greater than 6x6).
//
// This test verifies that the leak is no longer present, and checks for
// correctness. There is a break in method at N=25, so both sides are examined.
//
// A similar (though unreported) leak was present in the Determinant method;
// since this was also fixed, this test tests for correctness of determinant as
// well.
//
#include <iostream>
#include <math.h>
#include <malloc.h>
#include "CLHEP/Matrix/Matrix.h"
#include "CLHEP/Matrix/SymMatrix.h"
int test_inversion (int N) {
HepSymMatrix S(N,0);
for(int i=1;i<=N;++i) {
for(int j=1;j<=N;++j) {
if(i<=j) {
S (i,j) = (10.0*i+j)/10;
}
}
}
HepSymMatrix SS(N,0);
SS = S;
int ierr = 0;
SS.invert(ierr);
if (ierr) {
std::cout<<"SS.invert failed!!!! N = " << N <<
" ierr = "<< ierr <<std::endl;
return 2+100*N;
}
HepMatrix SI(N,N,0);
HepMatrix MS(N,N,0);
HepMatrix MSS(N,N,0);
MS = S;
MSS = SS;
SI = MSS*MS;
int i;
int j;
for(i=1;i<=N;++i) {
for(j=1;j<=N;++j) {
if(i!=j) {
if (fabs(SI(i,j)) > 1.0e-6) {
std::cout<<"SS.invert incorrect N = " << N <<
" error = "<< fabs(SI(i,j)) <<std::endl;
return 3+100*N;
}
}
if(i==j) {
if (fabs(1-SI(i,j)) > 1.0e-6) {
std::cout<<"SS.invert incorrect N = " << N <<
" error = "<< fabs(1-SI(i,j)) <<std::endl;
return 4+100*N;
}
}
}
}
#define DET_ALSO
#ifdef DET_ALSO
double detS = S.determinant();
// std::cout<<"Determinant N = " << N <<
// " = " << detS <<std::endl;
double detSS = SS.determinant();
// std::cout<<"Determinant Inverse N = " << N <<
// " = " << detSS <<std::endl;
if (fabs((detS-1.0/detSS)/detS) > 1.0e-6) {
std::cout<<"Determinant incorrect N = " << N <<
" error = " << fabs((detS-1.0/detSS)/detS) <<std::endl;
return 5+100*N;
}
#endif
return 0;
}
void heapAddresses ( double * &hNew,
double * &hMalloc,
double * &hNew10000,
double * &hMalloc80000 ) {
hNew = new double;
hMalloc = (double*) malloc(8);
hNew10000 = new double[10000];
hMalloc80000 = (double*) malloc(80000);
free (hMalloc80000);
delete[] hNew10000;
free (hMalloc);
delete hNew;
}
int checkHeap ( double * &hNew,
double * &hMalloc,
double * &hNew10000,
double * &hMalloc80000,
double * &xhNew,
double * &xhMalloc,
double * &xhNew10000,
double * &xhMalloc80000 ) {
int ret = 0;
if (hNew + 2000 < xhNew) {
std::cout<< "Leak:\n"
<< "xhNew - hNew = " << xhNew - hNew << "\n";
ret |= 1;
}
if (hMalloc + 2000 < xhMalloc) {
std::cout<< "Leak:\n"
<< "xhMalloc - hMalloc = " << xhMalloc - hMalloc << "\n";
ret |= 2;
}
if (hNew10000 + 2000 < xhNew10000) {
std::cout<< "Leak:\n"
<< "xhNew10000 - hNew10000 = " << xhNew10000 - hNew10000 << "\n";
ret |= 4;
}
if (hMalloc80000 + 2000 < xhMalloc80000) {
std::cout<< "Leak:\n"
<< "xhMalloc80000 - hMalloc80000 = " << xhMalloc80000 -hMalloc80000
<< "\n";
ret |= 8;
}
return ret;
}
int main(int, char **) {
int ret=0;
int rhp;
for (int i = 1; i < 50; i++) {
ret = test_inversion(i);
if (ret) return ret;
}
double *hNew, *hMalloc, *hNew10000, *hMalloc80000;
double *xhNew, *xhMalloc, *xhNew10000, *xhMalloc80000;
int n1 = 4000;
int n2 = 25;
heapAddresses ( hNew, hMalloc, hNew10000, hMalloc80000 );
for (int i=1; i<n1; i++) {
for (int j=1; j < n2; j++) {
ret = test_inversion(j);
if (ret) return ret;
}
}
heapAddresses ( xhNew, xhMalloc, xhNew10000, xhMalloc80000 );
rhp = checkHeap ( hNew, hMalloc, hNew10000, hMalloc80000,
xhNew, xhMalloc, xhNew10000, xhMalloc80000 );
if (rhp) std::cout << "Above Leak is after " << n1*n2 << " test inversions\n";
ret |= rhp;
heapAddresses ( hNew, hMalloc, hNew10000, hMalloc80000 );
for (int i=1; i<2; i++) {
for (int j=1; j < 20; j++) {
ret |= test_inversion(25+2*j);
if (ret) return ret;
}
}
heapAddresses ( xhNew, xhMalloc, xhNew10000, xhMalloc80000 );
rhp = checkHeap ( hNew, hMalloc, hNew10000, hMalloc80000,
xhNew, xhMalloc, xhNew10000, xhMalloc80000 );
if (rhp) std::cout << "Leak after big inversions\n";
ret |= rhp;
return ret;
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment