From 67ecd248591e16b16c5ddd31ceda2a3cdd9dffce Mon Sep 17 00:00:00 2001
From: Markus Blatt <mblatt@dune-project.org>
Date: Thu, 27 Sep 2007 06:23:49 +0000
Subject: [PATCH] Added and improved multiplicative Schwarz methods.
[[Imported from SVN: r820]]
---
istl/overlappingschwarz.hh | 51 +++++++++++++++++---
istl/test/laplacian.hh | 53 +++++++++++++++++----
istl/test/overlappingschwarztest.cc | 73 ++++++++++++++++++++++-------
3 files changed, 144 insertions(+), 33 deletions(-)
diff --git a/istl/overlappingschwarz.hh b/istl/overlappingschwarz.hh
index e1233ae9e..f35ff60b8 100644
--- a/istl/overlappingschwarz.hh
+++ b/istl/overlappingschwarz.hh
@@ -24,7 +24,7 @@ namespace Dune
/**
* @file
* @author Markus Blatt
- * @brief Contains an additive overlapping Schwarz preconditioner
+ * @brief Contains one level overlapping Schwarz preconditioners
*/
#ifdef HAVE_SUPERLU
@@ -104,6 +104,9 @@ namespace Dune
struct MultiplicativeSchwarzMode
{};
+ struct SymmetricMultiplicativeSchwarzMode
+ {};
+
namespace
{
template<typename T>
@@ -161,6 +164,12 @@ namespace Dune
typedef MultiplicativeAdder<X> Adder;
};
+ template<class X>
+ struct AdderSelector<SymmetricMultiplicativeSchwarzMode,X>
+ {
+ typedef MultiplicativeAdder<X> Adder;
+ };
+
template<typename T1, typename T2, bool forward>
struct ApplyHelper
{
@@ -217,10 +226,34 @@ namespace Dune
}
};
+ template<class T>
+ struct Applier
+ {
+ typedef T smoother;
+ typedef typename smoother::range_type range_type;
+
+ static void apply(smoother& sm, range_type& v, const range_type& b)
+ {
+ sm.template apply<true>(v, b);
+ }
+ };
+
+ template<class M, class X, class TA>
+ struct Applier<SeqOverlappingSchwarz<M,X,SymmetricMultiplicativeSchwarzMode, TA> >
+ {
+ typedef SeqOverlappingSchwarz<M,X,SymmetricMultiplicativeSchwarzMode, TA> smoother;
+ typedef typename smoother::range_type range_type;
+
+ static void apply(smoother& sm, range_type& v, const range_type& b)
+ {
+ sm.template apply<true>(v, b);
+ sm.template apply<false>(v, b);
+ }
+ };
}
/**
- * @brief Sequential additive overlapping Schwarz preconditioner
+ * @brief Sequential overlapping Schwarz preconditioner
*/
template<class M, class X, class TM=AdditiveSchwarzMode, class TA=std::allocator<X> >
class SeqOverlappingSchwarz
@@ -246,7 +279,7 @@ namespace Dune
* @brief The mode (additive or multiplicative) of the Schwarz
* method.
*
- * Either Schwarz::Additive or Schwarz::Multiplicative
+ * Either AdditiveSchwarzMode or MultiplicativeSchwarzMode
*/
typedef TM Mode;
@@ -342,6 +375,7 @@ namespace Dune
const BlockVector<FieldVector<T,n>,A>* x;
int i;
};
+
};
@@ -481,6 +515,7 @@ namespace Dune
for(DomainIterator d=iter->begin(); d != iter->end(); ++d)
subDomains[*d].insert(row);
+#ifdef DUNE_ISTL_WITH_CHECKING
int i=0;
typedef typename std::vector<std::set<size_type> >::const_iterator iterator;
for(iterator iter=subDomains.begin(); iter != subDomains.end(); ++iter) {
@@ -491,6 +526,7 @@ namespace Dune
}
std::cout<<std::endl;
}
+#endif
initialize(rowToDomain, mat);
}
@@ -502,14 +538,14 @@ namespace Dune
{
typedef typename subdomain_vector::const_iterator DomainIterator;
-#ifndef NDEBUG
+#ifdef DUNE_ISTL_WITH_CHECKING
int i=0;
for(DomainIterator d=sd.begin(); d != sd.end(); ++d,++i) {
//std::cout<<i<<": "<<d->size()<<std::endl;
assert(d->size()>0);
typedef typename DomainIterator::value_type::const_iterator entry_iterator;
- std::cout<<"domain "<<i++<<":";
+ std::cout<<"domain "<<i<<":";
for(entry_iterator entry = d->begin(); entry != d->end(); ++entry) {
std::cout<<" "<<*entry;
}
@@ -549,6 +585,7 @@ namespace Dune
++initializer, ++solver, ++domain) {
solver->mat.N_=domain->size();
solver->mat.M_=domain->size();
+ //solver->setVerbosity(true);
*initializer=SuperMatrixInitializer<matrix_type>(solver->mat);
}
@@ -578,8 +615,9 @@ namespace Dune
template<class M, class X, class TM, class TA>
void SeqOverlappingSchwarz<M,X,TM,TA>::apply(X& x, const X& b)
{
- this->template apply<true>(x, b);
+ Applier<SeqOverlappingSchwarz>::apply(*this, x, b);
}
+
template<class M, class X, class TM, class TA>
template<bool forward>
void SeqOverlappingSchwarz<M,X,TM,TA>::apply(X& x, const X& b)
@@ -600,7 +638,6 @@ namespace Dune
X v(x); // temporary for the update
v=0;
- //x=1;
typedef typename AdderSelector<TM,X>::Adder Adder;
for(iterator solver=ApplyHelper<solver_vector,subdomain_vector,forward>::begin(solvers);
diff --git a/istl/test/laplacian.hh b/istl/test/laplacian.hh
index ca2d5ebe3..270f2cf1a 100644
--- a/istl/test/laplacian.hh
+++ b/istl/test/laplacian.hh
@@ -53,19 +53,54 @@ void setupLaplacian(Dune::BCRSMatrix<B>& A, int N)
int x = i.index()%N; // x coordinate in the 2d field
int y = i.index()/N; // y coordinate in the 2d field
- i->operator[](i.index())=diagonal;
+ /* if(x==0 || x==N-1 || y==0||y==N-1){
- if(y>0)
- i->operator[](i.index()-N)=bone;
+ i->operator[](i.index())=1.0;
- if(y<N-1)
- i->operator[](i.index()+N)=bone;
+ if(y>0)
+ i->operator[](i.index()-N)=0;
- if(x>0)
- i->operator[](i.index()-1)=bone;
+ if(y<N-1)
+ i->operator[](i.index()+N)=0.0;
+
+ if(x>0)
+ i->operator[](i.index()-1)=0.0;
+
+ if(x < N-1)
+ i->operator[](i.index()+1)=0.0;
+
+ }else*/
+ {
+
+ i->operator[](i.index())=diagonal;
- if(x < N-1)
- i->operator[](i.index()+1)=bone;
+ if(y>0)
+ i->operator[](i.index()-N)=bone;
+
+ if(y<N-1)
+ i->operator[](i.index()+N)=bone;
+
+ if(x>0)
+ i->operator[](i.index()-1)=bone;
+
+ if(x < N-1)
+ i->operator[](i.index()+1)=bone;
+ }
+ }
+}
+template<int BS>
+void setBoundary(Dune::BlockVector<Dune::FieldVector<double,BS> >& lhs,
+ Dune::BlockVector<Dune::FieldVector<double,BS> >& rhs,
+ int N)
+{
+ for(int i=0; i < lhs.size(); ++i) {
+ int x = i/N;
+ int y = i%N;
+
+ if(x==0 || y ==0 || x==N-1 || y==N-1) {
+ double h = 1.0 / ((double) (N-1));
+ lhs[i]=rhs[i]=0; //((double)x)*((double)y)*h*h;
+ }
}
}
#endif
diff --git a/istl/test/overlappingschwarztest.cc b/istl/test/overlappingschwarztest.cc
index c662e4b51..107079aa9 100644
--- a/istl/test/overlappingschwarztest.cc
+++ b/istl/test/overlappingschwarztest.cc
@@ -1,18 +1,20 @@
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#include "config.h"
+#include <dune/istl/io.hh>
#include <dune/istl/bvector.hh>
#include <dune/istl/operators.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <laplacian.hh>
#include <dune/common/timer.hh>
+#include <dune/common/sllist.hh>
#include <dune/istl/overlappingschwarz.hh>
int main(int argc, char** argv)
{
const int BS=1;
- int N=100;
+ int N=4;
if(argc>1)
N = atoi(argv[1]);
@@ -30,17 +32,26 @@ int main(int argc, char** argv)
Vector b(N*N), x(N*N);
setupLaplacian(mat,N);
- b=1;
- x=0;
-
+ b=0;
+ x=100;
+ //setBoundary(x,b,N);
+ /*
+ for (BCRSMat::RowIterator row = mat.begin(); row != mat.end(); ++row)
+ (*row)[row.index()]+=row.index();
+ */
// create the subdomains
- int domainSize=4;
- int overlap = 1;
+ int domainSize=2;
+ if(argc>2)
+ domainSize = atoi(argv[2]);
+ int overlap = 0;
int domainsPerDim=(N+domainSize-1)/domainSize;
// set up the overlapping domains
- std::vector<std::set<BCRSMat::size_type> > domains(domainsPerDim*domainsPerDim);
+ typedef Dune::SeqOverlappingSchwarz<BCRSMat,Vector> Schwarz;
+ typedef Schwarz::subdomain_vector subdomain_vector;
+
+ subdomain_vector domains(domainsPerDim*domainsPerDim);
for(int j=0; j < N; ++j)
@@ -72,27 +83,55 @@ int main(int argc, char** argv)
domains[domain*domainsPerDim+xdomain].insert(j*N+i);
}
- typedef std::vector<std::set<BCRSMat::size_type> >::const_iterator iterator;
- int i=0;
- for(iterator iter=domains.begin(); iter != domains.end(); ++iter) {
- typedef std::set<BCRSMat::size_type>::const_iterator entry_iterator;
- std::cout<<"domain "<<i++<<":";
- for(entry_iterator entry = iter->begin(); entry != iter->end(); ++entry) {
- std::cout<<" "<<*entry;
+ typedef subdomain_vector::const_iterator iterator;
+
+ if(N<10) {
+ int i=0;
+ for(iterator iter=domains.begin(); iter != domains.end(); ++iter) {
+ typedef iterator::value_type::const_iterator entry_iterator;
+ std::cout<<"domain "<<i++<<":";
+ for(entry_iterator entry = iter->begin(); entry != iter->end(); ++entry) {
+ std::cout<<" "<<*entry;
+ }
+ std::cout<<std::endl;
}
- std::cout<<std::endl;
+ Dune::printmatrix(std::cout, mat, std::string("A"), std::string("A"));
+ Dune::printvector(std::cout, b, std::string("B"), std::string("B"));
+ Dune::printvector(std::cout, x, std::string("X"), std::string("X"));
}
- Dune::SeqOverlappingSchwarz<BCRSMat,Vector> prec(mat, domains);
+ Dune::SeqOverlappingSchwarz<BCRSMat,Vector> prec(mat, domains, 1);
Dune::Timer watch;
watch.reset();
- Dune::LoopSolver<Vector> solver(fop, prec, 10e-8,80,2);
+ Dune::LoopSolver<Vector> solver(fop, prec, 1e-2,100,2);
Dune::InverseOperatorResult res;
+
+ // b=0;
+ // x=100;
+ // setBoundary(x,b,N);
+ std::cout<<"Additive Schwarz"<<std::endl;
solver.apply(x,b, res);
+ Dune::SeqOverlappingSchwarz<BCRSMat,Vector,Dune::MultiplicativeSchwarzMode> prec1(mat,domains, 1);
+ Dune::LoopSolver<Vector> solver1(fop, prec1, 1e-2,100,2);
+
+ b=0;
+ x=100;
+ //setBoundary(x,b,N);
+ std::cout << "Multiplicative Schwarz"<<std::endl;
+
+ solver1.apply(x,b, res);
+
+ Dune::SeqSOR<BCRSMat,Vector,Vector> sor(mat, 1,1);
+ Dune::LoopSolver<Vector> solver2(fop, sor, 1e-2,100,2);
+ b=0;
+ x=100;
+ //setBoundary(x,b,N);
+ std::cout << "SOR"<<std::endl;
+ solver2.apply(x,b, res);
}
--
GitLab