diff --git a/dune/solvers/operators/lowrankoperator.hh b/dune/solvers/operators/lowrankoperator.hh index b47da905b7142d6df5dd3bbfb6e30347e4ab08e3..5f70e2fa0e3191710f021f54ea534b9d5bfcc260 100644 --- a/dune/solvers/operators/lowrankoperator.hh +++ b/dune/solvers/operators/lowrankoperator.hh @@ -21,26 +21,40 @@ class LowRankOperator //! export size type typedef typename LRFactorType::size_type size_type; + //! default constructor + LowRankOperator(): + factor_allocated_internally_(true) + { + lowRankFactor_ = new LowRankFactorType; + } + //! constructor taking the defining factor as argument LowRankOperator(LowRankFactorType& lrFactor): - lowRankFactor_(lrFactor) + lowRankFactor_(&lrFactor), + factor_allocated_internally_(false) {} + ~LowRankOperator() + { + if (factor_allocated_internally_) + delete lowRankFactor_; + } + //! b += Ax template <class LVectorType, class RVectorType> void umv(const LVectorType& x, RVectorType& b) const { - LVectorType temp(lowRankFactor_.N()); - lowRankFactor_.mv(x,temp); + LVectorType temp(lowRankFactor_->N()); + lowRankFactor_->mv(x,temp); - typename LowRankFactorType::RowIterator row_end = lowRankFactor_.end(); - typename LowRankFactorType::RowIterator row_it = lowRankFactor_.begin(); + typename LowRankFactorType::ConstRowIterator row_end = lowRankFactor_->end(); + typename LowRankFactorType::ConstRowIterator row_it = lowRankFactor_->begin(); for (; row_it!=row_end; ++row_it) { - typename LowRankFactorType::ColIterator col_end = row_it->end(); - typename LowRankFactorType::ColIterator col_it = row_it->begin(); + typename LowRankFactorType::ConstColIterator col_end = row_it->end(); + typename LowRankFactorType::ConstColIterator col_it = row_it->begin(); for (; col_it!=col_end; ++col_it) - lowRankFactor_[row_it.index()][col_it.index()].umv(temp[row_it.index()],b[col_it.index()]); + (*lowRankFactor_)[row_it.index()][col_it.index()].umv(temp[row_it.index()],b[col_it.index()]); } } @@ -48,17 +62,35 @@ class LowRankOperator template <class LVectorType, class RVectorType> void usmv(const field_type& alpha, const LVectorType& x, RVectorType& b) const { - LVectorType temp(lowRankFactor_.N()); - lowRankFactor_.mv(x,temp); + LVectorType temp(lowRankFactor_->N()); + lowRankFactor_->mv(x,temp); - typename LowRankFactorType::RowIterator row_end = lowRankFactor_.end(); - typename LowRankFactorType::RowIterator row_it = lowRankFactor_.begin(); + typename LowRankFactorType::ConstRowIterator row_end = lowRankFactor_->end(); + typename LowRankFactorType::ConstRowIterator row_it = lowRankFactor_->begin(); for (; row_it!=row_end; ++row_it) { - typename LowRankFactorType::ColIterator col_end = row_it->end(); - typename LowRankFactorType::ColIterator col_it = row_it->begin(); + typename LowRankFactorType::ConstColIterator col_end = row_it->end(); + typename LowRankFactorType::ConstColIterator col_it = row_it->begin(); for (; col_it!=col_end; ++col_it) - lowRankFactor_[row_it.index()][col_it.index()].usmv(alpha,temp[row_it.index()],b[col_it.index()]); + (*lowRankFactor_)[row_it.index()][col_it.index()].usmv(alpha,temp[row_it.index()],b[col_it.index()]); + } + } + + //! b -= Ax + template <class LVectorType, class RVectorType> + void mmv(const LVectorType& x, RVectorType& b) const + { + LVectorType temp(lowRankFactor_->N()); + lowRankFactor_->mv(x,temp); + + typename LowRankFactorType::ConstRowIterator row_end = lowRankFactor_->end(); + typename LowRankFactorType::ConstRowIterator row_it = lowRankFactor_->begin(); + for (; row_it!=row_end; ++row_it) + { + typename LowRankFactorType::ConstColIterator col_end = row_it->end(); + typename LowRankFactorType::ConstColIterator col_it = row_it->begin(); + for (; col_it!=col_end; ++col_it) + (*lowRankFactor_)[row_it.index()][col_it.index()].mmv(temp[row_it.index()],b[col_it.index()]); } } @@ -71,26 +103,33 @@ class LowRankOperator } //! returns number of rows - size_type N() + const size_type N() const { - return lowRankFactor_.M(); + return lowRankFactor_->M(); } //! returns number of columns - size_type M() + const size_type M() const { - return lowRankFactor_.M(); + return lowRankFactor_->M(); } //! returns the defining factor - const LowRankFactorType& getLowRankFactor() const + const LowRankFactorType& lowRankFactor() const { - return lowRankFactor_; + return *lowRankFactor_; } + //! returns the defining factor + LowRankFactorType& lowRankFactor() + { + return *lowRankFactor_; + } private: //! the defining factor - LowRankFactorType& lowRankFactor_; + LowRankFactorType* lowRankFactor_; + + const bool factor_allocated_internally_; };