"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "dune/istl/novlpschwarz.hh" between
dune-istl-2.6.0.tar.gz and dune-istl-2.7.0.tar.gz

About: dune-istl - DUNE (Distributed and Unified Numerics Environment) is a modular toolbox for solving partial differential equations (PDEs) with grid-based methods: DUNE Iterative Solver Template Library.

novlpschwarz.hh  (dune-istl-2.6.0):novlpschwarz.hh  (dune-istl-2.7.0)
skipping to change at line 92 skipping to change at line 92
typedef typename RIMap::iterator RIMapit; typedef typename RIMap::iterator RIMapit;
/** /**
* @brief constructor: just store a reference to a matrix. * @brief constructor: just store a reference to a matrix.
* *
* @param A The assembled matrix. * @param A The assembled matrix.
* @param com The communication object for syncing owner and copy * @param com The communication object for syncing owner and copy
* data points. (E.~g. OwnerOverlapCommunication ) * data points. (E.~g. OwnerOverlapCommunication )
*/ */
NonoverlappingSchwarzOperator (const matrix_type& A, const communication_typ e& com) NonoverlappingSchwarzOperator (const matrix_type& A, const communication_typ e& com)
: _A_(stackobject_to_shared_ptr(A)), communication(com), buildcomm(true)
{}
NonoverlappingSchwarzOperator (std::shared_ptr<const matrix_type> A, const c
ommunication_type& com)
: _A_(A), communication(com), buildcomm(true) : _A_(A), communication(com), buildcomm(true)
{} {}
//! apply operator to x: \f$ y = A(x) \f$ //! apply operator to x: \f$ y = A(x) \f$
virtual void apply (const X& x, Y& y) const virtual void apply (const X& x, Y& y) const
{ {
y = 0; y = 0;
novlp_op_apply(x,y,1); novlp_op_apply(x,y,1);
communication.addOwnerCopyToOwnerCopy(y,y); communication.addOwnerCopyToOwnerCopy(y,y);
} }
skipping to change at line 118 skipping to change at line 122
Y y1(y); Y y1(y);
y = 0; y = 0;
novlp_op_apply(x,y,alpha); novlp_op_apply(x,y,alpha);
communication.addOwnerCopyToOwnerCopy(y,y); communication.addOwnerCopyToOwnerCopy(y,y);
y += y1; y += y1;
} }
//! get matrix via * //! get matrix via *
virtual const matrix_type& getmat () const virtual const matrix_type& getmat () const
{ {
return _A_; return *_A_;
} }
void novlp_op_apply (const X& x, Y& y, field_type alpha) const void novlp_op_apply (const X& x, Y& y, field_type alpha) const
{ {
//get index sets //get index sets
const PIS& pis=communication.indexSet(); const PIS& pis=communication.indexSet();
const RI& ri = communication.remoteIndices(); const RI& ri = communication.remoteIndices();
// at the beginning make a multimap "bordercontribution". // at the beginning make a multimap "bordercontribution".
// process has i and j as border dofs but is not the owner // process has i and j as border dofs but is not the owner
skipping to change at line 152 skipping to change at line 156
} }
for (MM::iterator iter = bordercontribution.begin(); for (MM::iterator iter = bordercontribution.begin();
iter != bordercontribution.end(); ++iter) iter != bordercontribution.end(); ++iter)
bordercontribution.erase(iter); bordercontribution.erase(iter);
std::map<int,int> owner; //key: local index i, value: process, that owns i std::map<int,int> owner; //key: local index i, value: process, that owns i
RIMap rimap; RIMap rimap;
// for each local index make multimap rimap: // for each local index make multimap rimap:
// key: local index i, data: pair of process that knows i and pointer to RI entry // key: local index i, data: pair of process that knows i and pointer to RI entry
for (RowIterator i = _A_.begin(); i != _A_.end(); ++i) for (RowIterator i = _A_->begin(); i != _A_->end(); ++i)
if (mask[i.index()] == 0) if (mask[i.index()] == 0)
for (RIIterator remote = ri.begin(); remote != ri.end(); ++remote) { for (RIIterator remote = ri.begin(); remote != ri.end(); ++remote) {
RIL& ril = *(remote->second.first); RIL& ril = *(remote->second.first);
for (RILIterator rindex = ril.begin(); rindex != ril.end(); ++rind ex) for (RILIterator rindex = ril.begin(); rindex != ril.end(); ++rind ex)
if (rindex->attribute() != OwnerOverlapCopyAttributeSet::overlap ) if (rindex->attribute() != OwnerOverlapCopyAttributeSet::overlap )
if (rindex->localIndexPair().local().local() == i.index()) { if (rindex->localIndexPair().local().local() == i.index()) {
rimap.insert rimap.insert
(std::make_pair(i.index(), (std::make_pair(i.index(),
std::pair<int,RILIterator>(remote->first, rindex))); std::pair<int,RILIterator>(remote->first, rindex)));
if(rindex->attribute()==OwnerOverlapCopyAttributeSet::owner) if(rindex->attribute()==OwnerOverlapCopyAttributeSet::owner)
owner.insert(std::make_pair(i.index(),remote->first)); owner.insert(std::make_pair(i.index(),remote->first));
} }
} }
int iowner = 0; int iowner = 0;
for (RowIterator i = _A_.begin(); i != _A_.end(); ++i) { for (RowIterator i = _A_->begin(); i != _A_->end(); ++i) {
if (mask[i.index()] == 0) { if (mask[i.index()] == 0) {
std::map<int,int>::iterator it = owner.find(i.index()); std::map<int,int>::iterator it = owner.find(i.index());
iowner = it->second; iowner = it->second;
std::pair<RIMapit, RIMapit> foundiit = rimap.equal_range(i.index()); std::pair<RIMapit, RIMapit> foundiit = rimap.equal_range(i.index());
for (ColIterator j = _A_[i.index()].begin(); j != _A_[i.index()].end (); ++j) { for (ColIterator j = (*_A_)[i.index()].begin(); j != (*_A_)[i.index( )].end(); ++j) {
if (mask[j.index()] == 0) { if (mask[j.index()] == 0) {
bool flag = true; bool flag = true;
for (RIMapit foundi = foundiit.first; foundi != foundiit.second; ++foundi) { for (RIMapit foundi = foundiit.first; foundi != foundiit.second; ++foundi) {
std::pair<RIMapit, RIMapit> foundjit = rimap.equal_range(j.ind ex()); std::pair<RIMapit, RIMapit> foundjit = rimap.equal_range(j.ind ex());
for (RIMapit foundj = foundjit.first; foundj != foundjit.secon d; ++foundj) for (RIMapit foundj = foundjit.first; foundj != foundjit.secon d; ++foundj)
if (foundj->second.first == foundi->second.first) if (foundj->second.first == foundi->second.first)
if (foundj->second.second->attribute() == OwnerOverlapCopy AttributeSet::owner if (foundj->second.second->attribute() == OwnerOverlapCopy AttributeSet::owner
|| foundj->second.first == iowner || foundj->second.first == iowner
|| foundj->second.first < communication.communicator( ).rank()) { || foundj->second.first < communication.communicator( ).rank()) {
flag = false; flag = false;
skipping to change at line 206 skipping to change at line 210
if (flag==true) if (flag==true)
bordercontribution.insert(std::pair<int,int>(i.index(),j.index ())); bordercontribution.insert(std::pair<int,int>(i.index(),j.index ()));
} }
} }
} }
} }
buildcomm = false; buildcomm = false;
} }
//compute alpha*A*x nonoverlapping case //compute alpha*A*x nonoverlapping case
for (RowIterator i = _A_.begin(); i != _A_.end(); ++i) { for (RowIterator i = _A_->begin(); i != _A_->end(); ++i) {
if (mask[i.index()] == 0) { if (mask[i.index()] == 0) {
//dof doesn't belong to process but is border (not ghost) //dof doesn't belong to process but is border (not ghost)
for (ColIterator j = _A_[i.index()].begin(); j != _A_[i.index()].end() ; ++j) { for (ColIterator j = (*_A_)[i.index()].begin(); j != (*_A_)[i.index()] .end(); ++j) {
if (mask[j.index()] == 1) //j is owner => then sum entries if (mask[j.index()] == 1) //j is owner => then sum entries
(*j).usmv(alpha,x[j.index()],y[i.index()]); (*j).usmv(alpha,x[j.index()],y[i.index()]);
else if (mask[j.index()] == 0) { else if (mask[j.index()] == 0) {
std::pair<MM::iterator, MM::iterator> itp = std::pair<MM::iterator, MM::iterator> itp =
bordercontribution.equal_range(i.index()); bordercontribution.equal_range(i.index());
for (MM::iterator it = itp.first; it != itp.second; ++it) for (MM::iterator it = itp.first; it != itp.second; ++it)
if ((*it).second == (int)j.index()) if ((*it).second == (int)j.index())
(*j).usmv(alpha,x[j.index()],y[i.index()]); (*j).usmv(alpha,x[j.index()],y[i.index()]);
} }
} }
} }
else if (mask[i.index()] == 1) { else if (mask[i.index()] == 1) {
for (ColIterator j = _A_[i.index()].begin(); j != _A_[i.index()].end() ; ++j) for (ColIterator j = (*_A_)[i.index()].begin(); j != (*_A_)[i.index()] .end(); ++j)
if (mask[j.index()] != 2) if (mask[j.index()] != 2)
(*j).usmv(alpha,x[j.index()],y[i.index()]); (*j).usmv(alpha,x[j.index()],y[i.index()]);
} }
} }
} }
//! Category of the linear operator (see SolverCategory::Category) //! Category of the linear operator (see SolverCategory::Category)
virtual SolverCategory::Category category() const virtual SolverCategory::Category category() const
{ {
return SolverCategory::nonoverlapping; return SolverCategory::nonoverlapping;
} }
private: private:
const matrix_type&; _A_; std::shared_ptr<const matrix_type>; _A_;
const communication_type& communication; const communication_type& communication;
mutable bool buildcomm; mutable bool buildcomm;
mutable std::vector<double> mask; mutable std::vector<double> mask;
mutable std::multimap<int,int> bordercontribution; mutable std::multimap<int,int> bordercontribution;
}; };
/** @} */ /** @} */
namespace Amg namespace Amg
{ {
skipping to change at line 266 skipping to change at line 270
* @brief Nonoverlapping parallel preconditioner. * @brief Nonoverlapping parallel preconditioner.
* *
* This is essentially a wrapper that take a sequential * This is essentially a wrapper that take a sequential
* preconditoner. In each step the sequential preconditioner * preconditoner. In each step the sequential preconditioner
* is applied and then all owner data points are updated on * is applied and then all owner data points are updated on
* all other processes. * all other processes.
*/ */
template<class C, class P> template<class C, class P>
class NonoverlappingBlockPreconditioner class NonoverlappingBlockPreconditioner
: public Dune::Preconditioner<typename P::domain_type,typename P::range_type : public Preconditioner<typename P::domain_type,typename P::range_type> {
> { friend struct Amg::ConstructionTraits<NonoverlappingBlockPreconditioner<C,P>
friend class Amg::ConstructionTraits<NonoverlappingBlockPreconditioner<C,P> >;
>; using X = typename P::domain_type;
using Y = typename P::range_type;
public: public:
//! \brief The domain type of the preconditioner. //! \brief The domain type of the preconditioner.
typedef typename P::domain_type domain_type; typedef typename P::domain_type domain_type;
//! \brief The range type of the preconditioner. //! \brief The range type of the preconditioner.
typedef typename P::range_type range_type; typedef typename P::range_type range_type;
//! \brief The type of the communication object. //! \brief The type of the communication object.
typedef C communication_type; typedef C communication_type;
/*! \brief Constructor. /*! \brief Constructor.
constructor gets all parameters to operate the prec. constructor gets all parameters to operate the prec.
\param prec The sequential preconditioner. \param prec The sequential preconditioner.
\param c The communication object for syncing owner and copy \param c The communication object for syncing owner and copy
data points. (E.~g. OwnerOverlapCommunication ) data points. (E.~g. OwnerOverlapCommunication )
*/ */
NonoverlappingBlockPreconditioner (P& prec, const communication_type& c) /*! \brief Constructor.
: preconditioner(prec), communication(c)
{} constructor gets all parameters to operate the prec.
\param p The sequential preconditioner.
\param c The communication object for syncing overlap and copy
data points. (E.~g. OwnerOverlapCopyCommunication )
*/
NonoverlappingBlockPreconditioner (P& p, const communication_type& c)
: _preconditioner(stackobject_to_shared_ptr(p)), _communication(c)
{ }
/*! \brief Constructor.
constructor gets all parameters to operate the prec.
\param p The sequential preconditioner.
\param c The communication object for syncing overlap and copy
data points. (E.~g. OwnerOverlapCopyCommunication )
*/
NonoverlappingBlockPreconditioner (const std::shared_ptr<P>& p, const commun
ication_type& c)
: _preconditioner(p), _communication(c)
{ }
/*! /*!
\brief Prepare the preconditioner. \brief Prepare the preconditioner.
\copydoc Preconditioner::pre(domain_type&,range_type&) \copydoc Preconditioner::pre(domain_type&,range_type&)
*/ */
virtual void pre (domain_type& x, range_type& b) virtual void pre (domain_type& x, range_type& b)
{ {
preconditioner.pre(x,b); _preconditioner->pre(x,b);
} }
/*! /*!
\brief Apply the preconditioner \brief Apply the preconditioner
\copydoc Preconditioner::apply(domain_type&,const range_type&) \copydoc Preconditioner::apply(domain_type&,const range_type&)
*/ */
virtual void apply (domain_type& v, const range_type& d) virtual void apply (domain_type& v, const range_type& d)
{ {
// block preconditioner equivalent to WrappedPreconditioner from // block preconditioner equivalent to WrappedPreconditioner from
// pdelab/backend/ovlpistsolverbackend.hh, // pdelab/backend/ovlpistsolverbackend.hh,
// but not to BlockPreconditioner from schwarz.hh // but not to BlockPreconditioner from schwarz.hh
preconditioner.apply(v,d); _preconditioner->apply(v,d);
communication.addOwnerCopyToOwnerCopy(v,v); _communication.addOwnerCopyToOwnerCopy(v,v);
}
template<bool forward>
void apply (X& v, const Y& d)
{
_preconditioner->template apply<forward>(v,d);
_communication.addOwnerCopyToOwnerCopy(v,v);
} }
/*! /*!
\brief Clean up. \brief Clean up.
\copydoc Preconditioner::post(domain_type&) \copydoc Preconditioner::post(domain_type&)
*/ */
virtual void post (domain_type& x) virtual void post (domain_type& x)
{ {
preconditioner.post(x); _preconditioner->post(x);
} }
//! Category of the preconditioner (see SolverCategory::Category) //! Category of the preconditioner (see SolverCategory::Category)
virtual SolverCategory::Category category() const virtual SolverCategory::Category category() const
{ {
return SolverCategory::nonoverlapping; return SolverCategory::nonoverlapping;
} }
private: private:
//! \brief a sequential preconditioner //! \brief a sequential preconditioner
P& preconditioner; std::shared_ptr<P> _preconditioner;
//! \brief the communication object //! \brief the communication object
const communication_type& communication; const communication_type& _communication;
}; };
/** @} end documentation */ /** @} end documentation */
} // end namespace } // end namespace
#endif #endif
 End of changes. 16 change blocks. 
21 lines changed or deleted 53 lines changed or added

Home  |  About  |  Features  |  All  |  Newest  |  Dox  |  Diffs  |  RSS Feeds  |  Screenshots  |  Comments  |  Imprint  |  Privacy  |  HTTP(S)