Reference documentation for deal.II version 9.1.1
\(\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 template <typename number>
33 class Vector;
34 
35 
79 template <typename number>
81 {
82 public:
87 
91  Householder() = default;
92 
96  template <typename number2>
98 
104  template <typename number2>
105  void
107 
118  template <typename number2>
119  double
120  least_squares(Vector<number2> &dst, const Vector<number2> &src) const;
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 
147 private:
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 
167 template <typename number>
168 template <typename number2>
169 void
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 
225 template <typename number>
226 template <typename number2>
228 {
229  initialize(M);
230 }
231 
232 
233 
234 template <typename number>
235 template <typename number2>
236 double
237 Householder<number>::least_squares(Vector<number2> & dst,
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);
270 
271  // Compute solution
272  storage.backward(dst, *aux);
273 
274  return std::sqrt(sum);
275 }
276 
277 
278 
279 template <typename number>
280 template <typename number2>
281 double
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);
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 
332 template <typename number>
333 template <class VectorType>
334 void
335 Householder<number>::vmult(VectorType &dst, const VectorType &src) const
336 {
337  least_squares(dst, src);
338 }
339 
340 
341 template <typename number>
342 template <class VectorType>
343 void
344 Householder<number>::Tvmult(VectorType &, const VectorType &) const
345 {
346  Assert(false, ExcNotImplemented());
347 }
348 
349 
350 
351 #endif // DOXYGEN
352 
353 DEAL_II_NAMESPACE_CLOSE
354 
355 #endif
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:86
LinearAlgebra::distributed::Vector< Number > Vector
Matrix is diagonal.
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:157
SymmetricTensor< rank, dim, Number > sum(const SymmetricTensor< rank, dim, Number > &local, const MPI_Comm &mpi_communicator)
unsigned int global_dof_index
Definition: types.h:89
std::vector< number > diagonal
Definition: householder.h:152
void initialize(const FullMatrix< number2 > &A)
double least_squares(Vector< number2 > &dst, const Vector< number2 > &src) const
static ::ExceptionBase & ExcNotImplemented()
#define AssertIsFinite(number)
Definition: exceptions.h:1669
Householder()=default