Line data Source code
1 : /*
2 : * Copyright (c) 2010-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__CPU_ALGEBRA__CORE_SMOOTHERS__
34 : #define __H__UG__CPU_ALGEBRA__CORE_SMOOTHERS__
35 : ////////////////////////////////////////////////////////////////////////////////////////////////
36 :
37 : namespace ug
38 : {
39 :
40 : /// \addtogroup lib_algebra
41 : /// @{
42 :
43 : /////////////////////////////////////////////////////////////////////////////////////////////
44 : /// Gauss-Seidel-Iterations
45 : /**
46 : * Here, iteration schemes of gauss-seidel-type are described for solving the
47 : * linear equation
48 : *
49 : * \f$ A * x = b. A \in R^{nxn}, x \in R^n, b \in R^n \f$.
50 : *
51 : * Most of the common linear iteration-methods base on the decomposition of A into
52 : * its diagonal (D) and strict-upper(-U) and strict-lower part (-L),
53 : *
54 : * \f$ A = D - L - U \f$.
55 : *
56 : * Among others, W. Hackbusch ('Iterative Loesung grosser Gleichungssysteme'),
57 : * distinguishes three different forms for describing a linear iteration scheme.
58 : * The general 'first normal-form' of a linear iteration scheme takes the form
59 : *
60 : * \f$ x^{m+1} = M * x^m + N * b \f$,
61 : *
62 : * with some Matrices \f$ M \f$ and \f$ N \in R^{nxn} \f$. m denotes the iteration index.
63 : * The general 'second normal-form' of a linear iteration scheme takes the form
64 : *
65 : * \f$ x^{m+1} = x^m - N * (A * x^m - b) \f$.
66 : *
67 : * Those linear iteration schemes, which can be represented by the second normal-form
68 : * are the linear, consistent iteration schemes.
69 : *
70 : * Introducing the correction \f$ c{m+1} := x^{m+1} - x^m \f$ and the defect
71 : * \f$ d^m := b - A * x^m \f$ the second normal-form can be rewritten as
72 : *
73 : * \f$ c = N * d \f$.
74 : *
75 : * Below, methods for the (forward) Gauss-Seidel step, the backward Gauss-Seidel step and the symmetric
76 : * Gauss-Seidel step are implemented ('gs_step_LL', 'gs_step_UR' resp. 'sgs_step').
77 : * The matrices of the second normal-form are
78 : *
79 : * \f$ N = (D-L)^{-1}\f$ for the (forward) Gauss-Seidel step,
80 : * \f$ N = (D-U)^{-1}\f$ for the backward Gauss-Seidel step and
81 : * \f$ N = (D-u)^{-1}D(D-U)^{-1} \f$ for the symmetric Gauss-Seidel step.
82 : *
83 : * References:
84 : * <ul>
85 : * <li> W. Hackbusch. Iterative Loesung grosser Gleichungssysteme
86 : * </ul>
87 : *
88 : * \sa gs_step_LL, gs_step_UR, sgs_step
89 : */
90 :
91 :
92 : /////////////////////////////////////////////////////////////////////////////////////////////
93 : // gs_step_LL
94 : /** \brief Performs a forward gauss-seidel-step, that is, solve on the lower left of A.
95 : * Using gs_step_LL within a preconditioner-scheme leads to the fact that we get the correction
96 : * by successive inserting the already computed values of c in c = N * d, with c being the correction
97 : * and d being the defect. N denotes the matrix of the second normal-form of a linear iteration
98 : * scheme.
99 : *
100 : * \param A Matrix \f$A = D - L - U\f$
101 : * \param c Vector. \f$ c = N * d = (D-L)^{-1} * d \f$
102 : * \param d Vector d.
103 : * \sa gs_step_UR, sgs_step
104 : */
105 : template<typename Matrix_type, typename Vector_type>
106 0 : void gs_step_LL(const Matrix_type &A, Vector_type &c, const Vector_type &d, const number relaxFactor)
107 : {
108 : // gs LL has preconditioning matrix N = (D-L)^{-1}
109 :
110 : typedef typename Matrix_type::value_type matrix_block;
111 : typedef typename Matrix_type::const_row_iterator const_row_it;
112 : typename Vector_type::value_type s;
113 :
114 : const size_t sz = c.size();
115 0 : for (size_t i = 0; i < sz; ++i)
116 : {
117 0 : s = d[i];
118 :
119 : // loop over all lower left matrix entries.
120 : // Note: Here the corrections c, which have already been computed in previous loops (wrt. i),
121 : // are taken to compute the i-th correction. For example the correction of the second row
122 : // is computed by s[2] = (d[2] - A[2][1] * c[1]); and c[2] = s[2]/A[2][2];
123 : const const_row_it rowEnd = A.end_row(i);
124 : const_row_it it = A.begin_row(i);
125 0 : for(; it != rowEnd && it.index() < i; ++it)
126 : // s -= it.value() * c[it.index()];
127 0 : MatMultAdd(s, 1.0, s, -1.0, it.value(), c[it.index()]);
128 :
129 : // c[i] = relaxFactor * s/A(i,i)
130 0 : const matrix_block& A_ii = it.index() == i ? it.value() : matrix_block(0);
131 : InverseMatMult(c[i], relaxFactor, A_ii, s);
132 : }
133 0 : }
134 :
135 : /////////////////////////////////////////////////////////////////////////////////////////////
136 : // gs_step_UR
137 : /**
138 : * \brief Performs a backward gauss-seidel-step, that is, solve on the upper right of A.
139 : * Using gs_step_UR within a preconditioner-scheme leads to the fact that we get the correction
140 : * by successive inserting the already computed values of c in c = N * d, with c being the correction
141 : * and d being the defect. N denotes the matrix of the second normal-form of a linear iteration
142 : * scheme.
143 : *
144 : * \param A Matrix \f$A = D - L - U\f$
145 : * \param c will be \f$c = N * d = (D-U)^{-1} * d \f$
146 : * \param d the vector d.
147 : * \sa gs_step_LL, sgs_step
148 : */
149 : template<typename Matrix_type, typename Vector_type>
150 0 : void gs_step_UR(const Matrix_type &A, Vector_type &c, const Vector_type &d, const number relaxFactor)
151 : {
152 : // gs UR has preconditioning matrix N = (D-U)^{-1}
153 :
154 : typename Vector_type::value_type s;
155 :
156 0 : if(c.size() == 0) return;
157 0 : size_t i = c.size()-1;
158 : do
159 : {
160 0 : s = d[i];
161 : typename Matrix_type::const_row_iterator diag = A.get_connection(i, i);
162 :
163 : typename Matrix_type::const_row_iterator it = diag; ++it;
164 0 : for(; it != A.end_row(i); ++it)
165 : // s -= it.value() * x[it.index()];
166 0 : MatMultAdd(s, 1.0, s, -1.0, it.value(), c[it.index()]);
167 :
168 : // c[i] = relaxFactor * s/A(i,i)
169 : InverseMatMult(c[i], relaxFactor, diag.value(), s);
170 0 : } while(i-- != 0);
171 :
172 : }
173 :
174 : /////////////////////////////////////////////////////////////////////////////////////////////
175 : // sgs_step
176 : /**
177 : * \brief Performs a symmetric gauss-seidel step.
178 : * Using sgs_step within a preconditioner-scheme leads to the fact that we get the correction
179 : * by successive inserting the already computed values of c in c = N * d, with c being the correction
180 : * and d being the defect. N denotes the matrix of the second normal-form of a linear iteration
181 : * scheme.
182 : *
183 : * \param A Matrix \f$A = D - L - R\f$
184 : * \param c will be \f$c = N * d = (D-U)^{-1} D (D-L)^{-1} d \f$
185 : * \param d the vector d.
186 : * \sa gs_step_LL, gs_step_LL
187 : */
188 : template<typename Matrix_type, typename Vector_type>
189 0 : void sgs_step(const Matrix_type &A, Vector_type &c, const Vector_type &d, const number relaxFactor)
190 : {
191 : // sgs has preconditioning matrix N = (D-U)^{-1} D (D-L)^{-1}
192 :
193 : // c1 = (D-L)^{-1} d
194 0 : gs_step_LL(A, c, d, relaxFactor);
195 :
196 : // c2 = D c1
197 : typename Vector_type::value_type s;
198 0 : for(size_t i = 0; i<c.size(); i++)
199 : {
200 0 : s=c[i];
201 0 : MatMult(c[i], 1.0, A(i, i), s);
202 : }
203 :
204 : // c3 = (D-U)^{-1} c2
205 0 : gs_step_UR(A, c, c, relaxFactor);
206 0 : }
207 :
208 : /////////////////////////////////////////////////////////////////////////////////////////////
209 : // diag_step
210 : /**
211 : * \brief Performs a jacobi-step
212 : * \param A Matrix \f$A = D - L - R\f$
213 : * \param c will be \f$c = N * d = D^{-1} d \f$
214 : * \param d the vector d.
215 : * \param damp the damping factor
216 : */
217 : template<typename Matrix_type, typename Vector_type>
218 : void diag_step(const Matrix_type& A, Vector_type& c, const Vector_type& d, number damp)
219 : {
220 : //exit(3);
221 : UG_ASSERT(c.size() == d.size() && c.size() == A.num_rows(), c << ", " << d <<
222 : " and " << A << " need to have same size.");
223 :
224 : for(size_t i=0; i < c.size(); i++)
225 : // c[i] = damp * d[i]/A(i,i)
226 : InverseMatMult(c[i], damp, A(i,i), d[i]);
227 : }
228 :
229 :
230 : /// @}
231 : }
232 : #endif // __H__UG__CPU_ALGEBRA__CORE_SMOOTHERS__
|
