Added transpose ( .t() ) to SpMat class

SpMat.h

@@ -11,14 +11,6 @@
 //  for example inserting elements in a random order may be
 //  a bit slow.  
 //
-//  You will notice that some "obvious" functionality like e.g.
-//  transpose is missing. Instead I have supplied functionality
-//  like A.trans_mult(x) which is equivalent to A.t()*x in NEWMAT
-//  lingo. The reason for this is partly to supply an api that is
-//  consistent with that expected by IML++. But more so because I
-//  believe that when a matrix is unwieldily enough to warrant
-//  the use of a sparse matrix class, then one should not attempt
-//  to represent both the matrix and its transpose.
 //
 // Jesper Andersson, FMRIB Image Analysis Group
 //
@@ -134,6 +126,8 @@ public:
   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> t() const;                                                        // Returns transpose(*this). Avoid, if at all possible.
   
   friend class Accumulator<T>; 
 
@@ -146,18 +140,18 @@ public:
                                  unsigned int                          miter = 200,
 				 boost::shared_ptr<Preconditioner<T> > C = boost::shared_ptr<Preconditioner<T> >()) const;
 
-NEWMAT::ReturnMatrix SolveForx(const NEWMAT::ColumnVector&            b,
-			                 MatrixType                             type,
-			                 double                                 tol,
-			                 unsigned int                           miter,
-					 const NEWMAT::ColumnVector&            x_init) const;
+  NEWMAT::ReturnMatrix SolveForx(const NEWMAT::ColumnVector&            b,
+			         MatrixType                             type,
+			         double                                 tol,
+			         unsigned int                           miter,
+				 const NEWMAT::ColumnVector&            x_init) const;
 
-NEWMAT::ReturnMatrix SolveForx(const NEWMAT::ColumnVector&            b,
-			                 MatrixType                             type,
-			                 double                                 tol,
-			                 unsigned int                           miter,
-					 boost::shared_ptr<Preconditioner<T> >  C,
-					 const NEWMAT::ColumnVector&           x_init) const;
+  NEWMAT::ReturnMatrix SolveForx(const NEWMAT::ColumnVector&            b,
+			         MatrixType                             type,
+			         double                                 tol,
+			         unsigned int                           miter,
+				 boost::shared_ptr<Preconditioner<T> >  C,
+				 const NEWMAT::ColumnVector&            x_init) const;
 
 private:
   unsigned int                              _m;
@@ -461,6 +455,39 @@ NEWMAT::ReturnMatrix SpMat<T>::SolveForx(const NEWMAT::ColumnVector&
   x.Release();
   return(x);
 }
+
+/////////////////////////////////////////////////////////////////////
+//
+// Returns a sparse matrix that is the transpose of *this
+//
+/////////////////////////////////////////////////////////////////////
+
+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);
+    }
+    t_mat._nz += t_col.NO();
+  }
+  return(t_mat);  
+}
+
 /////////////////////////////////////////////////////////////////////
 //
 // Returns value at position i,j (one offset)