Merge sparse direct interface

This adds an interface to the sparse direct solver kernels that were added previously.
It consists of a generic Factorization type storing all factorizations based on triangular matrices,
a LU factorization factory invoking the kernels and a Direct solver type using a Factorization to
dispatch to a pair of triangular solvers.

Related PR: #1082

benchmark/solver/solver.cpp

@@ -78,7 +78,7 @@ DEFINE_string(solvers, "cg",
               "Supported values are: bicgstab, bicg, cb_gmres_keep, "
               "cb_gmres_reduce1, cb_gmres_reduce2, cb_gmres_integer, "
               "cb_gmres_ireduce1, cb_gmres_ireduce2, cg, cgs, fcg, gmres, idr, "
-              "lower_trs, upper_trs, overhead");
+              "lower_trs, upper_trs, symm_direct, overhead");
 
 DEFINE_uint32(
     nrhs, 1,
@@ -312,6 +312,13 @@ std::unique_ptr<gko::LinOpFactory> generate_solver(
         return gko::solver::UpperTrs<etype>::build()
             .with_num_rhs(FLAGS_nrhs)
             .on(exec);
+    } else if (description == "symm_direct") {
+        return gko::experimental::solver::Direct<etype, itype>::build()
+            .with_factorization(
+                gko::experimental::factorization::Lu<etype, itype>::build()
+                    .with_symmetric_sparsity(true)
+                    .on(exec))
+            .on(exec);
     } else if (description == "overhead") {
         return add_criteria_precond_finalize<gko::Overhead<etype>>(
             exec, precond, max_iters);

benchmark/sparse_blas/sparse_blas.cpp

@@ -231,8 +231,8 @@ class SpgemmOperation : public BenchmarkOperation {
 
     gko::size_type get_flops() const override
     {
-        auto host_mtx = Mtx::create(mtx_->get_executor());
-        auto host_mtx2 = Mtx::create(mtx_->get_executor());
+        auto host_mtx = Mtx::create(mtx_->get_executor()->get_master());
+        auto host_mtx2 = Mtx::create(mtx_->get_executor()->get_master());
         host_mtx->copy_from(mtx_);
         host_mtx2->copy_from(mtx2_.get());
         // count the individual products a_ik * b_kj

core/CMakeLists.txt

@@ -16,12 +16,15 @@ target_sources(ginkgo
     base/version.cpp
     distributed/partition.cpp
     factorization/elimination_forest.cpp
+    factorization/factorization.cpp
     factorization/ic.cpp
     factorization/ilu.cpp
+    factorization/lu.cpp
     factorization/par_ic.cpp
     factorization/par_ict.cpp
     factorization/par_ilu.cpp
     factorization/par_ilut.cpp
+    factorization/symbolic.cpp
     log/convergence.cpp
     log/logger.cpp
     log/performance_hint.cpp
@@ -51,6 +54,7 @@ target_sources(ginkgo
     solver/cb_gmres.cpp
     solver/cg.cpp
     solver/cgs.cpp
+    solver/direct.cpp
     solver/fcg.cpp
     solver/gmres.cpp
     solver/idr.cpp

core/factorization/factorization.cpp

@@ -0,0 +1,328 @@
+/*******************************<GINKGO LICENSE>******************************
+Copyright (c) 2017-2022, the Ginkgo authors
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions
+are met:
+
+1. Redistributions of source code must retain the above copyright
+notice, this list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright
+notice, this list of conditions and the following disclaimer in the
+documentation and/or other materials provided with the distribution.
+
+3. Neither the name of the copyright holder nor the names of its
+contributors may be used to endorse or promote products derived from
+this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
+IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
+PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+******************************<GINKGO LICENSE>*******************************/
+
+#include <ginkgo/core/factorization/factorization.hpp>
+
+
+#include <ginkgo/core/base/types.hpp>
+#include <ginkgo/core/matrix/csr.hpp>
+
+
+namespace gko {
+namespace experimental {
+namespace factorization {
+
+
+template <typename ValueType, typename IndexType>
+std::unique_ptr<Factorization<ValueType, IndexType>>
+Factorization<ValueType, IndexType>::unpack() const GKO_NOT_IMPLEMENTED;
+
+
+template <typename ValueType, typename IndexType>
+storage_type Factorization<ValueType, IndexType>::get_storage_type() const
+{
+    return storage_type_;
+}
+
+
+template <typename ValueType, typename IndexType>
+std::shared_ptr<const gko::matrix::Csr<ValueType, IndexType>>
+Factorization<ValueType, IndexType>::get_lower_factor() const
+{
+    switch (storage_type_) {
+    case storage_type::composition:
+    case storage_type::symm_composition:
+        GKO_ASSERT(factors_->get_operators().size() == 2 ||
+                   factors_->get_operators().size() == 3);
+        return as<matrix_type>(factors_->get_operators()[0]);
+    case storage_type::empty:
+    case storage_type::combined_lu:
+    case storage_type::combined_ldu:
+    case storage_type::symm_combined_cholesky:
+    case storage_type::symm_combined_ldl:
+    default:
+        return nullptr;
+    }
+}
+
+
+template <typename ValueType, typename IndexType>
+std::shared_ptr<const gko::matrix::Diagonal<ValueType>>
+Factorization<ValueType, IndexType>::get_diagonal() const
+{
+    switch (storage_type_) {
+    case storage_type::composition:
+    case storage_type::symm_composition:
+        if (factors_->get_operators().size() == 3) {
+            return as<diag_type>(factors_->get_operators()[1]);
+        } else {
+            return nullptr;
+        }
+    case storage_type::empty:
+    case storage_type::combined_lu:
+    case storage_type::combined_ldu:
+    case storage_type::symm_combined_cholesky:
+    case storage_type::symm_combined_ldl:
+    default:
+        return nullptr;
+    }
+}
+
+
+template <typename ValueType, typename IndexType>
+std::shared_ptr<const gko::matrix::Csr<ValueType, IndexType>>
+Factorization<ValueType, IndexType>::get_upper_factor() const
+{
+    switch (storage_type_) {
+    case storage_type::composition:
+    case storage_type::symm_composition:
+        GKO_ASSERT(factors_->get_operators().size() == 2 ||
+                   factors_->get_operators().size() == 3);
+        return as<matrix_type>(factors_->get_operators().back());
+    case storage_type::empty:
+    case storage_type::combined_lu:
+    case storage_type::combined_ldu:
+    case storage_type::symm_combined_cholesky:
+    case storage_type::symm_combined_ldl:
+    default:
+        return nullptr;
+    }
+}
+
+
+template <typename ValueType, typename IndexType>
+std::shared_ptr<const gko::matrix::Csr<ValueType, IndexType>>
+Factorization<ValueType, IndexType>::get_combined() const
+{
+    switch (storage_type_) {
+    case storage_type::combined_lu:
+    case storage_type::combined_ldu:
+    case storage_type::symm_combined_cholesky:
+    case storage_type::symm_combined_ldl:
+        GKO_ASSERT(factors_->get_operators().size() == 1);
+        return as<matrix_type>(factors_->get_operators()[0]);
+    case storage_type::empty:
+    case storage_type::composition:
+    case storage_type::symm_composition:
+    default:
+        return nullptr;
+    }
+}
+
+
+template <typename ValueType, typename IndexType>
+Factorization<ValueType, IndexType>::Factorization(const Factorization& fact)
+    : Factorization{fact.get_executor()}
+{
+    *this = fact;
+}
+
+
+template <typename ValueType, typename IndexType>
+Factorization<ValueType, IndexType>::Factorization(Factorization&& fact)
+    : Factorization{fact.get_executor()}
+{
+    *this = std::move(fact);
+}
+
+
+template <typename ValueType, typename IndexType>
+Factorization<ValueType, IndexType>&
+Factorization<ValueType, IndexType>::operator=(const Factorization& fact)
+{
+    if (this != &fact) {
+        EnableLinOp<Factorization<ValueType, IndexType>>::operator=(fact);
+        storage_type_ = fact.storage_type_;
+        *factors_ = *fact.factors_;
+    }
+    return *this;
+}
+
+
+template <typename ValueType, typename IndexType>
+Factorization<ValueType, IndexType>&
+Factorization<ValueType, IndexType>::operator=(Factorization&& fact)
+{
+    if (this != &fact) {
+        EnableLinOp<Factorization<ValueType, IndexType>>::operator=(
+            std::move(fact));
+        storage_type_ = std::exchange(fact.storage_type_, storage_type::empty);
+        factors_ =
+            std::exchange(fact.factors_, fact.factors_->create_default());
+        if (factors_->get_executor() != this->get_executor()) {
+            factors_ = factors_->clone(this->get_executor());
+        }
+    }
+    return *this;
+}
+
+
+template <typename ValueType, typename IndexType>
+Factorization<ValueType, IndexType>::Factorization(
+    std::shared_ptr<const Executor> exec)
+    : EnableLinOp<Factorization<ValueType, IndexType>>{exec},
+      storage_type_{storage_type::empty},
+      factors_{Composition<ValueType>::create(exec)}
+{}
+
+
+template <typename ValueType, typename IndexType>
+Factorization<ValueType, IndexType>::Factorization(
+    std::unique_ptr<Composition<ValueType>> factors, storage_type type)
+    : EnableLinOp<Factorization<ValueType, IndexType>>{factors->get_executor(),
+                                                       factors->get_size()},
+      storage_type_{type},
+      factors_{std::move(factors)}
+{}
+
+
+template <typename ValueType, typename IndexType>
+std::unique_ptr<Factorization<ValueType, IndexType>>
+Factorization<ValueType, IndexType>::create_from_composition(
+    std::unique_ptr<composition_type> composition)
+{
+    return std::unique_ptr<Factorization<ValueType, IndexType>>{
+        new Factorization<ValueType, IndexType>{std::move(composition),
+                                                storage_type::composition}};
+}
+
+
+template <typename ValueType, typename IndexType>
+std::unique_ptr<Factorization<ValueType, IndexType>>
+Factorization<ValueType, IndexType>::create_from_symm_composition(
+    std::unique_ptr<composition_type> composition)
+{
+    return std::unique_ptr<Factorization<ValueType, IndexType>>{
+        new Factorization<ValueType, IndexType>{
+            std::move(composition), storage_type::symm_composition}};
+}
+
+
+template <typename ValueType, typename IndexType>
+std::unique_ptr<Factorization<ValueType, IndexType>>
+Factorization<ValueType, IndexType>::create_from_combined_lu(
+    std::unique_ptr<matrix_type> combined)
+{
+    return std::unique_ptr<Factorization<ValueType, IndexType>>{
+        new Factorization<ValueType, IndexType>{
+            composition_type::create(gko::share(std::move(combined))),
+            storage_type::combined_lu}};
+}
+
+
+template <typename ValueType, typename IndexType>
+std::unique_ptr<Factorization<ValueType, IndexType>>
+Factorization<ValueType, IndexType>::create_from_combined_ldu(
+    std::unique_ptr<matrix_type> combined)
+{
+    return std::unique_ptr<Factorization<ValueType, IndexType>>{
+        new Factorization<ValueType, IndexType>{
+            composition_type::create(gko::share(std::move(combined))),
+            storage_type::combined_ldu}};
+}
+
+
+template <typename ValueType, typename IndexType>
+std::unique_ptr<Factorization<ValueType, IndexType>>
+Factorization<ValueType, IndexType>::create_from_combined_cholesky(
+    std::unique_ptr<matrix_type> combined)
+{
+    return std::unique_ptr<Factorization<ValueType, IndexType>>{
+        new Factorization<ValueType, IndexType>{
+            composition_type::create(gko::share(std::move(combined))),
+            storage_type::symm_combined_cholesky}};
+}
+
+
+template <typename ValueType, typename IndexType>
+std::unique_ptr<Factorization<ValueType, IndexType>>
+Factorization<ValueType, IndexType>::create_from_combined_ldl(
+    std::unique_ptr<matrix_type> combined)
+{
+    return std::unique_ptr<Factorization<ValueType, IndexType>>{
+        new Factorization<ValueType, IndexType>{
+            composition_type::create(gko::share(std::move(combined))),
+            storage_type::symm_combined_ldl}};
+}
+
+
+template <typename ValueType, typename IndexType>
+void Factorization<ValueType, IndexType>::apply_impl(const LinOp* b,
+                                                     LinOp* x) const
+{
+    switch (storage_type_) {
+    case storage_type::composition:
+    case storage_type::symm_composition:
+        factors_->apply(b, x);
+        break;
+    case storage_type::empty:
+    case storage_type::combined_lu:
+    case storage_type::combined_ldu:
+    case storage_type::symm_combined_cholesky:
+    case storage_type::symm_combined_ldl:
+    default:
+        GKO_NOT_SUPPORTED(storage_type_);
+    }
+}
+
+
+template <typename ValueType, typename IndexType>
+void Factorization<ValueType, IndexType>::apply_impl(const LinOp* alpha,
+                                                     const LinOp* b,
+                                                     const LinOp* beta,
+                                                     LinOp* x) const
+{
+    switch (storage_type_) {
+    case storage_type::composition:
+    case storage_type::symm_composition:
+        factors_->apply(alpha, b, beta, x);
+        break;
+    case storage_type::empty:
+    case storage_type::combined_lu:
+    case storage_type::combined_ldu:
+    case storage_type::symm_combined_cholesky:
+    case storage_type::symm_combined_ldl:
+    default:
+        GKO_NOT_SUPPORTED(storage_type_);
+    }
+}
+
+
+#define GKO_DECLARE_FACTORIZATION(ValueType, IndexType) \
+    class Factorization<ValueType, IndexType>
+
+GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(GKO_DECLARE_FACTORIZATION);
+
+
+}  // namespace factorization
+}  // namespace experimental
+}  // namespace gko