[Rivet-svn] r1990 - trunk/include/Rivet/Math

[blackhole at projects.hepforge.org]

Author: buckley
Date: Wed Nov  4 00:11:28 2009
New Revision: 1990

Log:
More matrix/vector improvements

Modified:
   trunk/include/Rivet/Math/MatrixN.hh
   trunk/include/Rivet/Math/VectorN.hh

Modified: trunk/include/Rivet/Math/MatrixN.hh
==============================================================================
--- trunk/include/Rivet/Math/MatrixN.hh	Tue Nov  3 17:25:47 2009	(r1989)
+++ trunk/include/Rivet/Math/MatrixN.hh	Wed Nov  4 00:11:28 2009	(r1990)
@@ -131,57 +131,61 @@
   //   return *this;
   // }
 
+  /// Calculate inverse
   Matrix<N> inverse() const {
     Matrix<N> tmp;
     tmp._matrix = _matrix.inverse();
     return tmp;
   }
 
+  /// Calculate determinant
   double det() const  {
     return _matrix.determinant();
   }
 
+  /// Calculate trace
   double trace() const {
-    double tr = 0.0;
-    for (size_t i = 0; i < N; ++i) {
-      tr += _matrix(i,i);
-    }
-    return tr;
-    // return _matrix.trace();
+    return _matrix.trace();
   }
 
+  /// Negate
   Matrix<N> operator-() const {
     Matrix<N> rtn;
     rtn._matrix = -_matrix;
     return rtn;
   }
 
+  /// Get dimensionality
   const size_t size() const {
     return N;
   }
 
-  const bool isZero() const {
+  /// Index-wise check for nullness, allowing for numerical precision
+  const bool isZero(double tolerance=1E-5) const {
     for (size_t i=0; i < N; ++i) {
       for (size_t j=0; j < N; ++j) {
-        if (! Rivet::isZero(_matrix(i,j)) ) return false;
+        if (! Rivet::isZero(_matrix(i,j), tolerance) ) return false;
       }
     }
     return true;
   }
 
+  /// Check for index-wise equality, allowing for numerical precision
   const bool isEqual(Matrix<N> other) const {
     for (size_t i=0; i < N; ++i) {
       for (size_t j=i; j < N; ++j) {
-        if (! Rivet::isZero(_matrix(i,j) - other._matrix(i,j)) ) return false;
+        if (! Rivet::fuzzyEquals(_matrix(i,j), other._matrix(i,j)) ) return false;
       }
     }
     return true;
   }
 
+  /// Check for symmetry under transposition
   const bool isSymm() const {
     return isEqual(this->transpose());
   }
 
+  /// Check that all off-diagonal elements are zero, allowing for numerical precision
   const bool isDiag() const {
     for (size_t i=0; i < N; ++i) {
       for (size_t j=0; j < N; ++j) {
@@ -356,6 +360,7 @@
 /////////////////////////////////
 
 
+/// Make string representation
 template <size_t N>
 inline string toString(const Matrix<N>& m) {
   ostringstream ss;
@@ -373,6 +378,7 @@
 }
 
 
+/// Stream out string representation
 template <size_t N>
 inline ostream& operator<<(std::ostream& out, const Matrix<N>& m) {
   out << toString(m);
@@ -380,6 +386,10 @@
 }
 
 
+/////////////////////////////////////////////////
+
+
+/// Compare two matrices by index, allowing for numerical precision
 template <size_t N>
 inline bool fuzzyEquals(const Matrix<N>& ma, const Matrix<N>& mb, double tolerance=1E-5) {
   for (size_t i = 0; i < N; ++i) {
@@ -393,6 +403,13 @@
 }
 
 
+/// External form of numerically safe nullness check
+template <size_t N>
+inline bool isZero(const Matrix<N>& m, double tolerance=1E-5) {
+  return m.isZero(tolerance);
+}
+
+
 }
 
 

Modified: trunk/include/Rivet/Math/VectorN.hh
==============================================================================
--- trunk/include/Rivet/Math/VectorN.hh	Tue Nov  3 17:25:47 2009	(r1989)
+++ trunk/include/Rivet/Math/VectorN.hh	Wed Nov  4 00:11:28 2009	(r1990)
@@ -49,6 +49,7 @@
       return get(index);
     }
 
+    /// Set indexed value
     Vector<N>& set(const size_t index, const double value) {
       if (index >= N) {
         throw std::runtime_error("Tried to access an invalid vector index.");
@@ -58,13 +59,15 @@
       return *this;
     }
 
+    /// Vector dimensionality
     const size_t size() const {
       return N;
     }
 
-    const bool isZero() const {
+    /// Check for nullness, allowing for numerical precision
+    const bool isZero(double tolerance=1E-5) const {
       for (size_t i=0; i < N; ++i) {
-        if (! Rivet::isZero(_vec[i]) ) return false;
+        if (! Rivet::isZero(_vec[i], tolerance) ) return false;
       }
       return true;
     }
@@ -154,6 +157,7 @@
   /////////////////////////////////////////////////
 
 
+  /// Make string representation
   template <size_t N>
   inline const string toString(const Vector<N>& v) {
     ostringstream out;
@@ -166,6 +170,7 @@
     return out.str();
   }
 
+  /// Stream out string representation
   template <size_t N>
   inline std::ostream& operator<<(std::ostream& out, const Vector<N>& v) {
     out << toString(v);
@@ -176,6 +181,7 @@
   /////////////////////////////////////////////////
 
 
+  /// Compare two vectors by index, allowing for numerical precision
   template <size_t N>
   inline bool fuzzyEquals(const Vector<N>& va, const Vector<N>& vb, double tolerance=1E-5) {
     for (size_t i = 0; i < N; ++i) {
@@ -187,6 +193,13 @@
   }
 
 
+  /// External form of numerically safe nullness check
+  template <size_t N>
+  inline bool isZero(const Vector<N>& v, double tolerance=1E-5) {
+    return v.isZero(tolerance);
+  }
+
+
 }
 
 #endif