diff --git a/core/CMakeLists.txt b/core/CMakeLists.txt
index f5388785d6a..55c7696c8c5 100644
--- a/core/CMakeLists.txt
+++ b/core/CMakeLists.txt
@@ -9,6 +9,7 @@ target_sources(ginkgo
     base/executor.cpp
     base/mtx_io.cpp
     base/perturbation.cpp
+    base/utils.cpp
     base/version.cpp
     factorization/ic.cpp
     factorization/ilu.cpp
diff --git a/core/base/utils.cpp b/core/base/utils.cpp
new file mode 100644
index 00000000000..642f20a7ff3
--- /dev/null
+++ b/core/base/utils.cpp
@@ -0,0 +1,63 @@
+/*******************************<GINKGO LICENSE>******************************
+Copyright (c) 2017-2021, the Ginkgo authors
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions
+are met:
+
+1. Redistributions of source code must retain the above copyright
+notice, this list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright
+notice, this list of conditions and the following disclaimer in the
+documentation and/or other materials provided with the distribution.
+
+3. Neither the name of the copyright holder nor the names of its
+contributors may be used to endorse or promote products derived from
+this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
+IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
+PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+******************************<GINKGO LICENSE>*******************************/
+
+#ifndef GKO_CORE_BASE_UTILS_HPP_
+#define GKO_CORE_BASE_UTILS_HPP_
+
+
+#include "utils.hpp"
+
+
+#include <memory>
+
+
+#include <ginkgo/core/base/executor.hpp>
+#include <ginkgo/core/matrix/dense.hpp>
+
+
+namespace gko {
+
+template <typename ValueType>
+std::unique_ptr<matrix::Dense<ValueType>> scalar_to_dense(
+    const ValueType val, std::shared_ptr<const Executor> exec)
+{
+    return initialize<matrix::Dense<ValueType>>({val}, exec);
+}
+#define GKO_DECLARE_SCALAR_TO_DENSE(_type)                 \
+    std::unique_ptr<matrix::Dense<_type>> scalar_to_dense( \
+        const _type, std::shared_ptr<const Executor>)
+GKO_INSTANTIATE_FOR_EACH_VALUE_TYPE(GKO_DECLARE_SCALAR_TO_DENSE);
+
+}  // namespace gko
+
+
+#endif  // GKO_CORE_BASE_UTILS_HPP_
diff --git a/core/test/base/CMakeLists.txt b/core/test/base/CMakeLists.txt
index 6f8c7291165..551e7985720 100644
--- a/core/test/base/CMakeLists.txt
+++ b/core/test/base/CMakeLists.txt
@@ -21,4 +21,5 @@ ginkgo_create_test(range_accessors)
 ginkgo_create_thread_test(sanitizers)
 ginkgo_create_test(types)
 ginkgo_create_test(utils)
+ginkgo_create_test(value_typed_lin_op)
 ginkgo_create_test(version)
diff --git a/core/test/base/value_typed_lin_op.cpp b/core/test/base/value_typed_lin_op.cpp
new file mode 100644
index 00000000000..9db3f5edb3c
--- /dev/null
+++ b/core/test/base/value_typed_lin_op.cpp
@@ -0,0 +1,176 @@
+/*******************************<GINKGO LICENSE>******************************
+Copyright (c) 2017-2021, the Ginkgo authors
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions
+are met:
+
+1. Redistributions of source code must retain the above copyright
+notice, this list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright
+notice, this list of conditions and the following disclaimer in the
+documentation and/or other materials provided with the distribution.
+
+3. Neither the name of the copyright holder nor the names of its
+contributors may be used to endorse or promote products derived from
+this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
+IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
+PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+******************************<GINKGO LICENSE>*******************************/
+
+#include <ginkgo/core/base/lin_op.hpp>
+
+
+#include <complex>
+#include <memory>
+#include <type_traits>
+
+
+#include <gtest/gtest.h>
+
+
+#include <ginkgo/core/base/math.hpp>
+#include <ginkgo/core/matrix/dense.hpp>
+
+
+#include "core/test/utils.hpp"
+
+
+namespace {
+
+
+template <typename T>
+T get_value_from_linop(const gko::LinOp *val);
+
+
+template <typename T>
+class DummyValueTypedLinOp
+    : public gko::EnableValueTypedLinOp<DummyValueTypedLinOp<T>, T>,
+      public gko::EnableCreateMethod<DummyValueTypedLinOp<T>> {
+    using Self = DummyValueTypedLinOp<T>;
+
+public:
+    explicit DummyValueTypedLinOp(std::shared_ptr<const gko::Executor> exec,
+                                  gko::dim<2> size = gko::dim<2>{})
+        : gko::EnableValueTypedLinOp<Self, T>(exec, size), value_()
+    {}
+    DummyValueTypedLinOp(std::shared_ptr<const gko::Executor> exec,
+                         gko::dim<2> size, T value)
+        : gko::EnableValueTypedLinOp<Self, T>(exec, size), value_(value)
+    {}
+
+    T get_value() const { return value_; }
+
+protected:
+    void apply_impl(const gko::LinOp *b, gko::LinOp *x) const override {}
+
+    void apply_impl(const gko::LinOp *alpha, const gko::LinOp *b,
+                    const gko::LinOp *beta, gko::LinOp *x) const override
+    {
+        T alpha_v = get_value_from_linop<T>(alpha);
+        T beta_v = get_value_from_linop<T>(beta);
+
+        gko::as<Self>(x)->value_ =
+            alpha_v * gko::as<Self>(this)->value_ * gko::as<Self>(b)->value_ +
+            beta_v * gko::as<Self>(x)->value_;
+    }
+
+    T value_;
+};
+
+
+template <typename T>
+T get_value_from_linop(const gko::LinOp *val)
+{
+    if (auto *dense = dynamic_cast<const gko::matrix::Dense<T> *>(val)) {
+        return dense->at(0, 0);
+    } else {
+        return (dynamic_cast<const DummyValueTypedLinOp<T> *>(val))
+            ->get_value();
+    }
+}
+
+
+template <typename T>
+class EnableValueTypedLinOp : public ::testing::Test {
+protected:
+    using dummy_type = DummyValueTypedLinOp<T>;
+    using value_type = T;
+
+    gko::dim<2> dim{1, 1};
+    std::shared_ptr<const gko::ReferenceExecutor> ref{
+        gko::ReferenceExecutor::create()};
+    std::unique_ptr<dummy_type> op{dummy_type::create(ref, dim, T{1.0})};
+    T alpha_v{2.0};
+    std::unique_ptr<dummy_type> alpha{dummy_type::create(ref, dim, alpha_v)};
+    T beta_v{3.0};
+    std::unique_ptr<dummy_type> beta{dummy_type::create(ref, dim, beta_v)};
+    std::unique_ptr<dummy_type> b{dummy_type::create(ref, dim, T{4.0})};
+    std::unique_ptr<dummy_type> x{dummy_type::create(ref, dim, T{5.0})};
+};
+
+
+TYPED_TEST_SUITE(EnableValueTypedLinOp, gko::test::ValueTypes);
+
+
+TYPED_TEST(EnableValueTypedLinOp, CanCallExtendedApplyImplLinopLinop)
+{
+    using value_type = typename TestFixture::value_type;
+    this->op->apply(lend(this->alpha), lend(this->b), lend(this->beta),
+                    lend(this->x));
+
+    ASSERT_EQ(this->x->get_value(), value_type{23.0});
+}
+
+
+TYPED_TEST(EnableValueTypedLinOp, CanCallExtendedApplyImplValueValue)
+{
+    auto reference_result = gko::clone(this->x);
+
+    this->op->apply(this->alpha_v, lend(this->b), this->beta_v, lend(this->x));
+    this->op->apply(lend(this->alpha), lend(this->b), lend(this->beta),
+                    lend(reference_result));
+
+    ASSERT_EQ(reference_result->get_value(), this->x->get_value());
+}
+
+
+TYPED_TEST(EnableValueTypedLinOp, CanCallExtendedApplyImplValueLinop)
+{
+    auto reference_result = gko::clone(this->x);
+
+    this->op->apply(this->alpha_v, lend(this->b), lend(this->beta),
+                    lend(this->x));
+    this->op->apply(lend(this->alpha), lend(this->b), lend(this->beta),
+                    lend(reference_result));
+
+    ASSERT_EQ(reference_result->get_value(), this->x->get_value());
+}
+
+
+TYPED_TEST(EnableValueTypedLinOp, CanCallExtendedApplyImplLinopValue)
+{
+    auto reference_result = gko::clone(this->x);
+
+    this->op->apply(lend(this->alpha), lend(this->b), this->beta_v,
+                    lend(this->x));
+    this->op->apply(lend(this->alpha), lend(this->b), lend(this->beta),
+                    lend(reference_result));
+
+    ASSERT_EQ(reference_result->get_value(), this->x->get_value());
+}
+
+
+}  // namespace
diff --git a/examples/adaptiveprecision-blockjacobi/adaptiveprecision-blockjacobi.cpp b/examples/adaptiveprecision-blockjacobi/adaptiveprecision-blockjacobi.cpp
index b3e672e6c1f..9dcee9b69b4 100644
--- a/examples/adaptiveprecision-blockjacobi/adaptiveprecision-blockjacobi.cpp
+++ b/examples/adaptiveprecision-blockjacobi/adaptiveprecision-blockjacobi.cpp
@@ -100,10 +100,8 @@ int main(int argc, char *argv[])
     b->copy_from(host_x.get());
 
     // Calculate initial residual by overwriting b
-    auto one = gko::initialize<vec>({1.0}, exec);
-    auto neg_one = gko::initialize<vec>({-1.0}, exec);
-    auto initres = gko::initialize<real_vec>({0.0}, exec);
-    A->apply(lend(one), lend(x), lend(neg_one), lend(b));
+    auto initres = gko::initialize<real_vec>({0.0}, exec->get_master());
+    A->apply(1.0, lend(x), -1.0, lend(b));
     b->compute_norm2(lend(initres));
 
     // copy b again
@@ -146,8 +144,8 @@ int main(int argc, char *argv[])
     time += std::chrono::duration_cast<std::chrono::nanoseconds>(toc - tic);
 
     // Calculate residual
-    auto res = gko::initialize<real_vec>({0.0}, exec);
-    A->apply(lend(one), lend(x), lend(neg_one), lend(b));
+    auto res = gko::initialize<real_vec>({0.0}, exec->get_master());
+    A->apply(1.0, lend(x), -1.0, lend(b));
     b->compute_norm2(lend(res));
     auto impl_res = gko::as<real_vec>(logger->get_implicit_sq_resnorm());
 
diff --git a/examples/cb-gmres/cb-gmres.cpp b/examples/cb-gmres/cb-gmres.cpp
index ffb52579132..d1ea569c355 100644
--- a/examples/cb-gmres/cb-gmres.cpp
+++ b/examples/cb-gmres/cb-gmres.cpp
@@ -138,7 +138,7 @@ int main(int argc, char *argv[])
         b_host->at(i, 0) =
             ValueType{1} / std::sqrt(static_cast<ValueType>(A_size[0]));
     }
