Add algorithm choice for triangular solvers

.gitlab-ci.yml

@@ -252,39 +252,41 @@ build/cuda101/nompi/clang/all/release/static:
     CUDA_ARCH: 35
 
 # clang-cuda with cuda 10.1 and friends
-build/clang-cuda101/openmpi/gcc/cuda/release/shared:
-  extends:
-    - .build_and_test_template
-    - .default_variables
-    - .quick_test_condition
-    - .use_gko-cuda101-openmpi-gnu8-llvm11-intel2019
-  variables:
-    CUDA_COMPILER: "clang++"
-    BUILD_OMP: "ON"
-    BUILD_CUDA: "ON"
-    BUILD_MPI: "ON"
-    MPI_AS_ROOT: "ON"
-    BUILD_HIP: "OFF"
-    BUILD_TYPE: "Release"
-    CUDA_ARCH: 61
-
-
-build/clang-cuda101/nompi/clang/cuda/debug/static:
-  extends:
-    - .build_and_test_template
-    - .default_variables
-    - .full_test_condition
-    - .use_gko-cuda101-openmpi-gnu8-llvm11-intel2019
-  variables:
-    C_COMPILER: "clang"
-    CXX_COMPILER: "clang++"
-    CUDA_COMPILER: "clang++"
-    BUILD_OMP: "ON"
-    BUILD_CUDA: "ON"
-    BUILD_TYPE: "Debug"
-    FAST_TESTS: "ON"
-    BUILD_SHARED_LIBS: "OFF"
-    CUDA_ARCH: 61
+#build/clang-cuda101/openmpi/clang/cuda/release/shared:
+#  extends:
+#    - .build_and_test_template
+#    - .default_variables
+#    - .quick_test_condition
+#    - .use_gko-cuda101-openmpi-gnu8-llvm13-intel2019
+#  variables:
+#    C_COMPILER: "clang"
+#    CXX_COMPILER: "clang++"
+#    CUDA_COMPILER: "clang++"
+#    BUILD_OMP: "ON"
+#    BUILD_CUDA: "ON"
+#    BUILD_MPI: "ON"
+#    MPI_AS_ROOT: "ON"
+#    BUILD_HIP: "OFF"
+#    BUILD_TYPE: "Release"
+#    CUDA_ARCH: 61
+
+
+#build/clang-cuda101/nompi/clang/cuda/debug/static:
+#  extends:
+#    - .build_and_test_template
+#    - .default_variables
+#    - .full_test_condition
+#    - .use_gko-cuda101-openmpi-gnu8-llvm13-intel2019
+#  variables:
+#    C_COMPILER: "clang"
+#    CXX_COMPILER: "clang++"
+#    CUDA_COMPILER: "clang++"
+#    BUILD_OMP: "ON"
+#    BUILD_CUDA: "ON"
+#    BUILD_TYPE: "Debug"
+#    FAST_TESTS: "ON"
+#    BUILD_SHARED_LIBS: "OFF"
+#    CUDA_ARCH: 61
 
 
 # cuda 10.2 and friends
@@ -931,7 +933,7 @@ threadsanitizer:
     - .default_variables
     - .deploy_condition
     - .before_script_template
-    - .use_gko-cuda101-openmpi-gnu8-llvm11-intel2019
+    - .use_gko-cuda101-openmpi-gnu8-llvm13-intel2019
   script:
     - LD_PRELOAD=/usr/local/lib/libomp.so
       CC=clang CXX=clang++
@@ -946,7 +948,7 @@ leaksanitizer:
     - .default_variables
     - .deploy_condition
     - .before_script_template
-    - .use_gko-cuda101-openmpi-gnu8-llvm11-intel2019
+    - .use_gko-cuda101-openmpi-gnu8-llvm13-intel2019
   script:
     - ctest -V -S cmake/CTestScript.cmake -DCTEST_BUILD_CONFIGURATION=LSAN
       -DCTEST_MEMORYCHECK_TYPE=LeakSanitizer
@@ -957,7 +959,7 @@ addresssanitizer:
     - .default_variables
     - .deploy_condition
     - .before_script_template
-    - .use_gko-cuda101-openmpi-gnu8-llvm11-intel2019
+    - .use_gko-cuda101-openmpi-gnu8-llvm13-intel2019
   script:
     - ctest -V -S cmake/CTestScript.cmake -DCTEST_BUILD_CONFIGURATION=ASAN
       -DCTEST_MEMORYCHECK_TYPE=AddressSanitizer
@@ -968,7 +970,7 @@ undefinedsanitizer:
     - .default_variables
     - .deploy_condition
     - .before_script_template
