Line data Source code
1 : /*
2 : * Copyright (c) 2014-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Martin Rupp
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 : #ifndef __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__GAUSS__
34 : #define __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__GAUSS__
35 :
36 : #include "lib_algebra/operator/interface/preconditioner.h"
37 : #include "lib_algebra/algebra_common/core_smoothers.h"
38 : #include "lib_algebra/algebra_common/sparsematrix_util.h"
39 : #ifdef UG_PARALLEL
40 : #include "lib_algebra/parallelization/parallelization.h"
41 : #include "lib_algebra/parallelization/parallel_matrix_overlap_impl.h"
42 : #endif
43 :
44 : #include "lib_algebra/algebra_common/sparsematrix_util.h"
45 : #include "lib_algebra/algebra_common/local_helper.h"
46 : #include "lib_algebra/small_algebra/small_matrix/print.h"
47 : #include "lib_algebra/operator/preconditioner/ilut.h"
48 : #include "common/debug_print.h"
49 : #include "common/progress.h"
50 : #include "lib_algebra/small_algebra/additional_math.h"
51 :
52 : #include "lib_algebra/common/graph/graph.h"
53 : #include "lib_algebra/algebra_common/sparse_vector.h"
54 :
55 :
56 :
57 : namespace ug{
58 :
59 : //////////////////// FROM AMG //////////////////////////////
60 :
61 : template<typename TMatrixType, typename TRowType>
62 0 : void CopyOffDiagEntries(const TMatrixType &A, size_t i, TRowType &row, bool enforceNew=false)
63 : {
64 : // copy matrix entries selected by rule
65 0 : for(typename TMatrixType::const_row_iterator connij = A.begin_row(i); connij != A.end_row(i); ++connij)
66 : {
67 : const size_t j=connij.index();
68 0 : row(connij.index()) = (i==j) ? 0.0 : connij.value();
69 : }
70 0 : }
71 :
72 : //! Adds 'strong negative connections' to graph.
73 : /*! Criterion: \f$ -a_{ij} \ge \epsilon \max_{a_{ik}<0} |a_{ik}| \f$ */
74 : class StrongNegativeConnectionsByBlockNorm
75 : {
76 : protected:
77 : double m_theta;
78 : public:
79 0 : StrongNegativeConnectionsByBlockNorm(double theta) : m_theta(theta){}
80 :
81 : template<typename TRowType>
82 0 : void find(const TRowType &Ci, size_t i, cgraph &graph)
83 : {
84 : double dmax = 0.0;
85 : typedef typename TRowType::const_iterator const_iterator;
86 :
87 0 : for(const_iterator conn = Ci.begin(); conn != Ci.end(); ++conn)
88 : {
89 0 : if(conn.index() == i) continue; // skip diag
90 0 : double d = BlockNorm(conn.value());
91 : // std::cout << d << std::endl;
92 0 : if(d > dmax) dmax = BlockNorm(conn.value());
93 : }
94 :
95 0 : for(const_iterator conn = Ci.begin(); conn != Ci.end(); ++conn)
96 : {
97 0 : if(conn.index() == i) continue; // skip diag
98 0 : if( BlockNorm(conn.value()) >= m_theta * dmax)
99 0 : graph.set_connection(i, conn.index());
100 : }
101 0 : }
102 : };
103 :
104 :
105 : template<typename matrix_type>
106 0 : void CreateStrongConnectionGraphForSystems(const matrix_type &A, cgraph &graph, double theta)
107 : {
108 : PROFILE_FUNC();
109 0 : graph.resize(A.num_rows());
110 :
111 0 : for(size_t i=0; i< A.num_rows(); i++)
112 : {
113 : // Skip isolated node
114 0 : if(A.is_isolated(i)) continue;
115 :
116 : // copy all off-diag entries
117 : SparseVector<typename matrix_type::value_type> Ri(A.num_cols());
118 0 : CopyOffDiagEntries(A, i, Ri);
119 :
120 : // Find strong connections
121 : StrongNegativeConnectionsByBlockNorm strong(theta);
122 0 : strong.find(Ri, i, graph);
123 : }
124 0 : }
125 :
126 :
127 : /////////////////////////////////
128 :
129 : /**
130 : * Calculates the gs correction by updating all unknowns in 'indices' at once
131 : * with the inverse of this block stored in AlocInv.
132 : * @param A a SparseMatrix
133 : * @param x solution to be GS-updated
134 : * @param b right hand side b
135 : * @param AlocInv cached inverse on a block/slice, calculated with GetSliceDenseInverse
136 : * @param indices the indices of the block
137 : * @param tmp a temporary vector
138 : */
139 : // we cache the inverse! computing the inverse "on the fly" seems to be a lot slower.
140 : template<typename TSparseMatrixType, typename TVectorType>
141 0 : void GetBlockGSCorrection(const TSparseMatrixType &A, TVectorType &x, TVectorType &b,
142 : DenseMatrix<VariableArray2<double> > &AlocInv,
143 : std::vector<size_t> &indices, DenseVector<VariableArray1<double> > &tmp, DenseVector<VariableArray1<double> > &tmp2)
144 : {
145 : // assume tmp is big enough
146 : size_t k;
147 :
148 : typedef typename TSparseMatrixType::const_row_iterator const_row_iterator;
149 : typedef typename TVectorType::value_type smallvec_type;
150 :
151 : k = 0;
152 0 : tmp.resize(indices.size()*GetSize(b[0]));
153 0 : tmp2.resize(indices.size()*GetSize(b[0]));
154 0 : for(size_t i=0; i<indices.size(); i++)
155 : {
156 0 : int j = indices[i];
157 0 : smallvec_type s = b[j];
158 : // calc b-Ax
159 0 : for(const_row_iterator it = A.begin_row(j); it != A.end_row(j); ++it)
160 0 : s -= it.value() * x[it.index()];
161 0 : for(size_t u=0; u<GetSize(s); u++)
162 0 : tmp[k++] = BlockRef(s, u);
163 : }
164 0 : MatMult(tmp2, 1.0, AlocInv, tmp);
165 :
166 : k=0;
167 0 : for(size_t i=0; i<indices.size(); i++)
168 : {
169 0 : smallvec_type &xi = x[indices[i]];
170 0 : for(size_t j=0; j<GetSize(xi); j++)
171 0 : BlockRef(xi, j) += tmp2[k++];
172 : }
173 0 : }
174 : template<typename TSparseMatrixType, typename TVectorType>
175 : void GetBlockGSCorrection(const TSparseMatrixType &A, TVectorType &x, TVectorType &b,
176 : DenseMatrix<VariableArray2<double> > &AlocInv,
177 : std::vector<size_t> &indices)
178 : {
179 : DenseVector<VariableArray1<double> > tmp;
180 : DenseVector<VariableArray1<double> > tmp2;
181 : GetBlockGSCorrection(A, x, b, AlocInv, indices, tmp, tmp2);
182 : }
183 :
184 : /**
185 : * Calculates the gs correction by updating all unknowns in 'indices' at once
186 : * with the inverse of this block implicitely in 'ilut'.
187 : * @param A a SparseMatrix
188 : * @param x solution to be GS-updated
189 : * @param b right hand side b
190 : * @param ilut cached ilut inverse type
191 : * @param indices the indices of the block
192 : * @param tmp a temporary vector
193 : * @param tmp2 another temporary vector
194 : */
195 : template<typename TSparseMatrixType, typename TVectorType>
196 0 : void GetBlockGSCorrectionILUT(const TSparseMatrixType &A, TVectorType &x, TVectorType &b,
197 : SmartPtr<ILUTPreconditioner<CPUAlgebra> > &ilut,
198 : std::vector<size_t> &indices, CPUAlgebra::vector_type &tmp, CPUAlgebra::vector_type &tmp2)
199 : {
200 : size_t blockSize = GetSize(b[0]);
201 : size_t k;
202 : typedef typename TSparseMatrixType::const_row_iterator const_row_iterator;
203 : typedef typename TVectorType::value_type smallvec_type;
204 :
205 : k = 0;
206 0 : tmp.resize(indices.size()*blockSize);
207 0 : tmp2.resize(indices.size()*blockSize);
208 0 : for(size_t i=0; i<indices.size(); i++)
209 : {
210 0 : int j = indices[i];
211 0 : smallvec_type s = b[j];
212 : // calc b-Ax
213 0 : for(const_row_iterator it = A.begin_row(j); it != A.end_row(j); ++it)
214 0 : s -= it.value() * x[it.index()];
215 0 : for(size_t u=0; u<GetSize(s); u++)
216 0 : tmp[k++] = BlockRef(s, u);
217 : }
218 0 : ilut->solve(tmp2, tmp);
219 :
220 : k=0;
221 0 : for(size_t i=0; i<indices.size(); i++)
222 : {
223 0 : smallvec_type &xi = x[indices[i]];
224 0 : for(size_t j=0; j<GetSize(xi); j++)
225 0 : BlockRef(xi, j) += tmp2[k++];
226 : }
227 0 : }
228 :
229 : /**
230 : * @param A a sparse matrix
231 : * @param indices map local -> global indices
232 : * @param AlocInv inverse on the indices to be used in @sa GetBlockGSCorrection
233 : */
234 : template<typename TSparseMatrixType>
235 0 : void GetSliceDenseInverse(const TSparseMatrixType &A, const std::vector<size_t> &indices,
236 : DenseMatrix<VariableArray2<double> > &AlocInv,
237 : DenseMatrix<VariableArray2<typename TSparseMatrixType::value_type> > &tmp)
238 : {
239 0 : tmp.resize(indices.size(), indices.size());
240 : GetLocalMatrix(A, tmp, &indices[0], &indices[0]);
241 0 : BlockMatrixToDoubleMatrix(AlocInv, tmp);
242 0 : Invert(AlocInv);
243 0 : }
244 :
245 : template<typename TSparseMatrixType>
246 : void GetSliceDenseInverse(const TSparseMatrixType &A, const std::vector<size_t> &indices,
247 : DenseMatrix<VariableArray2<double> > &AlocInv)
248 : {
249 : DenseMatrix<VariableArray2<typename TSparseMatrixType::value_type> > L;
250 : GetSliceDenseInverse(A, indices, AlocInv, L);
251 : }
252 :
253 : /**
254 : * @param A a sparse matrix
255 : * @param indices map local -> global indices
256 : * @param R sparse matrix slice in local indices
257 : */
258 : template<typename TSparseMatrixType>
259 0 : void GetSliceSparse(const TSparseMatrixType &A, const std::vector<size_t> &indices,
260 : CPUAlgebra::matrix_type &R)
261 : {
262 0 : size_t blockSize = GetRows(A(0,0));
263 : size_t numRows = indices.size();
264 : size_t numCols = indices.size();
265 :
266 0 : R.resize_and_clear(numRows*blockSize, numCols*blockSize);
267 0 : for(size_t ri=0; ri<indices.size(); ri++)
268 : {
269 0 : size_t r=indices[ri];
270 0 : for(size_t ci=0; ci<indices.size(); ci++)
271 : {
272 0 : size_t c = indices[ci];
273 0 : if(A.has_connection(r, c))
274 : {
275 0 : const typename TSparseMatrixType::value_type &m = A(r, c);
276 0 : for(size_t j1=0; j1<blockSize; j1++)
277 0 : for(size_t j2=0; j2<blockSize; j2++)
278 0 : R(ri*blockSize + j1, ci*blockSize + j2) = BlockRef(m, j1, j2);
279 : }
280 : }
281 : }
282 0 : R.defragment();
283 0 : }
284 :
285 :
286 : template<typename TAlgebra>
287 : class IBlockJacobiPreconditioner : public IPreconditioner<TAlgebra>
288 : {
289 : public:
290 : // Algebra type
291 : typedef TAlgebra algebra_type;
292 :
293 : // Vector type
294 : typedef typename TAlgebra::vector_type vector_type;
295 :
296 : // Matrix type
297 : typedef typename TAlgebra::matrix_type matrix_type;
298 :
299 : /// Matrix Operator type
300 : typedef typename IPreconditioner<TAlgebra>::matrix_operator_type matrix_operator_type;
301 :
302 0 : IBlockJacobiPreconditioner() {}
303 0 : IBlockJacobiPreconditioner(const IBlockJacobiPreconditioner &parent) : IPreconditioner<TAlgebra>(parent)
304 : { }
305 :
306 : protected:
307 : #ifdef UG_PARALLEL
308 : matrix_type A;
309 : #endif
310 0 : virtual bool preprocess(SmartPtr<MatrixOperator<matrix_type, vector_type> > pOp)
311 : {
312 :
313 0 : matrix_type &mat = *pOp;
314 :
315 : #ifdef UG_PARALLEL
316 : A = mat;
317 : MatAddSlaveRowsToMasterRowOverlap0(A);
318 :
319 : // set zero on slaves
320 : std::vector<IndexLayout::Element> vIndex;
321 : CollectUniqueElements(vIndex, mat.layouts()->slave());
322 : SetDirichletRow(A, vIndex);
323 : return block_preprocess(A);
324 : #else
325 0 : return block_preprocess(mat);
326 : #endif
327 : }
328 :
329 : // Preprocess routine
330 : virtual bool block_preprocess(matrix_type &A) = 0;
331 :
332 : // Postprocess routine
333 0 : virtual bool postprocess() {return true;}
334 0 : virtual bool supports_parallel() const { return true; }
335 :
336 : SmartPtr<vector_type> m_spDtmp;
337 :
338 : // Stepping routine
339 0 : virtual bool step(SmartPtr<MatrixOperator<matrix_type, vector_type> > pOp, vector_type& c, const vector_type& d)
340 : {
341 : #ifdef UG_PARALLEL
342 : if(!m_spDtmp.valid() || m_spDtmp->size() != d.size())
343 : m_spDtmp = d.clone();
344 : else
345 : (*m_spDtmp) = d;
346 : m_spDtmp->change_storage_type(PST_UNIQUE);
347 : bool b = block_step(A, c, *m_spDtmp);
348 : c.set_storage_type(PST_ADDITIVE);
349 : c.change_storage_type(PST_CONSISTENT);
350 : return b;
351 : #else
352 0 : return block_step(*pOp, c, d);
353 : #endif
354 : return true;
355 : }
356 :
357 : // Stepping routine
358 : virtual bool block_step(matrix_type &A, vector_type& c, const vector_type& d) = 0;
359 : };
360 :
361 : /**
362 : * BlockGaussSeidel
363 : * use the constructor or set_depth to set the depth
364 : * depth = 0 -> GS
365 : * depth = 1 -> Block i and neighbors of i
366 : * depth = 2 -> Block i and neighbors of i and their neighbors
367 : * depth = 3 ...
368 : */
369 : template <typename TAlgebra, bool backward, bool forward>
370 : class BlockGaussSeidel : public IBlockJacobiPreconditioner<TAlgebra>
371 : {
372 : public:
373 : // Algebra type
374 : typedef TAlgebra algebra_type;
375 :
376 : // Vector type
377 : typedef typename TAlgebra::vector_type vector_type;
378 :
379 : // Matrix type
380 : typedef typename TAlgebra::matrix_type matrix_type;
381 :
382 : /// Matrix Operator type
383 : typedef typename IPreconditioner<TAlgebra>::matrix_operator_type matrix_operator_type;
384 :
385 : /// Base type
386 : typedef IBlockJacobiPreconditioner<TAlgebra> base_type;
387 :
388 : typedef BlockGaussSeidel<algebra_type, backward, forward> this_type;
389 :
390 : protected:
391 : typedef typename matrix_type::value_type block_type;
392 :
393 : public:
394 : // Constructor
395 0 : BlockGaussSeidel() {
396 0 : m_depth = 1;
397 0 : };
398 :
399 0 : BlockGaussSeidel(int depth) {
400 0 : m_depth = depth;
401 0 : };
402 :
403 0 : BlockGaussSeidel(const this_type &parent) : base_type(parent)
404 : {
405 0 : m_depth = parent.m_depth;
406 : }
407 :
408 : // Clone
409 0 : virtual SmartPtr<ILinearIterator<vector_type> > clone()
410 : {
411 0 : return make_sp(new this_type(*this));
412 : }
413 :
414 : typedef typename matrix_type::value_type smallmat_type;
415 : typedef typename vector_type::value_type smallvec_type;
416 : typedef typename matrix_type::const_row_iterator const_row_iterator;
417 :
418 : std::vector< DenseMatrix<VariableArray2<double> > > AlocInv;
419 : size_t m_depth;
420 : std::vector<std::vector<size_t> > indices;
421 :
422 :
423 0 : void set_depth(size_t d)
424 : {
425 0 : m_depth = d;
426 0 : }
427 :
428 : protected:
429 : // Name of preconditioner
430 0 : virtual const char* name() const {return "BlockGaussSeidel";}
431 :
432 : // Preprocess routine
433 0 : virtual bool block_preprocess(matrix_type &A)
434 : {
435 : size_t N = A.num_rows();
436 : DenseMatrix<VariableArray2<smallmat_type> > tmpMat;
437 :
438 0 : AlocInv.clear();
439 0 : AlocInv.resize(N);
440 :
441 0 : indices.resize(N);
442 :
443 : size_t maxSize = 0;
444 0 : PROGRESS_START(prog, N, "BlockGaussSeidel: compute blocks");
445 0 : for(size_t i=0; i<N; i++)
446 : {
447 : PROGRESS_UPDATE(prog, i);
448 0 : std::vector<bool> bVisited(N, false);
449 : indices[i].clear();
450 0 : GetNeighborhood(A, i, m_depth, indices[i], bVisited);
451 : //indices[i].push_back(i);
452 :
453 0 : GetSliceDenseInverse(A, indices[i], AlocInv[i], tmpMat);
454 :
455 : //PrintVector(indices[i], "indices");
456 : // UG_LOG("AlocInv " << i << ":\n" << JuliaString(AlocInv[i], "Aloc") << "\n");
457 : // UG_LOG("AlocInv " << i << ":\n" << JuliaString(AlocInv[i], "AlocInv") << "\n");
458 :
459 0 : if(AlocInv[i].num_rows() > maxSize) maxSize = AlocInv[i].num_rows();
460 : }
461 0 : PROGRESS_FINISH(prog);
462 :
463 : UG_LOG("Max Size = " << maxSize << "\n");
464 0 : return true;
465 0 : }
466 :
467 : typedef typename matrix_type::const_row_iterator matrix_const_row_iterator;
468 : typedef typename matrix_type::row_iterator matrix_row_iterator;
469 :
470 : // Postprocess routine
471 0 : virtual bool postprocess() {return true;}
472 0 : virtual bool supports_parallel() const { return true; }
473 :
474 : // Stepping routine
475 0 : virtual bool block_step(matrix_type &A, vector_type& c, const vector_type& d)
476 : {
477 : PROFILE_BEGIN(BlockGaussSeidel_step);
478 : vector_type &x = c;
479 : x.set(0.0);
480 : vector_type b;
481 0 : b = d;
482 :
483 : DenseVector<VariableArray1<double> > tmp;
484 : DenseVector<VariableArray1<double> > tmp2;
485 :
486 : if(forward)
487 0 : for(size_t i=0; i<x.size(); i++)
488 : {
489 : // c = D^{-1}(b-Ax)
490 : // x = x + c
491 0 : GetBlockGSCorrection(A, x, b, AlocInv[i], indices[i], tmp, tmp2);
492 : // do_correction_implicit(A, x, b, indices[i]);
493 : }
494 : if(backward)
495 0 : for(size_t i=x.size()-1; ; i--)
496 : {
497 : // c = D^{-1}(b-Ax)
498 : // x = x + c
499 0 : GetBlockGSCorrection(A, x, b, AlocInv[i], indices[i], tmp, tmp2);
500 0 : if(i==0) break;
501 : // do_correction_implicit(A, x, b, indices[i]);
502 : }
503 :
504 : // Correction is always consistent
505 : #ifdef UG_PARALLEL
506 : c.set_storage_type(PST_CONSISTENT);
507 : #endif
508 0 : return true;
509 : }
510 :
511 0 : virtual std::string config_string() const
512 : {
513 0 : std::stringstream ss ; ss << "BlockGaussSeidel(depth = " << m_depth << ")";
514 0 : return ss.str();
515 0 : }
516 :
517 : };
518 :
519 :
520 : template <typename TAlgebra, bool backward, bool forward>
521 : class BlockGaussSeidelIterative : public IBlockJacobiPreconditioner<TAlgebra>
522 : {
523 : public:
524 : // Algebra type
525 : typedef TAlgebra algebra_type;
526 :
527 : // Vector type
528 : typedef typename TAlgebra::vector_type vector_type;
529 :
530 : // Matrix type
531 : typedef typename TAlgebra::matrix_type matrix_type;
532 :
533 : /// Matrix Operator type
534 : typedef typename IPreconditioner<TAlgebra>::matrix_operator_type matrix_operator_type;
535 :
536 : /// Base type
537 : typedef IBlockJacobiPreconditioner<TAlgebra> base_type;
538 :
539 : typedef BlockGaussSeidelIterative<algebra_type, backward, forward> this_type;
540 :
541 : protected:
542 : typedef typename matrix_type::value_type block_type;
543 : public:
544 : // Constructor
545 0 : BlockGaussSeidelIterative() {
546 0 : m_depth = 1;
547 0 : m_nu = 2;
548 0 : };
549 :
550 0 : BlockGaussSeidelIterative(int depth, int nu) {
551 0 : m_depth = depth;
552 0 : m_nu = nu;
553 0 : };
554 :
555 0 : BlockGaussSeidelIterative(const this_type &parent) : base_type(parent)
556 : {
557 0 : m_depth = parent.m_depth;
558 0 : m_nu = parent.m_nu;
559 : }
560 :
561 : // Clone
562 0 : virtual SmartPtr<ILinearIterator<vector_type> > clone()
563 : {
564 0 : return make_sp(new this_type(*this));
565 : }
566 :
567 : typedef typename matrix_type::value_type smallmat_type;
568 : typedef typename vector_type::value_type smallvec_type;
569 : typedef typename matrix_type::const_row_iterator const_row_iterator;
570 :
571 : std::vector< DenseMatrix<VariableArray2<double> > > AlocInv;
572 : size_t m_depth;
573 : std::vector<std::vector<size_t> > indices;
574 :
575 :
576 0 : void set_depth(size_t d)
577 : {
578 0 : m_depth = d;
579 0 : }
580 :
581 :
582 : protected:
583 : // Name of preconditioner
584 0 : virtual const char* name() const
585 : {
586 : if(backward&&forward) return "SymmetricBlockGaussSeidelIterative";
587 : else if(backward) return "BackwardBlockGaussSeidelIterative";
588 : return "BlockGaussSeidelIterative";
589 : }
590 :
591 :
592 :
593 : // Preprocess routine
594 0 : virtual bool block_preprocess(matrix_type &A)
595 : {
596 : size_t N = A.num_rows();
597 0 : indices.resize(N);
598 :
599 0 : size_t maxSize = 0;
600 0 : for(size_t i=0; i<N; i++)
601 : {
602 0 : std::vector<bool> bVisited(N, false);
603 : indices[i].clear();
604 0 : GetNeighborhood(A, i, m_depth, indices[i], bVisited);
605 0 : maxSize = std::max(indices[i].size(), maxSize);
606 : }
607 0 : UG_LOG("Max Size = " << maxSize << "\n");
608 0 : return true;
609 : }
610 :
611 : typedef typename matrix_type::const_row_iterator matrix_const_row_iterator;
612 : typedef typename matrix_type::row_iterator matrix_row_iterator;
613 :
614 : // Postprocess routine
615 0 : virtual bool postprocess() {return true;}
616 0 : virtual bool supports_parallel() const { return true; }
617 :
618 :
619 0 : void correct(size_t i, const matrix_type &A, vector_type& x, const vector_type& b)
620 : {
621 :
622 0 : smallvec_type s = b[i];
623 : // calc b-Ax
624 0 : for(const_row_iterator it = A.begin_row(i); it != A.end_row(i); ++it)
625 0 : s -= it.value() * x[it.index()];
626 : smallvec_type c;
627 0 : InverseMatMult(c, 1.0, A(i,i), s);
628 0 : x[i] += c;
629 :
630 0 : }
631 :
632 0 : void correct_forward(size_t i, matrix_type &A, vector_type& x, const vector_type& b)
633 : {
634 0 : for(size_t k=0; k<m_nu; k++)
635 0 : for(size_t j=0; j<indices[i].size(); j++)
636 0 : correct(indices[i][j], A, x, b);
637 0 : }
638 :
639 0 : void correct_backward(size_t i, matrix_type &A, vector_type& x, const vector_type& b)
640 : {
641 0 : for(size_t k=0; k<m_nu; k++)
642 0 : for(int j=(int)(indices[i].size())-1; j>=0 ; j--)
643 0 : correct(indices[i][j], A, x, b);
644 0 : }
645 :
646 : size_t m_nu;
647 :
648 : // Stepping routine
649 0 : virtual bool block_step(matrix_type &A, vector_type& c, const vector_type& d)
650 : {
651 : PROFILE_BEGIN(BlockGaussSeidelIterative_step);
652 : vector_type &x = c;
653 : x.set(0.0);
654 : vector_type b;
655 0 : b = d;
656 :
657 : if(forward)
658 0 : for(size_t i=0; i<x.size(); i++)
659 : {
660 0 : correct_forward(i, A, x, b);
661 : }
662 : if(backward)
663 0 : for(size_t i=x.size()-1; ; i--)
664 : {
665 0 : correct_backward(i, A, x, b);
666 0 : if(i==0) break;
667 : }
668 :
669 : // Correction is always consistent
670 : #ifdef UG_PARALLEL
671 : c.set_storage_type(PST_CONSISTENT);
672 : #endif
673 0 : return true;
674 : }
675 :
676 : public:
677 0 : void set_iterative_steps(size_t nu)
678 : {
679 0 : m_nu=nu;
680 0 : }
681 :
682 :
683 0 : virtual std::string config_string() const
684 : {
685 0 : std::stringstream ss ;
686 0 : if(backward&&forward) ss << "Symmetric";
687 0 : else if(backward) ss << "Backward";
688 0 : ss << "BlockGaussSeidelIterative(depth = " << m_depth << ", nu = " << m_nu << ")";
689 0 : return ss.str();
690 0 : }
691 :
692 : };
693 :
694 :
695 : /**
696 : * SparseBlockGaussSeidel
697 : * experimental version
698 : * a) can use bigger stencils since it uses SparseLU for solving blocks
699 : * b) tries to use some overlapping blocks (BlockGaussSeidel always uses N blocks)
700 : */
701 : template <typename TAlgebra, bool backward, bool forward>
702 : class SparseBlockGaussSeidel : public IBlockJacobiPreconditioner<TAlgebra>
703 : {
704 : public:
705 : // Algebra type
706 : typedef TAlgebra algebra_type;
707 :
708 : // Vector type
709 : typedef typename TAlgebra::vector_type vector_type;
710 :
711 : // Matrix type
712 : typedef typename TAlgebra::matrix_type matrix_type;
713 :
714 : /// Matrix Operator type
715 : typedef typename IPreconditioner<TAlgebra>::matrix_operator_type matrix_operator_type;
716 :
717 : /// Base type
718 : typedef IBlockJacobiPreconditioner<TAlgebra> base_type;
719 :
720 : protected:
721 : typedef typename matrix_type::value_type block_type;
722 : typedef SparseBlockGaussSeidel<TAlgebra, backward, forward> this_type;
723 :
724 : public:
725 : // Constructor
726 0 : SparseBlockGaussSeidel() {
727 0 : m_depth = 1;
728 0 : };
729 :
730 0 : SparseBlockGaussSeidel(int depth) {
731 0 : m_depth = depth;
732 0 : };
733 :
734 0 : SparseBlockGaussSeidel(const this_type &parent) : base_type(parent)
735 : {
736 0 : m_depth = parent.m_depth;
737 0 : }
738 :
739 : // Clone
740 0 : virtual SmartPtr<ILinearIterator<vector_type> > clone()
741 : {
742 0 : return make_sp(new this_type(*this));
743 : }
744 :
745 :
746 : typedef typename matrix_type::value_type smallmat_type;
747 : typedef typename vector_type::value_type smallvec_type;
748 : typedef typename matrix_type::const_row_iterator const_row_iterator;
749 :
750 : std::map<size_t, SmartPtr<ILUTPreconditioner<CPUAlgebra> > > m_ilut;
751 :
752 : size_t m_depth;
753 : std::vector<std::vector<size_t> > indices;
754 : std::map<size_t, SmartPtr<CPUAlgebra::matrix_type> > Aloc;
755 :
756 :
757 0 : void set_depth(size_t d)
758 : {
759 0 : m_depth = d;
760 0 : }
761 :
762 :
763 : protected:
764 : // Name of preconditioner
765 0 : virtual const char* name() const {return "SparseBlockGaussSeidel";}
766 :
767 : // Preprocess routine
768 0 : virtual bool block_preprocess(matrix_type &A)
769 : {
770 : size_t N = A.num_rows();
771 : DenseMatrix<VariableArray2<smallmat_type> > tmpMat;
772 :
773 :
774 : m_ilut.clear();
775 :
776 0 : indices.resize(N);
777 :
778 : size_t maxSize = 0;
779 0 : PROGRESS_START(prog, N, "SparseBlockGaussSeidel: compute blocks");
780 :
781 :
782 :
783 0 : std::vector<size_t> bVisited(N, 0);
784 0 : std::vector<bool> bVisited2(N, false);
785 :
786 : std::vector<size_t> levels;
787 :
788 0 : for(size_t i=0; i<N; i++)
789 : {
790 : // if(bVisited[i]>5) continue;
791 : PROGRESS_UPDATE(prog, i);
792 :
793 0 : indices[i].clear();
794 :
795 0 : GetNeighborhood(A, i, m_depth, indices[i], bVisited2);
796 0 : for(size_t j=0; j<indices[i].size(); j++)
797 0 : bVisited[indices[i][j]] ++;
798 : //for(size_t j=0; j<levels[m_depth-1]; j++)
799 : //bVisited2[indices[j]]=true;
800 :
801 0 : Aloc[i] = make_sp(new CPUAlgebra::matrix_type);
802 :
803 0 : GetSliceSparse(A, indices[i], *Aloc[i]);
804 :
805 :
806 0 : m_ilut[i] = make_sp(new ILUTPreconditioner<CPUAlgebra> (0.0));
807 0 : m_ilut[i]->preprocess_mat(*Aloc[i]);
808 :
809 :
810 :
811 0 : if(Aloc[i]->num_rows() > maxSize) maxSize = Aloc[i]->num_rows();
812 : }
813 0 : PROGRESS_FINISH(prog);
814 :
815 : UG_LOG("Max Size = " << maxSize << "\n");
816 0 : return true;
817 0 : }
818 :
819 : typedef typename matrix_type::const_row_iterator matrix_const_row_iterator;
820 : typedef typename matrix_type::row_iterator matrix_row_iterator;
821 :
822 : // Postprocess routine
823 0 : virtual bool postprocess() {return true;}
824 0 : virtual bool supports_parallel() const { return true; }
825 :
826 : // Stepping routine
827 0 : virtual bool block_step(matrix_type &A, vector_type& c, const vector_type& d)
828 : {
829 : vector_type &x = c;
830 : x.set(0.0);
831 : vector_type b;
832 0 : b = d;
833 :
834 : DenseMatrix<VariableArray2<double> > Adense;
835 : DenseMatrix<VariableArray2<smallmat_type> > Atmp;
836 : CPUAlgebra::vector_type tmp, tmp2;
837 0 : PROGRESS_START(prog, x.size(), "SparseBlockGaussSeidel: step");
838 : if(forward)
839 0 : for(size_t i=0; i<x.size(); i++)
840 : {
841 : PROGRESS_UPDATE(prog, i);
842 : // c = D^{-1}(b-Ax)
843 : // x = x + c
844 0 : if(indices[i].size() != 0)
845 0 : GetBlockGSCorrectionILUT(A, x, b, m_ilut[i], indices[i], tmp, tmp2);
846 : }
847 : if(backward)
848 0 : for(size_t i=x.size()-1; ; i--)
849 : {
850 0 : PROGRESS_UPDATE(prog, i);
851 : // c = D^{-1}(b-Ax)
852 : // x = x + c
853 0 : if(indices[i].size() != 0)
854 0 : GetBlockGSCorrectionILUT(A, x, b, m_ilut[i], indices[i], tmp, tmp2);
855 0 : if(i==0) break;
856 : }
857 0 : PROGRESS_FINISH(prog);
858 :
859 : // Correction is always consistent
860 : #ifdef UG_PARALLEL
861 : c.set_storage_type(PST_CONSISTENT);
862 : #endif
863 0 : return true;
864 0 : }
865 :
866 0 : virtual std::string config_string() const
867 : {
868 0 : std::stringstream ss ;
869 0 : if(backward&&forward) ss << "Symmetric";
870 0 : else if(backward) ss << "Backward";
871 0 : ss << "SparseBlockGaussSeidel(depth = " << m_depth << ")";
872 0 : return ss.str();
873 0 : }
874 :
875 : };
876 :
877 :
878 : /**
879 : * SparseBlockGaussSeidel
880 : * experimental version
881 : * a) can use bigger stencils since it uses SparseLU for solving blocks
882 : * b) tries to use some overlapping blocks (BlockGaussSeidel always uses N blocks)
883 : */
884 : template <typename TAlgebra, bool backward, bool forward>
885 : class SparseBlockGaussSeidel2 : public IBlockJacobiPreconditioner<TAlgebra>
886 : {
887 : public:
888 : // Algebra type
889 : typedef TAlgebra algebra_type;
890 :
891 : // Vector type
892 : typedef typename TAlgebra::vector_type vector_type;
893 :
894 : // Matrix type
895 : typedef typename TAlgebra::matrix_type matrix_type;
896 :
897 : /// Matrix Operator type
898 : typedef typename IPreconditioner<TAlgebra>::matrix_operator_type matrix_operator_type;
899 :
900 : /// Base type
901 : typedef IBlockJacobiPreconditioner<TAlgebra> base_type;
902 : typedef SparseBlockGaussSeidel2<TAlgebra, backward, forward> this_type;
903 :
904 : protected:
905 : typedef typename matrix_type::value_type block_type;
906 : public:
907 : // Constructor
908 0 : SparseBlockGaussSeidel2() {
909 0 : m_depth = 1;
910 0 : };
911 :
912 0 : SparseBlockGaussSeidel2(int depth) {
913 0 : m_depth = depth;
914 0 : };
915 :
916 0 : SparseBlockGaussSeidel2(const this_type &parent) : base_type(parent)
917 : {
918 0 : m_depth = parent.m_depth;
919 0 : }
920 :
921 : // Clone
922 0 : virtual SmartPtr<ILinearIterator<vector_type> > clone()
923 : {
924 0 : return make_sp(new this_type(*this));
925 : }
926 :
927 : typedef typename matrix_type::value_type smallmat_type;
928 : typedef typename vector_type::value_type smallvec_type;
929 : typedef typename matrix_type::const_row_iterator const_row_iterator;
930 :
931 : std::map<size_t, SmartPtr<ILUTPreconditioner<CPUAlgebra> > > m_ilut;
932 :
933 : size_t m_depth;
934 : std::vector<std::vector<size_t> > indices;
935 : std::map<size_t, SmartPtr<CPUAlgebra::matrix_type> > Aloc;
936 :
937 :
938 0 : void set_depth(size_t d)
939 : {
940 0 : m_depth = d;
941 0 : }
942 :
943 :
944 : protected:
945 : // Name of preconditioner
946 0 : virtual const char* name() const {return "SparseBlockGaussSeidel2";}
947 :
948 : // Preprocess routine
949 0 : virtual bool block_preprocess(matrix_type &A)
950 : {
951 : size_t N = A.num_rows();
952 : DenseMatrix<VariableArray2<smallmat_type> > tmpMat;
953 : m_ilut.clear();
954 :
955 0 : indices.resize(N);
956 :
957 : size_t maxSize = 0;
958 0 : PROGRESS_START(prog, N, "SparseBlockGaussSeidel: compute blocks");
959 :
960 : cgraph G;
961 : {
962 : cgraph graph;
963 0 : CreateStrongConnectionGraphForSystems(A, graph, 0.3);
964 0 : G.resize(N);
965 0 : for(size_t i=0; i<graph.size(); i++)
966 0 : for(cgraph::row_iterator it = graph.begin_row(i); it != graph.end_row(i); ++it)
967 : {
968 0 : G.set_connection(*it, i);
969 0 : G.set_connection(i, *it);
970 : }
971 : }
972 0 : std::vector<int> iComponent(N, -1);
973 : std::vector<std::vector<int> > components;
974 : // G.print();
975 0 : for(size_t i=0; i<G.size(); i++)
976 : {
977 0 : if(iComponent[i] != -1) continue;
978 0 : int myComponent = iComponent[i] = components.size();
979 0 : components.resize(components.size()+1);
980 0 : std::vector<int> &myComponents = components[myComponent];
981 0 : myComponents.push_back(i);
982 0 : for(size_t c=0; c<myComponents.size(); c++)
983 : {
984 0 : int j = myComponents[c];
985 0 : for(cgraph::row_iterator it = G.begin_row(j); it != G.end_row(j); ++it)
986 : {
987 0 : if(iComponent[*it] == myComponent) continue;
988 0 : myComponents.push_back(*it);
989 0 : UG_COND_THROW(iComponent[*it] != -1, i );
990 0 : iComponent[*it] = myComponent;
991 : }
992 : }
993 : }
994 0 : for(size_t i=0; i<N; i++)
995 : indices[i].clear();
996 0 : for(size_t c=0; c<components.size(); c++)
997 : {
998 0 : int i=components[c][0];
999 0 : indices[i].insert(indices[i].begin(), components[c].begin(), components[c].end());
1000 0 : PRINT_VECTOR(indices[i], i);
1001 :
1002 0 : Aloc[i] = make_sp(new CPUAlgebra::matrix_type);
1003 0 : GetSliceSparse(A, indices[i], *Aloc[i]);
1004 :
1005 0 : m_ilut[i] = make_sp(new ILUTPreconditioner<CPUAlgebra> (0.0));
1006 0 : m_ilut[i]->preprocess_mat(*Aloc[i]);
1007 0 : if(Aloc[i]->num_rows() > maxSize) maxSize = Aloc[i]->num_rows();
1008 : }
1009 0 : PROGRESS_FINISH(prog);
1010 :
1011 : UG_LOG("Max Size = " << maxSize << "\n");
1012 0 : return true;
1013 0 : }
1014 :
1015 : typedef typename matrix_type::const_row_iterator matrix_const_row_iterator;
1016 : typedef typename matrix_type::row_iterator matrix_row_iterator;
1017 :
1018 : // Postprocess routine
1019 0 : virtual bool postprocess() {return true;}
1020 0 : virtual bool supports_parallel() const { return true; }
1021 :
1022 : // Stepping routine
1023 0 : virtual bool block_step(matrix_type &A, vector_type& c, const vector_type& d)
1024 : {
1025 : vector_type &x = c;
1026 : x.set(0.0);
1027 : vector_type b;
1028 0 : b = d;
1029 :
1030 : DenseMatrix<VariableArray2<double> > Adense;
1031 : DenseMatrix<VariableArray2<smallmat_type> > Atmp;
1032 : CPUAlgebra::vector_type tmp, tmp2;
1033 0 : PROGRESS_START(prog, x.size(), "SparseBlockGaussSeidel: step");
1034 : if(forward)
1035 : for(size_t i=0; i<x.size(); i++)
1036 : {
1037 : PROGRESS_UPDATE(prog, i);
1038 : // c = D^{-1}(b-Ax)
1039 : // x = x + c
1040 : if(indices[i].size() != 0)
1041 : GetBlockGSCorrectionILUT(A, x, b, m_ilut[i], indices[i], tmp, tmp2);
1042 : }
1043 : if(backward)
1044 0 : for(size_t i=x.size()-1; ; i--)
1045 : {
1046 0 : PROGRESS_UPDATE(prog, i);
1047 : // c = D^{-1}(b-Ax)
1048 : // x = x + c
1049 0 : if(indices[i].size() != 0)
1050 0 : GetBlockGSCorrectionILUT(A, x, b, m_ilut[i], indices[i], tmp, tmp2);
1051 0 : if(i==0) break;
1052 : }
1053 0 : PROGRESS_FINISH(prog);
1054 :
1055 : // Correction is always consistent
1056 : #ifdef UG_PARALLEL
1057 : c.set_storage_type(PST_CONSISTENT);
1058 : #endif
1059 0 : return true;
1060 0 : }
1061 :
1062 0 : virtual std::string config_string() const
1063 : {
1064 0 : std::stringstream ss ;
1065 : if(backward&&forward) ss << "Symmetric";
1066 0 : else if(backward) ss << "Backward";
1067 0 : ss << "SparseBlockGaussSeidel(depth = " << m_depth << ")";
1068 0 : return ss.str();
1069 0 : }
1070 :
1071 : };
1072 :
1073 :
1074 : } // end namespace ug
1075 :
1076 : #endif // __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__GAUSS_SEIDEL__
|
