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

Added ILU(p) subdomain solver and documentation. 

[[Imported from SVN: r1263]]
parent 5c492b68
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dune/istl/ilusubdomainsolver.hh
+ 145
− 10
View file @ 87e8c2fb
......@@ -5,11 +5,30 @@
#include <map>
#include <dune/common/typetraits.hh>
#include "matrix.hh"
namespace Dune {
/**
* @file
* @brief Various local subdomain solvers based on ILU
* for SeqOverlappingSchwarz.
* @author Markus Blatt
*/
/**
* @addtogroup ISTL
* @{
*/
/**
* @brief base class encapsulating common algorithms of ILU0SubdomainSolver
* and ILUNSubdomainSolver.
* @tparam M The type of the matrix.
* @tparam X The type of the vector for the domain.
* @tparam X The type of the vector for the range.
*
*/
template<class M, class X, class Y>
class ILU0SubdomainSolver {
class ILUSubdomainSolver {
public:
//! \brief The matrix type the preconditioner is for.
typedef typename Dune::remove_const<M>::type matrix_type;
......@@ -17,26 +36,114 @@ namespace Dune {
typedef X domain_type;
//! \brief The range type of the preconditioner.
typedef Y range_type;
typedef typename X::field_type field_type;
/**
* @brief Apply the subdomain solver.
*
* On entry v=? and d=b-A(x) (although this might not be
* computed in that way. On exit v contains the update
*/
virtual void apply (X& v, const Y& d) =0;
virtual ~ILUSubdomainSolver()
{}
protected:
/**
* @brief Copy the local part of the global matrix to ILU.
* @param A The global matrix.
* @param rowset The global indices of the local problem.
*/
template<class S>
std::size_t copyToLocalMatrix(const M& A, S& rowset);
//! \brief The ILU0 decomposition of the matrix, or the local matrix
// for ILUN
matrix_type ILU;
};
/**
* @brief Exact subdomain solver using ILU(p) with appropriate p.
* @tparam M The type of the matrix.
* @tparam X The type of the vector for the domain.
* @tparam X The type of the vector for the range.
*/
template<class M, class X, class Y>
class ILU0SubdomainSolver
: public ILUSubdomainSolver<M,X,Y>{
public:
//! \brief The matrix type the preconditioner is for.
typedef typename Dune::remove_const<M>::type matrix_type;
typedef typename Dune::remove_const<M>::type rilu_type;
//! \brief The domain type of the preconditioner.
typedef X domain_type;
//! \brief The range type of the preconditioner.
typedef Y range_type;
/**
* @brief Apply the subdomain solver.
* @copydoc ILUSubdomainSolver::apply
*/
void apply (X& v, const Y& d)
{
bilu_backsolve(ILU,v,d);
bilu_backsolve(this->ILU,v,d);
}
/**
* @brief Set the data of the local problem.
*
* @param A The global matrix.
* @param rowset The global indices of the local problem.
* @tparam S The type of the set with the indices.
*/
template<class S>
void setSubMatrix(const M& A, S& rowset);
};
template<class M, class X, class Y>
class ILUNSubdomainSolver
: public ILUSubdomainSolver<M,X,Y>{
public:
//! \brief The matrix type the preconditioner is for.
typedef typename Dune::remove_const<M>::type matrix_type;
typedef typename Dune::remove_const<M>::type rilu_type;
//! \brief The domain type of the preconditioner.
typedef X domain_type;
//! \brief The range type of the preconditioner.
typedef Y range_type;
/**
* @brief Apply the subdomain solver.
* @copydoc ILUSubdomainSolver::apply
*/
void apply (X& v, const Y& d)
{
bilu_backsolve(RILU,v,d);
}
/**
* @brief Set the data of the local problem.
*
* @param A The global matrix.
* @param rowset The global indices of the local problem.
* @tparam S The type of the set with the indices.
*/
template<class S>
void setSubMatrix(const M& A, S& rowset);
private:
//! \brief The relaxation factor to use.
field_type _w;
//! \brief The ILU0 decomposition of the matrix.
matrix_type ILU;
/**
* @brief Storage for the ILUN decomposition.
*/
rilu_type RILU;
};
template<class M, class X, class Y>
template<class S>
void ILU0SubdomainSolver<M,X,Y>::setSubMatrix(const M& A, S& rowSet)
std::size_t ILUSubdomainSolver<M,X,Y>::copyToLocalMatrix(const M& A, S& rowSet)
{
// Calculate consecutive indices for local problem
// while perserving the ordering
......@@ -53,6 +160,9 @@ namespace Dune {
guess = indexMap.insert(guess,
std::make_pair(*rowIdx,localIndex));
// Build Matrix for local subproblem
ILU.setSize(rowSet.size(),rowSet.size());
// Build Matrix for local subproblem
ILU.setSize(rowSet.size(),rowSet.size());
ILU.setBuildMode(matrix_type::row_wise);
......@@ -61,19 +171,24 @@ namespace Dune {
typedef typename matrix_type::CreateIterator CIter;
CIter rowCreator = ILU.createbegin();
typedef typename S::const_iterator RIter;
std::size_t offset=0;
for(SIter rowIdx = rowSet.begin(), rowEnd=rowSet.end();
rowIdx!= rowEnd; ++rowIdx, ++rowCreator) {
// See wich row entries are in our subset and add them to
// the sparsity pattern
guess = indexMap.begin();
for(typename matrix_type::ConstColIterator col=A[*rowIdx].begin(),
endcol=A[*rowIdx].end(); col != endcol; ++col) {
// search for the entry in the row set
guess = indexMap.find(col.index());
if(guess!=indexMap.end())
if(guess!=indexMap.end()) {
// add local index to row
rowCreator.insert(guess->second);
offset=std::max(offset,(std::size_t)std::abs((int)(guess->second-rowCreator.index())));
}
}
}
// Insert the matrix values for the local problem
......@@ -95,9 +210,29 @@ namespace Dune {
}
}
}
bilu0_decomposition(ILU);
return offset;
}
template<class M, class X, class Y>
template<class S>
void ILU0SubdomainSolver<M,X,Y>::setSubMatrix(const M& A, S& rowSet)
{
copyToLocalMatrix(A,rowSet);
bilu0_decomposition(this->ILU);
}
template<class M, class X, class Y>
template<class S>
void ILUNSubdomainSolver<M,X,Y>::setSubMatrix(const M& A, S& rowSet)
{
std::size_t offset=copyToLocalMatrix(A,rowSet);
RILU.setSize(rowSet.size(),rowSet.size(), (this->ILU.nonzeroes()/rowSet.size()+2*offset)*rowSet.size());
RILU.setBuildMode(matrix_type::row_wise);
bilu_decomposition(this->ILU, (offset+1)/2, RILU);
}
/** @} */
} // end name space DUNE
......
dune/istl/overlappingschwarz.hh
+ 138
− 23
View file @ 87e8c2fb
......@@ -151,15 +151,34 @@ namespace Dune
*/
OverlappingAssigner(std::size_t maxlength, const BCRSMatrix<FieldMatrix<T,n,m>,A>& mat,
const range_type& b, range_type& x);
/**
* @brief Deallocates memory of the local vector.
* @warning memory is released by the destructor as this Functor
* is copied and the copy needs to still have the data.
*/
void deallocate();
/*
* @brief Resets the local index to zero.
*/
void resetIndexForNextDomain();
/**
* @brief Get the local left hand side.
* @return The local left hand side.
*/
field_type* lhs();
/**
* @brief Get the local left hand side.
* @return The local left hand side.
*/
field_type* rhs();
/**
* @brief relax the result.
* @param relax The relaxation parameter.
*/
void relaxResult(field_type relax);
/**
......@@ -168,22 +187,38 @@ namespace Dune
*/
void operator()(const size_type& domain);
/**
* @brief Assigns the block to the current local
* index.
* At the same time the local defect is calculated
* for the index and stored in the rhs.
* Afterwards the is incremented for the next block.
*/
void assignResult(block_type& res);
private:
/**
* @brief The global matrix for the defect calculation.
*/
const matrix_type* mat;
/** @brief The local right hand side. */
field_type* rhs_;
/** @brief The local left hand side. */
field_type* lhs_;
/** @brief The global right hand side for the defect calculation. */
const range_type* b;
/** @brief The global left hand side for adding the local update to. */
range_type* x;
/** @brief The current local index. */
std::size_t i;
/** @brief The maximum entries over all subdomains. */
std::size_t maxlength_;
};
#endif
template<class M, class X, class Y>
class OverlappingAssigner<ILU0SubdomainSolver<M,X,Y> >
class OverlappingAssignerILUBase
{
public:
typedef M matrix_type;
......@@ -200,18 +235,35 @@ namespace Dune
* @param b the global right hand side.
* @param x the global left hand side.
*/
OverlappingAssigner(std::size_t maxlength, const M& mat,
const Y& b, X& x);
OverlappingAssignerILUBase(std::size_t maxlength, const M& mat,
const Y& b, X& x);
/**
* @brief Deallocates memory of the local vector.
* @warning memory is released by the destructor as this Functor
* is copied and the copy needs to still have the data.
*/
void deallocate();
/*
* @brief Resets the local index to zero.
*/
void resetIndexForNextDomain();
/*
* @brief Resets the local index to zero.
*/
X& lhs();
/**
* @brief Get the local left hand side.
* @return The local left hand side.
*/
Y& rhs();
/**
* @brief relax the result.
* @param relax The relaxation parameter.
*/
void relaxResult(field_type relax);
/**
......@@ -220,16 +272,69 @@ namespace Dune
*/
void operator()(const size_type& domain);
/**
* @brief Assigns the block to the current local
* index.
* At the same time the local defect is calculated
* for the index and stored in the rhs.
* Afterwards the is incremented for the next block.
*/
void assignResult(block_type& res);
private:
/**
* @brief The global matrix for the defect calculation.
*/
const M* mat;
/** @brief The local right hand side. */
X* lhs_;
/** @brief The local left hand side. */
Y* rhs_;
/** @brief The global right hand side for the defect calculation. */
const Y* b;
/** @brief The global left hand side for adding the local update to. */
X* x;
/** @brief The maximum entries over all subdomains. */
size_type i;
};
// specialization for ILU0
template<class M, class X, class Y>
class OverlappingAssigner<ILU0SubdomainSolver<M,X,Y> >
: public OverlappingAssignerILUBase<M,X,Y>
{
public:
/**
* @brief Constructor.
* @param mat The global matrix.
* @param rhs storage for the local defect.
* @param b the global right hand side.
* @param x the global left hand side.
*/
OverlappingAssigner(std::size_t maxlength, const M& mat,
const Y& b, X& x)
: OverlappingAssignerILUBase<M,X,Y>(maxlength, mat,b,x)
{}
};
// specialization for ILUN
template<class M, class X, class Y>
class OverlappingAssigner<ILUNSubdomainSolver<M,X,Y> >
: public OverlappingAssignerILUBase<M,X,Y>
{
public:
/**
* @brief Constructor.
* @param mat The global matrix.
* @param rhs storage for the local defect.
* @param b the global right hand side.
* @param x the global left hand side.
*/
OverlappingAssigner(std::size_t maxlength, const M& mat,
const Y& b, X& x)
: OverlappingAssignerILUBase<M,X,Y>(maxlength, mat,b,x)
{}
};
template<typename S, typename T>
struct AdditiveAdder
......@@ -420,7 +525,7 @@ namespace Dune
#endif
template<class M,class X, class Y>
struct SeqOverlappingSchwarzAssembler<ILU0SubdomainSolver<M,X,Y> >
struct SeqOverlappingSchwarzAssemblerILUBase
{
typedef M matrix_type;
template<class RowToDomain, class Solvers, class SubDomains>
......@@ -429,6 +534,16 @@ namespace Dune
bool onTheFly);
};
template<class M,class X, class Y>
struct SeqOverlappingSchwarzAssembler<ILU0SubdomainSolver<M,X,Y> >
: public SeqOverlappingSchwarzAssemblerILUBase<M,X,Y>
{};
template<class M,class X, class Y>
struct SeqOverlappingSchwarzAssembler<ILUNSubdomainSolver<M,X,Y> >
: public SeqOverlappingSchwarzAssemblerILUBase<M,X,Y>
{};
/**
* @brief Sequential overlapping Schwarz preconditioner
*
......@@ -866,11 +981,11 @@ namespace Dune
template<class M,class X,class Y>
template<class RowToDomain, class Solvers, class SubDomains>
std::size_t SeqOverlappingSchwarzAssembler<ILU0SubdomainSolver<M,X,Y> >::assembleLocalProblems(const RowToDomain& rowToDomain,
const matrix_type& mat,
Solvers& solvers,
const SubDomains& subDomains,
bool onTheFly)
std::size_t SeqOverlappingSchwarzAssemblerILUBase<M,X,Y>::assembleLocalProblems(const RowToDomain& rowToDomain,
const matrix_type& mat,
Solvers& solvers,
const SubDomains& subDomains,
bool onTheFly)
{
typedef typename SubDomains::const_iterator DomainIterator;
typedef typename Solvers::iterator SolverIterator;
......@@ -1046,10 +1161,10 @@ namespace Dune
#endif
template<class M, class X, class Y>
OverlappingAssigner<ILU0SubdomainSolver<M,X,Y> >::OverlappingAssigner(std::size_t maxlength,
const M& mat_,
const Y& b_,
X& x_)
OverlappingAssignerILUBase<M,X,Y>::OverlappingAssignerILUBase(std::size_t maxlength,
const M& mat_,
const Y& b_,
X& x_)
: mat(&mat_),
b(&b_),
x(&x_), i(0)
......@@ -1059,14 +1174,14 @@ namespace Dune
}
template<class M, class X, class Y>
void OverlappingAssigner<ILU0SubdomainSolver<M,X,Y> >::deallocate()
void OverlappingAssignerILUBase<M,X,Y>::deallocate()
{
delete rhs_;
delete lhs_;
}
template<class M, class X, class Y>
void OverlappingAssigner<ILU0SubdomainSolver<M,X,Y> >::operator()(const size_type& domainIndex)
void OverlappingAssignerILUBase<M,X,Y>::operator()(const size_type& domainIndex)
{
(*rhs_)[i]=(*b)[domainIndex];
......@@ -1082,31 +1197,31 @@ namespace Dune
}
template<class M, class X, class Y>
void OverlappingAssigner<ILU0SubdomainSolver<M,X,Y> >::relaxResult(field_type relax)
void OverlappingAssignerILUBase<M,X,Y>::relaxResult(field_type relax)
{
(*lhs_)[i]*=relax;
}
template<class M, class X, class Y>
void OverlappingAssigner<ILU0SubdomainSolver<M,X,Y> >::assignResult(block_type& res)
void OverlappingAssignerILUBase<M,X,Y>::assignResult(block_type& res)
{
res+=(*lhs_)[i++];
}
template<class M, class X, class Y>
X& OverlappingAssigner<ILU0SubdomainSolver<M,X,Y> >::lhs()
X& OverlappingAssignerILUBase<M,X,Y>::lhs()
{
return *lhs_;
}
template<class M, class X, class Y>
Y& OverlappingAssigner<ILU0SubdomainSolver<M,X,Y> >::rhs()
Y& OverlappingAssignerILUBase<M,X,Y>::rhs()
{
return *rhs_;
}
template<class M, class X, class Y>
void OverlappingAssigner<ILU0SubdomainSolver<M,X,Y> >::resetIndexForNextDomain()
void OverlappingAssignerILUBase<M,X,Y>::resetIndexForNextDomain()
{
i=0;
}
......
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