-    - .use_gko-cuda101-openmpi-gnu8-llvm11-intel2019
+    - .use_gko-cuda101-openmpi-gnu8-llvm13-intel2019
   script:
     # the Gold linker is required because of a linker flag issues given by UBsan
     # in the Ubuntu setup we are using.
@@ -981,7 +983,7 @@ cudamemcheck:
     - .before_script_template
     - .default_variables
     - .deploy_condition
-  image: ginkgohub/cuda:101-openmpi-gnu8-llvm11-intel2019
+  image: ginkgohub/cuda:101-openmpi-gnu8-llvm13-intel2019
   tags:
     - private_ci
     - nvidia-gpu

.gitlab/image.yml

@@ -38,8 +38,8 @@
     - controller
     - cpu
 
-.use_gko-cuda101-openmpi-gnu8-llvm11-intel2019:
-  image: ginkgohub/cuda:101-openmpi-gnu8-llvm11-intel2019
+.use_gko-cuda101-openmpi-gnu8-llvm13-intel2019:
+  image: ginkgohub/cuda:101-openmpi-gnu8-llvm13-intel2019
   tags:
     - private_ci
     - nvidia-gpu

core/preconditioner/isai.cpp

@@ -44,8 +44,7 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 #include <ginkgo/core/matrix/csr.hpp>
 #include <ginkgo/core/preconditioner/jacobi.hpp>
 #include <ginkgo/core/solver/gmres.hpp>
-#include <ginkgo/core/solver/lower_trs.hpp>
-#include <ginkgo/core/solver/upper_trs.hpp>
+#include <ginkgo/core/solver/triangular.hpp>
 #include <ginkgo/core/stop/iteration.hpp>
 #include <ginkgo/core/stop/residual_norm.hpp>
 

core/solver/ir.cpp

@@ -89,6 +89,7 @@ Ir<ValueType>& Ir<ValueType>::operator=(const Ir& other)
         EnableLinOp<Ir>::operator=(other);
         EnableSolverBase<Ir>::operator=(other);
         EnableIterativeBase<Ir>::operator=(other);
+        this->parameters_ = other.parameters_;
         this->set_solver(other.get_solver());
         this->set_relaxation_factor(other.relaxation_factor_);
         parameters_ = other.parameters_;
@@ -104,6 +105,7 @@ Ir<ValueType>& Ir<ValueType>::operator=(Ir&& other)
         EnableLinOp<Ir>::operator=(std::move(other));
         EnableSolverBase<Ir>::operator=(std::move(other));
         EnableIterativeBase<Ir>::operator=(std::move(other));
+        this->parameters_ = std::exchange(other.parameters_, parameters_type{});
         this->set_solver(other.get_solver());
         this->set_relaxation_factor(other.relaxation_factor_);
         other.set_solver(nullptr);

core/solver/lower_trs.cpp

@@ -30,9 +30,6 @@ THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 ******************************<GINKGO LICENSE>*******************************/
 
-#include <ginkgo/core/solver/lower_trs.hpp>
-
-
 #include <ginkgo/core/base/array.hpp>
 #include <ginkgo/core/base/exception_helpers.hpp>
 #include <ginkgo/core/base/executor.hpp>
@@ -42,7 +39,7 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 #include <ginkgo/core/base/utils.hpp>
 #include <ginkgo/core/matrix/csr.hpp>
 #include <ginkgo/core/matrix/dense.hpp>
-#include <ginkgo/core/solver/upper_trs.hpp>
+#include <ginkgo/core/solver/triangular.hpp>
 
 
 #include "core/solver/lower_trs_kernels.hpp"
@@ -87,6 +84,7 @@ LowerTrs<ValueType, IndexType>& LowerTrs<ValueType, IndexType>::operator=(
     if (this != &other) {
         EnableLinOp<LowerTrs>::operator=(other);
         EnableSolverBase<LowerTrs, CsrMatrix>::operator=(other);
+        this->parameters_ = other.parameters_;
         this->generate();
     }
     return *this;
@@ -100,6 +98,7 @@ LowerTrs<ValueType, IndexType>& LowerTrs<ValueType, IndexType>::operator=(
     if (this != &other) {
         EnableLinOp<LowerTrs>::operator=(std::move(other));
         EnableSolverBase<LowerTrs, CsrMatrix>::operator=(std::move(other));
+        this->parameters_ = std::exchange(other.parameters_, parameters_type{});
         if (this->get_executor() == other.get_executor()) {
             this->solve_struct_ = std::exchange(other.solve_struct_, nullptr);
         } else {
@@ -136,7 +135,8 @@ void LowerTrs<ValueType, IndexType>::generate()
     if (this->get_system_matrix()) {
         this->get_executor()->run(lower_trs::make_generate(
             this->get_system_matrix().get(), this->solve_struct_,
-            this->get_parameters().unit_diagonal, parameters_.num_rhs));
+            this->get_parameters().unit_diagonal, parameters_.algorithm,
+            parameters_.num_rhs));
     }
 }
 
@@ -178,8 +178,8 @@ void LowerTrs<ValueType, IndexType>::apply_impl(const LinOp* b, LinOp* x) const
             }
             exec->run(lower_trs::make_solve(
                 lend(this->get_system_matrix()), lend(this->solve_struct_),
-                this->get_parameters().unit_diagonal, trans_b, trans_x, dense_b,
-                dense_x));
+                this->get_parameters().unit_diagonal, parameters_.algorithm,
+                trans_b, trans_x, dense_b, dense_x));
         },
         b, x);
 }

