16 #ifndef dealii_householder_h 17 #define dealii_householder_h 20 #include <deal.II/base/config.h> 22 #include <deal.II/lac/full_matrix.h> 23 #include <deal.II/lac/vector_memory.h> 28 DEAL_II_NAMESPACE_OPEN
32 template <
typename number>
79 template <
typename number>
96 template <
typename number2>
104 template <
typename number2>
118 template <
typename number2>
120 least_squares(Vector<number2> &dst,
const Vector<number2> &src)
const;
125 template <
typename number2>
134 template <
class VectorType>
136 vmult(VectorType &dst,
const VectorType &src)
const;
142 template <
class VectorType>
144 Tvmult(VectorType &dst,
const VectorType &src)
const;
167 template <
typename number>
168 template <
typename number2>
172 const size_type m = M.n_rows(), n = M.n_cols();
187 for (i = j; i < m; ++i)
191 if (std::fabs(sigma) < 1.e-15)
194 number2 s = (
storage(j, j) < 0) ? std::sqrt(sigma) : -std::sqrt(sigma);
196 number2 beta = std::sqrt(1. / (sigma - s *
storage(j, j)));
204 for (i = j + 1; i < m; ++i)
213 for (i = j + 1; i < m; ++i)
217 for (i = j + 1; i < m; ++i)
225 template <
typename number>
226 template <
typename number2>
234 template <
typename number>
235 template <
typename number2>
238 const Vector<number2> &src)
const 248 aux->reinit(src,
true);
257 number2 sum =
diagonal[j] * (*aux)(j);
259 sum += static_cast<number2>(
storage(i, j)) * (*aux)(i);
263 (*aux)(i) -= sum * static_cast<number2>(
storage(i, j));
268 sum += (*aux)(i) * (*aux)(i);
274 return std::sqrt(sum);
279 template <
typename number>
280 template <
typename number2>
293 aux->reinit(src,
true);
302 number2 sum =
diagonal[j] * (*aux)(j);
304 sum +=
storage(i, j) * (*aux)(i);
308 (*aux)(i) -= sum *
storage(i, j);
313 sum += (*aux)(i) * (*aux)(i);
319 Vector<number2> v_dst, v_aux;
328 return std::sqrt(sum);
332 template <
typename number>
333 template <
class VectorType>
341 template <
typename number>
342 template <
class VectorType>
353 DEAL_II_NAMESPACE_CLOSE
static::ExceptionBase & ExcEmptyMatrix()
#define AssertDimension(dim1, dim2)
void backward(Vector< number2 > &dst, const Vector< number2 > &src) const
types::global_dof_index size_type
void Tvmult(VectorType &dst, const VectorType &src) const
unsigned long long int global_dof_index
void reinit(const TableIndices< N > &new_size, const bool omit_default_initialization=false)
#define Assert(cond, exc)
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
FullMatrix< double > storage
std::vector< number > diagonal
void initialize(const FullMatrix< number2 > &A)
void vmult(VectorType &dst, const VectorType &src) const
static::ExceptionBase & ExcNotImplemented()
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)