16 #ifndef dealii__householder_h 17 #define dealii__householder_h 21 #include <deal.II/base/config.h> 22 #include <deal.II/lac/full_matrix.h> 23 #include <deal.II/lac/vector_memory.h> 27 DEAL_II_NAMESPACE_OPEN
31 template<
typename number>
class Vector;
55 template<
typename number>
72 template<
typename number2>
78 template<
typename number2>
92 template<
typename number2>
99 template<
typename number2>
107 template<
class VectorType>
108 void vmult (VectorType &dst,
const VectorType &src)
const;
110 template<
class VectorType>
111 void Tvmult (VectorType &dst,
const VectorType &src)
const;
128 template <
typename number>
134 template <
typename number>
135 template <
typename number2>
146 Assert (this->n_cols() <= this->n_rows(),
154 for (i=j ; i<
m ; ++i)
155 sigma += this->
el(i,j)*this->
el(i,j);
158 if (std::fabs(sigma) < 1.e-15)
return;
160 number2 s = (this->
el(j,j) < 0) ? std::sqrt(sigma) : -std::sqrt(sigma);
162 number2 beta = std::sqrt(1./(sigma-s*this->
el(j,j)));
170 for (i=j+1 ; i<
m ; ++i)
171 this->
el(i,j) *= beta;
179 for (i=j+1 ; i<
m ; ++i)
180 sum += this->
el(i,j)*this->
el(i,k);
183 for (i=j+1 ; i<
m ; ++i)
184 this->
el(i,k) -= sum*this->
el(i,j);
190 template <
typename number>
191 template <
typename number2>
198 template <
typename number>
199 template <
typename number2>
221 number2 sum =
diagonal[j]* (*aux)(j);
223 sum += static_cast<number2>(this->
el(i,j))*(*aux)(i);
227 (*aux)(i) -= sum*static_cast<number2>(this->
el(i,j));
232 sum += (*aux)(i) * (*aux)(i);
240 return std::sqrt(sum);
243 template <
typename number>
244 template <
typename number2>
266 number2 sum =
diagonal[j]* (*aux)(j);
268 sum += this->
el(i,j)*(*aux)(i);
272 (*aux)(i) -= sum*this->
el(i,j);
277 sum += (*aux)(i) * (*aux)(i);
294 return std::sqrt(sum);
298 template <
typename number>
299 template <
class VectorType>
307 template <
typename number>
308 template <
class VectorType>
319 DEAL_II_NAMESPACE_CLOSE
#define AssertDimension(dim1, dim2)
void backward(Vector< number2 > &dst, const Vector< number2 > &src) const
types::global_dof_index size_type
unsigned int global_dof_index
void reinit(const TableIndices< N > &new_size, const bool omit_default_initialization=false)
virtual VectorType * alloc()
#define Assert(cond, exc)
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
AlignedVector< T >::reference el(const TableIndices< N > &indices)
virtual void free(const VectorType *const)
void reinit(const unsigned int n_blocks, const size_type block_size=0, const bool omit_zeroing_entries=false)
std::vector< number > diagonal
void vmult(VectorType &dst, const VectorType &src) const
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
static::ExceptionBase & ExcNotImplemented()
void initialize(const FullMatrix< number2 > &)
AlignedVector< T >::size_type size_type
double least_squares(Vector< number2 > &dst, const Vector< number2 > &src) const
#define AssertIsFinite(number)
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)