ug4
Loading...
Searching...
No Matches
math_matrix_vector_functions_common_impl.hpp
Go to the documentation of this file.
1/*
2 * Copyright (c) 2009-2015: G-CSC, Goethe University Frankfurt
3 * Author: Andreas Vogel
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__LGMATH__MATRIX_VECTOR_FUNCTIONS_COMMON_IMPL__
34#define __H__LGMATH__MATRIX_VECTOR_FUNCTIONS_COMMON_IMPL__
35
36#include <cmath>
37#include "common/error.h"
38#include "math_matrix.h"
39#include "math_vector.h"
40
41namespace ug
42{
43
45// vOut = m * v
46template <typename vector_t_out, typename matrix_t, typename vector_t_in>
47inline
48void
49MatVecMult(vector_t_out& vOut, const matrix_t& m, const vector_t_in& v)
50{
51 assert(vector_t_out::Size == matrix_t::RowSize);
52 assert(vector_t_in::Size == matrix_t::ColSize);
53
54 typedef typename matrix_t::size_type size_type;
55 for(size_type i = 0; i < vOut.size(); ++i)
56 {
57 vOut[i] = 0.0;
58 for(size_type j = 0; j < v.size(); ++j)
59 {
60 vOut[i] += m(i,j) * v[j];
61 }
62 }
63}
64
66// vOut += m * v
67template <typename vector_t_out, typename matrix_t, typename vector_t_in>
68inline
69void
70MatVecMultAppend(vector_t_out& vOut, const matrix_t& m, const vector_t_in& v)
71{
72 assert(vector_t_out::Size == matrix_t::RowSize);
73 assert(vector_t_in::Size == matrix_t::ColSize);
74
75 typedef typename matrix_t::size_type size_type;
76 for(size_type i = 0; i < vOut.size(); ++i)
77 {
78 for(size_type j = 0; j < v.size(); ++j)
79 {
80 vOut[i] += m(i,j) * v[j];
81 }
82 }
83}
84
86// vOut += s * m * v
87template <typename vector_t_out, typename matrix_t, typename vector_t_in>
88inline
89void
90MatVecScaleMultAppend(vector_t_out& vOut, typename vector_t_out::value_type s, const matrix_t& m, const vector_t_in& v)
91{
92 assert(vector_t_out::Size == matrix_t::RowSize);
93 assert(vector_t_in::Size == matrix_t::ColSize);
94
95 typedef typename matrix_t::size_type size_type;
96 for(size_type i = 0; i < vOut.size(); ++i)
97 {
98 for(size_type j = 0; j < v.size(); ++j)
99 {
100 vOut[i] += s * m(i,j) * v[j];
101 }
102 }
103}
104
105
107// vOut = Transpose(m) * v
108template <typename vector_t_out, typename matrix_t, typename vector_t_in>
109inline
110void
111TransposedMatVecMult(vector_t_out& vOut, const matrix_t& m, const vector_t_in& v)
112{
113 assert(vector_t_out::Size == matrix_t::ColSize);
114 assert(vector_t_in::Size == matrix_t::RowSize);
115
116 typedef typename matrix_t::size_type size_type;
117 for(size_type i = 0; i < vOut.size(); ++i)
118 {
119 vOut[i] = 0.0;
120 for(size_type j = 0; j < v.size(); ++j)
121 {
122 vOut[i] += m(j,i) * v[j];
123 }
124 }
125}
126
128// vOut += Transpose(m) * v
129template <typename vector_t_out, typename matrix_t, typename vector_t_in>
130inline
131void
132TransposedMatVecMultAdd(vector_t_out& vOut, const matrix_t& m, const vector_t_in& v)
133{
134 assert(vector_t_out::Size == matrix_t::ColSize);
135 assert(vector_t_in::Size == matrix_t::RowSize);
136
137 typedef typename matrix_t::size_type size_type;
138 for(size_type i = 0; i < vOut.size(); ++i)
139 {
140 for(size_type j = 0; j < v.size(); ++j)
141 {
142 vOut[i] += m(j,i) * v[j];
143 }
144 }
145}
146
155template <typename matrix_t, typename vector_t>
156inline
157void
158GivensMatVecMult (matrix_t& A, vector_t& v)
159{
160 typedef typename matrix_t::size_type size_type;
161 typedef typename matrix_t::value_type value_type;
162
163 assert (vector_t::Size == matrix_t::RowSize);
164 assert (matrix_t::RowSize >= matrix_t::ColSize);
165
166 value_type s, c;
167 value_type d, x, y;
168
169 for (size_type i = 0; i < matrix_t::RowSize - 1; i++) // the 1st index of the elementar rotation
170 for (size_type k = matrix_t::RowSize - 1; k > i; k--) // the 2nd index of this rotation
171 {
172 d = A (i, i); x = A (k, i);
173
174 // Computation of the entries of the elementar transformation:
175 if (std::abs (d) > std::abs (x))
176 {
177 s = x / d;
178 c = 1.0 / std::sqrt (1 + s * s);
179 s *= c;
180 }
181 else if (x != 0) // to ensure that A (k, i) != 0; note that abs(x) >= abs(d) here!
182 {
183 c = d / x;
184 s = 1.0 / std::sqrt (1 + c * c);
185 c *= s;
186 }
187 else continue; // nothing to eliminate
188
189 // Multiplication of A by the elementar transformation:
190 for (size_type j = i; j < matrix_t::ColSize; j++)
191 {
192 x = A (i, j); y = A (k, j);
193 A (i, j) = c * x + s * y;
194 A (k, j) = c * y - s * x;
195 }
196
197 // Multiplication of v by the elementar transformation:
198 x = v [i]; y = v [k];
199 v [i] = c * x + s * y;
200 v [k] = c * y - s * x;
201 }
202}
203
220template <typename matrix_t, typename vector_t>
221inline
222void
223InvMatVecMult_byGivens (matrix_t& A, vector_t& v)
224{
225 typedef typename matrix_t::size_type size_type;
226
227// I. Multiply 'this' by the Givens rotation:
228 GivensMatVecMult (A, v);
229
230// II. Computation of the result:
231 size_type i = matrix_t::ColSize; // <= matrix_t::RowSize, i.e. we invert only the square block
232 do
233 {
234 i--;
235 for (size_type j = i + 1; j < matrix_t::ColSize; j++)
236 v [i] -= A (i, j) * v [j];
237 if (std::abs (A (i, i)) < SMALL * std::abs (v [i]))
238 UG_THROW ("InvMatVecMult_byGivens: Inverting singular matrix by the Givens rotations");
239 v [i] /= A (i, i);
240 }
241 while (i != 0);
242}
243
248template <typename matrix_t, typename vector_t>
249inline
250void
251OrthogProjectVec (vector_t& v, const matrix_t& A)
252{
253 typedef typename matrix_t::size_type size_type;
254 typedef typename matrix_t::value_type value_type;
255
256// I. Solve the least square problem:
257 matrix_t M = A; // we do not work with the original matrix; otherwise it would be destroyed
259
260// II. Compute the linear combination of the columns:
262 for (size_type i = 0; i < matrix_t::ColSize; i++) coeff [i] = v [i];
263 MatVecMult (v, A, coeff);
264}
265
266} // end of namespace ug
267
268#endif /* __H__LGMATH__LGMATH_MATRIX_VECTOR_FUNCTIONS_COMMON_IMPL__ */
parameterString s
a mathematical Vector with N entries.
Definition math_vector.h:97
#define UG_THROW(msg)
Definition error.h:57
void TransposedMatVecMultAdd(vector_t_out &vOut, const matrix_t &m, const vector_t_in &v)
Transposed Matrix - Vector Muliplication.
Definition math_matrix_vector_functions_common_impl.hpp:132
void MatVecMultAppend(vector_t_out &vOut, const matrix_t &m, const vector_t_in &v)
Matrix - Vector Multiplication adding to a second vector.
Definition math_matrix_vector_functions_common_impl.hpp:70
void InvMatVecMult_byGivens(matrix_t &A, vector_t &v)
Multiplication by the inverse using the Givens rotations.
Definition math_matrix_vector_functions_common_impl.hpp:223
void GivensMatVecMult(matrix_t &A, vector_t &v)
Multiplication by the Givens rotation of the QR-Decomposition.
Definition math_matrix_vector_functions_common_impl.hpp:158
void TransposedMatVecMult(vector_t_out &vOut, const matrix_t &m, const vector_t_in &v)
Transposed Matrix - Vector Muliplication.
Definition math_matrix_vector_functions_common_impl.hpp:111
void OrthogProjectVec(vector_t &v, const matrix_t &A)
Orthogonal projection.
Definition math_matrix_vector_functions_common_impl.hpp:251
void MatVecScaleMultAppend(vector_t_out &vOut, typename vector_t_out::value_type s, const matrix_t &m, const vector_t_in &v)
Matrix - Vector Multiplication added scaled to a second vector.
Definition math_matrix_vector_functions_common_impl.hpp:90
void MatVecMult(vector_t_out &vOut, const matrix_t &m, const vector_t_in &v)
Matrix - Vector Multiplication.
Definition math_matrix_vector_functions_common_impl.hpp:49
the ug namespace
const number SMALL
Definition math_constants.h:41
T value_type
Definition sparsematrix_interface.h:2