various coordinate performance improvements

casa/Quanta/Euler.cc

@@ -38,14 +38,44 @@ namespace casacore { //# NAMESPACE CASACORE - BEGIN
 
 // Euler class
 
+// simplistic Vector(3) cache to reduce allocation overhead for temporaries
+#if defined(AIPS_CXX11) && !defined(__APPLE__)
+thread_local size_t Euler::available = 0;
+thread_local Euler::DataArrays Euler::arrays[50];
+#endif
+
+Euler::DataArrays Euler::get_arrays()
+{
+#if defined(AIPS_CXX11) && !defined(__APPLE__)
+  if (available > 0) {
+    return arrays[--available];
+  }
+#endif
+  return std::make_pair(new Vector<Double>(3), new Vector<Int>(3));
+}
+
+void Euler::return_arrays(Euler::DataArrays array)
+{
+#if defined(AIPS_CXX11) && !defined(__APPLE__)
+  if (available < sizeof(arrays) / sizeof(arrays[0]) &&
+		  array.first->size() == 3 && array.second->size() == 3 &&
+		  array.first->nrefs() == 1 && array.second->nrefs() == 1) {
+    arrays[available++] = array;
+    return;
+  }
+#endif
+  delete array.first;
+  delete array.second;
+}
+
 //# Constructors
-Euler::Euler() : euler(3), axes(3) {
+Euler::Euler() : data(get_arrays()), euler(*data.first), axes(*data.second) {
     euler = Double(0.0);
     indgen(axes,1,1);
 }
 
 Euler::Euler(const Euler &other) : 
-euler(3), axes(3) {
+data(get_arrays()), euler(*data.first), axes(*data.second) {
     euler = other.euler;
     axes = other.axes;
 }
@@ -59,7 +89,7 @@ Euler &Euler::operator=(const Euler &other) {
 }
 
 Euler::Euler(Double in0, Double in1, Double in2) : 
-euler(3), axes(3){
+data(get_arrays()), euler(*data.first), axes(*data.second) {
     euler(0) = in0;
     euler(1) = in1;
     euler(2) = in2;
@@ -69,7 +99,7 @@ euler(3), axes(3){
 
 Euler::Euler(Double in0, uInt ax0, Double in1, uInt ax1, Double in2,
 	     uInt ax2) : 
-euler(3), axes(3) {
+data(get_arrays()), euler(*data.first), axes(*data.second) {
     DebugAssert(ax0 <= 3 && ax1 <=3 && ax2 <=3, AipsError);
     euler(0) = in0;
     euler(1) = in1;
@@ -80,31 +110,31 @@ euler(3), axes(3) {
 }
 
 Euler::Euler(const Quantity &in0) :
-euler(3), axes(3) {
+data(get_arrays()), euler(*data.first), axes(*data.second) {
     euler(0) = Euler::makeRad(in0);
     euler(1) = 0;
     euler(2) = 0;
     indgen(axes,1,1);
 }
 
 Euler::Euler(const Quantity &in0, const Quantity &in1) :
-euler(3), axes(3) {
+data(get_arrays()), euler(*data.first), axes(*data.second) {
     euler(0) = Euler::makeRad(in0);
     euler(1) = Euler::makeRad(in1);
     euler(2) = 0;
     indgen(axes,1,1);
 }
 
 Euler::Euler(const Quantity &in0, const Quantity &in1, const Quantity &in2) :
-euler(3), axes(3) {
+data(get_arrays()), euler(*data.first), axes(*data.second) {
     euler(0) = Euler::makeRad(in0);
     euler(1) = Euler::makeRad(in1);
     euler(2) = Euler::makeRad(in2);
     indgen(axes,1,1);
 }
 
 Euler::Euler(const Quantity &in0, uInt ax0) :
-euler(3), axes(3) {
+data(get_arrays()), euler(*data.first), axes(*data.second) {
     DebugAssert(ax0 <= 3, AipsError);
     euler(0) = Euler::makeRad(in0);
     euler(1) = 0;
@@ -114,7 +144,7 @@ euler(3), axes(3) {
     axes(2) = 0;
 }
 Euler::Euler(const Quantity &in0, uInt ax0, const Quantity &in1, uInt ax1) :
-euler(3), axes(3) {
+data(get_arrays()), euler(*data.first), axes(*data.second) {
     DebugAssert(ax0 <= 3 && ax1 <=3, AipsError);
     euler(0) = Euler::makeRad(in0);
     euler(1) = Euler::makeRad(in1);
@@ -125,7 +155,7 @@ euler(3), axes(3) {
 }
 Euler::Euler(const Quantity &in0, uInt ax0, const Quantity &in1, uInt ax1,
 	     const Quantity &in2, uInt ax2) :
-euler(3), axes(3) {
+data(get_arrays()), euler(*data.first), axes(*data.second) {
     DebugAssert(ax0 <= 3 && ax1 <=3 && ax2 <=3, AipsError);
     euler(0) = Euler::makeRad(in0);
     euler(1) = Euler::makeRad(in1);
@@ -136,7 +166,7 @@ euler(3), axes(3) {
 }
 
 Euler::Euler(const Quantum<Vector<Double> > &in) :
-euler(3), axes(3) {
+data(get_arrays()), euler(*data.first), axes(*data.second) {
     Int i;
     Vector<Double> tmp = Euler::makeRad(in);
     Int j=tmp.size(); j=min(j,3);
@@ -150,7 +180,7 @@ euler(3), axes(3) {
 }
 
 Euler::Euler(const Quantum<Vector<Double> > &in, const Vector<uInt> &ax) :
-euler(3), axes(3) {
+data(get_arrays()), euler(*data.first), axes(*data.second) {
     Vector<Double> tmp = Euler::makeRad(in);
     Int j=tmp.size(); j=min(j,3); Int i=ax.size(); j=min(j,i);
     for (i=0; i<j; i++) {
@@ -165,7 +195,9 @@ euler(3), axes(3) {
 }
 
 //# Destructor
-Euler::~Euler() {}
+Euler::~Euler() {
+  return_arrays(data);
+}
 
 //# Operators
 Euler Euler::operator-() const {

casa/Quanta/Euler.h

@@ -33,6 +33,7 @@
 #include <casacore/casa/aips.h>
 #include <casacore/casa/Arrays/Vector.h>
 #include <casacore/casa/Quanta/Quantum.h>
+#include <utility>
 
 namespace casacore { //# NAMESPACE CASACORE - BEGIN
 
@@ -189,17 +190,26 @@ class Euler
 
 private:
 //# Data
-// Actual vector with 3 Euler angles
-    Vector<Double> euler;
-// Axes
-    Vector<Int> axes;
+    typedef std::pair<Vector<Double> *, Vector<Int> *> DataArrays;
+// data container
+    DataArrays data;
+// vector with 3 Euler angles (data.first)
+    Vector<Double> & euler;
+// Axes (data.second)
+    Vector<Int> & axes;
 
 //# Private Member Functions
 // The makeRad functions check and convert the input Quantities to radians
 // <group>
     static Double makeRad(const Quantity &in);
     static Vector<Double> makeRad(const Quantum<Vector<Double> > &in);
 // </group>
+    DataArrays get_arrays();
+    void return_arrays(DataArrays array);
+#if defined(AIPS_CXX11) && !defined(__APPLE__)
+    static thread_local DataArrays arrays[50];
+    static thread_local size_t available;
+#endif
 };
 
 

casa/Quanta/MVEpoch.cc

@@ -40,6 +40,7 @@ namespace casacore { //# NAMESPACE CASACORE - BEGIN
 
 //# Constants
 const Double MVEpoch::secInDay(3600*24);
+const Unit MVEpoch::unitDay("d");
 
 //# Constructors
 MVEpoch::MVEpoch() :
@@ -178,7 +179,7 @@ Double MVEpoch::get() const {
 }
 
 Quantity MVEpoch::getTime() const {
-  return (Quantity(get(), "d"));
+  return (Quantity(get(), unitDay));
 }
 
 Quantity MVEpoch::getTime(const Unit &unit) const {
@@ -249,7 +250,7 @@ Bool MVEpoch::putValue(const Vector<Quantum<Double> > &in) {
 
 Double MVEpoch::makeDay(const Quantity &in) const {
   in.assure(UnitVal::TIME);
-  return in.get("d").getValue();
+  return in.get(unitDay).getValue();
 }
 
 void MVEpoch::addTime(Double in) {

casa/Quanta/MVEpoch.h

@@ -141,6 +141,7 @@ class MVEpoch : public MeasValue {
   //# General Member Functions
   // Constants
   static const Double secInDay;
+  static const Unit unitDay;
 
   // Tell me your type
   // <group>

casa/Quanta/MVPosition.cc

@@ -46,15 +46,44 @@ namespace casacore { //# NAMESPACE CASACORE - BEGIN
 const Double MVPosition::loLimit = 743.568;
 const Double MVPosition::hiLimit = 743.569;
 
+// simplistic vector(3) cache to reduce allocation overhead for temporaries
+#if defined(AIPS_CXX11) && !defined(__APPLE__)
+static thread_local size_t available = 0;
+static thread_local Vector<Double> * arrays[50];
+#endif
+
+Vector<Double> * get_array()
+{
+#if defined(AIPS_CXX11) && !defined(__APPLE__)
+  if (available > 0) {
+    return arrays[--available];
+  }
+#endif
+  return new Vector<Double>(3);
+}
+
+void return_array(Vector<Double> * array)
+{
+#if defined(AIPS_CXX11) && !defined(__APPLE__)
+  if (available < sizeof(arrays) / sizeof(arrays[0]) &&
+		  array->size() == 3 &&
+		  array->nrefs() == 1) {
+    arrays[available++] = array;
+    return;
+  }
+#endif
+  delete array;
+}
+
 //# Constructors
 MVPosition::MVPosition() :
-  xyz(3) {
+  xyz(*get_array()) {
     xyz = Double(0.0);
 }
 
 MVPosition::MVPosition(const MVPosition &other) : 
   MeasValue(),
-  xyz(3)
+  xyz(*get_array())
 {
   xyz = other.xyz;
 }
@@ -67,27 +96,27 @@ MVPosition &MVPosition::operator=(const MVPosition &other) {
 }
 
 MVPosition::MVPosition(Double in) :
-  xyz(3) {
+	xyz(*get_array()) {
     xyz = Double(0.0);
     xyz(2) = in;
   }
 
 MVPosition::MVPosition(const Quantity &l) :
-  xyz(3) {
+  xyz(*get_array()) {
     xyz = Double(0.0);
     l.assure(UnitVal::LENGTH);
     xyz(2) = l.getBaseValue();
   }
 
 MVPosition::MVPosition(Double in0, Double in1, Double in2) : 
-  xyz(3) {
+	xyz(*get_array()) {
     xyz(0) = in0;
     xyz(1) = in1;
     xyz(2) = in2;
   }
 
 MVPosition::MVPosition(const Quantity &l, Double angle0, Double angle1) : 
-  xyz(3) {
+  xyz(*get_array()) {
   Double loc = std::cos(angle1);
   xyz(0) = std::cos(angle0)*loc;
   xyz(1) = std::sin(angle0)*loc;
@@ -100,7 +129,7 @@ MVPosition::MVPosition(const Quantity &l, Double angle0, Double angle1) :
 
 MVPosition::MVPosition(const Quantity &l, const Quantity &angle0, 
 		       const Quantity &angle1) : 
-  xyz(3) {
+  xyz(*get_array()) {
   Double loc = (cos(angle1)).getValue();
   xyz(0) = ((cos(angle0)).getValue()) * loc;
   xyz(1) = ((sin(angle0)).getValue()) * loc;
@@ -113,7 +142,7 @@ MVPosition::MVPosition(const Quantity &l, const Quantity &angle0,
 }
 
 MVPosition::MVPosition(const Quantum<Vector<Double> > &angle) :
-  xyz(3) {
+  xyz(*get_array()) {
   uInt i; i = angle.getValue().nelements();
   if (i > 3 ) {
     throw (AipsError("Illegeal vector length in MVPosition constructor"));
@@ -139,7 +168,7 @@ MVPosition::MVPosition(const Quantum<Vector<Double> > &angle) :
 
 MVPosition::MVPosition(const Quantity &l, 
 		       const Quantum<Vector<Double> > &angle) :
-  xyz(3) {
+	xyz(*get_array()) {
     uInt i; i = angle.getValue().nelements();
     if (i > 3 ) {
       throw (AipsError("Illegal vector length in MVPosition constructor"));
@@ -166,7 +195,7 @@ MVPosition::MVPosition(const Quantity &l,
   }
 
 MVPosition::MVPosition(const Vector<Double> &other) :
-  xyz(3) {
+	xyz(*get_array()) {
     uInt i; i = other.nelements();
     if (i > 3 ) {
       throw (AipsError("Illegal vector length in MVPosition constructor"));
@@ -190,14 +219,17 @@ MVPosition::MVPosition(const Vector<Double> &other) :
   }
 
 MVPosition::MVPosition(const Vector<Quantity> &other) :
-  xyz(3) {
+  xyz(*get_array()) {
   if (!putValue(other)) {
     throw (AipsError("Illegal quantum vector in MVPosition constructor"));
   }
 }
 
 //# Destructor
-MVPosition::~MVPosition() {}
+MVPosition::~MVPosition()
+{
+  return_array(&xyz);
+}
 
 //# Operators
 Bool MVPosition::
@@ -273,14 +305,12 @@ MVPosition MVPosition::operator-(const MVPosition &right) const{
 }
 
 MVPosition &MVPosition::operator*=(const RotMatrix &right) {
-  MVPosition result;
-  for (Int i=0; i<3; i++) {
-    result(i) = 0;
-    for (Int j=0; j<3; j++) {
-      result(i) += xyz(j) * right(j,i);
-    }
-  }
-  *this = result;
+  Double x = xyz(0);
+  Double y = xyz(1);
+  Double z = xyz(2);
+  xyz(0) = x * right(0, 0) + y * right(1, 0) + z * right(2, 0);
+  xyz(1) = x * right(0, 1) + y * right(1, 1) + z * right(2, 1);
+  xyz(2) = x * right(0, 2) + y * right(1, 2) + z * right(2, 2);
   return *this;
 }
 

casa/Quanta/MVPosition.h

@@ -155,7 +155,7 @@ class MVPosition : public MeasValue {
   MVPosition &operator=(const MVPosition &other);
 
   // Destructor
-  ~MVPosition();
+  virtual ~MVPosition();
 
   //# Operators
   // Multiplication defined as in-product
@@ -279,7 +279,7 @@ class MVPosition : public MeasValue {
   Double getLat(Double ln) const;
   //# Data
   // Position vector (in m)
-  Vector<Double> xyz;
+  Vector<Double> & xyz;
 };
 
 //# Global functions