From 4e3dccb27af1015495d313f525e63ee5573cda42 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Tue, 20 May 2025 14:38:52 +0200 Subject: [PATCH] Make Householder class compile with complex numbers --- include/deal.II/lac/householder.h | 55 +++++++++++++++---------------- 1 file changed, 26 insertions(+), 29 deletions(-) diff --git a/include/deal.II/lac/householder.h b/include/deal.II/lac/householder.h index c0a64e3b70..0f5051b36c 100644 --- a/include/deal.II/lac/householder.h +++ b/include/deal.II/lac/householder.h @@ -115,7 +115,7 @@ public: * return. */ template - double + number2 least_squares(Vector &dst, const Vector &src) const; /** @@ -153,7 +153,7 @@ private: /** * Storage that is internally used for the Householder transformation. */ - FullMatrix storage; + FullMatrix storage; }; /** @} */ @@ -190,7 +190,11 @@ Householder::initialize(const FullMatrix &M) 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::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))); @@ -232,7 +236,7 @@ Householder::Householder(const FullMatrix &M) template template -double +number2 Householder::least_squares(Vector &dst, const Vector &src) const { @@ -242,10 +246,7 @@ Householder::least_squares(Vector &dst, const size_type m = storage.m(), n = storage.n(); - GrowingVectorMemory> mem; - typename VectorMemory>::Pointer aux(mem); - aux->reinit(src, true); - *aux = src; + Vector aux(src); // m > n, m = src.n, n = dst.n // Multiply Q_n ... Q_2 Q_1 src @@ -253,22 +254,22 @@ Householder::least_squares(Vector &dst, 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(storage(i, j)) * (*aux)(i); + sum += static_cast(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(storage(i, j)); + aux(i) -= sum * static_cast(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); } @@ -287,10 +288,9 @@ Householder::least_squares(BlockVector &dst, const size_type m = storage.m(), n = storage.n(); - GrowingVectorMemory> mem; - typename VectorMemory>::Pointer aux(mem); - aux->reinit(src, true); - *aux = src; + BlockVector aux; + aux.reinit(src, true); + aux = src; // m > n, m = src.n, n = dst.n // Multiply Q_n ... Q_2 Q_1 src @@ -298,30 +298,27 @@ Householder::least_squares(BlockVector &dst, 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 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); -- 2.39.5