core/solver/lower_trs_kernels.hpp

@@ -34,15 +34,13 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 #define GKO_CORE_SOLVER_LOWER_TRS_KERNELS_HPP_
 
 
-#include <ginkgo/core/solver/lower_trs.hpp>
-
-
 #include <memory>
 
 
 #include <ginkgo/core/base/executor.hpp>
 #include <ginkgo/core/matrix/csr.hpp>
 #include <ginkgo/core/matrix/dense.hpp>
+#include <ginkgo/core/solver/triangular.hpp>
 
 
 #include "core/base/kernel_declaration.hpp"
@@ -58,17 +56,19 @@ namespace lower_trs {
                                   bool& do_transpose)
 
 
-#define GKO_DECLARE_LOWER_TRS_GENERATE_KERNEL(_vtype, _itype)         \
-    void generate(std::shared_ptr<const DefaultExecutor> exec,        \
-                  const matrix::Csr<_vtype, _itype>* matrix,          \
-                  std::shared_ptr<solver::SolveStruct>& solve_struct, \
-                  bool unit_diag, const gko::size_type num_rhs)
+#define GKO_DECLARE_LOWER_TRS_GENERATE_KERNEL(_vtype, _itype)                 \
+    void generate(std::shared_ptr<const DefaultExecutor> exec,                \
+                  const matrix::Csr<_vtype, _itype>* matrix,                  \
+                  std::shared_ptr<solver::SolveStruct>& solve_struct,         \
+                  bool unit_diag, const solver::trisolve_algorithm algorithm, \
+                  const size_type num_rhs)
 
 
 #define GKO_DECLARE_LOWER_TRS_SOLVE_KERNEL(_vtype, _itype)                     \
     void solve(std::shared_ptr<const DefaultExecutor> exec,                    \
                const matrix::Csr<_vtype, _itype>* matrix,                      \
                const solver::SolveStruct* solve_struct, bool unit_diag,        \
+               const solver::trisolve_algorithm algorithm,                     \
                matrix::Dense<_vtype>* trans_b, matrix::Dense<_vtype>* trans_x, \
                const matrix::Dense<_vtype>* b, matrix::Dense<_vtype>* x)
 

core/solver/upper_trs.cpp

@@ -30,9 +30,6 @@ THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 ******************************<GINKGO LICENSE>*******************************/
 
-#include <ginkgo/core/solver/upper_trs.hpp>
-
-
 #include <ginkgo/core/base/array.hpp>
 #include <ginkgo/core/base/exception_helpers.hpp>
 #include <ginkgo/core/base/executor.hpp>
@@ -42,7 +39,7 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 #include <ginkgo/core/base/utils.hpp>
 #include <ginkgo/core/matrix/csr.hpp>
 #include <ginkgo/core/matrix/dense.hpp>
-#include <ginkgo/core/solver/lower_trs.hpp>
+#include <ginkgo/core/solver/triangular.hpp>
 
 
 #include "core/solver/upper_trs_kernels.hpp"
