Reference documentation for deal.II version 9.3.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\}}\)
householder.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2005 - 2020 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
33 template <typename number>
34 class Vector;
35 #endif
36 
78 template <typename number>
80 {
81 public:
86 
90  Householder() = default;
91 
95  template <typename number2>
97 
103  template <typename number2>
104  void
106 
117  template <typename number2>
118  double
119  least_squares(Vector<number2> &dst, const Vector<number2> &src) const;
120 
124  template <typename number2>
125  double
127  const BlockVector<number2> &src) const;
128 
133  template <class VectorType>
134  void
135  vmult(VectorType &dst, const VectorType &src) const;
136 
141  template <class VectorType>
142  void
143  Tvmult(VectorType &dst, const VectorType &src) const;
144 
145 
146 private:
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 
166 template <typename number>
167 template <typename number2>
168 void
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 
224 template <typename number>
225 template <typename number2>
227 {
228  initialize(M);
229 }
230 
231 
232 
233 template <typename number>
234 template <typename number2>
235 double
236 Householder<number>::least_squares(Vector<number2> & dst,
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 
278 template <typename number>
279 template <typename number2>
280 double
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 
331 template <typename number>
332 template <class VectorType>
333 void
334 Householder<number>::vmult(VectorType &dst, const VectorType &src) const
335 {
336  least_squares(dst, src);
337 }
338 
339 
340 template <typename number>
341 template <class VectorType>
342 void
344 {
345  Assert(false, ExcNotImplemented());
346 }
347 
348 
349 
350 #endif // DOXYGEN
351 
353 
354 #endif
size_type m() const
static ::ExceptionBase & ExcEmptyMatrix()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
void Tvmult(VectorType &dst, const VectorType &src) const
void backward(Vector< number2 > &dst, const Vector< number2 > &src) const
size_type n() const
T sum(const T &t, const MPI_Comm &mpi_communicator)
#define Assert(cond, exc)
Definition: exceptions.h:1465
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:395
void vmult(VectorType &dst, const VectorType &src) const
FullMatrix< double > storage
Definition: householder.h:156
Expression fabs(const Expression &x)
static const char A
unsigned int global_dof_index
Definition: types.h:76
std::vector< number > diagonal
Definition: householder.h:151
void initialize(const FullMatrix< number2 > &A)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:394
double least_squares(Vector< number2 > &dst, const Vector< number2 > &src) const
static ::ExceptionBase & ExcNotImplemented()
#define AssertIsFinite(number)
Definition: exceptions.h:1721
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)
Householder()=default