* return.
*/
template <typename number2>
- double
+ number2
least_squares(Vector<number2> &dst, const Vector<number2> &src) const;
/**
/**
* Storage that is internally used for the Householder transformation.
*/
- FullMatrix<double> storage;
+ FullMatrix<number> storage;
};
/** @} */
if (std::fabs(sigma) < 1.e-15)
return;
- number2 s = (storage(j, j) < 0) ? std::sqrt(sigma) : -std::sqrt(sigma);
+ number2 s;
+ if constexpr (numbers::NumberTraits<number2>::is_complex)
+ s = storage(j, j).real() < 0 ? std::sqrt(sigma) : -std::sqrt(sigma);
+ else
+ s = storage(j, j) < 0 ? std::sqrt(sigma) : -std::sqrt(sigma);
//
number2 beta = std::sqrt(1. / (sigma - s * storage(j, j)));
template <typename number>
template <typename number2>
-double
+number2
Householder<number>::least_squares(Vector<number2> &dst,
const Vector<number2> &src) const
{
const size_type m = storage.m(), n = storage.n();
- GrowingVectorMemory<Vector<number2>> mem;
- typename VectorMemory<Vector<number2>>::Pointer aux(mem);
- aux->reinit(src, true);
- *aux = src;
+ Vector<number2> aux(src);
// m > n, m = src.n, n = dst.n
// Multiply Q_n ... Q_2 Q_1 src
for (size_type j = 0; j < n; ++j)
{
// sum = v_i^T dst
- number2 sum = diagonal[j] * (*aux)(j);
+ number2 sum = diagonal[j] * aux(j);
for (size_type i = j + 1; i < m; ++i)
- sum += static_cast<number2>(storage(i, j)) * (*aux)(i);
+ sum += static_cast<number2>(storage(i, j)) * aux(i);
// dst -= v * sum
- (*aux)(j) -= sum * diagonal[j];
+ aux(j) -= sum * diagonal[j];
for (size_type i = j + 1; i < m; ++i)
- (*aux)(i) -= sum * static_cast<number2>(storage(i, j));
+ aux(i) -= sum * static_cast<number2>(storage(i, j));
}
// Compute norm of residual
number2 sum = 0.;
for (size_type i = n; i < m; ++i)
- sum += (*aux)(i) * (*aux)(i);
+ sum += aux(i) * aux(i);
AssertIsFinite(sum);
// Compute solution
- storage.backward(dst, *aux);
+ storage.backward(dst, aux);
return std::sqrt(sum);
}
const size_type m = storage.m(), n = storage.n();
- GrowingVectorMemory<BlockVector<number2>> mem;
- typename VectorMemory<BlockVector<number2>>::Pointer aux(mem);
- aux->reinit(src, true);
- *aux = src;
+ BlockVector<number2> aux;
+ aux.reinit(src, true);
+ aux = src;
// m > n, m = src.n, n = dst.n
// Multiply Q_n ... Q_2 Q_1 src
for (size_type j = 0; j < n; ++j)
{
// sum = v_i^T dst
- number2 sum = diagonal[j] * (*aux)(j);
+ number2 sum = diagonal[j] * aux(j);
for (size_type i = j + 1; i < m; ++i)
- sum += storage(i, j) * (*aux)(i);
+ sum += storage(i, j) * aux(i);
// dst -= v * sum
- (*aux)(j) -= sum * diagonal[j];
+ aux(j) -= sum * diagonal[j];
for (size_type i = j + 1; i < m; ++i)
- (*aux)(i) -= sum * storage(i, j);
+ aux(i) -= sum * storage(i, j);
}
// Compute norm of residual
number2 sum = 0.;
for (size_type i = n; i < m; ++i)
- sum += (*aux)(i) * (*aux)(i);
+ sum += *aux(i) * aux(i);
AssertIsFinite(sum);
- // backward works for
- // Vectors only, so copy
- // them before
+ // backward works for Vectors only, so copy them before
Vector<number2> v_dst, v_aux;
v_dst = dst;
- v_aux = *aux;
+ v_aux = aux;
// Compute solution
storage.backward(v_dst, v_aux);
- // copy the result back
- // to the BlockVector
+ // copy the result back to the BlockVector
dst = v_dst;
return std::sqrt(sum);