Support complex numbers in serial linear solvers, add BLAS zdot

src/DataStructures/DynamicMatrix.cpp

@@ -3,7 +3,10 @@
 
 #include "DataStructures/DynamicMatrix.hpp"
 
+#include <complex>
+
 #include "Options/ParseOptions.hpp"
+#include "Options/StdComplex.hpp"
 #include "Utilities/GenerateInstantiations.hpp"
 
 namespace DynamicMatrix_detail {
@@ -20,7 +23,7 @@ std::vector<std::vector<Type>> parse_to_vectors(
   template std::vector<std::vector<TYPE(data)>> parse_to_vectors( \
       const Options::Option& options);
 
-GENERATE_INSTANTIATIONS(INSTANTIATE, (double))
+GENERATE_INSTANTIATIONS(INSTANTIATE, (double, std::complex<double>))
 
 #undef INSTANTIATE
 #undef TYPE

src/DataStructures/DynamicVector.cpp

@@ -3,8 +3,11 @@
 
 #include "DataStructures/DynamicVector.hpp"
 
+#include <complex>
+
 #include "Options/Options.hpp"
 #include "Options/ParseOptions.hpp"
+#include "Options/StdComplex.hpp"
 #include "Utilities/GenerateInstantiations.hpp"
 
 namespace DynamicVector_detail {
@@ -20,7 +23,7 @@ std::vector<T> parse_to_vector(const Options::Option& options) {
   template std::vector<TYPE(data)> parse_to_vector( \
       const Options::Option& options);
 
-GENERATE_INSTANTIATIONS(INSTANTIATE, (double))
+GENERATE_INSTANTIATIONS(INSTANTIATE, (double, std::complex<double>))
 
 #undef INSTANTIATE
 #undef TYPE

src/Elliptic/SubdomainPreconditioners/MinusLaplacian.hpp

@@ -55,7 +55,7 @@ template <size_t Dim, typename OptionsGroup,
               ::LinearSolver::Serial::Registrars::Gmres<
                   ::LinearSolver::Schwarz::ElementCenteredSubdomainData<
                       Dim, tmpl::list<Poisson::Tags::Field>>>,
-              ::LinearSolver::Serial::Registrars::ExplicitInverse>>>
+              ::LinearSolver::Serial::Registrars::ExplicitInverse<double>>>>
 struct MinusLaplacian {
   template <typename LinearSolverRegistrars>
   using f = subdomain_preconditioners::MinusLaplacian<Dim, OptionsGroup, Solver,
@@ -97,7 +97,7 @@ template <size_t Dim, typename OptionsGroup,
               ::LinearSolver::Serial::Registrars::Gmres<
                   ::LinearSolver::Schwarz::ElementCenteredSubdomainData<
                       Dim, tmpl::list<Poisson::Tags::Field>>>,
-              ::LinearSolver::Serial::Registrars::ExplicitInverse>>,
+              ::LinearSolver::Serial::Registrars::ExplicitInverse<double>>>,
           typename LinearSolverRegistrars =
               tmpl::list<Registrars::MinusLaplacian<Dim, OptionsGroup, Solver>>>
 class MinusLaplacian

src/Elliptic/SubdomainPreconditioners/RegisterDerived.cpp

@@ -14,11 +14,12 @@
 namespace {
 template <size_t Dim>
 void register_derived_with_charm_impl() {
-  register_derived_classes_with_charm<LinearSolver::Serial::LinearSolver<
-      tmpl::list<::LinearSolver::Serial::Registrars::Gmres<
-                     ::LinearSolver::Schwarz::ElementCenteredSubdomainData<
-                         Dim, tmpl::list<Poisson::Tags::Field>>>,
-                 ::LinearSolver::Serial::Registrars::ExplicitInverse>>>();
+  register_derived_classes_with_charm<
+      LinearSolver::Serial::LinearSolver<tmpl::list<
+          ::LinearSolver::Serial::Registrars::Gmres<
+              ::LinearSolver::Schwarz::ElementCenteredSubdomainData<
+                  Dim, tmpl::list<Poisson::Tags::Field>>>,
+          ::LinearSolver::Serial::Registrars::ExplicitInverse<double>>>>();
 }
 }  // namespace
 

src/NumericalAlgorithms/LinearSolver/BuildMatrix.hpp

@@ -56,7 +56,7 @@ void build_matrix(const gsl::not_null<MatrixType*> matrix,
   // Re-using the iterators for all operator invocations
   auto result_iterator_begin = result_buffer->begin();
   auto result_iterator_end = result_buffer->end();
-  for (double& unit_vector_data : *operand_buffer) {
+  for (auto& unit_vector_data : *operand_buffer) {
     // Set a 1 at the unit vector location i
     unit_vector_data = 1.;
     // Invoke the operator on the unit vector

src/NumericalAlgorithms/LinearSolver/ExplicitInverse.hpp

@@ -30,13 +30,17 @@
 namespace LinearSolver::Serial {
 
 /// \cond
-template <typename LinearSolverRegistrars>
+template <typename ValueType, typename LinearSolverRegistrars>
 struct ExplicitInverse;
 /// \endcond
 
 namespace Registrars {
 /// Registers the `LinearSolver::Serial::ExplicitInverse` linear solver
-using ExplicitInverse = Registration::Registrar<Serial::ExplicitInverse>;
+template <typename ValueType>
+struct ExplicitInverse {
+  template <typename LinearSolverRegistrars>
+  using f = Serial::ExplicitInverse<ValueType, LinearSolverRegistrars>;
+};
 }  // namespace Registrars
 
 /*!
@@ -73,8 +77,9 @@ using ExplicitInverse = Registration::Registrar<Serial::ExplicitInverse>;
  *   subdomain problems only approximately, but possibly still sufficiently to
  *   provide effective preconditioning.
  */
-template <typename LinearSolverRegistrars =
-              tmpl::list<Registrars::ExplicitInverse>>
+template <typename ValueType,
+          typename LinearSolverRegistrars =
+              tmpl::list<Registrars::ExplicitInverse<ValueType>>>
 class ExplicitInverse : public LinearSolver<LinearSolverRegistrars> {
  private:
   using Base = LinearSolver<LinearSolverRegistrars>;
@@ -143,7 +148,7 @@ class ExplicitInverse : public LinearSolver<LinearSolverRegistrars> {
 
   /// The matrix representation of the solver. This matrix approximates the
   /// inverse of the subdomain operator.
-  const blaze::DynamicMatrix<double, blaze::columnMajor>&
+  const blaze::DynamicMatrix<ValueType, blaze::columnMajor>&
   matrix_representation() const {
     return inverse_;
   }
@@ -171,19 +176,20 @@ class ExplicitInverse : public LinearSolver<LinearSolverRegistrars> {
   // We currently store the matrix representation in a dense matrix because
   // Blaze doesn't support the inversion of sparse matrices (yet).
   // NOLINTNEXTLINE(spectre-mutable)
-  mutable blaze::DynamicMatrix<double, blaze::columnMajor> inverse_{};
+  mutable blaze::DynamicMatrix<ValueType, blaze::columnMajor> inverse_{};
 
   // Buffers to avoid re-allocating memory for applying the operator
   // NOLINTNEXTLINE(spectre-mutable)
-  mutable blaze::DynamicVector<double> source_workspace_{};
+  mutable blaze::DynamicVector<ValueType> source_workspace_{};
   // NOLINTNEXTLINE(spectre-mutable)
-  mutable blaze::DynamicVector<double> solution_workspace_{};
+  mutable blaze::DynamicVector<ValueType> solution_workspace_{};
 };
 
-template <typename LinearSolverRegistrars>
+template <typename ValueType, typename LinearSolverRegistrars>
 template <typename LinearOperator, typename VarsType, typename SourceType,
           typename... OperatorArgs>
-Convergence::HasConverged ExplicitInverse<LinearSolverRegistrars>::solve(
+Convergence::HasConverged
+ExplicitInverse<ValueType, LinearSolverRegistrars>::solve(
     const gsl::not_null<VarsType*> solution,
     const LinearOperator& linear_operator, const SourceType& source,
     const std::tuple<OperatorArgs...>& operator_args) const {
@@ -245,9 +251,11 @@ Convergence::HasConverged ExplicitInverse<LinearSolverRegistrars>::solve(
 }
 
 /// \cond
-template <typename LinearSolverRegistrars>
-// NOLINTNEXTLINE
-PUP::able::PUP_ID ExplicitInverse<LinearSolverRegistrars>::my_PUP_ID = 0;
+// NOLINTBEGIN
+template <typename ValueType, typename LinearSolverRegistrars>
+PUP::able::PUP_ID
+    ExplicitInverse<ValueType, LinearSolverRegistrars>::my_PUP_ID = 0;
+// NOLINTEND
 /// \endcond
 
 }  // namespace LinearSolver::Serial

src/NumericalAlgorithms/LinearSolver/Gmres.cpp

@@ -9,57 +9,65 @@
 #include <blaze/math/Submatrix.h>
 #include <blaze/math/Subvector.h>
 #include <blaze/math/lapack/trsv.h>
+#include <complex>
 
 #include "Utilities/ConstantExpressions.hpp"
+#include "Utilities/ErrorHandling/Assert.hpp"
+#include "Utilities/GenerateInstantiations.hpp"
 #include "Utilities/Gsl.hpp"
 
 namespace LinearSolver::gmres::detail {
 
 namespace {
 
-void update_givens_rotation(const gsl::not_null<double*> givens_cosine,
+template <typename ValueType>
+void update_givens_rotation(const gsl::not_null<ValueType*> givens_cosine,
                             const gsl::not_null<double*> givens_sine,
-                            const double rho, const double sigma) {
+                            const ValueType rho, const double sigma) {
   if (UNLIKELY(rho == 0.)) {
     *givens_cosine = 0.;
     *givens_sine = 1.;
   } else {
-    const double tmp = sqrt(square(rho) + square(sigma));
-    *givens_cosine = abs(rho) / tmp;
-    *givens_sine = *givens_cosine * sigma / rho;
+    const double tmp = sqrt(square(abs(rho)) + square(sigma));
+    *givens_cosine = rho / tmp;
+    *givens_sine = sigma / tmp;
   }
 }
 
-template <typename Arg>
+template <typename ValueType, typename Arg>
 void apply_and_update_givens_rotation(
     Arg argument,
     const gsl::not_null<blaze::DynamicVector<double>*> givens_sine_history,
-    const gsl::not_null<blaze::DynamicVector<double>*> givens_cosine_history,
+    const gsl::not_null<blaze::DynamicVector<ValueType>*> givens_cosine_history,
     const size_t iteration) {
   const size_t k = iteration + 1;
   for (size_t i = 0; i < k - 1; ++i) {
-    const double tmp = (*givens_cosine_history)[i] * argument[i] +
-                       (*givens_sine_history)[i] * argument[i + 1];
+    const ValueType tmp = (*givens_cosine_history)[i] * argument[i] +
+                          (*givens_sine_history)[i] * argument[i + 1];
     argument[i + 1] = (*givens_cosine_history)[i] * argument[i + 1] -
                       (*givens_sine_history)[i] * argument[i];
     argument[i] = tmp;
   }
+  // Note: argument[k] is real, since it is a normalization
+  ASSERT(equal_within_roundoff(std::imag(argument[k]), 0.),
+         "Normalization is not real: " << argument[k]);
   update_givens_rotation(make_not_null(&(*givens_cosine_history)[k - 1]),
                          make_not_null(&(*givens_sine_history)[k - 1]),
-                         argument[k - 1], argument[k]);
+                         argument[k - 1], std::real(argument[k]));
   argument[k - 1] = (*givens_cosine_history)[k - 1] * argument[k - 1] +
                     (*givens_sine_history)[k - 1] * argument[k];
   argument[k] = 0.;
 }
 
 }  // namespace
 
+template <typename ValueType>
 void solve_minimal_residual(
-    const gsl::not_null<blaze::DynamicMatrix<double>*>
+    const gsl::not_null<blaze::DynamicMatrix<ValueType>*>
         orthogonalization_history,
-    const gsl::not_null<blaze::DynamicVector<double>*> residual_history,
+    const gsl::not_null<blaze::DynamicVector<ValueType>*> residual_history,
     const gsl::not_null<blaze::DynamicVector<double>*> givens_sine_history,
-    const gsl::not_null<blaze::DynamicVector<double>*> givens_cosine_history,
+    const gsl::not_null<blaze::DynamicVector<ValueType>*> givens_cosine_history,
     const size_t iteration) {
   residual_history->resize(iteration + 2);
   givens_sine_history->resize(iteration + 1);
@@ -81,15 +89,33 @@ void solve_minimal_residual(
       (*givens_cosine_history)[iteration] * (*residual_history)[iteration];
 }
 
-blaze::DynamicVector<double> minimal_residual_vector(
-    const blaze::DynamicMatrix<double>& orthogonalization_history,
-    const blaze::DynamicVector<double>& residual_history) {
+template <typename ValueType>
+blaze::DynamicVector<ValueType> minimal_residual_vector(
+    const blaze::DynamicMatrix<ValueType>& orthogonalization_history,
+    const blaze::DynamicVector<ValueType>& residual_history) {
   const size_t length = orthogonalization_history.columns();
-  blaze::DynamicVector<double> minres =
+  blaze::DynamicVector<ValueType> minres =
       blaze::subvector(residual_history, 0, length);
   blaze::trsv(blaze::submatrix(orthogonalization_history, 0, 0, length, length),
               minres, 'U', 'N', 'N');
   return minres;
 }
 
+#define DTYPE(data) BOOST_PP_TUPLE_ELEM(0, data)
+
+#define INSTANTIATE(r, data)                                          \
+  template void solve_minimal_residual(                               \
+      gsl::not_null<blaze::DynamicMatrix<DTYPE(data)>*>,              \
+      gsl::not_null<blaze::DynamicVector<DTYPE(data)>*>,              \
+      gsl::not_null<blaze::DynamicVector<double>*>,                   \
+      gsl::not_null<blaze::DynamicVector<DTYPE(data)>*>, size_t);     \
+  template blaze::DynamicVector<DTYPE(data)> minimal_residual_vector( \
+      const blaze::DynamicMatrix<DTYPE(data)>&,                       \
+      const blaze::DynamicVector<DTYPE(data)>&);
+
+GENERATE_INSTANTIATIONS(INSTANTIATE, (double, std::complex<double>))
+
+#undef DTYPE
+#undef INSTANTIATE
+
 }  // namespace LinearSolver::gmres::detail