-    auto b_norm = gko::initialize<real_vec>({0.0}, exec);
+    auto b_norm = gko::initialize<real_vec>({0.0}, exec->get_master());
     b_host->compute_norm2(lend(b_norm));
     auto b = clone(exec, lend(b_host));
 
@@ -194,25 +194,20 @@ int main(int argc, char *argv[])
 
     // To measure if your solution has actually converged, the error of the
     // solution is measured.
-    // one, neg_one are objects that represent the numbers which allow for a
-    // uniform interface when computing on any device. To compute the residual,
-    // the (advanced) apply method is used.
-    auto one = gko::initialize<vec>({1.0}, exec);
-    auto neg_one = gko::initialize<vec>({-1.0}, exec);
-
-    auto res_norm_keep = gko::initialize<real_vec>({0.0}, exec);
-    auto res_norm_reduce = gko::initialize<real_vec>({0.0}, exec);
+    // To compute the residual, the (advanced) apply method is used.
+    auto res_norm_keep = gko::initialize<real_vec>({0.0}, exec->get_master());
+    auto res_norm_reduce = gko::initialize<real_vec>({0.0}, exec->get_master());
     auto tmp = gko::clone(gko::lend(b));
 
     // tmp = Ax - tmp
-    A->apply(lend(one), lend(x_keep), lend(neg_one), lend(tmp));
+    A->apply(1.0, lend(x_keep), -1.0, lend(tmp));
     tmp->compute_norm2(lend(res_norm_keep));
 
     std::cout << "\nResidual norm without compression:\n";
     write(std::cout, lend(res_norm_keep));
 
     tmp->copy_from(lend(b));
-    A->apply(lend(one), lend(x_reduce), lend(neg_one), lend(tmp));
+    A->apply(1.0, lend(x_reduce), -1.0, lend(tmp));
     tmp->compute_norm2(lend(res_norm_reduce));
 
     std::cout << "\nResidual norm with compression:\n";
diff --git a/examples/custom-logger/custom-logger.cpp b/examples/custom-logger/custom-logger.cpp
index 39fcb440eb1..77b1d930afa 100644
--- a/examples/custom-logger/custom-logger.cpp
+++ b/examples/custom-logger/custom-logger.cpp
@@ -343,18 +343,14 @@ int main(int argc, char *argv[])
     logger->write();
 
     // To measure if your solution has actually converged, you can measure
-    // the error of the solution. one, neg_one are objects that represent
-    // the numbers which allow for a uniform interface when computing on any
-    // device. To compute the residual, all you need to do is call the apply
-    // method, which in this case is an spmv and equivalent to the LAPACK
-    // z_spmv routine. Finally, you compute the euclidean 2-norm with the
-    // compute_norm2 function.
-    auto one = gko::initialize<vec>({1.0}, exec);
-    auto neg_one = gko::initialize<vec>({-1.0}, exec);
-    auto res = gko::initialize<real_vec>({0.0}, exec);
-    A->apply(gko::lend(one), gko::lend(x), gko::lend(neg_one), gko::lend(b));
-    b->compute_norm2(gko::lend(res));
+    // the error of the solution. To compute the residual, all you need to do is
+    // call the apply method, which in this case is an spmv and equivalent to
+    // the LAPACK z_spmv routine. Finally, you compute the euclidean 2-norm with
+    // the compute_norm2 function.
+    auto res = gko::initialize<real_vec>({0.0}, exec->get_master());
+    A->apply(1.0, gko::lend(x), -1.0, gko::lend(b));
+    b->compute_norm2(lend(res));
 
     std::cout << "Residual norm sqrt(r^T r):\n";
-    write(std::cout, gko::lend(res));
+    write(std::cout, lend(res));
 }
diff --git a/examples/custom-stopping-criterion/custom-stopping-criterion.cpp b/examples/custom-stopping-criterion/custom-stopping-criterion.cpp
index 7c5074009bc..f209d3611b9 100644
--- a/examples/custom-stopping-criterion/custom-stopping-criterion.cpp
+++ b/examples/custom-stopping-criterion/custom-stopping-criterion.cpp
@@ -127,10 +127,8 @@ void run_solver(volatile bool *stop_iteration_process,
     write(std::cout, lend(x));
 
     // Calculate residual
-    auto one = gko::initialize<vec>({1.0}, exec);
-    auto neg_one = gko::initialize<vec>({-1.0}, exec);
-    auto res = gko::initialize<real_vec>({0.0}, exec);
-    A->apply(lend(one), lend(x), lend(neg_one), lend(b));
+    auto res = gko::initialize<real_vec>({0.0}, exec->get_master());
+    A->apply(1.0, lend(x), -1.0, lend(b));
     b->compute_norm2(lend(res));
 
     std::cout << "Residual norm sqrt(r^T r): \n";
diff --git a/examples/heat-equation/heat-equation.cpp b/examples/heat-equation/heat-equation.cpp
index 60ca7f53eb3..212ee45a9e3 100644
--- a/examples/heat-equation/heat-equation.cpp
+++ b/examples/heat-equation/heat-equation.cpp
@@ -158,7 +158,7 @@ int main(int argc, char *argv[])
     // create output vector with initial guess for
     auto out_vector = in_vector->clone();
     // create scalar for source update
-    auto tau_source_scalar = gko::initialize<vec>({source_scale * tau}, exec);
+    auto tau_source_scalar = source_scale * tau;
     // create stencil matrix as shared_ptr for solver
     auto stencil_matrix = gko::share(mtx::create(exec));
     // assemble matrix
@@ -212,7 +212,7 @@ int main(int argc, char *argv[])
                     ->get_const_values());
         }
         // add heat source contribution
-        in_vector->add_scaled(gko::lend(tau_source_scalar), gko::lend(source));
+        in_vector->add_scaled(tau_source_scalar, gko::lend(source));
         // execute Euler step
         solver->apply(gko::lend(in_vector), gko::lend(out_vector));
         // swap input and output
diff --git a/examples/ilu-preconditioned-solver/ilu-preconditioned-solver.cpp b/examples/ilu-preconditioned-solver/ilu-preconditioned-solver.cpp
index 1b9e7b5f53b..559a95d9a7b 100644
--- a/examples/ilu-preconditioned-solver/ilu-preconditioned-solver.cpp
+++ b/examples/ilu-preconditioned-solver/ilu-preconditioned-solver.cpp
@@ -134,12 +134,10 @@ int main(int argc, char *argv[])
     write(std::cout, gko::lend(x));
 
     // Calculate residual
-    auto one = gko::initialize<vec>({1.0}, exec);
-    auto neg_one = gko::initialize<vec>({-1.0}, exec);
-    auto res = gko::initialize<real_vec>({0.0}, exec);
-    A->apply(gko::lend(one), gko::lend(x), gko::lend(neg_one), gko::lend(b));
-    b->compute_norm2(gko::lend(res));
+    auto res = gko::initialize<real_vec>({0.0}, exec->get_master());
+    A->apply(1.0, gko::lend(x), -1.0, gko::lend(b));
+    b->compute_norm2(lend(res));
 
     std::cout << "Residual norm sqrt(r^T r):\n";
-    write(std::cout, gko::lend(res));
+    write(std::cout, lend(res));
 }
diff --git a/examples/ir-ilu-preconditioned-solver/ir-ilu-preconditioned-solver.cpp b/examples/ir-ilu-preconditioned-solver/ir-ilu-preconditioned-solver.cpp
index aab893c6efa..458d16adf2b 100644
--- a/examples/ir-ilu-preconditioned-solver/ir-ilu-preconditioned-solver.cpp
+++ b/examples/ir-ilu-preconditioned-solver/ir-ilu-preconditioned-solver.cpp
@@ -185,16 +185,14 @@ int main(int argc, char *argv[])
     write(std::cout, gko::lend(x));
 
     // Calculate residual
-    auto one = gko::initialize<vec>({1.0}, exec);
-    auto neg_one = gko::initialize<vec>({-1.0}, exec);
-    auto res = gko::initialize<real_vec>({0.0}, exec);
-    A->apply(gko::lend(one), gko::lend(x), gko::lend(neg_one), gko::lend(b));
-    b->compute_norm2(gko::lend(res));
+    auto res = gko::initialize<real_vec>({0.0}, exec->get_master());
+    A->apply(1.0, gko::lend(x), -1.0, gko::lend(b));
+    b->compute_norm2(lend(res));
 
     std::cout << "GMRES iteration count:     " << logger->get_num_iterations()
               << "\n";
     std::cout << "GMRES execution time [ms]: "
               << static_cast<double>(time.count()) / 100000000.0 << "\n";
     std::cout << "Residual norm sqrt(r^T r):\n";
-    write(std::cout, gko::lend(res));
+    write(std::cout, lend(res));
 }
diff --git a/examples/iterative-refinement/iterative-refinement.cpp b/examples/iterative-refinement/iterative-refinement.cpp
index 20a3514316e..b26045c0d74 100644
--- a/examples/iterative-refinement/iterative-refinement.cpp
+++ b/examples/iterative-refinement/iterative-refinement.cpp
@@ -101,10 +101,8 @@ int main(int argc, char *argv[])
     b->copy_from(host_x.get());
 
     // Calculate initial residual by overwriting b
-    auto one = gko::initialize<vec>({1.0}, exec);
-    auto neg_one = gko::initialize<vec>({-1.0}, exec);
-    auto initres = gko::initialize<real_vec>({0.0}, exec);
-    A->apply(lend(one), lend(x), lend(neg_one), lend(b));
+    auto initres = gko::initialize<real_vec>({0.0}, exec->get_master());
+    A->apply(1.0, lend(x), -1.0, lend(b));
     b->compute_norm2(lend(initres));
 
     // copy b again
@@ -148,8 +146,8 @@ int main(int argc, char *argv[])
     time += std::chrono::duration_cast<std::chrono::nanoseconds>(toc - tic);
 
     // Calculate residual
-    auto res = gko::initialize<real_vec>({0.0}, exec);
-    A->apply(lend(one), lend(x), lend(neg_one), lend(b));
+    auto res = gko::initialize<real_vec>({0.0}, exec->get_master());
+    A->apply(1.0, lend(x), -1.0, lend(b));
     b->compute_norm2(lend(res));
 
     std::cout << "Initial residual norm sqrt(r^T r):\n";
