Skip to content
Snippets Groups Projects
Commit 9348a968 authored by Markus Blatt's avatar Markus Blatt
Browse files

[Merge][Feature] from speedup_amg_build branch 

[[Imported from SVN: r1900]]
parents a388fca7 1192d2f5 57a0e068
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 1166 additions and 262 deletions
dune/istl/paamg/aggregates.hh
+ 282
− 228
View file @ 9348a968
...@@ -114,6 +114,129 @@ namespace Dune ...@@ -114,6 +114,129 @@ namespace Dune
return os; return os;
} }
/**
* @brief Dependency policy for symmetric matrices.
*
* We assume that not only the sparsity pattern is symmetric
* but also the entries (a_ij=aji). If that is not the case
* the resulting dependency graph might be unsymmetric.
*
* \tparam M The type of the matrix
* \tparam N The type of the metric that turns matrix blocks into
* field values
*/
template<class M, class N>
class SymmetricMatrixDependency : public Dune::Amg::Parameters
{
public:
/**
* @brief The matrix type we build the dependency of.
*/
typedef M Matrix;
/**
* @brief The norm to use for examining the matrix entries.
*/
typedef N Norm;
/**
* @brief Constant Row iterator of the matrix.
*/
typedef typename Matrix::row_type Row;
/**
* @brief Constant column iterator of the matrix.
*/
typedef typename Matrix::ConstColIterator ColIter;
void init(const Matrix* matrix);
void initRow(const Row& row, int index);
void examine(const ColIter& col);
template<class G>
void examine(G& graph, const typename G::EdgeIterator& edge, const ColIter& col);
bool isIsolated();
SymmetricMatrixDependency(const Parameters& parms)
: Parameters(parms)
{}
SymmetricMatrixDependency()
: Parameters()
{}
protected:
/** @brief The matrix we work on. */
const Matrix* matrix_;
/** @brief The current max value.*/
typename Matrix::field_type maxValue_;
/** @brief The functor for calculating the norm. */
Norm norm_;
/** @brief index of the currently evaluated row. */
int row_;
/** @brief The norm of the current diagonal. */
typename Matrix::field_type diagonal_;
std::vector<typename Matrix::field_type> vals_;
typename std::vector<typename Matrix::field_type>::iterator valIter_;
};
template<class M, class N>
inline void SymmetricMatrixDependency<M,N>::init(const Matrix* matrix)
{
matrix_ = matrix;
}
template<class M, class N>
inline void SymmetricMatrixDependency<M,N>::initRow(const Row& row, int index)
{
vals_.assign(row.size(), 0.0);
assert(vals_.size()==row.size());
valIter_=vals_.begin();
maxValue_ = std::min(- std::numeric_limits<typename Matrix::field_type>::max(), std::numeric_limits<typename Matrix::field_type>::min());
diagonal_=norm_(row[index]);
row_ = index;
}
template<class M, class N>
inline void SymmetricMatrixDependency<M,N>::examine(const ColIter& col)
{
// skip positive offdiagonals if norm preserves sign of them.
typename Matrix::field_type eij = norm_(*col);
if(!N::is_sign_preserving || eij<0) // || eji<0)
{
*valIter_ = eij/diagonal_*eij/norm_(matrix_->operator[](col.index())[col.index()]);
maxValue_ = std::max(maxValue_, *valIter_);
}else
*valIter_ =0;
++valIter_;
}
template<class M, class N>
template<class G>
inline void SymmetricMatrixDependency<M,N>::examine(G& graph, const typename G::EdgeIterator& edge, const ColIter& col)
{
if(*valIter_ > alpha() * maxValue_) {
edge.properties().setDepends();
edge.properties().setInfluences();
}
++valIter_;
}
template<class M, class N>
inline bool SymmetricMatrixDependency<M,N>::isIsolated()
{
if(diagonal_==0)
DUNE_THROW(Dune::ISTLError, "No diagonal entry for row "<<row_<<".");
valIter_=vals_.begin();
return maxValue_ < beta();
}
/** /**
* @brief Dependency policy for symmetric matrices. * @brief Dependency policy for symmetric matrices.
*/ */
...@@ -618,13 +741,6 @@ namespace Dune ...@@ -618,13 +741,6 @@ namespace Dune
*/ */
typedef PoolAllocator<Vertex,100> Allocator; typedef PoolAllocator<Vertex,100> Allocator;
/**
* @brief The type of a single linked list of vertex
* descriptors.
*/
typedef SLList<Vertex,Allocator> VertexList;
/** /**
* @brief The type of a single linked list of vertex * @brief The type of a single linked list of vertex
* descriptors. * descriptors.
...@@ -632,7 +748,7 @@ namespace Dune ...@@ -632,7 +748,7 @@ namespace Dune
typedef S VertexSet; typedef S VertexSet;
/** @brief Const iterator over a vertex list. */ /** @brief Const iterator over a vertex list. */
typedef typename VertexList::const_iterator const_iterator; typedef typename VertexSet::const_iterator const_iterator;
/** /**
* @brief Type of the mapping of aggregate members onto distance spheres. * @brief Type of the mapping of aggregate members onto distance spheres.
...@@ -646,8 +762,13 @@ namespace Dune ...@@ -646,8 +762,13 @@ namespace Dune
* @param connectivity The set of vertices connected to the aggregate. * @param connectivity The set of vertices connected to the aggregate.
* distance spheres. * distance spheres.
*/ */
Aggregate(const MatrixGraph& graph, AggregatesMap<Vertex>& aggregates, Aggregate(MatrixGraph& graph, AggregatesMap<Vertex>& aggregates,
VertexSet& connectivity); VertexSet& connectivity, std::vector<Vertex>& front_);
void invalidate()
{
--id_;
}
/** /**
* @brief Reconstruct the aggregat from an seed node. * @brief Reconstruct the aggregat from an seed node.
...@@ -667,6 +788,7 @@ namespace Dune ...@@ -667,6 +788,7 @@ namespace Dune
*/ */
void add(const Vertex& vertex); void add(const Vertex& vertex);
void add(std::vector<Vertex>& vertex);
/** /**
* @brief Clear the aggregate. * @brief Clear the aggregate.
*/ */
...@@ -675,11 +797,11 @@ namespace Dune ...@@ -675,11 +797,11 @@ namespace Dune
/** /**
* @brief Get the size of the aggregate. * @brief Get the size of the aggregate.
*/ */
typename VertexList::size_type size(); typename VertexSet::size_type size();
/** /**
* @brief Get tne number of connections to other aggregates. * @brief Get tne number of connections to other aggregates.
*/ */
typename VertexList::size_type connectSize(); typename VertexSet::size_type connectSize();
/** /**
* @brief Get the id identifying the aggregate. * @brief Get the id identifying the aggregate.
...@@ -696,7 +818,7 @@ namespace Dune ...@@ -696,7 +818,7 @@ namespace Dune
/** /**
* @brief The vertices of the aggregate. * @brief The vertices of the aggregate.
*/ */
VertexList vertices_; VertexSet vertices_;
/** /**
* @brief The number of the currently referenced * @brief The number of the currently referenced
...@@ -707,7 +829,7 @@ namespace Dune ...@@ -707,7 +829,7 @@ namespace Dune
/** /**
* @brief The matrix graph the aggregates live on. * @brief The matrix graph the aggregates live on.
*/ */
const MatrixGraph& graph_; MatrixGraph& graph_;
/** /**
* @brief The aggregate mapping we build. * @brief The aggregate mapping we build.
...@@ -718,6 +840,8 @@ namespace Dune ...@@ -718,6 +840,8 @@ namespace Dune
* @brief The connections to other aggregates. * @brief The connections to other aggregates.
*/ */
VertexSet& connected_; VertexSet& connected_;
std::vector<Vertex>& front_;
}; };
/** /**
...@@ -806,7 +930,7 @@ namespace Dune ...@@ -806,7 +930,7 @@ namespace Dune
/** /**
* @brief The vertices of the current aggregate front. * @brief The vertices of the current aggregate front.
*/ */
VertexList front_; std::vector<Vertex> front_;
/** /**
* @brief The set of connected vertices. * @brief The set of connected vertices.
...@@ -830,11 +954,9 @@ namespace Dune ...@@ -830,11 +954,9 @@ namespace Dune
const Aggregator<G>& aggregatesBuilder, const Aggregator<G>& aggregatesBuilder,
const AggregatesMap<Vertex>& aggregates); const AggregatesMap<Vertex>& aggregates);
~Stack(); ~Stack();
bool push(const Vertex& v);
void fill();
Vertex pop(); Vertex pop();
private: private:
enum { N = 256000 }; enum { N = 1300000 };
/** @brief The graph we work on. */ /** @brief The graph we work on. */
const MatrixGraph& graph_; const MatrixGraph& graph_;
...@@ -844,15 +966,14 @@ namespace Dune ...@@ -844,15 +966,14 @@ namespace Dune
const AggregatesMap<Vertex>& aggregates_; const AggregatesMap<Vertex>& aggregates_;
/** @brief The current size. */ /** @brief The current size. */
int size_; int size_;
int maxSize_; Vertex maxSize_;
/** @brief The index of the top element. */ /** @brief The index of the top element. */
int head_; typename MatrixGraph::ConstVertexIterator begin_;
int filled_; typename MatrixGraph::ConstVertexIterator end_;
/** @brief The values on the stack. */ /** @brief The values on the stack. */
Vertex* vals_; Vertex* vals_;
void localPush(const Vertex& v);
}; };
friend class Stack; friend class Stack;
...@@ -1096,25 +1217,17 @@ namespace Dune ...@@ -1096,25 +1217,17 @@ namespace Dune
* @param front The list to store the front vertices in. * @param front The list to store the front vertices in.
* @param graph The matrix graph we work on. * @param graph The matrix graph we work on.
*/ */
FrontMarker(VertexList& front, MatrixGraph& graph); FrontMarker(std::vector<Vertex>& front, MatrixGraph& graph);
void operator()(const typename MatrixGraph::ConstEdgeIterator& edge); void operator()(const typename MatrixGraph::ConstEdgeIterator& edge);
private: private:
/** @brief The list of front vertices. */ /** @brief The list of front vertices. */
VertexList& front_; std::vector<Vertex>& front_;
/** @brief The matrix graph we work on. */ /** @brief The matrix graph we work on. */
MatrixGraph& graph_; MatrixGraph& graph_;
}; };
/**
* @brief Mark the front of the current aggregate.
*
* The front are the direct (unaggregated) neighbours of
* the aggregate vertices.
*/
void markFront(const AggregatesMap<Vertex>& aggregates);
/** /**
* @brief Unmarks all front vertices. * @brief Unmarks all front vertices.
*/ */
...@@ -1179,14 +1292,6 @@ namespace Dune ...@@ -1179,14 +1292,6 @@ namespace Dune
*/ */
bool admissible(const Vertex& vertex, const AggregateDescriptor& aggregate, const AggregatesMap<Vertex>& aggregates) const; bool admissible(const Vertex& vertex, const AggregateDescriptor& aggregate, const AggregatesMap<Vertex>& aggregates) const;
/**
* @brief Push the neighbours of the current aggregate on the stack.
*
* @param stack The stack to push them on.
* @param isolated If true only isolated vertices are push onto the stack.
*/
void seedFromFront(Stack& stack, bool isolated);
/** /**
* @brief The maximum distance of the vertex to any vertex in the * @brief The maximum distance of the vertex to any vertex in the
* current aggregate. * current aggregate.
...@@ -1331,19 +1436,20 @@ namespace Dune ...@@ -1331,19 +1436,20 @@ namespace Dune
} }
template<class G,class S> template<class G,class S>
Aggregate<G,S>::Aggregate(const MatrixGraph& graph, AggregatesMap<Vertex>& aggregates, Aggregate<G,S>::Aggregate(MatrixGraph& graph, AggregatesMap<Vertex>& aggregates,
VertexSet& connected) VertexSet& connected, std::vector<Vertex>& front)
: vertices_(), id_(-1), graph_(graph), aggregates_(aggregates), : vertices_(), id_(-1), graph_(graph), aggregates_(aggregates),
connected_(connected) connected_(connected), front_(front)
{} {}
template<class G,class S> template<class G,class S>
void Aggregate<G,S>::reconstruct(const Vertex& vertex) void Aggregate<G,S>::reconstruct(const Vertex& vertex)
{ {
vertices_.push_back(vertex); /*
typedef typename VertexList::const_iterator iterator; vertices_.push_back(vertex);
iterator begin = vertices_.begin(); typedef typename VertexList::const_iterator iterator;
iterator end = vertices_.end(); iterator begin = vertices_.begin();
iterator end = vertices_.end();*/
throw "Not yet implemented"; throw "Not yet implemented";
while(begin!=end) { while(begin!=end) {
...@@ -1360,7 +1466,7 @@ namespace Dune ...@@ -1360,7 +1466,7 @@ namespace Dune
vertices_.clear(); vertices_.clear();
connected_.insert(vertex); connected_.insert(vertex);
dvverb << " Inserting "<<vertex<<" size="<<connected_.size(); dvverb << " Inserting "<<vertex<<" size="<<connected_.size();
id_ = vertex; ++id_ ;
add(vertex); add(vertex);
} }
...@@ -1368,8 +1474,11 @@ namespace Dune ...@@ -1368,8 +1474,11 @@ namespace Dune
template<class G,class S> template<class G,class S>
inline void Aggregate<G,S>::add(const Vertex& vertex) inline void Aggregate<G,S>::add(const Vertex& vertex)
{ {
vertices_.push_back(vertex); vertices_.insert(vertex);
aggregates_[vertex]=id_; aggregates_[vertex]=id_;
if(front_.size())
front_.erase(std::lower_bound(front_.begin(), front_.end(), vertex));
typedef typename MatrixGraph::ConstEdgeIterator iterator; typedef typename MatrixGraph::ConstEdgeIterator iterator;
const iterator end = graph_.endEdges(vertex); const iterator end = graph_.endEdges(vertex);
...@@ -1377,8 +1486,58 @@ namespace Dune ...@@ -1377,8 +1486,58 @@ namespace Dune
dvverb << " Inserting "<<aggregates_[edge.target()]; dvverb << " Inserting "<<aggregates_[edge.target()];
connected_.insert(aggregates_[edge.target()]); connected_.insert(aggregates_[edge.target()]);
dvverb <<" size="<<connected_.size(); dvverb <<" size="<<connected_.size();
if(aggregates_[edge.target()]==AggregatesMap<Vertex>::UNAGGREGATED &&
!graph_.getVertexProperties(edge.target()).front())
{
front_.push_back(edge.target());
graph_.getVertexProperties(edge.target()).setFront();
}
} }
dvverb <<std::endl; dvverb <<std::endl;
std::sort(front_.begin(), front_.end());
}
template<class G,class S>
inline void Aggregate<G,S>::add(std::vector<Vertex>& vertices)
{
std::size_t oldsize = vertices_.size();
typedef typename std::vector<Vertex>::iterator Iterator;
typedef typename VertexSet::iterator SIterator;
SIterator pos=vertices_.begin();
std::vector<Vertex> newFront;
newFront.reserve(front_.capacity());
std::set_difference(front_.begin(), front_.end(), vertices.begin(), vertices.end(),
std::back_inserter(newFront));
front_=newFront;
for(Iterator vertex=vertices.begin(); vertex != vertices.end(); ++vertex)
{
pos=vertices_.insert(pos,*vertex);
vertices_.insert(*vertex);
graph_.getVertexProperties(*vertex).resetFront(); // Not a front node any more.
aggregates_[*vertex]=id_;
typedef typename MatrixGraph::ConstEdgeIterator iterator;
const iterator end = graph_.endEdges(*vertex);
for(iterator edge = graph_.beginEdges(*vertex); edge != end; ++edge) {
dvverb << " Inserting "<<aggregates_[edge.target()];
connected_.insert(aggregates_[edge.target()]);
if(aggregates_[edge.target()]==AggregatesMap<Vertex>::UNAGGREGATED &&
!graph_.getVertexProperties(edge.target()).front())
{
front_.push_back(edge.target());
graph_.getVertexProperties(edge.target()).setFront();
}
dvverb <<" size="<<connected_.size();
}
dvverb <<std::endl;
}
std::sort(front_.begin(), front_.end());
assert(oldsize+vertices.size()==vertices_.size());
} }
template<class G,class S> template<class G,class S>
inline void Aggregate<G,S>::clear() inline void Aggregate<G,S>::clear()
...@@ -1389,14 +1548,14 @@ namespace Dune ...@@ -1389,14 +1548,14 @@ namespace Dune
} }
template<class G,class S> template<class G,class S>
inline typename Aggregate<G,S>::VertexList::size_type inline typename Aggregate<G,S>::VertexSet::size_type
Aggregate<G,S>::size() Aggregate<G,S>::size()
{ {
return vertices_.size(); return vertices_.size();
} }
template<class G,class S> template<class G,class S>
inline typename Aggregate<G,S>::VertexList::size_type inline typename Aggregate<G,S>::VertexSet::size_type
Aggregate<G,S>::connectSize() Aggregate<G,S>::connectSize()
{ {
return connected_.size(); return connected_.size();
...@@ -1787,6 +1946,7 @@ namespace Dune ...@@ -1787,6 +1946,7 @@ namespace Dune
template<class G> template<class G>
std::size_t Aggregator<G>::distance(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates) std::size_t Aggregator<G>::distance(const Vertex& vertex, const AggregatesMap<Vertex>& aggregates)
{ {
return 0;
typename PropertyMapTypeSelector<VertexVisitedTag,G>::Type visitedMap = get(VertexVisitedTag(), *graph_); typename PropertyMapTypeSelector<VertexVisitedTag,G>::Type visitedMap = get(VertexVisitedTag(), *graph_);
VertexList vlist; VertexList vlist;
typename AggregatesMap<Vertex>::DummyEdgeVisitor dummy; typename AggregatesMap<Vertex>::DummyEdgeVisitor dummy;
...@@ -1796,7 +1956,7 @@ namespace Dune ...@@ -1796,7 +1956,7 @@ namespace Dune
} }
template<class G> template<class G>
inline Aggregator<G>::FrontMarker::FrontMarker(VertexList& front, MatrixGraph& graph) inline Aggregator<G>::FrontMarker::FrontMarker(std::vector<Vertex>& front, MatrixGraph& graph)
: front_(front), graph_(graph) : front_(front), graph_(graph)
{} {}
...@@ -1811,19 +1971,6 @@ namespace Dune ...@@ -1811,19 +1971,6 @@ namespace Dune
} }
} }
template<class G>
void Aggregator<G>::markFront(const AggregatesMap<Vertex>& aggregates)
{
assert(front_.size()==0);
FrontMarker frontBuilder(front_, *graph_);
typedef typename Aggregate<G,VertexSet>::const_iterator Iterator;
for(Iterator vertex=aggregate_->begin(); vertex != aggregate_->end(); ++vertex)
visitAggregateNeighbours(*vertex, AggregatesMap<Vertex>::UNAGGREGATED, aggregates, frontBuilder);
}
template<class G> template<class G>
inline bool Aggregator<G>::admissible(const Vertex& vertex, const AggregateDescriptor& aggregate, const AggregatesMap<Vertex>& aggregates) const inline bool Aggregator<G>::admissible(const Vertex& vertex, const AggregateDescriptor& aggregate, const AggregatesMap<Vertex>& aggregates) const
{ {
...@@ -1860,7 +2007,9 @@ namespace Dune ...@@ -1860,7 +2007,9 @@ namespace Dune
} }
} }
if(found) if(found)
{
return true; return true;
}
} }
} }
} }
...@@ -1893,7 +2042,9 @@ namespace Dune ...@@ -1893,7 +2042,9 @@ namespace Dune
} }
} }
if(found) if(found)
{
return true; return true;
}
} }
} }
} }
...@@ -1904,7 +2055,7 @@ namespace Dune ...@@ -1904,7 +2055,7 @@ namespace Dune
template<class G> template<class G>
void Aggregator<G>::unmarkFront() void Aggregator<G>::unmarkFront()
{ {
typedef typename VertexList::const_iterator Iterator; typedef typename std::vector<Vertex>::const_iterator Iterator;
for(Iterator vertex=front_.begin(); vertex != front_.end(); ++vertex) for(Iterator vertex=front_.begin(); vertex != front_.end(); ++vertex)
graph_->getVertexProperties(*vertex).resetFront(); graph_->getVertexProperties(*vertex).resetFront();
...@@ -2002,9 +2153,8 @@ namespace Dune ...@@ -2002,9 +2153,8 @@ namespace Dune
std::size_t maxFrontNeighbours=0; std::size_t maxFrontNeighbours=0;
Vertex candidate=AggregatesMap<Vertex>::UNAGGREGATED; Vertex candidate=AggregatesMap<Vertex>::UNAGGREGATED;
unmarkFront();
markFront(aggregates); typedef typename std::vector<Vertex>::const_iterator Iterator;
typedef typename VertexList::const_iterator Iterator;
for(Iterator vertex = front_.begin(); vertex != front_.end(); ++vertex) { for(Iterator vertex = front_.begin(); vertex != front_.end(); ++vertex) {
if(distance(*vertex, aggregates)>c.maxDistance()) if(distance(*vertex, aggregates)>c.maxDistance())
...@@ -2045,16 +2195,16 @@ namespace Dune ...@@ -2045,16 +2195,16 @@ namespace Dune
template<class C> template<class C>
void Aggregator<G>::growAggregate(const Vertex& seed, const AggregatesMap<Vertex>& aggregates, const C& c) void Aggregator<G>::growAggregate(const Vertex& seed, const AggregatesMap<Vertex>& aggregates, const C& c)
{ {
while(aggregate_->size() < c.minAggregateSize()) {
std::size_t distance_ =0;
while(aggregate_->size() < c.minAggregateSize()&& distance_<c.maxDistance()) {
int maxTwoCons=0, maxOneCons=0, maxNeighbours=-1; int maxTwoCons=0, maxOneCons=0, maxNeighbours=-1;
double maxCon=-1; double maxCon=-1;
Vertex candidate = AggregatesMap<Vertex>::UNAGGREGATED; std::vector<Vertex> candidates;
candidates.reserve(30);
unmarkFront();
markFront(aggregates);
typedef typename VertexList::const_iterator Iterator; typedef typename std::vector<Vertex>::const_iterator Iterator;
for(Iterator vertex = front_.begin(); vertex != front_.end(); ++vertex) { for(Iterator vertex = front_.begin(); vertex != front_.end(); ++vertex) {
// Only nonisolated nodes are considered // Only nonisolated nodes are considered
...@@ -2072,38 +2222,32 @@ namespace Dune ...@@ -2072,38 +2222,32 @@ namespace Dune
if(neighbours > maxNeighbours) { if(neighbours > maxNeighbours) {
maxNeighbours = neighbours; maxNeighbours = neighbours;
candidates.clear();
std::size_t distance_ = distance(*vertex, aggregates); candidates.push_back(*vertex);
}else{
if(c.maxDistance() >= distance_) { candidates.push_back(*vertex);
candidate = *vertex;
}
} }
}else if( con > maxCon) { }else if( con > maxCon) {
maxCon = con; maxCon = con;
maxNeighbours = noFrontNeighbours(*vertex); maxNeighbours = noFrontNeighbours(*vertex);
std::size_t distance_ = distance(*vertex, aggregates); candidates.clear();
candidates.push_back(*vertex);
if(c.maxDistance() >= distance_) {
candidate = *vertex;
}
} }
}else if(twoWayCons > maxTwoCons) { }else if(twoWayCons > maxTwoCons) {
maxTwoCons = twoWayCons; maxTwoCons = twoWayCons;
maxCon = connectivity(*vertex, aggregates); maxCon = connectivity(*vertex, aggregates);
maxNeighbours = noFrontNeighbours(*vertex); maxNeighbours = noFrontNeighbours(*vertex);
std::size_t distance_ = distance(*vertex, aggregates); candidates.clear();
candidates.push_back(*vertex);
if(c.maxDistance() >= distance_) {
candidate = *vertex;
}
// two way connections preceed // two way connections preceed
maxOneCons = std::numeric_limits<int>::max(); maxOneCons = std::numeric_limits<int>::max();
} }
if(twoWayCons > 0) if(twoWayCons > 0)
{
continue; // THis is a two-way node, skip tests for one way nodes continue; // THis is a two-way node, skip tests for one way nodes
}
/* The one way case */ /* The one way case */
int oneWayCons = oneWayConnections(*vertex, aggregate_->id(), aggregates); int oneWayCons = oneWayConnections(*vertex, aggregate_->id(), aggregates);
...@@ -2122,37 +2266,36 @@ namespace Dune ...@@ -2122,37 +2266,36 @@ namespace Dune
if(neighbours > maxNeighbours) { if(neighbours > maxNeighbours) {
maxNeighbours = neighbours; maxNeighbours = neighbours;
std::size_t distance_ = distance(*vertex, aggregates); candidates.clear();
candidates.push_back(*vertex);
if(c.maxDistance() >= distance_) { }else{
candidate = *vertex; if(neighbours==maxNeighbours)
{
candidates.push_back(*vertex);
} }
} }
}else if( con > maxCon) { }else if( con > maxCon) {
maxCon = con; maxCon = con;
maxNeighbours = noFrontNeighbours(*vertex); maxNeighbours = noFrontNeighbours(*vertex);
std::size_t distance_ = distance(*vertex, aggregates); candidates.clear();
if(c.maxDistance() >= distance_) { candidates.push_back(*vertex);
candidate = *vertex;
}
} }
}else if(oneWayCons > maxOneCons) { }else if(oneWayCons > maxOneCons) {
maxOneCons = oneWayCons; maxOneCons = oneWayCons;
maxCon = connectivity(*vertex, aggregates); maxCon = connectivity(*vertex, aggregates);
maxNeighbours = noFrontNeighbours(*vertex); maxNeighbours = noFrontNeighbours(*vertex);
std::size_t distance_ = distance(*vertex, aggregates); candidates.clear();
candidates.push_back(*vertex);
if(c.maxDistance() >= distance_) {
candidate = *vertex;
}
} }
} }
if(candidate == AggregatesMap<Vertex>::UNAGGREGATED) if(!candidates.size())
break; // No more candidates found break; // No more candidates found
distance_=distance(seed, aggregates);
aggregate_->add(candidate); candidates.resize(std::min(candidates.size(), c.maxAggregateSize()-
aggregate_->size()));
aggregate_->add(candidates);
} }
} }
...@@ -2175,7 +2318,7 @@ namespace Dune ...@@ -2175,7 +2318,7 @@ namespace Dune
graph_ = &graph; graph_ = &graph;
aggregate_ = new Aggregate<G,VertexSet>(graph, aggregates, connected_); aggregate_ = new Aggregate<G,VertexSet>(graph, aggregates, connected_, front_);
Timer watch; Timer watch;
watch.reset(); watch.reset();
...@@ -2199,8 +2342,7 @@ namespace Dune ...@@ -2199,8 +2342,7 @@ namespace Dune
// Debugging output // Debugging output
if((noAggregates+1)%10000 == 0) if((noAggregates+1)%10000 == 0)
Dune::dverb<<"c"; Dune::dverb<<"c";
unmarkFront();
aggregate_->seed(seed);
if(graph.getVertexProperties(seed).excludedBorder()) { if(graph.getVertexProperties(seed).excludedBorder()) {
aggregates[seed]=AggregatesMap<Vertex>::ISOLATED; aggregates[seed]=AggregatesMap<Vertex>::ISOLATED;
...@@ -2215,21 +2357,22 @@ namespace Dune ...@@ -2215,21 +2357,22 @@ namespace Dune
++skippedAggregates; ++skippedAggregates;
// skip rest as no agglomeration is done. // skip rest as no agglomeration is done.
continue; continue;
}else }else{
aggregate_->seed(seed);
growIsolatedAggregate(seed, aggregates, c); growIsolatedAggregate(seed, aggregates, c);
}else }
}else{
aggregate_->seed(seed);
growAggregate(seed, aggregates, c); growAggregate(seed, aggregates, c);
}
/* The rounding step. */ /* The rounding step. */
while(!(graph.getVertexProperties(seed).isolated()) && aggregate_->size() < c.maxAggregateSize()) { while(!(graph.getVertexProperties(seed).isolated()) && aggregate_->size() < c.maxAggregateSize()) {
unmarkFront(); std::vector<Vertex> candidates;
markFront(aggregates); candidates.reserve(30);
Vertex candidate = AggregatesMap<Vertex>::UNAGGREGATED;
typedef typename VertexList::const_iterator Iterator; typedef typename std::vector<Vertex>::const_iterator Iterator;
for(Iterator vertex = front_.begin(); vertex != front_.end(); ++vertex) { for(Iterator vertex = front_.begin(); vertex != front_.end(); ++vertex) {
...@@ -2252,13 +2395,15 @@ namespace Dune ...@@ -2252,13 +2395,15 @@ namespace Dune
if(distance(*vertex, aggregates) > c.maxDistance()) if(distance(*vertex, aggregates) > c.maxDistance())
continue; // Distance too far continue; // Distance too far
candidate = *vertex; candidates.push_back(*vertex);
break; break;
} }
if(candidate == AggregatesMap<Vertex>::UNAGGREGATED) break; // no more candidates found. if(!candidates.size()) break; // no more candidates found.
aggregate_->add(candidate); candidates.resize(std::min(candidates.size(), c.maxAggregateSize()-
aggregate_->size()));
aggregate_->add(candidates);
} }
...@@ -2270,7 +2415,9 @@ namespace Dune ...@@ -2270,7 +2415,9 @@ namespace Dune
if(mergedNeighbour != AggregatesMap<Vertex>::UNAGGREGATED) { if(mergedNeighbour != AggregatesMap<Vertex>::UNAGGREGATED) {
// assign vertex to the neighbouring cluster // assign vertex to the neighbouring cluster
aggregates[seed] = aggregates[mergedNeighbour]; aggregates[seed] = aggregates[mergedNeighbour];
aggregate_->invalidate();
}else{ }else{
++avg;
minA=std::min(minA,static_cast<std::size_t>(1)); minA=std::min(minA,static_cast<std::size_t>(1));
maxA=std::max(maxA,static_cast<std::size_t>(1)); maxA=std::max(maxA,static_cast<std::size_t>(1));
++oneAggregates; ++oneAggregates;
...@@ -2293,10 +2440,7 @@ namespace Dune ...@@ -2293,10 +2440,7 @@ namespace Dune
else else
++conAggregates; ++conAggregates;
} }
unmarkFront();
markFront(aggregates);
seedFromFront(stack_, graph.getVertexProperties(seed).isolated());
unmarkFront();
} }
Dune::dinfo<<"connected aggregates: "<<conAggregates; Dune::dinfo<<"connected aggregates: "<<conAggregates;
...@@ -2311,128 +2455,38 @@ namespace Dune ...@@ -2311,128 +2455,38 @@ namespace Dune
oneAggregates,skippedAggregates); oneAggregates,skippedAggregates);
} }
template<class G>
inline void Aggregator<G>::seedFromFront(Stack& stack_, bool isolated)
{
typedef typename VertexList::const_iterator Iterator;
Iterator end= front_.end();
int count=0;
for(Iterator vertex=front_.begin(); vertex != end; ++vertex,++count)
if(graph_->getVertexProperties(*vertex).isolated()==isolated)
stack_.push(*vertex);
/*
if(MINIMAL_DEBUG_LEVEL<=2 && count==0 && !isolated)
Dune::dverb<< " no vertices pushed for nonisolated aggregate!"<<std::endl;
*/
}
template<class G> template<class G>
Aggregator<G>::Stack::Stack(const MatrixGraph& graph, const Aggregator<G>& aggregatesBuilder, Aggregator<G>::Stack::Stack(const MatrixGraph& graph, const Aggregator<G>& aggregatesBuilder,
const AggregatesMap<Vertex>& aggregates) const AggregatesMap<Vertex>& aggregates)
: graph_(graph), aggregatesBuilder_(aggregatesBuilder), aggregates_(aggregates), size_(0), maxSize_(0), head_(0), filled_(0) : graph_(graph), aggregatesBuilder_(aggregatesBuilder), aggregates_(aggregates), begin_(graph.begin()), end_(graph.end())
{ {
vals_ = new Vertex[N]; //vals_ = new Vertex[N];
} }
template<class G> template<class G>
Aggregator<G>::Stack::~Stack() Aggregator<G>::Stack::~Stack()
{ {
Dune::dverb << "Max stack size was "<<maxSize_<<" filled="<<filled_<<std::endl; //Dune::dverb << "Max stack size was "<<maxSize_<<" filled="<<filled_<<std::endl;
delete[] vals_; //delete[] vals_;
} }
template<class G> template<class G>
const typename Aggregator<G>::Vertex Aggregator<G>::Stack::NullEntry const typename Aggregator<G>::Vertex Aggregator<G>::Stack::NullEntry
= std::numeric_limits<typename G::VertexDescriptor>::max(); = std::numeric_limits<typename G::VertexDescriptor>::max();
template<class G>
inline bool Aggregator<G>::Stack::push(const Vertex & v)
{
if(aggregates_[v] == AggregatesMap<Vertex>::UNAGGREGATED) {
localPush(v);
return true;
}else
return false;
}
template<class G>
inline void Aggregator<G>::Stack::localPush(const Vertex & v)
{
vals_[head_] = v;
size_ = std::min<int>(size_+1, N);
head_ = (head_+N+1)%N;
}
template<class G>
void Aggregator<G>::Stack::fill()
{
int isolated = 0, connected=0;
int isoumin, umin;
filled_++;
head_ = size_ = 0;
isoumin = umin = std::numeric_limits<int>::max();
typedef typename MatrixGraph::ConstVertexIterator Iterator;
const Iterator end = graph_.end();
for(Iterator vertex = graph_.begin(); vertex != end; ++vertex) {
// Skip already aggregated vertices
if(aggregates_[*vertex] != AggregatesMap<Vertex>::UNAGGREGATED)
continue;
if(vertex.properties().isolated()) {
isoumin = std::min(isoumin, aggregatesBuilder_.unusedNeighbours(*vertex, aggregates_));
isolated++;
}else{
umin = std::min(umin, aggregatesBuilder_.unusedNeighbours(*vertex, aggregates_));
connected++;
}
}
if(connected + isolated == 0)
// No unaggregated vertices.
return;
if(connected > 0) {
// Connected vertices have higher priority.
for(Iterator vertex = graph_.begin(); vertex != end; ++vertex)
if(aggregates_[*vertex] == AggregatesMap<Vertex>::UNAGGREGATED && !vertex.properties().isolated()
&& aggregatesBuilder_.unusedNeighbours(*vertex, aggregates_) == umin)
localPush(*vertex);
}else{
for(Iterator vertex = graph_.begin(); vertex != end; ++vertex)
if(aggregates_[*vertex] == AggregatesMap<Vertex>::UNAGGREGATED && vertex.properties().isolated()
&& aggregatesBuilder_.unusedNeighbours(*vertex, aggregates_) == isoumin)
localPush(*vertex);
}
maxSize_ = std::max(size_, maxSize_);
}
template<class G> template<class G>
inline typename G::VertexDescriptor Aggregator<G>::Stack::pop() inline typename G::VertexDescriptor Aggregator<G>::Stack::pop()
{ {
while(size_>0) { for(; begin_!=end_ && aggregates_[*begin_] != AggregatesMap<Vertex>::UNAGGREGATED; ++begin_) ;
head_ = (head_ + N -1) % N;
size_--; if(begin_!=end_)
Vertex v = vals_[head_]; {
if(aggregates_[v]==AggregatesMap<Vertex>::UNAGGREGATED) typename G::VertexDescriptor current=*begin_;
return v; ++begin_;
} return current;
// Stack is empty try to fill it }else
fill(); return NullEntry;
// try again
while(size_>0) {
head_ = (head_ + N -1) % N;
size_--;
Vertex v = vals_[head_];
if(aggregates_[v]==AggregatesMap<Vertex>::UNAGGREGATED)
return v;
}
return NullEntry;
} }
#endif // DOXYGEN #endif // DOXYGEN
......
dune/istl/paamg/fastamg.hh 0 → 100644
+ 627
− 0
View file @ 9348a968
This diff is collapsed.
dune/istl/paamg/fastamgsmoother.hh 0 → 100644
+ 96
− 0
View file @ 9348a968
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_ISTL_FASTAMGSMOOTHER_HH
#define DUNE_ISTL_FASTAMGSMOOTHER_HH
namespace Dune
{
namespace Amg
{
template<std::size_t level>
struct GaussSeidelPresmoothDefect {
template<typename M, typename X, typename Y>
static void apply(const M& A, X& x, Y& d,
const Y& b)
{
typedef typename M::ConstRowIterator RowIterator;
typedef typename M::ConstColIterator ColIterator;
typedef typename Y::block_type YBlock;
typedef typename X::block_type XBlock;
typename Y::iterator dIter=d.begin();
typename Y::const_iterator bIter=b.begin();
typename X::iterator xIter=x.begin();
for(RowIterator row=A.begin(), end=A.end(); row != end;
++row, ++dIter, ++xIter, ++bIter)
{
ColIterator endCol=(*row).end();
ColIterator col=(*row).begin();
*dIter = *bIter;
for (; col.index()<row.index(); ++col)
(*col).mmv(x[col.index()],*dIter); // rhs -= sum_{j<i} a_ij * xnew_j
assert(row.index()==col.index());
ColIterator diag=col; // upper diagonal matrix not needed as x was 0 before.
// Not recursive yet. Just solve with the diagonal
diag->solve(*xIter,*dIter);
*dIter=0; //as r=v
// Update residual for the symmetric case
for(col=(*row).begin(); col.index()<row.index(); ++col)
col->mmv(*xIter, d[col.index()]); //d_j-=A_ij x_i
}
}
};
template<std::size_t level>
struct GaussSeidelPostsmoothDefect {
template<typename M, typename X, typename Y>
static void apply(const M& A, X& x, Y& d,
const Y& b)
{
typedef typename M::ConstRowIterator RowIterator;
typedef typename M::ConstColIterator ColIterator;
typedef typename Y::block_type YBlock;
typedef typename X::block_type XBlock;
typename Y::iterator dIter=d.beforeEnd();
typename X::iterator xIter=x.beforeEnd();
typename Y::const_iterator bIter=b.beforeEnd();
for(RowIterator row=A.beforeEnd(), end=A.beforeBegin(); row != end;
--row, --dIter, --xIter, --bIter)
{
ColIterator endCol=(*row).beforeBegin();
ColIterator col=(*row).beforeEnd();
*dIter = *bIter;
for (; col.index()>row.index(); --col)
(*col).mmv(x[col.index()],*dIter); // rhs -= sum_{i>j} a_ij * xnew_j
assert(row.index()==col.index());
ColIterator diag=col;
YBlock v=*dIter;
// upper diagonal matrix
for (--col; col!=endCol; --col)
(*col).mmv(x[col.index()],v); // v -= sum_{j<i} a_ij * xold_j
// Not recursive yet. Just solve with the diagonal
diag->solve(*xIter,v);
*dIter-=v;
// Update residual for the symmetric case
// Skip residual computation as it is not needed.
//for(col=(*row).begin();col.index()<row.index(); ++col)
//col.mmv(*xIter, d[col.index()]); //d_j-=A_ij x_i
}
}
};
} // end namespace Amg
} // end namespace Dune
#endif
dune/istl/paamg/hierarchy.hh
+ 5
− 4
View file @ 9348a968
...@@ -869,8 +869,8 @@ namespace Dune ...@@ -869,8 +869,8 @@ namespace Dune
*(get<1>(graphs)), *(get<1>(graphs)),
visitedMap, visitedMap,
*aggregatesMap, *aggregatesMap,
*infoLevel); *infoLevel,
noAggregates);
GraphCreator::free(graphs); GraphCreator::free(graphs);
if(criterion.debugLevel()>2) { if(criterion.debugLevel()>2) {
...@@ -909,7 +909,8 @@ namespace Dune ...@@ -909,7 +909,8 @@ namespace Dune
*aggregatesMap, *aggregatesMap,
aggregates, aggregates,
OverlapFlags()); OverlapFlags());
dverb<<"Building of sparsity pattern took "<<watch.elapsed()<<std::endl;
watch.reset();
info->freeGlobalLookup(); info->freeGlobalLookup();
delete get<0>(graphs); delete get<0>(graphs);
...@@ -917,7 +918,7 @@ namespace Dune ...@@ -917,7 +918,7 @@ namespace Dune
if(criterion.debugLevel()>2) { if(criterion.debugLevel()>2) {
if(rank==0) if(rank==0)
std::cout<<"Calculation of Galerkin product took "<<watch.elapsed()<<" seconds."<<std::endl; std::cout<<"Calculation entries of Galerkin product took "<<watch.elapsed()<<" seconds."<<std::endl;
} }
BIGINT nonzeros = countNonZeros(*coarseMatrix); BIGINT nonzeros = countNonZeros(*coarseMatrix);
......
dune/istl/paamg/indicescoarsener.hh
+ 10
− 25
View file @ 9348a968
...@@ -90,7 +90,8 @@ namespace Dune ...@@ -90,7 +90,8 @@ namespace Dune
Graph& fineGraph, Graph& fineGraph,
VM& visitedMap, VM& visitedMap,
AggregatesMap<typename Graph::VertexDescriptor>& aggregates, AggregatesMap<typename Graph::VertexDescriptor>& aggregates,
ParallelInformation& coarseInfo); ParallelInformation& coarseInfo,
typename Graph::VertexDescriptor noAggregates);
private: private:
template<typename G, typename I> template<typename G, typename I>
...@@ -212,11 +213,12 @@ namespace Dune ...@@ -212,11 +213,12 @@ namespace Dune
public: public:
template<typename Graph, typename VM> template<typename Graph, typename VM>
static typename Graph::VertexDescriptor static typename Graph::VertexDescriptor
coarsen(const SequentialInformation& fineInfo, coarsen(const SequentialInformation & fineInfo,
Graph& fineGraph, Graph& fineGraph,
VM& visitedMap, VM& visitedMap,
AggregatesMap<typename Graph::VertexDescriptor>& aggregates, AggregatesMap<typename Graph::VertexDescriptor>& aggregates,
SequentialInformation& coarseInfo); SequentialInformation& coarseInfo,
typename Graph::VertexDescriptor noAggregates);
}; };
#if HAVE_MPI #if HAVE_MPI
...@@ -227,7 +229,8 @@ namespace Dune ...@@ -227,7 +229,8 @@ namespace Dune
Graph& fineGraph, Graph& fineGraph,
VM& visitedMap, VM& visitedMap,
AggregatesMap<typename Graph::VertexDescriptor>& aggregates, AggregatesMap<typename Graph::VertexDescriptor>& aggregates,
ParallelInformation& coarseInfo) ParallelInformation& coarseInfo,
typename Graph::VertexDescriptor noAggregates)
{ {
ParallelAggregateRenumberer<Graph,typename ParallelInformation::GlobalLookupIndexSet> renumberer(aggregates, fineInfo.globalLookup()); ParallelAggregateRenumberer<Graph,typename ParallelInformation::GlobalLookupIndexSet> renumberer(aggregates, fineInfo.globalLookup());
buildCoarseIndexSet(fineInfo, fineGraph, visitedMap, aggregates, buildCoarseIndexSet(fineInfo, fineGraph, visitedMap, aggregates,
...@@ -387,28 +390,10 @@ namespace Dune ...@@ -387,28 +390,10 @@ namespace Dune
Graph& fineGraph, Graph& fineGraph,
VM& visitedMap, VM& visitedMap,
AggregatesMap<typename Graph::VertexDescriptor>& aggregates, AggregatesMap<typename Graph::VertexDescriptor>& aggregates,
SequentialInformation& coarseInfo) SequentialInformation& coarseInfo,
typename Graph::VertexDescriptor noAggregates)
{ {
typedef typename Graph::VertexDescriptor Vertex; return noAggregates;
AggregateRenumberer<Graph> renumberer(aggregates);
typedef typename Graph::VertexIterator Iterator;
for(Iterator vertex=fineGraph.begin(), endVertex=fineGraph.end();
vertex != endVertex; ++vertex)
if(aggregates[*vertex]!=AggregatesMap<Vertex>::ISOLATED &&
!get(visitedMap, *vertex)) {
aggregates.template breadthFirstSearch<false>(*vertex, aggregates[*vertex],
fineGraph, renumberer, visitedMap);
aggregates[*vertex] = renumberer;
++renumberer;
}
for(Iterator vertex=fineGraph.begin(), endVertex=fineGraph.end();
vertex != endVertex; ++vertex)
put(visitedMap, *vertex, false);
return renumberer;
} }
} //namespace Amg } //namespace Amg
......
dune/istl/paamg/test/Makefile.am
+ 9
− 1
View file @ 9348a968
...@@ -5,7 +5,7 @@ if MPI ...@@ -5,7 +5,7 @@ if MPI
TESTPROGS = galerkintest hierarchytest pamgtest transfertest pamg_comm_repart_test TESTPROGS = galerkintest hierarchytest pamgtest transfertest pamg_comm_repart_test
endif endif
NORMALTESTS = kamgtest amgtest graphtest $(MPITESTS) NORMALTESTS = kamgtest amgtest fastamg graphtest $(MPITESTS)
# which tests to run # which tests to run
TESTS = $(NORMALTESTS) $(TESTPROGS) TESTS = $(NORMALTESTS) $(TESTPROGS)
...@@ -51,6 +51,14 @@ amgtest_LDADD = \ ...@@ -51,6 +51,14 @@ amgtest_LDADD = \
$(SUPERLU_LIBS) \ $(SUPERLU_LIBS) \
$(LDADD) $(LDADD)
fastamg_SOURCES = fastamg.cc
fastamg_CPPFLAGS = $(AM_CPPFLAGS) $(SUPERLU_CPPFLAGS)
fastamg_LDFLAGS = $(AM_LDFLAGS) $(SUPERLU_LDFLAGS)
fastamg_LDADD = \
$(SUPERLU_LIBS) \
$(LDADD)
kamgtest_SOURCES = kamgtest.cc kamgtest_SOURCES = kamgtest.cc
kamgtest_CPPFLAGS = $(AM_CPPFLAGS) $(SUPERLU_CPPFLAGS) kamgtest_CPPFLAGS = $(AM_CPPFLAGS) $(SUPERLU_CPPFLAGS)
kamgtest_LDFLAGS = $(AM_LDFLAGS) $(SUPERLU_LDFLAGS) kamgtest_LDFLAGS = $(AM_LDFLAGS) $(SUPERLU_LDFLAGS)
......
dune/istl/paamg/test/amgtest.cc
+ 6
− 4
View file @ 9348a968
...@@ -94,10 +94,12 @@ void testAMG(int N, int coarsenTarget, int ml) ...@@ -94,10 +94,12 @@ void testAMG(int N, int coarsenTarget, int ml)
watch.reset(); watch.reset();
Operator fop(mat); Operator fop(mat);
typedef Dune::Amg::Dependency<BCRSMat,Dune::Amg::FirstDiagonal> //typedef Dune::Amg::CoarsenCriterion<Dune::Amg::UnSymmetricCriterion<BCRSMat,Dune::Amg::FirstDiagonal> >
Criterion; // Criterion;
typedef Dune::SeqSSOR<BCRSMat,Vector,Vector> Smoother; typedef Dune::Amg::AggregationCriterion<Dune::Amg::SymmetricMatrixDependency<BCRSMat,Dune::Amg::FirstDiagonal> > CriterionBase;
//typedef Dune::SeqSOR<BCRSMat,Vector,Vector> Smoother; typedef Dune::Amg::CoarsenCriterion<CriterionBase> Criterion;
//typedef Dune::SeqSSOR<BCRSMat,Vector,Vector> Smoother;
typedef Dune::SeqSOR<BCRSMat,Vector,Vector> Smoother;
//typedef Dune::SeqJac<BCRSMat,Vector,Vector> Smoother; //typedef Dune::SeqJac<BCRSMat,Vector,Vector> Smoother;
//typedef Dune::SeqOverlappingSchwarz<BCRSMat,Vector,Dune::MultiplicativeSchwarzMode> Smoother; //typedef Dune::SeqOverlappingSchwarz<BCRSMat,Vector,Dune::MultiplicativeSchwarzMode> Smoother;
//typedef Dune::SeqOverlappingSchwarz<BCRSMat,Vector,Dune::SymmetricMultiplicativeSchwarzMode> Smoother; //typedef Dune::SeqOverlappingSchwarz<BCRSMat,Vector,Dune::SymmetricMultiplicativeSchwarzMode> Smoother;
......
dune/istl/paamg/test/fastamg.cc 0 → 100644
+ 131
− 0
View file @ 9348a968
// -*- 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/common/timer.hh>
#include <dune/common/parallel/indexset.hh>
#include <dune/common/parallel/collectivecommunication.hh>
#include <dune/istl/paamg/fastamg.hh>
#include <dune/istl/paamg/pinfo.hh>
#include <dune/istl/solvers.hh>
#include <cstdlib>
#include <ctime>
template<class M, class V>
void randomize(const M& mat, V& b)
{
V x=b;
srand((unsigned)std::clock());
typedef typename V::iterator iterator;
for(iterator i=x.begin(); i != x.end(); ++i)
*i=(rand() / (RAND_MAX + 1.0));
mat.mv(static_cast<const V&>(x), b);
}
template <int BS>
void testAMG(int N, int coarsenTarget, int ml)
{
std::cout<<"N="<<N<<" coarsenTarget="<<coarsenTarget<<" maxlevel="<<ml<<std::endl;
typedef Dune::ParallelIndexSet<int,LocalIndex,512> ParallelIndexSet;
ParallelIndexSet indices;
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::CollectiveCommunication<void*> Comm;
int n;
Comm c;
BCRSMat mat = setupAnisotropic2d<BS,double>(N, indices, c, &n, 1);
Vector b(mat.N()), x(mat.M());
b=0;
x=100;
setBoundary(x, b, N);
x=0;
randomize(mat, b);
if(N<6) {
Dune::printmatrix(std::cout, mat, "A", "row");
Dune::printvector(std::cout, x, "x", "row");
}
Dune::Timer watch;
watch.reset();
Operator fop(mat);
typedef Dune::Amg::AggregationCriterion<Dune::Amg::SymmetricMatrixDependency<BCRSMat,Dune::Amg::FirstDiagonal> > CriterionBase;
typedef Dune::Amg::CoarsenCriterion<CriterionBase> Criterion;
Criterion criterion(15,coarsenTarget);
criterion.setDefaultValuesIsotropic(2);
criterion.setAlpha(.67);
criterion.setBeta(1.0e-4);
criterion.setMaxLevel(ml);
criterion.setSkipIsolated(false);
Dune::SeqScalarProduct<Vector> sp;
typedef Dune::Amg::FastAMG<Operator,Vector> AMG;
Dune::Amg::Parameters parms;
AMG amg(fop, criterion, parms);
double buildtime = watch.elapsed();
std::cout<<"Building hierarchy took "<<buildtime<<" seconds"<<std::endl;
Dune::GeneralizedPCGSolver<Vector> amgCG(fop,amg,1e-6,80,2);
//Dune::LoopSolver<Vector> amgCG(fop, amg, 1e-4, 10000, 2);
watch.reset();
Dune::InverseOperatorResult r;
amgCG.apply(x,b,r);
double solvetime = watch.elapsed();
std::cout<<"AMG solving took "<<solvetime<<" seconds"<<std::endl;
std::cout<<"AMG building took "<<(buildtime/r.elapsed*r.iterations)<<" iterations"<<std::endl;
std::cout<<"AMG building together with solving took "<<buildtime+solvetime<<std::endl;
/*
watch.reset();
cg.apply(x,b,r);
std::cout<<"CG solving took "<<watch.elapsed()<<" seconds"<<std::endl;
*/
}
int main(int argc, char** argv)
{
int N=100;
int coarsenTarget=1200;
int ml=10;
if(argc>1)
N = atoi(argv[1]);
if(argc>2)
coarsenTarget = atoi(argv[2]);
if(argc>3)
ml = atoi(argv[3]);
testAMG<1>(N, coarsenTarget, ml);
testAMG<2>(N, coarsenTarget, ml);
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment