ug4
core_smoothers.h
Go to the documentation of this file.
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__
36 
37 namespace ug
38 {
39 
42 
44 // gs_step_LL
52 template<typename Matrix_type, typename Vector_type>
53 bool gs_step_LL(const Matrix_type &A, Vector_type &x, const Vector_type &b)
54 {
55  // gs LL has preconditioning matrix
56  // (D-L)^{-1}
57 
58  typename Vector_type::value_type s;
59 
60  for(size_t i=0; i < x.size(); i++)
61  {
62  s = b[i];
63 
64  for(typename Matrix_type::const_row_iterator it = A.begin_row(i); it != A.end_row(i) && it.index() < i; ++it)
65  // s -= it.value() * x[it.index()];
66  MatMultAdd(s, 1.0, s, -1.0, it.value(), x[it.index()]);
67 
68  // x[i] = s/A(i,i)
69  InverseMatMult(x[i], 1.0, A(i,i), s);
70  }
71 
72  return true;
73 }
74 
76 // gs_step_UR
84 template<typename Matrix_type, typename Vector_type>
85 bool gs_step_UR(const Matrix_type &A, Vector_type &x, const Vector_type &b)
86 {
87  // gs UR has preconditioning matrix
88  // (D-U)^{-1}
89  typename Vector_type::value_type s;
90 
91  if(x.size() == 0) return true;
92  size_t i = x.size()-1;
93  do
94  {
95  s = b[i];
96  typename Matrix_type::const_row_iterator diag = A.get_connection(i, i);
97 
98  typename Matrix_type::const_row_iterator it = diag; ++it;
99  for(; it != A.end_row(i); ++it)
100  // s -= it.value() * x[it.index()];
101  MatMultAdd(s, 1.0, s, -1.0, it.value(), x[it.index()]);
102 
103  // x[i] = s/A(i,i)
104  InverseMatMult(x[i], 1.0, diag.value(), s);
105  } while(i-- != 0);
106 
107  return true;
108 }
109 
111 // sgs_step
119 template<typename Matrix_type, typename Vector_type>
120 bool sgs_step(const Matrix_type &A, Vector_type &x, const Vector_type &b)
121 {
122  // sgs has preconditioning matrix
123  // (D-U)^{-1} D (D-L)^{-1}
124 
125  // x1 = (D-L)^{-1} b
126  gs_step_LL(A, x, b);
127 
128  // x2 = D x1
129  for(size_t i = 0; i<x.size(); i++)
130  MatMult(x[i], 1.0, A(i, i), x[i]);
131 
132  // x3 = (D-U)^{-1} x2
133  gs_step_UR(A, x, x);
134 
135  return true;
136 }
137 
139 // diag_step
147 template<typename Matrix_type, typename Vector_type>
148 bool diag_step(const Matrix_type& A, Vector_type& x, const Vector_type& b, number damp)
149 {
150  //exit(3);
151  UG_ASSERT(x.size() == b.size() && x.size() == A.num_rows(), x << ", " << b << " and " << A << " need to have same size.");
152 
153  for(size_t i=0; i < x.size(); i++)
154  // x[i] = damp * b[i]/A(i,i)
155  InverseMatMult(x[i], damp, A(i,i), b[i]);
156 
157  return true;
158 }
159 
161 }
162 #endif // __H__UG__CPU_ALGEBRA__CORE_SMOOTHERS__
parameterString s
void gs_step_UR(const Matrix_type &A, Vector_type &c, const Vector_type &d, const number relaxFactor)
Performs a backward gauss-seidel-step, that is, solve on the upper right of A. Using gs_step_UR withi...
Definition: core_smoothers.h:150
void gs_step_LL(const Matrix_type &A, Vector_type &c, const Vector_type &d, const number relaxFactor)
Gauss-Seidel-Iterations.
Definition: core_smoothers.h:106
void diag_step(const Matrix_type &A, Vector_type &c, const Vector_type &d, number damp)
Performs a jacobi-step.
Definition: core_smoothers.h:218
void sgs_step(const Matrix_type &A, Vector_type &c, const Vector_type &d, const number relaxFactor)
Performs a symmetric gauss-seidel step. Using sgs_step within a preconditioner-scheme leads to the fa...
Definition: core_smoothers.h:189
#define UG_ASSERT(expr, msg)
Definition: assert.h:70
double number
Definition: types.h:124
the ug namespace
bool MatMultAdd(vector_t &dest, const number &alpha1, const vector_t &v1, const number &beta1, const matrix_t &A1, const vector_t &w1)
calculates dest = alpha1*v1 + beta1 * A1 *w1;
Definition: operations_mat.h:68
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
T value_type
Definition: sparsematrix_interface.h:2