diff --git a/istl/paamg/amg.hh b/istl/paamg/amg.hh
index a8eab6068d8f13fb49b4db6b7f3b2b358108e4c7..34ac5667c933a3036d3d12be97fd648cf0d5b838 100644
--- a/istl/paamg/amg.hh
+++ b/istl/paamg/amg.hh
@@ -1,9 +1,10 @@
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
-// $Id : $
+// $Id$
#ifndef DUNE_AMG_AMG_HH
#define DUNE_AMG_AMG_HH
+#include <memory>
#include <dune/istl/paamg/smoother.hh>
#include <dune/istl/paamg/transfer.hh>
#include <dune/istl/paamg/hierarchy.hh>
@@ -23,7 +24,7 @@ namespace Dune
* @author Markus Blatt
* @brief The AMG preconditioner.
*/
- template<class M, class X, class Y, class CS, class S, class A>
+ template<class M, class X, class Y, class CS, class S, class A=std::allocator<X> >
class AMG : public Preconditioner<X,Y>
{
public:
@@ -65,14 +66,16 @@ namespace Dune
void post(Domain& x);
+ private:
+ /** @brief Multigrid cycle on a level. */
void mgc(typename Hierarchy<Smoother,A>::Iterator& smoother,
- const typename MatrixHierarchy::Iterator& matrix,
+ typename MatrixHierarchy::ParallelMatrixHierarchy::ConstIterator& matrix,
+ typename MatrixHierarchy::AggregatesMapList::const_iterator& aggregates,
typename Hierarchy<Domain,A>::Iterator& lhs,
typename Hierarchy<Domain,A>::Iterator& update,
typename Hierarchy<Range,A>::Iterator& rhs,
typename Hierarchy<Range,A>::Iterator& defect);
- private:
/** @brief The matrix we solve. */
const MatrixHierarchy* matrices_;
/** @brief The arguments to construct the smoother */
@@ -80,7 +83,7 @@ namespace Dune
/** @brief The hierarchy of the smoothers. */
Hierarchy<Smoother,A> smoothers_;
/** @brief The solver of the coarsest level. */
- const CoarseSolver* solver_;
+ CoarseSolver* solver_;
/** @brief The right hand side of our problem. */
Hierarchy<Range,A>* rhs_;
/** @brief The current defect. */
@@ -89,28 +92,24 @@ namespace Dune
Hierarchy<Domain,A>* lhs_;
/** @brief The vectors to store the update in. */
Hierarchy<Domain,A>* update_;
+ /** @brief Gamma, 1 for V-cycle and 2 for W-cycle. */
std::size_t gamma_;
+ /** @brief The number of pre and postsmoothing steps. */
std::size_t steps_;
+ std::size_t level;
};
template<class M, class X, class Y, class CS, class S, class A>
AMG<M,X,Y,CS,S,A>::AMG(const MatrixHierarchy& matrices, CoarseSolver& coarseSolver,
const SmootherArgs& smootherArgs,
std::size_t gamma, std::size_t smoothingSteps)
- : matrices_(&matrices), solver_(coarseSolver), smootherArgs_(smootherArgs),
- smoothers_(), solver_(0), rhs_(), defect_(), lhs_(), update_(), gamma_(gamma),
+ : matrices_(&matrices), smootherArgs_(smootherArgs),
+ smoothers_(), solver_(&coarseSolver), gamma_(gamma),
steps_(smoothingSteps)
{
- assert(matrices_.built());
- // Initialize smoother on the finest level if necessary
- if(smoothers_.levels()==0) {
- typename SmootherTraits<S>::Arguments args;
- args.setMatrix(matrices_.finest()->matrix());
- args.setArgs(smootherArgs_);
- smoothers_.addCoarser(args);
- }
+ assert(matrices_->isBuilt());
// build the necessary smoother hierarchies
- matrices_.coarsenSmoother(smoothers_, smootherArgs_);
+ matrices_->coarsenSmoother(smoothers_, smootherArgs_);
}
@@ -118,17 +117,17 @@ namespace Dune
void AMG<M,X,Y,CS,S,A>::pre(X& x, Y& b)
{
Y* copy = new Y(b);
- rhs_ = new Hierarchy<Range,A>(copy);
+ rhs_ = new Hierarchy<Range,A>(*copy);
copy = new Y(b);
- defect_ = new Hierarchy<Range,A>(copy);
+ defect_ = new Hierarchy<Range,A>(*copy);
copy = new X(x);
- lhs_ = new Hierarchy<Domain,A>(copy);
+ lhs_ = new Hierarchy<Domain,A>(*copy);
copy = new X(x);
- update_ = new Hierarchy<Domain,A>(copy);
- matrices_.coarsenVector(*rhs_);
- matrices_.coarsenVector(*defect_);
- matrices_.coarsenVector(*lhs_);
- matrices_.coarsenVector(*update_);
+ update_ = new Hierarchy<Domain,A>(*copy);
+ matrices_->coarsenVector(*rhs_);
+ matrices_->coarsenVector(*defect_);
+ matrices_->coarsenVector(*lhs_);
+ matrices_->coarsenVector(*update_);
// Preprocess all smoothers
typedef typename Hierarchy<Smoother,A>::Iterator Iterator;
@@ -136,13 +135,13 @@ namespace Dune
typedef typename Hierarchy<Domain,A>::Iterator DIterator;
Iterator coarsest = smoothers_.coarsest();
Iterator smoother = smoothers_.finest();
- RIterator rhs = rhs_.finest();
- DIterator lhs = lhs_.finest();
+ RIterator rhs = rhs_->finest();
+ DIterator lhs = lhs_->finest();
- smoother.pre(*lhs,*rhs);
+ smoother->pre(*lhs,*rhs);
for(++smoother, ++lhs, ++rhs; smoother != coarsest; ++smoother, ++lhs, ++rhs)
- smoother.pre(*lhs,*rhs);
+ smoother->pre(*lhs,*rhs);
}
template<class M, class X, class Y, class CS, class S, class A>
@@ -151,32 +150,36 @@ namespace Dune
*(lhs_->finest())=v;
*(rhs_->finest())=d;
- //Create the iterator over the levels
+ typename Hierarchy<Smoother,A>::Iterator smoother = smoothers_.finest();
+ typename MatrixHierarchy::ParallelMatrixHierarchy::ConstIterator matrix = matrices_->matrices().finest();
+ typename MatrixHierarchy::AggregatesMapList::const_iterator aggregates = matrices_->aggregatesMaps().begin();
+ typename Hierarchy<Domain,A>::Iterator lhs = lhs_->finest();
+ typename Hierarchy<Domain,A>::Iterator update = update_->finest();
+ typename Hierarchy<Range,A>::Iterator rhs = rhs_->finest();
+ typename Hierarchy<Range,A>::Iterator defect = defect_->finest();
- // Coarsest level
- typename Hierarchy<Domain>::Iterator coarsest =
- mgc(smoothers_.finest, matrices_.finest(), lhs_.finest(),
- update_.finest(), rhs_.finest(), defect_.finest());
+ level=0;
+ mgc(smoother, matrix, aggregates, lhs, update, rhs, defect);
}
template<class M, class X, class Y, class CS, class S, class A>
void AMG<M,X,Y,CS,S,A>::mgc(typename Hierarchy<Smoother,A>::Iterator& smoother,
- const typename MatrixHierarchy::Iterator& matrix,
+ typename MatrixHierarchy::ParallelMatrixHierarchy::ConstIterator& matrix,
+ typename MatrixHierarchy::AggregatesMapList::const_iterator& aggregates,
typename Hierarchy<Domain,A>::Iterator& lhs,
typename Hierarchy<Domain,A>::Iterator& update,
typename Hierarchy<Range,A>::Iterator& rhs,
- typename Hierarchy<Range,A>::Iterator& defect)
- {
- if(matrix == matrices_.coarsest()) {
+ typename Hierarchy<Range,A>::Iterator& defect){
+ if(matrix == matrices_->matrices().coarsest()) {
// Solve directly
InverseOperatorResult res;
- solver_.apply(*lhs, rhs, res);
+ solver_->apply(*lhs, *rhs, res);
}else{
// presmoothing
- for(int i=0; i < steps_; ++i) {
+ for(std::size_t i=0; i < steps_; ++i) {
// calculate defect
*defect = *rhs;
- matrix().matrix().mmv(*lhs, *defect);
+ matrix->matrix().mmv(*lhs, *defect);
// smooth
smoother->apply(*update, *rhs);
@@ -187,43 +190,46 @@ namespace Dune
// calculate defect
*defect = *rhs;
- matrix().matrix().mmv(*lhs, *defect);
+ matrix->matrix().mmv(*lhs, *defect);
//restrict defect to coarse level right hand side.
++rhs;
- Transfer<typename MatrixHierarchy::AggregatesMap::Vertex,Range>::restrict (matrix->aggregates(),
- rhs,
- defect);
+ Transfer<typename MatrixHierarchy::AggregatesMap::AggregateDescriptor,Range>
+ ::restrict (*(*aggregates), *rhs, *defect);
+
// prepare coarse system
- ++lhs=0;
+ ++lhs;
++update;
++defect;
++smoother;
++matrix;
+ ++aggregates;
+ ++level;
+ *lhs=0;
// next level
- mgc(smoother, matrix, lhs, update, rhs, defect);
-
+ mgc(smoother, matrix, aggregates, lhs, update, rhs, defect);
+ --level;
//prolongate (coarse x is the new update)
--matrix;
+ --aggregates;
--update;
- Transfer<typename MatrixHierarchy::AggregatesMap::Vertex,Range>::prolongate(matrix->aggregates(),
- update,
- lhs,
- 0.8);
+ Transfer<typename MatrixHierarchy::AggregatesMap::AggregateDescriptor,Range>
+ ::prolongate(*(*aggregates), *update, *lhs, 0.8);
+
--lhs;
--rhs;
--defect;
--smoother;
// postsmoothing
- for(int i=0; i < steps_; ++i) {
+ for(std::size_t i=0; i < steps_; ++i) {
// calculate defect
*defect = *rhs;
- matrix().matrix().mmv(*lhs, *defect);
+ matrix->matrix().mmv(*lhs, *defect);
// smooth
smoother->apply(*update, *defect);
@@ -243,20 +249,20 @@ namespace Dune
typedef typename Hierarchy<Domain,A>::Iterator DIterator;
Iterator coarsest = smoothers_.coarsest();
Iterator smoother = smoothers_.finest();
- DIterator lhs = lhs_.finest();
+ DIterator lhs = lhs_->finest();
- smoother.post(*lhs);
+ smoother->post(*lhs);
for(++smoother; smoother != coarsest; ++smoother, ++lhs)
- smoother.post(*lhs);
+ smoother->post(*lhs);
- delete &(*update_.finest());
+ delete &(*update_->finest());
delete update_;
- delete &(*lhs_.finest());
+ delete &(*lhs_->finest());
delete lhs_;
- delete &(*defect_.finest());
+ delete &(*defect_->finest());
delete defect_;
- delete &(*rhs_.finest());
+ delete &(*rhs_->finest());
delete rhs_;
}
} // end namespace Amg
diff --git a/istl/paamg/hierarchy.hh b/istl/paamg/hierarchy.hh
index 9b2a9560c3fae0c167e0c95b37be4cbca9ee67a1..7a9f8382b77b79d9bd5b7c83c71be92f61c27632 100644
--- a/istl/paamg/hierarchy.hh
+++ b/istl/paamg/hierarchy.hh
@@ -297,7 +297,11 @@ namespace Dune
typedef A Allocator;
/** @brief The type of the aggregates map we use. */
typedef AggregatesMap<typename MatrixGraph<M>::VertexDescriptor> AggregatesMap;
+ /** @brief The type of the parallel matrix hierarchy. */
+ typedef Hierarchy<ParallelMatrix,Allocator> ParallelMatrixHierarchy;
+ /** @brief The type of the aggregates maps list. */
+ typedef std::list<AggregatesMap*,Allocator> AggregatesMapList;
/**
* @brief Constructor
* @param fineMatrix The matrix to coarsen.
@@ -324,8 +328,8 @@ namespace Dune
template<class V, class TA>
void coarsenVector(Hierarchy<BlockVector<V,TA> >& hierarchy) const;
- template<class S>
- void coarsenSmoother(Hierarchy<S>& smoothers,
+ template<class S, class TA>
+ void coarsenSmoother(Hierarchy<S,TA>& smoothers,
const typename SmootherTraits<S>::Arguments& args) const;
/**
@@ -340,15 +344,16 @@ namespace Dune
*/
bool isBuilt() const;
+ const ParallelMatrixHierarchy& matrices() const;
+
+ const AggregatesMapList& aggregatesMaps() const;
+
private:
typedef typename ConstructionTraits<ParallelMatrix>::Arguments MatrixArgs;
- /** @brief The type of the aggregates maps list. */
- typedef std::list<AggregatesMap*,Allocator> AggregatesMapList;
/** @brief The list of aggregates maps. */
- AggregatesMapList aggregatesMaps;
- typedef Hierarchy<ParallelMatrix,Allocator> MMatrixHierarchy;
+ AggregatesMapList aggregatesMaps_;
/** @brief The hierarchy of parallel matrices. */
- MMatrixHierarchy matrices_;
+ ParallelMatrixHierarchy matrices_;
typedef Hierarchy<Interface,Allocator> InterfaceHierarchy;
/** @brief The hierarchy of interfaces. */
InterfaceHierarchy interfaces_;
@@ -361,7 +366,7 @@ namespace Dune
template<class T>
bool coarsenTargetReached(const T& crit,
- const typename MMatrixHierarchy::Iterator& matrix);
+ const typename ParallelMatrixHierarchy::Iterator& matrix);
};
template<class T>
@@ -428,7 +433,7 @@ namespace Dune
template<typename T>
inline bool
MatrixHierarchy<M,IS,O,A>::coarsenTargetReached(const T& crit,
- const typename MMatrixHierarchy::Iterator& matrix)
+ const typename ParallelMatrixHierarchy::Iterator& matrix)
{
int nodes = matrix->matrix().N();
int totalNodes;
@@ -444,7 +449,7 @@ namespace Dune
{
GalerkinProduct productBuilder;
int procs;
- typedef typename MMatrixHierarchy::Iterator MatIterator;
+ typedef typename ParallelMatrixHierarchy::Iterator MatIterator;
typedef typename InterfaceHierarchy::Iterator InterfaceIterator;
typedef typename CommunicatorHierarchy::Iterator CommIterator;
@@ -479,7 +484,7 @@ namespace Dune
PropertiesGraph pg(sg, IdentityMap(), sg.getEdgeIndexMap());
AggregatesMap* aggregatesMap=new AggregatesMap(pg.maxVertex());
- aggregatesMaps.push_back(aggregatesMap);
+ aggregatesMaps_.push_back(aggregatesMap);
int noAggregates = aggregatesMap->buildAggregates(mlevel->matrix(), pg, criterion);
@@ -536,16 +541,31 @@ namespace Dune
matrices_.addCoarser(MatrixArgs(*coarseMatrix, *coarseIndices, *coarseRemote));
}
built_=true;
+ AggregatesMap* aggregatesMap=new AggregatesMap(0);
+ aggregatesMaps_.push_back(aggregatesMap);
+ }
+ template<class M, class IS, class R, class I>
+ const typename MatrixHierarchy<M,IS,R,I>::ParallelMatrixHierarchy&
+ MatrixHierarchy<M,IS,R,I>::matrices() const
+ {
+ return matrices_;
}
+
+ template<class M, class IS, class R, class I>
+ const typename MatrixHierarchy<M,IS,R,I>::AggregatesMapList&
+ MatrixHierarchy<M,IS,R,I>::aggregatesMaps() const
+ {
+ return aggregatesMaps_;
+ }
template<class M, class IS, class R, class I>
MatrixHierarchy<M,IS,R,I>::~MatrixHierarchy()
{
typedef typename AggregatesMapList::reverse_iterator AggregatesMapIterator;
- typedef typename MMatrixHierarchy::Iterator Iterator;
+ typedef typename ParallelMatrixHierarchy::Iterator Iterator;
- AggregatesMapIterator amap = aggregatesMaps.rbegin();
+ AggregatesMapIterator amap = aggregatesMaps_.rbegin();
for(Iterator level=matrices_.coarsest(), finest=matrices_.finest(); level != finest; --level, ++amap) {
delete *amap;
@@ -567,7 +587,7 @@ namespace Dune
void MatrixHierarchy<M,IS,R,I>::coarsenVector(Hierarchy<BlockVector<V,TA> >& hierarchy) const
{
assert(hierarchy.levels()==1);
- typedef typename MMatrixHierarchy::ConstIterator Iterator;
+ typedef typename ParallelMatrixHierarchy::ConstIterator Iterator;
Iterator coarsest = matrices_.coarsest();
int level=0;
std::cout<<"Level "<<level<<" has "<<matrices_.finest()->matrix().N()<<" unknows!"<<std::endl;
@@ -581,16 +601,18 @@ namespace Dune
}
template<class M, class IS, class R, class I>
- template<class S>
- void MatrixHierarchy<M,IS,R,I>::coarsenSmoother(Hierarchy<S>& smoothers,
+ template<class S, class TA>
+ void MatrixHierarchy<M,IS,R,I>::coarsenSmoother(Hierarchy<S,TA>& smoothers,
const typename SmootherTraits<S>::Arguments& sargs) const
{
- assert(smoothers.levels()==1);
- typedef typename MMatrixHierarchy::ConstIterator Iterator;
+ assert(smoothers.levels()==0);
+ typedef typename ParallelMatrixHierarchy::ConstIterator Iterator;
typename ConstructionTraits<S>::Arguments cargs;
cargs.setArgs(sargs);
Iterator coarsest = matrices_.coarsest();
+ int level=0;
for(Iterator matrix = matrices_.finest(); matrix != coarsest; ++matrix) {
+ std::cout<<"level "<<level++<<" "<<&matrix->matrix();
cargs.setMatrix(matrix->matrix());
smoothers.addCoarser(cargs);
}
@@ -600,9 +622,9 @@ namespace Dune
void MatrixHierarchy<M,IS,R,I>::recalculateGalerkin()
{
typedef typename AggregatesMapList::iterator AggregatesMapIterator;
- typedef typename MMatrixHierarchy::Iterator Iterator;
+ typedef typename ParallelMatrixHierarchy::Iterator Iterator;
- AggregatesMapIterator amap = aggregatesMaps.begin();
+ AggregatesMapIterator amap = aggregatesMaps_.begin();
GalerkinProduct productBuilder;
for(Iterator level = matrices_.finest(), coarsest=matrices_.coarsest(); level!=coarsest; ++amap) {
diff --git a/istl/paamg/smoother.hh b/istl/paamg/smoother.hh
index 06d96c39de54efce9e97fcb6e9a8e473e5d93e92..39a0e3aad7ad5126e86e52b53d4781d4c2018def 100644
--- a/istl/paamg/smoother.hh
+++ b/istl/paamg/smoother.hh
@@ -3,6 +3,7 @@
#ifndef DUNE_AMGSMOOTHER_HH
#define DUNE_AMGSMOOTHER_HH
+#include <iostream>
#include <dune/istl/paamg/construction.hh>
#include <dune/istl/preconditioners.hh>
namespace Dune
@@ -24,13 +25,13 @@ namespace Dune
/**
* @brief The default class for the smoother arguments.
*/
- template<class M>
+ template<class T>
struct DefaultSmootherArgs
{
/**
* @brief The type of matrix the smoother is for.
*/
- typedef M Matrix;
+ typedef typename T::matrix_type Matrix;
/**
* @brief The numbe of iterations to perform.
@@ -40,7 +41,7 @@ namespace Dune
/**
* @brief The relaxation factor to use.
*/
- typename M::field_type relaxationFactor;
+ typename Matrix::field_type relaxationFactor;
/**
* @brief Default constructor.
@@ -76,7 +77,9 @@ namespace Dune
public:
void setMatrix(const Matrix& matrix)
{
+ std::cout<<"Setting matrix "<<&matrix<<" size="<<matrix.N()<<std::endl;
matrix_=&matrix;
+ std::cout<<"Matrix "<<matrix_<<" size="<<matrix_->N()<<std::endl;
}
void setArgs(const DefaultSmootherArgs<T>& args)
diff --git a/istl/paamg/test/Makefile.am b/istl/paamg/test/Makefile.am
index ebd38715ec90f615eb6fecfc03a41f6aaee7865f..23e190a3d5143ca7b270e98a21200de25346dfa6 100644
--- a/istl/paamg/test/Makefile.am
+++ b/istl/paamg/test/Makefile.am
@@ -3,7 +3,7 @@
#if MPI
MPITESTS = graphtest
TESTPROGS = galerkintest hierarchytest amgtest
- TESTRUNS = run-galerkin run-hierarchy
+ TESTRUNS = run-galerkin run-hierarchy run-amg
#endif
NORMALTESTS = $(MPITESTS)
diff --git a/istl/paamg/test/amgtest.cc b/istl/paamg/test/amgtest.cc
index 825f82b7e4422d1b3b2452a955d3b7281513e457..d99ada7c5e198c90c2ea614199d7f0870b17f3bd 100644
--- a/istl/paamg/test/amgtest.cc
+++ b/istl/paamg/test/amgtest.cc
@@ -1,7 +1,74 @@
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#include "config.h"
+#include "anisotropic.hh"
#include <dune/istl/paamg/amg.hh>
+#include <dune/istl/indexset.hh>
-int main()
-{}
+int main(int argc, char** argv)
+{
+ MPI_Init(&argc, &argv);
+
+ const int BS=1, N=8;
+
+ int procs, rank;
+ MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+ MPI_Comm_size(MPI_COMM_WORLD, &procs);
+
+ IndexSet indices;
+ typedef Dune::FieldMatrix<double,BS,BS> MatrixBlock;
+ typedef Dune::BCRSMatrix<MatrixBlock> BCRSMat;
+ typedef Dune::FieldVector<double,BS> VectorBlock;
+ typedef Dune::BlockVector<VectorBlock> Vector;
+
+ int n;
+
+ BCRSMat mat = setupAnisotropic2d<N,BS>(indices, &n);
+ Vector b(indices.size()), x(indices.size());
+
+ b=1;
+ x=100;
+
+ RemoteIndices remoteIndices(indices,indices,MPI_COMM_WORLD);
+ remoteIndices.rebuild<false>();
+
+ typedef Dune::Interface<IndexSet> Interface;
+
+ Interface interface;
+
+ typedef Dune::EnumItem<GridFlag,overlap> OverlapFlags;
+ typedef Dune::Amg::MatrixHierarchy<BCRSMat,IndexSet,OverlapFlags> MHierarchy;
+ typedef Dune::Amg::Hierarchy<Vector> VHierarchy;
+
+ interface.build(remoteIndices, Dune::NegateSet<OverlapFlags>(), OverlapFlags());
+
+ MHierarchy hierarchy(mat, indices, remoteIndices, interface);
+ VHierarchy vh(b);
+
+ typedef Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<BCRSMat,Dune::Amg::FirstDiagonal> >
+ Criterion;
+
+ Criterion criterion(10,4);
+
+ hierarchy.build(criterion);
+ hierarchy.coarsenVector(vh);
+
+ typedef Dune::SeqSSOR<BCRSMat,Vector,Vector> Smoother;
+ typedef Dune::Amg::SmootherTraits<Smoother>::Arguments SmootherArgs;
+ SmootherArgs smootherArgs;
+
+ Dune::MatrixAdapter<BCRSMat,Vector,Vector> op(hierarchy.matrices().coarsest()->matrix());
+ Dune::SeqSSOR<BCRSMat,Vector,Vector> ssor(hierarchy.matrices().coarsest()->matrix(),1,1.0);
+ typedef Dune::LoopSolver<Vector> Solver;
+ Solver solver(op,ssor,1E-6,8000,2);
+ Dune::SeqScalarProduct<Vector> sp;
+ typedef Dune::Amg::AMG<MHierarchy,Vector,Vector,Solver,Smoother> AMG;
+
+ AMG amg(hierarchy, solver, smootherArgs, 1, 1);
+
+ amg.pre(x, b);
+ amg.apply(x,b);
+ amg.post(x);
+
+ MPI_Finalize();
+}
diff --git a/istl/paamg/test/run-amg b/istl/paamg/test/run-amg
new file mode 100755
index 0000000000000000000000000000000000000000..8ef12937f0b05796b2cab62c59324e9f68bb7ce9
--- /dev/null
+++ b/istl/paamg/test/run-amg
@@ -0,0 +1,17 @@
+#!/bin/sh
+# $Id$
+
+# run lamboot if possible
+if which lamboot; then
+ lamboot 2>&1 > /dev/null
+fi
+
+# start test and abort on fail
+set -e
+mpirun -np 1 amgtest
+set +e
+
+# possibly stop LAM again
+if which lamhalt; then
+ lamhalt 2>&1 > /dev/null
+fi
diff --git a/istl/paamg/transfer.hh b/istl/paamg/transfer.hh
index b14d2e8985b7768c7e8872a94f8dcee50494b0d7..df07786b89538d4ef75cbca8b8a5db4bf84c5e80 100644
--- a/istl/paamg/transfer.hh
+++ b/istl/paamg/transfer.hh
@@ -28,20 +28,20 @@ namespace Dune
typedef V1 Vertex;
typedef V2 Vector;
- void prolongate(const AggregatesMap<Vertex>& aggregates, const Vector & coarse, Vector fine,
- typename Vector::field_type damp) const;
- void restrict (const AggregatesMap<Vertex>& aggregates, Vector& coarse, const Vector fine) const;
+ static void prolongate(const AggregatesMap<Vertex>& aggregates, const Vector & coarse, Vector fine,
+ typename Vector::field_type damp);
+ static void restrict (const AggregatesMap<Vertex>& aggregates, Vector& coarse, const Vector fine);
};
template<class V1, class V2>
void Transfer<V1,V2>::prolongate(const AggregatesMap<Vertex>& aggregates, const Vector& coarse,
- Vector fine, typename Vector::field_type damp) const
+ Vector fine, typename Vector::field_type damp)
{
DUNE_THROW(NotImplemented, "There is no secialization available for this type of vector!");
}
template<class V1, class V2>
void Transfer<V1,V2>::restrict (const AggregatesMap<Vertex>& aggregates, Vector& coarse,
- const Vector fine) const
+ const Vector fine)
{
DUNE_THROW(NotImplemented, "There is no secialization available for this type of vector!");
}
@@ -52,36 +52,36 @@ namespace Dune
public:
typedef V Vertex;
typedef BlockVector<B> Vector;
- void prolongate(const AggregatesMap<Vertex>& aggregates, const Vector & coarse, Vector fine,
- typename Vector::field_type damp) const;
+ static void prolongate(const AggregatesMap<Vertex>& aggregates, const Vector & coarse, Vector fine,
+ typename Vector::field_type damp);
- void restrict (const AggregatesMap<Vertex>& aggregates, Vector& coarse, const Vector fine) const;
+ static void restrict (const AggregatesMap<Vertex>& aggregates, Vector& coarse, const Vector fine);
};
template<class V, class B>
- void Transfer<V,BlockVector<B> >::prolongate(const AggregatesMap<Vertex>& aggregates,
- const Vector& coarse, Vector fine,
- typename Vector::field_type damp) const
+ inline void Transfer<V,BlockVector<B> >::prolongate(const AggregatesMap<Vertex>& aggregates,
+ const Vector& coarse, Vector fine,
+ typename Vector::field_type damp)
{
typedef typename Vector::iterator Iterator;
Iterator end = fine.end();
for(Iterator block=fine.begin(); block != end; ++block) {
- block = coarse[aggregates[block.index()]];
- block *= damp;
+ *block = coarse[aggregates[block.index()]];
+ *block *= damp;
}
}
template<class V, class B>
- void Transfer<V,BlockVector<B> >::restrict (const AggregatesMap<Vertex>& aggregates,
- Vector& coarse,
- const Vector fine) const
+ inline void Transfer<V,BlockVector<B> >::restrict (const AggregatesMap<Vertex>& aggregates,
+ Vector& coarse,
+ const Vector fine)
{
// Set coarse vector to zero
coarse=0;
- typedef typename Vector::iterator Iterator;
+ typedef typename Vector::const_iterator Iterator;
Iterator end = fine.end();
for(Iterator block=fine.begin(); block != end; ++block)
coarse[aggregates[block.index()]] += *block;