@@ -87,6 +84,7 @@ UpperTrs<ValueType, IndexType>& UpperTrs<ValueType, IndexType>::operator=(
     if (this != &other) {
         EnableLinOp<UpperTrs>::operator=(other);
         EnableSolverBase<UpperTrs, CsrMatrix>::operator=(other);
+        this->parameters_ = other.parameters_;
         this->generate();
     }
     return *this;
@@ -100,6 +98,7 @@ UpperTrs<ValueType, IndexType>& UpperTrs<ValueType, IndexType>::operator=(
     if (this != &other) {
         EnableLinOp<UpperTrs>::operator=(std::move(other));
         EnableSolverBase<UpperTrs, CsrMatrix>::operator=(std::move(other));
+        this->parameters_ = std::exchange(other.parameters_, parameters_type{});
         if (this->get_executor() == other.get_executor()) {
             this->solve_struct_ = std::exchange(other.solve_struct_, nullptr);
         } else {
@@ -136,7 +135,8 @@ void UpperTrs<ValueType, IndexType>::generate()
     if (this->get_system_matrix()) {
         this->get_executor()->run(upper_trs::make_generate(
             this->get_system_matrix().get(), this->solve_struct_,
-            this->get_parameters().unit_diagonal, parameters_.num_rhs));
+            this->get_parameters().unit_diagonal, parameters_.algorithm,
+            parameters_.num_rhs));
     }
 }
 
@@ -178,8 +178,8 @@ void UpperTrs<ValueType, IndexType>::apply_impl(const LinOp* b, LinOp* x) const
             }
             exec->run(upper_trs::make_solve(
                 lend(this->get_system_matrix()), lend(this->solve_struct_),
-                this->get_parameters().unit_diagonal, trans_b, trans_x, dense_b,
-                dense_x));
+                this->get_parameters().unit_diagonal, parameters_.algorithm,
+                trans_b, trans_x, dense_b, dense_x));
         },
         b, x);
 }

core/solver/upper_trs_kernels.hpp

@@ -34,15 +34,13 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 #define GKO_CORE_SOLVER_UPPER_TRS_KERNELS_HPP_
 
 
-#include <ginkgo/core/solver/upper_trs.hpp>
-
-
 #include <memory>
 
 
 #include <ginkgo/core/base/executor.hpp>
 #include <ginkgo/core/matrix/csr.hpp>
 #include <ginkgo/core/matrix/dense.hpp>
+#include <ginkgo/core/solver/triangular.hpp>
 
 
 #include "core/base/kernel_declaration.hpp"
@@ -58,17 +56,19 @@ namespace upper_trs {
                                   bool& do_transpose)
 
 
-#define GKO_DECLARE_UPPER_TRS_GENERATE_KERNEL(_vtype, _itype)              \
-    void generate(std::shared_ptr<const DefaultExecutor> exec,             \
-                  const matrix::Csr<_vtype, _itype>* matrix,               \
-                  std::shared_ptr<gko::solver::SolveStruct>& solve_struct, \
-                  bool unit_diag, const gko::size_type num_rhs)
+#define GKO_DECLARE_UPPER_TRS_GENERATE_KERNEL(_vtype, _itype)                 \
+    void generate(std::shared_ptr<const DefaultExecutor> exec,                \
+                  const matrix::Csr<_vtype, _itype>* matrix,                  \
+                  std::shared_ptr<solver::SolveStruct>& solve_struct,         \
+                  bool unit_diag, const solver::trisolve_algorithm algorithm, \
+                  const size_type num_rhs)
 
 
 #define GKO_DECLARE_UPPER_TRS_SOLVE_KERNEL(_vtype, _itype)                     \
     void solve(std::shared_ptr<const DefaultExecutor> exec,                    \
                const matrix::Csr<_vtype, _itype>* matrix,                      \
                const solver::SolveStruct* solve_struct, bool unit_diag,        \
+               const solver::trisolve_algorithm algorithm,                     \
                matrix::Dense<_vtype>* trans_b, matrix::Dense<_vtype>* trans_x, \
                const matrix::Dense<_vtype>* b, matrix::Dense<_vtype>* x)
 

core/test/solver/lower_trs.cpp

@@ -30,16 +30,14 @@ THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 ******************************<GINKGO LICENSE>*******************************/
 
-#include <ginkgo/core/solver/lower_trs.hpp>
-
-
 #include <memory>
 
 
 #include <gtest/gtest.h>
 
 
 #include <ginkgo/core/base/executor.hpp>
+#include <ginkgo/core/solver/triangular.hpp>
 
 
 #include "core/test/utils.hpp"

core/test/solver/upper_trs.cpp

@@ -30,16 +30,14 @@ THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 ******************************<GINKGO LICENSE>*******************************/
 
-#include <ginkgo/core/solver/upper_trs.hpp>
-
-
 #include <memory>
 
 
 #include <gtest/gtest.h>
 
 
 #include <ginkgo/core/base/executor.hpp>
+#include <ginkgo/core/solver/triangular.hpp>
 
 
 #include "core/test/utils.hpp"

cuda/solver/common_trs_kernels.cuh

@@ -315,7 +315,7 @@ void solve_kernel(std::shared_ptr<const CudaExecutor> exec,
                   const matrix::Dense<ValueType>* b,
                   matrix::Dense<ValueType>* x)
 {
-    if (matrix->get_size()[0] == 0) {
+    if (matrix->get_size()[0] == 0 || b->get_size()[1] == 0) {
         return;
     }
     using vec = matrix::Dense<ValueType>;