Reference documentation for deal.II version 8.1.0
householder.h
1 // ---------------------------------------------------------------------
2 // @f$Id: householder.h 30036 2013-07-18 16:55:32Z maier @f$
3 //
4 // Copyright (C) 2005 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 #ifndef __deal2__householder_h
18 #define __deal2__householder_h
19 
20 
21 #include <cmath>
22 #include <deal.II/base/config.h>
23 #include <deal.II/lac/full_matrix.h>
24 #include <deal.II/lac/vector_memory.h>
25 
26 #include <vector>
27 
29 
30 
31 // forward declarations
32 template<typename number> class Vector;
33 
34 
55 template<typename number>
56 class Householder : private FullMatrix<number>
57 {
58 public:
63 
67  Householder ();
68 
73  template<typename number2>
75 
80  template<typename number2>
81  void
83 
101  template<typename number2>
102  double least_squares (Vector<number2> &dst,
103  const Vector<number2> &src) const;
104 
109  template<typename number2>
110  double least_squares (BlockVector<number2> &dst,
111  const BlockVector<number2> &src) const;
112 
118  template<class VECTOR>
119  void vmult (VECTOR &dst, const VECTOR &src) const;
120 
121  template<class VECTOR>
122  void Tvmult (VECTOR &dst, const VECTOR &src) const;
123 
124 
125 private:
131  std::vector<number> diagonal;
132 };
133 
136 #ifndef DOXYGEN
137 /*-------------------------Inline functions -------------------------------*/
138 
139 // QR-transformation cf. Stoer 1 4.8.2 (p. 191)
140 
141 template <typename number>
143 {}
144 
145 
146 
147 template <typename number>
148 template <typename number2>
149 void
151 {
152  const size_type m = M.n_rows(), n = M.n_cols();
153  this->reinit(m, n);
154  this->fill(M);
155  Assert (!this->empty(), typename FullMatrix<number2>::ExcEmptyMatrix());
156  diagonal.resize(m);
157 
158  // m > n, src.n() = m
159  Assert (this->n_cols() <= this->n_rows(),
160  ExcDimensionMismatch(this->n_cols(), this->n_rows()));
161 
162  for (size_type j=0 ; j<n ; ++j)
163  {
164  number2 sigma = 0;
165  size_type i;
166  // sigma = ||v||^2
167  for (i=j ; i<m ; ++i)
168  sigma += this->el(i,j)*this->el(i,j);
169  // We are ready if the column is
170  // empty. Are we?
171  if (std::fabs(sigma) < 1.e-15) return;
172 
173  number2 s = (this->el(j,j) < 0) ? std::sqrt(sigma) : -std::sqrt(sigma);
174  //
175  number2 beta = std::sqrt(1./(sigma-s*this->el(j,j)));
176 
177  // Make column j the Householder
178  // vector, store first entry in
179  // diagonal
180  diagonal[j] = beta*(this->el(j,j) - s);
181  this->el(j,j) = s;
182 
183  for (i=j+1 ; i<m ; ++i)
184  this->el(i,j) *= beta;
185 
186 
187  // For all subsequent columns do
188  // the Householder reflexion
189  for (size_type k=j+1 ; k<n ; ++k)
190  {
191  number2 sum = diagonal[j]*this->el(j,k);
192  for (i=j+1 ; i<m ; ++i)
193  sum += this->el(i,j)*this->el(i,k);
194 
195  this->el(j,k) -= sum*this->diagonal[j];
196  for (i=j+1 ; i<m ; ++i)
197  this->el(i,k) -= sum*this->el(i,j);
198  }
199  }
200 }
201 
202 
203 template <typename number>
204 template <typename number2>
206 {
207  initialize(M);
208 }
209 
210 
211 template <typename number>
212 template <typename number2>
213 double
215  const Vector<number2> &src) const
216 {
217  Assert (!this->empty(), typename FullMatrix<number2>::ExcEmptyMatrix());
218  AssertDimension(dst.size(), this->n());
219  AssertDimension(src.size(), this->m());
220 
221  const size_type m = this->m(), n = this->n();
222 
224  Vector<number2> *aux = mem.alloc();
225  aux->reinit(src, true);
226  *aux = src;
227  // m > n, m = src.n, n = dst.n
228 
229  // Multiply Q_n ... Q_2 Q_1 src
230  // Where Q_i = I-v_i v_i^T
231  for (size_type j=0; j<n; ++j)
232  {
233  // sum = v_i^T dst
234  number2 sum = diagonal[j]* (*aux)(j);
235  for (size_type i=j+1 ; i<m ; ++i)
236  sum += static_cast<number2>(this->el(i,j))*(*aux)(i);
237  // dst -= v * sum
238  (*aux)(j) -= sum*diagonal[j];
239  for (size_type i=j+1 ; i<m ; ++i)
240  (*aux)(i) -= sum*static_cast<number2>(this->el(i,j));
241  }
242  // Compute norm of residual
243  number2 sum = 0.;
244  for (size_type i=n ; i<m ; ++i)
245  sum += (*aux)(i) * (*aux)(i);
247 
248  // Compute solution
249  this->backward(dst, *aux);
250 
251  mem.free(aux);
252 
253  return std::sqrt(sum);
254 }
255 
256 template <typename number>
257 template <typename number2>
258 double
260  const BlockVector<number2> &src) const
261 {
262  Assert (!this->empty(), typename FullMatrix<number2>::ExcEmptyMatrix());
263  AssertDimension(dst.size(), this->n());
264  AssertDimension(src.size(), this->m());
265 
266  const size_type m = this->m(), n = this->n();
267 
269  BlockVector<number2> *aux = mem.alloc();
270  aux->reinit(src, true);
271  *aux = src;
272  // m > n, m = src.n, n = dst.n
273 
274  // Multiply Q_n ... Q_2 Q_1 src
275  // Where Q_i = I-v_i v_i^T
276  for (size_type j=0; j<n; ++j)
277  {
278  // sum = v_i^T dst
279  number2 sum = diagonal[j]* (*aux)(j);
280  for (size_type i=j+1 ; i<m ; ++i)
281  sum += this->el(i,j)*(*aux)(i);
282  // dst -= v * sum
283  (*aux)(j) -= sum*diagonal[j];
284  for (size_type i=j+1 ; i<m ; ++i)
285  (*aux)(i) -= sum*this->el(i,j);
286  }
287  // Compute norm of residual
288  number2 sum = 0.;
289  for (size_type i=n ; i<m ; ++i)
290  sum += (*aux)(i) * (*aux)(i);
292 
293  //backward works for
294  //Vectors only, so copy
295  //them before
296  Vector<number2> v_dst, v_aux;
297  v_dst = dst;
298  v_aux = *aux;
299  // Compute solution
300  this->backward(v_dst, v_aux);
301  //copy the result back
302  //to the BlockVector
303  dst = v_dst;
304 
305  mem.free(aux);
306 
307  return std::sqrt(sum);
308 }
309 
310 
311 template <typename number>
312 template <class VECTOR>
313 void
314 Householder<number>::vmult (VECTOR &dst, const VECTOR &src) const
315 {
316  least_squares (dst, src);
317 }
318 
319 
320 template <typename number>
321 template <class VECTOR>
322 void
323 Householder<number>::Tvmult (VECTOR &, const VECTOR &) const
324 {
325  Assert(false, ExcNotImplemented());
326 }
327 
328 
329 
330 #endif // DOXYGEN
331 
332 DEAL_II_NAMESPACE_CLOSE
333 
334 #endif
335 
void reinit(const TableIndices< N > &new_size, const bool fast=false)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:858
void backward(Vector< number2 > &dst, const Vector< number2 > &src) const
virtual VECTOR * alloc()
bool is_finite(const double x)
void reinit(const unsigned int num_blocks, const size_type block_size=0, const bool fast=false)
types::global_dof_index size_type
Definition: householder.h:62
size_type n() const
unsigned int global_dof_index
Definition: types.h:100
virtual void reinit(const size_type N, const bool fast=false)
#define Assert(cond, exc)
Definition: exceptions.h:299
std::size_t size() const
void vmult(VECTOR &dst, const VECTOR &src) const
size_type m() const
std::vector< number > diagonal
Definition: householder.h:131
::ExceptionBase & ExcNumberNotFinite()
::ExceptionBase & ExcNotImplemented()
::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
void initialize(const FullMatrix< number2 > &)
virtual void free(const VECTOR *const)
std::vector< number >::reference el(const TableIndices< N > &indices)
double least_squares(Vector< number2 > &dst, const Vector< number2 > &src) 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)