Merge Make GMRES work with custom Dense matrices

Add new functionalities and capabilities to Dense

+ Give `create_submatrix` a virtual behavior.
+ Add a new `create_with_type_of` to allow virtual
creation of Dense matrices.
+ Use `create_with_type_of` in GMRES to so that
it works with custom Dense matrices.

Related PR: #718

core/solver/gmres.cpp

@@ -107,9 +107,10 @@ void Gmres<ValueType>::apply_impl(const LinOp *b, LinOp *x) const
     auto dense_b = as<const Vector>(b);
     auto dense_x = as<Vector>(x);
     auto residual = Vector::create_with_config_of(dense_b);
-    auto krylov_bases = Vector::create(
-        exec, dim<2>{system_matrix_->get_size()[1] * (krylov_dim_ + 1),
-                     dense_b->get_size()[1]});
+    auto krylov_bases = Vector::create_with_type_of(
+        dense_b, exec,
+        dim<2>{system_matrix_->get_size()[1] * (krylov_dim_ + 1),
+               dense_b->get_size()[1]});
     std::shared_ptr<matrix::Dense<ValueType>> preconditioned_vector =
         Vector::create_with_config_of(dense_b);
     auto hessenberg = Vector::create(

include/ginkgo/core/matrix/dense.hpp

@@ -161,6 +161,41 @@ class Dense
         return (*other).create_with_same_config();
     }
 
+    /**
+     * Creates a Dense matrix with the same type and executor as another Dense
+     * matrix but 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 Dense matrix with the type of other.
+     */
+    static std::unique_ptr<Dense> create_with_type_of(
+        const Dense *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 Dense*, 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<Dense> create_with_type_of(
+        const Dense *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);
+    }
+
     friend class Dense<next_precision<ValueType>>;
 
     void convert_to(Dense<next_precision<ValueType>> *result) const override;
@@ -510,19 +545,7 @@ class Dense
                                             const span &columns,
                                             const size_type stride)
     {
-        row_major_range range_this{this->get_values(), this->get_size()[0],
-                                   this->get_size()[1], this->get_stride()};
-        auto range_result = range_this(rows, columns);
-        // TODO: can result in HUGE padding - which will be copied with the
-        // vector
-        return Dense::create(
-            this->get_executor(),
-            dim<2>{range_result.length(0), range_result.length(1)},
-            Array<ValueType>::view(
-                this->get_executor(),
-                range_result.length(0) * range_this.length(1) - columns.begin,
-                range_result->data),
-            stride);
+        return this->create_submatrix_impl(rows, columns, stride);
     }
 
     /**
@@ -650,6 +673,20 @@ class Dense
                              this->get_stride());
     }
 
+    /**
+     * Creates a Dense matrix with the same type as the callers matrix.
+     *
+     * @param size  size of the matrix
+     *
+     * @returns a Dense matrix with the same type as the caller.
+     */
+    virtual std::unique_ptr<Dense> create_with_type_of_impl(
+        std::shared_ptr<const Executor> exec, const dim<2> &size,
+        size_type stride) const
+    {
+        return Dense::create(exec, size, stride);
+    }
+
     /**
      * @copydoc scale(const LinOp *)
      *
@@ -682,6 +719,32 @@ class Dense
      */
     virtual void compute_norm2_impl(LinOp *result) const;
 
+    /**
+     * @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<Dense> create_submatrix_impl(const span &rows,
+                                                         const span &columns,
+                                                         const size_type stride)
+    {
+        row_major_range range_this{this->get_values(), this->get_size()[0],
+                                   this->get_size()[1], this->get_stride()};
+        auto range_result = range_this(rows, columns);
+        // TODO: can result in HUGE padding - which will be copied with the
+        // vector
+        return Dense::create(
+            this->get_executor(),
+            dim<2>{range_result.length(0), range_result.length(1)},
+            Array<ValueType>::view(
+                this->get_executor(),
+                range_result.length(0) * range_this.length(1) - columns.begin,
+                range_result->data),
+            stride);
+    }
+
     void apply_impl(const LinOp *b, LinOp *x) const override;
 
     void apply_impl(const LinOp *alpha, const LinOp *b, const LinOp *beta,