16 #ifndef dealii_householder_h
17 #define dealii_householder_h
33 template <
typename number>
78 template <
typename number>
95 template <
typename number2>
103 template <
typename number2>
117 template <
typename number2>
124 template <
typename number2>
133 template <
class VectorType>
135 vmult(VectorType &dst,
const VectorType &src)
const;
141 template <
class VectorType>
143 Tvmult(VectorType &dst,
const VectorType &src)
const;
166 template <
typename number>
167 template <
typename number2>
171 const size_type m = M.n_rows(), n = M.n_cols();
172 storage.reinit(m, n);
178 Assert(storage.n_cols() <= storage.n_rows(),
186 for (i = j; i < m; ++i)
187 sigma += storage(i, j) * storage(i, j);
193 number2 s = (storage(j, j) < 0) ? std::sqrt(sigma) : -std::sqrt(sigma);
195 number2 beta = std::sqrt(1. / (sigma - s * storage(j, j)));
200 diagonal[j] = beta * (storage(j, j) - s);
203 for (i = j + 1; i < m; ++i)
204 storage(i, j) *= beta;
212 for (i = j + 1; i < m; ++i)
213 sum += storage(i, j) * storage(i, k);
216 for (i = j + 1; i < m; ++i)
217 storage(i, k) -=
sum * storage(i, j);
224 template <
typename number>
225 template <
typename number2>
233 template <
typename number>
234 template <
typename number2>
243 const size_type m = storage.m(), n = storage.n();
247 aux->reinit(src,
true);
258 sum +=
static_cast<number2
>(storage(i, j)) * (*aux)(i);
262 (*aux)(i) -=
sum *
static_cast<number2
>(storage(i, j));
267 sum += (*aux)(i) * (*aux)(i);
271 storage.backward(dst, *aux);
278 template <
typename number>
279 template <
typename number2>
288 const size_type m = storage.m(), n = storage.n();
292 aux->reinit(src,
true);
303 sum += storage(i, j) * (*aux)(i);
307 (*aux)(i) -=
sum * storage(i, j);
312 sum += (*aux)(i) * (*aux)(i);
322 storage.backward(v_dst, v_aux);
331 template <
typename number>
332 template <
class VectorType>
336 least_squares(dst, src);
340 template <
typename number>
341 template <
class VectorType>
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)
SymmetricTensor< rank, dim, Number > sum(const SymmetricTensor< rank, dim, Number > &local, const MPI_Comm &mpi_communicator)
#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)
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)
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 > &)