Iterative solvers via lis (#145)

.ci-ignore

@@ -6,3 +6,4 @@ Toolbox/metis-src
 Toolbox/mumps-src
 Toolbox/superlu-src
 Toolbox/superlumt-src
+Toolbox/lis-src

.github/.codecov.yml

@@ -19,6 +19,7 @@ ignore:
   - "Toolbox/arpack-src"
   - "Toolbox/feast-src"
   - "Toolbox/lapack-ext"
+  - "Toolbox/lis-src"
   - "Toolbox/metis-src"
   - "Toolbox/mumps-src"
   - "Toolbox/superlu-src"

.gitignore

@@ -43,6 +43,7 @@ $RECYCLE.BIN/
 /MSVC/suanPan/superlumt/x64/*
 /MSVC/suanPan/solver.pardiso/x64/*
 /MSVC/suanPan/tester.pardiso/x64/*
+/MSVC/suanPan/lis/x64/*
 venv
 .cache
 

CHANGELOG.md

@@ -6,6 +6,10 @@
 2. The `SPIKE` solver may occasionally hang, to circumvent, try to avoid mixing compilers from different vendors.
 3. The `MUMPS` solver throws segfaults with large matrices.
 
+## version 3.1
+
+1. iterative solvers by the Lis library
+
 ## version 3.0
 
 1. add experimental `MAGMA` based GPU sparse solver [#123](https://github.com/TLCFEM/suanPan/pull/123)

CMakeLists.txt

@@ -76,8 +76,9 @@ add_subdirectory(Solver)
 add_subdirectory(Step)
 add_subdirectory(Toolbox)
 add_subdirectory(UnitTest)
+add_subdirectory(Toolbox/lis-src)
 
-target_link_libraries(${PROJECT_NAME} Element Material Section Solver)
+target_link_libraries(${PROJECT_NAME} Element Material Section Solver lis)
 
 if (FORTRAN_STATUS)
     add_subdirectory(Toolbox/amd-src)

Domain/Factory.hpp

@@ -66,7 +66,8 @@ enum class SolverType {
     CUDA,
     PARDISO,
     FGMRES,
-    MAGMA
+    MAGMA,
+    LIS
 };
 
 template<sp_d T> class Factory final {
@@ -1376,6 +1377,7 @@ template<sp_d T> unique_ptr<MetaMat<T>> Factory<T>::get_basic_container() {
         return std::make_unique<SymmPackMat<T>>(n_size);
     case StorageScheme::SPARSE:
         if(SolverType::MUMPS == solver) return std::make_unique<SparseMatMUMPS<T>>(n_size, n_size, n_elem);
+        if(SolverType::LIS == solver) return std::make_unique<SparseMatLis<T>>(n_size, n_size, n_elem);
         if(SolverType::SUPERLU == solver) return std::make_unique<SparseMatSuperLU<T>>(n_size, n_size, n_elem);
 #ifdef SUANPAN_MKL
         if(SolverType::PARDISO == solver) return std::make_unique<SparseMatPARDISO<T>>(n_size, n_size, n_elem);

Domain/MetaMat/MetaMat

@@ -11,11 +11,12 @@
 #include "FullMatCUDA.hpp"
 #include "SparseMatCUDA.hpp"
 #include "SparseMatFGMRES.hpp"
-#include "SparseMatMUMPS.hpp"
+#include "SparseMatLis.hpp"
+#include "SparseMatMAGMA.hpp"
 #include "SparseMatMPI.hpp"
+#include "SparseMatMUMPS.hpp"
 #include "SparseMatPARDISO.hpp"
 #include "SparseMatSuperLU.hpp"
-#include "SparseMatMAGMA.hpp"
 #include "SymmPackMat.hpp"
 
 // my sparse implementation

Domain/MetaMat/SolverSetting.hpp

@@ -46,6 +46,7 @@ template<sp_d data_t> struct SolverSetting {
     IterativeSolver iterative_solver = IterativeSolver::NONE;
     PreconditionerType preconditioner_type = PreconditionerType::JACOBI;
     Preconditioner<data_t>* preconditioner = nullptr;
+    string lis_options{};
 };
 
 #endif

Domain/MetaMat/SparseMatLis.hpp

@@ -0,0 +1,114 @@
+/*******************************************************************************
+ * Copyright (C) 2017-2023 Theodore Chang
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ ******************************************************************************/
+/**
+ * @class SparseMatLis
+ * @brief A SparseMatLis class that holds matrices.
+ *
+ * @author tlc
+ * @date 15/06/2023
+ * @version 0.1.0
+ * @file SparseMatLis.hpp
+ * @addtogroup MetaMat
+ * @{
+ */
+
+#ifndef SPARSEMATLIS_HPP
+#define SPARSEMATLIS_HPP
+
+#include <lis/lislib.h>
+#include "SparseMat.hpp"
+#include "csr_form.hpp"
+
+template<sp_d T> class SparseMatLis final : public SparseMat<T> {
+    LIS_SOLVER solver = nullptr;
+
+protected:
+    using SparseMat<T>::direct_solve;
+    using SparseMat<T>::setting;
+
+    int direct_solve(Mat<T>&, const Mat<T>&) override;
+
+public:
+    SparseMatLis(const uword in_row, const uword in_col, const uword in_elem = 0)
+        : SparseMat<T>(in_row, in_col, in_elem) {
+        lis_solver_create(&solver);
+        lis_solver_set_option("-i fgmres -p ilu", solver);
+    }
+
+    SparseMatLis(const SparseMatLis& other)
+        : SparseMat<T>(other) {
+        lis_solver_create(&solver);
+        lis_solver_set_option("-i fgmres -p ilu", solver);
+    }
+
+    SparseMatLis(SparseMatLis&&) noexcept = delete;
+    SparseMatLis& operator=(const SparseMatLis&) = delete;
+    SparseMatLis& operator=(SparseMatLis&&) noexcept = delete;
+
+    ~SparseMatLis() override { lis_solver_destroy(solver); }
+
+    unique_ptr<MetaMat<T>> make_copy() override { return std::make_unique<SparseMatLis>(*this); }
+};
+
+template<sp_d T> int SparseMatLis<T>::direct_solve(Mat<T>& X, const Mat<T>& B) {
+    X.set_size(B.n_rows, B.n_cols);
+
+    csr_form<double, LIS_INT> csr_mat(this->triplet_mat, SparseBase::ZERO, true);
+
+    const sp_i auto n = csr_mat.n_rows;
+    const sp_i auto nnz = csr_mat.n_elem;
+
+    LIS_MATRIX A;
+    LIS_VECTOR b, x;
+
+    lis_matrix_create(0, &A);
+    lis_matrix_set_size(A, n, 0);
+    lis_matrix_set_csr(nnz, csr_mat.row_mem(), csr_mat.col_mem(), csr_mat.val_mem(), A);
+    lis_matrix_assemble(A);
+
+    lis_vector_create(0, &b);
+    lis_vector_create(0, &x);
+    lis_vector_set_size(b, n, 0);
+    lis_vector_set_size(x, n, 0);
+
+    lis_solver_set_option(setting.lis_options.c_str(), solver);
+
+    for(uword I = 0; I < B.n_cols; ++I) {
+        // ReSharper disable CppCStyleCast
+        lis_vector_set(b, (double*)B.colptr(I));
+        lis_vector_set(x, (double*)X.colptr(I));
+        // ReSharper restore CppCStyleCast
+
+        lis_solve(A, b, x, solver);
+    }
+
+    A->ptr = nullptr;
+    A->index = nullptr;
+    A->value = nullptr;
+    b->value = nullptr;
+    x->value = nullptr;
+
+    lis_matrix_destroy(A);
+    lis_vector_destroy(b);
+    lis_vector_destroy(x);
+
+    return SUANPAN_SUCCESS;
+}
+
+#endif
+
+//! @}

Example/Solver/lis.supan

@@ -0,0 +1,45 @@
+node 1 0 0
+node 2 .5 0
+node 3 .5 .5
+node 4 0 .5
+
+material Elastic2D 1 12 .1 1E-4
+
+element CP3 1 1 2 3 1 1
+element CP3 2 1 3 4 1 1
+
+fix2 1 1 1 4
+fix2 2 2 1
+
+cload 1 0 1 2 2 3
+
+step static 1
+set ini_step_size 1E-1
+set fixed_step_size 1
+set sparse_mat 1
+set system_solver lis -print 2 -i 17 -p ilu -ilu_fill 1
+
+converger RelIncreDisp 1 1E-8 20 1
+
+analyze
+
+# Node 2:
+# Coordinate:
+# 0.5000 0
+# Displacement:
+# 0.1498 0.6417
+# Resistance:
+# -1.1102e-16 1.0000e+00
+# 
+# Node 3:
+# Coordinate:
+# 0.5000 0.5000
+# Displacement:
+# -0.1665 0.6081
+# Resistance:
+# 0 1.0000
+peek node 2 3
+
+reset
+clear
+exit