Reference documentation for deal.II version GIT relicensing-245-g36f19064f7 2024-03-29 07:20:02+00:00
\(\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// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2005 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_householder_h
16#define dealii_householder_h
17
18
19#include <deal.II/base/config.h>
20
23
24#include <cmath>
25#include <vector>
26
28
29
30// forward declarations
31#ifndef DOXYGEN
32template <typename number>
33class Vector;
34#endif
35
78template <typename number>
80{
81public:
86
90 Householder() = default;
91
95 template <typename number2>
97
103 template <typename number2>
104 void
106
117 template <typename number2>
118 double
120
124 template <typename number2>
125 double
127 const BlockVector<number2> &src) const;
128
133 template <typename VectorType>
134 void
135 vmult(VectorType &dst, const VectorType &src) const;
136
141 template <typename VectorType>
142 void
143 Tvmult(VectorType &dst, const VectorType &src) const;
144
145
146private:
151 std::vector<number> diagonal;
152
157};
158
161#ifndef DOXYGEN
162/*-------------------------Inline functions -------------------------------*/
163
164// QR-transformation cf. Stoer 1 4.8.2 (p. 191)
165
166template <typename number>
167template <typename number2>
168void
170{
171 const size_type m = M.n_rows(), n = M.n_cols();
172 storage.reinit(m, n);
173 storage.fill(M);
174 Assert(!storage.empty(), typename FullMatrix<number2>::ExcEmptyMatrix());
175 diagonal.resize(m);
176
177 // m > n, src.n() = m
178 Assert(storage.n_cols() <= storage.n_rows(),
179 ExcDimensionMismatch(storage.n_cols(), storage.n_rows()));
180
181 for (size_type j = 0; j < n; ++j)
182 {
183 number2 sigma = 0;
184 size_type i;
185 // sigma = ||v||^2
186 for (i = j; i < m; ++i)
187 sigma += storage(i, j) * storage(i, j);
188 // We are ready if the column is
189 // empty. Are we?
190 if (std::fabs(sigma) < 1.e-15)
191 return;
192
193 number2 s = (storage(j, j) < 0) ? std::sqrt(sigma) : -std::sqrt(sigma);
194 //
195 number2 beta = std::sqrt(1. / (sigma - s * storage(j, j)));
196
197 // Make column j the Householder
198 // vector, store first entry in
199 // diagonal
200 diagonal[j] = beta * (storage(j, j) - s);
201 storage(j, j) = s;
202
203 for (i = j + 1; i < m; ++i)
204 storage(i, j) *= beta;
205
206
207 // For all subsequent columns do
208 // the Householder reflection
209 for (size_type k = j + 1; k < n; ++k)
210 {
211 number2 sum = diagonal[j] * storage(j, k);
212 for (i = j + 1; i < m; ++i)
213 sum += storage(i, j) * storage(i, k);
214
215 storage(j, k) -= sum * this->diagonal[j];
216 for (i = j + 1; i < m; ++i)
217 storage(i, k) -= sum * storage(i, j);
218 }
219 }
220}
221
222
223
224template <typename number>
225template <typename number2>
227{
228 initialize(M);
229}
230
231
232
233template <typename number>
234template <typename number2>
235double
237 const Vector<number2> &src) const
238{
239 Assert(!storage.empty(), typename FullMatrix<number2>::ExcEmptyMatrix());
240 AssertDimension(dst.size(), storage.n());
241 AssertDimension(src.size(), storage.m());
242
243 const size_type m = storage.m(), n = storage.n();
244
246 typename VectorMemory<Vector<number2>>::Pointer aux(mem);
247 aux->reinit(src, true);
248 *aux = src;
249 // m > n, m = src.n, n = dst.n
250
251 // Multiply Q_n ... Q_2 Q_1 src
252 // Where Q_i = I - v_i v_i^T
253 for (size_type j = 0; j < n; ++j)
254 {
255 // sum = v_i^T dst
256 number2 sum = diagonal[j] * (*aux)(j);
257 for (size_type i = j + 1; i < m; ++i)
258 sum += static_cast<number2>(storage(i, j)) * (*aux)(i);
259 // dst -= v * sum
260 (*aux)(j) -= sum * diagonal[j];
261 for (size_type i = j + 1; i < m; ++i)
262 (*aux)(i) -= sum * static_cast<number2>(storage(i, j));
263 }
264 // Compute norm of residual
265 number2 sum = 0.;
266 for (size_type i = n; i < m; ++i)
267 sum += (*aux)(i) * (*aux)(i);
268 AssertIsFinite(sum);
269
270 // Compute solution
271 storage.backward(dst, *aux);
272
273 return std::sqrt(sum);
274}
275
276
277
278template <typename number>
279template <typename number2>
280double
282 const BlockVector<number2> &src) const
283{
284 Assert(!storage.empty(), typename FullMatrix<number2>::ExcEmptyMatrix());
285 AssertDimension(dst.size(), storage.n());
286 AssertDimension(src.size(), storage.m());
287
288 const size_type m = storage.m(), n = storage.n();
289
291 typename VectorMemory<BlockVector<number2>>::Pointer aux(mem);
292 aux->reinit(src, true);
293 *aux = src;
294 // m > n, m = src.n, n = dst.n
295
296 // Multiply Q_n ... Q_2 Q_1 src
297 // Where Q_i = I-v_i v_i^T
298 for (size_type j = 0; j < n; ++j)
299 {
300 // sum = v_i^T dst
301 number2 sum = diagonal[j] * (*aux)(j);
302 for (size_type i = j + 1; i < m; ++i)
303 sum += storage(i, j) * (*aux)(i);
304 // dst -= v * sum
305 (*aux)(j) -= sum * diagonal[j];
306 for (size_type i = j + 1; i < m; ++i)
307 (*aux)(i) -= sum * storage(i, j);
308 }
309 // Compute norm of residual
310 number2 sum = 0.;
311 for (size_type i = n; i < m; ++i)
312 sum += (*aux)(i) * (*aux)(i);
313 AssertIsFinite(sum);
314
315 // backward works for
316 // Vectors only, so copy
317 // them before
318 Vector<number2> v_dst, v_aux;
319 v_dst = dst;
320 v_aux = *aux;
321 // Compute solution
322 storage.backward(v_dst, v_aux);
323 // copy the result back
324 // to the BlockVector
325 dst = v_dst;
326
327 return std::sqrt(sum);
328}
329
330
331template <typename number>
332template <typename VectorType>
333void
334Householder<number>::vmult(VectorType &dst, const VectorType &src) const
335{
336 least_squares(dst, src);
337}
338
339
340template <typename number>
341template <typename VectorType>
342void
343Householder<number>::Tvmult(VectorType &, const VectorType &) const
344{
346}
347
348
349
350#endif // DOXYGEN
351
353
354#endif
virtual size_type size() const override
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)
virtual size_type size() const override
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
#define Assert(cond, exc)
#define AssertIsFinite(number)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define DEAL_II_NOT_IMPLEMENTED()
@ 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:81