diff --git a/examples/mixed-precision-ir/mixed-precision-ir.cpp b/examples/mixed-precision-ir/mixed-precision-ir.cpp
index 499924fcfe9..58d6e6f8ad5 100644
--- a/examples/mixed-precision-ir/mixed-precision-ir.cpp
+++ b/examples/mixed-precision-ir/mixed-precision-ir.cpp
@@ -108,11 +108,9 @@ int main(int argc, char *argv[])
     b->copy_from(host_x.get());
 
     // Calculate initial residual by overwriting b
-    auto one = gko::initialize<vec>({1.0}, exec);
-    auto neg_one = gko::initialize<vec>({-1.0}, exec);
-    auto initres_vec = gko::initialize<real_vec>({0.0}, exec);
-    A->apply(lend(one), lend(x), lend(neg_one), lend(b));
-    b->compute_norm2(lend(initres_vec));
+    auto initres = gko::initialize<real_vec>({0.0}, exec->get_master());
+    A->apply(1.0, lend(x), -1.0, lend(b));
+    b->compute_norm2(lend(initres));
 
     // Build lower-precision system matrix and residual
     auto solver_A = solver_mtx::create(exec);
@@ -139,8 +137,7 @@ int main(int argc, char *argv[])
     // Solve system
     exec->synchronize();
     std::chrono::nanoseconds time(0);
-    auto res_vec = gko::initialize<real_vec>({0.0}, exec);
-    auto initres = exec->copy_val_to_host(initres_vec->get_const_values());
+    auto res = gko::initialize<real_vec>({0.0}, exec->get_master());
     auto inner_solution = solver_vec::create(exec);
     auto outer_delta = vec::create(exec);
     auto tic = std::chrono::steady_clock::now();
@@ -150,11 +147,11 @@ int main(int argc, char *argv[])
 
         // convert residual to inner precision
         outer_residual->convert_to(lend(inner_residual));
-        outer_residual->compute_norm2(lend(res_vec));
-        auto res = exec->copy_val_to_host(res_vec->get_const_values());
+        outer_residual->compute_norm2(lend(res));
 
         // break if we exceed the number of iterations or have converged
-        if (iter > max_outer_iters || res / initres < outer_reduction_factor) {
+        if (iter > max_outer_iters ||
+            res->at(0) / initres->at(0) < outer_reduction_factor) {
             break;
         }
 
@@ -168,24 +165,24 @@ int main(int argc, char *argv[])
         inner_solution->convert_to(lend(outer_delta));
 
         // x = x + inner_solution
-        x->add_scaled(lend(one), lend(outer_delta));
+        x->add_scaled(1.0, lend(outer_delta));
 
         // residual = b - A * x
         outer_residual->copy_from(lend(b));
-        A->apply(lend(neg_one), lend(x), lend(one), lend(outer_residual));
+        A->apply(-1.0, lend(x), 1.0, lend(outer_residual));
     }
 
     auto toc = std::chrono::steady_clock::now();
     time += std::chrono::duration_cast<std::chrono::nanoseconds>(toc - tic);
 
     // Calculate residual
-    A->apply(lend(one), lend(x), lend(neg_one), lend(b));
-    b->compute_norm2(lend(res_vec));
+    A->apply(1.0, lend(x), -1.0, lend(b));
+    b->compute_norm2(lend(res));
 
     std::cout << "Initial residual norm sqrt(r^T r):\n";
-    write(std::cout, lend(initres_vec));
+    write(std::cout, lend(initres));
     std::cout << "Final residual norm sqrt(r^T r):\n";
-    write(std::cout, lend(res_vec));
+    write(std::cout, lend(res));
 
     // Print solver statistics
     std::cout << "MPIR iteration count:     " << iter << std::endl;
diff --git a/examples/mixed-spmv/mixed-spmv.cpp b/examples/mixed-spmv/mixed-spmv.cpp
index 7a6e6a8b542..f709dd2a757 100644
--- a/examples/mixed-spmv/mixed-spmv.cpp
+++ b/examples/mixed-spmv/mixed-spmv.cpp
@@ -240,12 +240,10 @@ int main(int argc, char *argv[])
 
 
     // To measure error of result.
-    // neg_one is an object that represent the number -1.0 which allows for a
-    // uniform interface when computing on any device. To compute the residual,
-    // all you need to do is call the add_scaled method, which in this case is
-    // an axpy and equivalent to the LAPACK axpy routine. Finally, you compute
-    // the euclidean 2-norm with the compute_norm2 function.
-    auto neg_one = gko::initialize<hp_vec>({-1.0}, exec);
+    // To compute the residual, all you need to do is call the add_scaled
+    // method, which in this case is an axpy and equivalent to the LAPACK axpy
+    // routine. Finally, you compute the euclidean 2-norm with the compute_norm2
+    // function.
     auto hp_x_norm = gko::initialize<real_vec>({0.0}, exec->get_master());
     auto lp_diff_norm = gko::initialize<real_vec>({0.0}, exec->get_master());
     auto hplp_diff_norm = gko::initialize<real_vec>({0.0}, exec->get_master());
@@ -257,31 +255,31 @@ int main(int argc, char *argv[])
     auto lphp_diff = hp_x->clone();
 
     hp_x->compute_norm2(lend(hp_x_norm));
-    lp_diff->add_scaled(lend(neg_one), lend(lp_x));
+    lp_diff->add_scaled(-1.0, lend(lp_x));
     lp_diff->compute_norm2(lend(lp_diff_norm));
-    hplp_diff->add_scaled(lend(neg_one), lend(hplp_x));
+    hplp_diff->add_scaled(-1.0, lend(hplp_x));
     hplp_diff->compute_norm2(lend(hplp_diff_norm));
-    lplp_diff->add_scaled(lend(neg_one), lend(lplp_x));
+    lplp_diff->add_scaled(-1.0, lend(lplp_x));
     lplp_diff->compute_norm2(lend(lplp_diff_norm));
-    lphp_diff->add_scaled(lend(neg_one), lend(lphp_x));
+    lphp_diff->add_scaled(-1.0, lend(lphp_x));
     lphp_diff->compute_norm2(lend(lphp_diff_norm));
     exec->synchronize();
 
     std::cout.precision(10);
     std::cout << std::scientific;
     std::cout << "High Precision time(s): " << hp_sec << std::endl;
-    std::cout << "High Precision result norm: " << hp_x_norm->at(0)
+    std::cout << "High Precision result norm: " << hp_x_norm->get_values()[0]
               << std::endl;
     std::cout << "Low Precision time(s): " << lp_sec << std::endl;
     std::cout << "Low Precision relative error: "
-              << lp_diff_norm->at(0) / hp_x_norm->at(0) << "\n";
+              << lp_diff_norm->get_values()[0] / hp_x_norm->get_values()[0]
+              << "\n";
     std::cout << "Hp * Lp -> Hp time(s): " << hplp_sec << std::endl;
     std::cout << "Hp * Lp -> Hp relative error: "
-              << hplp_diff_norm->at(0) / hp_x_norm->at(0) << "\n";
+              << hplp_diff_norm->get_values()[0] / hp_x_norm->get_values()[0]
+              << "\n";
     std::cout << "Lp * Lp -> Hp time(s): " << lplp_sec << std::endl;
     std::cout << "Lp * Lp -> Hp relative error: "
-              << lplp_diff_norm->at(0) / hp_x_norm->at(0) << "\n";
-    std::cout << "Lp * Hp -> Hp time(s): " << lplp_sec << std::endl;
-    std::cout << "Lp * Hp -> Hp relative error: "
-              << lphp_diff_norm->at(0) / hp_x_norm->at(0) << "\n";
+              << lplp_diff_norm->get_values()[0] / hp_x_norm->get_values()[0]
+              << "\n";
 }
diff --git a/examples/papi-logging/papi-logging.cpp b/examples/papi-logging/papi-logging.cpp
index 6aff2772980..c05747d21e3 100644
--- a/examples/papi-logging/papi-logging.cpp
+++ b/examples/papi-logging/papi-logging.cpp
@@ -213,10 +213,8 @@ int main(int argc, char *argv[])
     write(std::cout, lend(x));
 
     // Calculate residual
-    auto one = gko::initialize<vec>({1.0}, exec);
-    auto neg_one = gko::initialize<vec>({-1.0}, exec);
-    auto res = gko::initialize<real_vec>({0.0}, exec);
-    A->apply(lend(one), lend(x), lend(neg_one), lend(b));
+    auto res = gko::initialize<real_vec>({0.0}, exec->get_master());
+    A->apply(1.0, lend(x), -1.0, lend(b));
     b->compute_norm2(lend(res));
 
     std::cout << "Residual norm sqrt(r^T r): \n";
diff --git a/examples/preconditioned-solver/preconditioned-solver.cpp b/examples/preconditioned-solver/preconditioned-solver.cpp
index ce01f0c0cc5..214977c864a 100644
--- a/examples/preconditioned-solver/preconditioned-solver.cpp
+++ b/examples/preconditioned-solver/preconditioned-solver.cpp
@@ -116,10 +116,8 @@ int main(int argc, char *argv[])
     write(std::cout, lend(x));
 
     // Calculate residual
-    auto one = gko::initialize<vec>({1.0}, exec);
-    auto neg_one = gko::initialize<vec>({-1.0}, exec);
-    auto res = gko::initialize<real_vec>({0.0}, exec);
-    A->apply(lend(one), lend(x), lend(neg_one), lend(b));
+    auto res = gko::initialize<real_vec>({0.0}, exec->get_master());
+    A->apply(1.0, lend(x), -1.0, lend(b));
     b->compute_norm2(lend(res));
 
     std::cout << "Residual norm sqrt(r^T r):\n";
diff --git a/examples/simple-solver-logging/simple-solver-logging.cpp b/examples/simple-solver-logging/simple-solver-logging.cpp
index a7f5bbaca84..3541f70b280 100644
--- a/examples/simple-solver-logging/simple-solver-logging.cpp
+++ b/examples/simple-solver-logging/simple-solver-logging.cpp
@@ -187,10 +187,8 @@ int main(int argc, char *argv[])
     write(std::cout, lend(x));
 
     // Calculate residual
-    auto one = gko::initialize<vec>({1.0}, exec);
-    auto neg_one = gko::initialize<vec>({-1.0}, exec);
-    auto res = gko::initialize<real_vec>({0.0}, exec);
-    A->apply(lend(one), lend(x), lend(neg_one), lend(b));
+    auto res = gko::initialize<real_vec>({0.0}, exec->get_master());
+    A->apply(1.0, lend(x), -1.0, lend(b));
     b->compute_norm2(lend(res));
 
     std::cout << "Residual norm sqrt(r^T r):\n";
