17 #ifndef __deal2__householder_h
18 #define __deal2__householder_h
22 #include <deal.II/base/config.h>
23 #include <deal.II/lac/full_matrix.h>
24 #include <deal.II/lac/vector_memory.h>
32 template<
typename number>
class Vector;
55 template<
typename number>
73 template<
typename number2>
80 template<
typename number2>
101 template<
typename number2>
109 template<
typename number2>
118 template<
class VECTOR>
119 void vmult (VECTOR &dst,
const VECTOR &src)
const;
121 template<
class VECTOR>
122 void Tvmult (VECTOR &dst,
const VECTOR &src)
const;
141 template <
typename number>
147 template <
typename number>
148 template <
typename number2>
152 const size_type m = M.n_rows(), n = M.n_cols();
159 Assert (this->n_cols() <= this->n_rows(),
167 for (i=j ; i<m ; ++i)
168 sigma += this->el(i,j)*this->el(i,j);
171 if (std::fabs(sigma) < 1.e-15)
return;
173 number2 s = (this->el(j,j) < 0) ? std::sqrt(sigma) : -std::sqrt(sigma);
175 number2 beta = std::sqrt(1./(sigma-s*this->el(j,j)));
180 diagonal[j] = beta*(this->el(j,j) - s);
183 for (i=j+1 ; i<m ; ++i)
184 this->el(i,j) *= beta;
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);
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);
203 template <
typename number>
204 template <
typename number2>
211 template <
typename number>
212 template <
typename number2>
221 const size_type m = this->m(), n = this->n();
234 number2 sum = diagonal[j]* (*aux)(j);
236 sum += static_cast<number2>(this->el(i,j))*(*aux)(i);
238 (*aux)(j) -= sum*diagonal[j];
240 (*aux)(i) -= sum*static_cast<number2>(this->el(i,j));
245 sum += (*aux)(i) * (*aux)(i);
249 this->backward(dst, *aux);
253 return std::sqrt(sum);
256 template <
typename number>
257 template <
typename number2>
266 const size_type m = this->m(), n = this->n();
279 number2 sum = diagonal[j]* (*aux)(j);
281 sum += this->el(i,j)*(*aux)(i);
283 (*aux)(j) -= sum*diagonal[j];
285 (*aux)(i) -= sum*this->el(i,j);
290 sum += (*aux)(i) * (*aux)(i);
300 this->backward(v_dst, v_aux);
307 return std::sqrt(sum);
311 template <
typename number>
312 template <
class VECTOR>
316 least_squares (dst, src);
320 template <
typename number>
321 template <
class VECTOR>
332 DEAL_II_NAMESPACE_CLOSE
#define AssertDimension(dim1, dim2)
virtual void reinit(const size_type N, const bool fast=false)
bool is_finite(const double x)
types::global_dof_index size_type
void reinit(const unsigned int num_blocks, const size_type block_size=0, const bool fast=false)
T sum(const T &t, const MPI_Comm &mpi_communicator)
unsigned int global_dof_index
#define Assert(cond, exc)
void vmult(VECTOR &dst, const VECTOR &src) const
std::vector< number > diagonal
::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)
double least_squares(Vector< number2 > &dst, const Vector< number2 > &src) const