Erweiterungen und Aufräumarbeiten am LowRankOperator

[[Imported from SVN: r5253]]

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
         {

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)