diff --git a/core/matrix/multivector.cpp b/core/matrix/multivector.cpp
index a71c0a57227..97e465dc155 100644
--- a/core/matrix/multivector.cpp
+++ b/core/matrix/multivector.cpp
@@ -14,6 +14,7 @@
 #include <ginkgo/core/base/math.hpp>
 #include <ginkgo/core/base/precision_dispatch.hpp>
 #include <ginkgo/core/base/temporary_clone.hpp>
+#include <ginkgo/core/base/temporary_conversion.hpp>
 #include <ginkgo/core/base/utils.hpp>
 #include <ginkgo/core/matrix/dense.hpp>
 #include <ginkgo/core/matrix/diagonal.hpp>
@@ -85,381 +86,8 @@ GKO_REGISTER_OPERATION(get_imag, multivector::get_imag);
 
 
 template <typename ValueType>
-void MultiVector<ValueType>::apply_impl(const LinOp* b,
-                                        LinOp* x) const GKO_NOT_IMPLEMENTED;
-
-
-template <typename ValueType>
-void MultiVector<ValueType>::apply_impl(const LinOp* alpha, const LinOp* b,
-                                        const LinOp* beta,
-                                        LinOp* x) const GKO_NOT_IMPLEMENTED;
-
-template <typename ValueType>
-void MultiVector<ValueType>::fill(const ValueType value)
-{
-    this->get_executor()->run(
-        multivector::make_fill(this->get_device_view(), value));
-}
-
-
-template <typename ValueType>
-void MultiVector<ValueType>::scale(ptr_param<const LinOp> alpha)
-{
-    auto exec = this->get_executor();
-    this->scale_impl(make_temporary_clone(exec, alpha).get());
-}
-
-
-template <typename ValueType>
-void MultiVector<ValueType>::inv_scale(ptr_param<const LinOp> alpha)
-{
-    auto exec = this->get_executor();
-    this->inv_scale_impl(make_temporary_clone(exec, alpha).get());
-}
-
-
-template <typename ValueType>
-void MultiVector<ValueType>::add_scaled(ptr_param<const LinOp> alpha,
-                                        ptr_param<const LinOp> b)
-{
-    auto exec = this->get_executor();
-    this->add_scaled_impl(make_temporary_clone(exec, alpha).get(),
-                          make_temporary_clone(exec, b).get());
-}
-
-
-template <typename ValueType>
-void MultiVector<ValueType>::sub_scaled(ptr_param<const LinOp> alpha,
-                                        ptr_param<const LinOp> b)
-{
-    auto exec = this->get_executor();
-    this->sub_scaled_impl(make_temporary_clone(exec, alpha).get(),
-                          make_temporary_clone(exec, b).get());
-}
-
-
-template <typename ValueType>
-void MultiVector<ValueType>::compute_dot(ptr_param<const LinOp> b,
-                                         ptr_param<LinOp> result) const
-{
-    auto exec = this->get_executor();
-    this->compute_dot_impl(make_temporary_clone(exec, b).get(),
-                           make_temporary_output_clone(exec, result).get());
-}
-
-
-template <typename ValueType>
-void MultiVector<ValueType>::compute_conj_dot(ptr_param<const LinOp> b,
-                                              ptr_param<LinOp> result) const
-{
-    auto exec = this->get_executor();
-    this->compute_conj_dot_impl(
-        make_temporary_clone(exec, b).get(),
-        make_temporary_output_clone(exec, result).get());
-}
-
-
-template <typename ValueType>
-void MultiVector<ValueType>::compute_norm2(ptr_param<LinOp> result) const
-{
-    auto exec = this->get_executor();
-    this->compute_norm2_impl(make_temporary_output_clone(exec, result).get());
-}
-
-
-template <typename ValueType>
-void MultiVector<ValueType>::compute_norm1(ptr_param<LinOp> result) const
-{
-    auto exec = this->get_executor();
-    this->compute_norm1_impl(make_temporary_output_clone(exec, result).get());
-}
-
-
-template <typename ValueType>
-void MultiVector<ValueType>::compute_squared_norm2(
-    ptr_param<LinOp> result) const
-{
-    auto exec = this->get_executor();
-    this->compute_squared_norm2_impl(
-        make_temporary_output_clone(exec, result).get());
-}
-
-
-template <typename ValueType>
-void MultiVector<ValueType>::inv_scale_impl(const LinOp* alpha)
-{
-    GKO_ASSERT_EQUAL_ROWS(alpha, dim<2>(1, 1));
-    if (alpha->get_size()[1] != 1) {
-        // different alpha for each column
-        GKO_ASSERT_EQUAL_COLS(this, alpha);
-    }
-    auto exec = this->get_executor();
-    // if alpha is real (convertible to real) and ValueType complex
-    if (dynamic_cast<const ConvertibleTo<MultiVector<>>*>(alpha) &&
-        is_complex<ValueType>()) {
-        // use the real-complex kernel
-        exec->run(multivector::make_inv_scale(
-            make_temporary_conversion<remove_complex<ValueType>>(alpha)
-                ->get_const_device_view(),
-            dynamic_cast<complex_type*>(this)->get_device_view()));
-        // this last cast is a no-op for complex value type and the branch is
-        // never taken for real value type
-    } else {
-        // otherwise: use the normal kernel
-        exec->run(multivector::make_inv_scale(
-            make_temporary_conversion<ValueType>(alpha)
-                ->get_const_device_view(),
-            this->get_device_view()));
-    }
-}
-
-
-template <typename ValueType>
-void MultiVector<ValueType>::scale_impl(const LinOp* alpha)
-{
-    GKO_ASSERT_EQUAL_ROWS(alpha, dim<2>(1, 1));
-    if (alpha->get_size()[1] != 1) {
-        // different alpha for each column
-        GKO_ASSERT_EQUAL_COLS(this, alpha);
-    }
-    auto exec = this->get_executor();
-    // if alpha is real (convertible to real) and ValueType complex
-    if (dynamic_cast<const ConvertibleTo<MultiVector<>>*>(alpha) &&
-        is_complex<ValueType>()) {
-        // use the real-complex kernel
-        exec->run(multivector::make_scale(
-            make_temporary_conversion<remove_complex<ValueType>>(alpha)
-                ->get_const_device_view(),
-            dynamic_cast<complex_type*>(this)->get_device_view()));
-        // this last cast is a no-op for complex value type and the branch is
-        // never taken for real value type
-    } else {
-        // otherwise: use the normal kernel
-        exec->run(
-            multivector::make_scale(make_temporary_conversion<ValueType>(alpha)
-                                        ->get_const_device_view(),
-                                    this->get_device_view()));
-    }
-}
-
-
-template <typename ValueType>
-void MultiVector<ValueType>::add_scaled_impl(const LinOp* alpha, const LinOp* b)
-{
-    GKO_ASSERT_EQUAL_ROWS(alpha, dim<2>(1, 1));
-    if (alpha->get_size()[1] != 1) {
-        // different alpha for each column
-        GKO_ASSERT_EQUAL_COLS(this, alpha);
-    }
-    GKO_ASSERT_EQUAL_DIMENSIONS(this, b);
-    auto exec = this->get_executor();
-
-    // if alpha is real and value type complex
-    if (dynamic_cast<const ConvertibleTo<MultiVector<>>*>(alpha) &&
-        is_complex<ValueType>()) {
-        exec->run(multivector::make_add_scaled(
-            make_temporary_conversion<remove_complex<ValueType>>(alpha)
-                ->get_const_device_view(),
-            make_temporary_conversion<to_complex<ValueType>>(b)
-                ->get_const_device_view(),
-            dynamic_cast<complex_type*>(this)->get_device_view()));
-    } else {
-        exec->run(multivector::make_add_scaled(
-            make_temporary_conversion<ValueType>(alpha)
-                ->get_const_device_view(),
-            make_temporary_conversion<ValueType>(b)->get_const_device_view(),
-            this->get_device_view()));
-    }
-}
-
-
-template <typename ValueType>
-void MultiVector<ValueType>::sub_scaled_impl(const LinOp* alpha, const LinOp* b)
-{
-    GKO_ASSERT_EQUAL_ROWS(alpha, dim<2>(1, 1));
-    if (alpha->get_size()[1] != 1) {
-        // different alpha for each column
-        GKO_ASSERT_EQUAL_COLS(this, alpha);
-    }
-    GKO_ASSERT_EQUAL_DIMENSIONS(this, b);
-    auto exec = this->get_executor();
-
-    if (dynamic_cast<const ConvertibleTo<MultiVector<>>*>(alpha) &&
-        is_complex<ValueType>()) {
-        exec->run(multivector::make_sub_scaled(
-            make_temporary_conversion<remove_complex<ValueType>>(alpha)
-                ->get_const_device_view(),
-            make_temporary_conversion<to_complex<ValueType>>(b)
-                ->get_const_device_view(),
-            dynamic_cast<complex_type*>(this)->get_device_view()));
-    } else {
-        exec->run(multivector::make_sub_scaled(
-            make_temporary_conversion<ValueType>(alpha)
-                ->get_const_device_view(),
-            make_temporary_conversion<ValueType>(b)->get_const_device_view(),
-            this->get_device_view()));
-    }
-}
-
-
-template <typename ValueType>
-void MultiVector<ValueType>::compute_dot(ptr_param<const LinOp> b,
-                                         ptr_param<LinOp> result,
-                                         array<char>& tmp) const
-{
-    GKO_ASSERT_EQUAL_DIMENSIONS(this, b);
-    GKO_ASSERT_EQUAL_DIMENSIONS(result, dim<2>(1, this->get_size()[1]));
-    auto exec = this->get_executor();
-    if (tmp.get_executor() != exec) {
-        tmp.clear();
-        tmp.set_executor(exec);
-    }
-    auto local_b = make_temporary_clone(exec, b);
-    auto local_res = make_temporary_clone(exec, result);
-    auto dense_b = make_temporary_conversion<ValueType>(local_b.get());
-    auto dense_res = make_temporary_conversion<ValueType>(local_res.get());
-    exec->run(multivector::make_compute_dot(this->get_const_device_view(),
-                                            dense_b->get_const_device_view(),
-                                            dense_res->get_device_view(), tmp));
-}
-
-
-template <typename ValueType>
-void MultiVector<ValueType>::compute_dot_impl(const LinOp* b,
-                                              LinOp* result) const
-{
-    GKO_ASSERT_EQUAL_DIMENSIONS(this, b);
-    GKO_ASSERT_EQUAL_DIMENSIONS(result, dim<2>(1, this->get_size()[1]));
-    auto exec = this->get_executor();
-    auto dense_b = make_temporary_conversion<ValueType>(b);
-    auto dense_res = make_temporary_conversion<ValueType>(result);
-    array<char> tmp{exec};
-    exec->run(multivector::make_compute_dot(this->get_const_device_view(),
-                                            dense_b->get_const_device_view(),
-                                            dense_res->get_device_view(), tmp));
-}
-
-
-template <typename ValueType>
-void MultiVector<ValueType>::compute_conj_dot(ptr_param<const LinOp> b,
-                                              ptr_param<LinOp> result,
-                                              array<char>& tmp) const
-{
-    GKO_ASSERT_EQUAL_DIMENSIONS(this, b);
-    GKO_ASSERT_EQUAL_DIMENSIONS(result, dim<2>(1, this->get_size()[1]));
-    auto exec = this->get_executor();
-    if (tmp.get_executor() != exec) {
-        tmp.clear();
-        tmp.set_executor(exec);
-    }
-    auto local_b = make_temporary_clone(exec, b);
-    auto local_res = make_temporary_clone(exec, result);
-    auto dense_b = make_temporary_conversion<ValueType>(local_b.get());
-    auto dense_res = make_temporary_conversion<ValueType>(local_res.get());
-    exec->run(multivector::make_compute_conj_dot(
-        this->get_const_device_view(), dense_b->get_const_device_view(),
-        dense_res->get_device_view(), tmp));
-}
-
-
-template <typename ValueType>
-void MultiVector<ValueType>::compute_conj_dot_impl(const LinOp* b,
-                                                   LinOp* result) const
-{
-    GKO_ASSERT_EQUAL_DIMENSIONS(this, b);
-    GKO_ASSERT_EQUAL_DIMENSIONS(result, dim<2>(1, this->get_size()[1]));
-    auto exec = this->get_executor();
-    auto dense_b = make_temporary_conversion<ValueType>(b);
-    auto dense_res = make_temporary_conversion<ValueType>(result);
-    array<char> tmp{exec};
-    exec->run(multivector::make_compute_conj_dot(
-        this->get_const_device_view(), dense_b->get_const_device_view(),
-        dense_res->get_device_view(), tmp));
-}
-
-
-template <typename ValueType>
-void MultiVector<ValueType>::compute_norm2(ptr_param<LinOp> result,
-                                           array<char>& tmp) const
-{
-    GKO_ASSERT_EQUAL_DIMENSIONS(result, dim<2>(1, this->get_size()[1]));
-    auto exec = this->get_executor();
-    if (tmp.get_executor() != exec) {
-        tmp.clear();
-        tmp.set_executor(exec);
-    }
-    auto local_result = make_temporary_clone(exec, result);
-    auto dense_res = make_temporary_conversion<remove_complex<ValueType>>(
-        local_result.get());
-    exec->run(multivector::make_compute_norm2(
-        this->get_const_device_view(), dense_res->get_device_view(), tmp));
-}
-
-
-template <typename ValueType>
-void MultiVector<ValueType>::compute_norm2_impl(LinOp* result) const
-{
-    GKO_ASSERT_EQUAL_DIMENSIONS(result, dim<2>(1, this->get_size()[1]));
-    auto exec = this->get_executor();
-    auto dense_res =
-        make_temporary_conversion<remove_complex<ValueType>>(result);
-    array<char> tmp{exec};
-    exec->run(multivector::make_compute_norm2(
-        this->get_const_device_view(), dense_res->get_device_view(), tmp));
-}
-
-
-template <typename ValueType>
-void MultiVector<ValueType>::compute_norm1(ptr_param<LinOp> result,
-                                           array<char>& tmp) const
-{
-    GKO_ASSERT_EQUAL_DIMENSIONS(result, dim<2>(1, this->get_size()[1]));
-    auto exec = this->get_executor();
-    if (tmp.get_executor() != exec) {
-        tmp.clear();
-        tmp.set_executor(exec);
-    }
-    auto local_result = make_temporary_clone(exec, result);
-    auto dense_res = make_temporary_conversion<remove_complex<ValueType>>(
-        local_result.get());
-    exec->run(multivector::make_compute_norm1(
-        this->get_const_device_view(), dense_res->get_device_view(), tmp));
-}
-
-
-template <typename ValueType>
-void MultiVector<ValueType>::compute_norm1_impl(LinOp* result) const
-{
-    GKO_ASSERT_EQUAL_DIMENSIONS(result, dim<2>(1, this->get_size()[1]));
-    auto exec = this->get_executor();
-    auto dense_res =
-        make_temporary_conversion<remove_complex<ValueType>>(result);
-    array<char> tmp{exec};
-    exec->run(multivector::make_compute_norm1(
-        this->get_const_device_view(), dense_res->get_device_view(), tmp));
-}
-
-
-template <typename ValueType>
-void MultiVector<ValueType>::compute_squared_norm2(ptr_param<LinOp> result,
-                                                   array<char>& tmp) const
-{
-    GKO_ASSERT_EQUAL_DIMENSIONS(result, dim<2>(1, this->get_size()[1]));
-    auto exec = this->get_executor();
-    if (tmp.get_executor() != exec) {
-        tmp.clear();
-        tmp.set_executor(exec);
-    }
-    auto local_result = make_temporary_clone(exec, result);
-    auto dense_res = make_temporary_conversion<remove_complex<ValueType>>(
-        local_result.get());
-    exec->run(multivector::make_compute_squared_norm2(
-        this->get_const_device_view(), dense_res->get_device_view(), tmp));
-}
-
-
-template <typename ValueType>
-void MultiVector<ValueType>::compute_mean(ptr_param<LinOp> result) const
+void MultiVector<ValueType>::compute_mean(
+    ptr_param<AbstractMultiVector> result) const
 {
     auto exec = this->get_executor();
     this->compute_mean_impl(make_temporary_output_clone(exec, result).get());
@@ -467,7 +95,7 @@ void MultiVector<ValueType>::compute_mean(ptr_param<LinOp> result) const
 
 
 template <typename ValueType>
-void MultiVector<ValueType>::compute_mean(ptr_param<LinOp> result,
+void MultiVector<ValueType>::compute_mean(ptr_param<AbstractMultiVector> result,
                                           array<char>& tmp) const
 {
     GKO_ASSERT_EQUAL_COLS(result, this);
@@ -476,24 +104,15 @@ void MultiVector<ValueType>::compute_mean(ptr_param<LinOp> result,
         tmp.clear();
         tmp.set_executor(exec);
     }
-    auto dense_res = make_temporary_conversion<ValueType>(result);
+    auto dense_res = as<MultiVector>(result->as_precision(this));
     exec->run(multivector::make_compute_mean(
         this->get_const_device_view(), dense_res->get_device_view(), tmp));
 }
 
 
 template <typename ValueType>
-void MultiVector<ValueType>::compute_squared_norm2_impl(LinOp* result) const
-{
-    auto exec = this->get_executor();
-    array<char> tmp{exec};
-    this->compute_squared_norm2(make_temporary_output_clone(exec, result).get(),
-                                tmp);
-}
-
-
-template <typename ValueType>
-void MultiVector<ValueType>::compute_mean_impl(LinOp* result) const
+void MultiVector<ValueType>::compute_mean_impl(
+    AbstractMultiVector* result) const
 {
     auto exec = this->get_executor();
     array<char> tmp{exec};
@@ -507,7 +126,7 @@ MultiVector<ValueType>& MultiVector<ValueType>::operator=(
 {
     if (&other != this) {
         auto old_size = this->get_size();
-        LinOp::operator=(other);
+        AbstractMultiVector::operator=(other);
         // NOTE: keep this consistent with resize(...)
         if (old_size != other.get_size()) {
             this->stride_ = this->get_size()[1];
@@ -537,7 +156,7 @@ MultiVector<ValueType>& MultiVector<ValueType>::operator=(
     MultiVector<ValueType>&& other)
 {
     if (&other != this) {
-        LinOp::operator=(std::move(other));
+        AbstractMultiVector::operator=(std::move(other));
         values_ = std::move(other.values_);
         stride_ = std::exchange(other.stride_, 0);
     }
@@ -561,6 +180,16 @@ MultiVector<ValueType>::MultiVector(MultiVector<ValueType>&& other)
 }
 
 
+template <typename ValueType>
+std::unique_ptr<MultiVector<ValueType>>
+MultiVector<ValueType>::create_with_type_of(
+    ptr_param<const MultiVector> other, std::shared_ptr<const Executor> exec,
+    const dim<2>& size, size_type stride)
+{
+    return other->create_with_type_of_impl(exec, size, stride);
+}
+
+
 template <typename ValueType>
 void MultiVector<ValueType>::convert_to(
     MultiVector<next_precision<ValueType>>* result) const
@@ -767,7 +396,8 @@ void MultiVector<ValueType>::write(mat_data32& data) const
 
 
 template <typename ValueType>
-std::unique_ptr<LinOp> MultiVector<ValueType>::transpose() const
+std::unique_ptr<MultiVector<ValueType>> MultiVector<ValueType>::transpose()
+    const
 {
     auto result = MultiVector::create(this->get_executor(),
                                       gko::transpose(this->get_size()));
@@ -777,7 +407,8 @@ std::unique_ptr<LinOp> MultiVector<ValueType>::transpose() const
 
 
 template <typename ValueType>
-std::unique_ptr<LinOp> MultiVector<ValueType>::conj_transpose() const
+std::unique_ptr<MultiVector<ValueType>> MultiVector<ValueType>::conj_transpose()
+    const
 {
     auto result = MultiVector::create(this->get_executor(),
                                       gko::transpose(this->get_size()));
@@ -1021,7 +652,7 @@ void MultiVector<ValueType>::row_gather_impl(
 
 
 template <typename ValueType>
-std::unique_ptr<LinOp> MultiVector<ValueType>::permute(
+std::unique_ptr<MultiVector<ValueType>> MultiVector<ValueType>::permute(
     const array<int32>* permutation_indices) const
 {
     auto result = MultiVector::create(this->get_executor(), this->get_size());
@@ -1031,7 +662,7 @@ std::unique_ptr<LinOp> MultiVector<ValueType>::permute(
 
 
 template <typename ValueType>
-std::unique_ptr<LinOp> MultiVector<ValueType>::permute(
+std::unique_ptr<MultiVector<ValueType>> MultiVector<ValueType>::permute(
     const array<int64>* permutation_indices) const
 {
     auto result = MultiVector::create(this->get_executor(), this->get_size());
@@ -1152,7 +783,7 @@ void MultiVector<ValueType>::permute(
 
 
 template <typename ValueType>
-std::unique_ptr<LinOp> MultiVector<ValueType>::inverse_permute(
+std::unique_ptr<MultiVector<ValueType>> MultiVector<ValueType>::inverse_permute(
     const array<int32>* permutation_indices) const
 {
     auto result = MultiVector::create(this->get_executor(), this->get_size());
@@ -1162,7 +793,7 @@ std::unique_ptr<LinOp> MultiVector<ValueType>::inverse_permute(
 
 
 template <typename ValueType>
-std::unique_ptr<LinOp> MultiVector<ValueType>::inverse_permute(
+std::unique_ptr<MultiVector<ValueType>> MultiVector<ValueType>::inverse_permute(
     const array<int64>* permutation_indices) const
 {
     auto result = MultiVector::create(this->get_executor(), this->get_size());
@@ -1192,7 +823,7 @@ void MultiVector<ValueType>::inverse_permute(
 
 
 template <typename ValueType>
-std::unique_ptr<LinOp> MultiVector<ValueType>::row_permute(
+std::unique_ptr<MultiVector<ValueType>> MultiVector<ValueType>::row_permute(
     const array<int32>* permutation_indices) const
 {
     auto result = MultiVector::create(this->get_executor(), this->get_size());
@@ -1202,7 +833,7 @@ std::unique_ptr<LinOp> MultiVector<ValueType>::row_permute(
 
 
 template <typename ValueType>
-std::unique_ptr<LinOp> MultiVector<ValueType>::row_permute(
+std::unique_ptr<MultiVector<ValueType>> MultiVector<ValueType>::row_permute(
     const array<int64>* permutation_indices) const
 {
     auto result = MultiVector::create(this->get_executor(), this->get_size());
@@ -1258,13 +889,14 @@ namespace {
 
 
 template <typename ValueType, typename Function>
-void gather_mixed_real_complex(Function fn, LinOp* out)
+void gather_mixed_real_complex(Function fn, AbstractMultiVector* out)
 {
 #ifdef GINKGO_MIXED_PRECISION
     run<matrix::MultiVector, ValueType, next_precision<ValueType>,
         next_precision<ValueType, 2>, next_precision<ValueType, 3>>(out, fn);
 #else
-    precision_dispatch<ValueType>(fn, out);
+    fn(as<MultiVector<ValueType>>(out->as_precision(precision_v<ValueType>))
+           .get());
 #endif
 }
 
@@ -1273,8 +905,9 @@ void gather_mixed_real_complex(Function fn, LinOp* out)
 
 
 template <typename ValueType>
-void MultiVector<ValueType>::row_gather(const array<int32>* row_idxs,
-                                        ptr_param<LinOp> row_collection) const
+void MultiVector<ValueType>::row_gather(
+    const array<int32>* row_idxs,
+    ptr_param<AbstractMultiVector> row_collection) const
 {
     gather_mixed_real_complex<ValueType>(
         [&](auto dense) { this->row_gather_impl(row_idxs, dense); },
@@ -1283,8 +916,9 @@ void MultiVector<ValueType>::row_gather(const array<int32>* row_idxs,
 
 
 template <typename ValueType>
-void MultiVector<ValueType>::row_gather(const array<int64>* row_idxs,
-                                        ptr_param<LinOp> row_collection) const
+void MultiVector<ValueType>::row_gather(
+    const array<int64>* row_idxs,
+    ptr_param<AbstractMultiVector> row_collection) const
 {
     gather_mixed_real_complex<ValueType>(
         [&](auto dense) { this->row_gather_impl(row_idxs, dense); },
@@ -1293,10 +927,9 @@ void MultiVector<ValueType>::row_gather(const array<int64>* row_idxs,
 
 
 template <typename ValueType>
-void MultiVector<ValueType>::row_gather(ptr_param<const LinOp> alpha,
-                                        const array<int32>* gather_indices,
-                                        ptr_param<const LinOp> beta,
-                                        ptr_param<LinOp> out) const
+void MultiVector<ValueType>::row_gather(
+    ptr_param<const LinOp> alpha, const array<int32>* gather_indices,
+    ptr_param<const LinOp> beta, ptr_param<AbstractMultiVector> out) const
 {
     auto dense_alpha = make_temporary_conversion<ValueType>(alpha);
     auto dense_beta = make_temporary_conversion<ValueType>(beta);
@@ -1311,10 +944,9 @@ void MultiVector<ValueType>::row_gather(ptr_param<const LinOp> alpha,
 }
 
 template <typename ValueType>
-void MultiVector<ValueType>::row_gather(ptr_param<const LinOp> alpha,
-                                        const array<int64>* gather_indices,
-                                        ptr_param<const LinOp> beta,
-                                        ptr_param<LinOp> out) const
+void MultiVector<ValueType>::row_gather(
+    ptr_param<const LinOp> alpha, const array<int64>* gather_indices,
+    ptr_param<const LinOp> beta, ptr_param<AbstractMultiVector> out) const
 {
     auto dense_alpha = make_temporary_conversion<ValueType>(alpha);
     auto dense_beta = make_temporary_conversion<ValueType>(beta);
@@ -1330,7 +962,7 @@ void MultiVector<ValueType>::row_gather(ptr_param<const LinOp> alpha,
 
 
 template <typename ValueType>
-std::unique_ptr<LinOp> MultiVector<ValueType>::column_permute(
+std::unique_ptr<MultiVector<ValueType>> MultiVector<ValueType>::column_permute(
     const array<int32>* permutation_indices) const
 {
     auto result = MultiVector::create(this->get_executor(), this->get_size());
@@ -1340,7 +972,7 @@ std::unique_ptr<LinOp> MultiVector<ValueType>::column_permute(
 
 
 template <typename ValueType>
-std::unique_ptr<LinOp> MultiVector<ValueType>::column_permute(
+std::unique_ptr<MultiVector<ValueType>> MultiVector<ValueType>::column_permute(
     const array<int64>* permutation_indices) const
 {
     auto result = MultiVector::create(this->get_executor(), this->get_size());
@@ -1370,7 +1002,8 @@ void MultiVector<ValueType>::column_permute(
 
 
 template <typename ValueType>
-std::unique_ptr<LinOp> MultiVector<ValueType>::inverse_row_permute(
+std::unique_ptr<MultiVector<ValueType>>
+MultiVector<ValueType>::inverse_row_permute(
     const array<int32>* permutation_indices) const
 {
     auto result = MultiVector::create(this->get_executor(), this->get_size());
@@ -1380,7 +1013,8 @@ std::unique_ptr<LinOp> MultiVector<ValueType>::inverse_row_permute(
 
 
 template <typename ValueType>
-std::unique_ptr<LinOp> MultiVector<ValueType>::inverse_row_permute(
+std::unique_ptr<MultiVector<ValueType>>
+MultiVector<ValueType>::inverse_row_permute(
     const array<int64>* permutation_indices) const
 {
     auto result = MultiVector::create(this->get_executor(), this->get_size());
@@ -1410,7 +1044,8 @@ void MultiVector<ValueType>::inverse_row_permute(
 
 
 template <typename ValueType>
-std::unique_ptr<LinOp> MultiVector<ValueType>::inverse_column_permute(
+std::unique_ptr<MultiVector<ValueType>>
+MultiVector<ValueType>::inverse_column_permute(
     const array<int32>* permutation_indices) const
 {
     auto result = MultiVector::create(this->get_executor(), this->get_size());
@@ -1420,7 +1055,8 @@ std::unique_ptr<LinOp> MultiVector<ValueType>::inverse_column_permute(
 
 
 template <typename ValueType>
-std::unique_ptr<LinOp> MultiVector<ValueType>::inverse_column_permute(
+std::unique_ptr<MultiVector<ValueType>>
+MultiVector<ValueType>::inverse_column_permute(
     const array<int64>* permutation_indices) const
 {
     auto result = MultiVector::create(this->get_executor(), this->get_size());
@@ -1536,16 +1172,212 @@ void MultiVector<ValueType>::scale_permute(
 
 
 template <typename ValueType>
-void MultiVector<ValueType>::compute_absolute_inplace()
+typename MultiVector<ValueType>::device_view
+MultiVector<ValueType>::get_device_view()
+{
+    return device_view{this->get_size(), this->get_stride(),
+                       this->get_values()};
+}
+
+
+template <typename ValueType>
+typename MultiVector<ValueType>::const_device_view
+MultiVector<ValueType>::get_const_device_view() const
+{
+    return const_device_view{this->get_size(), this->get_stride(),
+                             this->get_const_values()};
+}
+
+
+template <typename ValueType>
+template <typename OtherValueType>
+temporary_conversion<MultiVector<OtherValueType>>
+MultiVector<ValueType>::as_precision()
+{
+    if constexpr (is_complex<ValueType>() == is_complex<OtherValueType>()) {
+        // The value types are either both real or both complex
+        return temporary_conversion<MultiVector<OtherValueType>>::create(this);
+    } else if constexpr (is_complex<ValueType>() &&
+                         std::is_same_v<to_complex<OtherValueType>,
+                                        ValueType>) {
+        // The value type of this is complex and the other value type is the
+        // corresponding real value type (std::complex<double> vs double).
+        // This conversion is allowed, since the real view of this can be used
+        return temporary_conversion<MultiVector<OtherValueType>>::create(
+            this->create_real_view());
+    } else {
+        // Conversions from real to complex (or vice versa) that change the
+        // precision are not allowed.
+        GKO_NOT_IMPLEMENTED;
+    }
+}
+
+#define GKO_DECLARE_MULTIVECTOR_AS_PRECISION(ValueType, OtherValueType) \
+    auto MultiVector<ValueType>::as_precision()                         \
+        ->temporary_conversion<MultiVector<OtherValueType>>
+#define GKO_DECLARE_MULTIVECTOR_AS_PRECISION_same(ValueType) \
+    GKO_DECLARE_MULTIVECTOR_AS_PRECISION(ValueType, ValueType)
+GKO_INSTANTIATE_FOR_EACH_VALUE_CONVERSION(GKO_DECLARE_MULTIVECTOR_AS_PRECISION);
+GKO_INSTANTIATE_FOR_EACH_VALUE_TYPE(GKO_DECLARE_MULTIVECTOR_AS_PRECISION_same);
+
+
+template <typename ValueType>
+template <typename OtherValueType>
+temporary_conversion<const MultiVector<OtherValueType>>
+MultiVector<ValueType>::as_precision() const
+{
+    if constexpr (is_complex<ValueType>() == is_complex<OtherValueType>()) {
+        // The value types are either both real or both complex
+        return temporary_conversion<const MultiVector<OtherValueType>>::create(
+            this);
+    } else if constexpr (is_complex<ValueType>() &&
+                         std::is_same_v<to_complex<OtherValueType>,
+                                        ValueType>) {
+        // The value type of this is complex and the other value type is the
+        // corresponding real value type (std::complex<double> vs double).
+        // This conversion is allowed, since the real view of this can be used
+        return gko::detail::
+            temporary_conversion<const MultiVector<OtherValueType>>::create(
+                this->create_real_view());
+    } else {
+        // Conversions from real to complex (or vice versa) that change the
+        // precision are not allowed.
+        GKO_NOT_IMPLEMENTED;
+    }
+}
+
+#define GKO_DECLARE_MULTIVECTOR_CONST_AS_PRECISION(ValueType, OtherValueType) \
+    auto MultiVector<ValueType>::as_precision()                               \
+        const->temporary_conversion<const MultiVector<OtherValueType>>
+#define GKO_DECLARE_MULTIVECTOR_CONST_AS_PRECISION_same(ValueType) \
+    GKO_DECLARE_MULTIVECTOR_CONST_AS_PRECISION(ValueType, ValueType)
+GKO_INSTANTIATE_FOR_EACH_VALUE_CONVERSION(
+    GKO_DECLARE_MULTIVECTOR_CONST_AS_PRECISION);
+GKO_INSTANTIATE_FOR_EACH_VALUE_TYPE(
+    GKO_DECLARE_MULTIVECTOR_CONST_AS_PRECISION_same);
+
+
+template <typename ValueType>
+void MultiVector<ValueType>::compute_absolute_inplace_impl()
 {
     this->get_executor()->run(
         multivector::make_inplace_absolute_dense(this->get_device_view()));
 }
 
 
+template <typename ValueType>
+std::unique_ptr<MultiVector<ValueType>>
+MultiVector<ValueType>::create_with_same_config_impl() const
+{
+    return MultiVector::create(this->get_executor(), this->get_size(),
+                               this->get_stride());
+}
+
+
+template <typename ValueType>
+std::unique_ptr<MultiVector<ValueType>>
+MultiVector<ValueType>::create_with_type_of_impl(
+    std::shared_ptr<const Executor> exec, const dim<2>& global_size,
+    const dim<2>& local_size, size_type stride) const
+{
+    GKO_ASSERT_EQUAL_DIMENSIONS(global_size, local_size);
+    return create_with_type_of_impl(std::move(exec), global_size, stride);
+}
+
+
+template <typename ValueType>
+std::unique_ptr<MultiVector<ValueType>>
+MultiVector<ValueType>::create_subview_impl(local_span rows, local_span columns)
+{
+    return create_subview_impl(rows, columns,
+                               dim<2>(rows.length(), columns.length()));
+}
+
+
+template <typename ValueType>
+std::unique_ptr<const MultiVector<ValueType>>
+MultiVector<ValueType>::create_subview_impl(local_span rows,
+                                            local_span columns) const
+{
+    return const_cast<MultiVector&>(*this).create_subview(rows, columns);
+}
+
+
+template <typename ValueType>
+std::unique_ptr<MultiVector<ValueType>>
+MultiVector<ValueType>::create_subview_impl(local_span rows, local_span columns,
+                                            dim<2> global_size)
+{
+    dim<2> actual_size{rows.length(), columns.length()};
+    GKO_ASSERT_EQUAL_DIMENSIONS(actual_size, global_size);
+
+    row_major_range range_this{this->get_values(), this->get_size()[0],
+                               this->get_size()[1], this->get_stride()};
+    auto sub_range = range_this(rows, columns);
+    size_type storage_size =
+        rows.length() > 0 ? sub_range.length(1) +
+                                (sub_range.length(0) - 1) * this->get_stride()
+                          : 0;
+    return MultiVector::create(
+        this->get_executor(), dim<2>{sub_range.length(0), sub_range.length(1)},
+        make_array_view(this->get_executor(), storage_size, sub_range->data),
+        this->get_stride());
+}
+
+
+template <typename ValueType>
+std::unique_ptr<const MultiVector<ValueType>>
+MultiVector<ValueType>::create_subview_impl(local_span rows, local_span columns,
+                                            dim<2> global_size) const
+{
+    dim<2> actual_size{rows.length(), columns.length()};
+    GKO_ASSERT_EQUAL_DIMENSIONS(actual_size, global_size);
+    return const_cast<MultiVector&>(*this).create_subview(rows, columns);
+}
+
+
+template <typename ValueType>
+std::unique_ptr<const typename MultiVector<ValueType>::real_type>
+MultiVector<ValueType>::create_real_view_impl() const
+{
+    const auto num_rows = this->get_size()[0];
+    constexpr bool complex = is_complex<ValueType>();
+    const auto num_cols =
+        complex ? 2 * this->get_size()[1] : this->get_size()[1];
+    const auto stride = complex ? 2 * this->get_stride() : this->get_stride();
+
+    return MultiVector<remove_complex<ValueType>>::create_const(
+        this->get_executor(), dim<2>{num_rows, num_cols},
+        make_const_array_view(
+            this->get_executor(), num_rows * stride,
+            reinterpret_cast<const remove_complex<ValueType>*>(
+                this->get_const_values())),
+        stride);
+}
+
+
+template <typename ValueType>
+std::unique_ptr<typename MultiVector<ValueType>::real_type>
+MultiVector<ValueType>::create_real_view_impl()
+{
+    const auto num_rows = this->get_size()[0];
+    constexpr bool complex = is_complex<ValueType>();
+    const auto num_cols =
+        complex ? 2 * this->get_size()[1] : this->get_size()[1];
+    const auto stride = complex ? 2 * this->get_stride() : this->get_stride();
+
+    return MultiVector<remove_complex<ValueType>>::create(
+        this->get_executor(), dim<2>{num_rows, num_cols},
+        make_array_view(
+            this->get_executor(), num_rows * stride,
+            reinterpret_cast<remove_complex<ValueType>*>(this->get_values())),
+        stride);
+}
+
+
 template <typename ValueType>
 std::unique_ptr<typename MultiVector<ValueType>::absolute_type>
-MultiVector<ValueType>::compute_absolute() const
+MultiVector<ValueType>::compute_absolute_impl() const
 {
     // do not inherit the stride
     auto result = absolute_type::create(this->get_executor(), this->get_size());
@@ -1555,21 +1387,19 @@ MultiVector<ValueType>::compute_absolute() const
 
 
 template <typename ValueType>
-void MultiVector<ValueType>::compute_absolute(
-    ptr_param<absolute_type> output) const
+void MultiVector<ValueType>::compute_absolute_impl(absolute_type* result) const
 {
-    GKO_ASSERT_EQUAL_DIMENSIONS(this, output);
     auto exec = this->get_executor();
 
     exec->run(multivector::make_outplace_absolute_dense(
         this->get_const_device_view(),
-        make_temporary_output_clone(exec, output)->get_device_view()));
+        make_temporary_output_clone(exec, result)->get_device_view()));
 }
 
 
 template <typename ValueType>
 std::unique_ptr<typename MultiVector<ValueType>::complex_type>
-MultiVector<ValueType>::make_complex() const
+MultiVector<ValueType>::make_complex_impl() const
 {
     auto result = complex_type::create(this->get_executor(), this->get_size());
     this->make_complex(result);
@@ -1578,133 +1408,209 @@ MultiVector<ValueType>::make_complex() const
 
 
 template <typename ValueType>
-void MultiVector<ValueType>::make_complex(ptr_param<complex_type> result) const
+std::unique_ptr<typename MultiVector<ValueType>::real_type>
+MultiVector<ValueType>::get_real_impl() const
 {
-    GKO_ASSERT_EQUAL_DIMENSIONS(this, result);
-    auto exec = this->get_executor();
-
-    exec->run(multivector::make_make_complex(
-        this->get_const_device_view(),
-        make_temporary_output_clone(exec, result)->get_device_view()));
+    auto result = real_type::create(this->get_executor(), this->get_size());
+    this->get_real(result);
+    return result;
 }
 
 
 template <typename ValueType>
 std::unique_ptr<typename MultiVector<ValueType>::real_type>
-MultiVector<ValueType>::get_real() const
+MultiVector<ValueType>::get_imag_impl() const
 {
     auto result = real_type::create(this->get_executor(), this->get_size());
-    this->get_real(result);
+    this->get_imag(result);
     return result;
 }
 
 
 template <typename ValueType>
-void MultiVector<ValueType>::get_real(ptr_param<real_type> result) const
+void MultiVector<ValueType>::make_complex_impl(complex_type* result) const
 {
-    GKO_ASSERT_EQUAL_DIMENSIONS(this, result);
     auto exec = this->get_executor();
 
-    exec->run(multivector::make_get_real(
+    exec->run(multivector::make_make_complex(
         this->get_const_device_view(),
         make_temporary_output_clone(exec, result)->get_device_view()));
 }
 
 
 template <typename ValueType>
-std::unique_ptr<typename MultiVector<ValueType>::real_type>
-MultiVector<ValueType>::get_imag() const
+void MultiVector<ValueType>::get_real_impl(real_type* result) const
 {
-    auto result = real_type::create(this->get_executor(), this->get_size());
-    this->get_imag(result);
-    return result;
+    auto exec = this->get_executor();
+
+    exec->run(multivector::make_get_real(this->get_const_device_view(),
+                                         result->get_device_view()));
 }
 
 
 template <typename ValueType>
-void MultiVector<ValueType>::get_imag(ptr_param<real_type> result) const
+void MultiVector<ValueType>::get_imag_impl(real_type* result) const
 {
-    GKO_ASSERT_EQUAL_DIMENSIONS(this, result);
     auto exec = this->get_executor();
 
-    exec->run(multivector::make_get_imag(
-        this->get_const_device_view(),
-        make_temporary_output_clone(exec, result)->get_device_view()));
+    exec->run(multivector::make_get_imag(this->get_const_device_view(),
+                                         result->get_device_view()));
 }
 
 
 template <typename ValueType>
-auto MultiVector<ValueType>::get_device_view() -> device_view
+void MultiVector<ValueType>::fill_impl(value_type value)
 {
-    return device_view{this->get_size(), this->get_stride(),
-                       this->get_values()};
-};
+    this->get_executor()->run(
+        multivector::make_fill(this->get_device_view(), value));
+}
 
 
 template <typename ValueType>
-auto MultiVector<ValueType>::get_const_device_view() const -> const_device_view
+void MultiVector<ValueType>::scale_impl(scaling_param<value_type> alpha)
 {
-    return const_device_view{this->get_size(), this->get_stride(),
-                             this->get_const_values()};
-};
+    std::visit(
+        [this](auto alpha_v) {
+            auto exec = this->get_executor();
+            exec->run(multivector::make_scale(alpha_v->get_const_device_view(),
+                                              this->get_device_view()));
+        },
+        alpha.variant);
+}
 
 
 template <typename ValueType>
-std::unique_ptr<typename MultiVector<ValueType>::real_type>
-MultiVector<ValueType>::create_real_view()
+void MultiVector<ValueType>::inv_scale_impl(scaling_param<value_type> alpha)
 {
-    const auto num_rows = this->get_size()[0];
-    constexpr bool complex = is_complex<ValueType>();
-    const auto num_cols =
-        complex ? 2 * this->get_size()[1] : this->get_size()[1];
-    const auto stride = complex ? 2 * this->get_stride() : this->get_stride();
+    std::visit(
+        [this](auto alpha_v) {
+            auto exec = this->get_executor();
+            exec->run(multivector::make_inv_scale(
+                alpha_v->get_const_device_view(), this->get_device_view()));
+        },
+        alpha.variant);
+}
 
-    return MultiVector<remove_complex<ValueType>>::create(
-        this->get_executor(), dim<2>{num_rows, num_cols},
-        make_array_view(
-            this->get_executor(), num_rows * stride,
-            reinterpret_cast<remove_complex<ValueType>*>(this->get_values())),
-        stride);
+
+template <typename ValueType>
+void MultiVector<ValueType>::add_scaled_impl(scaling_param<value_type> alpha,
+                                             const MultiVector* b)
+{
+    std::visit(
+        [this, b](auto alpha_v) {
+            auto exec = this->get_executor();
+            exec->run(multivector::make_add_scaled(
+                alpha_v->get_const_device_view(), b->get_const_device_view(),
+                this->get_device_view()));
+        },
+        alpha.variant);
 }
 
 
 template <typename ValueType>
-std::unique_ptr<const typename MultiVector<ValueType>::real_type>
-MultiVector<ValueType>::create_real_view() const
+void MultiVector<ValueType>::sub_scaled_impl(scaling_param<value_type> alpha,
+                                             const MultiVector* b)
 {
-    const auto num_rows = this->get_size()[0];
-    constexpr bool complex = is_complex<ValueType>();
-    const auto num_cols =
-        complex ? 2 * this->get_size()[1] : this->get_size()[1];
-    const auto stride = complex ? 2 * this->get_stride() : this->get_stride();
+    std::visit(
+        [this, b](auto alpha_v) {
+            auto exec = this->get_executor();
 
-    return MultiVector<remove_complex<ValueType>>::create_const(
-        this->get_executor(), dim<2>{num_rows, num_cols},
-        make_const_array_view(
-            this->get_executor(), num_rows * stride,
-            reinterpret_cast<const remove_complex<ValueType>*>(
-                this->get_const_values())),
-        stride);
+            exec->run(multivector::make_sub_scaled(
+                alpha_v->get_const_device_view(), b->get_const_device_view(),
+                this->get_device_view()));
+        },
+        alpha.variant);
 }
 
 
 template <typename ValueType>
-std::unique_ptr<MultiVector<ValueType>>
-MultiVector<ValueType>::create_submatrix_impl(const span& rows,
-                                              const span& columns,
-                                              const size_type stride)
+void MultiVector<ValueType>::compute_dot_impl(
+    const MultiVector* b, matrix::MultiVector<value_type>* result,
+    array<char>& tmp) const
 {
-    row_major_range range_this{this->get_values(), this->get_size()[0],
-                               this->get_size()[1], this->get_stride()};
-    auto sub_range = range_this(rows, columns);
-    size_type storage_size =
-        rows.length() > 0 ? sub_range.length(1) +
-                                (sub_range.length(0) - 1) * this->get_stride()
-                          : 0;
-    return MultiVector::create(
-        this->get_executor(), dim<2>{sub_range.length(0), sub_range.length(1)},
-        make_array_view(this->get_executor(), storage_size, sub_range->data),
-        stride);
+    auto exec = this->get_executor();
+    if (tmp.get_executor() != exec) {
+        tmp.clear();
+        tmp.set_executor(exec);
+    }
+    exec->run(multivector::make_compute_dot(this->get_const_device_view(),
+                                            b->get_const_device_view(),
+                                            result->get_device_view(), tmp));
+}
+
+
+template <typename ValueType>
+void MultiVector<ValueType>::compute_conj_dot_impl(
+    const MultiVector* b, matrix::MultiVector<value_type>* result,
+    array<char>& tmp) const
+{
+    auto exec = this->get_executor();
+    if (tmp.get_executor() != exec) {
+        tmp.clear();
+        tmp.set_executor(exec);
+    }
+    exec->run(multivector::make_compute_conj_dot(
+        this->get_const_device_view(), b->get_const_device_view(),
+        result->get_device_view(), tmp));
+}
+
+
+template <typename ValueType>
+void MultiVector<ValueType>::compute_norm2_impl(norm_type* result,
+                                                array<char>& tmp) const
+{
+    auto exec = this->get_executor();
+    if (tmp.get_executor() != exec) {
+        tmp.clear();
+        tmp.set_executor(exec);
+    }
+    exec->run(multivector::make_compute_norm2(this->get_const_device_view(),
+                                              result->get_device_view(), tmp));
+}
+
+
+template <typename ValueType>
+void MultiVector<ValueType>::compute_squared_norm2_impl(norm_type* result,
+                                                        array<char>& tmp) const
+{
+    auto exec = this->get_executor();
+    if (tmp.get_executor() != exec) {
+        tmp.clear();
+        tmp.set_executor(exec);
+    }
+    exec->run(multivector::make_compute_squared_norm2(
+        this->get_const_device_view(), result->get_device_view(), tmp));
+}
+
+
+template <typename ValueType>
+void MultiVector<ValueType>::compute_norm1_impl(norm_type* result,
+                                                array<char>& tmp) const
+{
+    auto exec = this->get_executor();
+    if (tmp.get_executor() != exec) {
+        tmp.clear();
+        tmp.set_executor(exec);
+    }
+    exec->run(multivector::make_compute_norm1(this->get_const_device_view(),
+                                              result->get_device_view(), tmp));
+}
+
+
+template <typename ValueType>
+AbstractMultiVector::device_view<typename MultiVector<ValueType>::value_type>
+MultiVector<ValueType>::get_local_device_view_impl()
+{
+    return this->get_device_view();
+}
+
+
+template <typename ValueType>
+AbstractMultiVector::device_view<
+    const typename MultiVector<ValueType>::value_type>
+MultiVector<ValueType>::get_const_local_device_view_impl() const
+{
+    return this->get_const_device_view();
 }
 
 
@@ -1760,7 +1666,7 @@ std::unique_ptr<Dense<ValueType>> MultiVector<ValueType>::as_dense_view()
 template <typename ValueType>
 MultiVector<ValueType>::MultiVector(std::shared_ptr<const Executor> exec,
                                     const dim<2>& size, size_type stride)
-    : LinOp(exec, size, precision_v<ValueType>),
+    : EnableMultiVector<MultiVector>(exec, size),
       stride_(stride == 0 ? size[1] : stride),
       values_(exec, size[0] * stride_)
 {}
@@ -1770,7 +1676,7 @@ template <typename ValueType>
 MultiVector<ValueType>::MultiVector(std::shared_ptr<const Executor> exec,
                                     const dim<2>& size,
                                     array<value_type> values, size_type stride)
-    : LinOp(exec, size, precision_v<ValueType>),
+    : EnableMultiVector<MultiVector>(exec, size),
       stride_{stride},
       values_{exec, std::move(values)}
 {
diff --git a/core/test/utils/dummy_vector.hpp b/core/test/utils/dummy_vector.hpp
index 391283f20ef..57a68d1339a 100644
--- a/core/test/utils/dummy_vector.hpp
+++ b/core/test/utils/dummy_vector.hpp
@@ -28,11 +28,7 @@ class AbstractDummyVector : public gko::AbstractMultiVector {
 protected:
     Cloneable* copy_from_impl(const Cloneable* other) override
         GKO_NOT_IMPLEMENTED;
-    Cloneable* copy_from_impl(std::unique_ptr<Cloneable> other) override
-        GKO_NOT_IMPLEMENTED;
     Cloneable* move_from_impl(Cloneable* other) override GKO_NOT_IMPLEMENTED;
-    Cloneable* move_from_impl(std::unique_ptr<Cloneable> other) override
-        GKO_NOT_IMPLEMENTED;
     [[nodiscard]] std::unique_ptr<Cloneable> clone_impl(
         std::shared_ptr<const gko::Executor> exec) const override
         GKO_NOT_IMPLEMENTED;
diff --git a/include/ginkgo/core/base/polymorphic_object.hpp b/include/ginkgo/core/base/polymorphic_object.hpp
index b4d9541e31e..24a2d6e1a14 100644
--- a/include/ginkgo/core/base/polymorphic_object.hpp
+++ b/include/ginkgo/core/base/polymorphic_object.hpp
@@ -336,12 +336,8 @@ class Cloneable {
 protected:
     virtual Cloneable* copy_from_impl(const Cloneable* other) = 0;
 
-    virtual Cloneable* copy_from_impl(std::unique_ptr<Cloneable> other) = 0;
-
     virtual Cloneable* move_from_impl(Cloneable* other) = 0;
 
-    virtual Cloneable* move_from_impl(std::unique_ptr<Cloneable> other) = 0;
-
     [[nodiscard]] virtual std::unique_ptr<Cloneable> clone_impl(
         std::shared_ptr<const Executor> exec) const = 0;
 
@@ -431,20 +427,6 @@ class EnableCloneable : public ConvertibleTo<ConcreteType>, public Cloneable {
         return this;
     }
 
-    Cloneable* copy_from_impl(std::unique_ptr<Cloneable> other) override
-    {
-        self()->template log<log::Logger::polymorphic_object_copy_started>(
-            self()->get_executor().get(),
-            dynamic_cast<const PolymorphicObject*>(other.get()),
-            dynamic_cast<const PolymorphicObject*>(this));
-        as<ConvertibleTo<ConcreteType>>(other.get())->convert_to(self());
-        self()->template log<log::Logger::polymorphic_object_copy_completed>(
-            self()->get_executor().get(),
-            dynamic_cast<const PolymorphicObject*>(other.get()),
-            dynamic_cast<const PolymorphicObject*>(this));
-        return this;
-    }
-
     Cloneable* move_from_impl(Cloneable* other) override
     {
         self()->template log<log::Logger::polymorphic_object_move_started>(
@@ -459,20 +441,6 @@ class EnableCloneable : public ConvertibleTo<ConcreteType>, public Cloneable {
         return this;
     }
 
-    Cloneable* move_from_impl(std::unique_ptr<Cloneable> other) override
-    {
-        self()->template log<log::Logger::polymorphic_object_move_started>(
-            self()->get_executor().get(),
-            dynamic_cast<const PolymorphicObject*>(other.get()),
-            dynamic_cast<const PolymorphicObject*>(this));
-        as<ConvertibleTo<ConcreteType>>(other.get())->move_to(self());
-        self()->template log<log::Logger::polymorphic_object_move_completed>(
-            self()->get_executor().get(),
-            dynamic_cast<const PolymorphicObject*>(other.get()),
-            dynamic_cast<const PolymorphicObject*>(this));
-        return this;
-    }
-
     [[nodiscard]] std::unique_ptr<Cloneable> clone_impl(
         std::shared_ptr<const Executor> exec) const override
     {
diff --git a/include/ginkgo/core/matrix/multivector.hpp b/include/ginkgo/core/matrix/multivector.hpp
index 5b4922a2572..a5655e245ad 100644
--- a/include/ginkgo/core/matrix/multivector.hpp
+++ b/include/ginkgo/core/matrix/multivector.hpp
@@ -13,6 +13,7 @@
 #include <ginkgo/core/base/exception_helpers.hpp>
 #include <ginkgo/core/base/executor.hpp>
 #include <ginkgo/core/base/lin_op.hpp>
+#include <ginkgo/core/base/multivector_mixin.hpp>
 #include <ginkgo/core/base/range_accessors.hpp>
 #include <ginkgo/core/base/types.hpp>
 #include <ginkgo/core/base/utils.hpp>
@@ -20,6 +21,8 @@
 #include <ginkgo/core/matrix/permutation.hpp>
 #include <ginkgo/core/matrix/scaled_permutation.hpp>
 
+#include "ginkgo/core/base/multivector_mixin.hpp"
+
 
 namespace gko {
 namespace experimental {
@@ -67,8 +70,7 @@ class Dense;
  */
 template <typename ValueType = default_precision>
 class MultiVector
-    : public LinOp,
-      public EnableCloneable<MultiVector<ValueType>>,
+    : public EnableMultiVector<MultiVector<ValueType>>,
       public ConvertibleTo<MultiVector<next_precision<ValueType>>>,
 #if GINKGO_ENABLE_HALF || GINKGO_ENABLE_BFLOAT16
       public ConvertibleTo<MultiVector<next_precision<ValueType, 2>>>,
@@ -80,21 +82,17 @@ class MultiVector
       public ReadableFromMatrixData<ValueType, int32>,
       public ReadableFromMatrixData<ValueType, int64>,
       public WritableToMatrixData<ValueType, int32>,
-      public WritableToMatrixData<ValueType, int64>,
-      public Transposable,
-      public Permutable<int32>,
-      public Permutable<int64>,
-      public EnableAbsoluteComputation<remove_complex<MultiVector<ValueType>>> {
+      public WritableToMatrixData<ValueType, int64> {
     friend class Dense<ValueType>;
     friend class MultiVector<to_complex<ValueType>>;
-    friend class EnableCloneable<MultiVector>;
+    friend class EnableMultiVector<MultiVector>;
     friend class experimental::distributed::Vector<ValueType>;
     friend class experimental::distributed::detail::VectorCache<ValueType>;
     GKO_ASSERT_SUPPORTED_VALUE_TYPE;
 
 public:
-    using EnableCloneable<MultiVector>::convert_to;
-    using EnableCloneable<MultiVector>::move_to;
+    using EnableMultiVector<MultiVector>::convert_to;
+    using EnableMultiVector<MultiVector>::move_to;
     using ConvertibleTo<MultiVector<next_precision<ValueType>>>::convert_to;
     using ConvertibleTo<MultiVector<next_precision<ValueType>>>::move_to;
     using ConvertibleTo<Dense<ValueType>>::convert_to;
@@ -102,118 +100,28 @@ class MultiVector
     using ReadableFromMatrixData<ValueType, int32>::read;
     using ReadableFromMatrixData<ValueType, int64>::read;
 
-    using value_type = ValueType;
+    using value_type = typename EnableMultiVector<MultiVector>::value_type;
     using index_type = int64;
     using transposed_type = MultiVector<value_type>;
     using mat_data64 = matrix_data<value_type, int64>;
     using mat_data32 = matrix_data<value_type, int32>;
     using device_mat_data64 = device_matrix_data<value_type, int64>;
     using device_mat_data32 = device_matrix_data<value_type, int32>;
-    using absolute_type = remove_complex<MultiVector>;
-    using real_type = absolute_type;
-    using complex_type = to_complex<MultiVector>;
-    using device_view = matrix::view::dense<value_type>;
-    using const_device_view = matrix::view::dense<const value_type>;
+    using absolute_type =
+        typename EnableMultiVector<MultiVector>::absolute_type;
+    using real_type = typename EnableMultiVector<MultiVector>::real_type;
+    using complex_type = typename EnableMultiVector<MultiVector>::complex_type;
+    using norm_type = typename EnableMultiVector<MultiVector>::norm_type;
+    using device_view = typename EnableMultiVector<MultiVector>::device_view;
+    using const_device_view =
+        typename EnableMultiVector<MultiVector>::const_device_view;
 
     using row_major_range = gko::range<gko::accessor::row_major<ValueType, 2>>;
 
-    /**
-     * Creates a MultiVector with the same size and stride as another
-     * MultiVector matrix.
-     *
-     * @param other  The other matrix whose configuration needs to copied.
-     */
-    static std::unique_ptr<MultiVector> create_with_config_of(
-        ptr_param<const MultiVector> other)
-    {
-        // De-referencing `other` before calling the functions (instead of
-        // using operator `->`) is currently required to be compatible with
-        // CUDA 10.1.
-        // Otherwise, it results in a compile error.
-        return (*other).create_with_same_config();
-    }
-
-    /**
-     * Creates a MultiVector with the same type as another MultiVector
-     * matrix but on a different executor and with a different size.
-     *
-     * @param other  The other matrix whose type we target.
-     * @param exec  The executor of the new matrix.
-     * @param size  The size of the new matrix.
-     * @param stride  The stride of the new matrix.
-     *
-     * @returns a MultiVector with the type of other.
-     */
-    static std::unique_ptr<MultiVector> create_with_type_of(
-        ptr_param<const MultiVector> other,
-        std::shared_ptr<const Executor> exec, const dim<2>& size = dim<2>{})
-    {
-        // See create_with_config_of()
-        return (*other).create_with_type_of_impl(exec, size, size[1]);
-    }
-
-    /**
-     * @copydoc create_with_type_of(const MultiVector*, std::shared_ptr<const
-     * Executor>, const dim<2>)
-     *
-     * @param stride  The stride of the new matrix.
-     *
-     * @note This is an overload which allows full parameter specification.
-     */
-    static std::unique_ptr<MultiVector> create_with_type_of(
+    [[nodiscard]] static std::unique_ptr<MultiVector> create_with_type_of(
         ptr_param<const MultiVector> other,
         std::shared_ptr<const Executor> exec, const dim<2>& size,
-        size_type stride)
-    {
-        // See create_with_config_of()
-        return (*other).create_with_type_of_impl(exec, size, stride);
-    }
-
-    /**
-     * @copydoc create_with_type_of(const MultiVector*, std::shared_ptr<const
-     * Executor>, const dim<2>)
-     *
-     * @param local_size  Unused
-     * @param stride  The stride of the new matrix.
-     *
-     * @note This is an overload to stay consistent with
-     *       gko::experimental::distributed::Vector
-     */
-    static std::unique_ptr<MultiVector> create_with_type_of(
-        ptr_param<const MultiVector> other,
-        std::shared_ptr<const Executor> exec, const dim<2>& size,
-        const dim<2>& local_size, size_type stride)
-    {
-        // See create_with_config_of()
-        return (*other).create_with_type_of_impl(exec, size, stride);
-    }
-
-    /**
-     * Creates a MultiVector, where the underlying array is a view of another
-     * MultiVector' array.
-     *
-     * @param other  The other matrix on which to create the view
-     *
-     * @return  A MultiVector that is a view of other
-     */
-    static std::unique_ptr<MultiVector> create_view_of(
-        ptr_param<MultiVector> other)
-    {
-        return other->create_view_of_impl();
-    }
-
-    /**
-     * Creates a immutable MultiVector, where the underlying array is a view of
-     * another MultiVector' array.
-     *
-     * @param other  The other matrix on which to create the view
-     * @return  A immutable MultiVector that is a view of other
-     */
-    static std::unique_ptr<const MultiVector> create_const_view_of(
-        ptr_param<const MultiVector> other)
-    {
-        return other->create_const_view_of_impl();
-    }
+        size_type stride);
 
     friend class MultiVector<previous_precision<ValueType>>;
 
@@ -264,9 +172,9 @@ class MultiVector
 
     void write(mat_data32& data) const override;
 
-    std::unique_ptr<LinOp> transpose() const override;
+    std::unique_ptr<MultiVector> transpose() const;
 
-    std::unique_ptr<LinOp> conj_transpose() const override;
+    std::unique_ptr<MultiVector> conj_transpose() const;
 
     /**
      * Writes the transposed matrix into the given output matrix.
@@ -284,13 +192,6 @@ class MultiVector
      */
     void conj_transpose(ptr_param<MultiVector> output) const;
 
-    /**
-     * Fill the dense matrix with a given value.
-     *
-     * @param value  the value to be filled
-     */
-    void fill(const ValueType value);
-
     /**
      * Creates a permuted copy $A'$ of this matrix $A$ with the given
      * permutation $P$. By default, this computes a symmetric permutation
@@ -466,11 +367,11 @@ class MultiVector
             column_permutation,
         ptr_param<MultiVector> output, bool invert = false) const;
 
-    std::unique_ptr<LinOp> permute(
-        const array<int32>* permutation_indices) const override;
+    std::unique_ptr<MultiVector> permute(
+        const array<int32>* permutation_indices) const;
 
-    std::unique_ptr<LinOp> permute(
-        const array<int64>* permutation_indices) const override;
+    std::unique_ptr<MultiVector> permute(
+        const array<int64>* permutation_indices) const;
 
     /**
      * Writes the symmetrically permuted matrix into the given output matrix.
@@ -490,11 +391,11 @@ class MultiVector
     void permute(const array<int64>* permutation_indices,
                  ptr_param<MultiVector> output) const;
 
-    std::unique_ptr<LinOp> inverse_permute(
-        const array<int32>* permutation_indices) const override;
+    std::unique_ptr<MultiVector> inverse_permute(
+        const array<int32>* permutation_indices) const;
 
-    std::unique_ptr<LinOp> inverse_permute(
-        const array<int64>* permutation_indices) const override;
+    std::unique_ptr<MultiVector> inverse_permute(
+        const array<int64>* permutation_indices) const;
 
     /**
      * Writes the inverse symmetrically permuted matrix into the given output
@@ -515,11 +416,11 @@ class MultiVector
     void inverse_permute(const array<int64>* permutation_indices,
                          ptr_param<MultiVector> output) const;
 
-    std::unique_ptr<LinOp> row_permute(
-        const array<int32>* permutation_indices) const override;
+    std::unique_ptr<MultiVector> row_permute(
+        const array<int32>* permutation_indices) const;
 
-    std::unique_ptr<LinOp> row_permute(
-        const array<int64>* permutation_indices) const override;
+    std::unique_ptr<MultiVector> row_permute(
+        const array<int64>* permutation_indices) const;
 
     /**
      * Writes the row-permuted matrix into the given output matrix.
@@ -572,13 +473,13 @@ class MultiVector
      *                        matrix and `gather_indices->get_size()` rows.
      */
     void row_gather(const array<int32>* gather_indices,
-                    ptr_param<LinOp> row_collection) const;
+                    ptr_param<AbstractMultiVector> row_collection) const;
 
     /**
      * @copydoc row_gather(const array<int32>*, LinOp*) const
      */
     void row_gather(const array<int64>* gather_indices,
-                    ptr_param<LinOp> row_collection) const;
+                    ptr_param<AbstractMultiVector> row_collection) const;
 
     /**
      * Copies the given rows from this matrix into `row_collection` with scaling
@@ -596,7 +497,7 @@ class MultiVector
     void row_gather(ptr_param<const LinOp> alpha,
                     const array<int32>* gather_indices,
                     ptr_param<const LinOp> beta,
-                    ptr_param<LinOp> row_collection) const;
+                    ptr_param<AbstractMultiVector> row_collection) const;
 
     /**
      * @copydoc row_gather(const LinOp*, const array<int32>*, const LinOp*,
@@ -605,13 +506,13 @@ class MultiVector
     void row_gather(ptr_param<const LinOp> alpha,
                     const array<int64>* gather_indices,
                     ptr_param<const LinOp> beta,
-                    ptr_param<LinOp> row_collection) const;
+                    ptr_param<AbstractMultiVector> row_collection) const;
 
-    std::unique_ptr<LinOp> column_permute(
-        const array<int32>* permutation_indices) const override;
+    std::unique_ptr<MultiVector> column_permute(
+        const array<int32>* permutation_indices) const;
 
-    std::unique_ptr<LinOp> column_permute(
-        const array<int64>* permutation_indices) const override;
+    std::unique_ptr<MultiVector> column_permute(
+        const array<int64>* permutation_indices) const;
 
     /**
      * Writes the column-permuted matrix into the given output matrix.
@@ -631,11 +532,11 @@ class MultiVector
     void column_permute(const array<int64>* permutation_indices,
                         ptr_param<MultiVector> output) const;
 
-    std::unique_ptr<LinOp> inverse_row_permute(
-        const array<int32>* permutation_indices) const override;
+    std::unique_ptr<MultiVector> inverse_row_permute(
+        const array<int32>* permutation_indices) const;
 
-    std::unique_ptr<LinOp> inverse_row_permute(
-        const array<int64>* permutation_indices) const override;
+    std::unique_ptr<MultiVector> inverse_row_permute(
+        const array<int64>* permutation_indices) const;
 
     /**
      * Writes the inverse row-permuted matrix into the given output matrix.
@@ -656,11 +557,11 @@ class MultiVector
     void inverse_row_permute(const array<int64>* permutation_indices,
                              ptr_param<MultiVector> output) const;
 
-    std::unique_ptr<LinOp> inverse_column_permute(
-        const array<int32>* permutation_indices) const override;
+    std::unique_ptr<MultiVector> inverse_column_permute(
+        const array<int32>* permutation_indices) const;
 
-    std::unique_ptr<LinOp> inverse_column_permute(
-        const array<int64>* permutation_indices) const override;
+    std::unique_ptr<MultiVector> inverse_column_permute(
+        const array<int64>* permutation_indices) const;
 
     /**
      * Writes the inverse column-permuted matrix into the given output matrix.
@@ -681,55 +582,6 @@ class MultiVector
     void inverse_column_permute(const array<int64>* permutation_indices,
                                 ptr_param<MultiVector> output) const;
 
-    std::unique_ptr<absolute_type> compute_absolute() const override;
-
-    /**
-     * Writes the absolute values of this matrix into an existing matrix.
-     *
-     * @param output  The output matrix. Its size must match the size of this
-     *                matrix.
-     * @see MultiVector::compute_absolute()
-     */
-    void compute_absolute(ptr_param<absolute_type> output) const;
-
-    void compute_absolute_inplace() override;
-
-    /**
-     * Creates a complex copy of the original matrix. If the original matrix
-     * was real, the imaginary part of the result will be zero.
-     */
-    std::unique_ptr<complex_type> make_complex() const;
-
-    /**
-     * Writes a complex copy of the original matrix to a given complex matrix.
-     * If the original matrix was real, the imaginary part of the result will
-     * be zero.
-     */
-    void make_complex(ptr_param<complex_type> result) const;
-
-    /**
-     * Creates a new real matrix and extracts the real part of the original
-     * matrix into that.
-     */
-    std::unique_ptr<real_type> get_real() const;
-
-    /**
-     * Extracts the real part of the original matrix into a given real matrix.
-     */
-    void get_real(ptr_param<real_type> result) const;
-
-    /**
-     * Creates a new real matrix and extracts the imaginary part of the
-     * original matrix into that.
-     */
-    std::unique_ptr<real_type> get_imag() const;
-
-    /**
-     * Extracts the imaginary part of the original matrix into a given real
-     * matrix.
-     */
-    void get_imag(ptr_param<real_type> result) const;
-
     /**
      * Returns a pointer to the array of values of the matrix.
      *
@@ -820,166 +672,6 @@ class MultiVector
         return values_.get_const_data()[linearize_index(idx)];
     }
 
-    /**
-     * Scales the matrix with a scalar (aka: BLAS scal).
-     *
-     * @param alpha  If alpha is 1x1 MultiVector is scaled
-     *               by alpha. If it is a MultiVector row vector of values,
-     *               then i-th column of the matrix is scaled with the i-th
-     *               element of alpha (the number of columns of alpha has to
-     *               match the number of columns of the matrix).
-     */
-    void scale(ptr_param<const LinOp> alpha);
-
-    /**
-     * Scales the matrix with the inverse of a scalar.
-     *
-     * @param alpha  If alpha is 1x1 MultiVector is scaled
-     *               by 1 / alpha. If it is a MultiVector row vector of values,
-     *               then i-th column of the matrix is scaled with the inverse
-     *               of the i-th element of alpha (the number of columns of
-     *               alpha has to match the number of columns of the matrix).
-     */
-    void inv_scale(ptr_param<const LinOp> alpha);
-
-    /**
-     * Adds `b` scaled by `alpha` to the matrix (aka: BLAS axpy).
-     *
-     * @param alpha  If alpha is 1x1 MultiVector is scaled
-     *               by alpha. If it is a MultiVector row vector of values,
-     *               then i-th column of the matrix is scaled with the i-th
-     *               element of alpha (the number of columns of alpha has to
-     *               match the number of columns of the matrix).
-     * @param b  a matrix of the same dimension as this
-     */
-    void add_scaled(ptr_param<const LinOp> alpha, ptr_param<const LinOp> b);
-
-    /**
-     * Subtracts `b` scaled by `alpha` from the matrix (aka: BLAS axpy).
-     *
-     * @param alpha  If alpha is 1x1 MultiVector, b is scaled
-     *               by alpha. If it is a MultiVector row vector of values,
-     *               then i-th column of b is scaled with the i-th
-     *               element of alpha (the number of columns of alpha has to
-     *               match the number of columns of the matrix).
-     * @param b  a matrix of the same dimension as this
-     */
-    void sub_scaled(ptr_param<const LinOp> alpha, ptr_param<const LinOp> b);
-
-    /**
-     * Computes the column-wise dot product of this matrix and `b`.
-     *
-     * @param b  a MultiVector of same dimension as this
-     * @param result  a MultiVector row vector, used to store the dot product
-     *                (the number of column in the vector must match the number
-     *                of columns of this)
-     */
-    void compute_dot(ptr_param<const LinOp> b, ptr_param<LinOp> result) const;
-
-    /**
-     * Computes the column-wise dot product of this matrix and `b`.
-     *
-     * @param b  a MultiVector of same dimension as this
-     * @param result  a MultiVector row vector, used to store the dot product
-     *                (the number of column in the vector must match the number
-     *                of columns of this)
-     * @param tmp  the temporary storage to use for partial sums during the
-     *             reduction computation. It may be resized and/or reset to the
-     *             correct executor.
-     */
-    void compute_dot(ptr_param<const LinOp> b, ptr_param<LinOp> result,
-                     array<char>& tmp) const;
-
-    /**
-     * Computes the column-wise dot product of `conj(this matrix)` and `b`.
-     *
-     * @param b  a MultiVector of same dimension as this
-     * @param result  a MultiVector row vector, used to store the dot product
-     *                (the number of column in the vector must match the number
-     *                of columns of this)
-     */
-    void compute_conj_dot(ptr_param<const LinOp> b,
-                          ptr_param<LinOp> result) const;
-
-    /**
-     * Computes the column-wise dot product of `conj(this matrix)` and `b`.
-     *
-     * @param b  a MultiVector of same dimension as this
-     * @param result  a MultiVector row vector, used to store the dot product
-     *                (the number of column in the vector must match the number
-     *                of columns of this)
-     * @param tmp  the temporary storage to use for partial sums during the
-     *             reduction computation. It may be resized and/or reset to the
-     *             correct executor.
-     */
-    void compute_conj_dot(ptr_param<const LinOp> b, ptr_param<LinOp> result,
-                          array<char>& tmp) const;
-
-    /**
-     * Computes the column-wise Euclidean (L^2) norm of this matrix.
-     *
-     * @param result  a MultiVector row vector, used to store the norm
-     *                (the number of columns in the vector must match the number
-     *                of columns of this)
-     */
-    void compute_norm2(ptr_param<LinOp> result) const;
-
-    /**
-     * Computes the column-wise Euclidean (L^2) norm of this matrix.
-     *
-     * @param result  a MultiVector row vector, used to store the norm
-     *                (the number of columns in the vector must match the
-     *                number of columns of this)
-     * @param tmp  the temporary storage to use for partial sums during the
-     *             reduction computation. It may be resized and/or reset to the
-     *             correct executor.
-     */
-    void compute_norm2(ptr_param<LinOp> result, array<char>& tmp) const;
-
-    /**
-     * Computes the column-wise (L^1) norm of this matrix.
-     *
-     * @param result  a MultiVector row vector, used to store the norm
-     *                (the number of columns in the vector must match the number
-     *                of columns of this)
-     */
-    void compute_norm1(ptr_param<LinOp> result) const;
-
-    /**
-     * Computes the column-wise (L^1) norm of this matrix.
-     *
-     * @param result  a MultiVector row vector, used to store the norm
-     *                (the number of columns in the vector must match the
-     *                number of columns of this)
-     * @param tmp  the temporary storage to use for partial sums during the
-     *             reduction computation. It may be resized and/or reset to the
-     *             correct executor.
-     */
-    void compute_norm1(ptr_param<LinOp> result, array<char>& tmp) const;
-
-    /**
-     * Computes the square of the column-wise Euclidean (L^2) norm of this
-     * matrix.
-     *
-     * @param result  a MultiVector row vector, used to store the norm
-     *                (the number of columns in the vector must match the number
-     *                of columns of this)
-     */
-    void compute_squared_norm2(ptr_param<LinOp> result) const;
-
-    /**
-     * Computes the square of the column-wise Euclidean (L^2) norm of this
-     * matrix.
-     *
-     * @param result  a MultiVector row vector, used to store the norm
-     *                (the number of columns in the vector must match the
-     *                number of columns of this)
-     * @param tmp  the temporary storage to use for partial sums during the
-     *             reduction computation. It may be resized and/or reset to the
-     *             correct executor.
-     */
-    void compute_squared_norm2(ptr_param<LinOp> result, array<char>& tmp) const;
-
     /**
      * Computes the column-wise arithmetic mean of this matrix.
      *
@@ -987,7 +679,7 @@ class MultiVector
      *                (the number of columns in the vector must match the number
      *                of columns of this)
      */
-    void compute_mean(ptr_param<LinOp> result) const;
+    void compute_mean(ptr_param<AbstractMultiVector> result) const;
 
     /**
      * Computes the column-wise arithmetic mean of this matrix.
@@ -999,68 +691,8 @@ class MultiVector
      *             reduction computation. It may be resized and/or reset to the
      *             correct executor.
      */
-    void compute_mean(ptr_param<LinOp> result, array<char>& tmp) const;
-
-    /**
-     * Create a submatrix from the original matrix.
-     * Warning: defining stride for this create_submatrix method might cause
-     * wrong memory access. Better use the create_submatrix(rows, columns)
-     * method instead.
-     *
-     * @param rows     row span
-     * @param columns  column span
-     * @param stride   stride of the new submatrix.
-     */
-    std::unique_ptr<MultiVector> create_submatrix(const span& rows,
-                                                  const span& columns,
-                                                  const size_type stride)
-    {
-        return this->create_submatrix_impl(rows, columns, stride);
-    }
-
-    /**
-     * Create a submatrix from the original matrix.
-     *
-     * @param rows     row span
-     * @param columns  column span
-     */
-    std::unique_ptr<MultiVector> create_submatrix(const span& rows,
-                                                  const span& columns)
-    {
-        return create_submatrix(rows, columns, this->get_stride());
-    }
-
-
-    /**
-     * Create a submatrix from the original matrix.
-     *
-     * @param rows  row span
-     * @param columns  column span
-     * @param size  size of the submatrix (only used for consistency with
-     *              distributed::Vector)
-     */
-    std::unique_ptr<MultiVector> create_submatrix(const local_span& rows,
-                                                  const local_span& columns,
-                                                  dim<2> size)
-    {
-        dim<2> deduced_size{rows.length(), columns.length()};
-        GKO_ASSERT_EQUAL_DIMENSIONS(deduced_size, size);
-        return create_submatrix(rows, columns, this->get_stride());
-    }
-
-    /**
-     * Create a real view of the (potentially) complex original matrix.
-     * If the original matrix is real, nothing changes. If the original matrix
-     * is complex, the result is created by viewing the complex matrix with as
-     * real with a reinterpret_cast with twice the number of columns and
-     * double the stride.
-     */
-    std::unique_ptr<real_type> create_real_view();
-
-    /**
-     * @copydoc create_real_view()
-     */
-    std::unique_ptr<const real_type> create_real_view() const;
+    void compute_mean(ptr_param<AbstractMultiVector> result,
+                      array<char>& tmp) const;
 
     /**
      * Creates an uninitialized MultiVector of the specified size.
@@ -1135,6 +767,14 @@ class MultiVector
 
     [[nodiscard]] std::unique_ptr<Dense<ValueType>> as_dense_view();
 
+    template <typename OtherValueType>
+    [[nodiscard]] temporary_conversion<MultiVector<OtherValueType>>
+    as_precision();
+
+    template <typename OtherValueType>
+    [[nodiscard]] temporary_conversion<const MultiVector<OtherValueType>>
+    as_precision() const;
+
     /**
      * Copy-assigns a MultiVector. Preserves the executor, reallocates
      * the matrix with minimal stride if the dimensions don't match, then copies
@@ -1168,19 +808,6 @@ class MultiVector
     MultiVector(std::shared_ptr<const Executor> exec, const dim<2>& size,
                 array<value_type> values, size_type stride);
 
-    /**
-     * Creates a MultiVector with the same size and stride as the callers
-     * matrix.
-     *
-     * @returns a MultiVector with the same size and stride as the
-     * caller.
-     */
-    virtual std::unique_ptr<MultiVector> create_with_same_config() const
-    {
-        return MultiVector::create(this->get_executor(), this->get_size(),
-                                   this->get_stride());
-    }
-
     /**
      * Creates a MultiVector with the same type as the callers matrix.
      *
@@ -1195,114 +822,10 @@ class MultiVector
         return MultiVector::create(exec, size, stride);
     }
 
-    /**
-     * Creates a MultiVector where the underlying array is a view of
-     * this' array.
-     *
-     * @return  A MultiVector that is a view of this.
-     */
-    virtual std::unique_ptr<MultiVector> create_view_of_impl()
-    {
-        auto exec = this->get_executor();
-        return MultiVector::create(
-            exec, this->get_size(),
-            gko::make_array_view(exec, this->get_num_stored_elements(),
-                                 this->get_values()),
-            this->get_stride());
-    }
-
-    /**
-     * Creates a immutable MultiVector where the underlying array is a
-     * view of this' array.
-     *
-     * @return  A immutable MultiVector that is a view of this.
-     */
-    virtual std::unique_ptr<const MultiVector> create_const_view_of_impl() const
-    {
-        auto exec = this->get_executor();
-        return MultiVector::create_const(
-            exec, this->get_size(),
-            gko::make_const_array_view(exec, this->get_num_stored_elements(),
-                                       this->get_const_values()),
-            this->get_stride());
-    }
-
-    /**
-     * @copydoc scale(const LinOp *)
-     *
-     * @deprecated  This function will be removed in the future,
-     *              we will instead always use Ginkgo's implementation.
-     */
-    virtual void scale_impl(const LinOp* alpha);
-
-    /**
-     * @copydoc inv_scale(const LinOp *)
-     *
-     * @deprecated  This function will be removed in the future,
-     *              we will instead always use Ginkgo's implementation.
-     */
-    virtual void inv_scale_impl(const LinOp* alpha);
-
-    /**
-     * @copydoc add_scaled(const LinOp *, const LinOp *)
-     *
-     * @deprecated  This function will be removed in the future,
-     *              we will instead always use Ginkgo's implementation.
-     */
-    virtual void add_scaled_impl(const LinOp* alpha, const LinOp* b);
-
-    /**
-     * @copydoc sub_scaled(const LinOp *, const LinOp *)
-     *
-     * @deprecated  This function will be removed in the future,
-     *              we will instead always use Ginkgo's implementation.
-     */
-    virtual void sub_scaled_impl(const LinOp* alpha, const LinOp* b);
-
-    /**
-     * @copydoc compute_dot(const LinOp*, LinOp*) const
-     *
-     * @deprecated  This function will be removed in the future,
-     *              we will instead always use Ginkgo's implementation.
-     */
-    virtual void compute_dot_impl(const LinOp* b, LinOp* result) const;
-
-    /**
-     * @copydoc compute_conj_dot(const LinOp*, LinOp*) const
-     *
-     * @deprecated  This function will be removed in the future,
-     *              we will instead always use Ginkgo's implementation.
-     */
-    virtual void compute_conj_dot_impl(const LinOp* b, LinOp* result) const;
-
-    /**
-     * @copydoc compute_norm2(LinOp*) const
-     *
-     * @deprecated  This function will be removed in the future,
-     *              we will instead always use Ginkgo's implementation.
-     */
-    virtual void compute_norm2_impl(LinOp* result) const;
-
-    /**
-     * @copydoc compute_norm1(LinOp*) const
-     *
-     * @deprecated  This function will be removed in the future,
-     *              we will instead always use Ginkgo's implementation.
-     */
-    virtual void compute_norm1_impl(LinOp* result) const;
-
-    /**
-     * @copydoc compute_squared_norm2(LinOp*) const
-     *
-     * @deprecated  This function will be removed in the future,
-     *              we will instead always use Ginkgo's implementation.
-     */
-    virtual void compute_squared_norm2_impl(LinOp* result) const;
-
     /**
      * @copydoc compute_mean(LinOp*) const
      */
-    virtual void compute_mean_impl(LinOp* result) const;
+    virtual void compute_mean_impl(AbstractMultiVector* result) const;
 
     /**
      * Resizes the matrix to the given size.
@@ -1314,21 +837,6 @@ class MultiVector
      */
     void resize(gko::dim<2> new_size);
 
-    /**
-     * @copydoc create_submatrix(const span, const span, const size_type)
-     *
-     * @note  Other implementations of dense should override this function
-     *        instead of create_submatrix(const span, const span, const
-     *        size_type).
-     */
-    virtual std::unique_ptr<MultiVector> create_submatrix_impl(
-        const span& rows, const span& columns, const size_type stride);
-
-    void apply_impl(const LinOp* b, LinOp* x) const override;
-
-    void apply_impl(const LinOp* alpha, const LinOp* b, const LinOp* beta,
-                    LinOp* x) const override;
-
     size_type linearize_index(size_type row, size_type col) const noexcept
     {
         return row * stride_ + col;
@@ -1370,6 +878,82 @@ class MultiVector
                          const MultiVector<ValueType>* beta,
                          MultiVector<OutputType>* row_collection) const;
 
+    void compute_absolute_inplace_impl() override;
+
+    [[nodiscard]] std::unique_ptr<MultiVector> create_with_same_config_impl()
+        const override;
+
+    [[nodiscard]] std::unique_ptr<MultiVector> create_with_type_of_impl(
+        std::shared_ptr<const Executor> exec, const dim<2>& global_size,
+        const dim<2>& local_size, size_type stride) const override;
+
+    [[nodiscard]] std::unique_ptr<MultiVector> create_subview_impl(
+        local_span rows, local_span columns) override;
+
+    [[nodiscard]] std::unique_ptr<const MultiVector> create_subview_impl(
+        local_span rows, local_span columns) const override;
+
+    [[nodiscard]] std::unique_ptr<MultiVector> create_subview_impl(
+        local_span rows, local_span columns, dim<2> global_size) override;
+
+    [[nodiscard]] std::unique_ptr<const MultiVector> create_subview_impl(
+        local_span rows, local_span columns, dim<2> global_size) const override;
+
+    [[nodiscard]] std::unique_ptr<const real_type> create_real_view_impl()
+        const override;
+
+    [[nodiscard]] std::unique_ptr<real_type> create_real_view_impl() override;
+
+    [[nodiscard]] std::unique_ptr<absolute_type> compute_absolute_impl()
+        const override;
+
+    void compute_absolute_impl(absolute_type* result) const override;
+
+    [[nodiscard]] std::unique_ptr<complex_type> make_complex_impl()
+        const override;
+
+    [[nodiscard]] std::unique_ptr<real_type> get_real_impl() const override;
+
+    [[nodiscard]] std::unique_ptr<real_type> get_imag_impl() const override;
+
+    void make_complex_impl(complex_type* result) const override;
+
+    void get_real_impl(real_type* result) const override;
+
+    void get_imag_impl(real_type* result) const override;
+
+    void fill_impl(value_type value) override;
+
+    void scale_impl(scaling_param<value_type> alpha) override;
+
+    void inv_scale_impl(scaling_param<value_type> alpha) override;
+
+    void add_scaled_impl(scaling_param<value_type> alpha,
+                         const MultiVector* b) override;
+
+    void sub_scaled_impl(scaling_param<value_type> alpha,
+                         const MultiVector* b) override;
+
+    void compute_dot_impl(const MultiVector* b,
+                          matrix::MultiVector<value_type>* result,
+                          array<char>& tmp) const override;
+
+    void compute_conj_dot_impl(const MultiVector* b,
+                               matrix::MultiVector<value_type>* result,
+                               array<char>& tmp) const override;
+
+    void compute_norm2_impl(norm_type* result, array<char>& tmp) const override;
+
+    void compute_squared_norm2_impl(norm_type* result,
+                                    array<char>& tmp) const override;
+
+    void compute_norm1_impl(norm_type* result, array<char>& tmp) const override;
+    AbstractMultiVector::device_view<value_type> get_local_device_view_impl()
+        override;
+
+    AbstractMultiVector::device_view<const value_type>
+    get_const_local_device_view_impl() const override;
+
 private:
     size_type stride_;
     array<value_type> values_;
@@ -1413,8 +997,8 @@ std::unique_ptr<
     matrix::MultiVector<typename detail::pointee<VecPtr>::value_type>>
 make_dense_view(VecPtr&& vector)
 {
-    using value_type = typename detail::pointee<VecPtr>::value_type;
-    return matrix::MultiVector<value_type>::create_view_of(vector);
+    return vector->create_subview({0, vector->get_size()[0]},
+                                  {0, vector->get_size()[1]});
 }
 
 
@@ -1430,8 +1014,8 @@ std::unique_ptr<
     const matrix::MultiVector<typename detail::pointee<VecPtr>::value_type>>
 make_const_dense_view(VecPtr&& vector)
 {
-    using value_type = typename detail::pointee<VecPtr>::value_type;
-    return matrix::MultiVector<value_type>::create_const_view_of(vector);
+    return vector->create_subview({0, vector->get_size()[0]},
+                                  {0, vector->get_size()[1]});
 }
 
 
diff --git a/test/matrix/multivector_kernels.cpp b/test/matrix/multivector_kernels.cpp
index 97a5332e3f1..7a6889ba752 100644
--- a/test/matrix/multivector_kernels.cpp
+++ b/test/matrix/multivector_kernels.cpp
@@ -413,10 +413,10 @@ TEST_F(MultiVector, IsTransposable)
 TEST_F(MultiVector, IsTransposableIntoMultiVectorCrossExecutor)
 {
     set_up_apply_data();
-    auto row_span = gko::span{0, x->get_size()[0] - 2};
-    auto col_span = gko::span{0, x->get_size()[1] - 2};
-    auto sub_x = x->create_submatrix(row_span, col_span);
-    auto sub_dx = dx->create_submatrix(row_span, col_span);
+    auto row_span = gko::local_span{0, x->get_size()[0] - 2};
+    auto col_span = gko::local_span{0, x->get_size()[1] - 2};
+    auto sub_x = x->create_subview(row_span, col_span);
+    auto sub_dx = dx->create_subview(row_span, col_span);
     // create the target matrices on another executor to
     // force temporary clone
     auto trans = Mtx::create(ref, gko::transpose(sub_x->get_size()));
@@ -445,10 +445,10 @@ TEST_F(MultiVector, IsConjugateTransposable)
 TEST_F(MultiVector, IsConjugateTransposableIntoMultiVectorCrossExecutor)
 {
     set_up_apply_data();
-    auto row_span = gko::span{0, c_x->get_size()[0] - 2};
-    auto col_span = gko::span{0, c_x->get_size()[1] - 2};
-    auto sub_x = c_x->create_submatrix(row_span, col_span);
-    auto sub_dx = dc_x->create_submatrix(row_span, col_span);
+    auto row_span = gko::local_span{0, c_x->get_size()[0] - 2};
+    auto col_span = gko::local_span{0, c_x->get_size()[1] - 2};
+    auto sub_x = c_x->create_subview(row_span, col_span);
+    auto sub_dx = dc_x->create_subview(row_span, col_span);
     // create the target matrices on another executor to
     // force temporary clone
     auto trans = ComplexMtx::create(ref, gko::transpose(sub_x->get_size()));
@@ -927,10 +927,10 @@ TEST_F(MultiVector, CanGatherRows)
 TEST_F(MultiVector, CanGatherRowsIntoMultiVectorCrossExecutor)
 {
     set_up_apply_data();
-    auto row_span = gko::span{0, x->get_size()[0]};
-    auto col_span = gko::span{0, x->get_size()[1] - 2};
-    auto sub_x = x->create_submatrix(row_span, col_span);
-    auto sub_dx = dx->create_submatrix(row_span, col_span);
+    auto row_span = gko::local_span{0, x->get_size()[0]};
+    auto col_span = gko::local_span{0, x->get_size()[1] - 2};
+    auto sub_x = x->create_subview(row_span, col_span);
+    auto sub_dx = dx->create_subview(row_span, col_span);
     auto gather_size =
         gko::dim<2>{rgather_idxs->get_size(), sub_x->get_size()[1]};
     auto r_gather = Mtx::create(ref, gather_size);
@@ -947,10 +947,10 @@ TEST_F(MultiVector, CanGatherRowsIntoMultiVectorCrossExecutor)
 TEST_F(MultiVector, CanAdvancedGatherRowsIntoMultiVectorCrossExecutor)
 {
     set_up_apply_data();
-    auto row_span = gko::span{0, x->get_size()[0]};
-    auto col_span = gko::span{0, x->get_size()[1] - 2};
-    auto sub_x = x->create_submatrix(row_span, col_span);
-    auto sub_dx = dx->create_submatrix(row_span, col_span);
+    auto row_span = gko::local_span{0, x->get_size()[0]};
+    auto col_span = gko::local_span{0, x->get_size()[1] - 2};
+    auto sub_x = x->create_subview(row_span, col_span);
+    auto sub_dx = dx->create_subview(row_span, col_span);
     auto gather_size =
         gko::dim<2>{rgather_idxs->get_size(), sub_x->get_size()[1]};
     auto r_gather = gen_mtx<Mtx>(gather_size[0], gather_size[1]);
@@ -968,10 +968,10 @@ TEST_F(MultiVector, CanAdvancedGatherRowsIntoMultiVectorCrossExecutor)
 TEST_F(MultiVector, CanGatherRowsIntoMixedMultiVectorCrossExecutor)
 {
     set_up_apply_data();
-    auto row_span = gko::span{0, x->get_size()[0]};
-    auto col_span = gko::span{0, x->get_size()[1] - 2};
-    auto sub_x = x->create_submatrix(row_span, col_span);
-    auto sub_dx = dx->create_submatrix(row_span, col_span);
+    auto row_span = gko::local_span{0, x->get_size()[0]};
+    auto col_span = gko::local_span{0, x->get_size()[1] - 2};
+    auto sub_x = x->create_subview(row_span, col_span);
+    auto sub_dx = dx->create_subview(row_span, col_span);
     auto gather_size =
         gko::dim<2>{rgather_idxs->get_size(), sub_x->get_size()[1]};
     auto r_gather = MixedMtx::create(ref, gather_size);
@@ -989,10 +989,10 @@ TEST_F(MultiVector, CanGatherRowsIntoMixedMultiVectorCrossExecutor)
 TEST_F(MultiVector, CanAdvancedGatherRowsIntoMixedMultiVectorCrossExecutor)
 {
     set_up_apply_data();
-    auto row_span = gko::span{0, x->get_size()[0]};
-    auto col_span = gko::span{0, x->get_size()[1] - 2};
-    auto sub_x = x->create_submatrix(row_span, col_span);
-    auto sub_dx = dx->create_submatrix(row_span, col_span);
+    auto row_span = gko::local_span{0, x->get_size()[0]};
+    auto col_span = gko::local_span{0, x->get_size()[1] - 2};
+    auto sub_x = x->create_subview(row_span, col_span);
+    auto sub_dx = dx->create_subview(row_span, col_span);
     auto gather_size =
         gko::dim<2>{rgather_idxs->get_size(), sub_x->get_size()[1]};
     auto r_gather = gen_mtx<MixedMtx>(gather_size[0], gather_size[1]);