Loading...
Searching...
No Matches
ilusubdomainsolver.hh
Go to the documentation of this file.
// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_ISTL_ILUSUBDOMAINSOLVER_HH
#define DUNE_ISTL_ILUSUBDOMAINSOLVER_HH
#include <map>
#include <dune/common/typetraits.hh>
#include <dune/istl/preconditioners.hh>
#include "matrix.hh"
#include <cmath>
#include <cstdlib>
namespace Dune {
template<class M, class X, class Y>
class ILUSubdomainSolver {
public:
typedef typename std::remove_const<M>::type matrix_type;
typedef X domain_type;
typedef Y range_type;
virtual void apply (X& v, const Y& d) =0;
virtual ~ILUSubdomainSolver()
{}
protected:
template<class S>
std::size_t copyToLocalMatrix(const M& A, S& rowset);
// for ILUN
matrix_type ILU;
};
template<class M, class X, class Y>
class ILU0SubdomainSolver
: public ILUSubdomainSolver<M,X,Y>{
public:
typedef typename std::remove_const<M>::type matrix_type;
typedef typename std::remove_const<M>::type rilu_type;
typedef X domain_type;
typedef Y range_type;
void apply (X& v, const Y& d)
{
ILU::blockILUBacksolve(this->ILU,v,d);
}
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:
typedef typename std::remove_const<M>::type matrix_type;
typedef typename std::remove_const<M>::type rilu_type;
typedef X domain_type;
typedef Y range_type;
void apply (X& v, const Y& d)
{
ILU::blockILUBacksolve(RILU,v,d);
}
template<class S>
void setSubMatrix(const M& A, S& rowset);
private:
rilu_type RILU;
};
template<class M, class X, class Y>
template<class S>
std::size_t ILUSubdomainSolver<M,X,Y>::copyToLocalMatrix(const M& A, S& rowSet)
{
// Calculate consecutive indices for local problem
// while preserving the ordering
typedef typename M::size_type size_type;
typedef std::map<typename S::value_type,size_type> IndexMap;
typedef typename IndexMap::iterator IMIter;
IndexMap indexMap;
IMIter guess = indexMap.begin();
size_type localIndex=0;
typedef typename S::const_iterator SIter;
for(SIter rowIdx = rowSet.begin(), rowEnd=rowSet.end();
rowIdx!= rowEnd; ++rowIdx, ++localIndex)
guess = indexMap.insert(guess,
std::make_pair(*rowIdx,localIndex));
// Build Matrix for local subproblem
ILU.setSize(rowSet.size(),rowSet.size());
ILU.setBuildMode(matrix_type::row_wise);
// Create sparsity pattern
typedef typename matrix_type::CreateIterator CIter;
CIter rowCreator = ILU.createbegin();
std::size_t offset=0;
for(SIter rowIdx = rowSet.begin(), rowEnd=rowSet.end();
rowIdx!= rowEnd; ++rowIdx, ++rowCreator) {
// See which 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()) {
// 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
typename matrix_type::iterator iluRow=ILU.begin();
for(SIter rowIdx = rowSet.begin(), rowEnd=rowSet.end();
rowIdx!= rowEnd; ++rowIdx, ++iluRow) {
// See which row entries are in our subset and add them to
// the sparsity pattern
typename matrix_type::ColIterator localCol=iluRow->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()) {
// set local value
(*localCol)=(*col);
++localCol;
}
}
}
return offset;
}
template<class M, class X, class Y>
template<class S>
void ILU0SubdomainSolver<M,X,Y>::setSubMatrix(const M& A, S& rowSet)
{
this->copyToLocalMatrix(A,rowSet);
ILU::blockILU0Decomposition(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(), (1+2*offset)*rowSet.size());
RILU.setBuildMode(matrix_type::row_wise);
ILU::blockILUDecomposition(this->ILU, (offset+1)/2, RILU);
}
} // end name space DUNE
#endif
preconditioners.hh
Define general preconditioner interface.
matrix.hh
A dynamic dense block matrix class.
Dune::ILUSubdomainSolver::copyToLocalMatrix
std::size_t copyToLocalMatrix(const M &A, S &rowset)
Copy the local part of the global matrix to ILU.
Definition ilusubdomainsolver.hh:151
Dune::ILUNSubdomainSolver::setSubMatrix
void setSubMatrix(const M &A, S &rowset)
Set the data of the local problem.
Definition ilusubdomainsolver.hh:229
Dune::ILU0SubdomainSolver::setSubMatrix
void setSubMatrix(const M &A, S &rowset)
Set the data of the local problem.
Definition ilusubdomainsolver.hh:221
Dune
Definition allocator.hh:11
Dune::col
MatrixSparsityPatternGatherScatter< M, I >::ColIter MatrixSparsityPatternGatherScatter< M, I >::col
Definition matrixredistribute.hh:591
Dune::ILU::blockILUBacksolve
void blockILUBacksolve(const M &A, X &v, const Y &d)
LU backsolve with stored inverse.
Definition ilu.hh:103
Dune::ILU::blockILU0Decomposition
void blockILU0Decomposition(M &A)
compute ILU decomposition of A. A is overwritten by its decomposition
Definition ilu.hh:35
Dune::ILU::blockILUDecomposition
void blockILUDecomposition(const M &A, int n, M &ILU)
Definition ilu.hh:176
Dune::size_type
Dune::BCRSMatrix::begin
Iterator begin()
Get iterator to first row.
Definition bcrsmatrix.hh:672
Dune::ILUSubdomainSolver
base class encapsulating common algorithms of ILU0SubdomainSolver and ILUNSubdomainSolver.
Definition ilusubdomainsolver.hh:36
Dune::ILUSubdomainSolver::ILU
matrix_type ILU
The ILU0 decomposition of the matrix, or the local matrix.
Definition ilusubdomainsolver.hh:67
Dune::ILUSubdomainSolver::domain_type
X domain_type
The domain type of the preconditioner.
Definition ilusubdomainsolver.hh:41
Dune::ILUSubdomainSolver::~ILUSubdomainSolver
virtual ~ILUSubdomainSolver()
Definition ilusubdomainsolver.hh:53
Dune::ILUSubdomainSolver::range_type
Y range_type
The range type of the preconditioner.
Definition ilusubdomainsolver.hh:43
Dune::ILUSubdomainSolver::matrix_type
std::remove_const< M >::type matrix_type
The matrix type the preconditioner is for.
Definition ilusubdomainsolver.hh:39
Dune::ILUSubdomainSolver::apply
virtual void apply(X &v, const Y &d)=0
Apply the subdomain solver.
Dune::ILU0SubdomainSolver
Exact subdomain solver using ILU(p) with appropriate p.
Definition ilusubdomainsolver.hh:78
Dune::ILU0SubdomainSolver::domain_type
X domain_type
The domain type of the preconditioner.
Definition ilusubdomainsolver.hh:84
Dune::ILU0SubdomainSolver::range_type
Y range_type
The range type of the preconditioner.
Definition ilusubdomainsolver.hh:86
Dune::ILU0SubdomainSolver::rilu_type
std::remove_const< M >::type rilu_type
Definition ilusubdomainsolver.hh:82
Dune::ILU0SubdomainSolver::matrix_type
std::remove_const< M >::type matrix_type
The matrix type the preconditioner is for.
Definition ilusubdomainsolver.hh:81
Dune::ILU0SubdomainSolver::apply
void apply(X &v, const Y &d)
Apply the subdomain solver.
Definition ilusubdomainsolver.hh:93
Dune::ILUNSubdomainSolver
Definition ilusubdomainsolver.hh:111
Dune::ILUNSubdomainSolver::domain_type
X domain_type
The domain type of the preconditioner.
Definition ilusubdomainsolver.hh:117
Dune::ILUNSubdomainSolver::matrix_type
std::remove_const< M >::type matrix_type
The matrix type the preconditioner is for.
Definition ilusubdomainsolver.hh:114
Dune::ILUNSubdomainSolver::rilu_type
std::remove_const< M >::type rilu_type
Definition ilusubdomainsolver.hh:115
Dune::ILUNSubdomainSolver::apply
void apply(X &v, const Y &d)
Apply the subdomain solver.
Definition ilusubdomainsolver.hh:125
Dune::ILUNSubdomainSolver::range_type
Y range_type
The range type of the preconditioner.
Definition ilusubdomainsolver.hh:119