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

Make transfer work for std::vector. 

Project coarsest aggregates to finest level for visualization.

[[Imported from SVN: r1035]]
parent 93d268b0
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
istl/paamg/hierarchy.hh
+ 74
− 3
View file @ f7ebac08
......@@ -25,6 +25,7 @@
#include <dune/istl/paamg/globalaggregates.hh>
#include <dune/istl/paamg/construction.hh>
#include <dune/istl/paamg/smoother.hh>
#include <dune/istl/paamg/transfer.hh>
namespace Dune
{
......@@ -124,7 +125,7 @@ namespace Dune
* @brief Iterator over the levels in the hierarchy.
*
* operator++() moves to the next coarser level in the hierarchy.
* while operator--() moved to the next finer level in the hierarchy.
* while operator--() moves to the next finer level in the hierarchy.
*/
template<class C, class T1>
class LevelIterator
......@@ -386,6 +387,18 @@ namespace Dune
return prolongDamp_;
}
/**
* @brief Get the mapping of fine level unknowns to coarse level
* aggregates.
*
* For each fine level unknown i the correcponding data[i] is the
* aggregate it belongs to on the coarsest level.
*
* @param[out] data The mapping of fine level unknowns to coarse level
* aggregates.
*/
void getCoarsestAggregatesOnFinest(std::vector<std::size_t>& data);
private:
typedef typename ConstructionTraits<MatrixOperator>::Arguments MatrixArgs;
/** @brief The list of aggregates maps. */
......@@ -815,6 +828,65 @@ namespace Dune
return parallelInformation_;
}
template<class M, class IS, class A>
void MatrixHierarchy<M,IS,A>::getCoarsestAggregatesOnFinest(std::vector<std::size_t>& data)
{
int levels=aggregatesMaps().size();
int maxlevels=parallelInformation_.finest()->communicator().max(levels);
std::size_t size=(*(aggregatesMaps().begin()))->noVertices();
// We need an auxiliary vector for the consecutive prolongation.
std::vector<std::size_t> tmp;
std::vector<std::size_t> *coarse, *fine;
// make sure the allocated space suffices.
tmp.reserve(size);
data.reserve(size);
// Correctly assign coarse and fine for the first prolongation such that
// we end up in data in the end.
if(levels%2==0) {
coarse=&tmp;
fine=&data;
}else{
coarse=&data;
fine=&tmp;
}
// Number the unknowns on the coarsest level consecutively for each process.
if(levels==maxlevels) {
const AggregatesMap& map = *(*(++aggregatesMaps().rbegin()));
std::size_t m=0;
for(typename AggregatesMap::const_iterator iter = map.begin(); iter != map.end(); ++iter)
m=std::max(*iter,m);
coarse->resize(m+1);
std::size_t i=0;
for(typename std::vector<std::size_t>::iterator iter=coarse->begin(); iter != coarse->end();
++iter, ++i)
*iter=i;
}
typename ParallelInformationHierarchy::Iterator pinfo = parallelInformation().coarsest();
--pinfo;
int l=levels;
// Now consecutively project the numbers to the finest level.
for(typename AggregatesMapList::const_reverse_iterator aggregates=++aggregatesMaps().rbegin();
aggregates != aggregatesMaps().rend(); ++aggregates,--levels) {
fine->resize((*aggregates)->noVertices());
fine->assign(fine->size(), 0);
Transfer<typename AggregatesMap::AggregateDescriptor, std::vector<std::size_t>, ParallelInformation>
::prolongate(*(*aggregates), *coarse, *fine, static_cast<std::size_t>(1), *pinfo);
--pinfo;
std::swap(coarse, fine);
}
// Assertion to check that we really projected to data on the last step.
assert(coarse==&data);
}
template<class M, class IS, class A>
const typename MatrixHierarchy<M,IS,A>::AggregatesMapList&
MatrixHierarchy<M,IS,A>::aggregatesMaps() const
......@@ -958,14 +1030,13 @@ namespace Dune
coarsest_->element_ = ConstructionTraits<MemberType>::construct(args);
finest_ = coarsest_;
coarsest_->finer_ = 0;
coarsest_->redistributed_ = 0;
}else{
coarsest_->coarser_ = allocator_.allocate(1,0);
coarsest_->coarser_->finer_ = coarsest_;
coarsest_ = coarsest_->coarser_;
coarsest_->element_ = ConstructionTraits<MemberType>::construct(args);
finest_->redistributed_ = 0;
}
coarsest_->redistributed_ = 0;
coarsest_->coarser_=0;
++levels_;
}
......
istl/paamg/test/hierarchytest.cc
+ 4
− 1
View file @ f7ebac08
......@@ -15,7 +15,7 @@ int main(int argc, char** argv)
{
MPI_Init(&argc, &argv);
const int BS=1, N=100;
const int BS=1, N=10;
int procs, rank;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
......@@ -71,6 +71,9 @@ int main(int argc, char** argv)
hierarchy.recalculateGalerkin(OverlapFlags());
std::vector<std::size_t> data;
hierarchy.getCoarsestAggregatesOnFinest(data);
MPI_Finalize();
......
istl/paamg/transfer.hh
+ 65
− 52
View file @ f7ebac08
......@@ -32,21 +32,23 @@ namespace Dune
typedef V1 Vertex;
typedef V2 Vector;
template<typename T1>
static void prolongate(const AggregatesMap<Vertex>& aggregates, Vector& coarse, Vector& fine,
typename Vector::field_type damp);
T1 damp);
static void restrict (const AggregatesMap<Vertex>& aggregates, Vector& coarse, const Vector & fine,
T& comm);
};
template<class V,class B>
class Transfer<V,BlockVector<B>, SequentialInformation>
template<class V,class V1>
class Transfer<V,V1, SequentialInformation>
{
public:
typedef V Vertex;
typedef BlockVector<B> Vector;
typedef V1 Vector;
template<typename T1>
static void prolongate(const AggregatesMap<Vertex>& aggregates, Vector& coarse, Vector& fine,
typename Vector::field_type damp,
T1 damp,
const SequentialInformation& comm=SequentialInformation());
......@@ -56,27 +58,29 @@ namespace Dune
#if HAVE_MPI
template<class V,class B, class T>
class Transfer<V,BlockVector<B>,ParallelInformation<T> >
template<class V,class V1, class T>
class Transfer<V,V1,ParallelInformation<T> >
{
public:
typedef V Vertex;
typedef BlockVector<B> Vector;
typedef V1 Vector;
template<typename T1>
static void prolongate(const AggregatesMap<Vertex>& aggregates, Vector& coarse, Vector& fine,
typename Vector::field_type damp);
T1 damp);
static void restrict (const AggregatesMap<Vertex>& aggregates, Vector& coarse, const Vector & fine,
ParallelInformation<T>& comm);
};
template<class V,class B, class T1, class T2>
class Transfer<V,BlockVector<B>,OwnerOverlapCopyCommunication<T1,T2> >
template<class V,class V1, class T1, class T2>
class Transfer<V,V1,OwnerOverlapCopyCommunication<T1,T2> >
{
public:
typedef V Vertex;
typedef BlockVector<B> Vector;
typedef V1 Vector;
template<typename T3>
static void prolongate(const AggregatesMap<Vertex>& aggregates, Vector& coarse, Vector& fine,
typename Vector::field_type damp, OwnerOverlapCopyCommunication<T1,T2>& comm);
T3 damp, OwnerOverlapCopyCommunication<T1,T2>& comm);
static void restrict (const AggregatesMap<Vertex>& aggregates, Vector& coarse, const Vector & fine,
OwnerOverlapCopyCommunication<T1,T2>& comm);
......@@ -84,78 +88,87 @@ namespace Dune
#endif
template<class V, class B>
inline void Transfer<V,BlockVector<B>,SequentialInformation>::prolongate(const AggregatesMap<Vertex>& aggregates,
Vector& coarse, Vector& fine,
typename Vector::field_type damp,
const SequentialInformation& comm)
template<class V, class V1>
template<typename T>
inline void Transfer<V,V1,SequentialInformation>::prolongate(const AggregatesMap<Vertex>& aggregates,
Vector& coarse, Vector& fine,
T damp,
const SequentialInformation& comm)
{
typedef typename Vector::iterator Iterator;
Iterator end = fine.end();
coarse *= damp;
Iterator end = coarse.end();
Iterator begin= coarse.begin();
for(; begin!=end; ++begin)
*begin*=damp;
end=fine.end();
begin=fine.begin();
for(Iterator block=fine.begin(); block != end; ++block) {
const Vertex& vertex = aggregates[block.index()];
for(Iterator block=begin; block != end; ++block) {
std::ptrdiff_t index=block-begin;
const Vertex& vertex = aggregates[index];
if(vertex != AggregatesMap<Vertex>::ISOLATED)
*block += coarse[aggregates[block.index()]];
*block += coarse[aggregates[index]];
}
}
template<class V, class B>
inline void Transfer<V,BlockVector<B>,SequentialInformation>::restrict (const AggregatesMap<Vertex>& aggregates,
Vector& coarse,
const Vector & fine,
const SequentialInformation & comm)
template<class V, class V1>
inline void Transfer<V,V1,SequentialInformation>::restrict (const AggregatesMap<Vertex>& aggregates,
Vector& coarse,
const Vector & fine,
const SequentialInformation & comm)
{
// Set coarse vector to zero
coarse=0;
typedef typename Vector::const_iterator Iterator;
Iterator end = fine.end();
for(Iterator block=fine.begin(); block != end; ++block) {
const Vertex& vertex = aggregates[block.index()];
Iterator begin=fine.begin();
for(Iterator block=begin; block != end; ++block) {
const Vertex& vertex = aggregates[block-begin];
if(vertex != AggregatesMap<Vertex>::ISOLATED)
coarse[vertex] += *block;
}
}
#if HAVE_MPI
template<class V, class B, class T>
inline void Transfer<V,BlockVector<B>,ParallelInformation<T> >::prolongate(const AggregatesMap<Vertex>& aggregates,
Vector& coarse, Vector& fine,
typename Vector::field_type damp)
template<class V, class V1, class T>
template<typename T1>
inline void Transfer<V,V1,ParallelInformation<T> >::prolongate(const AggregatesMap<Vertex>& aggregates,
Vector& coarse, Vector& fine,
T1 damp)
{
Transfer<V,BlockVector<B>,SequentialInformation>::prolongate(aggregates, coarse, fine, damp);
Transfer<V,V1,SequentialInformation>::prolongate(aggregates, coarse, fine, damp);
}
template<class V, class B, class T>
inline void Transfer<V,BlockVector<B>,ParallelInformation<T> >::restrict (const AggregatesMap<Vertex>& aggregates,
Vector& coarse, const Vector & fine,
ParallelInformation<T>& comm)
template<class V, class V1, class T>
inline void Transfer<V,V1,ParallelInformation<T> >::restrict (const AggregatesMap<Vertex>& aggregates,
Vector& coarse, const Vector & fine,
ParallelInformation<T>& comm)
{
Transfer<V,BlockVector<B>,SequentialInformation>::restrict (aggregates, coarse, fine, SequentialInformation());
Transfer<V,V1,SequentialInformation>::restrict (aggregates, coarse, fine, SequentialInformation());
// We need this here to avoid it in the smoothers on the coarse level.
// There (in the preconditioner d is const.
comm.project(coarse);
}
template<class V, class B, class T1, class T2>
inline void Transfer<V,BlockVector<B>,OwnerOverlapCopyCommunication<T1,T2> >::prolongate(const AggregatesMap<Vertex>& aggregates,
Vector& coarse, Vector& fine,
typename Vector::field_type damp,
OwnerOverlapCopyCommunication<T1,T2>& comm)
template<class V, class V1, class T1, class T2>
template<typename T3>
inline void Transfer<V,V1,OwnerOverlapCopyCommunication<T1,T2> >::prolongate(const AggregatesMap<Vertex>& aggregates,
Vector& coarse, Vector& fine,
T3 damp,
OwnerOverlapCopyCommunication<T1,T2>& comm)
{
Transfer<V,BlockVector<B>,SequentialInformation>::prolongate(aggregates, coarse, fine, damp);
Transfer<V,V1,SequentialInformation>::prolongate(aggregates, coarse, fine, damp);
}
template<class V, class B, class T1, class T2>
inline void Transfer<V,BlockVector<B>,OwnerOverlapCopyCommunication<T1,T2> >::restrict (const AggregatesMap<Vertex>& aggregates,
Vector& coarse, const Vector & fine,
OwnerOverlapCopyCommunication<T1,T2>& comm)
template<class V, class V1, class T1, class T2>
inline void Transfer<V,V1,OwnerOverlapCopyCommunication<T1,T2> >::restrict (const AggregatesMap<Vertex>& aggregates,
Vector& coarse, const Vector & fine,
OwnerOverlapCopyCommunication<T1,T2>& comm)
{
Transfer<V,BlockVector<B>,SequentialInformation>::restrict (aggregates, coarse, fine, SequentialInformation());
Transfer<V,V1,SequentialInformation>::restrict (aggregates, coarse, fine, SequentialInformation());
// We need this here to avoid it in the smoothers on the coarse level.
// There (in the preconditioner d is const.
comm.project(coarse);
......
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