ug4
Loading...
Searching...
No Matches
proj_gauss_seidel_impl.h
Go to the documentation of this file.
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
38namespace ug{
39
41template<typename Matrix_type, typename Vector_type>
42void forward_gs_step(Vector_type& c, const Matrix_type& A, const Vector_type& d,
43 const size_t i, const number relaxFactor)
44{
45 typename Vector_type::value_type s = d[i];
46
47 for(typename Matrix_type::const_row_iterator it = A.begin_row(i);
48 it != A.end_row(i) && it.index() < i; ++it)
49 // s -= it.value() * x[it.index()];
50 MatMultAdd(s, 1.0, s, -1.0, it.value(), c[it.index()]);
51
52 // c[i] = relaxFactor * s / A(i,i)
53 InverseMatMult(c[i], relaxFactor, A(i,i), s);
54}
55
56template<typename Matrix_type, typename Vector_type>
57void backward_gs_step(Vector_type& c, const Matrix_type& A, const Vector_type& d,
58 const size_t i, const number relaxFactor)
59{
60 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 for(; it != A.end_row(i); ++it)
66 // s -= it.value() * x[it.index()];
67 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}
72
73
74template <typename TDomain, typename TAlgebra>
75void
77step(const matrix_type& A, vector_type& c, const vector_type& d, const number relax)
78{
79 for(size_t i = 0; i < c.size(); i++)
80 {
81 forward_gs_step(c, A, d, i, relax);
82
83 // project correction on the subspace defined by the obstacle constraints
84 this->project_correction(c[i], i);
85 }
86}
87
88template <typename TDomain, typename TAlgebra>
89void
91step(const matrix_type& A, vector_type& c, const vector_type& d, const number relax)
92{
93 if(c.size() == 0) return;
94 size_t i = c.size()-1;
95 do
96 {
97 backward_gs_step(c, A, d, i, relax);
98
99 // project correction on the subspace defined by the obstacle constraints
100 this->project_correction(c[i], i);
101
102 } while(i-- != 0);
103}
104
105template <typename TDomain, typename TAlgebra>
106void
108step(const matrix_type& A, vector_type& c, const vector_type& d, const number relax)
109{
110 for(size_t i = 0; i < c.size(); i++)
111 {
112 // 1. perform a forward GaussSeidel step
113 // c1 = (D-L)^{-1} d
114 forward_gs_step(c, A, d, i, relax);
115
116 // 2. c2 = D c1
117 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 backward_gs_step(c, A, c, i, relax);
122
123 // project correction on the subspace defined by the obstacle constraints
124 this->project_correction(c[i], i);
125 }
126}
127
128} // end namespace ug
129
130#endif /* __H__UG__LIB_ALGEBRA__OPERATOR__PRECONDITIONER__PROJECTED_GAUSS_SEIDEL__PROJ_GAUSS_SEIDEL_IMPL__ */
parameterString s
algebra_type::vector_type vector_type
Vector type.
Definition proj_gauss_seidel.h:129
algebra_type::matrix_type matrix_type
Matrix type.
Definition proj_gauss_seidel.h:126
virtual void step(const matrix_type &mat, vector_type &c, const vector_type &d, const number relax)
computes a new correction c = B*d and projects on the underlying constraint
Definition proj_gauss_seidel_impl.h:91
algebra_type::matrix_type matrix_type
Matrix type.
Definition proj_gauss_seidel.h:88
algebra_type::vector_type vector_type
Vector type.
Definition proj_gauss_seidel.h:91
virtual void step(const matrix_type &mat, vector_type &c, const vector_type &d, const number relax)
computes a new correction c = B*d and projects on the underlying constraint
Definition proj_gauss_seidel_impl.h:77
virtual void step(const matrix_type &mat, vector_type &c, const vector_type &d, const number relax)
computes a new correction c = B*d and projects on the underlying constraint
Definition proj_gauss_seidel_impl.h:108
algebra_type::matrix_type matrix_type
Matrix type.
Definition proj_gauss_seidel.h:166
algebra_type::vector_type vector_type
Vector type.
Definition proj_gauss_seidel.h:169
double number
Definition types.h:124
the ug namespace
void backward_gs_step(Vector_type &c, const Matrix_type &A, const Vector_type &d, const size_t i, const number relaxFactor)
Definition proj_gauss_seidel_impl.h:57
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
void forward_gs_step(Vector_type &c, const Matrix_type &A, const Vector_type &d, const size_t i, const number relaxFactor)
commmon GaussSeidel-step-calls for a single index 'i'
Definition proj_gauss_seidel_impl.h:42
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