Changed, and massively speeded up, implementation of TransMult and TransMultSelf methods

SpMat.h

@@ -135,13 +135,13 @@ public:
   const NEWMAT::ReturnMatrix TransMult(const NEWMAT::ColumnVector& x) const {
     return(trans_mult(x));                                                         // Duplication for compatibility with IML++
   }
-  const SpMat<T> TransMult(const SpMat<T>& B) const;                               // Multiplication of transpose(*this) with sparse matrix B
+  const SpMat<T> TransMult(const SpMat<T>& B) const { return(this->t()*B); }       // Multiplication of transpose(*this) with sparse matrix B
   SpMat<T>& operator*=(double s);                                                  // Multiplication of self with scalar
   SpMat<T> operator-(const SpMat<T>& M) const {return(SpMat<T>(M) *= -1.0);}       // Unary minus
   SpMat<T>& operator|=(const SpMat<T>& rh);                                        // Hor concat to right
   SpMat<T>& operator&=(const SpMat<T>& bh);                                        // Vert concat below
 
-  const SpMat<T> TransMultSelf() const {return(TransMult(*this));}                 // Returns transpose(*this)*(*this)
+  const SpMat<T> TransMultSelf() const {return(this->t() * *this);}                 // Returns transpose(*this)*(*this)
 
   const SpMat<T> t() const;                                                        // Returns transpose(*this). 
   
@@ -745,9 +745,11 @@ const NEWMAT::ReturnMatrix SpMat<T>::operator*(const NEWMAT::ColumnVector& x) co
 /////////////////////////////////////////////////////////////////////
 //
 // Multiply transpose with sparse matrix B returning matrix C (C = A'*B)
+// This implementation was found to be inferior to C = A.t()*B, so
+// it has been replaced by that for now and the code below commented out.
 //
 /////////////////////////////////////////////////////////////////////
-
+/*
 template<class T>
 const SpMat<T> SpMat<T>::TransMult(const SpMat<T>& B) const
 {
@@ -785,7 +787,7 @@ const SpMat<T> SpMat<T>::TransMult(const SpMat<T>& B) const
  
   return(C);
 }
-
+*/
 /////////////////////////////////////////////////////////////////////
 //
 // Multiply transpose with vector x returning vector b (b = A'*x)