15#ifndef dealii_householder_h
16#define dealii_householder_h
32template <
typename number>
78template <
typename number>
95 template <
typename number2>
103 template <
typename number2>
117 template <
typename number2>
124 template <
typename number2>
133 template <
typename VectorType>
135 vmult(VectorType &dst,
const VectorType &src)
const;
141 template <
typename VectorType>
143 Tvmult(VectorType &dst,
const VectorType &src)
const;
166template <
typename number>
167template <
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(),
181 for (size_type
j = 0;
j < n; ++
j)
186 for (i =
j; i < m; ++i)
187 sigma += storage(i,
j) * storage(i,
j);
190 if (std::fabs(
sigma) < 1.e-15)
200 diagonal[
j] =
beta * (storage(
j,
j) - s);
203 for (i =
j + 1; i < m; ++i)
204 storage(i,
j) *=
beta;
209 for (size_type
k =
j + 1;
k < n; ++
k)
212 for (i =
j + 1; i < m; ++i)
213 sum += storage(i,
j) * storage(i,
k);
215 storage(
j,
k) -=
sum * this->diagonal[
j];
216 for (i =
j + 1; i < m; ++i)
217 storage(i,
k) -=
sum * storage(i,
j);
224template <
typename number>
225template <
typename number2>
233template <
typename number>
234template <
typename number2>
243 const size_type m = storage.m(), n = storage.n();
247 aux->reinit(src,
true);
253 for (size_type
j = 0;
j < n; ++
j)
257 for (size_type i =
j + 1; i < m; ++i)
258 sum +=
static_cast<number2>(storage(i,
j)) * (*aux)(i);
260 (*aux)(
j) -=
sum * diagonal[
j];
261 for (size_type i =
j + 1; i < m; ++i)
262 (*aux)(i) -=
sum *
static_cast<number2>(storage(i,
j));
266 for (size_type i = n; i < m; ++i)
267 sum += (*aux)(i) * (*aux)(i);
271 storage.backward(dst, *aux);
278template <
typename number>
279template <
typename number2>
288 const size_type m = storage.m(), n = storage.n();
292 aux->reinit(src,
true);
298 for (size_type
j = 0;
j < n; ++
j)
302 for (size_type i =
j + 1; i < m; ++i)
303 sum += storage(i,
j) * (*aux)(i);
305 (*aux)(
j) -=
sum * diagonal[
j];
306 for (size_type i =
j + 1; i < m; ++i)
307 (*aux)(i) -=
sum * storage(i,
j);
311 for (size_type i = n; i < m; ++i)
312 sum += (*aux)(i) * (*aux)(i);
331template <
typename number>
332template <
typename VectorType>
336 least_squares(dst, src);
340template <
typename number>
341template <
typename VectorType>
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
FullMatrix< double > storage
Householder(const FullMatrix< number2 > &A)
void initialize(const FullMatrix< number2 > &A)
Tensor< rank, dim, Number > sum(const Tensor< rank, dim, Number > &local, const MPI_Comm mpi_communicator)
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
#define Assert(cond, exc)
#define AssertIsFinite(number)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define DEAL_II_NOT_IMPLEMENTED()
types::global_dof_index size_type
@ diagonal
Matrix is diagonal.
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
unsigned int global_dof_index