Line data Source code
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 :
41 : namespace ug
42 : {
43 :
44 : /// Matrix - Vector Multiplication
45 : // vOut = m * v
46 : template <typename vector_t_out, typename matrix_t, typename vector_t_in>
47 : inline
48 : void
49 : MatVecMult(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 412128 : for(size_type i = 0; i < vOut.size(); ++i)
56 : {
57 274752 : vOut[i] = 0.0;
58 824256 : for(size_type j = 0; j < v.size(); ++j)
59 : {
60 549504 : vOut[i] += m(i,j) * v[j];
61 : }
62 : }
63 : }
64 :
65 : /// Matrix - Vector Multiplication adding to a second matrix
66 : // vOut += m * v
67 : template <typename vector_t_out, typename matrix_t, typename vector_t_in>
68 : inline
69 : void
70 : MatVecMultAppend(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 0 : for(size_type i = 0; i < vOut.size(); ++i)
77 : {
78 0 : for(size_type j = 0; j < v.size(); ++j)
79 : {
80 0 : vOut[i] += m(i,j) * v[j];
81 : }
82 : }
83 : }
84 :
85 : /// Matrix - Vector Multiplication added scaled to a second vector
86 : // vOut += s * m * v
87 : template <typename vector_t_out, typename matrix_t, typename vector_t_in>
88 : inline
89 : void
90 : MatVecScaleMultAppend(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 0 : for(size_type i = 0; i < vOut.size(); ++i)
97 : {
98 0 : for(size_type j = 0; j < v.size(); ++j)
99 : {
100 0 : vOut[i] += s * m(i,j) * v[j];
101 : }
102 : }
103 : }
104 :
105 :
106 : /// Transposed Matrix - Vector Muliplication
107 : // vOut = Transpose(m) * v
108 : template <typename vector_t_out, typename matrix_t, typename vector_t_in>
109 : inline
110 : void
111 : TransposedMatVecMult(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 0 : for(size_type i = 0; i < vOut.size(); ++i)
118 : {
119 0 : vOut[i] = 0.0;
120 0 : for(size_type j = 0; j < v.size(); ++j)
121 : {
122 0 : vOut[i] += m(j,i) * v[j];
123 : }
124 : }
125 : }
126 :
127 : /// Transposed Matrix - Vector Muliplication
128 : // vOut += Transpose(m) * v
129 : template <typename vector_t_out, typename matrix_t, typename vector_t_in>
130 : inline
131 : void
132 : TransposedMatVecMultAdd(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 :
147 : /**
148 : * Multiplication of a vector v by the Givens rotation transforming a given
149 : * matrix A to an upper triangular form R. After the call A stores the upper
150 : * triangular form R. Thus, this function computes the QR-decomposition
151 : * \f$ A = Q R \f$ with \f$ Q^T = Q^{-1} \f$ and performs the multiplication
152 : * \f$ v \gets Q^{-1} v \f$. Note that A is not necessarily a square matrix,
153 : * but it should have more lines than columns.
154 : */
155 : template <typename matrix_t, typename vector_t>
156 : inline
157 : void
158 0 : GivensMatVecMult (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 0 : for (size_type i = 0; i < matrix_t::RowSize - 1; i++) // the 1st index of the elementar rotation
170 0 : for (size_type k = matrix_t::RowSize - 1; k > i; k--) // the 2nd index of this rotation
171 : {
172 0 : d = A (i, i); x = A (k, i);
173 :
174 : // Computation of the entries of the elementar transformation:
175 0 : if (std::abs (d) > std::abs (x))
176 : {
177 0 : s = x / d;
178 0 : c = 1.0 / std::sqrt (1 + s * s);
179 0 : s *= c;
180 : }
181 0 : else if (x != 0) // to ensure that A (k, i) != 0; note that abs(x) >= abs(d) here!
182 : {
183 0 : c = d / x;
184 0 : s = 1.0 / std::sqrt (1 + c * c);
185 0 : c *= s;
186 : }
187 0 : else continue; // nothing to eliminate
188 :
189 : // Multiplication of A by the elementar transformation:
190 0 : for (size_type j = i; j < matrix_t::ColSize; j++)
191 : {
192 0 : x = A (i, j); y = A (k, j);
193 0 : A (i, j) = c * x + s * y;
194 0 : A (k, j) = c * y - s * x;
195 : }
196 :
197 : // Multiplication of v by the elementar transformation:
198 0 : x = v [i]; y = v [k];
199 0 : v [i] = c * x + s * y;
200 0 : v [k] = c * y - s * x;
201 : }
202 0 : }
203 :
204 : /**
205 : * Multiplication of a vector v by \f$ A^{-1} \f$ fr a given matrix A
206 : * (\f$ v := A^{-1} v \f$) using the QR-decomposition based on the Givens
207 : * rotations. Note that A is not necessarily a square matrix, but it should
208 : * have more lines than columns. In the latter case, the result is stored
209 : * in v[0]...v[A.num_cols() - 1], whereas the Euclidean norm of the rest
210 : * is the Euclidean distance between the original v and its projection to
211 : * the space spanned by the columns of the matrix. If the matrix is singular
212 : * (i.e. its columns are linearly dependent) then an exception is thrown.
213 : *
214 : * This function can be used in the least square method and for the orthogonal
215 : * projection of v to the space spanned by the columns of (the original) A.
216 : *
217 : * Remark: After the call, A stores the upper triangular form of the
218 : * QR-decomposition, i.e. the original matrix is destroyed.
219 : */
220 : template <typename matrix_t, typename vector_t>
221 : inline
222 : void
223 0 : InvMatVecMult_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 0 : 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 0 : i--;
235 0 : for (size_type j = i + 1; j < matrix_t::ColSize; j++)
236 0 : v [i] -= A (i, j) * v [j];
237 0 : if (std::abs (A (i, i)) < SMALL * std::abs (v [i]))
238 0 : UG_THROW ("InvMatVecMult_byGivens: Inverting singular matrix by the Givens rotations");
239 0 : v [i] /= A (i, i);
240 : }
241 0 : while (i != 0);
242 0 : }
243 :
244 : /**
245 : * Orthogonal projection of a given vector v to the space spanned by the columns
246 : * of a given matrix A. The projection is written to v.
247 : */
248 : template <typename matrix_t, typename vector_t>
249 : inline
250 : void
251 0 : OrthogProjectVec (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
258 0 : InvMatVecMult_byGivens (M, v);
259 :
260 : // II. Compute the linear combination of the columns:
261 : MathVector<matrix_t::ColSize, value_type> coeff;
262 0 : for (size_type i = 0; i < matrix_t::ColSize; i++) coeff [i] = v [i];
263 : MatVecMult (v, A, coeff);
264 0 : }
265 :
266 : } // end of namespace ug
267 :
268 : #endif /* __H__LGMATH__LGMATH_MATRIX_VECTOR_FUNCTIONS_COMMON_IMPL__ */
|