diff --git a/examples/simple-solver/simple-solver.cpp b/examples/simple-solver/simple-solver.cpp
index bf8095fc373..d925a18d049 100644
--- a/examples/simple-solver/simple-solver.cpp
+++ b/examples/simple-solver/simple-solver.cpp
@@ -157,15 +157,12 @@ int main(int argc, char *argv[])
 
     // To measure if your solution has actually converged, you can measure the
     // error of the solution.
-    // one, neg_one are objects that represent the numbers which allow for a
-    // uniform interface when computing on any device. To compute the residual,
-    // all you need to do is call the apply method, which in this case is an
-    // spmv and equivalent to the LAPACK z_spmv routine. Finally, you compute
-    // the euclidean 2-norm with the compute_norm2 function.
-    auto one = gko::initialize<vec>({1.0}, exec);
-    auto neg_one = gko::initialize<vec>({-1.0}, exec);
-    auto res = gko::initialize<real_vec>({0.0}, exec);
-    A->apply(lend(one), lend(x), lend(neg_one), lend(b));
+    // To compute the residual, all you need to do is call the apply method,
+    // which in this case is an spmv and equivalent to the LAPACK z_spmv
+    // routine. Finally, you compute the euclidean 2-norm with the compute_norm2
+    // function.
+    auto res = gko::initialize<real_vec>({0.0}, exec->get_master());
+    A->apply(1.0, lend(x), -1.0, lend(b));
     b->compute_norm2(lend(res));
 
     std::cout << "Residual norm sqrt(r^T r):\n";
diff --git a/include/ginkgo/core/base/lin_op.hpp b/include/ginkgo/core/base/lin_op.hpp
index cde3840c30f..d53df7e8af2 100644
--- a/include/ginkgo/core/base/lin_op.hpp
+++ b/include/ginkgo/core/base/lin_op.hpp
@@ -827,6 +827,66 @@ class EnableLinOp
     GKO_ENABLE_SELF(ConcreteLinOp);
 };
 
+/**
+ * The EnableValueTypeLinOp mixin adds several overloads to the EnableLinOp
+ * mixin. The added overloads allow using scalar arguments for the
+ * LinOp::apply() with four parameters. For more details see @ref EnableLinOp
+ */
+template <typename ConcreteLinOp, typename ValueType,
+          typename PolymorphicBase = LinOp>
+class EnableValueTypedLinOp
+    : public EnableLinOp<ConcreteLinOp, PolymorphicBase> {
+public:
+    using value_type = ValueType;
+
+    using EnableLinOp<ConcreteLinOp, PolymorphicBase>::EnableLinOp;
+    using EnableLinOp<ConcreteLinOp, PolymorphicBase>::apply;
+
+    const ConcreteLinOp *apply(const ValueType alpha, const LinOp *b,
+                               const ValueType beta, LinOp *x) const
+    {
+        auto dense_alpha = scalar_to_dense(alpha, this->get_executor());
+        auto dense_beta = scalar_to_dense(beta, this->get_executor());
+        return this->apply(dense_alpha.get(), b, dense_beta.get(), x);
+    }
+
+    const ConcreteLinOp *apply(const ValueType alpha, const LinOp *b,
+                               const LinOp *beta, LinOp *x) const
+    {
+        auto dense_alpha = scalar_to_dense(alpha, this->get_executor());
+        return this->apply(dense_alpha.get(), b, beta, x);
+    }
+
+    const ConcreteLinOp *apply(const LinOp *alpha, const LinOp *b,
+                               const ValueType beta, LinOp *x) const
+    {
+        auto dense_beta = scalar_to_dense(beta, this->get_executor());
+        return this->apply(alpha, b, dense_beta.get(), x);
+    }
+
+    ConcreteLinOp *apply(const ValueType alpha, const LinOp *b,
+                         const ValueType beta, LinOp *x)
+    {
+        auto dense_alpha = scalar_to_dense(alpha, this->get_executor());
+        auto dense_beta = scalar_to_dense(beta, this->get_executor());
+        return this->apply(dense_alpha.get(), b, dense_beta.get(), x);
+    }
+
+    ConcreteLinOp *apply(const ValueType alpha, const LinOp *b,
+                         const LinOp *beta, LinOp *x)
+    {
+        auto dense_alpha = scalar_to_dense(alpha, this->get_executor());
+        return this->apply(dense_alpha.get(), b, beta, x);
+    }
+
+    ConcreteLinOp *apply(const LinOp *alpha, const LinOp *b,
+                         const ValueType beta, LinOp *x)
+    {
+        auto dense_beta = scalar_to_dense(beta, this->get_executor());
+        return this->apply(alpha, b, dense_beta.get(), x);
+    }
+};
+
 
 /**
  * This is an alias for the EnableDefaultFactory mixin, which correctly sets the
diff --git a/include/ginkgo/core/base/mtx_io.hpp b/include/ginkgo/core/base/mtx_io.hpp
index 3fde2de108d..545a312bf01 100644
--- a/include/ginkgo/core/base/mtx_io.hpp
+++ b/include/ginkgo/core/base/mtx_io.hpp
@@ -120,10 +120,10 @@ inline std::unique_ptr<MatrixType> read(StreamType &&is, MatrixArgs &&... args)
 
 
 /**
- * Reads a matrix stored in matrix market format from an input stream.
+ * Writes a matrix into an output in the matrix market format.
  *
- * @tparam MatrixType  a ReadableFromMatrixData LinOp type used to store the
- *                     matrix once it's been read from disk.
+ * @tparam MatrixType  a WritableToMatrixData LinOp type used to store the
+ *                     matrix.
  * @tparam StreamType  type of stream used to write the data to
  *
  * @param os  output stream where the data is to be written
diff --git a/include/ginkgo/core/base/utils_helper.hpp b/include/ginkgo/core/base/utils_helper.hpp
index 51efb3d5ca7..fa1b906345e 100644
--- a/include/ginkgo/core/base/utils_helper.hpp
+++ b/include/ginkgo/core/base/utils_helper.hpp
@@ -432,6 +432,15 @@ class null_deleter<T[]> {
 };
 
 
+namespace matrix {
+template <typename ValueType>
+class Dense;
+}
+
+template <typename ValueType>
+std::unique_ptr<matrix::Dense<ValueType>> scalar_to_dense(
+    ValueType val, std::shared_ptr<const Executor> exec);
+
 }  // namespace gko
 
 
diff --git a/include/ginkgo/core/matrix/coo.hpp b/include/ginkgo/core/matrix/coo.hpp
index f0872f32787..b627c78ea80 100644
--- a/include/ginkgo/core/matrix/coo.hpp
+++ b/include/ginkgo/core/matrix/coo.hpp
@@ -74,7 +74,7 @@ class CooBuilder;
  * @ingroup LinOp
  */
 template <typename ValueType = default_precision, typename IndexType = int32>
-class Coo : public EnableLinOp<Coo<ValueType, IndexType>>,
+class Coo : public EnableValueTypedLinOp<Coo<ValueType, IndexType>, ValueType>,
             public EnableCreateMethod<Coo<ValueType, IndexType>>,
             public ConvertibleTo<Coo<next_precision<ValueType>, IndexType>>,
             public ConvertibleTo<Csr<ValueType, IndexType>>,
@@ -268,7 +268,7 @@ class Coo : public EnableLinOp<Coo<ValueType, IndexType>>,
      */
     Coo(std::shared_ptr<const Executor> exec, const dim<2> &size = dim<2>{},
         size_type num_nonzeros = {})
-        : EnableLinOp<Coo>(exec, size),
+        : EnableValueTypedLinOp<Coo, ValueType>(exec, size),
           values_(exec, num_nonzeros),
           col_idxs_(exec, num_nonzeros),
           row_idxs_(exec, num_nonzeros)
@@ -298,7 +298,7 @@ class Coo : public EnableLinOp<Coo<ValueType, IndexType>>,
               typename RowIdxsArray>
     Coo(std::shared_ptr<const Executor> exec, const dim<2> &size,
         ValuesArray &&values, ColIdxsArray &&col_idxs, RowIdxsArray &&row_idxs)
-        : EnableLinOp<Coo>(exec, size),
+        : EnableValueTypedLinOp<Coo, ValueType>(exec, size),
           values_{exec, std::forward<ValuesArray>(values)},
           col_idxs_{exec, std::forward<ColIdxsArray>(col_idxs)},
           row_idxs_{exec, std::forward<RowIdxsArray>(row_idxs)}
diff --git a/include/ginkgo/core/matrix/csr.hpp b/include/ginkgo/core/matrix/csr.hpp
index eded7da0434..b9a58b741c9 100644
--- a/include/ginkgo/core/matrix/csr.hpp
+++ b/include/ginkgo/core/matrix/csr.hpp
@@ -117,7 +117,7 @@ void strategy_rebuild_helper(Csr<ValueType, IndexType> *result);
  * @ingroup LinOp
  */
 template <typename ValueType = default_precision, typename IndexType = int32>
-class Csr : public EnableLinOp<Csr<ValueType, IndexType>>,
+class Csr : public EnableValueTypedLinOp<Csr<ValueType, IndexType>, ValueType>,
             public EnableCreateMethod<Csr<ValueType, IndexType>>,
             public ConvertibleTo<Csr<next_precision<ValueType>, IndexType>>,
             public ConvertibleTo<Dense<ValueType>>,
@@ -871,7 +871,7 @@ class Csr : public EnableLinOp<Csr<ValueType, IndexType>>,
     Csr(std::shared_ptr<const Executor> exec, const dim<2> &size = dim<2>{},
         size_type num_nonzeros = {},
         std::shared_ptr<strategy_type> strategy = std::make_shared<sparselib>())
-        : EnableLinOp<Csr>(exec, size),
+        : EnableValueTypedLinOp<Csr, ValueType>(exec, size),
           values_(exec, num_nonzeros),
           col_idxs_(exec, num_nonzeros),
           row_ptrs_(exec, size[0] + 1),
@@ -904,7 +904,7 @@ class Csr : public EnableLinOp<Csr<ValueType, IndexType>>,
     Csr(std::shared_ptr<const Executor> exec, const dim<2> &size,
         ValuesArray &&values, ColIdxsArray &&col_idxs, RowPtrsArray &&row_ptrs,
         std::shared_ptr<strategy_type> strategy = std::make_shared<sparselib>())
-        : EnableLinOp<Csr>(exec, size),
+        : EnableValueTypedLinOp<Csr, ValueType>(exec, size),
           values_{exec, std::forward<ValuesArray>(values)},
           col_idxs_{exec, std::forward<ColIdxsArray>(col_idxs)},
           row_ptrs_{exec, std::forward<RowPtrsArray>(row_ptrs)},
