diff --git a/src/ParallelAlgorithms/LinearSolver/Gmres/ElementActions.hpp b/src/ParallelAlgorithms/LinearSolver/Gmres/ElementActions.hpp
index 56a19f32bdf96..ef0541c8e2b07 100644
--- a/src/ParallelAlgorithms/LinearSolver/Gmres/ElementActions.hpp
+++ b/src/ParallelAlgorithms/LinearSolver/Gmres/ElementActions.hpp
@@ -33,6 +33,7 @@
 #include "Utilities/PrettyType.hpp"
 #include "Utilities/Requires.hpp"
 #include "Utilities/TMPL.hpp"
+#include "Utilities/TypeTraits/GetFundamentalType.hpp"
 
 /// \cond
 namespace tuples {
@@ -300,6 +301,8 @@ struct PerformStep {
       db::add_tag_prefix<LinearSolver::Tags::Operand, fields_tag>;
   using preconditioned_operand_tag =
       db::add_tag_prefix<LinearSolver::Tags::Preconditioned, operand_tag>;
+  using ValueType =
+      tt::get_complex_or_fundamental_type_t<typename fields_tag::type>;
 
  public:
   using const_global_cache_tags =
@@ -374,7 +377,7 @@ struct PerformStep {
         Parallel::ReductionData<
             Parallel::ReductionDatum<size_t, funcl::AssertEqual<>>,
             Parallel::ReductionDatum<size_t, funcl::AssertEqual<>>,
-            Parallel::ReductionDatum<double, funcl::Plus<>>>{
+            Parallel::ReductionDatum<ValueType, funcl::Plus<>>>{
             get<Convergence::Tags::IterationId<OptionsGroup>>(box),
             get<orthogonalization_iteration_id_tag>(box),
             inner_product(get<basis_history_tag>(box)[0],
@@ -400,9 +403,12 @@ struct OrthogonalizeOperand {
           Convergence::Tags::IterationId<OptionsGroup>>;
   using basis_history_tag =
       LinearSolver::Tags::KrylovSubspaceBasis<operand_tag>;
+  using ValueType =
+      tt::get_complex_or_fundamental_type_t<typename fields_tag::type>;
 
  public:
-  using inbox_tags = tmpl::list<Tags::Orthogonalization<OptionsGroup>>;
+  using inbox_tags =
+      tmpl::list<Tags::Orthogonalization<OptionsGroup, ValueType>>;
 
   template <typename DbTagsList, typename... InboxTags, typename Metavariables,
             typename ArrayIndex, typename ActionList,
@@ -414,12 +420,13 @@ struct OrthogonalizeOperand {
       const ParallelComponent* const /*meta*/) {
     const size_t iteration_id =
         db::get<Convergence::Tags::IterationId<OptionsGroup>>(box);
-    auto& inbox = get<Tags::Orthogonalization<OptionsGroup>>(inboxes);
+    auto& inbox =
+        get<Tags::Orthogonalization<OptionsGroup, ValueType>>(inboxes);
     if (inbox.find(iteration_id) == inbox.end()) {
       return {Parallel::AlgorithmExecution::Retry, std::nullopt};
     }
 
-    const double orthogonalization =
+    const ValueType orthogonalization =
         std::move(inbox.extract(iteration_id).mapped());
 
     db::mutate<operand_tag, orthogonalization_iteration_id_tag>(
@@ -437,7 +444,7 @@ struct OrthogonalizeOperand {
         get<orthogonalization_iteration_id_tag>(box);
     const bool orthogonalization_complete =
         next_orthogonalization_iteration_id == iteration_id;
-    const double local_orthogonalization =
+    const ValueType local_orthogonalization =
         inner_product(orthogonalization_complete
                           ? get<operand_tag>(box)
                           : gsl::at(get<basis_history_tag>(box),
@@ -451,7 +458,7 @@ struct OrthogonalizeOperand {
         Parallel::ReductionData<
             Parallel::ReductionDatum<size_t, funcl::AssertEqual<>>,
             Parallel::ReductionDatum<size_t, funcl::AssertEqual<>>,
-            Parallel::ReductionDatum<double, funcl::Plus<>>>{
+            Parallel::ReductionDatum<ValueType, funcl::Plus<>>>{
             iteration_id, next_orthogonalization_iteration_id,
             local_orthogonalization},
         Parallel::get_parallel_component<ParallelComponent>(cache)[array_index],
@@ -483,11 +490,14 @@ struct NormalizeOperandAndUpdateField {
   using preconditioned_basis_history_tag =
       LinearSolver::Tags::KrylovSubspaceBasis<std::conditional_t<
           Preconditioned, preconditioned_operand_tag, operand_tag>>;
+  using ValueType =
+      tt::get_complex_or_fundamental_type_t<typename fields_tag::type>;
 
  public:
   using const_global_cache_tags =
       tmpl::list<logging::Tags::Verbosity<OptionsGroup>>;
-  using inbox_tags = tmpl::list<Tags::FinalOrthogonalization<OptionsGroup>>;
+  using inbox_tags =
+      tmpl::list<Tags::FinalOrthogonalization<OptionsGroup, ValueType>>;
 
   template <typename DbTagsList, typename... InboxTags, typename Metavariables,
             typename ArrayIndex, typename ActionList,
@@ -499,7 +509,8 @@ struct NormalizeOperandAndUpdateField {
       const ParallelComponent* const /*meta*/) {
     const size_t iteration_id =
         db::get<Convergence::Tags::IterationId<OptionsGroup>>(box);
-    auto& inbox = get<Tags::FinalOrthogonalization<OptionsGroup>>(inboxes);
+    auto& inbox =
+        get<Tags::FinalOrthogonalization<OptionsGroup, ValueType>>(inboxes);
     if (inbox.find(iteration_id) == inbox.end()) {
       return {Parallel::AlgorithmExecution::Retry, std::nullopt};
     }
diff --git a/src/ParallelAlgorithms/LinearSolver/Gmres/ResidualMonitorActions.hpp b/src/ParallelAlgorithms/LinearSolver/Gmres/ResidualMonitorActions.hpp
index 52669e4103157..08d786e6b4fde 100644
--- a/src/ParallelAlgorithms/LinearSolver/Gmres/ResidualMonitorActions.hpp
+++ b/src/ParallelAlgorithms/LinearSolver/Gmres/ResidualMonitorActions.hpp
@@ -14,6 +14,7 @@
 #include "IO/Logging/Tags.hpp"
 #include "IO/Logging/Verbosity.hpp"
 #include "NumericalAlgorithms/Convergence/Tags.hpp"
+#include "NumericalAlgorithms/LinearSolver/InnerProduct.hpp"
 #include "Parallel/GlobalCache.hpp"
 #include "Parallel/Invoke.hpp"
 #include "Parallel/Printf/Printf.hpp"
@@ -106,6 +107,8 @@ struct StoreOrthogonalization {
       ::Tags::Previous<residual_magnitude_tag>;
   using orthogonalization_history_tag =
       LinearSolver::Tags::OrthogonalizationHistory<fields_tag>;
+  using ValueType =
+      tt::get_complex_or_fundamental_type_t<typename fields_tag::type>;
 
  public:
   template <typename ParallelComponent, typename DbTagsList,
@@ -116,7 +119,7 @@ struct StoreOrthogonalization {
                     const ArrayIndex& /*array_index*/,
                     const size_t iteration_id,
                     const size_t orthogonalization_iteration_id,
-                    const double orthogonalization) {
+                    const ValueType orthogonalization) {
     if (UNLIKELY(orthogonalization_iteration_id == 0)) {
       // Append a row and a column to the orthogonalization history. Zero the
       // entries that won't be set during the orthogonalization procedure below.
@@ -143,19 +146,21 @@ struct StoreOrthogonalization {
           },
           make_not_null(&box));
 
-      Parallel::receive_data<Tags::Orthogonalization<OptionsGroup>>(
+      Parallel::receive_data<Tags::Orthogonalization<OptionsGroup, ValueType>>(
           Parallel::get_parallel_component<BroadcastTarget>(cache),
           iteration_id, orthogonalization);
       return;
     }
 
     // At this point, the orthogonalization procedure is complete.
+    ASSERT(equal_within_roundoff(imag(orthogonalization), 0.0),
+           "Normalization is not real: " << orthogonalization);
+    const double normalization = sqrt(real(orthogonalization));
     db::mutate<orthogonalization_history_tag>(
-        [orthogonalization, iteration_id,
+        [normalization, iteration_id,
          orthogonalization_iteration_id](const auto orthogonalization_history) {
           (*orthogonalization_history)(orthogonalization_iteration_id,
-                                       iteration_id - 1) =
-              sqrt(orthogonalization);
+                                       iteration_id - 1) = normalization;
         },
         make_not_null(&box));
 
@@ -164,18 +169,19 @@ struct StoreOrthogonalization {
     const auto& orthogonalization_history =
         get<orthogonalization_history_tag>(box);
     const auto num_rows = orthogonalization_iteration_id + 1;
-    blaze::DynamicMatrix<double> qr_Q;
-    blaze::DynamicMatrix<double> qr_R;
+    blaze::DynamicMatrix<ValueType> qr_Q;
+    blaze::DynamicMatrix<ValueType> qr_R;
     blaze::qr(orthogonalization_history, qr_Q, qr_R);
     // Compute the residual vector from the QR decomposition
     blaze::DynamicVector<double> beta(num_rows, 0.);
     const double initial_residual_magnitude =
         get<initial_residual_magnitude_tag>(box);
     beta[0] = initial_residual_magnitude;
-    blaze::DynamicVector<double> minres =
-        blaze::inv(qr_R) * blaze::trans(qr_Q) * beta;
-    const double residual_magnitude =
-        blaze::length(beta - orthogonalization_history * minres);
+    blaze::DynamicVector<ValueType> minres =
+        blaze::inv(qr_R) * blaze::ctrans(qr_Q) * beta;
+    blaze::DynamicVector<ValueType> res =
+        beta - orthogonalization_history * minres;
+    const double residual_magnitude = sqrt(magnitude_square(res));
 
     // At this point, the iteration is complete. We proceed with observing,
     // logging and checking convergence before broadcasting back to the
@@ -233,9 +239,10 @@ struct StoreOrthogonalization {
       }
     }
 
-    Parallel::receive_data<Tags::FinalOrthogonalization<OptionsGroup>>(
+    Parallel::receive_data<
+        Tags::FinalOrthogonalization<OptionsGroup, ValueType>>(
         Parallel::get_parallel_component<BroadcastTarget>(cache), iteration_id,
-        std::make_tuple(sqrt(orthogonalization), std::move(minres),
+        std::make_tuple(normalization, std::move(minres),
                         // NOLINTNEXTLINE(performance-move-const-arg)
                         std::move(has_converged)));
   }
diff --git a/src/ParallelAlgorithms/LinearSolver/Gmres/Tags/InboxTags.hpp b/src/ParallelAlgorithms/LinearSolver/Gmres/Tags/InboxTags.hpp
index d39d68e579a62..1aa7516038843 100644
--- a/src/ParallelAlgorithms/LinearSolver/Gmres/Tags/InboxTags.hpp
+++ b/src/ParallelAlgorithms/LinearSolver/Gmres/Tags/InboxTags.hpp
@@ -22,19 +22,20 @@ struct InitialOrthogonalization
       std::map<temporal_id, std::tuple<double, Convergence::HasConverged>>;
 };
 
-template <typename OptionsGroup>
-struct Orthogonalization
-    : Parallel::InboxInserters::Value<Orthogonalization<OptionsGroup>> {
+template <typename OptionsGroup, typename ValueType>
+struct Orthogonalization : Parallel::InboxInserters::Value<
+                               Orthogonalization<OptionsGroup, ValueType>> {
   using temporal_id = size_t;
-  using type = std::map<temporal_id, double>;
+  using type = std::map<temporal_id, ValueType>;
 };
 
-template <typename OptionsGroup>
+template <typename OptionsGroup, typename ValueType>
 struct FinalOrthogonalization
-    : Parallel::InboxInserters::Value<FinalOrthogonalization<OptionsGroup>> {
+    : Parallel::InboxInserters::Value<
+          FinalOrthogonalization<OptionsGroup, ValueType>> {
   using temporal_id = size_t;
   using type =
-      std::map<temporal_id, std::tuple<double, blaze::DynamicVector<double>,
+      std::map<temporal_id, std::tuple<double, blaze::DynamicVector<ValueType>,
                                        Convergence::HasConverged>>;
 };
 
diff --git a/src/ParallelAlgorithms/LinearSolver/Tags.hpp b/src/ParallelAlgorithms/LinearSolver/Tags.hpp
index 6a7da4a09ff8e..4a65f0e1b8d90 100644
--- a/src/ParallelAlgorithms/LinearSolver/Tags.hpp
+++ b/src/ParallelAlgorithms/LinearSolver/Tags.hpp
@@ -17,6 +17,7 @@
 #include "DataStructures/DataBox/TagName.hpp"
 #include "DataStructures/DynamicMatrix.hpp"
 #include "Utilities/Gsl.hpp"
+#include "Utilities/TypeTraits/GetFundamentalType.hpp"
 
 /*!
  * \ingroup LinearSolverGroup
@@ -137,6 +138,7 @@ struct Orthogonalization : db::PrefixTag, db::SimpleTag {
  */
 template <typename Tag>
 struct OrthogonalizationHistory : db::PrefixTag, db::SimpleTag {
+  using ValueType = tt::get_complex_or_fundamental_type_t<typename Tag::type>;
   static std::string name() {
     // Add "Linear" prefix to abbreviate the namespace for uniqueness
     return "LinearOrthogonalizationHistory(" + db::tag_name<Tag>() + ")";
diff --git a/tests/Unit/Helpers/ParallelAlgorithms/LinearSolver/LinearSolverAlgorithmTestHelpers.hpp b/tests/Unit/Helpers/ParallelAlgorithms/LinearSolver/LinearSolverAlgorithmTestHelpers.hpp
index 7a911dc48c2d0..39e4ffe41c0b2 100644
--- a/tests/Unit/Helpers/ParallelAlgorithms/LinearSolver/LinearSolverAlgorithmTestHelpers.hpp
+++ b/tests/Unit/Helpers/ParallelAlgorithms/LinearSolver/LinearSolverAlgorithmTestHelpers.hpp
@@ -46,25 +46,30 @@
 #include "Utilities/Requires.hpp"
 #include "Utilities/TMPL.hpp"
 #include "Utilities/TaggedTuple.hpp"
+#include "Utilities/TypeTraits/GetFundamentalType.hpp"
 
 namespace LinearSolverAlgorithmTestHelpers {
 
 namespace OptionTags {
+template <typename ValueType>
 struct LinearOperator {
   static constexpr Options::String help = "The linear operator A to invert.";
-  using type = blaze::DynamicMatrix<double>;
+  using type = blaze::DynamicMatrix<ValueType>;
 };
+template <typename ValueType>
 struct Source {
   static constexpr Options::String help = "The source b in the equation Ax=b.";
-  using type = blaze::DynamicVector<double>;
+  using type = blaze::DynamicVector<ValueType>;
 };
+template <typename ValueType>
 struct InitialGuess {
   static constexpr Options::String help = "The initial guess for the vector x.";
-  using type = blaze::DynamicVector<double>;
+  using type = blaze::DynamicVector<ValueType>;
 };
+template <typename ValueType>
 struct ExpectedResult {
   static constexpr Options::String help = "The solution x in the equation Ax=b";
-  using type = blaze::DynamicVector<double>;
+  using type = blaze::DynamicVector<ValueType>;
 };
 struct ExpectedConvergenceReason {
   static std::string name() { return "ConvergenceReason"; }
@@ -73,9 +78,10 @@ struct ExpectedConvergenceReason {
 };
 }  // namespace OptionTags
 
+template <typename ValueType>
 struct LinearOperator : db::SimpleTag {
-  using type = blaze::DynamicMatrix<double>;
-  using option_tags = tmpl::list<OptionTags::LinearOperator>;
+  using type = blaze::DynamicMatrix<ValueType>;
+  using option_tags = tmpl::list<OptionTags::LinearOperator<ValueType>>;
 
   static constexpr bool pass_metavariables = false;
   static type create_from_options(const type& linear_operator) {
@@ -83,17 +89,19 @@ struct LinearOperator : db::SimpleTag {
   }
 };
 
+template <typename ValueType>
 struct Source : db::SimpleTag {
-  using type = blaze::DynamicVector<double>;
-  using option_tags = tmpl::list<OptionTags::Source>;
+  using type = blaze::DynamicVector<ValueType>;
+  using option_tags = tmpl::list<OptionTags::Source<ValueType>>;
 
   static constexpr bool pass_metavariables = false;
   static type create_from_options(const type& source) { return source; }
 };
 
+template <typename ValueType>
 struct InitialGuess : db::SimpleTag {
-  using type = blaze::DynamicVector<double>;
-  using option_tags = tmpl::list<OptionTags::InitialGuess>;
+  using type = blaze::DynamicVector<ValueType>;
+  using option_tags = tmpl::list<OptionTags::InitialGuess<ValueType>>;
 
   static constexpr bool pass_metavariables = false;
   static type create_from_options(const type& initial_guess) {
@@ -101,9 +109,10 @@ struct InitialGuess : db::SimpleTag {
   }
 };
 
+template <typename ValueType>
 struct ExpectedResult : db::SimpleTag {
-  using type = blaze::DynamicVector<double>;
-  using option_tags = tmpl::list<OptionTags::ExpectedResult>;
+  using type = blaze::DynamicVector<ValueType>;
+  using option_tags = tmpl::list<OptionTags::ExpectedResult<ValueType>>;
 
   static constexpr bool pass_metavariables = false;
   static type create_from_options(const type& expected_result) {
@@ -122,11 +131,13 @@ struct ExpectedConvergenceReason : db::SimpleTag {
 };
 
 // The vector `x` we want to solve for
+template <typename ValueType>
 struct VectorTag : db::SimpleTag {
-  using type = blaze::DynamicVector<double>;
+  using type = blaze::DynamicVector<ValueType>;
 };
 
-using fields_tag = VectorTag;
+template <typename ValueType>
+using fields_tag = VectorTag<ValueType>;
 
 template <typename OperandTag>
 struct ComputeOperatorAction {
@@ -142,23 +153,26 @@ struct ComputeOperatorAction {
       const ActionList /*meta*/,
       // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
       const ParallelComponent* const /*meta*/) {
+    using ValueType =
+        tt::get_complex_or_fundamental_type_t<typename OperandTag::type>;
     db::mutate<LinearSolver::Tags::OperatorAppliedTo<OperandTag>>(
-        [](const gsl::not_null<blaze::DynamicVector<double>*>
+        [](const gsl::not_null<blaze::DynamicVector<ValueType>*>
                operator_applied_to_operand,
-           const blaze::DynamicMatrix<double>& linear_operator,
-           const blaze::DynamicVector<double>& operand) {
+           const blaze::DynamicMatrix<ValueType>& linear_operator,
+           const blaze::DynamicVector<ValueType>& operand) {
           *operator_applied_to_operand = linear_operator * operand;
         },
-        make_not_null(&box), get<LinearOperator>(box), get<OperandTag>(box));
+        make_not_null(&box), get<LinearOperator<ValueType>>(box),
+        get<OperandTag>(box));
     return {Parallel::AlgorithmExecution::Continue, std::nullopt};
   }
 };
 
 // Checks for the correct solution after the algorithm has terminated.
-template <typename OptionsGroup>
+template <typename ValueType, typename OptionsGroup>
 struct TestResult {
   using const_global_cache_tags =
-      tmpl::list<ExpectedResult, ExpectedConvergenceReason>;
+      tmpl::list<ExpectedResult<ValueType>, ExpectedConvergenceReason>;
 
   template <typename DbTagsList, typename... InboxTags, typename Metavariables,
             typename ActionList, typename ParallelComponent>
@@ -177,17 +191,26 @@ struct TestResult {
     SPECTRE_PARALLEL_REQUIRE(has_converged);
     SPECTRE_PARALLEL_REQUIRE(has_converged.reason() ==
                              get<ExpectedConvergenceReason>(box));
-    const auto& result = get<VectorTag>(box);
-    const auto& expected_result = get<ExpectedResult>(box);
+    const auto& result = get<VectorTag<ValueType>>(box);
+    const auto& expected_result = get<ExpectedResult<ValueType>>(box);
     for (size_t i = 0; i < expected_result.size(); i++) {
-      SPECTRE_PARALLEL_REQUIRE(result[i] == approx(expected_result[i]));
+      if constexpr (std::is_same_v<ValueType, std::complex<double>>) {
+        SPECTRE_PARALLEL_REQUIRE(real(result[i]) ==
+                                 approx(real(expected_result[i])));
+        SPECTRE_PARALLEL_REQUIRE(imag(result[i]) ==
+                                 approx(imag(expected_result[i])));
+      } else {
+        SPECTRE_PARALLEL_REQUIRE(result[i] == approx(expected_result[i]));
+      }
     }
     return {Parallel::AlgorithmExecution::Pause, std::nullopt};
   }
 };
 
+template <typename ValueType>
 struct InitializeElement {
-  using simple_tags = tmpl::list<VectorTag, ::Tags::FixedSource<VectorTag>>;
+  using simple_tags = tmpl::list<VectorTag<ValueType>,
+                                 ::Tags::FixedSource<VectorTag<ValueType>>>;
   template <typename DbTagsList, typename... InboxTags, typename Metavariables,
             typename ActionList, typename ParallelComponent>
   static Parallel::iterable_action_return_t apply(
@@ -197,7 +220,8 @@ struct InitializeElement {
       const int /*array_index*/, const ActionList /*meta*/,
       const ParallelComponent* const /*meta*/) {
     Initialization::mutate_assign<simple_tags>(
-        make_not_null(&box), get<InitialGuess>(box), get<Source>(box));
+        make_not_null(&box), get<InitialGuess<ValueType>>(box),
+        get<Source<ValueType>>(box));
     return {Parallel::AlgorithmExecution::Continue, std::nullopt};
   }
 };
@@ -262,6 +286,8 @@ struct ElementArray {
   using metavariables = Metavariables;
   using linear_solver = typename Metavariables::linear_solver;
   using preconditioner = typename Metavariables::preconditioner;
+  using ValueType = tt::get_complex_or_fundamental_type_t<
+      typename linear_solver::fields_tag::type>;
 
   // In each step of the algorithm we must provide A(p). The linear solver then
   // takes care of updating x and p, as well as the internal variables r, its
@@ -270,9 +296,9 @@ struct ElementArray {
   using phase_dependent_action_list = tmpl::list<
       Parallel::PhaseActions<
           Parallel::Phase::Initialization,
-          tmpl::list<InitializeElement,
+          tmpl::list<InitializeElement<ValueType>,
                      typename linear_solver::initialize_element,
-                     ComputeOperatorAction<fields_tag>,
+                     ComputeOperatorAction<fields_tag<ValueType>>,
                      detail::init_preconditioner<preconditioner>,
                      Parallel::Actions::TerminatePhase>>,
       Parallel::PhaseActions<
@@ -289,12 +315,14 @@ struct ElementArray {
               Parallel::Actions::TerminatePhase>>,
       Parallel::PhaseActions<
           Parallel::Phase::Testing,
-          tmpl::list<TestResult<typename linear_solver::options_group>>>>;
+          tmpl::list<
+              TestResult<ValueType, typename linear_solver::options_group>>>>;
   /// [action_list]
   using simple_tags_from_options = Parallel::get_simple_tags_from_options<
       Parallel::get_initialization_actions_list<phase_dependent_action_list>>;
   using const_global_cache_tags =
-      tmpl::list<LinearOperator, Source, InitialGuess, ExpectedResult>;
+      tmpl::list<LinearOperator<ValueType>, Source<ValueType>,
+                 InitialGuess<ValueType>, ExpectedResult<ValueType>>;
   using array_allocation_tags = tmpl::list<>;
 
   static void allocate_array(
diff --git a/tests/Unit/ParallelAlgorithms/LinearSolver/ConjugateGradient/Test_ConjugateGradientAlgorithm.cpp b/tests/Unit/ParallelAlgorithms/LinearSolver/ConjugateGradient/Test_ConjugateGradientAlgorithm.cpp
index 096968a2e6316..959a3e7f0db5f 100644
--- a/tests/Unit/ParallelAlgorithms/LinearSolver/ConjugateGradient/Test_ConjugateGradientAlgorithm.cpp
+++ b/tests/Unit/ParallelAlgorithms/LinearSolver/ConjugateGradient/Test_ConjugateGradientAlgorithm.cpp
@@ -25,9 +25,8 @@ struct Metavariables {
   static constexpr const char* const help{
       "Test the conjugate gradient linear solver algorithm"};
 
-  using linear_solver =
-      LinearSolver::cg::ConjugateGradient<Metavariables, helpers::fields_tag,
-                                          SerialCg>;
+  using linear_solver = LinearSolver::cg::ConjugateGradient<
+      Metavariables, helpers::fields_tag<double>, SerialCg>;
   using preconditioner = void;
 
   using component_list = helpers::component_list<Metavariables>;
diff --git a/tests/Unit/ParallelAlgorithms/LinearSolver/Gmres/CMakeLists.txt b/tests/Unit/ParallelAlgorithms/LinearSolver/Gmres/CMakeLists.txt
index 3b96d130e9932..52404dc2f2a8f 100644
--- a/tests/Unit/ParallelAlgorithms/LinearSolver/Gmres/CMakeLists.txt
+++ b/tests/Unit/ParallelAlgorithms/LinearSolver/Gmres/CMakeLists.txt
@@ -28,6 +28,13 @@ target_link_libraries(
   "Test_GmresAlgorithm"
   PRIVATE
   "${INTEGRATION_TEST_LINK_LIBRARIES}")
+add_standalone_test(
+  "Integration.LinearSolver.ComplexGmresAlgorithm"
+  INPUT_FILE "Test_ComplexGmresAlgorithm.yaml")
+target_link_libraries(
+  "Test_ComplexGmresAlgorithm"
+  PRIVATE
+  "${INTEGRATION_TEST_LINK_LIBRARIES}")
 add_standalone_test(
   "Integration.LinearSolver.GmresPreconditionedAlgorithm"
   INPUT_FILE "Test_GmresPreconditionedAlgorithm.yaml")
diff --git a/tests/Unit/ParallelAlgorithms/LinearSolver/Gmres/Test_ComplexGmresAlgorithm.cpp b/tests/Unit/ParallelAlgorithms/LinearSolver/Gmres/Test_ComplexGmresAlgorithm.cpp
new file mode 100644
index 0000000000000..d9ae50fa7096d
--- /dev/null
+++ b/tests/Unit/ParallelAlgorithms/LinearSolver/Gmres/Test_ComplexGmresAlgorithm.cpp
@@ -0,0 +1,49 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include <vector>
+
+#include "Helpers/ParallelAlgorithms/LinearSolver/LinearSolverAlgorithmTestHelpers.hpp"
+#include "Parallel/CharmMain.tpp"
+#include "ParallelAlgorithms/LinearSolver/Gmres/Gmres.hpp"
+#include "Utilities/TMPL.hpp"
+
+namespace PUP {
+class er;
+}  // namespace PUP
+
+namespace helpers = LinearSolverAlgorithmTestHelpers;
+
+namespace {
+
+struct SerialGmres {
+  static constexpr Options::String help =
+      "Options for the iterative linear solver";
+};
+
+struct Metavariables {
+  static constexpr const char* const help{
+      "Test the GMRES linear solver algorithm"};
+
+  using linear_solver =
+      LinearSolver::gmres::Gmres<Metavariables,
+                                 helpers::fields_tag<std::complex<double>>,
+                                 SerialGmres, false>;
+  using preconditioner = void;
+
+  using component_list = helpers::component_list<Metavariables>;
+  using observed_reduction_data_tags =
+      helpers::observed_reduction_data_tags<Metavariables>;
+  static constexpr bool ignore_unrecognized_command_line_options = false;
+  static constexpr auto default_phase_order = helpers::default_phase_order;
+
+  // NOLINTNEXTLINE(google-runtime-references)
+  void pup(PUP::er& /*p*/) {}
+};
+
+}  // namespace
+
+extern "C" void CkRegisterMainModule() {
+  Parallel::charmxx::register_main_module<Metavariables>();
+  Parallel::charmxx::register_init_node_and_proc({}, {});
+}
diff --git a/tests/Unit/ParallelAlgorithms/LinearSolver/Gmres/Test_ComplexGmresAlgorithm.yaml b/tests/Unit/ParallelAlgorithms/LinearSolver/Gmres/Test_ComplexGmresAlgorithm.yaml
new file mode 100644
index 0000000000000..6d5b01dacc0c4
--- /dev/null
+++ b/tests/Unit/ParallelAlgorithms/LinearSolver/Gmres/Test_ComplexGmresAlgorithm.yaml
@@ -0,0 +1,29 @@
+# Distributed under the MIT License.
+# See LICENSE.txt for details.
+
+---
+---
+
+LinearOperator:
+  - [[1, 2], [2, -1]]
+  - [[3, 4], [4, 1]]
+Source: [[1, 1], [2, -3]]
+InitialGuess: [[0, 0], [0, 0]]
+ExpectedResult: [[0.45, -1.4], [-1.2, 0.15]]
+
+Observers:
+  VolumeFileName: "Test_ComplexGmresAlgorithm_Volume"
+  ReductionFileName: "Test_ComplexGmresAlgorithm_Reductions"
+
+SerialGmres:
+  ConvergenceCriteria:
+    MaxIterations: 2
+    AbsoluteResidual: 1e-14
+    RelativeResidual: 0
+  Verbosity: Verbose
+
+ConvergenceReason: AbsoluteResidual
+
+ResourceInfo:
+  AvoidGlobalProc0: false
+  Singletons: Auto
diff --git a/tests/Unit/ParallelAlgorithms/LinearSolver/Gmres/Test_ElementActions.cpp b/tests/Unit/ParallelAlgorithms/LinearSolver/Gmres/Test_ElementActions.cpp
index a74d82af93231..94104306c296f 100644
--- a/tests/Unit/ParallelAlgorithms/LinearSolver/Gmres/Test_ElementActions.cpp
+++ b/tests/Unit/ParallelAlgorithms/LinearSolver/Gmres/Test_ElementActions.cpp
@@ -206,9 +206,9 @@ void test_element_actions() {
         REQUIRE_FALSE(ActionTesting::next_action_if_ready<element_array>(
             make_not_null(&runner), 0));
         auto& inbox = ActionTesting::get_inbox_tag<
-            element_array, LinearSolver::gmres::detail::Tags::
-                               FinalOrthogonalization<DummyOptionsGroup>>(
-            make_not_null(&runner), 0);
+            element_array,
+            LinearSolver::gmres::detail::Tags::FinalOrthogonalization<
+                DummyOptionsGroup, double>>(make_not_null(&runner), 0);
         const double normalization = 4.;
         const blaze::DynamicVector<double> minres{2., 4.};
         CAPTURE(has_converged);
diff --git a/tests/Unit/ParallelAlgorithms/LinearSolver/Gmres/Test_GmresAlgorithm.cpp b/tests/Unit/ParallelAlgorithms/LinearSolver/Gmres/Test_GmresAlgorithm.cpp
index 6b9f5009253f1..ce1ef47f92e3e 100644
--- a/tests/Unit/ParallelAlgorithms/LinearSolver/Gmres/Test_GmresAlgorithm.cpp
+++ b/tests/Unit/ParallelAlgorithms/LinearSolver/Gmres/Test_GmresAlgorithm.cpp
@@ -26,7 +26,7 @@ struct Metavariables {
       "Test the GMRES linear solver algorithm"};
 
   using linear_solver =
-      LinearSolver::gmres::Gmres<Metavariables, helpers::fields_tag,
+      LinearSolver::gmres::Gmres<Metavariables, helpers::fields_tag<double>,
                                  SerialGmres, false>;
   using preconditioner = void;
 
diff --git a/tests/Unit/ParallelAlgorithms/LinearSolver/Gmres/Test_GmresPreconditionedAlgorithm.cpp b/tests/Unit/ParallelAlgorithms/LinearSolver/Gmres/Test_GmresPreconditionedAlgorithm.cpp
index 4ce042f7cf928..35c8fcc26a2b4 100644
--- a/tests/Unit/ParallelAlgorithms/LinearSolver/Gmres/Test_GmresPreconditionedAlgorithm.cpp
+++ b/tests/Unit/ParallelAlgorithms/LinearSolver/Gmres/Test_GmresPreconditionedAlgorithm.cpp
@@ -31,7 +31,7 @@ struct Metavariables {
       "Test the preconditioned GMRES linear solver algorithm"};
 
   using linear_solver =
-      LinearSolver::gmres::Gmres<Metavariables, helpers::fields_tag,
+      LinearSolver::gmres::Gmres<Metavariables, helpers::fields_tag<double>,
                                  SerialGmres, true>;
   using preconditioner = LinearSolver::Richardson::Richardson<
       typename linear_solver::operand_tag, Preconditioner,
diff --git a/tests/Unit/ParallelAlgorithms/LinearSolver/Gmres/Test_ResidualMonitorActions.cpp b/tests/Unit/ParallelAlgorithms/LinearSolver/Gmres/Test_ResidualMonitorActions.cpp
index 458da2a09656c..c12566b671664 100644
--- a/tests/Unit/ParallelAlgorithms/LinearSolver/Gmres/Test_ResidualMonitorActions.cpp
+++ b/tests/Unit/ParallelAlgorithms/LinearSolver/Gmres/Test_ResidualMonitorActions.cpp
@@ -79,12 +79,13 @@ struct MockElementArray {
   using array_index = int;
   using phase_dependent_action_list = tmpl::list<
       Parallel::PhaseActions<Parallel::Phase::Initialization, tmpl::list<>>>;
-  using inbox_tags = tmpl::list<
-      LinearSolver::gmres::detail::Tags::InitialOrthogonalization<
-          TestLinearSolver>,
-      LinearSolver::gmres::detail::Tags::Orthogonalization<TestLinearSolver>,
-      LinearSolver::gmres::detail::Tags::FinalOrthogonalization<
-          TestLinearSolver>>;
+  using inbox_tags =
+      tmpl::list<LinearSolver::gmres::detail::Tags::InitialOrthogonalization<
+                     TestLinearSolver>,
+                 LinearSolver::gmres::detail::Tags::Orthogonalization<
+                     TestLinearSolver, double>,
+                 LinearSolver::gmres::detail::Tags::FinalOrthogonalization<
+                     TestLinearSolver, double>>;
 };
 
 struct Metavariables {
@@ -206,7 +207,7 @@ SPECTRE_TEST_CASE(
     // Test element state
     CHECK(get_element_inbox_tag(
               LinearSolver::gmres::detail::Tags::Orthogonalization<
-                  TestLinearSolver>{})
+                  TestLinearSolver, double>{})
               .at(1) == 2.);
   }
 
@@ -239,7 +240,7 @@ SPECTRE_TEST_CASE(
     const auto& element_inbox =
         get_element_inbox_tag(
             LinearSolver::gmres::detail::Tags::FinalOrthogonalization<
-                TestLinearSolver>{})
+                TestLinearSolver, double>{})
             .at(1);
     // beta = [2., 0.]
     // minres = inv(qr_R(H)) * trans(qr_Q(H)) * beta = [0.4615384615384615]
@@ -286,7 +287,7 @@ SPECTRE_TEST_CASE(
     const auto& element_inbox =
         get_element_inbox_tag(
             LinearSolver::gmres::detail::Tags::FinalOrthogonalization<
-                TestLinearSolver>{})
+                TestLinearSolver, double>{})
             .at(1);
     // beta = [2., 0.]
     // minres = inv(qr_R(H)) * trans(qr_Q(H)) * beta = [2.]
@@ -335,7 +336,7 @@ SPECTRE_TEST_CASE(
     const auto& element_inbox =
         get_element_inbox_tag(
             LinearSolver::gmres::detail::Tags::FinalOrthogonalization<
-                TestLinearSolver>{})
+                TestLinearSolver, double>{})
             .at(2);
     // beta = [1., 0., 0.]
     // minres = inv(qr_R(H)) * trans(qr_Q(H)) * beta = [0.13178295, 0.03100775]
@@ -373,7 +374,7 @@ SPECTRE_TEST_CASE(
     const auto& element_inbox =
         get_element_inbox_tag(
             LinearSolver::gmres::detail::Tags::FinalOrthogonalization<
-                TestLinearSolver>{})
+                TestLinearSolver, double>{})
             .at(1);
     // beta = [2., 0.]
     // minres = inv(qr_R(H)) * trans(qr_Q(H)) * beta = [0.6]
@@ -411,7 +412,7 @@ SPECTRE_TEST_CASE(
     const auto& element_inbox =
         get_element_inbox_tag(
             LinearSolver::gmres::detail::Tags::FinalOrthogonalization<
-                TestLinearSolver>{})
+                TestLinearSolver, double>{})
             .at(1);
     const auto& minres = get<1>(element_inbox);
     CHECK(minres.size() == 1);
diff --git a/tests/Unit/ParallelAlgorithms/LinearSolver/Richardson/Test_RichardsonAlgorithm.cpp b/tests/Unit/ParallelAlgorithms/LinearSolver/Richardson/Test_RichardsonAlgorithm.cpp
index 4fc5237fae009..4f5cd79477ef1 100644
--- a/tests/Unit/ParallelAlgorithms/LinearSolver/Richardson/Test_RichardsonAlgorithm.cpp
+++ b/tests/Unit/ParallelAlgorithms/LinearSolver/Richardson/Test_RichardsonAlgorithm.cpp
@@ -26,7 +26,7 @@ struct Metavariables {
       "Test the Richardson linear solver algorithm"};
 
   using linear_solver =
-      LinearSolver::Richardson::Richardson<helpers::fields_tag,
+      LinearSolver::Richardson::Richardson<helpers::fields_tag<double>,
                                            SerialRichardson>;
   using preconditioner = void;
 
diff --git a/tests/Unit/ParallelAlgorithms/NonlinearSolver/NewtonRaphson/Test_NewtonRaphsonAlgorithm.cpp b/tests/Unit/ParallelAlgorithms/NonlinearSolver/NewtonRaphson/Test_NewtonRaphsonAlgorithm.cpp
index 5236ff2d89bb5..93144f42b5d27 100644
--- a/tests/Unit/ParallelAlgorithms/NonlinearSolver/NewtonRaphson/Test_NewtonRaphsonAlgorithm.cpp
+++ b/tests/Unit/ParallelAlgorithms/NonlinearSolver/NewtonRaphson/Test_NewtonRaphsonAlgorithm.cpp
@@ -71,7 +71,7 @@ struct Metavariables {
       "Test the Newton-Raphson nonlinear solver algorithm"};
 
   using nonlinear_solver = NonlinearSolver::newton_raphson::NewtonRaphson<
-      Metavariables, helpers::fields_tag, NonlinearSolverGroup>;
+      Metavariables, helpers::fields_tag<double>, NonlinearSolverGroup>;
   using linear_solver = LinearSolver::gmres::Gmres<
       Metavariables, typename nonlinear_solver::linear_solver_fields_tag,
       LinearSolverGroup, false,