Changed the implementation of transpose ( .t() ). Massive speed up

SpMat.h

@@ -142,7 +142,7 @@ public:
 
   const SpMat<T> TransMultSelf() const {return(TransMult(*this));}                 // Returns transpose(*this)*(*this)
 
-  const SpMat<T> t() const;                                                        // Returns transpose(*this). Avoid, if at all possible.
+  const SpMat<T> t() const;                                                        // Returns transpose(*this). 
   
   friend class Accumulator<T>; 
 
@@ -632,26 +632,27 @@ template<class T>
 const SpMat<T> SpMat<T>::t() const
 {
   SpMat<T>         t_mat(_n,_m);
-  Accumulator<T>   t_col(_n);
-  for (unsigned int new_col=0; new_col<_m; new_col++) {    // For all columns of transposed matrix
-    t_col.Reset();
-    for (unsigned int old_col=0; old_col<_n; old_col++) {  // Search old colums for row-index corresponding to new_col
-      int pos = 0;
-      if (found(_ri[old_col],new_col,pos)) {
-	t_col(old_col) = _val[old_col][pos];
-      }
-    }
-    t_mat._ri[new_col].resize(t_col.NO());  
-    t_mat._val[new_col].resize(t_col.NO());
-    std::vector<unsigned int>&   t_mat_ri = t_mat._ri[new_col];
-    std::vector<T>&              t_mat_val = t_mat._val[new_col];
-    for (unsigned int i=0; i<t_col.NO(); i++) {
-      t_mat_ri[i] = t_col.ri(i);
-      t_mat_val[i] = t_col.val(i);
+  // First make a list of number of elements in each row of *this (columns of t_mat)
+  std::vector<unsigned int> no_per_row(_m,static_cast<unsigned int>(0));
+  for (unsigned int i=0; i<_n; i++) {
+    for (unsigned int j=0; j<_ri[i].size(); j++) no_per_row[_ri[i][j]]++;
+  }
+  // A second pass to allocate storage
+  for (unsigned int i=0; i<_m; i++) {
+    t_mat._ri[i].resize(no_per_row[i]);
+    t_mat._val[i].resize(no_per_row[i]);
+  }
+  // A third pass to populate
+  no_per_row.assign(no_per_row.size(),static_cast<unsigned int>(0)); // Reset to use as counter
+  for (unsigned int i=0; i<_n; i++) {
+    for (unsigned int j=0; j<_ri[i].size(); j++) {
+      t_mat._ri[_ri[i][j]][no_per_row[_ri[i][j]]] = i;
+      t_mat._val[_ri[i][j]][no_per_row[_ri[i][j]]++] = _val[i][j];
     }
-    t_mat._nz += t_col.NO();
   }
-  return(t_mat);  
+  t_mat._nz = _nz;
+
+  return(t_mat);
 }
 
 /////////////////////////////////////////////////////////////////////