ginkgo-project · MarcelKoch · Jul 16, 2025 · Jul 9, 2025 · Jul 11, 2025
diff --git a/benchmark/blas/blas.cpp b/benchmark/blas/blas.cpp
@@ -1,4 +1,4 @@
-// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors
+// SPDX-FileCopyrightText: 2017 - 2025 The Ginkgo authors
 //
 // SPDX-License-Identifier: BSD-3-Clause
 
@@ -32,6 +32,12 @@ std::map<std::string, std::function<std::unique_ptr<BenchmarkOperation>(
                  exec, Generator{}, dims.n, dims.r, dims.stride_x,
                  dims.stride_y, false);
          }},
+        {"sub_scaled",
+         [](std::shared_ptr<const gko::Executor> exec, dimensions dims) {
+             return std::make_unique<SubScaledOperation<Generator>>(
+                 exec, Generator{}, dims.n, dims.r, dims.stride_x,
+                 dims.stride_y, false);
+         }},
         {"multiaxpy",
          [](std::shared_ptr<const gko::Executor> exec, dimensions dims) {
              return std::make_unique<AxpyOperation<Generator>>(

diff --git a/benchmark/blas/blas_common.hpp b/benchmark/blas/blas_common.hpp
@@ -28,6 +28,7 @@ DEFINE_string(
     "BLAS algorithms:\n"
     "   copy (y = x),\n"
     "   axpy (y = y + a * x),\n"
+    "   sub_scaled (y = y - a * x),\n"
     "   multiaxpy (like axpy, but a has one entry per column),\n"
     "   scal (y = a * y),\n"
     "   multiscal (like scal, but a has one entry per column),\n"
@@ -136,6 +137,46 @@ class AxpyOperation : public BenchmarkOperation {
 };
 
 
+template <typename Generator>
+class SubScaledOperation : public BenchmarkOperation {
+public:
+    SubScaledOperation(std::shared_ptr<const gko::Executor> exec,
+                       const Generator& generator, gko::size_type rows,
+                       gko::size_type cols, gko::size_type stride_in,
+                       gko::size_type stride_out, bool multi)
+    {
+        auto size = gko::dim<2>{rows, cols};
+        alpha_ = gko::matrix::Dense<etype>::create(
+            exec, gko::dim<2>{1, multi ? cols : 1});
+        x_ = generator.create_multi_vector_strided(
+            exec, size, generator.create_default_local_size(size), stride_in);
+        y_ = generator.create_multi_vector_strided(
+            exec, size, generator.create_default_local_size(size), stride_out);
+        alpha_->fill(1);
+        as_vector<Generator>(x_)->fill(1);
+    }
+
+    gko::size_type get_flops() const override
+    {
+        return y_->get_size()[0] * y_->get_size()[1] * 2;
+    }
+
+    gko::size_type get_memory() const override
+    {
+        return y_->get_size()[0] * y_->get_size()[1] * sizeof(etype) * 3;
+    }
+
+    void prepare() override { as_vector<Generator>(y_)->fill(1); }
+
+    void run() override { as_vector<Generator>(y_)->sub_scaled(alpha_, x_); }
+
+private:
+    std::unique_ptr<gko::matrix::Dense<etype>> alpha_;
+    std::unique_ptr<gko::LinOp> x_;
+    std::unique_ptr<gko::LinOp> y_;
+};
+
+
 template <typename Generator>
 class ScalOperation : public BenchmarkOperation {
 public:

diff --git a/benchmark/solver/solver_common.hpp b/benchmark/solver/solver_common.hpp
@@ -53,6 +53,9 @@ DEFINE_uint32(gcr_restart, 100,
 DEFINE_uint32(gmres_restart, 100,
               "Maximum dimension of the Krylov space to use in GMRES");
 
+DEFINE_string(gmres_ortho_method, "mgs",
+              "The orthogonalization method to use in GMRES.");
+
 DEFINE_uint32(idr_subspace_dim, 2,
               "What dimension of the subspace to use in IDR");
 
@@ -200,9 +203,23 @@ std::unique_ptr<gko::LinOpFactory> generate_solver(
                 .with_kappa(static_cast<rc_etype>(FLAGS_idr_kappa)),
             exec, precond, max_iters);
     } else if (description == "gmres") {
+        gko::solver::gmres::ortho_method ortho_method;
+        if (FLAGS_gmres_ortho_method == "mgs") {
+            ortho_method = gko::solver::gmres::ortho_method::mgs;
+        } else if (FLAGS_gmres_ortho_method == "cgs") {
+            ortho_method = gko::solver::gmres::ortho_method::cgs;
+        } else if (FLAGS_gmres_ortho_method == "cgs2") {
+            ortho_method = gko::solver::gmres::ortho_method::cgs2;
+        } else {
+            throw std::range_error(
+                std::string(
+                    "GMRES doesn't support the orthogonalization method <") +
+                FLAGS_gmres_ortho_method + ">!");
+        }
         return add_criteria_precond_finalize(
-            gko::solver::Gmres<etype>::build().with_krylov_dim(
-                FLAGS_gmres_restart),
+            gko::solver::Gmres<etype>::build()
+                .with_krylov_dim(FLAGS_gmres_restart)
+                .with_ortho_method(ortho_method),
             exec, precond, max_iters);
     } else if (description == "minres") {
         return add_criteria_precond_finalize<gko::solver::Minres<etype>>(