Reference documentation for deal.II version Git 840703ecf3 2019-10-14 17:15:20 -0400
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
householder.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2005 - 2018 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 
22 #include <deal.II/lac/full_matrix.h>
23 #include <deal.II/lac/vector_memory.h>
24 
25 #include <cmath>
26 #include <vector>
27 
28 DEAL_II_NAMESPACE_OPEN
29 
30 
31 // forward declarations
32 #ifndef DOXYGEN
33 template <typename number>
34 class Vector;
35 #endif
36 
80 template <typename number>
82 {
83 public:
88 
92  Householder() = default;
93 
97  template <typename number2>
99 
105  template <typename number2>
106  void
108 
119  template <typename number2>
120  double
121  least_squares(Vector<number2> &dst, const Vector<number2> &src) const;
122 
126  template <typename number2>
127  double
129  const BlockVector<number2> &src) const;
130 
135  template <class VectorType>
136  void
137  vmult(VectorType &dst, const VectorType &src) const;
138 
143  template <class VectorType>
144  void
145  Tvmult(VectorType &dst, const VectorType &src) const;
146 
147 
148 private:
153  std::vector<number> diagonal;
154 
159 };
160 
163 #ifndef DOXYGEN
164 /*-------------------------Inline functions -------------------------------*/
165 
166 // QR-transformation cf. Stoer 1 4.8.2 (p. 191)
167 
168 template <typename number>
169 template <typename number2>
170 void
172 {
173  const size_type m = M.n_rows(), n = M.n_cols();
174  storage.reinit(m, n);
175  storage.fill(M);
177  diagonal.resize(m);
178 
179  // m > n, src.n() = m
180  Assert(storage.n_cols() <= storage.n_rows(),
181  ExcDimensionMismatch(storage.n_cols(), storage.n_rows()));
182 
183  for (size_type j = 0; j < n; ++j)
184  {
185  number2 sigma = 0;
186  size_type i;
187  // sigma = ||v||^2
188  for (i = j; i < m; ++i)
189  sigma += storage(i, j) * storage(i, j);
190  // We are ready if the column is
191  // empty. Are we?
192  if (std::fabs(sigma) < 1.e-15)
193  return;
194 
195  number2 s = (storage(j, j) < 0) ? std::sqrt(sigma) : -std::sqrt(sigma);
196  //
197  number2 beta = std::sqrt(1. / (sigma - s * storage(j, j)));
198 
199  // Make column j the Householder
200  // vector, store first entry in
201  // diagonal
202  diagonal[j] = beta * (storage(j, j) - s);
203  storage(j, j) = s;
204 
205  for (i = j + 1; i < m; ++i)
206  storage(i, j) *= beta;
207 
208 
209  // For all subsequent columns do
210  // the Householder reflection
211  for (size_type k = j + 1; k < n; ++k)
212  {
213  number2 sum = diagonal[j] * storage(j, k);
214  for (i = j + 1; i < m; ++i)
215  sum += storage(i, j) * storage(i, k);
216 
217  storage(j, k) -= sum * this->diagonal[j];
218  for (i = j + 1; i < m; ++i)
219  storage(i, k) -= sum * storage(i, j);
220  }
221  }
222 }
223 
224 
225 
226 template <typename number>
227 template <typename number2>
229 {
230  initialize(M);
231 }
232 
233 
234 
235 template <typename number>
236 template <typename number2>
237 double
238 Householder<number>::least_squares(Vector<number2> & dst,
239  const Vector<number2> &src) const
240 {
242  AssertDimension(dst.size(), storage.n());
243  AssertDimension(src.size(), storage.m());
244 
245  const size_type m = storage.m(), n = storage.n();
246 
248  typename VectorMemory<Vector<number2>>::Pointer aux(mem);
249  aux->reinit(src, true);
250  *aux = src;
251  // m > n, m = src.n, n = dst.n
252 
253  // Multiply Q_n ... Q_2 Q_1 src
254  // Where Q_i = I - v_i v_i^T
255  for (size_type j = 0; j < n; ++j)
256  {
257  // sum = v_i^T dst
258  number2 sum = diagonal[j] * (*aux)(j);
259  for (size_type i = j + 1; i < m; ++i)
260  sum += static_cast<number2>(storage(i, j)) * (*aux)(i);
261  // dst -= v * sum
262  (*aux)(j) -= sum * diagonal[j];
263  for (size_type i = j + 1; i < m; ++i)
264  (*aux)(i) -= sum * static_cast<number2>(storage(i, j));
265  }
266  // Compute norm of residual
267  number2 sum = 0.;
268  for (size_type i = n; i < m; ++i)
269  sum += (*aux)(i) * (*aux)(i);
270  AssertIsFinite(sum);
271 
272  // Compute solution
273  storage.backward(dst, *aux);
274 
275  return std::sqrt(sum);
276 }
277 
278 
279 
280 template <typename number>
281 template <typename number2>
282 double
284  const BlockVector<number2> &src) const
285 {
287  AssertDimension(dst.size(), storage.n());
288  AssertDimension(src.size(), storage.m());
289 
290  const size_type m = storage.m(), n = storage.n();
291 
293  typename VectorMemory<BlockVector<number2>>::Pointer aux(mem);
294  aux->reinit(src, true);
295  *aux = src;
296  // m > n, m = src.n, n = dst.n
297 
298  // Multiply Q_n ... Q_2 Q_1 src
299  // Where Q_i = I-v_i v_i^T
300  for (size_type j = 0; j < n; ++j)
301  {
302  // sum = v_i^T dst
303  number2 sum = diagonal[j] * (*aux)(j);
304  for (size_type i = j + 1; i < m; ++i)
305  sum += storage(i, j) * (*aux)(i);
306  // dst -= v * sum
307  (*aux)(j) -= sum * diagonal[j];
308  for (size_type i = j + 1; i < m; ++i)
309  (*aux)(i) -= sum * storage(i, j);
310  }
311  // Compute norm of residual
312  number2 sum = 0.;
313  for (size_type i = n; i < m; ++i)
314  sum += (*aux)(i) * (*aux)(i);
315  AssertIsFinite(sum);
316 
317  // backward works for
318  // Vectors only, so copy
319  // them before
320  Vector<number2> v_dst, v_aux;
321  v_dst = dst;
322  v_aux = *aux;
323  // Compute solution
324  storage.backward(v_dst, v_aux);
325  // copy the result back
326  // to the BlockVector
327  dst = v_dst;
328 
329  return std::sqrt(sum);
330 }
331 
332 
333 template <typename number>
334 template <class VectorType>
335 void
336 Householder<number>::vmult(VectorType &dst, const VectorType &src) const
337 {
338  least_squares(dst, src);
339 }
340 
341 
342 template <typename number>
343 template <class VectorType>
344 void
345 Householder<number>::Tvmult(VectorType &, const VectorType &) const
346 {
347  Assert(false, ExcNotImplemented());
348 }
349 
350 
351 
352 #endif // DOXYGEN
353 
354 DEAL_II_NAMESPACE_CLOSE
355 
356 #endif
size_type m() const
static ::ExceptionBase & ExcEmptyMatrix()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1567
void Tvmult(VectorType &dst, const VectorType &src) const
types::global_dof_index size_type
Definition: householder.h:87
LinearAlgebra::distributed::Vector< Number > Vector
void backward(Vector< number2 > &dst, const Vector< number2 > &src) const
size_type n() const
void reinit(const TableIndices< N > &new_size, const bool omit_default_initialization=false)
#define Assert(cond, exc)
Definition: exceptions.h:1407
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
void vmult(VectorType &dst, const VectorType &src) const
FullMatrix< double > storage
Definition: householder.h:158
unsigned int global_dof_index
Definition: types.h:89
std::vector< number > diagonal
Definition: householder.h:153
void initialize(const FullMatrix< number2 > &A)
double least_squares(Vector< number2 > &dst, const Vector< number2 > &src) const
static ::ExceptionBase & ExcNotImplemented()
bool empty() const
#define AssertIsFinite(number)
Definition: exceptions.h:1669
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