From 53ec4c3e9088d3acaa6a0358c792358796ecde5b Mon Sep 17 00:00:00 2001 From: Jesper Andersson <jesper@fmrib.ox.ac.uk> Date: Mon, 10 Sep 2007 14:54:53 +0000 Subject: [PATCH] Added transpose ( .t() ) to SpMat class --- SpMat.h | 65 ++++++++++++++++++++++++++++++++++++++++----------------- 1 file changed, 46 insertions(+), 19 deletions(-) diff --git a/SpMat.h b/SpMat.h index 57823ea..32a37a7 100644 --- a/SpMat.h +++ b/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) -- GitLab