From 65aaf74a26f8d3d2b2b19ca4f13e8797dc897922 Mon Sep 17 00:00:00 2001
From: Christian Engwer <christi@dune-project.org>
Date: Fri, 14 Dec 2007 10:17:35 +0000
Subject: [PATCH] added test for bug #333

[[Imported from SVN: r841]]
---
 istl/overlappingschwarz.hh          | 18 ++++---
 istl/test/.gitignore                |  1 +
 istl/test/overlappingschwarztest.cc | 79 +++++++++++++++++++++--------
 3 files changed, 68 insertions(+), 30 deletions(-)

diff --git a/istl/overlappingschwarz.hh b/istl/overlappingschwarz.hh
index f35ff60b8..8fa771c5b 100644
--- a/istl/overlappingschwarz.hh
+++ b/istl/overlappingschwarz.hh
@@ -301,6 +301,8 @@ namespace Dune
     /** @brief The vector type containing the subdomain to row index mapping. */
     typedef std::vector<subdomain_type,allocator> subdomain_vector;
 
+    /** @brief The vector type containing the row index to subdomain mapping. */
+    typedef std::vector<SLList<size_type,TA>,TA> rowtodomain_vector;
 
     enum {
       //! \brief The category the precondtioner is part of.
@@ -322,7 +324,7 @@ namespace Dune
      * @param mat The matrix to precondition.
      * @param rowToDomain The mapping of the rows onto the domains.
      */
-    SeqOverlappingSchwarz(const matrix_type& mat, const std::vector<SLList<size_type,TA>,TA>& rowToDomain,
+    SeqOverlappingSchwarz(const matrix_type& mat, const rowtodomain_vector& rowToDomain,
                           field_type relaxationFactor=1);
 
     /*!
@@ -357,7 +359,7 @@ namespace Dune
 
     int maxlength;
 
-    void initialize(const std::vector<SLList<size_type,TA>,TA>& rowToDomain, const matrix_type& mat);
+    void initialize(const rowtodomain_vector& rowToDomain, const matrix_type& mat);
 
     template<typename T>
     struct Assigner
@@ -485,11 +487,11 @@ namespace Dune
   }
 
   template<class M, class X, class TM, class TA>
-  SeqOverlappingSchwarz<M,X,TM,TA>::SeqOverlappingSchwarz(const matrix_type& mat_, const std::vector<SLList<size_type,TA>,TA>& rowToDomain,
+  SeqOverlappingSchwarz<M,X,TM,TA>::SeqOverlappingSchwarz(const matrix_type& mat_, const rowtodomain_vector& rowToDomain,
                                                           field_type relaxationFactor)
     : mat(mat_), relax(relaxationFactor)
   {
-    typedef typename std::vector<SLList<size_type,TA>,TA>::const_iterator RowDomainIterator;
+    typedef typename rowtodomain_vector::const_iterator RowDomainIterator;
     typedef typename SLList<size_type,TA>::const_iterator DomainIterator;
 #ifdef DUNE_ISTL_WITH_CHECKING
     assert(rowToDomain.size()==mat.N());
@@ -555,7 +557,7 @@ namespace Dune
 #endif
 
     // Create a row to subdomain mapping
-    std::vector<SLList<size_type,TA>,TA> rowToDomain(mat.N());
+    rowtodomain_vector rowToDomain(mat.N());
 
     int domainId=0;
 
@@ -569,7 +571,7 @@ namespace Dune
   }
 
   template<class M, class X, class TM, class TA>
-  void SeqOverlappingSchwarz<M,X,TM,TA>::initialize(const std::vector<SLList<size_type,TA>,TA>& rowToDomain, const matrix_type& mat)
+  void SeqOverlappingSchwarz<M,X,TM,TA>::initialize(const rowtodomain_vector& rowToDomain, const matrix_type& mat)
   {
     typedef typename std::vector<SuperMatrixInitializer<matrix_type> >::iterator InitializerIterator;
     typedef typename subdomain_vector::const_iterator DomainIterator;
@@ -591,7 +593,7 @@ namespace Dune
 
     // Set up the supermatrices according to the subdomains
     typedef OverlappingSchwarzInitializer<std::vector<SuperMatrixInitializer<matrix_type> >,
-        std::vector<SLList<size_type,TA>,TA> > Initializer;
+        rowtodomain_vector > Initializer;
 
     Initializer initializer(initializers, rowToDomain);
     copyToSuperMatrix(initializer, mat);
@@ -675,7 +677,7 @@ namespace Dune
 
     // loop over all Matrix row entries and calculate defect.
     typedef typename M::ConstColIterator col_iterator;
-    typedef typename  subdomain_type::const_iterator domain_iterator;
+    typedef typename subdomain_type::const_iterator domain_iterator;
 
     for(col_iterator col=(*mat)[domainIndex].begin(); col!=(*mat)[domainIndex].end(); ++col) {
       typename X::block_type tmp;
diff --git a/istl/test/.gitignore b/istl/test/.gitignore
index 5574e7820..5370ee5f3 100644
--- a/istl/test/.gitignore
+++ b/istl/test/.gitignore
@@ -13,6 +13,7 @@ gmon.out
 vectorcommtest
 matrixtest
 matrixiteratortest
+overlappingschwarztest
 bcrsbuildtest
 superlutest
 mv
diff --git a/istl/test/overlappingschwarztest.cc b/istl/test/overlappingschwarztest.cc
index 107079aa9..02b9bf445 100644
--- a/istl/test/overlappingschwarztest.cc
+++ b/istl/test/overlappingschwarztest.cc
@@ -10,6 +10,7 @@
 #include <dune/common/timer.hh>
 #include <dune/common/sllist.hh>
 #include <dune/istl/overlappingschwarz.hh>
+
 int main(int argc, char** argv)
 {
 
@@ -24,12 +25,12 @@ int main(int argc, char** argv)
   typedef Dune::FieldMatrix<double,BS,BS> MatrixBlock;
   typedef Dune::BCRSMatrix<MatrixBlock> BCRSMat;
   typedef Dune::FieldVector<double,BS> VectorBlock;
-  typedef Dune::BlockVector<VectorBlock> Vector;
-  typedef Dune::MatrixAdapter<BCRSMat,Vector,Vector> Operator;
+  typedef Dune::BlockVector<VectorBlock> BVector;
+  typedef Dune::MatrixAdapter<BCRSMat,BVector,BVector> Operator;
 
   BCRSMat mat;
   Operator fop(mat);
-  Vector b(N*N), x(N*N);
+  BVector b(N*N), x(N*N);
 
   setupLaplacian(mat,N);
   b=0;
@@ -48,11 +49,14 @@ int main(int argc, char** argv)
   int domainsPerDim=(N+domainSize-1)/domainSize;
 
   // set up the overlapping domains
-  typedef Dune::SeqOverlappingSchwarz<BCRSMat,Vector> Schwarz;
+  typedef Dune::SeqOverlappingSchwarz<BCRSMat,BVector> Schwarz;
   typedef Schwarz::subdomain_vector subdomain_vector;
 
   subdomain_vector domains(domainsPerDim*domainsPerDim);
 
+  // set up the rowToDomain vector
+  typedef Schwarz::rowtodomain_vector rowtodomain_vector;
+  rowtodomain_vector rowToDomain(N*N);
 
   for(int j=0; j < N; ++j)
     for(int i=0; i < N; ++i)
@@ -60,27 +64,41 @@ int main(int argc, char** argv)
       int xdomain = i/domainSize;
       int ydomain = j/domainSize;
       int mainDomain=ydomain*domainsPerDim+xdomain;
-      domains[mainDomain].insert(j*N+i);
+      int id=j*N+i;
+      domains[mainDomain].insert(id);
+      rowToDomain[id].push_back(mainDomain);
 
       // check left domain
       int domain = (i-overlap)/domainSize;
       if(domain>=0 && domain<domainsPerDim)
-        domains[ydomain*domainsPerDim+domain].insert(j*N+i);
+      {
+        domains[ydomain*domainsPerDim+domain].insert(id);
+        rowToDomain[id].push_back(ydomain*domainsPerDim+domain);
+      }
 
       //check right domain
       domain = (i+overlap)/domainSize;
       if(domain>=0 && domain<domainsPerDim)
-        domains[ydomain*domainsPerDim+domain].insert(j*N+i);
+      {
+        domains[ydomain*domainsPerDim+domain].insert(id);
+        rowToDomain[id].push_back(ydomain*domainsPerDim+domain);
+      }
 
       // check lower domain
       domain = (j-overlap)/domainSize;
       if(domain>=0 && domain<domainsPerDim)
-        domains[domain*domainsPerDim+xdomain].insert(j*N+i);
+      {
+        domains[domain*domainsPerDim+xdomain].insert(id);
+        rowToDomain[id].push_back(domain*domainsPerDim+xdomain);
+      }
 
       //check right domain
       domain = (j+overlap)/domainSize;
       if(domain>=0 && domain<domainsPerDim)
-        domains[domain*domainsPerDim+xdomain].insert(j*N+i);
+      {
+        domains[domain*domainsPerDim+xdomain].insert(id);
+        rowToDomain[id].push_back(domain*domainsPerDim+xdomain);
+      }
     }
 
   typedef subdomain_vector::const_iterator iterator;
@@ -100,38 +118,55 @@ int main(int argc, char** argv)
     Dune::printvector(std::cout, x, std::string("X"), std::string("X"));
   }
 
-  Dune::SeqOverlappingSchwarz<BCRSMat,Vector> prec(mat, domains, 1);
-
-
   Dune::Timer watch;
 
   watch.reset();
 
-  Dune::LoopSolver<Vector> solver(fop, prec, 1e-2,100,2);
   Dune::InverseOperatorResult res;
 
+  std::cout<<"Additive Schwarz (domains vector)"<<std::endl;
 
   //  b=0;
   //  x=100;
   //  setBoundary(x,b,N);
-  std::cout<<"Additive Schwarz"<<std::endl;
-  solver.apply(x,b, res);
+  Dune::SeqOverlappingSchwarz<BCRSMat,BVector> prec0(mat, domains, 1);
+  Dune::LoopSolver<BVector> solver0(fop, prec0, 1e-2,100,2);
+  solver0.apply(x,b, res);
+
 
-  Dune::SeqOverlappingSchwarz<BCRSMat,Vector,Dune::MultiplicativeSchwarzMode> prec1(mat,domains, 1);
-  Dune::LoopSolver<Vector> solver1(fop, prec1, 1e-2,100,2);
+  std::cout << "Multiplicative Schwarz (domains vector)"<<std::endl;
 
   b=0;
   x=100;
   //setBoundary(x,b,N);
-  std::cout << "Multiplicative Schwarz"<<std::endl;
-
+  Dune::SeqOverlappingSchwarz<BCRSMat,BVector,Dune::MultiplicativeSchwarzMode> prec1(mat, domains, 1);
+  Dune::LoopSolver<BVector> solver1(fop, prec1, 1e-2,100,2);
   solver1.apply(x,b, res);
 
-  Dune::SeqSOR<BCRSMat,Vector,Vector> sor(mat, 1,1);
-  Dune::LoopSolver<Vector> solver2(fop, sor, 1e-2,100,2);
+  std::cout<<"Additive Schwarz (rowToDomain vector)"<<std::endl;
+
+  b=0;
+  x=100;
+  //  setBoundary(x,b,N);
+  Dune::SeqOverlappingSchwarz<BCRSMat,BVector> prec2(mat, rowToDomain, 1);
+  Dune::LoopSolver<BVector> solver2(fop, prec2, 1e-2,100,2);
+  solver2.apply(x,b, res);
+
+  std::cout << "Multiplicative Schwarz (rowToDomain vector)"<<std::endl;
+
   b=0;
   x=100;
   //setBoundary(x,b,N);
+  Dune::SeqOverlappingSchwarz<BCRSMat,BVector,Dune::MultiplicativeSchwarzMode> prec3(mat, rowToDomain, 1);
+  Dune::LoopSolver<BVector> solver3(fop, prec3, 1e-2,100,2);
+  solver3.apply(x,b, res);
+
   std::cout << "SOR"<<std::endl;
-  solver2.apply(x,b, res);
+
+  b=0;
+  x=100;
+  //setBoundary(x,b,N);
+  Dune::SeqSOR<BCRSMat,BVector,BVector> sor(mat, 1,1);
+  Dune::LoopSolver<BVector> solver4(fop, sor, 1e-2,100,2);
+  solver4.apply(x,b, res);
 }
-- 
GitLab