Commit 8b9d504b authored by kraus's avatar kraus
Browse files

properly allocate and release memory

parent 559b6f79
......@@ -83,14 +83,18 @@ inline BMultipoleField::Pair BMultipoleField::Pair::operator-() const
// ------------------------------------------------------------------------
BMultipoleField::BMultipoleField():
pairs(0),
pairs(NULL),
itsOrder(0)
{}
BMultipoleField::BMultipoleField(const BMultipoleField &rhs):
pairs(new Pair[rhs.itsOrder]),
pairs(NULL),
itsOrder(rhs.itsOrder) {
if (itsOrder > 0)
pairs = new Pair[itsOrder];
for(int i = 0; i < itsOrder; i++) {
pairs[i] = rhs.pairs[i];
}
......@@ -104,9 +108,10 @@ BMultipoleField::~BMultipoleField() {
BMultipoleField &BMultipoleField::operator=(const BMultipoleField &rhs) {
if(&rhs != this) {
Pair *temp = new Pair[rhs.itsOrder];
if(pairs != 0) delete [] pairs;
pairs = temp;
if(pairs != NULL) delete[] pairs;
if (rhs.itsOrder > 0)
pairs = new Pair[rhs.itsOrder];
itsOrder = rhs.itsOrder;
for(int i = 0; i < itsOrder; i++) {
......
......@@ -105,6 +105,8 @@ CommMPI::CommMPI(int& argc , char**& argv, int procs, bool mpiinit, MPI_Comm mpi
// restore original executable name without absolute path
strcpy(argv[0],execname);
free(execname);
// duplicate the MPI_COMM_WORLD communicator, so that we can use
// a communicator that will not conflict with other users of MPI_COMM_WORLD
MPI_Comm_dup(mpicomm, &communicator);
......
......@@ -55,6 +55,8 @@ extern Inform *gmsg;
#define DISTDBG1
#define noDISTDBG2
#define SMALLESTCUTOFF 1e-12
SymTenzor<double, 6> getUnit6x6() {
SymTenzor<double, 6> unit6x6;
for (unsigned int i = 0; i < 6u; ++ i) {
......@@ -2387,8 +2389,8 @@ void Distribution::GenerateFlattopT(size_t numberOfParticles) {
}
}
if (randGenStandard)
gsl_rng_free(randGenStandard);
gsl_rng_free(randGenStandard);
if (quasiRandGen2D != NULL)
gsl_qrng_free(quasiRandGen2D);
......@@ -2474,14 +2476,13 @@ void Distribution::GenerateFlattopZ(size_t numberOfParticles) {
pzDist_m.push_back(pz);
}
}
if (randGenStandard)
gsl_rng_free(randGenStandard);
if (quasiRandGen1D != NULL)
gsl_qrng_free(quasiRandGen1D);
gsl_rng_free(randGenStandard);
if (quasiRandGen2D != NULL)
if (quasiRandGen1D != NULL) {
gsl_qrng_free(quasiRandGen1D);
gsl_qrng_free(quasiRandGen2D);
}
}
void Distribution::GenerateGaussZ(size_t numberOfParticles) {
......@@ -2489,15 +2490,15 @@ void Distribution::GenerateGaussZ(size_t numberOfParticles) {
gsl_rng_env_setup();
gsl_rng *randGen = gsl_rng_alloc(gsl_rng_default);
gsl_matrix * m = gsl_matrix_alloc (6, 6);
gsl_matrix * corMat = gsl_matrix_alloc (6, 6);
gsl_vector * rx = gsl_vector_alloc(6);
gsl_vector * ry = gsl_vector_alloc(6);
for (unsigned int i = 0; i < 6; ++ i) {
gsl_matrix_set(m, i, i, correlationMatrix_m(i, i));
gsl_matrix_set(corMat, i, i, correlationMatrix_m(i, i));
for (unsigned int j = 0; j < i; ++ j) {
gsl_matrix_set(m, i, j, correlationMatrix_m(i, j));
gsl_matrix_set(m, j, i, correlationMatrix_m(i, j));
gsl_matrix_set(corMat, i, j, correlationMatrix_m(i, j));
gsl_matrix_set(corMat, j, i, correlationMatrix_m(i, j));
}
}
......@@ -2508,15 +2509,15 @@ void Distribution::GenerateGaussZ(size_t numberOfParticles) {
for (int j = 0; j < 6; j++) {
if (j==0)
*gmsg << "* r(" << std::setprecision(1) << i << ", : ) = "
<< std::setprecision(4) << std::setw(10) << gsl_matrix_get (m, i, j);
<< std::setprecision(4) << std::setw(10) << gsl_matrix_get (corMat, i, j);
else
*gmsg << " & " << std::setprecision(4) << std::setw(10) << gsl_matrix_get (m, i, j);
*gmsg << " & " << std::setprecision(4) << std::setw(10) << gsl_matrix_get (corMat, i, j);
}
*gmsg << " \\\\" << endl;
}
#endif
int errcode = gsl_linalg_cholesky_decomp(m);
int errcode = gsl_linalg_cholesky_decomp(corMat);
if (errcode == GSL_EDOM) {
throw OpalException("Distribution::GenerateGaussZ",
......@@ -2525,7 +2526,7 @@ void Distribution::GenerateGaussZ(size_t numberOfParticles) {
// so make the upper part zero.
for (int i = 0; i < 5; ++ i) {
for (int j = i+1; j < 6 ; ++ j) {
gsl_matrix_set (m, i, j, 0.0);
gsl_matrix_set (corMat, i, j, 0.0);
}
}
......@@ -2536,9 +2537,9 @@ void Distribution::GenerateGaussZ(size_t numberOfParticles) {
for (int j = 0; j < 6; j++) {
if (j==0)
*gmsg << "* r(" << std::setprecision(1) << i << ", : ) = "
<< std::setprecision(4) << std::setw(10) << gsl_matrix_get (m, i, j);
<< std::setprecision(4) << std::setw(10) << gsl_matrix_get (corMat, i, j);
else
*gmsg << " & " << std::setprecision(4) << std::setw(10) << gsl_matrix_get (m, i, j);
*gmsg << " & " << std::setprecision(4) << std::setw(10) << gsl_matrix_get (corMat, i, j);
}
*gmsg << " \\\\" << endl;
}
......@@ -2579,7 +2580,7 @@ void Distribution::GenerateGaussZ(size_t numberOfParticles) {
saveProcessor = 0;
if (Ippl::myNode() == saveProcessor) {
if (gsl_blas_dgemv(CblasNoTrans,1.0,m,rx,0.0,ry)) {
if (gsl_blas_dgemv(CblasNoTrans,1.0,corMat,rx,0.0,ry)) {
throw OpalException("Distribution::GenerateGaussZ",
"oops... something wrong with GSL matvec");
}
......@@ -2592,15 +2593,10 @@ void Distribution::GenerateGaussZ(size_t numberOfParticles) {
}
}
/*
//std::for_each(v.rbegin(), v.rend(), [&](int n) { sum_of_elements += n; });
double pxm = std::accumulate(pxDist_m.begin(), pxDist_m.end(), 0.0);
double pym = std::accumulate(pyDist_m.begin(), pyDist_m.end(), 0.0);
double pzm = std::accumulate(pzDist_m.begin(), pzDist_m.end(), 0.0);
std::cout << "bega= " << std::sqrt(pxm*pxm + pym*pym + pzm*pzm) << std::endl;
*/
if (randGen)
gsl_rng_free(randGen);
gsl_rng_free(randGen);
gsl_vector_free(rx);
gsl_vector_free(ry);
gsl_matrix_free(corMat);
}
void Distribution::GenerateLongFlattopT(size_t numberOfParticles) {
......@@ -2762,15 +2758,12 @@ void Distribution::GenerateLongFlattopT(size_t numberOfParticles) {
}
gsl_rng_free(randGenGSL);
gsl_rng_free(randGen);
if (randGen)
gsl_rng_free(randGen);
if (quasiRandGen1D != NULL)
if (quasiRandGen1D != NULL) {
gsl_qrng_free(quasiRandGen1D);
if (quasiRandGen2D != NULL)
gsl_qrng_free(quasiRandGen2D);
}
}
void Distribution::GenerateTransverseGauss(size_t numberOfParticles) {
......@@ -2778,19 +2771,19 @@ void Distribution::GenerateTransverseGauss(size_t numberOfParticles) {
// Generate coordinates.
gsl_rng_env_setup();
gsl_rng *randGen = gsl_rng_alloc(gsl_rng_default);
gsl_matrix * m = gsl_matrix_alloc (4, 4);
gsl_matrix * corMat = gsl_matrix_alloc (4, 4);
gsl_vector * rx = gsl_vector_alloc(4);
gsl_vector * ry = gsl_vector_alloc(4);
for (unsigned int i = 0; i < 4; ++ i) {
gsl_matrix_set(m, i, i, correlationMatrix_m(i, i));
gsl_matrix_set(corMat, i, i, correlationMatrix_m(i, i));
for (unsigned int j = 0; j < i; ++ j) {
gsl_matrix_set(m, i, j, correlationMatrix_m(i, j));
gsl_matrix_set(m, j, i, correlationMatrix_m(i, j));
gsl_matrix_set(corMat, i, j, correlationMatrix_m(i, j));
gsl_matrix_set(corMat, j, i, correlationMatrix_m(i, j));
}
}
int errcode = gsl_linalg_cholesky_decomp(m);
int errcode = gsl_linalg_cholesky_decomp(corMat);
if (errcode == GSL_EDOM) {
throw OpalException("Distribution::GenerateTransverseGauss",
......@@ -2799,7 +2792,7 @@ void Distribution::GenerateTransverseGauss(size_t numberOfParticles) {
for (int i = 0; i < 3; ++ i) {
for (int j = i+1; j < 4 ; ++ j) {
gsl_matrix_set (m, i, j, 0.0);
gsl_matrix_set (corMat, i, j, 0.0);
}
}
......@@ -2819,28 +2812,28 @@ void Distribution::GenerateTransverseGauss(size_t numberOfParticles) {
gsl_vector_set(rx, 2, gsl_ran_gaussian (randGen,1.0));
gsl_vector_set(rx, 3, gsl_ran_gaussian (randGen,1.0));
gsl_blas_dgemv(CblasNoTrans, 1.0, m, rx, 0.0, ry);
gsl_blas_dgemv(CblasNoTrans, 1.0, corMat, rx, 0.0, ry);
x = sigmaR_m[0] * gsl_vector_get(ry, 0);
px = sigmaP_m[0] * gsl_vector_get(ry, 1);
y = sigmaR_m[1] * gsl_vector_get(ry, 2);
py = sigmaP_m[1] * gsl_vector_get(ry, 3);
bool xAndYOk = false;
if (cutoffR_m[0] <= 0.0 && cutoffR_m[1] <= 0.0)
if (cutoffR_m[0] < SMALLESTCUTOFF && cutoffR_m[1] < SMALLESTCUTOFF)
xAndYOk = true;
else if (cutoffR_m[0] <= 0.0)
else if (cutoffR_m[0] < SMALLESTCUTOFF)
xAndYOk = (std::abs(y) <= sigmaR_m[1] * cutoffR_m[1]);
else if (cutoffR_m[1] <= 0.0)
else if (cutoffR_m[1] < SMALLESTCUTOFF)
xAndYOk = (std::abs(x) <= sigmaR_m[1] * cutoffR_m[0]);
else
xAndYOk = (pow(x / (sigmaR_m[0] * cutoffR_m[0]), 2.0) + pow(y / (sigmaR_m[1] * cutoffR_m[1]), 2.0) <= 1.0);
bool pxAndPyOk = false;
if (cutoffP_m[0] <= 0.0 && cutoffP_m[1] <= 0.0)
if (cutoffP_m[0] < SMALLESTCUTOFF && cutoffP_m[1] < SMALLESTCUTOFF)
pxAndPyOk = true;
else if (cutoffP_m[0] <= 0.0)
else if (cutoffP_m[0] < SMALLESTCUTOFF)
pxAndPyOk = (std::abs(py) <= sigmaP_m[1] * cutoffP_m[1]);
else if (cutoffP_m[1] <= 0.0)
else if (cutoffP_m[1] < SMALLESTCUTOFF)
pxAndPyOk = (std::abs(px) <= sigmaP_m[1] * cutoffP_m[0]);
else
pxAndPyOk = (pow(px / (sigmaP_m[0] * cutoffP_m[0]), 2.0) + pow(py / (sigmaP_m[1] * cutoffP_m[1]), 2.0) <= 1.0);
......@@ -2861,8 +2854,11 @@ void Distribution::GenerateTransverseGauss(size_t numberOfParticles) {
pyDist_m.push_back(py);
}
}
if (randGen)
gsl_rng_free(randGen);
gsl_rng_free(randGen);
gsl_matrix_free(corMat);
gsl_vector_free(rx);
gsl_vector_free(ry);
}
void Distribution::InitializeBeam(PartBunch &beam) {
......@@ -4644,4 +4640,4 @@ void Distribution::WriteOutFileInjection() {
reduce(numberOfParticles, numberOfParticles, OpAddAssign());
}
}
}
\ No newline at end of file
}
Markdown is supported
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