Reference documentation for deal.II version 9.5.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
householder.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2005 - 2022 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16#ifndef dealii_householder_h
17#define dealii_householder_h
18
19
20#include <deal.II/base/config.h>
21
24
25#include <cmath>
26#include <vector>
27
29
30
31// forward declarations
32#ifndef DOXYGEN
33template <typename number>
34class Vector;
35#endif
36
79template <typename number>
81{
82public:
87
91 Householder() = default;
92
96 template <typename number2>
98
104 template <typename number2>
105 void
107
118 template <typename number2>
119 double
121
125 template <typename number2>
126 double
128 const BlockVector<number2> &src) const;
129
134 template <class VectorType>
135 void
136 vmult(VectorType &dst, const VectorType &src) const;
137
142 template <class VectorType>
143 void
144 Tvmult(VectorType &dst, const VectorType &src) const;
145
146
147private:
152 std::vector<number> diagonal;
153
158};
159
162#ifndef DOXYGEN
163/*-------------------------Inline functions -------------------------------*/
164
165// QR-transformation cf. Stoer 1 4.8.2 (p. 191)
166
167template <typename number>
168template <typename number2>
169void
171{
172 const size_type m = M.n_rows(), n = M.n_cols();
173 storage.reinit(m, n);
174 storage.fill(M);
175 Assert(!storage.empty(), typename FullMatrix<number2>::ExcEmptyMatrix());
176 diagonal.resize(m);
177
178 // m > n, src.n() = m
179 Assert(storage.n_cols() <= storage.n_rows(),
180 ExcDimensionMismatch(storage.n_cols(), storage.n_rows()));
181
182 for (size_type j = 0; j < n; ++j)
183 {
184 number2 sigma = 0;
185 size_type i;
186 // sigma = ||v||^2
187 for (i = j; i < m; ++i)
188 sigma += storage(i, j) * storage(i, j);
189 // We are ready if the column is
190 // empty. Are we?
191 if (std::fabs(sigma) < 1.e-15)
192 return;
193
194 number2 s = (storage(j, j) < 0) ? std::sqrt(sigma) : -std::sqrt(sigma);
195 //
196 number2 beta = std::sqrt(1. / (sigma - s * storage(j, j)));
197
198 // Make column j the Householder
199 // vector, store first entry in
200 // diagonal
201 diagonal[j] = beta * (storage(j, j) - s);
202 storage(j, j) = s;
203
204 for (i = j + 1; i < m; ++i)
205 storage(i, j) *= beta;
206
207
208 // For all subsequent columns do
209 // the Householder reflection
210 for (size_type k = j + 1; k < n; ++k)
211 {
212 number2 sum = diagonal[j] * storage(j, k);
213 for (i = j + 1; i < m; ++i)
214 sum += storage(i, j) * storage(i, k);
215
216 storage(j, k) -= sum * this->diagonal[j];
217 for (i = j + 1; i < m; ++i)
218 storage(i, k) -= sum * storage(i, j);
219 }
220 }
221}
222
223
224
225template <typename number>
226template <typename number2>
228{
229 initialize(M);
230}
231
232
233
234template <typename number>
235template <typename number2>
236double
238 const Vector<number2> &src) const
239{
240 Assert(!storage.empty(), typename FullMatrix<number2>::ExcEmptyMatrix());
241 AssertDimension(dst.size(), storage.n());
242 AssertDimension(src.size(), storage.m());
243
244 const size_type m = storage.m(), n = storage.n();
245
247 typename VectorMemory<Vector<number2>>::Pointer aux(mem);
248 aux->reinit(src, true);
249 *aux = src;
250 // m > n, m = src.n, n = dst.n
251
252 // Multiply Q_n ... Q_2 Q_1 src
253 // Where Q_i = I - v_i v_i^T
254 for (size_type j = 0; j < n; ++j)
255 {
256 // sum = v_i^T dst
257 number2 sum = diagonal[j] * (*aux)(j);
258 for (size_type i = j + 1; i < m; ++i)
259 sum += static_cast<number2>(storage(i, j)) * (*aux)(i);
260 // dst -= v * sum
261 (*aux)(j) -= sum * diagonal[j];
262 for (size_type i = j + 1; i < m; ++i)
263 (*aux)(i) -= sum * static_cast<number2>(storage(i, j));
264 }
265 // Compute norm of residual
266 number2 sum = 0.;
267 for (size_type i = n; i < m; ++i)
268 sum += (*aux)(i) * (*aux)(i);
269 AssertIsFinite(sum);
270
271 // Compute solution
272 storage.backward(dst, *aux);
273
274 return std::sqrt(sum);
275}
276
277
278
279template <typename number>
280template <typename number2>
281double
283 const BlockVector<number2> &src) const
284{
285 Assert(!storage.empty(), typename FullMatrix<number2>::ExcEmptyMatrix());
286 AssertDimension(dst.size(), storage.n());
287 AssertDimension(src.size(), storage.m());
288
289 const size_type m = storage.m(), n = storage.n();
290
292 typename VectorMemory<BlockVector<number2>>::Pointer aux(mem);
293 aux->reinit(src, true);
294 *aux = src;
295 // m > n, m = src.n, n = dst.n
296
297 // Multiply Q_n ... Q_2 Q_1 src
298 // Where Q_i = I-v_i v_i^T
299 for (size_type j = 0; j < n; ++j)
300 {
301 // sum = v_i^T dst
302 number2 sum = diagonal[j] * (*aux)(j);
303 for (size_type i = j + 1; i < m; ++i)
304 sum += storage(i, j) * (*aux)(i);
305 // dst -= v * sum
306 (*aux)(j) -= sum * diagonal[j];
307 for (size_type i = j + 1; i < m; ++i)
308 (*aux)(i) -= sum * storage(i, j);
309 }
310 // Compute norm of residual
311 number2 sum = 0.;
312 for (size_type i = n; i < m; ++i)
313 sum += (*aux)(i) * (*aux)(i);
314 AssertIsFinite(sum);
315
316 // backward works for
317 // Vectors only, so copy
318 // them before
319 Vector<number2> v_dst, v_aux;
320 v_dst = dst;
321 v_aux = *aux;
322 // Compute solution
323 storage.backward(v_dst, v_aux);
324 // copy the result back
325 // to the BlockVector
326 dst = v_dst;
327
328 return std::sqrt(sum);
329}
330
331
332template <typename number>
333template <class VectorType>
334void
335Householder<number>::vmult(VectorType &dst, const VectorType &src) const
336{
337 least_squares(dst, src);
338}
339
340
341template <typename number>
342template <class VectorType>
343void
344Householder<number>::Tvmult(VectorType &, const VectorType &) const
345{
346 Assert(false, ExcNotImplemented());
347}
348
349
350
351#endif // DOXYGEN
352
354
355#endif
std::size_t size() const
void fill(const FullMatrix< number2 > &src, const size_type dst_offset_i=0, const size_type dst_offset_j=0, const size_type src_offset_i=0, const size_type src_offset_j=0)
void vmult(VectorType &dst, const VectorType &src) const
double least_squares(BlockVector< number2 > &dst, const BlockVector< number2 > &src) const
void Tvmult(VectorType &dst, const VectorType &src) const
double least_squares(Vector< number2 > &dst, const Vector< number2 > &src) const
std::vector< number > diagonal
FullMatrix< double > storage
Householder()=default
Householder(const FullMatrix< number2 > &A)
void initialize(const FullMatrix< number2 > &A)
size_type size() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertIsFinite(number)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
@ diagonal
Matrix is diagonal.
types::global_dof_index size_type
T sum(const T &t, const MPI_Comm mpi_communicator)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
unsigned int global_dof_index
Definition types.h:82