diff --git a/ippl/src/Utility/IpplMemoryUsage.cpp b/ippl/src/Utility/IpplMemoryUsage.cpp index 4520beb24b4390310cc028794a7b4559abce16df..4d97ff2ef5968b8d76df90f08831009cbdc293d3 100644 --- a/ippl/src/Utility/IpplMemoryUsage.cpp +++ b/ippl/src/Utility/IpplMemoryUsage.cpp @@ -40,9 +40,11 @@ IpplMemoryUsage::IpplMemoryUsage(Unit unit, bool reset) case Unit::BIT: conversion_factor_m = 8.0e3; unit_m = "bit"; + break; case Unit::B: conversion_factor_m = 1.0e3; unit_m = "B"; + break; case Unit::KB: conversion_factor_m = 1.0; unit_m = "kB"; diff --git a/optimizer/Expression/FromFile.h b/optimizer/Expression/FromFile.h index b7ebc10c91803f85bbc5bbde348915e2cadb571e..9d142352c8d51eac0b71c3cad8873101ac28441e 100644 --- a/optimizer/Expression/FromFile.h +++ b/optimizer/Expression/FromFile.h @@ -41,7 +41,7 @@ struct FromFile { try { readValues(); - } catch(OptPilotException e) { + } catch(const OptPilotException& e) { return boost::make_tuple(0.0, false); } diff --git a/src/Algorithms/CMakeLists.txt b/src/Algorithms/CMakeLists.txt index 54a1c7c43f93852b080d0443fef071b9314dbdc8..ebe21b7f4c87d1121ce38d73c3b214b5de43c5d2 100644 --- a/src/Algorithms/CMakeLists.txt +++ b/src/Algorithms/CMakeLists.txt @@ -19,7 +19,6 @@ set (_SRCS bet/EnvelopeSlice.cpp bet/EnvelopeBunch.cpp bet/profile.cpp - bet/math/sort.cpp bet/math/integrate.cpp bet/math/interpol.cpp bet/math/rk.cpp @@ -66,7 +65,6 @@ set (HDRS bet/math/rk.h bet/math/root.h bet/math/savgol.h - bet/math/sort.h bet/math/svdfit.h ) diff --git a/src/Algorithms/bet/EnvelopeBunch.cpp b/src/Algorithms/bet/EnvelopeBunch.cpp index 8fe9bf5aacaa4c21b25d6ee28b10d156635d5ec3..a7e067a074ab8b5f03c5fea4cba0b334de8d596e 100644 --- a/src/Algorithms/bet/EnvelopeBunch.cpp +++ b/src/Algorithms/bet/EnvelopeBunch.cpp @@ -13,7 +13,6 @@ #include <mpi.h> #include "Algorithms/bet/math/root.h" // root finding routines -#include "Algorithms/bet/math/sort.h" // sorting routines #include "Algorithms/bet/math/linfit.h" // linear fitting routines #include "Algorithms/bet/math/savgol.h" // savgol smoothing routine #include "Algorithms/bet/math/rk.h" // Runge-Kutta Integration @@ -721,16 +720,12 @@ void EnvelopeBunch::calcI() { already_called = 1; std::vector<double> z1(numSlices_m, 0.0); - std::vector<double> b(numSlices_m, 0.0); + std::vector<double> b (numSlices_m, 0.0); double bSum = 0.0; double dz2Sum = 0.0; int n1 = 0; int my_start = 0, my_end = 0; - for(int i = 0; i < numSlices_m; i++) { - z1[i] = 0.0; - b[i] = 0.0; - } for(int i = 0; i < numSlices_m; i++) { if(b_m[i] > 0.0) { if((unsigned int) i == mySliceStartOffset_m) @@ -759,8 +754,20 @@ void EnvelopeBunch::calcI() { double sigma_dz = sqrt(dz2Sum / (n1 - 1)); double beta = bSum / n1; - //sort z1 according to beta's - sort2(&z1.front(), &b.front(), n1); + //sort beta's according to z1 with temporary vector of pairs + std::vector<std::pair<double,double>> z1_b(numSlices_m); + for (int i=0; i<numSlices_m; i++) { + z1_b[i] = std::make_pair(z1[i],b[i]); + } + + std::sort(z1_b.begin(), z1_b.end(), + [](const std::pair<double,double> &left, const std::pair<double,double> &right) + {return left.first < right.first;}); + + for (int i=0; i<numSlices_m; i++) { + z1[i] = z1_b[i].first; + b [i] = z1_b[i].second; + } double q = Q_m > 0.0 ? Q_m / numSlices_m : Physics::q_e; double dz = 0.0; diff --git a/src/Algorithms/bet/math/sort.cpp b/src/Algorithms/bet/math/sort.cpp deleted file mode 100644 index 1a9c5aba62013144db5a5bc354ba2827eaff6214..0000000000000000000000000000000000000000 --- a/src/Algorithms/bet/math/sort.cpp +++ /dev/null @@ -1,387 +0,0 @@ -/* Sort.C - Sorting routines - - NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) - - Revision history - Date Description Programmer - ------------ -------------------------------------------- ------------ - 10/03/2006 Copied from "numerical recipies" Rene Bakker - - Last Revision: - $Id: sort.C 85 2007-05-03 15:55:00Z bakker $ -*/ - -#include <stdlib.h> -#include <stdio.h> -#include <math.h> - -#include "Algorithms/bet/BetError.h" -#include "Algorithms/bet/math/sort.h" - - -/* ivector - allocate an array of integers [0..n-1] and check memory, - reallocate old vector if applicable */ -static int *ivector(int n, int *old = NULL) { - int *b; - - b = (int *) malloc(sizeof(int) * n); - if(!b) { - //writeBetError(errModeAll,errFatal,"Insufficient memory malloc %d bytes (sort.C)",sizeof(int)*n); - writeBetError("Insufficient memory malloc in sort.C"); - } - - if(old) { - int - i, - iMax = sizeof(old) / sizeof(int); - - for(i = 0; (i < iMax) && (i < n); i++) b[i] = old[i]; - free(old); - } - return b; -} /* ivector */ - - -#define SWAP(a,b) temp=(a);(a)=(b);(b)=temp -#define ISWAP(ia,ib) itemp=(ia);(ia)=(ib);(ib)=itemp -#define PSWAP(ap,bp) tempp=(ap);(ap)=(bp);(bp)=tempp -#define M 7 - -/* sort() - Sorts an array arr[0..n-1] into ascending numerical order using the Quicksort - algorithm. n is input; arr is replaced on output by its sorted rearrangement. -*/ -void sortN(double *x, double **y, int n) { - int - nstack = 50; - - int - i, ir = n - 1, - j, k, l = 0; - int - *istack, jstack = 0; - double - a, temp, - *bp = NULL, *tempp; - - istack = ivector(nstack); - - for(;;) { - if(ir - l < M) { - for(j = l + 1; j <= ir; j++) { - a = x[j]; - if(y) bp = y[j]; - for(i = j - 1; i >= 0; i--) { - if(x[i] <= a) break; - x[i+1] = x[i]; - if(y) y[i+1] = y[i]; - } - x[i+1] = a; - if(y) y[i+1] = bp; - } - if(!jstack) { - free(istack); - return; - } - ir = istack[jstack]; - l = istack[jstack-1]; - jstack -= 2; - } else { - k = (l + ir) >> 1; - SWAP(x[k], x[l+1]); - if(y) {PSWAP(y[k], y[l+1]);} - if(x[l+1] > x[ir]) { - SWAP(x[l+1], x[ir]); - if(y) {PSWAP(y[l+1], y[ir]);} - } - if(x[l] > x[ir]) { - SWAP(x[l], x[ir]); - if(y) {PSWAP(y[l], y[ir]);} - } - if(x[l+1] > x[l]) { - SWAP(x[l+1], x[l]); - if(y) {PSWAP(y[l+1], y[l]);} - } - i = l + 1; - j = ir; - a = x[l]; - if(y) bp = y[l]; - for(;;) { - do i++; - while(x[i] < a); - do j--; - while(x[j] > a); - if(j < i) break; - SWAP(x[i], x[j]); - if(y) {PSWAP(y[i], y[j]);} - } - x[l] = x[j]; - x[j] = a; - if(y) { - y[l] = y[j]; - y[j] = bp; - } - jstack += 2; - if(jstack > nstack - 1) { // make stack bigger - nstack += 50; - istack = ivector(nstack, istack); - } - if(ir - i + 1 >= j - l) { - istack[jstack] = ir; - istack[jstack-1] = i; - ir = j - 1; - } else { - istack[jstack] = j - 1; - istack[jstack-1] = l; - l = i; - } - } - } -} /* sortN() */ - -void sort2(double *x, double *y, int n) { - int - nstack = 50; - int - i, ir = n - 1, - j, k, l = 0; - int - *istack, jstack = 0; - double - a, b = 0.0, temp; - - istack = ivector(nstack); - - for(;;) { - if(ir - l < M) { - for(j = l + 1; j <= ir; j++) { - a = x[j]; - b = y[j]; - for(i = j - 1; i >= 0; i--) { - if(x[i] <= a) break; - x[i+1] = x[i]; - y[i+1] = y[i]; - } - x[i+1] = a; - y[i+1] = b; - } - if(!jstack) { - free(istack); - return; - } - ir = istack[jstack]; - l = istack[jstack-1]; - jstack -= 2; - } else { - k = (l + ir) >> 1; - SWAP(x[k], x[l+1]); - SWAP(y[k], y[l+1]); - if(x[l+1] > x[ir]) { - SWAP(x[l+1], x[ir]); - SWAP(y[l+1], y[ir]); - } - if(x[l] > x[ir]) { - SWAP(x[l], x[ir]); - SWAP(y[l], y[ir]); - } - if(x[l+1] > x[l]) { - SWAP(x[l+1], x[l]); - SWAP(y[l+1], y[l]); - } - i = l + 1; - j = ir; - a = x[l]; - if(y) b = y[l]; - for(;;) { - do i++; - while(x[i] < a); - do j--; - while(x[j] > a); - if(j < i) break; - SWAP(x[i], x[j]); - SWAP(y[i], y[j]); - } - x[l] = x[j]; - x[j] = a; - y[l] = y[j]; - y[j] = b; - jstack += 2; - if(jstack > nstack - 1) { // make stack bigger - nstack += 50; - istack = ivector(nstack, istack); - } - if(ir - i + 1 >= j - l) { - istack[jstack] = ir; - istack[jstack-1] = i; - ir = j - 1; - } else { - istack[jstack] = j - 1; - istack[jstack-1] = l; - l = i; - } - } - } -} /* sort2() */ - -void isort2(double *x, int *y, int n) { - int - nstack = 50; - - int - i, ir = n - 1, - j, k, l = 0; - int - ib = -1, itemp, - *istack, jstack = 0; - double - a, temp; - - istack = ivector(nstack); - - for(;;) { - if(ir - l < M) { - for(j = l + 1; j <= ir; j++) { - a = x[j]; - if(y) ib = y[j]; - for(i = j - 1; i >= 0; i--) { - if(x[i] <= a) break; - x[i+1] = x[i]; - if(y) y[i+1] = y[i]; - } - x[i+1] = a; - if(y) y[i+1] = ib; - } - if(!jstack) { - free(istack); - return; - } - ir = istack[jstack]; - l = istack[jstack-1]; - jstack -= 2; - } else { - k = (l + ir) >> 1; - SWAP(x[k], x[l+1]); - if(y) {ISWAP(y[k], y[l+1]);} - if(x[l+1] > x[ir]) { - SWAP(x[l+1], x[ir]); - if(y) {ISWAP(y[l+1], y[ir]);} - } - if(x[l] > x[ir]) { - SWAP(x[l], x[ir]); - if(y) {ISWAP(y[l], y[ir]);} - } - if(x[l+1] > x[l]) { - SWAP(x[l+1], x[l]); - if(y) {ISWAP(y[l+1], y[l]);} - } - i = l + 1; - j = ir; - a = x[l]; - if(y) ib = y[l]; - for(;;) { - do i++; - while(x[i] < a); - do j--; - while(x[j] > a); - if(j < i) break; - SWAP(x[i], x[j]); - if(y) {ISWAP(y[i], y[j]);} - } - x[l] = x[j]; - x[j] = a; - if(y) { - y[l] = y[j]; - y[j] = ib; - } - jstack += 2; - if(jstack > nstack - 1) { // make stack bigger - nstack += 50; - istack = ivector(nstack, istack); - } - if(ir - i + 1 >= j - l) { - istack[jstack] = ir; - istack[jstack-1] = i; - ir = j - 1; - } else { - istack[jstack] = j - 1; - istack[jstack-1] = l; - l = i; - } - } - } -} /* isort2() */ - -void sort1(double *x, int n) { - int nstack = 50; - - int - i, ir = n - 1, - j, k, l = 0; - int jstack = 0, *istack; - double a, temp; - - istack = ivector(nstack); - for(;;) { - if(ir - l < M) { - for(j = l + 1; j <= ir; j++) { - a = x[j]; - for(i = j - 1; i >= 0; i--) { - if(x[i] <= a) break; - x[i+1] = x[i]; - } - x[i+1] = a; - } - if(jstack == 0) break; - ir = istack[jstack--]; - l = istack[jstack--]; - } else { - k = (l + ir) >> 1; - SWAP(x[k], x[l+1]); - if(x[l+1] > x[ir]) { - SWAP(x[l+1], x[ir]); - } - if(x[l] > x[ir]) { - SWAP(x[l], x[ir]); - } - if(x[l+1] > x[l]) { - SWAP(x[l+1], x[l]); - } - i = l + 1; - j = ir; - a = x[l]; - for(;;) { - do i++; - while(x[i] < a); - do j--; - while(x[j] > a); - if(j < i) break; - SWAP(x[i], x[j]); - } - x[l] = x[j]; - x[j] = a; - jstack += 2; - if(jstack > nstack - 1) { // make stack bigger - nstack += 50; - istack = ivector(nstack, istack); - } - if(ir - i + 1 >= j - l) { - istack[jstack] = ir; - istack[jstack-1] = i; - ir = j - 1; - } else { - istack[jstack] = j - 1; - istack[jstack-1] = l; - l = i; - } - } - } - free(istack); -} - -#undef M -#undef SWAP -#undef ISWAP -#undef PSWAP - diff --git a/src/Algorithms/bet/math/sort.h b/src/Algorithms/bet/math/sort.h deleted file mode 100644 index ad0b2b2c9abec4fd5e4e41cd1ae5640805dbd599..0000000000000000000000000000000000000000 --- a/src/Algorithms/bet/math/sort.h +++ /dev/null @@ -1,42 +0,0 @@ -/* sort.h - sorting routines - - Project: Beam Envelope Tracker (BET) - - Revision history - Date Description Programmer - ------------ -------------------------------------------- -------------- - 09-03-06 Created Rene Bakker - - Last Revision: - $Id: sort.h 85 2007-05-03 15:55:00Z bakker $ -*/ - - -#ifndef _SORT_DEF -#define _SORT_DEF - -/* sort() - Sorts an array arr[0..n-1] into ascending numerical order using the Quicksort - algorithm. n is input; arr is replaced on output by its sorted rearrangement. -*/ -void sort1( - double *, // array to sort - int); // number of elements in array - -void sort2( - double *, // array to sort - double *, // array to sort along - int); // number of elements in array - -void isort2( - double *, // array to sort - int *, // index array to sort along - int); // number of elements in array - -void sortN( - double *, // array to sort - double **, // matrix to sort with x - int); // number of elements in array - -#endif diff --git a/src/Algorithms/bet/profile.cpp b/src/Algorithms/bet/profile.cpp index 95c5c9e2f28f95c6cddf3cf316376e8d7aa266fc..614f535e59497405b8f7cbde1562edc23adec692 100644 --- a/src/Algorithms/bet/profile.cpp +++ b/src/Algorithms/bet/profile.cpp @@ -14,11 +14,11 @@ #include "Ippl.h" #include <iostream> +#include <algorithm> #include <cmath> #include <stdlib.h> #include <string.h> -#include "Algorithms/bet/math/sort.h" #include "Algorithms/bet/math/interpol.h" #include "Algorithms/bet/math/integrate.h" #include "Algorithms/bet/profile.h" @@ -42,25 +42,14 @@ static double f3(double x) { Profile::Profile(double v) { n = 0; - x = NULL; - y = NULL; - y2 = NULL; sf = 1.0; yMin = v; yMax = v; } -Profile::Profile(double *_x, double *_y, int _n) { - n = _n; - - x = (double *) malloc(sizeof(double) * n); - y = (double *) malloc(sizeof(double) * n); - if((x == NULL) || (y == NULL)) - std::cout << "Profile::Profile: Insuffient memory for Profile(n = " << n << " )" << std::endl; - memcpy(x, _x, sizeof(double)*n); - memcpy(y, _y, sizeof(double)*n); - +Profile::Profile(double *_x, double *_y, int _n) : + n(_n), x(_x, _x + _n), y(_y, _y + _n) { create(); } /* Profile::Profile() */ @@ -82,12 +71,9 @@ Profile::Profile(char *fname, double eps) { fclose(f); n = i; - x = (double *) malloc(sizeof(double) * n); - y = (double *) malloc(sizeof(double) * n); + x.resize(n); + y.resize(n); - if((x == NULL) || (y == NULL)) { - std::cout << "Profile::Profile: Insuffient memory for Profile(n = " << n << " )" << std::endl; - } // read all values f = fopen(fname, "r"); for(i = 0; i < n; i++) { @@ -117,12 +103,6 @@ Profile::Profile(char *fname, double eps) { create(); } /* Profile::Profile() */ -Profile::~Profile() { - if(x) free(x); - if(y) free(y); - if(y2) free(y2); -} /* Profile::~Profile() */ - void Profile::create() { int i, j, k; @@ -131,18 +111,30 @@ void Profile::create() { sf = 1.0; - y2 = (double *) malloc(sizeof(double) * n); + y2.resize(n); yMax = y[0]; yMin = yMax; - if(y2 == NULL) { - std::cout << "Profile::Profile: Insuffient memory for Profile y2 (n = " << n << " )" << std::endl; - } for(i = 0; i < n; i++) { if(y[i] < yMin) yMin = y[i]; if(y[i] > yMax) yMax = y[i]; } - sort2(x, y, n); + + //sort y according to x + std::vector<std::pair<double,double>> xy(n); + for (int i=0; i<n; i++) { + xy[i] = std::make_pair(x[i],y[i]); + } + + std::sort(xy.begin(), xy.end(), + [](const std::pair<double,double> &left, const std::pair<double,double> &right) + {return left.first < right.first;}); + + for (int i=0; i<n; i++) { + x[i] = xy[i].first; + y[i] = xy[i].second; + } + /* remove points with identical x-value take the average if the situation does occur */ @@ -163,21 +155,21 @@ void Profile::create() { } n = i; - spline(x, y, n, y2); + spline(&x[0], &y[0], n, &y2[0]); } /* Profile::create() */ double Profile::get(double xa, Interpol_type tp) { double val = 0.0; - if(x) { + if(x.empty()==false) { if(xa < x[0]) val = 0.0; else if(xa > x[n-1]) val = 0.0; else switch(tp) { case itype_lin : - lsplint(x, y, y2, n, xa, &val); + lsplint(&x[0], &y[0], &y2[0], n, xa, &val); break; default : - lsplint(x, y, y2, n, xa, &val); + lsplint(&x[0], &y[0], &y2[0], n, xa, &val); break; } } @@ -262,11 +254,11 @@ double Profile::max() { } double Profile::xMax() { - return (x ? x[n-1] : 0.0); + return ((x.empty()==false) ? x[n-1] : 0.0); } double Profile::xMin() { - return (x ? x[0] : 0.0); + return ((x.empty()==false) ? x[0] : 0.0); } double Profile::Leff() { @@ -274,7 +266,7 @@ double Profile::Leff() { ym = fabs((fabs(yMin) > fabs(yMax)) ? yMin : yMax); cProfile = this; - return (((x == NULL) || (x[n-1] == x[0]) || (ym == 0.0)) ? 0.0 : + return ((x.empty() || (x[n-1] == x[0]) || (ym == 0.0)) ? 0.0 : fabs(qromb(f1, x[0], x[n-1]) / ym)); } @@ -283,7 +275,7 @@ double Profile::Leff2() { ym = pow((fabs(yMin) > fabs(yMax)) ? yMin : yMax, 2); cProfile = this; - return (((x == NULL) || (x[n-1] == x[0]) || (ym == 0.0)) ? 0.0 : + return ((x.empty() || (x[n-1] == x[0]) || (ym == 0.0)) ? 0.0 : fabs(qromb(f2, x[0], x[n-1]) / ym)); } @@ -292,6 +284,6 @@ double Profile::Labs() { ym = fabs((fabs(yMin) > fabs(yMax)) ? yMin : yMax); cProfile = this; - return (((x == NULL) || (x[n-1] == x[0]) || (ym == 0.0)) ? 0.0 : + return ((x.empty() || (x[n-1] == x[0]) || (ym == 0.0)) ? 0.0 : fabs(qromb(f3, x[0], x[n-1]) / ym)); } \ No newline at end of file diff --git a/src/Algorithms/bet/profile.h b/src/Algorithms/bet/profile.h index 47e4f034755b2a3fdb676e20f2898c746439f347..79ecd304eb23757d39f098023fbbd30f76af9d4a 100644 --- a/src/Algorithms/bet/profile.h +++ b/src/Algorithms/bet/profile.h @@ -19,6 +19,8 @@ #include <stdio.h> +#include <vector> + enum Interpol_type { itype_spline, // spline interpolation itype_lin, // linearl interpolation @@ -30,8 +32,8 @@ class Profile { double yMax, // maximum of array yMin, // minimum of array - sf, // scaling factor (1.0 by default) - *x, *y, *y2; + sf; // scaling factor (1.0 by default) + std::vector<double> x, y, y2; private: void create(); // general creator routine @@ -45,7 +47,6 @@ public: Profile( // creator from file char *, // filename double = 0.0); // cutoff value - ~Profile(); void normalize(); // set max of profile to 1.0 void scale(double); // scale the amplitude diff --git a/src/Structure/BoundaryGeometry.cpp b/src/Structure/BoundaryGeometry.cpp index 0a2ce21f9a97b3b2da94b781169cd7fb935488a0..829302bc546650e7dd9672ec03fb17cfcc1a8f12 100644 --- a/src/Structure/BoundaryGeometry.cpp +++ b/src/Structure/BoundaryGeometry.cpp @@ -964,10 +964,12 @@ BoundaryGeometry::intersectLineTriangle ( if (r < 0.0) { // ray goes away from triangle return 0; // => no intersect } + break; case SEGMENT: if (r < 0 || 1.0 < r) { // intersection on extended return 0; // segment } + break; case LINE: break; }; diff --git a/src/Tables/ThreadAll.cpp b/src/Tables/ThreadAll.cpp index bc9f32a16657b8dd969175bcb40931ab06e7919e..337b9c0d55aa6c5b572671338b67e4625504f4a4 100644 --- a/src/Tables/ThreadAll.cpp +++ b/src/Tables/ThreadAll.cpp @@ -131,7 +131,7 @@ void ThreadAll::execute() { } else { ++iter; } - } catch(DomainError) { + } catch(const DomainError &) { // Domain error detected; // attempt to correct in the plane which has the larger deviation. FVector<double, 6> orbit = itsTracker->getOrbit(); diff --git a/src/Tables/ThreadBpm.cpp b/src/Tables/ThreadBpm.cpp index cff826efbd2e42ae6ef3c6dd69a8955c392be44d..a074f2ace53a25e240eb48705181db618f5c47ba 100644 --- a/src/Tables/ThreadBpm.cpp +++ b/src/Tables/ThreadBpm.cpp @@ -132,7 +132,7 @@ void ThreadBpm::execute() { // Go to next element. ++iter; } - } catch(DomainError) { + } catch(const DomainError &) { break; } }