diff --git a/include/ginkgo/core/matrix/dense.hpp b/include/ginkgo/core/matrix/dense.hpp
index 95ca95e720b..bf28d51e9f7 100644
--- a/include/ginkgo/core/matrix/dense.hpp
+++ b/include/ginkgo/core/matrix/dense.hpp
@@ -90,7 +90,7 @@ class SparsityCsr;
  */
 template <typename ValueType = default_precision>
 class Dense
-    : public EnableLinOp<Dense<ValueType>>,
+    : public EnableValueTypedLinOp<Dense<ValueType>, ValueType>,
       public EnableCreateMethod<Dense<ValueType>>,
       public ConvertibleTo<Dense<next_precision<ValueType>>>,
       public ConvertibleTo<Coo<ValueType, int32>>,
@@ -632,6 +632,15 @@ class Dense
         this->scale_impl(make_temporary_clone(exec, alpha).get());
     }
 
+    /**
+     * Overload for scalar alpha, see @ref scale
+     */
+    void scale(const ValueType alpha)
+    {
+        auto dense_alpha = scalar_to_dense(alpha, this->get_executor());
+        this->scale(dense_alpha.get());
+    }
+
     /**
      * Adds `b` scaled by `alpha` to the matrix (aka: BLAS axpy).
      *
@@ -649,6 +658,15 @@ class Dense
                               make_temporary_clone(exec, b).get());
     }
 
+    /**
+     * Overload for scalar alpha, see @ref add_scaled
+     */
+    void add_scaled(const ValueType alpha, const LinOp *b)
+    {
+        auto dense_alpha = scalar_to_dense(alpha, this->get_executor());
+        this->add_scaled(dense_alpha.get(), b);
+    }
+
     /**
      * Computes the column-wise dot product of this matrix and `b`.
      *
@@ -793,7 +811,7 @@ class Dense
      */
     Dense(std::shared_ptr<const Executor> exec, const dim<2> &size,
           size_type stride)
-        : EnableLinOp<Dense>(exec, size),
+        : EnableValueTypedLinOp<Dense, ValueType>(exec, size),
           values_(exec, size[0] * stride),
           stride_(stride)
     {}
@@ -817,7 +835,7 @@ class Dense
     template <typename ValuesArray>
     Dense(std::shared_ptr<const Executor> exec, const dim<2> &size,
           ValuesArray &&values, size_type stride)
