diff --git a/sparsesmat.h b/sparsesmat.h index 1cd19e5..25ca17d 100644 --- a/sparsesmat.h +++ b/sparsesmat.h @@ -87,7 +87,8 @@ public: void gemm(const T beta, const SparseSMat &a, const char transa, const SparseSMat &b, const char transb, const T alpha); //this := alpha*op( A )*op( B ) + beta*this !!!NOT A GENERAL ROUTINE, JUST FOR THE CASES WHEN THE RESULT STAYS SYMMETRIC inline void add(const SPMatindex n, const SPMatindex m, const T elem, const bool both=true); inline unsigned long long length() {return simplify();}; - void transposeme() const {}; + void transposeme() const {laerror("in-place transposition not necessary/implemented for SparseSMat");}; + SparseSMat transpose(bool conj=false) const; //if we store a non-symmetric matrix there void get(int fd, bool dimen, bool transp); void put(int fd, bool dimen, bool transp) const; int nrows() const {return nn;} @@ -404,11 +405,20 @@ std::istream& operator>>(std::istream &s, SparseSMat &x) return s; } +template +SparseSMat SparseSMat::transpose(bool conj) const +{ +SparseSMat r(nn); +typename SparseSMat::iterator p(*this); +for(; p.notend(); ++p) r.add(p->col, p->row, (conj?LA_traits::conjugate(p->elem):p->elem), false); +return r; +} + //Cholesky decomposition, pivoted, positive semidefinite, not in place //it is NOT checked that the input matrix is symmetric/hermitean //result.transpose(true)*result reproduces the original matrix -//Due to pivoting the result is upper triangular only before permutation +//Due to pivoting the result is upper triangular only before applying final permutation // template SparseSMat SparseSMat::cholesky(void) const @@ -427,6 +437,8 @@ for(int i=0; i