16 #ifndef dealii_householder_h
17 #define dealii_householder_h
33 template <
typename number>
79 template <
typename number>
96 template <
typename number2>
104 template <
typename number2>
118 template <
typename number2>
125 template <
typename number2>
134 template <
typename VectorType>
136 vmult(VectorType &dst,
const VectorType &src)
const;
142 template <
typename 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();
173 storage.reinit(m, n);
179 Assert(storage.n_cols() <= storage.n_rows(),
187 for (i = j; i < m; ++i)
188 sigma += storage(i, j) * storage(i, j);
194 number2 s = (storage(j, j) < 0) ? std::sqrt(sigma) : -std::sqrt(sigma);
196 number2 beta = std::sqrt(1. / (sigma - s * storage(j, j)));
201 diagonal[j] = beta * (storage(j, j) - s);
204 for (i = j + 1; i < m; ++i)
205 storage(i, j) *= beta;
213 for (i = j + 1; i < m; ++i)
214 sum += storage(i, j) * storage(i, k);
217 for (i = j + 1; i < m; ++i)
218 storage(i, k) -=
sum * storage(i, j);
225 template <
typename number>
226 template <
typename number2>
234 template <
typename number>
235 template <
typename number2>
244 const size_type m = storage.m(), n = storage.n();
248 aux->reinit(src,
true);
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);
272 storage.backward(dst, *aux);
279 template <
typename number>
280 template <
typename number2>
289 const size_type m = storage.m(), n = storage.n();
293 aux->reinit(src,
true);
304 sum += storage(i, j) * (*aux)(i);
308 (*aux)(i) -=
sum * storage(i, j);
313 sum += (*aux)(i) * (*aux)(i);
323 storage.backward(v_dst, v_aux);
332 template <
typename number>
333 template <
typename VectorType>
337 least_squares(dst, src);
341 template <
typename number>
342 template <
typename VectorType>
virtual size_type size() const override
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)
void vmult(VectorType &dst, const VectorType &src) const
double least_squares(BlockVector< number2 > &dst, const BlockVector< number2 > &src) const
void Tvmult(VectorType &dst, const VectorType &src) const
double least_squares(Vector< number2 > &dst, const Vector< number2 > &src) const
std::vector< number > diagonal
types::global_dof_index size_type
FullMatrix< double > storage
Householder(const FullMatrix< number2 > &A)
void initialize(const FullMatrix< number2 > &A)
SymmetricTensor< rank, dim, Number > sum(const SymmetricTensor< rank, dim, Number > &local, const MPI_Comm mpi_communicator)
virtual size_type size() const override
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
static ::ExceptionBase & ExcNotImplemented()
#define AssertIsFinite(number)
#define AssertDimension(dim1, dim2)
Expression fabs(const Expression &x)
@ diagonal
Matrix is diagonal.
types::global_dof_index size_type
unsigned int global_dof_index
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)