diff --git a/dune/solvers/operators/lowrankoperator.hh b/dune/solvers/operators/lowrankoperator.hh index 916112d3df30821ca9b0366cd38fc2fc40caa85a..b47da905b7142d6df5dd3bbfb6e30347e4ab08e3 100644 --- a/dune/solvers/operators/lowrankoperator.hh +++ b/dune/solvers/operators/lowrankoperator.hh @@ -2,7 +2,7 @@ #define LOWRANKOPERATOR_HH -/** \brief Class that behaves like a filled in low-rank matrix defined via a factor \f$ M\in\mathbb R^{k\times n}, k<<n \f$ +/** \brief Class that behaves like a filled in square low-rank matrix defined via a factor \f$ M\in\mathbb R^{k\times n}, k<<n \f$ * * The full matrix is given by \f$ A=M^TM\f$ * @@ -14,6 +14,12 @@ class LowRankOperator public: //! export the type of the defining factor typedef LRFactorType LowRankFactorType; + //! export block type + typedef typename LRFactorType::block_type block_type; + //! export field type + typedef typename LRFactorType::field_type field_type; + //! export size type + typedef typename LRFactorType::size_type size_type; //! constructor taking the defining factor as argument LowRankOperator(LowRankFactorType& lrFactor): @@ -27,10 +33,33 @@ class LowRankOperator LVectorType temp(lowRankFactor_.N()); lowRankFactor_.mv(x,temp); - for (size_t i=0; i<b.size(); ++i) - for (size_t j=0; j<lowRankFactor_.N(); ++j) - if (lowRankFactor_.exists(j,i)) - lowRankFactor_[j][i].umv(temp[j],b[i]); + typename LowRankFactorType::RowIterator row_end = lowRankFactor_.end(); + typename LowRankFactorType::RowIterator 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(); + for (; col_it!=col_end; ++col_it) + lowRankFactor_[row_it.index()][col_it.index()].umv(temp[row_it.index()],b[col_it.index()]); + } + } + + //! b += alpha*Ax + 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); + + typename LowRankFactorType::RowIterator row_end = lowRankFactor_.end(); + typename LowRankFactorType::RowIterator 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(); + for (; col_it!=col_end; ++col_it) + lowRankFactor_[row_it.index()][col_it.index()].usmv(alpha,temp[row_it.index()],b[col_it.index()]); + } } //! b = Ax @@ -41,6 +70,18 @@ class LowRankOperator umv(x,b); } + //! returns number of rows + size_type N() + { + return lowRankFactor_.M(); + } + + //! returns number of columns + size_type M() + { + return lowRankFactor_.M(); + } + //! returns the defining factor const LowRankFactorType& getLowRankFactor() const { diff --git a/dune/solvers/test/lowrankoperatortest.cc b/dune/solvers/test/lowrankoperatortest.cc index 8a456e293617020edd781391ea39de788354b27b..0b3f3c4de8fbe64cb120dcb3a9c0270945498621 100644 --- a/dune/solvers/test/lowrankoperatortest.cc +++ b/dune/solvers/test/lowrankoperatortest.cc @@ -80,6 +80,8 @@ bool check() x[i] = (1.0*rand()) / RAND_MAX; b = ref_vec = 0.0; + double alpha = (1.0*rand())/RAND_MAX; + LowRankFactorType lr_factor(m,n); std::cout << "Checking LowRankOperator<" << Dune::className(lr_factor) << " > ... "; @@ -142,6 +144,20 @@ bool check() break; } + full_op.usmv(alpha,x,ref_vec); + + lr_op.usmv(alpha,x,b); + for (size_t i=0; i<b.size(); ++i) + if ((b[i] - ref_vec[i]).two_norm()/b[i].two_norm()>1e-12) + { + std::cout << "LowRankOperator::usmv does not yield correct result. Difference = " << (b[i] - ref_vec[i]).two_norm() << std::endl; + passed = false; + for (size_t j=0; j<b.size(); ++j) + std::cout << b[j] << " "<<ref_vec[j] << std::endl; + std::cout << std::endl; + break; + } + LowRankFactorType lr_factor_reborn = lr_op.getLowRankFactor(); if (passed)