Line data Source code
1 : /*
2 : * Copyright (c) 2010-2015: G-CSC, Goethe University Frankfurt
3 : * Authors: Andreas Vogel, Arne Nägel
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 :
34 : #ifndef __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__ILU__
35 : #define __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__ILU__
36 :
37 : #include <limits>
38 : #include "common/error.h"
39 : #ifndef NDEBUG
40 : #include "common/stopwatch.h"
41 : #endif
42 : #include "common/util/smart_pointer.h"
43 : #include "lib_algebra/operator/interface/preconditioner.h"
44 :
45 : #ifdef UG_PARALLEL
46 : #include "pcl/pcl_util.h"
47 : #include "lib_algebra/parallelization/parallelization_util.h"
48 : #include "lib_algebra/parallelization/matrix_overlap.h"
49 : #include "lib_algebra/parallelization/overlap_writer.h"
50 : #endif
51 :
52 : #include "lib_algebra/ordering_strategies/algorithms/IOrderingAlgorithm.h"
53 : #include "lib_algebra/ordering_strategies/algorithms/native_cuthill_mckee.h" // for backward compatibility
54 :
55 : #include "lib_algebra/algebra_common/permutation_util.h"
56 :
57 : namespace ug{
58 :
59 :
60 : // ILU(0) solver, i.e. static pattern ILU w/ P=P(A)
61 : // (cf. Y Saad, Iterative methods for Sparse Linear Systems, p. 270)
62 : template<typename Matrix_type>
63 : bool FactorizeILU(Matrix_type &A)
64 : {
65 : PROFILE_FUNC_GROUP("algebra ILU");
66 : typedef typename Matrix_type::row_iterator row_iterator;
67 : typedef typename Matrix_type::value_type block_type;
68 :
69 : // for all rows
70 : for(size_t i=1; i < A.num_rows(); i++)
71 : {
72 : // eliminate all entries A(i, k) with k<i with rows A(k, .) and k<i
73 : const row_iterator rowEnd = A.end_row(i);
74 : for(row_iterator it_k = A.begin_row(i); it_k != rowEnd && (it_k.index() < i); ++it_k)
75 : {
76 : const size_t k = it_k.index();
77 : block_type &a_ik = it_k.value();
78 :
79 : // add row k to row i by A(i, .) -= A(k,.) A(i,k) / A(k,k)
80 : // so that A(i,k) is zero.
81 : // save A(i,k)/A(k,k) in A(i,k)
82 : if(fabs(BlockNorm(A(k,k))) < 1e-15*BlockNorm(A(i,k)))
83 : UG_THROW("Diag is Zero for k="<<k<<", cannot factorize ILU.");
84 :
85 : a_ik /= A(k,k);
86 :
87 : row_iterator it_j = it_k;
88 : for(++it_j; it_j != rowEnd; ++it_j)
89 : {
90 : const size_t j = it_j.index();
91 : block_type& a_ij = it_j.value();
92 : bool bFound;
93 : row_iterator p = A.get_connection(k,j, bFound);
94 : if(bFound)
95 : {
96 : const block_type a_kj = p.value();
97 : a_ij -= a_ik *a_kj;
98 : }
99 : }
100 : }
101 : }
102 :
103 : return true;
104 : }
105 :
106 : // ILU(0)-beta solver, i.e.
107 : // -> static pattern ILU w/ P=P(A)
108 : // -> Fill-in is computed and lumped onto the diagonal
109 : template<typename Matrix_type>
110 0 : bool FactorizeILUBeta(Matrix_type &A, number beta)
111 : {
112 : PROFILE_FUNC_GROUP("algebra ILUBeta");
113 : typedef typename Matrix_type::row_iterator row_iterator;
114 : typedef typename Matrix_type::value_type block_type;
115 :
116 : // for all rows i=1:n do
117 0 : for(size_t i=1; i < A.num_rows(); i++)
118 : {
119 : block_type &Aii = A(i,i);
120 0 : block_type Nii(Aii); Nii*=0.0;
121 :
122 : // for k=1:(i-1) do
123 : // eliminate all entries A(i, k) with k<i with rows A(k, .) and k<i
124 : const row_iterator it_iEnd = A.end_row(i);
125 0 : for(row_iterator it_ik = A.begin_row(i); it_ik != it_iEnd && (it_ik.index() < i); ++it_ik)
126 : {
127 :
128 : // add row k to row i by A(i, .) -= [A(i,k) / A(k,k)] A(k,.)
129 : // such that A(i,k) is zero.
130 :
131 : // 1) Contribution to L part:
132 : // store A(i,k)/A(k,k) in A(i,k)
133 : const size_t k = it_ik.index();
134 : block_type &a_ik = it_ik.value();
135 0 : a_ik /= A(k,k);
136 :
137 : // 2) Contribution to U part:
138 : // compute contributions from row k for j=k:N
139 : const row_iterator it_kEnd = A.end_row(k);
140 0 : for (row_iterator it_kj=A.begin_row(k); it_kj != it_kEnd ;++it_kj)
141 : {
142 : const size_t j = it_kj.index();
143 0 : if (j<=k) continue; // index j belongs L part?
144 :
145 : // this index j belongs U part
146 : const block_type& a_kj = it_kj.value();
147 :
148 : bool aijFound;
149 0 : row_iterator pij = A.get_connection(i,j, aijFound);
150 0 : if(aijFound) {
151 : // entry belongs to pattern
152 : // -> proceed with standard elimination
153 : block_type& a_ij = pij.value();
154 0 : a_ij -= a_ik *a_kj ;
155 :
156 : } else {
157 : // entry DOES NOT belong to pattern
158 : // -> we lump it onto the diagonal
159 : // TODO : non square matrices!!!
160 0 : Nii -= a_ik * a_kj;
161 : }
162 :
163 : }
164 : }
165 :
166 : // add fill-in to diagonal
167 : AddMult(Aii, beta, Nii);
168 : }
169 :
170 0 : return true;
171 : }
172 :
173 : template<typename Matrix_type>
174 0 : bool FactorizeILUSorted(Matrix_type &A, const number eps = 1e-50)
175 : {
176 : PROFILE_FUNC_GROUP("algebra ILU");
177 : typedef typename Matrix_type::row_iterator row_iterator;
178 : typedef typename Matrix_type::value_type block_type;
179 :
180 : // for all rows
181 0 : for(size_t i=1; i < A.num_rows(); i++)
182 : {
183 :
184 : // eliminate all entries A(i, k) with k<i with rows A(k, .) and k<i
185 0 : for(row_iterator it_k = A.begin_row(i); it_k != A.end_row(i) && (it_k.index() < i); ++it_k)
186 : {
187 : const size_t k = it_k.index();
188 : block_type &a_ik = it_k.value();
189 : block_type &a_kk = A(k,k);
190 :
191 : // add row k to row i by A(i, .) -= A(k,.) A(i,k) / A(k,k)
192 : // so that A(i,k) is zero.
193 : // safe A(i,k)/A(k,k) in A(i,k)
194 0 : if(fabs(BlockNorm(A(k,k))) < eps * BlockNorm(A(i,k)))
195 0 : UG_THROW("ILU: Blocknorm of diagonal is near-zero for k="<<k<<
196 : " with eps: "<< eps <<", ||A_kk||="<<fabs(BlockNorm(A(k,k)))
197 : <<", ||A_ik||="<<BlockNorm(A(i,k)));
198 :
199 0 : try {a_ik /= a_kk;}
200 0 : UG_CATCH_THROW("Failed to calculate A_ik /= A_kk "
201 : "with i = " << i << " and k = " << k << ".");
202 :
203 :
204 : typename Matrix_type::row_iterator it_ij = it_k; // of row i
205 : ++it_ij; // skip a_ik
206 : typename Matrix_type::row_iterator it_kj = A.begin_row(k); // of row k
207 :
208 0 : while(it_ij != A.end_row(i) && it_kj != A.end_row(k))
209 : {
210 0 : if(it_ij.index() > it_kj.index())
211 : ++it_kj;
212 0 : else if(it_ij.index() < it_kj.index())
213 : ++it_ij;
214 : else
215 : {
216 : block_type &a_ij = it_ij.value();
217 : const block_type &a_kj = it_kj.value();
218 0 : a_ij -= a_ik * a_kj;
219 : ++it_kj; ++it_ij;
220 : }
221 : }
222 : }
223 : }
224 :
225 0 : return true;
226 : }
227 :
228 :
229 : // solve x = L^-1 b
230 : // Returns true on success, or false on issues that lead to some changes in the solution
231 : // (the solution is computed unless no exceptions are thrown)
232 : template<typename Matrix_type, typename Vector_type>
233 0 : bool invert_L(const Matrix_type &A, Vector_type &x, const Vector_type &b)
234 : {
235 : PROFILE_FUNC_GROUP("algebra ILU");
236 : typedef typename Matrix_type::const_row_iterator const_row_iterator;
237 :
238 : typename Vector_type::value_type s;
239 0 : for(size_t i=0; i < x.size(); i++)
240 : {
241 0 : s = b[i];
242 0 : for(const_row_iterator it = A.begin_row(i); it != A.end_row(i); ++it)
243 : {
244 0 : if(it.index() >= i) continue;
245 0 : MatMultAdd(s, 1.0, s, -1.0, it.value(), x[it.index()]);
246 : }
247 0 : x[i] = s;
248 : }
249 :
250 0 : return true;
251 : }
252 :
253 : // solve x = U^-1 * b
254 : // Returns true on success, or false on issues that lead to some changes in the solution
255 : // (the solution is computed unless no exceptions are thrown)
256 : template<typename Matrix_type, typename Vector_type>
257 0 : bool invert_U(const Matrix_type &A, Vector_type &x, const Vector_type &b,
258 : const number eps = 1e-8)
259 : {
260 : PROFILE_FUNC_GROUP("algebra ILU");
261 : typedef typename Matrix_type::const_row_iterator const_row_iterator;
262 :
263 : typename Vector_type::value_type s;
264 :
265 : bool result = true;
266 :
267 : // last row diagonal U entry might be close to zero with corresponding close to zero rhs
268 : // when solving Navier Stokes system, therefore handle separately
269 0 : if(x.size() > 0)
270 : {
271 0 : size_t i=x.size()-1;
272 0 : s = b[i];
273 :
274 : // check if diag part is significantly smaller than rhs
275 : // This may happen when matrix is indefinite with one eigenvalue
276 : // zero. In that case, the factorization on the last row is
277 : // nearly zero due to round-off errors. In order to allow ill-
278 : // scaled matrices (i.e. small matrix entries row-wise) this
279 : // is compared to the rhs, that is small in this case as well.
280 : //TODO: Note that this may happen for problems with naturally
281 : // non-zero kernels, e.g. for the Stokes equation. One should
282 : // probably suppress this message in those cases but set the
283 : // rhs to 0.
284 0 : if (BlockNorm(A(i,i)) <= eps * BlockNorm(s))
285 : {
286 0 : UG_LOG("ILU Warning: Near-zero last diagonal entry "
287 : "with norm "<<BlockNorm(A(i,i))<<" in U "
288 : "for non-near-zero rhs entry with norm "
289 : << BlockNorm(s) << ". Setting rhs to zero.\n"
290 : "NOTE: Reduce 'eps' using e.g. ILU::set_inversion_eps(...) "
291 : "to avoid this warning. Current eps: " << eps << ".\n")
292 : // set correction to zero
293 0 : x[i] = 0;
294 : result = false;
295 : } else {
296 : // c[i] = s/uii;
297 0 : InverseMatMult(x[i], 1.0, A(i,i), s);
298 : }
299 : }
300 0 : if(x.size() <= 1) return result;
301 :
302 : // handle all other rows
303 0 : for(size_t i = x.size()-2; ; --i)
304 : {
305 0 : s = b[i];
306 0 : for(const_row_iterator it = A.begin_row(i); it != A.end_row(i); ++it)
307 : {
308 0 : if(it.index() <= i) continue;
309 : // s -= it.value() * x[it.index()];
310 0 : MatMultAdd(s, 1.0, s, -1.0, it.value(), x[it.index()]);
311 :
312 : }
313 : // x[i] = s/A(i,i);
314 0 : InverseMatMult(x[i], 1.0, A(i,i), s);
315 0 : if(i == 0) break;
316 : }
317 :
318 : return result;
319 : }
320 :
321 : /// ILU / ILU(beta) preconditioner
322 : template <typename TAlgebra>
323 : class ILU : public IPreconditioner<TAlgebra>
324 : {
325 : public:
326 : /// Algebra type
327 : typedef TAlgebra algebra_type;
328 :
329 : /// Vector type
330 : typedef typename TAlgebra::vector_type vector_type;
331 :
332 : /// Matrix type
333 : typedef typename TAlgebra::matrix_type matrix_type;
334 :
335 : /// Matrix Operator type
336 : typedef typename IPreconditioner<TAlgebra>::matrix_operator_type matrix_operator_type;
337 :
338 : /// Base type
339 : typedef IPreconditioner<TAlgebra> base_type;
340 :
341 : /// Ordering type
342 : typedef std::vector<size_t> ordering_container_type;
343 : typedef IOrderingAlgorithm<TAlgebra, ordering_container_type> ordering_algo_type;
344 :
345 : protected:
346 : using base_type::set_debug;
347 : using base_type::debug_writer;
348 : using base_type::write_debug;
349 : using base_type::print_debugger_message;
350 :
351 : public:
352 : // Constructor
353 0 : ILU (double beta=0.0) :
354 0 : m_beta(beta),
355 0 : m_sortEps(1.e-50),
356 0 : m_invEps(1.e-8),
357 0 : m_bDisablePreprocessing(false),
358 0 : m_useConsistentInterfaces(false),
359 0 : m_useOverlap(false),
360 : m_spOrderingAlgo(SPNULL),
361 0 : m_bSortIsIdentity(false),
362 0 : m_u(nullptr)
363 0 : {};
364 :
365 : /// clone constructor
366 0 : ILU (const ILU<TAlgebra> &parent) :
367 : base_type(parent),
368 0 : m_beta(parent.m_beta),
369 0 : m_sortEps(parent.m_sortEps),
370 0 : m_invEps(parent.m_invEps),
371 0 : m_bDisablePreprocessing(parent.m_bDisablePreprocessing),
372 0 : m_useConsistentInterfaces(parent.m_useConsistentInterfaces),
373 0 : m_useOverlap(parent.m_useOverlap),
374 : m_spOrderingAlgo(parent.m_spOrderingAlgo),
375 0 : m_bSortIsIdentity(false),
376 0 : m_u(nullptr)
377 0 : {}
378 :
379 : /// Clone
380 0 : virtual SmartPtr<ILinearIterator<vector_type> > clone()
381 : {
382 0 : return make_sp(new ILU<algebra_type>(*this));
383 : }
384 :
385 : /// Destructor
386 0 : virtual ~ILU(){}
387 :
388 : /// returns if parallel solving is supported
389 0 : virtual bool supports_parallel() const {return true;}
390 :
391 : /// set factor for \f$ ILU_{\beta} \f$
392 0 : void set_beta(double beta) {m_beta = beta;}
393 :
394 : /// sets an ordering algorithm
395 0 : void set_ordering_algorithm(SmartPtr<ordering_algo_type> ordering_algo){
396 0 : m_spOrderingAlgo = ordering_algo;
397 0 : }
398 :
399 : /// set cuthill-mckee sort on/off
400 0 : void set_sort(bool b)
401 : {
402 0 : if(b){
403 0 : m_spOrderingAlgo = make_sp(new NativeCuthillMcKeeOrdering<TAlgebra, ordering_container_type>());
404 : }
405 : else{
406 0 : m_spOrderingAlgo = SPNULL;
407 : }
408 :
409 : UG_LOG("\nILU: please use 'set_ordering_algorithm(..)' in the future\n");
410 0 : }
411 :
412 : /// disable preprocessing (if underlying matrix has not changed)
413 0 : void set_disable_preprocessing(bool bDisable) {m_bDisablePreprocessing = bDisable;}
414 :
415 : /// sets the smallest allowed value for sorted factorization
416 0 : void set_sort_eps(number eps) {m_sortEps = eps;}
417 :
418 : /// sets the smallest allowed value for the Aii/Bi quotient
419 0 : void set_inversion_eps(number eps) {m_invEps = eps;}
420 :
421 : /// enables consistent interfaces.
422 : /** Connections between coefficients which lie in the same parallel interface
423 : * are made consistent between processes.*/
424 0 : void enable_consistent_interfaces (bool enable) {m_useConsistentInterfaces = enable;}
425 :
426 0 : void enable_overlap (bool enable) {m_useOverlap = enable;}
427 :
428 : protected:
429 : // Name of preconditioner
430 0 : virtual const char* name() const {return "ILU";}
431 :
432 0 : void apply_ordering()
433 : {
434 0 : if (!m_spOrderingAlgo.valid())
435 : return;
436 :
437 0 : if (m_useOverlap)
438 0 : UG_THROW ("ILU: Ordering for overlap has not been implemented yet.");
439 :
440 : #ifndef NDEBUG
441 : double start = get_clock_s();
442 : #endif
443 0 : if (m_u)
444 0 : m_spOrderingAlgo->init(&m_ILU, *m_u);
445 : else
446 0 : m_spOrderingAlgo->init(&m_ILU);
447 :
448 0 : m_spOrderingAlgo->compute();
449 : #ifndef NDEBUG
450 : double end = get_clock_s();
451 : UG_LOG("ILU: ordering took " << end-start << " seconds\n");
452 : #endif
453 0 : m_ordering = m_spOrderingAlgo->ordering();
454 :
455 0 : m_bSortIsIdentity = GetInversePermutation(m_ordering, m_old_ordering);
456 :
457 0 : if (!m_bSortIsIdentity)
458 : {
459 0 : matrix_type tmp;
460 0 : tmp = m_ILU;
461 0 : SetMatrixAsPermutation(m_ILU, tmp, m_ordering);
462 0 : }
463 : }
464 :
465 : protected:
466 0 : virtual bool init(SmartPtr<ILinearOperator<vector_type> > J,
467 : const vector_type& u)
468 : {
469 : // cast to matrix based operator
470 0 : SmartPtr<MatrixOperator<matrix_type, vector_type> > pOp =
471 : J.template cast_dynamic<MatrixOperator<matrix_type, vector_type> >();
472 :
473 : // Check that matrix if of correct type
474 0 : if(pOp.invalid())
475 0 : UG_THROW(name() << "::init': Passed Operator is "
476 : "not based on matrix. This Preconditioner can only "
477 : "handle matrix-based operators.");
478 :
479 0 : m_u = &u;
480 :
481 : // forward request to matrix based implementation
482 0 : return base_type::init(pOp);
483 : }
484 :
485 0 : bool init(SmartPtr<ILinearOperator<vector_type> > L)
486 : {
487 : // cast to matrix based operator
488 0 : SmartPtr<MatrixOperator<matrix_type, vector_type> > pOp =
489 : L.template cast_dynamic<MatrixOperator<matrix_type, vector_type> >();
490 :
491 : // Check that matrix if of correct type
492 0 : if(pOp.invalid())
493 0 : UG_THROW(name() << "::init': Passed Operator is "
494 : "not based on matrix. This Preconditioner can only "
495 : "handle matrix-based operators.");
496 :
497 0 : m_u = NULL;
498 :
499 : // forward request to matrix based implementation
500 0 : return base_type::init(pOp);
501 : }
502 :
503 : bool init(SmartPtr<MatrixOperator<matrix_type, vector_type> > Op)
504 : {
505 : m_u = NULL;
506 :
507 : return base_type::init(Op);
508 : }
509 :
510 : protected:
511 :
512 : // Preprocess routine
513 0 : virtual bool preprocess(SmartPtr<MatrixOperator<matrix_type, vector_type> > pOp)
514 : {
515 : // do not do a thing if preprocessing disabled
516 0 : if (m_bDisablePreprocessing) return true;
517 :
518 0 : matrix_type &mat = *pOp;
519 : PROFILE_BEGIN_GROUP(ILU_preprocess, "algebra ILU");
520 : // Debug output of matrices
521 : #ifdef UG_PARALLEL
522 : write_overlap_debug(mat, "ILU_prep_01_A_BeforeMakeUnique");
523 : #else
524 0 : write_debug(mat, "ILU_PreProcess_orig_A");
525 : #endif
526 :
527 0 : m_ILU = mat;
528 :
529 : #ifdef UG_PARALLEL
530 : if(m_useOverlap){
531 : CreateOverlap(m_ILU);
532 : m_oD.set_layouts(m_ILU.layouts());
533 : m_oC.set_layouts(m_ILU.layouts());
534 : m_oD.resize(m_ILU.num_rows(), false);
535 : m_oC.resize(m_ILU.num_rows(), false);
536 :
537 : if(debug_writer().valid()){
538 : m_overlapWriter = make_sp(new OverlapWriter<TAlgebra>());
539 : m_overlapWriter->init (*m_ILU.layouts(),
540 : *debug_writer(),
541 : m_ILU.num_rows());
542 : }
543 : }
544 : else if(m_useConsistentInterfaces){
545 : MatMakeConsistentOverlap0(m_ILU);
546 : }
547 : else {
548 : MatAddSlaveRowsToMasterRowOverlap0(m_ILU);
549 : // set dirichlet rows on slaves
550 : std::vector<IndexLayout::Element> vIndex;
551 : CollectUniqueElements(vIndex, m_ILU.layouts()->slave());
552 : SetDirichletRow(m_ILU, vIndex);
553 : }
554 : #endif
555 :
556 0 : m_h.resize(m_ILU.num_cols());
557 :
558 : #ifdef UG_PARALLEL
559 : write_overlap_debug(m_ILU, "ILU_prep_02_A_AfterMakeUnique");
560 : #endif
561 :
562 0 : apply_ordering();
563 :
564 : // Debug output of matrices
565 : #ifdef UG_PARALLEL
566 : write_overlap_debug(m_ILU, "ILU_prep_03_A_BeforeFactorize");
567 : #else
568 0 : write_debug(m_ILU, "ILU_PreProcess_U_BeforeFactor");
569 : #endif
570 :
571 :
572 : // Compute ILU Factorization
573 0 : if (m_beta!=0.0) FactorizeILUBeta(m_ILU, m_beta);
574 0 : else if(matrix_type::rows_sorted) FactorizeILUSorted(m_ILU, m_sortEps);
575 : else FactorizeILU(m_ILU);
576 0 : m_ILU.defragment();
577 :
578 : // Debug output of matrices
579 : #ifdef UG_PARALLEL
580 : write_overlap_debug(m_ILU, "ILU_prep_04_A_AfterFactorize");
581 : #else
582 0 : write_debug(m_ILU, "ILU_PreProcess_U_AfterFactor");
583 : #endif
584 :
585 : // we're done
586 0 : return true;
587 : }
588 :
589 :
590 0 : void applyLU(vector_type &c, const vector_type &d, vector_type &tmp)
591 : {
592 :
593 0 : if(m_spOrderingAlgo.invalid() || m_bSortIsIdentity)
594 : {
595 : // apply iterator: c = LU^{-1}*d
596 0 : if(! invert_L(m_ILU, tmp, d)) // h := L^-1 d
597 : print_debugger_message("ILU: There were issues at inverting L\n");
598 0 : if(! invert_U(m_ILU, c, tmp, m_invEps)) // c := U^-1 h = (LU)^-1 d
599 : print_debugger_message("ILU: There were issues at inverting U\n");
600 : }
601 : ///*
602 : else
603 : {
604 : // we save one vector here by renaming
605 0 : SetVectorAsPermutation(tmp, d, m_ordering);
606 0 : if(! invert_L(m_ILU, c, tmp)) // c = L^{-1} d
607 : print_debugger_message("ILU: There were issues at inverting L (after permutation)\n");
608 0 : if(! invert_U(m_ILU, tmp, c, m_invEps)) // tmp = (LU)^{-1} d
609 : print_debugger_message("ILU: There were issues at inverting U (after permutation)\n");
610 0 : SetVectorAsPermutation(c, tmp, m_old_ordering);
611 : }
612 : //*/
613 0 : }
614 :
615 : // Stepping routine
616 0 : virtual bool step(SmartPtr<MatrixOperator<matrix_type, vector_type> > pOp,
617 : vector_type& c,
618 : const vector_type& d)
619 : {
620 : PROFILE_BEGIN_GROUP(ILU_step, "algebra ILU");
621 :
622 :
623 : // \todo: introduce damping
624 : #ifdef UG_PARALLEL
625 : // for debug output (only for application is written)
626 : static bool first = true;
627 :
628 : if(first) write_overlap_debug(d, "ILU_step_1_d");
629 : if(m_useOverlap){
630 : for(size_t i = 0; i < d.size(); ++i)
631 : m_oD[i] = d[i];
632 : for(size_t i = d.size(); i < m_oD.size(); ++i)
633 : m_oD[i] = 0;
634 : m_oD.set_storage_type(PST_ADDITIVE);
635 : m_oD.change_storage_type(PST_CONSISTENT);
636 :
637 : if(first) write_overlap_debug(m_oD, "ILU_step_2_oD_consistent");
638 :
639 : applyLU(m_oC, m_oD, m_h);
640 :
641 : for(size_t i = 0; i < c.size(); ++i)
642 : c[i] = m_oC[i];
643 : SetLayoutValues(&c, c.layouts()->slave(), 0);
644 : c.set_storage_type(PST_UNIQUE);
645 : }
646 : else if(m_useConsistentInterfaces){
647 : // make defect consistent
648 : SmartPtr<vector_type> spDtmp = d.clone();
649 : spDtmp->change_storage_type(PST_CONSISTENT);
650 : applyLU(c, *spDtmp, m_h);
651 :
652 : // declare c unique to enforce that only master correction is used
653 : // when it is made consistent below
654 : c.set_storage_type(PST_UNIQUE);
655 : }
656 : else{
657 : // make defect unique
658 : SmartPtr<vector_type> spDtmp = d.clone();
659 : // SetVectorAsPermutation(*spDtmp, d, m_ordering);
660 :
661 : spDtmp->change_storage_type(PST_UNIQUE);
662 : if(first) write_debug(*spDtmp, "ILU_step_2_d_unique");
663 : applyLU(c, *spDtmp, m_h);
664 : c.set_storage_type(PST_ADDITIVE);
665 : }
666 :
667 : // write debug
668 : if(first) write_overlap_debug(c, "ILU_step_3_c");
669 :
670 : c.change_storage_type(PST_CONSISTENT);
671 :
672 : // write debug
673 : if(first) {write_overlap_debug(c, "ILU_step_4_c_consistent"); first = false;}
674 :
675 : #else
676 0 : write_debug(d, "ILU_step_d");
677 0 : applyLU(c, d, m_h);
678 0 : write_debug(c, "ILU_step_c");
679 : #endif
680 :
681 : // we're done
682 0 : return true;
683 : }
684 :
685 : /// Postprocess routine
686 0 : virtual bool postprocess() {return true;}
687 :
688 : private:
689 : #ifdef UG_PARALLEL
690 : template <class T> void write_overlap_debug(const T& t, std::string name)
691 : {
692 : if(debug_writer().valid()){
693 : if(m_useOverlap && m_overlapWriter.valid() && t.layouts()->overlap_enabled())
694 : m_overlapWriter->write(t, name);
695 : else
696 : write_debug(t, name.c_str());
697 : }
698 : }
699 : #endif
700 :
701 : protected:
702 : /// storage for factorization
703 : matrix_type m_ILU;
704 :
705 : /// help vector
706 : vector_type m_h;
707 :
708 : /// for overlaps only
709 : vector_type m_oD;
710 : vector_type m_oC;
711 : #ifdef UG_PARALLEL
712 : SmartPtr<OverlapWriter<TAlgebra> > m_overlapWriter;
713 : #endif
714 :
715 : /// factor for ILU-beta
716 : number m_beta;
717 :
718 : /// smallest allowed value for sorted factorization
719 : number m_sortEps;
720 :
721 : /// smallest allowed value for the Aii/Bi quotient
722 : number m_invEps;
723 :
724 : /// whether or not to disable preprocessing
725 : bool m_bDisablePreprocessing;
726 :
727 : bool m_useConsistentInterfaces;
728 : bool m_useOverlap;
729 :
730 : /// for ordering algorithms
731 : SmartPtr<ordering_algo_type> m_spOrderingAlgo;
732 : ordering_container_type m_ordering, m_old_ordering;
733 : std::vector<size_t> m_newIndex, m_oldIndex;
734 : bool m_bSortIsIdentity;
735 :
736 : const vector_type* m_u;
737 : };
738 :
739 : } // end namespace ug
740 :
741 : #endif
|
