Line data Source code
1 : /*
2 : * Copyright (c) 2013-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Raphael Prohl
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_ALGEBRA__OPERATOR__PRECONDITIONER__PROJECTED_GAUSS_SEIDEL__PROJ_GAUSS_SEIDEL_IMPL__
34 : #define __H__UG__LIB_ALGEBRA__OPERATOR__PRECONDITIONER__PROJECTED_GAUSS_SEIDEL__PROJ_GAUSS_SEIDEL_IMPL__
35 :
36 : #include "proj_gauss_seidel.h"
37 :
38 : namespace ug{
39 :
40 : /// commmon GaussSeidel-step-calls for a single index 'i'
41 : template<typename Matrix_type, typename Vector_type>
42 0 : void forward_gs_step(Vector_type& c, const Matrix_type& A, const Vector_type& d,
43 : const size_t i, const number relaxFactor)
44 : {
45 0 : typename Vector_type::value_type s = d[i];
46 :
47 : for(typename Matrix_type::const_row_iterator it = A.begin_row(i);
48 0 : it != A.end_row(i) && it.index() < i; ++it)
49 : // s -= it.value() * x[it.index()];
50 0 : MatMultAdd(s, 1.0, s, -1.0, it.value(), c[it.index()]);
51 :
52 : // c[i] = relaxFactor * s / A(i,i)
53 0 : InverseMatMult(c[i], relaxFactor, A(i,i), s);
54 0 : }
55 :
56 : template<typename Matrix_type, typename Vector_type>
57 0 : void backward_gs_step(Vector_type& c, const Matrix_type& A, const Vector_type& d,
58 : const size_t i, const number relaxFactor)
59 : {
60 0 : typename Vector_type::value_type s = d[i];
61 :
62 : typename Matrix_type::const_row_iterator diag = A.get_connection(i, i);
63 : typename Matrix_type::const_row_iterator it = diag; ++it;
64 :
65 0 : for(; it != A.end_row(i); ++it)
66 : // s -= it.value() * x[it.index()];
67 0 : MatMultAdd(s, 1.0, s, -1.0, it.value(), c[it.index()]);
68 :
69 : // c[i] = relaxFactor * s/A(i,i)
70 : InverseMatMult(c[i], relaxFactor, diag.value(), s);
71 0 : }
72 :
73 :
74 : template <typename TDomain, typename TAlgebra>
75 : void
76 0 : ProjGaussSeidel<TDomain,TAlgebra>::
77 : step(const matrix_type& A, vector_type& c, const vector_type& d, const number relax)
78 : {
79 0 : for(size_t i = 0; i < c.size(); i++)
80 : {
81 0 : forward_gs_step(c, A, d, i, relax);
82 :
83 : // project correction on the subspace defined by the obstacle constraints
84 0 : this->project_correction(c[i], i);
85 : }
86 0 : }
87 :
88 : template <typename TDomain, typename TAlgebra>
89 : void
90 0 : ProjBackwardGaussSeidel<TDomain,TAlgebra>::
91 : step(const matrix_type& A, vector_type& c, const vector_type& d, const number relax)
92 : {
93 0 : if(c.size() == 0) return;
94 0 : size_t i = c.size()-1;
95 : do
96 : {
97 0 : backward_gs_step(c, A, d, i, relax);
98 :
99 : // project correction on the subspace defined by the obstacle constraints
100 0 : this->project_correction(c[i], i);
101 :
102 0 : } while(i-- != 0);
103 : }
104 :
105 : template <typename TDomain, typename TAlgebra>
106 : void
107 0 : ProjSymmetricGaussSeidel<TDomain,TAlgebra>::
108 : step(const matrix_type& A, vector_type& c, const vector_type& d, const number relax)
109 : {
110 0 : for(size_t i = 0; i < c.size(); i++)
111 : {
112 : // 1. perform a forward GaussSeidel step
113 : // c1 = (D-L)^{-1} d
114 0 : forward_gs_step(c, A, d, i, relax);
115 :
116 : // 2. c2 = D c1
117 0 : MatMult(c[i], 1.0, A(i, i), c[i]);
118 :
119 : // 3. perform a backward GaussSeidel step
120 : // c3 = (D-U)^{-1} c2
121 0 : backward_gs_step(c, A, c, i, relax);
122 :
123 : // project correction on the subspace defined by the obstacle constraints
124 0 : this->project_correction(c[i], i);
125 : }
126 0 : }
127 :
128 : } // end namespace ug
129 :
130 : #endif /* __H__UG__LIB_ALGEBRA__OPERATOR__PRECONDITIONER__PROJECTED_GAUSS_SEIDEL__PROJ_GAUSS_SEIDEL_IMPL__ */
|