-        : EnableLinOp<Dense>(exec, size),
+        : EnableValueTypedLinOp<Dense, ValueType>(exec, size),
           values_{exec, std::forward<ValuesArray>(values)},
           stride_{stride}
     {
diff --git a/include/ginkgo/core/matrix/diagonal.hpp b/include/ginkgo/core/matrix/diagonal.hpp
index ddd4aa73bf0..f7b89a94066 100644
--- a/include/ginkgo/core/matrix/diagonal.hpp
+++ b/include/ginkgo/core/matrix/diagonal.hpp
@@ -66,7 +66,7 @@ class Dense;
  */
 template <typename ValueType = default_precision>
 class Diagonal
-    : public EnableLinOp<Diagonal<ValueType>>,
+    : public EnableValueTypedLinOp<Diagonal<ValueType>, ValueType>,
       public EnableCreateMethod<Diagonal<ValueType>>,
       public ConvertibleTo<Csr<ValueType, int32>>,
       public ConvertibleTo<Csr<ValueType, int64>>,
@@ -169,7 +169,8 @@ class Diagonal
      * @param size  size of the matrix
      */
     Diagonal(std::shared_ptr<const Executor> exec, size_type size)
-        : EnableLinOp<Diagonal>(exec, dim<2>{size}), values_(exec, size)
+        : EnableValueTypedLinOp<Diagonal, ValueType>(exec, dim<2>{size}),
+          values_(exec, size)
     {}
 
     /**
@@ -189,7 +190,7 @@ class Diagonal
     template <typename ValuesArray>
     Diagonal(std::shared_ptr<const Executor> exec, const size_type size,
              ValuesArray &&values)
-        : EnableLinOp<Diagonal>(exec, dim<2>(size)),
+        : EnableValueTypedLinOp<Diagonal, ValueType>(exec, dim<2>(size)),
           values_{exec, std::forward<ValuesArray>(values)}
     {
         GKO_ENSURE_IN_BOUNDS(size - 1, values_.get_num_elems());
diff --git a/include/ginkgo/core/matrix/ell.hpp b/include/ginkgo/core/matrix/ell.hpp
index ea24770fbea..6e0e88fd5d6 100644
--- a/include/ginkgo/core/matrix/ell.hpp
+++ b/include/ginkgo/core/matrix/ell.hpp
@@ -68,7 +68,7 @@ class Csr;
  * @ingroup LinOp
  */
 template <typename ValueType = default_precision, typename IndexType = int32>
-class Ell : public EnableLinOp<Ell<ValueType, IndexType>>,
+class Ell : public EnableValueTypedLinOp<Ell<ValueType, IndexType>, ValueType>,
             public EnableCreateMethod<Ell<ValueType, IndexType>>,
             public ConvertibleTo<Ell<next_precision<ValueType>, IndexType>>,
             public ConvertibleTo<Dense<ValueType>>,
@@ -269,7 +269,7 @@ class Ell : public EnableLinOp<Ell<ValueType, IndexType>>,
      */
     Ell(std::shared_ptr<const Executor> exec, const dim<2> &size,
         size_type num_stored_elements_per_row, size_type stride)
-        : EnableLinOp<Ell>(exec, size),
+        : EnableValueTypedLinOp<Ell, ValueType>(exec, size),
           values_(exec, stride * num_stored_elements_per_row),
           col_idxs_(exec, stride * num_stored_elements_per_row),
           num_stored_elements_per_row_(num_stored_elements_per_row),
@@ -301,7 +301,7 @@ class Ell : public EnableLinOp<Ell<ValueType, IndexType>>,
     Ell(std::shared_ptr<const Executor> exec, const dim<2> &size,
         ValuesArray &&values, ColIdxsArray &&col_idxs,
         size_type num_stored_elements_per_row, size_type stride)
-        : EnableLinOp<Ell>(exec, size),
+        : EnableValueTypedLinOp<Ell, ValueType>(exec, size),
           values_{exec, std::forward<ValuesArray>(values)},
           col_idxs_{exec, std::forward<ColIdxsArray>(col_idxs)},
           num_stored_elements_per_row_{num_stored_elements_per_row},
diff --git a/include/ginkgo/core/matrix/fbcsr.hpp b/include/ginkgo/core/matrix/fbcsr.hpp
index 38cac854bc5..60a600eb9dd 100644
--- a/include/ginkgo/core/matrix/fbcsr.hpp
+++ b/include/ginkgo/core/matrix/fbcsr.hpp
@@ -124,18 +124,19 @@ inline IndexType get_num_blocks(const int block_size, const IndexType size)
  * @ingroup LinOp
  */
 template <typename ValueType = default_precision, typename IndexType = int32>
-class Fbcsr : public EnableLinOp<Fbcsr<ValueType, IndexType>>,
-              public EnableCreateMethod<Fbcsr<ValueType, IndexType>>,
-              public ConvertibleTo<Fbcsr<next_precision<ValueType>, IndexType>>,
-              public ConvertibleTo<Dense<ValueType>>,
-              public ConvertibleTo<Csr<ValueType, IndexType>>,
-              public ConvertibleTo<SparsityCsr<ValueType, IndexType>>,
-              public DiagonalExtractable<ValueType>,
-              public ReadableFromMatrixData<ValueType, IndexType>,
-              public WritableToMatrixData<ValueType, IndexType>,
-              public Transposable,
-              public EnableAbsoluteComputation<
-                  remove_complex<Fbcsr<ValueType, IndexType>>> {
+class Fbcsr
+    : public EnableValueTypedLinOp<Fbcsr<ValueType, IndexType>, ValueType>,
+      public EnableCreateMethod<Fbcsr<ValueType, IndexType>>,
+      public ConvertibleTo<Fbcsr<next_precision<ValueType>, IndexType>>,
+      public ConvertibleTo<Dense<ValueType>>,
+      public ConvertibleTo<Csr<ValueType, IndexType>>,
+      public ConvertibleTo<SparsityCsr<ValueType, IndexType>>,
+      public DiagonalExtractable<ValueType>,
+      public ReadableFromMatrixData<ValueType, IndexType>,
+      public WritableToMatrixData<ValueType, IndexType>,
+      public Transposable,
+      public EnableAbsoluteComputation<
+          remove_complex<Fbcsr<ValueType, IndexType>>> {
     friend class EnableCreateMethod<Fbcsr>;
     friend class EnablePolymorphicObject<Fbcsr, LinOp>;
     friend class Dense<ValueType>;
@@ -342,7 +343,7 @@ class Fbcsr : public EnableLinOp<Fbcsr<ValueType, IndexType>>,
      */
     Fbcsr(std::shared_ptr<const Executor> exec, const dim<2> &size,
           size_type num_nonzeros, int block_size)
-        : EnableLinOp<Fbcsr>(exec, size),
+        : EnableValueTypedLinOp<Fbcsr, ValueType>(exec, size),
           bs_{block_size},
           nbcols_{static_cast<index_type>(
               detail::get_num_blocks(block_size, size[1]))},
@@ -378,7 +379,7 @@ class Fbcsr : public EnableLinOp<Fbcsr<ValueType, IndexType>>,
     Fbcsr(std::shared_ptr<const Executor> exec, const dim<2> &size,
           int block_size, ValuesArray &&values, ColIdxsArray &&col_idxs,
           RowPtrsArray &&row_ptrs)
-        : EnableLinOp<Fbcsr>(exec, size),
+        : EnableValueTypedLinOp<Fbcsr, ValueType>(exec, size),
           bs_{block_size},
           nbcols_{static_cast<index_type>(
               detail::get_num_blocks(block_size, size[1]))},
diff --git a/include/ginkgo/core/matrix/hybrid.hpp b/include/ginkgo/core/matrix/hybrid.hpp
index 7511000fade..1130addd850 100644
--- a/include/ginkgo/core/matrix/hybrid.hpp
+++ b/include/ginkgo/core/matrix/hybrid.hpp
@@ -69,7 +69,7 @@ class Csr;
  */
 template <typename ValueType = default_precision, typename IndexType = int32>
 class Hybrid
-    : public EnableLinOp<Hybrid<ValueType, IndexType>>,
+    : public EnableValueTypedLinOp<Hybrid<ValueType, IndexType>, ValueType>,
       public EnableCreateMethod<Hybrid<ValueType, IndexType>>,
       public ConvertibleTo<Hybrid<next_precision<ValueType>, IndexType>>,
       public ConvertibleTo<Dense<ValueType>>,
@@ -726,7 +726,7 @@ class Hybrid
         size_type num_stored_elements_per_row, size_type stride,
         size_type num_nonzeros = {},
         std::shared_ptr<strategy_type> strategy = std::make_shared<automatic>())
-        : EnableLinOp<Hybrid>(exec, size),
+        : EnableValueTypedLinOp<Hybrid, ValueType>(exec, size),
           ell_(std::move(ell_type::create(
               exec, size, num_stored_elements_per_row, stride))),
           coo_(std::move(coo_type::create(exec, size, num_nonzeros))),
diff --git a/include/ginkgo/core/matrix/identity.hpp b/include/ginkgo/core/matrix/identity.hpp
index 5c5315eedce..ef5ca405c5b 100644
--- a/include/ginkgo/core/matrix/identity.hpp
+++ b/include/ginkgo/core/matrix/identity.hpp
@@ -60,7 +60,7 @@ namespace matrix {
  * @ingroup LinOp
  */
 template <typename ValueType = default_precision>
-class Identity : public EnableLinOp<Identity<ValueType>>,
+class Identity : public EnableValueTypedLinOp<Identity<ValueType>, ValueType>,
                  public EnableCreateMethod<Identity<ValueType>>,
                  public Transposable {
     friend class EnablePolymorphicObject<Identity, LinOp>;
@@ -85,7 +85,7 @@ class Identity : public EnableLinOp<Identity<ValueType>>,
      * @param exec  Executor associated to the matrix
      */
     explicit Identity(std::shared_ptr<const Executor> exec)
-        : EnableLinOp<Identity>(exec)
+        : EnableValueTypedLinOp<Identity, ValueType>(exec)
     {}
 
     /**
@@ -94,7 +94,7 @@ class Identity : public EnableLinOp<Identity<ValueType>>,
      * @param size  size of the matrix (must be square)
      */
     Identity(std::shared_ptr<const Executor> exec, dim<2> size)
-        : EnableLinOp<Identity>(exec, size)
+        : EnableValueTypedLinOp<Identity, ValueType>(exec, size)
     {
         GKO_ASSERT_IS_SQUARE_MATRIX(this);
     }
@@ -105,7 +105,7 @@ class Identity : public EnableLinOp<Identity<ValueType>>,
      * @param size  size of the matrix
      */
     Identity(std::shared_ptr<const Executor> exec, size_type size)
-        : EnableLinOp<Identity>(exec, dim<2>{size})
+        : EnableValueTypedLinOp<Identity, ValueType>(exec, dim<2>{size})
     {}
 
     void apply_impl(const LinOp *b, LinOp *x) const override;
diff --git a/include/ginkgo/core/matrix/sellp.hpp b/include/ginkgo/core/matrix/sellp.hpp
index d2dff8199bb..ac43d959169 100644
--- a/include/ginkgo/core/matrix/sellp.hpp
+++ b/include/ginkgo/core/matrix/sellp.hpp
@@ -65,16 +65,17 @@ class Csr;
  * @ingroup LinOp
  */
 template <typename ValueType = default_precision, typename IndexType = int32>
-class Sellp : public EnableLinOp<Sellp<ValueType, IndexType>>,
-              public EnableCreateMethod<Sellp<ValueType, IndexType>>,
-              public ConvertibleTo<Sellp<next_precision<ValueType>, IndexType>>,
-              public ConvertibleTo<Dense<ValueType>>,
-              public ConvertibleTo<Csr<ValueType, IndexType>>,
-              public DiagonalExtractable<ValueType>,
-              public ReadableFromMatrixData<ValueType, IndexType>,
-              public WritableToMatrixData<ValueType, IndexType>,
-              public EnableAbsoluteComputation<
-                  remove_complex<Sellp<ValueType, IndexType>>> {
+class Sellp
+    : public EnableValueTypedLinOp<Sellp<ValueType, IndexType>, ValueType>,
+      public EnableCreateMethod<Sellp<ValueType, IndexType>>,
+      public ConvertibleTo<Sellp<next_precision<ValueType>, IndexType>>,
+      public ConvertibleTo<Dense<ValueType>>,
+      public ConvertibleTo<Csr<ValueType, IndexType>>,
+      public DiagonalExtractable<ValueType>,
+      public ReadableFromMatrixData<ValueType, IndexType>,
+      public WritableToMatrixData<ValueType, IndexType>,
+      public EnableAbsoluteComputation<
+          remove_complex<Sellp<ValueType, IndexType>>> {
     friend class EnableCreateMethod<Sellp>;
     friend class EnablePolymorphicObject<Sellp, LinOp>;
     friend class Dense<ValueType>;
@@ -322,7 +323,7 @@ class Sellp : public EnableLinOp<Sellp<ValueType, IndexType>>,
      */
     Sellp(std::shared_ptr<const Executor> exec, const dim<2> &size,
           size_type slice_size, size_type stride_factor, size_type total_cols)
-        : EnableLinOp<Sellp>(exec, size),
+        : EnableValueTypedLinOp<Sellp, ValueType>(exec, size),
           values_(exec, slice_size * total_cols),
           col_idxs_(exec, slice_size * total_cols),
           slice_lengths_(exec, ceildiv(size[0], slice_size)),
diff --git a/include/ginkgo/core/multigrid/amgx_pgm.hpp b/include/ginkgo/core/multigrid/amgx_pgm.hpp
index 5c5b36ebf7c..18d7cb752b5 100644
--- a/include/ginkgo/core/multigrid/amgx_pgm.hpp
+++ b/include/ginkgo/core/multigrid/amgx_pgm.hpp
@@ -72,8 +72,9 @@ namespace multigrid {
  * @ingroup LinOp
  */
 template <typename ValueType = default_precision, typename IndexType = int32>
-class AmgxPgm : public EnableLinOp<AmgxPgm<ValueType, IndexType>>,
-                public EnableMultigridLevel<ValueType> {
+class AmgxPgm
+    : public EnableValueTypedLinOp<AmgxPgm<ValueType, IndexType>, ValueType>,
+      public EnableMultigridLevel<ValueType> {
     friend class EnableLinOp<AmgxPgm>;
     friend class EnablePolymorphicObject<AmgxPgm, LinOp>;
 
@@ -156,13 +157,13 @@ class AmgxPgm : public EnableLinOp<AmgxPgm<ValueType, IndexType>>,
     }
 
     explicit AmgxPgm(std::shared_ptr<const Executor> exec)
-        : EnableLinOp<AmgxPgm>(std::move(exec))
+        : EnableValueTypedLinOp<AmgxPgm, ValueType>(std::move(exec))
     {}
 
     explicit AmgxPgm(const Factory *factory,
                      std::shared_ptr<const LinOp> system_matrix)
-        : EnableLinOp<AmgxPgm>(factory->get_executor(),
-                               system_matrix->get_size()),
+        : EnableValueTypedLinOp<AmgxPgm, ValueType>(factory->get_executor(),
+                                                    system_matrix->get_size()),
           EnableMultigridLevel<ValueType>(system_matrix),
           parameters_{factory->get_parameters()},
           system_matrix_{system_matrix},
diff --git a/include/ginkgo/core/preconditioner/ic.hpp b/include/ginkgo/core/preconditioner/ic.hpp
index dcc35a9b317..9494ddd8ec5 100644
--- a/include/ginkgo/core/preconditioner/ic.hpp
+++ b/include/ginkgo/core/preconditioner/ic.hpp
@@ -105,7 +105,9 @@ namespace preconditioner {
  * @ingroup LinOp
  */
 template <typename LSolverType = solver::LowerTrs<>, typename IndexType = int32>
-class Ic : public EnableLinOp<Ic<LSolverType, IndexType>>, public Transposable {
+class Ic : public EnableValueTypedLinOp<Ic<LSolverType, IndexType>,
+                                        typename LSolverType::value_type>,
+           public Transposable {
     friend class EnableLinOp<Ic>;
     friend class EnablePolymorphicObject<Ic, LinOp>;
 
@@ -218,11 +220,12 @@ class Ic : public EnableLinOp<Ic<LSolverType, IndexType>>, public Transposable {
     }
 
     explicit Ic(std::shared_ptr<const Executor> exec)
-        : EnableLinOp<Ic>(std::move(exec))
+        : EnableValueTypedLinOp<Ic, value_type>(std::move(exec))
     {}
 
     explicit Ic(const Factory *factory, std::shared_ptr<const LinOp> lin_op)
-        : EnableLinOp<Ic>(factory->get_executor(), lin_op->get_size()),
+        : EnableValueTypedLinOp<Ic, value_type>(factory->get_executor(),
+                                                lin_op->get_size()),
           parameters_{factory->get_parameters()}
     {
         auto comp =
diff --git a/include/ginkgo/core/preconditioner/ilu.hpp b/include/ginkgo/core/preconditioner/ilu.hpp
index 9eae5447af3..969f5b46348 100644
--- a/include/ginkgo/core/preconditioner/ilu.hpp
+++ b/include/ginkgo/core/preconditioner/ilu.hpp
@@ -112,8 +112,9 @@ namespace preconditioner {
 template <typename LSolverType = solver::LowerTrs<>,
           typename USolverType = solver::UpperTrs<>, bool ReverseApply = false,
           typename IndexType = int32>
-class Ilu : public EnableLinOp<
-                Ilu<LSolverType, USolverType, ReverseApply, IndexType>>,
+class Ilu : public EnableValueTypedLinOp<
+                Ilu<LSolverType, USolverType, ReverseApply, IndexType>,
+                typename LSolverType::value_type>,
             public Transposable {
     friend class EnableLinOp<Ilu>;
     friend class EnablePolymorphicObject<Ilu, LinOp>;
@@ -250,11 +251,12 @@ class Ilu : public EnableLinOp<
     }
 
     explicit Ilu(std::shared_ptr<const Executor> exec)
-        : EnableLinOp<Ilu>(std::move(exec))
+        : EnableValueTypedLinOp<Ilu, value_type>(std::move(exec))
     {}
 
     explicit Ilu(const Factory *factory, std::shared_ptr<const LinOp> lin_op)
-        : EnableLinOp<Ilu>(factory->get_executor(), lin_op->get_size()),
+        : EnableValueTypedLinOp<Ilu, value_type>(factory->get_executor(),
+                                                 lin_op->get_size()),
           parameters_{factory->get_parameters()}
     {
         auto comp =
diff --git a/include/ginkgo/core/preconditioner/isai.hpp b/include/ginkgo/core/preconditioner/isai.hpp
index d3ff87d7b10..d54fc663546 100644
--- a/include/ginkgo/core/preconditioner/isai.hpp
+++ b/include/ginkgo/core/preconditioner/isai.hpp
@@ -103,7 +103,8 @@ enum struct isai_type { lower, upper, general, spd };
  * @ingroup LinOp
  */
 template <isai_type IsaiType, typename ValueType, typename IndexType>
-class Isai : public EnableLinOp<Isai<IsaiType, ValueType, IndexType>>,
+class Isai : public EnableValueTypedLinOp<Isai<IsaiType, ValueType, IndexType>,
+                                          ValueType>,
              public Transposable {
     friend class EnableLinOp<Isai>;
     friend class EnablePolymorphicObject<Isai, LinOp>;
@@ -197,7 +198,7 @@ class Isai : public EnableLinOp<Isai<IsaiType, ValueType, IndexType>>,
 
 protected:
     explicit Isai(std::shared_ptr<const Executor> exec)
-        : EnableLinOp<Isai>(std::move(exec))
+        : EnableValueTypedLinOp<Isai, ValueType>(std::move(exec))
     {}
 
     /**
@@ -208,7 +209,8 @@ class Isai : public EnableLinOp<Isai<IsaiType, ValueType, IndexType>>,
      */
     explicit Isai(const Factory *factory,
                   std::shared_ptr<const LinOp> system_matrix)
-        : EnableLinOp<Isai>(factory->get_executor(), system_matrix->get_size()),
+        : EnableValueTypedLinOp<Isai, ValueType>(factory->get_executor(),
+                                                 system_matrix->get_size()),
           parameters_{factory->get_parameters()}
     {
         const auto skip_sorting = parameters_.skip_sorting;
diff --git a/include/ginkgo/core/preconditioner/jacobi.hpp b/include/ginkgo/core/preconditioner/jacobi.hpp
index dd05296abde..af039a5c8ea 100644
--- a/include/ginkgo/core/preconditioner/jacobi.hpp
+++ b/include/ginkgo/core/preconditioner/jacobi.hpp
@@ -206,10 +206,11 @@ struct block_interleaved_storage_scheme {
  * @ingroup LinOp
  */
 template <typename ValueType = default_precision, typename IndexType = int32>
-class Jacobi : public EnableLinOp<Jacobi<ValueType, IndexType>>,
-               public ConvertibleTo<matrix::Dense<ValueType>>,
-               public WritableToMatrixData<ValueType, IndexType>,
-               public Transposable {
+class Jacobi
+    : public EnableValueTypedLinOp<Jacobi<ValueType, IndexType>, ValueType>,
+      public ConvertibleTo<matrix::Dense<ValueType>>,
+      public WritableToMatrixData<ValueType, IndexType>,
+      public Transposable {
     friend class EnableLinOp<Jacobi>;
     friend class EnablePolymorphicObject<Jacobi, LinOp>;
 
@@ -499,7 +500,7 @@ class Jacobi : public EnableLinOp<Jacobi<ValueType, IndexType>>,
      * @param exec  the executor this object is assigned to
      */
     explicit Jacobi(std::shared_ptr<const Executor> exec)
-        : EnableLinOp<Jacobi>(exec),
+        : EnableValueTypedLinOp<Jacobi, ValueType>(exec),
           num_blocks_{},
           blocks_(exec),
           conditioning_(exec)
@@ -517,8 +518,9 @@ class Jacobi : public EnableLinOp<Jacobi<ValueType, IndexType>>,
      */
     explicit Jacobi(const Factory *factory,
                     std::shared_ptr<const LinOp> system_matrix)
-        : EnableLinOp<Jacobi>(factory->get_executor(),
-                              gko::transpose(system_matrix->get_size())),
+        : EnableValueTypedLinOp<Jacobi, ValueType>(
+              factory->get_executor(),
+              gko::transpose(system_matrix->get_size())),
           parameters_{factory->get_parameters()},
           storage_scheme_{this->compute_storage_scheme(
               parameters_.max_block_size, parameters_.max_block_stride)},
diff --git a/include/ginkgo/core/solver/bicg.hpp b/include/ginkgo/core/solver/bicg.hpp
index da57e30e732..25cc4194628 100644
--- a/include/ginkgo/core/solver/bicg.hpp
+++ b/include/ginkgo/core/solver/bicg.hpp
@@ -76,7 +76,7 @@ namespace solver {
  * @ingroup LinOp
  */
 template <typename ValueType = default_precision>
-class Bicg : public EnableLinOp<Bicg<ValueType>>,
+class Bicg : public EnableValueTypedLinOp<Bicg<ValueType>, ValueType>,
              public Preconditionable,
              public Transposable {
     friend class EnableLinOp<Bicg>;
@@ -163,13 +163,14 @@ class Bicg : public EnableLinOp<Bicg<ValueType>>,
                     LinOp *x) const override;
 
     explicit Bicg(std::shared_ptr<const Executor> exec)
-        : EnableLinOp<Bicg>(std::move(exec))
+        : EnableValueTypedLinOp<Bicg, ValueType>(std::move(exec))
     {}
 
     explicit Bicg(const Factory *factory,
                   std::shared_ptr<const LinOp> system_matrix)
-        : EnableLinOp<Bicg>(factory->get_executor(),
-                            gko::transpose(system_matrix->get_size())),
+        : EnableValueTypedLinOp<Bicg, ValueType>(
+              factory->get_executor(),
+              gko::transpose(system_matrix->get_size())),
           parameters_{factory->get_parameters()},
           system_matrix_{std::move(system_matrix)}
     {
diff --git a/include/ginkgo/core/solver/bicgstab.hpp b/include/ginkgo/core/solver/bicgstab.hpp
index d842a4fbd31..30020e5f163 100644
--- a/include/ginkgo/core/solver/bicgstab.hpp
+++ b/include/ginkgo/core/solver/bicgstab.hpp
@@ -74,7 +74,7 @@ namespace solver {
  * @ingroup LinOp
  */
 template <typename ValueType = default_precision>
-class Bicgstab : public EnableLinOp<Bicgstab<ValueType>>,
+class Bicgstab : public EnableValueTypedLinOp<Bicgstab<ValueType>, ValueType>,
                  public Preconditionable,
                  public Transposable {
     friend class EnableLinOp<Bicgstab>;
@@ -161,13 +161,14 @@ class Bicgstab : public EnableLinOp<Bicgstab<ValueType>>,
                     LinOp *x) const override;
 
     explicit Bicgstab(std::shared_ptr<const Executor> exec)
-        : EnableLinOp<Bicgstab>(std::move(exec))
+        : EnableValueTypedLinOp<Bicgstab, ValueType>(std::move(exec))
     {}
 
     explicit Bicgstab(const Factory *factory,
                       std::shared_ptr<const LinOp> system_matrix)
-        : EnableLinOp<Bicgstab>(factory->get_executor(),
-                                gko::transpose(system_matrix->get_size())),
+        : EnableValueTypedLinOp<Bicgstab, ValueType>(
+              factory->get_executor(),
+              gko::transpose(system_matrix->get_size())),
           parameters_{factory->get_parameters()},
           system_matrix_{std::move(system_matrix)}
     {
diff --git a/include/ginkgo/core/solver/cb_gmres.hpp b/include/ginkgo/core/solver/cb_gmres.hpp
index 4ee9fe18711..25d236c1324 100644
--- a/include/ginkgo/core/solver/cb_gmres.hpp
+++ b/include/ginkgo/core/solver/cb_gmres.hpp
@@ -119,7 +119,7 @@ enum class storage_precision {
  * @ingroup LinOp
  */
 template <typename ValueType = default_precision>
-class CbGmres : public EnableLinOp<CbGmres<ValueType>>,
+class CbGmres : public EnableValueTypedLinOp<CbGmres<ValueType>, ValueType>,
                 public Preconditionable {
     friend class EnableLinOp<CbGmres>;
     friend class EnablePolymorphicObject<CbGmres, LinOp>;
@@ -206,13 +206,13 @@ class CbGmres : public EnableLinOp<CbGmres<ValueType>>,
                     LinOp *x) const override;
 
     explicit CbGmres(std::shared_ptr<const Executor> exec)
-        : EnableLinOp<CbGmres>(std::move(exec))
+        : EnableValueTypedLinOp<CbGmres, ValueType>(std::move(exec))
     {}
 
     explicit CbGmres(const Factory *factory,
                      std::shared_ptr<const LinOp> system_matrix)
-        : EnableLinOp<CbGmres>(factory->get_executor(),
-                               transpose(system_matrix->get_size())),
+        : EnableValueTypedLinOp<CbGmres, ValueType>(
+              factory->get_executor(), transpose(system_matrix->get_size())),
           parameters_{factory->get_parameters()},
           system_matrix_{std::move(system_matrix)}
     {
diff --git a/include/ginkgo/core/solver/cg.hpp b/include/ginkgo/core/solver/cg.hpp
index d84df60906d..4f093e56356 100644
--- a/include/ginkgo/core/solver/cg.hpp
+++ b/include/ginkgo/core/solver/cg.hpp
@@ -70,7 +70,7 @@ namespace solver {
  * @ingroup LinOp
  */
 template <typename ValueType = default_precision>
-class Cg : public EnableLinOp<Cg<ValueType>>,
+class Cg : public EnableValueTypedLinOp<Cg<ValueType>, ValueType>,
            public Preconditionable,
            public Transposable {
     friend class EnableLinOp<Cg>;
@@ -157,13 +157,14 @@ class Cg : public EnableLinOp<Cg<ValueType>>,
                     LinOp *x) const override;
 
     explicit Cg(std::shared_ptr<const Executor> exec)
-        : EnableLinOp<Cg>(std::move(exec))
+        : EnableValueTypedLinOp<Cg, ValueType>(std::move(exec))
     {}
 
     explicit Cg(const Factory *factory,
                 std::shared_ptr<const LinOp> system_matrix)
-        : EnableLinOp<Cg>(factory->get_executor(),
-                          gko::transpose(system_matrix->get_size())),
+        : EnableValueTypedLinOp<Cg, ValueType>(
+              factory->get_executor(),
+              gko::transpose(system_matrix->get_size())),
           parameters_{factory->get_parameters()},
           system_matrix_{std::move(system_matrix)}
     {
diff --git a/include/ginkgo/core/solver/cgs.hpp b/include/ginkgo/core/solver/cgs.hpp
index 14282bb6c8f..9aeef372160 100644
--- a/include/ginkgo/core/solver/cgs.hpp
+++ b/include/ginkgo/core/solver/cgs.hpp
@@ -67,7 +67,7 @@ namespace solver {
  * @ingroup LinOp
  */
 template <typename ValueType = default_precision>
-class Cgs : public EnableLinOp<Cgs<ValueType>>,
+class Cgs : public EnableValueTypedLinOp<Cgs<ValueType>, ValueType>,
             public Preconditionable,
             public Transposable {
     friend class EnableLinOp<Cgs>;
@@ -154,13 +154,14 @@ class Cgs : public EnableLinOp<Cgs<ValueType>>,
                     LinOp *x) const override;
 
     explicit Cgs(std::shared_ptr<const Executor> exec)
-        : EnableLinOp<Cgs>(std::move(exec))
+        : EnableValueTypedLinOp<Cgs, ValueType>(std::move(exec))
     {}
 
     explicit Cgs(const Factory *factory,
                  std::shared_ptr<const LinOp> system_matrix)
-        : EnableLinOp<Cgs>(factory->get_executor(),
-                           gko::transpose(system_matrix->get_size())),
+        : EnableValueTypedLinOp<Cgs, ValueType>(
+              factory->get_executor(),
+              gko::transpose(system_matrix->get_size())),
           parameters_{factory->get_parameters()},
           system_matrix_{std::move(system_matrix)}
     {
diff --git a/include/ginkgo/core/solver/fcg.hpp b/include/ginkgo/core/solver/fcg.hpp
index 55ded85f667..7c88ca99121 100644
--- a/include/ginkgo/core/solver/fcg.hpp
+++ b/include/ginkgo/core/solver/fcg.hpp
@@ -75,7 +75,7 @@ namespace solver {
  * @ingroup LinOp
  */
 template <typename ValueType = default_precision>
-class Fcg : public EnableLinOp<Fcg<ValueType>>,
+class Fcg : public EnableValueTypedLinOp<Fcg<ValueType>, ValueType>,
             public Preconditionable,
             public Transposable {
     friend class EnableLinOp<Fcg>;
@@ -162,13 +162,14 @@ class Fcg : public EnableLinOp<Fcg<ValueType>>,
                     LinOp *x) const override;
 
     explicit Fcg(std::shared_ptr<const Executor> exec)
-        : EnableLinOp<Fcg>(std::move(exec))
+        : EnableValueTypedLinOp<Fcg, ValueType>(std::move(exec))
     {}
 
     explicit Fcg(const Factory *factory,
                  std::shared_ptr<const LinOp> system_matrix)
-        : EnableLinOp<Fcg>(factory->get_executor(),
-                           gko::transpose(system_matrix->get_size())),
+        : EnableValueTypedLinOp<Fcg, ValueType>(
+              factory->get_executor(),
+              gko::transpose(system_matrix->get_size())),
           parameters_{factory->get_parameters()},
           system_matrix_{std::move(system_matrix)}
     {
diff --git a/include/ginkgo/core/solver/gmres.hpp b/include/ginkgo/core/solver/gmres.hpp
index a24e42392d3..17a02fa1bcd 100644
--- a/include/ginkgo/core/solver/gmres.hpp
+++ b/include/ginkgo/core/solver/gmres.hpp
@@ -70,7 +70,7 @@ constexpr size_type default_krylov_dim = 100u;
  * @ingroup LinOp
  */
 template <typename ValueType = default_precision>
-class Gmres : public EnableLinOp<Gmres<ValueType>>,
+class Gmres : public EnableValueTypedLinOp<Gmres<ValueType>, ValueType>,
               public Preconditionable,
               public Transposable {
     friend class EnableLinOp<Gmres>;
@@ -176,13 +176,14 @@ class Gmres : public EnableLinOp<Gmres<ValueType>>,
                     LinOp *x) const override;
 
     explicit Gmres(std::shared_ptr<const Executor> exec)
-        : EnableLinOp<Gmres>(std::move(exec))
+        : EnableValueTypedLinOp<Gmres, ValueType>(std::move(exec))
     {}
 
     explicit Gmres(const Factory *factory,
                    std::shared_ptr<const LinOp> system_matrix)
-        : EnableLinOp<Gmres>(factory->get_executor(),
-                             gko::transpose(system_matrix->get_size())),
+        : EnableValueTypedLinOp<Gmres, ValueType>(
+              factory->get_executor(),
+              gko::transpose(system_matrix->get_size())),
           parameters_{factory->get_parameters()},
           system_matrix_{std::move(system_matrix)}
     {
diff --git a/include/ginkgo/core/solver/idr.hpp b/include/ginkgo/core/solver/idr.hpp
index 45551eea415..876730a56ca 100644
--- a/include/ginkgo/core/solver/idr.hpp
+++ b/include/ginkgo/core/solver/idr.hpp
@@ -79,7 +79,7 @@ namespace solver {
  * @ingroup LinOp
  */
 template <typename ValueType = default_precision>
-class Idr : public EnableLinOp<Idr<ValueType>>,
+class Idr : public EnableValueTypedLinOp<Idr<ValueType>, ValueType>,
             public Preconditionable,
             public Transposable {
     friend class EnableLinOp<Idr>;
@@ -255,13 +255,14 @@ class Idr : public EnableLinOp<Idr<ValueType>>,
                  matrix::Dense<SubspaceType> *dense_x) const;
 
     explicit Idr(std::shared_ptr<const Executor> exec)
-        : EnableLinOp<Idr>(std::move(exec))
+        : EnableValueTypedLinOp<Idr, ValueType>(std::move(exec))
     {}
 
     explicit Idr(const Factory *factory,
                  std::shared_ptr<const LinOp> system_matrix)
-        : EnableLinOp<Idr>(factory->get_executor(),
-                           gko::transpose(system_matrix->get_size())),
+        : EnableValueTypedLinOp<Idr, ValueType>(
+              factory->get_executor(),
+              gko::transpose(system_matrix->get_size())),
           parameters_{factory->get_parameters()},
           system_matrix_{std::move(system_matrix)}
     {
diff --git a/include/ginkgo/core/solver/ir.hpp b/include/ginkgo/core/solver/ir.hpp
index 16ca3e1d3fa..88436796dfc 100644
--- a/include/ginkgo/core/solver/ir.hpp
+++ b/include/ginkgo/core/solver/ir.hpp
@@ -103,7 +103,8 @@ namespace solver {
  * @ingroup LinOp
  */
 template <typename ValueType = default_precision>
-class Ir : public EnableLinOp<Ir<ValueType>>, public Transposable {
+class Ir : public EnableValueTypedLinOp<Ir<ValueType>, ValueType>,
+           public Transposable {
     friend class EnableLinOp<Ir>;
     friend class EnablePolymorphicObject<Ir, LinOp>;
 
@@ -212,13 +213,14 @@ class Ir : public EnableLinOp<Ir<ValueType>>, public Transposable {
                     LinOp *x) const override;
 
     explicit Ir(std::shared_ptr<const Executor> exec)
-        : EnableLinOp<Ir>(std::move(exec))
+        : EnableValueTypedLinOp<Ir, ValueType>(std::move(exec))
     {}
 
     explicit Ir(const Factory *factory,
                 std::shared_ptr<const LinOp> system_matrix)
-        : EnableLinOp<Ir>(factory->get_executor(),
-                          gko::transpose(system_matrix->get_size())),
+        : EnableValueTypedLinOp<Ir, ValueType>(
+              factory->get_executor(),
+              gko::transpose(system_matrix->get_size())),
           parameters_{factory->get_parameters()},
           system_matrix_{std::move(system_matrix)}
     {
diff --git a/include/ginkgo/core/solver/lower_trs.hpp b/include/ginkgo/core/solver/lower_trs.hpp
index c3b95702fd0..ac7105f7769 100644
--- a/include/ginkgo/core/solver/lower_trs.hpp
+++ b/include/ginkgo/core/solver/lower_trs.hpp
@@ -80,8 +80,9 @@ class UpperTrs;
  * @ingroup LinOp
  */
 template <typename ValueType = default_precision, typename IndexType = int32>
-class LowerTrs : public EnableLinOp<LowerTrs<ValueType, IndexType>>,
-                 public Transposable {
+class LowerTrs
+    : public EnableValueTypedLinOp<LowerTrs<ValueType, IndexType>, ValueType>,
+      public Transposable {
     friend class EnableLinOp<LowerTrs>;
     friend class EnablePolymorphicObject<LowerTrs, LinOp>;
     friend class UpperTrs<ValueType, IndexType>;
@@ -137,13 +138,14 @@ class LowerTrs : public EnableLinOp<LowerTrs<ValueType, IndexType>>,
     void generate();
 
     explicit LowerTrs(std::shared_ptr<const Executor> exec)
-        : EnableLinOp<LowerTrs>(std::move(exec))
+        : EnableValueTypedLinOp<LowerTrs, ValueType>(std::move(exec))
     {}
 
     explicit LowerTrs(const Factory *factory,
                       std::shared_ptr<const LinOp> system_matrix)
-        : EnableLinOp<LowerTrs>(factory->get_executor(),
-                                gko::transpose(system_matrix->get_size())),
+        : EnableValueTypedLinOp<LowerTrs, ValueType>(
+              factory->get_executor(),
+              gko::transpose(system_matrix->get_size())),
           parameters_{factory->get_parameters()},
           system_matrix_{}
     {
diff --git a/include/ginkgo/core/solver/upper_trs.hpp b/include/ginkgo/core/solver/upper_trs.hpp
index dd82b5df6e4..8dd0165faac 100644
--- a/include/ginkgo/core/solver/upper_trs.hpp
+++ b/include/ginkgo/core/solver/upper_trs.hpp
@@ -80,8 +80,9 @@ class LowerTrs;
  * @ingroup LinOp
  */
 template <typename ValueType = default_precision, typename IndexType = int32>
-class UpperTrs : public EnableLinOp<UpperTrs<ValueType, IndexType>>,
-                 public Transposable {
+class UpperTrs
+    : public EnableValueTypedLinOp<UpperTrs<ValueType, IndexType>, ValueType>,
+      public Transposable {
     friend class EnableLinOp<UpperTrs>;
     friend class EnablePolymorphicObject<UpperTrs, LinOp>;
     friend class LowerTrs<ValueType, IndexType>;
@@ -137,13 +138,14 @@ class UpperTrs : public EnableLinOp<UpperTrs<ValueType, IndexType>>,
     void generate();
 
     explicit UpperTrs(std::shared_ptr<const Executor> exec)
-        : EnableLinOp<UpperTrs>(std::move(exec))
+        : EnableValueTypedLinOp<UpperTrs, ValueType>(std::move(exec))
     {}
 
     explicit UpperTrs(const Factory *factory,
                       std::shared_ptr<const LinOp> system_matrix)
-        : EnableLinOp<UpperTrs>(factory->get_executor(),
-                                gko::transpose(system_matrix->get_size())),
+        : EnableValueTypedLinOp<UpperTrs, ValueType>(
+              factory->get_executor(),
+              gko::transpose(system_matrix->get_size())),
           parameters_{factory->get_parameters()},
           system_matrix_{}
     {