16 #ifndef dealii_householder_h
17 #define dealii_householder_h
33 template <
typename number>
80 template <
typename number>
97 template <
typename number2>
105 template <
typename number2>
119 template <
typename number2>
121 least_squares(Vector<number2> &dst,
const Vector<number2> &src)
const;
126 template <
typename number2>
135 template <
class VectorType>
143 template <
class VectorType>
168 template <
typename number>
169 template <
typename number2>
173 const size_type m = M.n_rows(), n = M.n_cols();
180 Assert(storage.n_cols() <= storage.n_rows(),
188 for (i = j; i < m; ++i)
189 sigma += storage(i, j) * storage(i, j);
195 number2 s = (storage(j, j) < 0) ? std::sqrt(sigma) : -std::sqrt(sigma);
197 number2 beta = std::sqrt(1. / (sigma - s * storage(j, j)));
202 diagonal[j] = beta * (storage(j, j) - s);
205 for (i = j + 1; i < m; ++i)
206 storage(i, j) *= beta;
214 for (i = j + 1; i < m; ++i)
215 sum += storage(i, j) * storage(i, k);
218 for (i = j + 1; i < m; ++i)
219 storage(i, k) -=
sum * storage(i, j);
226 template <
typename number>
227 template <
typename number2>
235 template <
typename number>
236 template <
typename number2>
239 const Vector<number2> &src)
const
245 const size_type m = storage.m(), n = storage.n();
249 aux->reinit(src,
true);
260 sum +=
static_cast<number2
>(storage(i, j)) * (*aux)(i);
264 (*aux)(i) -=
sum *
static_cast<number2
>(storage(i, j));
269 sum += (*aux)(i) * (*aux)(i);
273 storage.backward(dst, *aux);
280 template <
typename number>
281 template <
typename number2>
290 const size_type m = storage.m(), n = storage.n();
294 aux->reinit(src,
true);
305 sum += storage(i, j) * (*aux)(i);
309 (*aux)(i) -=
sum * storage(i, j);
314 sum += (*aux)(i) * (*aux)(i);
320 Vector<number2> v_dst, v_aux;
324 storage.backward(v_dst, v_aux);
333 template <
typename number>
334 template <
class VectorType>
338 least_squares(dst, src);
342 template <
typename number>
343 template <
class VectorType>