DUNE-FEM (unstable)

ilusubdomainsolver.hh
Go to the documentation of this file.
1// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
2// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
3// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
4// vi: set et ts=4 sw=2 sts=2:
5#ifndef DUNE_ISTL_ILUSUBDOMAINSOLVER_HH
6#define DUNE_ISTL_ILUSUBDOMAINSOLVER_HH
7
8#include <map>
11#include "matrix.hh"
12#include <cmath>
13#include <cstdlib>
14
15namespace Dune {
16
35 template<class M, class X, class Y>
37 public:
39 typedef typename std::remove_const<M>::type matrix_type;
41 typedef X domain_type;
43 typedef Y range_type;
44
51 virtual void apply (X& v, const Y& d) =0;
52
53 virtual ~ILUSubdomainSolver()
54 {}
55
56 protected:
62 template<class S>
63 std::size_t copyToLocalMatrix(const M& A, S& rowset);
64
66 // for ILUN
68 };
69
76 template<class M, class X, class Y>
78 : public ILUSubdomainSolver<M,X,Y>{
79 public:
81 typedef typename std::remove_const<M>::type matrix_type;
82 typedef typename std::remove_const<M>::type rilu_type;
84 typedef X domain_type;
86 typedef Y range_type;
87
88
93 void apply (X& v, const Y& d)
94 {
95 ILU::blockILUBacksolve(this->ILU,v,d);
96 }
104 template<class S>
105 void setSubMatrix(const M& A, S& rowset);
106
107 };
108
109 template<class M, class X, class Y>
110 class ILUNSubdomainSolver
111 : public ILUSubdomainSolver<M,X,Y>{
112 public:
114 typedef typename std::remove_const<M>::type matrix_type;
115 typedef typename std::remove_const<M>::type rilu_type;
117 typedef X domain_type;
119 typedef Y range_type;
120
125 void apply (X& v, const Y& d)
126 {
127 ILU::blockILUBacksolve(RILU,v,d);
128 }
129
137 template<class S>
138 void setSubMatrix(const M& A, S& rowset);
139
140 private:
144 rilu_type RILU;
145 };
146
147
148
149 template<class M, class X, class Y>
150 template<class S>
151 std::size_t ILUSubdomainSolver<M,X,Y>::copyToLocalMatrix(const M& A, S& rowSet)
152 {
153 // Calculate consecutive indices for local problem
154 // while preserving the ordering
155 typedef typename M::size_type size_type;
156 typedef std::map<typename S::value_type,size_type> IndexMap;
157 typedef typename IndexMap::iterator IMIter;
158 IndexMap indexMap;
159 IMIter guess = indexMap.begin();
160 size_type localIndex=0;
161
162 typedef typename S::const_iterator SIter;
163 for(SIter rowIdx = rowSet.begin(), rowEnd=rowSet.end();
164 rowIdx!= rowEnd; ++rowIdx, ++localIndex)
165 guess = indexMap.insert(guess,
166 std::make_pair(*rowIdx,localIndex));
167
168
169 // Build Matrix for local subproblem
170 ILU.setSize(rowSet.size(),rowSet.size());
171 ILU.setBuildMode(matrix_type::row_wise);
172
173 // Create sparsity pattern
174 typedef typename matrix_type::CreateIterator CIter;
175 CIter rowCreator = ILU.createbegin();
176 std::size_t offset=0;
177 for(SIter rowIdx = rowSet.begin(), rowEnd=rowSet.end();
178 rowIdx!= rowEnd; ++rowIdx, ++rowCreator) {
179 // See which row entries are in our subset and add them to
180 // the sparsity pattern
181 guess = indexMap.begin();
182
183 for(typename matrix_type::ConstColIterator col=A[*rowIdx].begin(),
184 endcol=A[*rowIdx].end(); col != endcol; ++col) {
185 // search for the entry in the row set
186 guess = indexMap.find(col.index());
187 if(guess!=indexMap.end()) {
188 // add local index to row
189 rowCreator.insert(guess->second);
190 offset=std::max(offset,(std::size_t)std::abs((int)(guess->second-rowCreator.index())));
191 }
192 }
193
194 }
195
196 // Insert the matrix values for the local problem
197 typename matrix_type::iterator iluRow=ILU.begin();
198
199 for(SIter rowIdx = rowSet.begin(), rowEnd=rowSet.end();
200 rowIdx!= rowEnd; ++rowIdx, ++iluRow) {
201 // See which row entries are in our subset and add them to
202 // the sparsity pattern
203 typename matrix_type::ColIterator localCol=iluRow->begin();
204 for(typename matrix_type::ConstColIterator col=A[*rowIdx].begin(),
205 endcol=A[*rowIdx].end(); col != endcol; ++col) {
206 // search for the entry in the row set
207 guess = indexMap.find(col.index());
208 if(guess!=indexMap.end()) {
209 // set local value
210 (*localCol)=(*col);
211 ++localCol;
212 }
213 }
214 }
215 return offset;
216 }
217
218
219 template<class M, class X, class Y>
220 template<class S>
222 {
223 this->copyToLocalMatrix(A,rowSet);
224 ILU::blockILU0Decomposition(this->ILU);
225 }
226
227 template<class M, class X, class Y>
228 template<class S>
229 void ILUNSubdomainSolver<M,X,Y>::setSubMatrix(const M& A, S& rowSet)
230 {
231 std::size_t offset=copyToLocalMatrix(A,rowSet);
232 RILU.setSize(rowSet.size(),rowSet.size(), (1+2*offset)*rowSet.size());
233 RILU.setBuildMode(matrix_type::row_wise);
234 ILU::blockILUDecomposition(this->ILU, (offset+1)/2, RILU);
235 }
236
238} // end name space DUNE
239
240
241#endif
Exact subdomain solver using ILU(p) with appropriate p.
Definition: ilusubdomainsolver.hh:78
X domain_type
The domain type of the preconditioner.
Definition: ilusubdomainsolver.hh:84
Y range_type
The range type of the preconditioner.
Definition: ilusubdomainsolver.hh:86
std::remove_const< M >::type matrix_type
The matrix type the preconditioner is for.
Definition: ilusubdomainsolver.hh:81
void apply(X &v, const Y &d)
Apply the subdomain solver.
Definition: ilusubdomainsolver.hh:93
base class encapsulating common algorithms of ILU0SubdomainSolver and ILUNSubdomainSolver.
Definition: ilusubdomainsolver.hh:36
matrix_type ILU
The ILU0 decomposition of the matrix, or the local matrix.
Definition: ilusubdomainsolver.hh:67
X domain_type
The domain type of the preconditioner.
Definition: ilusubdomainsolver.hh:41
Y range_type
The range type of the preconditioner.
Definition: ilusubdomainsolver.hh:43
std::remove_const< M >::type matrix_type
The matrix type the preconditioner is for.
Definition: ilusubdomainsolver.hh:39
virtual void apply(X &v, const Y &d)=0
Apply the subdomain solver.
std::size_t copyToLocalMatrix(const M &A, S &rowset)
Copy the local part of the global matrix to ILU.
Definition: ilusubdomainsolver.hh:151
void setSubMatrix(const M &A, S &rowset)
Set the data of the local problem.
Definition: ilusubdomainsolver.hh:221
A dynamic dense block matrix class.
Dune namespace.
Definition: alignedallocator.hh:13
Define general preconditioner interface.
Traits for type conversions and type information.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 21, 23:30, 2024)