ug4
Loading...
Searching...
No Matches
block_gauss_seidel.h
Go to the documentation of this file.
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
39#ifdef UG_PARALLEL
42#endif
43
48#include "common/debug_print.h"
49#include "common/progress.h"
51
54
55
56
57namespace ug{
58
60
61template<typename TMatrixType, typename TRowType>
62void CopyOffDiagEntries(const TMatrixType &A, size_t i, TRowType &row, bool enforceNew=false)
63{
64 // copy matrix entries selected by rule
65 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 row(connij.index()) = (i==j) ? 0.0 : connij.value();
69 }
70}
71
73
75{
76protected:
77 double m_theta;
78public:
80
81 template<typename TRowType>
82 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 for(const_iterator conn = Ci.begin(); conn != Ci.end(); ++conn)
88 {
89 if(conn.index() == i) continue; // skip diag
90 double d = BlockNorm(conn.value());
91// std::cout << d << std::endl;
92 if(d > dmax) dmax = BlockNorm(conn.value());
93 }
94
95 for(const_iterator conn = Ci.begin(); conn != Ci.end(); ++conn)
96 {
97 if(conn.index() == i) continue; // skip diag
98 if( BlockNorm(conn.value()) >= m_theta * dmax)
99 graph.set_connection(i, conn.index());
100 }
101 }
102};
103
104
105template<typename matrix_type>
106void CreateStrongConnectionGraphForSystems(const matrix_type &A, cgraph &graph, double theta)
107{
108 PROFILE_FUNC();
109 graph.resize(A.num_rows());
110
111 for(size_t i=0; i< A.num_rows(); i++)
112 {
113 // Skip isolated node
114 if(A.is_isolated(i)) continue;
115
116 // copy all off-diag entries
118 CopyOffDiagEntries(A, i, Ri);
119
120 // Find strong connections
122 strong.find(Ri, i, graph);
123 }
124}
125
126
128
139// we cache the inverse! computing the inverse "on the fly" seems to be a lot slower.
140template<typename TSparseMatrixType, typename TVectorType>
141void GetBlockGSCorrection(const TSparseMatrixType &A, TVectorType &x, TVectorType &b,
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 tmp.resize(indices.size()*GetSize(b[0]));
153 tmp2.resize(indices.size()*GetSize(b[0]));
154 for(size_t i=0; i<indices.size(); i++)
155 {
156 int j = indices[i];
157 smallvec_type s = b[j];
158 // calc b-Ax
159 for(const_row_iterator it = A.begin_row(j); it != A.end_row(j); ++it)
160 s -= it.value() * x[it.index()];
161 for(size_t u=0; u<GetSize(s); u++)
162 tmp[k++] = BlockRef(s, u);
163 }
164 MatMult(tmp2, 1.0, AlocInv, tmp);
165
166 k=0;
167 for(size_t i=0; i<indices.size(); i++)
168 {
169 smallvec_type &xi = x[indices[i]];
170 for(size_t j=0; j<GetSize(xi); j++)
171 BlockRef(xi, j) += tmp2[k++];
172 }
173}
174template<typename TSparseMatrixType, typename TVectorType>
175void GetBlockGSCorrection(const TSparseMatrixType &A, TVectorType &x, TVectorType &b,
177 std::vector<size_t> &indices)
178{
181 GetBlockGSCorrection(A, x, b, AlocInv, indices, tmp, tmp2);
182}
183
195template<typename TSparseMatrixType, typename TVectorType>
196void GetBlockGSCorrectionILUT(const TSparseMatrixType &A, TVectorType &x, TVectorType &b,
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 tmp.resize(indices.size()*blockSize);
207 tmp2.resize(indices.size()*blockSize);
208 for(size_t i=0; i<indices.size(); i++)
209 {
210 int j = indices[i];
211 smallvec_type s = b[j];
212 // calc b-Ax
213 for(const_row_iterator it = A.begin_row(j); it != A.end_row(j); ++it)
214 s -= it.value() * x[it.index()];
215 for(size_t u=0; u<GetSize(s); u++)
216 tmp[k++] = BlockRef(s, u);
217 }
218 ilut->solve(tmp2, tmp);
219
220 k=0;
221 for(size_t i=0; i<indices.size(); i++)
222 {
223 smallvec_type &xi = x[indices[i]];
224 for(size_t j=0; j<GetSize(xi); j++)
225 BlockRef(xi, j) += tmp2[k++];
226 }
227}
228
234template<typename TSparseMatrixType>
235void GetSliceDenseInverse(const TSparseMatrixType &A, const std::vector<size_t> &indices,
238{
239 tmp.resize(indices.size(), indices.size());
240 GetLocalMatrix(A, tmp, &indices[0], &indices[0]);
241 BlockMatrixToDoubleMatrix(AlocInv, tmp);
242 Invert(AlocInv);
243}
244
245template<typename TSparseMatrixType>
246void GetSliceDenseInverse(const TSparseMatrixType &A, const std::vector<size_t> &indices,
248{
250 GetSliceDenseInverse(A, indices, AlocInv, L);
251}
252
258template<typename TSparseMatrixType>
259void GetSliceSparse(const TSparseMatrixType &A, const std::vector<size_t> &indices,
261{
262 size_t blockSize = GetRows(A(0,0));
263 size_t numRows = indices.size();
264 size_t numCols = indices.size();
265
266 R.resize_and_clear(numRows*blockSize, numCols*blockSize);
267 for(size_t ri=0; ri<indices.size(); ri++)
268 {
269 size_t r=indices[ri];
270 for(size_t ci=0; ci<indices.size(); ci++)
271 {
272 size_t c = indices[ci];
273 if(A.has_connection(r, c))
274 {
275 const typename TSparseMatrixType::value_type &m = A(r, c);
276 for(size_t j1=0; j1<blockSize; j1++)
277 for(size_t j2=0; j2<blockSize; j2++)
278 R(ri*blockSize + j1, ci*blockSize + j2) = BlockRef(m, j1, j2);
279 }
280 }
281 }
282 R.defragment();
283}
284
285
286template<typename 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
301
305
306 protected:
307#ifdef UG_PARALLEL
309#endif
311 {
312
313 matrix_type &mat = *pOp;
314
315#ifdef UG_PARALLEL
316 A = mat;
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 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 virtual bool postprocess() {return true;}
334 virtual bool supports_parallel() const { return true; }
335
337
338 // Stepping routine
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 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
369template <typename TAlgebra, bool backward, bool forward>
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
384
387
389
390 protected:
391 typedef typename matrix_type::value_type block_type;
392
393 public:
394 // Constructor
396 m_depth = 1;
397 };
398
399 BlockGaussSeidel(int depth) {
400 m_depth = depth;
401 };
402
403 BlockGaussSeidel(const this_type &parent) : base_type(parent)
404 {
405 m_depth = parent.m_depth;
406 }
407
408 // Clone
410 {
411 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 void set_depth(size_t d)
424 {
425 m_depth = d;
426 }
427
428 protected:
429 // Name of preconditioner
430 virtual const char* name() const {return "BlockGaussSeidel";}
431
432 // Preprocess routine
434 {
435 size_t N = A.num_rows();
437
438 AlocInv.clear();
439 AlocInv.resize(N);
440
441 indices.resize(N);
442
443 size_t maxSize = 0;
444 PROGRESS_START(prog, N, "BlockGaussSeidel: compute blocks");
445 for(size_t i=0; i<N; i++)
446 {
447 PROGRESS_UPDATE(prog, i);
448 std::vector<bool> bVisited(N, false);
449 indices[i].clear();
450 GetNeighborhood(A, i, m_depth, indices[i], bVisited);
451 //indices[i].push_back(i);
452
453 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 if(AlocInv[i].num_rows() > maxSize) maxSize = AlocInv[i].num_rows();
460 }
461 PROGRESS_FINISH(prog);
462
463 UG_LOG("Max Size = " << maxSize << "\n");
464 return true;
465 }
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 virtual bool postprocess() {return true;}
472 virtual bool supports_parallel() const { return true; }
473
474 // Stepping routine
475 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 b = d;
482
485
486 if(forward)
487 for(size_t i=0; i<x.size(); i++)
488 {
489 // c = D^{-1}(b-Ax)
490 // x = x + c
491 GetBlockGSCorrection(A, x, b, AlocInv[i], indices[i], tmp, tmp2);
492 // do_correction_implicit(A, x, b, indices[i]);
493 }
494 if(backward)
495 for(size_t i=x.size()-1; ; i--)
496 {
497 // c = D^{-1}(b-Ax)
498 // x = x + c
499 GetBlockGSCorrection(A, x, b, AlocInv[i], indices[i], tmp, tmp2);
500 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 return true;
509 }
510
511 virtual std::string config_string() const
512 {
513 std::stringstream ss ; ss << "BlockGaussSeidel(depth = " << m_depth << ")";
514 return ss.str();
515 }
516
517};
518
519
520template <typename TAlgebra, bool backward, bool forward>
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
535
538
540
541 protected:
542 typedef typename matrix_type::value_type block_type;
543 public:
544 // Constructor
546 m_depth = 1;
547 m_nu = 2;
548 };
549
550 BlockGaussSeidelIterative(int depth, int nu) {
551 m_depth = depth;
552 m_nu = nu;
553 };
554
556 {
557 m_depth = parent.m_depth;
558 m_nu = parent.m_nu;
559 }
560
561 // Clone
563 {
564 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 void set_depth(size_t d)
577 {
578 m_depth = d;
579 }
580
581
582 protected:
583 // Name of preconditioner
584 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
595 {
596 size_t N = A.num_rows();
597 indices.resize(N);
598
599 size_t maxSize = 0;
600 for(size_t i=0; i<N; i++)
601 {
602 std::vector<bool> bVisited(N, false);
603 indices[i].clear();
604 GetNeighborhood(A, i, m_depth, indices[i], bVisited);
605 maxSize = std::max(indices[i].size(), maxSize);
606 }
607 UG_LOG("Max Size = " << maxSize << "\n");
608 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 virtual bool postprocess() {return true;}
616 virtual bool supports_parallel() const { return true; }
617
618
619 void correct(size_t i, const matrix_type &A, vector_type& x, const vector_type& b)
620 {
621
622 smallvec_type s = b[i];
623 // calc b-Ax
624 for(const_row_iterator it = A.begin_row(i); it != A.end_row(i); ++it)
625 s -= it.value() * x[it.index()];
627 InverseMatMult(c, 1.0, A(i,i), s);
628 x[i] += c;
629
630 }
631
632 void correct_forward(size_t i, matrix_type &A, vector_type& x, const vector_type& b)
633 {
634 for(size_t k=0; k<m_nu; k++)
635 for(size_t j=0; j<indices[i].size(); j++)
636 correct(indices[i][j], A, x, b);
637 }
638
639 void correct_backward(size_t i, matrix_type &A, vector_type& x, const vector_type& b)
640 {
641 for(size_t k=0; k<m_nu; k++)
642 for(int j=(int)(indices[i].size())-1; j>=0 ; j--)
643 correct(indices[i][j], A, x, b);
644 }
645
646 size_t m_nu;
647
648 // Stepping routine
649 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 b = d;
656
657 if(forward)
658 for(size_t i=0; i<x.size(); i++)
659 {
660 correct_forward(i, A, x, b);
661 }
662 if(backward)
663 for(size_t i=x.size()-1; ; i--)
664 {
665 correct_backward(i, A, x, b);
666 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 return true;
674 }
675
676 public:
677 void set_iterative_steps(size_t nu)
678 {
679 m_nu=nu;
680 }
681
682
683 virtual std::string config_string() const
684 {
685 std::stringstream ss ;
686 if(backward&&forward) ss << "Symmetric";
687 else if(backward) ss << "Backward";
688 ss << "BlockGaussSeidelIterative(depth = " << m_depth << ", nu = " << m_nu << ")";
689 return ss.str();
690 }
691
692};
693
694
701template <typename TAlgebra, bool backward, bool forward>
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
716
719
720 protected:
721 typedef typename matrix_type::value_type block_type;
723
724 public:
725 // Constructor
727 m_depth = 1;
728 };
729
731 m_depth = depth;
732 };
733
735 {
736 m_depth = parent.m_depth;
737 }
738
739 // Clone
741 {
742 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 void set_depth(size_t d)
758 {
759 m_depth = d;
760 }
761
762
763 protected:
764 // Name of preconditioner
765 virtual const char* name() const {return "SparseBlockGaussSeidel";}
766
767 // Preprocess routine
769 {
770 size_t N = A.num_rows();
772
773
774 m_ilut.clear();
775
776 indices.resize(N);
777
778 size_t maxSize = 0;
779 PROGRESS_START(prog, N, "SparseBlockGaussSeidel: compute blocks");
780
781
782
783 std::vector<size_t> bVisited(N, 0);
784 std::vector<bool> bVisited2(N, false);
785
786 std::vector<size_t> levels;
787
788 for(size_t i=0; i<N; i++)
789 {
790// if(bVisited[i]>5) continue;
791 PROGRESS_UPDATE(prog, i);
792
793 indices[i].clear();
794
795 GetNeighborhood(A, i, m_depth, indices[i], bVisited2);
796 for(size_t j=0; j<indices[i].size(); j++)
797 bVisited[indices[i][j]] ++;
798 //for(size_t j=0; j<levels[m_depth-1]; j++)
799 //bVisited2[indices[j]]=true;
800
802
803 GetSliceSparse(A, indices[i], *Aloc[i]);
804
805
807 m_ilut[i]->preprocess_mat(*Aloc[i]);
808
809
810
811 if(Aloc[i]->num_rows() > maxSize) maxSize = Aloc[i]->num_rows();
812 }
813 PROGRESS_FINISH(prog);
814
815 UG_LOG("Max Size = " << maxSize << "\n");
816 return true;
817 }
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 virtual bool postprocess() {return true;}
824 virtual bool supports_parallel() const { return true; }
825
826 // Stepping routine
827 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 b = d;
833
836 CPUAlgebra::vector_type tmp, tmp2;
837 PROGRESS_START(prog, x.size(), "SparseBlockGaussSeidel: step");
838 if(forward)
839 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 if(indices[i].size() != 0)
845 GetBlockGSCorrectionILUT(A, x, b, m_ilut[i], indices[i], tmp, tmp2);
846 }
847 if(backward)
848 for(size_t i=x.size()-1; ; i--)
849 {
850 PROGRESS_UPDATE(prog, i);
851 // c = D^{-1}(b-Ax)
852 // x = x + c
853 if(indices[i].size() != 0)
854 GetBlockGSCorrectionILUT(A, x, b, m_ilut[i], indices[i], tmp, tmp2);
855 if(i==0) break;
856 }
857 PROGRESS_FINISH(prog);
858
859 // Correction is always consistent
860 #ifdef UG_PARALLEL
861 c.set_storage_type(PST_CONSISTENT);
862 #endif
863 return true;
864 }
865
866 virtual std::string config_string() const
867 {
868 std::stringstream ss ;
869 if(backward&&forward) ss << "Symmetric";
870 else if(backward) ss << "Backward";
871 ss << "SparseBlockGaussSeidel(depth = " << m_depth << ")";
872 return ss.str();
873 }
874
875};
876
877
884template <typename TAlgebra, bool backward, bool forward>
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
899
903
904 protected:
905 typedef typename matrix_type::value_type block_type;
906 public:
907 // Constructor
909 m_depth = 1;
910 };
911
913 m_depth = depth;
914 };
915
917 {
918 m_depth = parent.m_depth;
919 }
920
921 // Clone
923 {
924 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 void set_depth(size_t d)
939 {
940 m_depth = d;
941 }
942
943
944 protected:
945 // Name of preconditioner
946 virtual const char* name() const {return "SparseBlockGaussSeidel2";}
947
948 // Preprocess routine
950 {
951 size_t N = A.num_rows();
953 m_ilut.clear();
954
955 indices.resize(N);
956
957 size_t maxSize = 0;
958 PROGRESS_START(prog, N, "SparseBlockGaussSeidel: compute blocks");
959
960 cgraph G;
961 {
962 cgraph graph;
964 G.resize(N);
965 for(size_t i=0; i<graph.size(); i++)
966 for(cgraph::row_iterator it = graph.begin_row(i); it != graph.end_row(i); ++it)
967 {
968 G.set_connection(*it, i);
969 G.set_connection(i, *it);
970 }
971 }
972 std::vector<int> iComponent(N, -1);
973 std::vector<std::vector<int> > components;
974// G.print();
975 for(size_t i=0; i<G.size(); i++)
976 {
977 if(iComponent[i] != -1) continue;
978 int myComponent = iComponent[i] = components.size();
979 components.resize(components.size()+1);
980 std::vector<int> &myComponents = components[myComponent];
981 myComponents.push_back(i);
982 for(size_t c=0; c<myComponents.size(); c++)
983 {
984 int j = myComponents[c];
985 for(cgraph::row_iterator it = G.begin_row(j); it != G.end_row(j); ++it)
986 {
987 if(iComponent[*it] == myComponent) continue;
988 myComponents.push_back(*it);
989 UG_COND_THROW(iComponent[*it] != -1, i );
990 iComponent[*it] = myComponent;
991 }
992 }
993 }
994 for(size_t i=0; i<N; i++)
995 indices[i].clear();
996 for(size_t c=0; c<components.size(); c++)
997 {
998 int i=components[c][0];
999 indices[i].insert(indices[i].begin(), components[c].begin(), components[c].end());
1000 PRINT_VECTOR(indices[i], i);
1001
1003 GetSliceSparse(A, indices[i], *Aloc[i]);
1004
1006 m_ilut[i]->preprocess_mat(*Aloc[i]);
1007 if(Aloc[i]->num_rows() > maxSize) maxSize = Aloc[i]->num_rows();
1008 }
1009 PROGRESS_FINISH(prog);
1010
1011 UG_LOG("Max Size = " << maxSize << "\n");
1012 return true;
1013 }
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 virtual bool postprocess() {return true;}
1020 virtual bool supports_parallel() const { return true; }
1021
1022 // Stepping routine
1023 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 b = d;
1029
1032 CPUAlgebra::vector_type tmp, tmp2;
1033 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 for(size_t i=x.size()-1; ; i--)
1045 {
1046 PROGRESS_UPDATE(prog, i);
1047 // c = D^{-1}(b-Ax)
1048 // x = x + c
1049 if(indices[i].size() != 0)
1050 GetBlockGSCorrectionILUT(A, x, b, m_ilut[i], indices[i], tmp, tmp2);
1051 if(i==0) break;
1052 }
1053 PROGRESS_FINISH(prog);
1054
1055 // Correction is always consistent
1056 #ifdef UG_PARALLEL
1057 c.set_storage_type(PST_CONSISTENT);
1058 #endif
1059 return true;
1060 }
1061
1062 virtual std::string config_string() const
1063 {
1064 std::stringstream ss ;
1065 if(backward&&forward) ss << "Symmetric";
1066 else if(backward) ss << "Backward";
1067 ss << "SparseBlockGaussSeidel(depth = " << m_depth << ")";
1068 return ss.str();
1069 }
1070
1071};
1072
1073
1074} // end namespace ug
1075
1076#endif // __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__GAUSS_SEIDEL__
parameterString s
Definition smart_pointer.h:108
bool valid() const
returns true if the pointer is valid, false if not.
Definition smart_pointer.h:206
Definition block_gauss_seidel.h:371
vector_type::value_type smallvec_type
Definition block_gauss_seidel.h:415
void set_depth(size_t d)
Definition block_gauss_seidel.h:423
std::vector< DenseMatrix< VariableArray2< double > > > AlocInv
Definition block_gauss_seidel.h:418
virtual bool postprocess()
cleans the operator
Definition block_gauss_seidel.h:471
std::vector< std::vector< size_t > > indices
Definition block_gauss_seidel.h:420
matrix_type::value_type smallmat_type
Definition block_gauss_seidel.h:414
virtual SmartPtr< ILinearIterator< vector_type > > clone()
clone
Definition block_gauss_seidel.h:409
virtual bool supports_parallel() const
returns if parallel solving is supported
Definition block_gauss_seidel.h:472
BlockGaussSeidel()
Definition block_gauss_seidel.h:395
virtual const char * name() const
returns the name of iterator
Definition block_gauss_seidel.h:430
matrix_type::const_row_iterator const_row_iterator
Definition block_gauss_seidel.h:416
size_t m_depth
Definition block_gauss_seidel.h:419
virtual bool block_preprocess(matrix_type &A)
Definition block_gauss_seidel.h:433
virtual std::string config_string() const
Definition block_gauss_seidel.h:511
IBlockJacobiPreconditioner< TAlgebra > base_type
Base type.
Definition block_gauss_seidel.h:386
BlockGaussSeidel(int depth)
Definition block_gauss_seidel.h:399
matrix_type::const_row_iterator matrix_const_row_iterator
Definition block_gauss_seidel.h:467
BlockGaussSeidel< algebra_type, backward, forward > this_type
Definition block_gauss_seidel.h:388
TAlgebra::matrix_type matrix_type
Definition block_gauss_seidel.h:380
matrix_type::row_iterator matrix_row_iterator
Definition block_gauss_seidel.h:468
IPreconditioner< TAlgebra >::matrix_operator_type matrix_operator_type
Matrix Operator type.
Definition block_gauss_seidel.h:383
matrix_type::value_type block_type
Definition block_gauss_seidel.h:391
TAlgebra algebra_type
Definition block_gauss_seidel.h:374
TAlgebra::vector_type vector_type
Definition block_gauss_seidel.h:377
BlockGaussSeidel(const this_type &parent)
Definition block_gauss_seidel.h:403
virtual bool block_step(matrix_type &A, vector_type &c, const vector_type &d)
Definition block_gauss_seidel.h:475
Definition block_gauss_seidel.h:522
matrix_type::value_type block_type
Definition block_gauss_seidel.h:542
virtual std::string config_string() const
Definition block_gauss_seidel.h:683
matrix_type::value_type smallmat_type
Definition block_gauss_seidel.h:567
void set_iterative_steps(size_t nu)
Definition block_gauss_seidel.h:677
virtual bool postprocess()
cleans the operator
Definition block_gauss_seidel.h:615
virtual const char * name() const
returns the name of iterator
Definition block_gauss_seidel.h:584
void correct_backward(size_t i, matrix_type &A, vector_type &x, const vector_type &b)
Definition block_gauss_seidel.h:639
vector_type::value_type smallvec_type
Definition block_gauss_seidel.h:568
virtual bool supports_parallel() const
returns if parallel solving is supported
Definition block_gauss_seidel.h:616
std::vector< DenseMatrix< VariableArray2< double > > > AlocInv
Definition block_gauss_seidel.h:571
void correct_forward(size_t i, matrix_type &A, vector_type &x, const vector_type &b)
Definition block_gauss_seidel.h:632
matrix_type::const_row_iterator const_row_iterator
Definition block_gauss_seidel.h:569
virtual bool block_step(matrix_type &A, vector_type &c, const vector_type &d)
Definition block_gauss_seidel.h:649
IBlockJacobiPreconditioner< TAlgebra > base_type
Base type.
Definition block_gauss_seidel.h:537
size_t m_depth
Definition block_gauss_seidel.h:572
std::vector< std::vector< size_t > > indices
Definition block_gauss_seidel.h:573
IPreconditioner< TAlgebra >::matrix_operator_type matrix_operator_type
Matrix Operator type.
Definition block_gauss_seidel.h:534
void correct(size_t i, const matrix_type &A, vector_type &x, const vector_type &b)
Definition block_gauss_seidel.h:619
void set_depth(size_t d)
Definition block_gauss_seidel.h:576
TAlgebra::matrix_type matrix_type
Definition block_gauss_seidel.h:531
BlockGaussSeidelIterative(int depth, int nu)
Definition block_gauss_seidel.h:550
TAlgebra::vector_type vector_type
Definition block_gauss_seidel.h:528
virtual bool block_preprocess(matrix_type &A)
Definition block_gauss_seidel.h:594
TAlgebra algebra_type
Definition block_gauss_seidel.h:525
matrix_type::row_iterator matrix_row_iterator
Definition block_gauss_seidel.h:612
BlockGaussSeidelIterative< algebra_type, backward, forward > this_type
Definition block_gauss_seidel.h:539
matrix_type::const_row_iterator matrix_const_row_iterator
Definition block_gauss_seidel.h:611
size_t m_nu
Definition block_gauss_seidel.h:646
virtual SmartPtr< ILinearIterator< vector_type > > clone()
clone
Definition block_gauss_seidel.h:562
BlockGaussSeidelIterative(const this_type &parent)
Definition block_gauss_seidel.h:555
BlockGaussSeidelIterative()
Definition block_gauss_seidel.h:545
Definition densematrix.h:57
Definition densevector.h:101
Definition block_gauss_seidel.h:288
IBlockJacobiPreconditioner(const IBlockJacobiPreconditioner &parent)
Definition block_gauss_seidel.h:303
IBlockJacobiPreconditioner()
Definition block_gauss_seidel.h:302
virtual bool supports_parallel() const
returns if parallel solving is supported
Definition block_gauss_seidel.h:334
SmartPtr< vector_type > m_spDtmp
Definition block_gauss_seidel.h:336
TAlgebra::vector_type vector_type
Definition block_gauss_seidel.h:294
TAlgebra algebra_type
Definition block_gauss_seidel.h:291
virtual bool preprocess(SmartPtr< MatrixOperator< matrix_type, vector_type > > pOp)
initializes the preconditioner
Definition block_gauss_seidel.h:310
matrix_type A
Definition block_gauss_seidel.h:308
virtual bool step(SmartPtr< MatrixOperator< matrix_type, vector_type > > pOp, vector_type &c, const vector_type &d)
computes a new correction c = B*d
Definition block_gauss_seidel.h:339
virtual bool block_preprocess(matrix_type &A)=0
TAlgebra::matrix_type matrix_type
Definition block_gauss_seidel.h:297
IPreconditioner< TAlgebra >::matrix_operator_type matrix_operator_type
Matrix Operator type.
Definition block_gauss_seidel.h:300
virtual bool block_step(matrix_type &A, vector_type &c, const vector_type &d)=0
virtual bool postprocess()
cleans the operator
Definition block_gauss_seidel.h:333
Definition ilut.h:25
describes a linear iterator that is based on a matrix operator
Definition preconditioner.h:103
Wrapper for sequential matrices to handle them in parallel.
Definition parallel_matrix.h:65
Definition block_gauss_seidel.h:886
virtual bool postprocess()
cleans the operator
Definition block_gauss_seidel.h:1019
SparseBlockGaussSeidel2(int depth)
Definition block_gauss_seidel.h:912
TAlgebra::vector_type vector_type
Definition block_gauss_seidel.h:892
void set_depth(size_t d)
Definition block_gauss_seidel.h:938
SparseBlockGaussSeidel2()
Definition block_gauss_seidel.h:908
matrix_type::value_type smallmat_type
Definition block_gauss_seidel.h:927
matrix_type::const_row_iterator matrix_const_row_iterator
Definition block_gauss_seidel.h:1015
virtual bool block_preprocess(matrix_type &A)
Definition block_gauss_seidel.h:949
IPreconditioner< TAlgebra >::matrix_operator_type matrix_operator_type
Matrix Operator type.
Definition block_gauss_seidel.h:898
vector_type::value_type smallvec_type
Definition block_gauss_seidel.h:928
matrix_type::const_row_iterator const_row_iterator
Definition block_gauss_seidel.h:929
SparseBlockGaussSeidel2< TAlgebra, backward, forward > this_type
Definition block_gauss_seidel.h:902
virtual std::string config_string() const
Definition block_gauss_seidel.h:1062
std::map< size_t, SmartPtr< ILUTPreconditioner< CPUAlgebra > > > m_ilut
Definition block_gauss_seidel.h:931
virtual SmartPtr< ILinearIterator< vector_type > > clone()
clone
Definition block_gauss_seidel.h:922
TAlgebra::matrix_type matrix_type
Definition block_gauss_seidel.h:895
SparseBlockGaussSeidel2(const this_type &parent)
Definition block_gauss_seidel.h:916
std::map< size_t, SmartPtr< CPUAlgebra::matrix_type > > Aloc
Definition block_gauss_seidel.h:935
virtual bool supports_parallel() const
returns if parallel solving is supported
Definition block_gauss_seidel.h:1020
size_t m_depth
Definition block_gauss_seidel.h:933
std::vector< std::vector< size_t > > indices
Definition block_gauss_seidel.h:934
matrix_type::row_iterator matrix_row_iterator
Definition block_gauss_seidel.h:1016
virtual bool block_step(matrix_type &A, vector_type &c, const vector_type &d)
Definition block_gauss_seidel.h:1023
matrix_type::value_type block_type
Definition block_gauss_seidel.h:905
IBlockJacobiPreconditioner< TAlgebra > base_type
Base type.
Definition block_gauss_seidel.h:901
virtual const char * name() const
returns the name of iterator
Definition block_gauss_seidel.h:946
TAlgebra algebra_type
Definition block_gauss_seidel.h:889
Definition block_gauss_seidel.h:703
std::map< size_t, SmartPtr< CPUAlgebra::matrix_type > > Aloc
Definition block_gauss_seidel.h:754
virtual bool supports_parallel() const
returns if parallel solving is supported
Definition block_gauss_seidel.h:824
SparseBlockGaussSeidel()
Definition block_gauss_seidel.h:726
SparseBlockGaussSeidel(const this_type &parent)
Definition block_gauss_seidel.h:734
virtual SmartPtr< ILinearIterator< vector_type > > clone()
clone
Definition block_gauss_seidel.h:740
virtual const char * name() const
returns the name of iterator
Definition block_gauss_seidel.h:765
SparseBlockGaussSeidel< TAlgebra, backward, forward > this_type
Definition block_gauss_seidel.h:722
SparseBlockGaussSeidel(int depth)
Definition block_gauss_seidel.h:730
virtual bool block_preprocess(matrix_type &A)
Definition block_gauss_seidel.h:768
vector_type::value_type smallvec_type
Definition block_gauss_seidel.h:747
void set_depth(size_t d)
Definition block_gauss_seidel.h:757
matrix_type::row_iterator matrix_row_iterator
Definition block_gauss_seidel.h:820
size_t m_depth
Definition block_gauss_seidel.h:752
TAlgebra::matrix_type matrix_type
Definition block_gauss_seidel.h:712
matrix_type::value_type smallmat_type
Definition block_gauss_seidel.h:746
matrix_type::const_row_iterator matrix_const_row_iterator
Definition block_gauss_seidel.h:819
IBlockJacobiPreconditioner< TAlgebra > base_type
Base type.
Definition block_gauss_seidel.h:718
virtual std::string config_string() const
Definition block_gauss_seidel.h:866
matrix_type::const_row_iterator const_row_iterator
Definition block_gauss_seidel.h:748
std::map< size_t, SmartPtr< ILUTPreconditioner< CPUAlgebra > > > m_ilut
Definition block_gauss_seidel.h:750
matrix_type::value_type block_type
Definition block_gauss_seidel.h:721
TAlgebra algebra_type
Definition block_gauss_seidel.h:706
virtual bool postprocess()
cleans the operator
Definition block_gauss_seidel.h:823
IPreconditioner< TAlgebra >::matrix_operator_type matrix_operator_type
Matrix Operator type.
Definition block_gauss_seidel.h:715
virtual bool block_step(matrix_type &A, vector_type &c, const vector_type &d)
Definition block_gauss_seidel.h:827
TAlgebra::vector_type vector_type
Definition block_gauss_seidel.h:709
std::vector< std::vector< size_t > > indices
Definition block_gauss_seidel.h:753
Definition sparse_vector.h:41
Adds 'strong negative connections' to graph.
Definition block_gauss_seidel.h:75
double m_theta
Definition block_gauss_seidel.h:77
StrongNegativeConnectionsByBlockNorm(double theta)
Definition block_gauss_seidel.h:79
void find(const TRowType &Ci, size_t i, cgraph &graph)
Definition block_gauss_seidel.h:82
Definition variable_array.h:59
Definition variable_array.h:139
void resize(size_t newSize, bool bCopyValues=true)
Definition vector.h:109
Definition new_graph.h:69
row_iterator begin_row(size_t row)
Definition new_graph.h:168
row_iterator end_row(size_t row)
Definition new_graph.h:175
size_t size() const
Definition new_graph.h:287
size_t * row_iterator
Definition new_graph.h:72
void resize(size_t n)
Definition new_graph.h:88
void set_connection(size_t from, size_t to)
set a connection from "from" to "to" if not already there
Definition new_graph.h:148
#define PRINT_VECTOR(v, msg)
Definition debug_print.h:53
void MatAddSlaveRowsToMasterRowOverlap0(TMatrix &mat)
Generates a set of unique global algebra ids.
Definition parallelization_util.h:106
@ PST_CONSISTENT
Definition parallel_storage_type.h:68
@ PST_UNIQUE
Definition parallel_storage_type.h:70
@ PST_ADDITIVE
Definition parallel_storage_type.h:69
void GetNeighborhood(const TSparseMatrix &A, size_t node, size_t depth, std::vector< size_t > &indices, std::vector< bool > &bVisited, bool bResetVisitedFlags=true)
Definition sparsematrix_util.h:495
void SetDirichletRow(TSparseMatrix &A, size_t i, size_t alpha)
Definition sparsematrix_util.h:796
#define UG_LOG(msg)
Definition log.h:367
#define UG_COND_THROW(cond, msg)
UG_COND_THROW(cond, msg) : performs a UG_THROW(msg) if cond == true.
Definition error.h:61
the ug namespace
void GetSliceDenseInverse(const TSparseMatrixType &A, const std::vector< size_t > &indices, DenseMatrix< VariableArray2< double > > &AlocInv, DenseMatrix< VariableArray2< typename TSparseMatrixType::value_type > > &tmp)
Definition block_gauss_seidel.h:235
bool Invert(T &m)
double & BlockRef(T &vec, size_t i)
Definition blocks.h:66
void BlockMatrixToDoubleMatrix(DenseMatrix< VariableArray2< double > > &Ad, TMatrix &Ab)
Definition additional_math.h:193
void CopyOffDiagEntries(const TMatrixType &A, size_t i, TRowType &row, bool enforceNew=false)
Definition block_gauss_seidel.h:62
bool GetLocalMatrix(const TGlobalMatrix &mat, TLocalMatrix &localMat)
Definition local_helper.h:191
void GetBlockGSCorrectionILUT(const TSparseMatrixType &A, TVectorType &x, TVectorType &b, SmartPtr< ILUTPreconditioner< CPUAlgebra > > &ilut, std::vector< size_t > &indices, CPUAlgebra::vector_type &tmp, CPUAlgebra::vector_type &tmp2)
Definition block_gauss_seidel.h:196
size_t GetRows(const T &t)
bool MatMult(vector_t &dest, const number &beta1, const matrix_t &A1, const vector_t &w1)
calculates dest = beta1 * A1;
Definition operations_mat.h:59
bool InverseMatMult(number &dest, const double &beta, const TMat &mat, const TVec &vec)
you can implement this function with GetInverse and MatMult
size_t GetSize(const T &t)
double BlockNorm(const TYPE &v)
Definition blocks.h:57
void GetSliceSparse(const TSparseMatrixType &A, const std::vector< size_t > &indices, CPUAlgebra::matrix_type &R)
Definition block_gauss_seidel.h:259
void CreateStrongConnectionGraphForSystems(const matrix_type &A, cgraph &graph, double theta)
Definition block_gauss_seidel.h:106
void GetBlockGSCorrection(const TSparseMatrixType &A, TVectorType &x, TVectorType &b, DenseMatrix< VariableArray2< double > > &AlocInv, std::vector< size_t > &indices, DenseVector< VariableArray1< double > > &tmp, DenseVector< VariableArray1< double > > &tmp2)
Definition block_gauss_seidel.h:141
#define PROFILE_BEGIN(name)
Definition profiler.h:254
#define PROFILE_FUNC()
Definition profiler.h:257
#define PROGRESS_FINISH(progVarName)
Definition progress.h:118
#define PROGRESS_START(progVarName, dSize, msg)
Definition progress.h:111
#define PROGRESS_UPDATE(progVarName, d)
Definition progress.h:117
SmartPtr< T, FreePolicy > make_sp(T *inst)
returns a SmartPtr for the passed raw pointer
Definition smart_pointer.h:836
size_t num_rows() const
Definition sparsematrix_interface.h:38