Add interface for Intel DSS

tbox/_src/tboxw.h

@@ -25,6 +25,7 @@
 #include "wCacheQuad4.h"
 #include "wCacheTetra4.h"
 #include "wCacheTri3.h"
+#include "wDss.h"
 #include "wElement.h"
 #include "wFct0.h"
 #include "wFct1.h"

tbox/_src/tboxw.i

@@ -156,6 +156,8 @@ namespace std {
 %include "wLinearSolver.h"
 %shared_ptr(tbox::Pardiso);
 %include "wPardiso.h"
+%shared_ptr(tbox::Dss);
+%include "wDss.h"
 %shared_ptr(tbox::Mumps);
 %include "wMumps.h"
 %shared_ptr(tbox::SparseLu);

tbox/src/tbox.h

@@ -99,6 +99,7 @@ class Linesearch;
 // Linear solvers
 class LinearSolver;
 class Pardiso;
+class Dss;
 class Mumps;
 class SparseLu;
 class Gmres;

tbox/src/wDss.cpp

+/*
+ * Copyright 2020 University of Liège
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#include "wDss.h"
+#ifdef USE_MKL
+
+using namespace tbox;
+
+Dss::Dss() : LinearSolver()
+{
+    handle = NULL;
+    MKL_INT opt = MKL_DSS_MSG_LVL_INFO + MKL_DSS_TERM_LVL_FATAL + MKL_DSS_ZERO_BASED_INDEXING; // continue exec until non-recoverable (fatal) error
+    print(dss_create(handle, opt));
+}
+Dss::~Dss()
+{
+    if (handle)
+    {
+        MKL_INT opt = MKL_DSS_DEFAULTS;
+        print(dss_delete(handle, opt));
+        handle = NULL;
+    }
+}
+
+void Dss::analyze(Eigen::SparseMatrix<double> const &A)
+{
+    // Check that matrix A is square
+    if (A.rows() != A.cols())
+        throw std::runtime_error("tbox::Dss: Matrix A is not square!\n");
+    // Convert to CRS
+    Eigen::SparseMatrix<double, Eigen::RowMajor> AA = A;
+    if (!AA.isCompressed())
+        AA.makeCompressed();
+
+    // Build data structure
+    MKL_INT opt = MKL_DSS_NON_SYMMETRIC;        // consider MKL_DSS_SYMMETRIC and MKL_DSS_SYMMETRIC_STRUCTURE
+    MKL_INT n = AA.rows();                      // matrix size (number of rows and columns)
+    MKL_INT nz = AA.nonZeros();                 // number of non-zero entries
+    const MKL_INT *rStart = AA.outerIndexPtr(); // starting index of non-zero entry for a row
+    const MKL_INT *cIndex = AA.innerIndexPtr(); // index of column of non-zero entry
+    print(dss_define_structure(handle, opt, rStart, n, n, cIndex, nz));
+
+    // Define reordering
+    opt = MKL_DSS_AUTO_ORDER;
+    print(dss_reorder(handle, opt, 0));
+}
+void Dss::factorize(Eigen::SparseMatrix<double> const &A)
+{
+    Eigen::SparseMatrix<double, Eigen::RowMajor> AA = A;
+    MKL_INT opt = MKL_DSS_INDEFINITE;
+    const double *val = AA.valuePtr();
+    print(dss_factor_real(handle, opt, val));
+}
+void Dss::solve(Eigen::Map<Eigen::VectorXd> const &b, Eigen::Map<Eigen::VectorXd> &x)
+{
+    MKL_INT opt = MKL_DSS_DEFAULTS;
+    MKL_INT n = 1; // number of rhs
+    const double *rhs = b.data();
+    double *sol = x.data();
+    print(dss_solve_real(handle, opt, rhs, n, sol));
+}
+void Dss::compute(Eigen::SparseMatrix<double> const &A, Eigen::Map<Eigen::VectorXd> const &b, Eigen::Map<Eigen::VectorXd> &x)
+{
+    this->analyze(A);
+    this->factorize(A);
+    this->solve(b, x);
+}
+
+/**
+ * @brief Print status code
+ */
+void Dss::print(int cod)
+{
+    if (cod != MKL_DSS_SUCCESS)
+    {
+        std::stringstream err;
+        err << "tbox::Dss: DSS returned code " << cod << "\n";
+        throw std::runtime_error(err.str());
+    }
+}
+
+void Dss::write(std::ostream &out) const
+{
+    out << "DSS (MKL)" << std::endl;
+}
+
+#endif // USE_MKL

tbox/src/wDss.h

+/*
+ * Copyright 2020 University of Liège
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#ifndef WDSS_H
+#define WDSS_H
+
+#include "tbox.h"
+#ifdef USE_MKL
+
+#include "wLinearSolver.h"
+#include <mkl_dss.h>
+#include <Eigen/Sparse>
+#include <limits>
+
+namespace tbox
+{
+/**
+ * @brief Intel DSS (PARDISO) interface
+ * @authors Adrien Crovato
+ */
+class TBOX_API Dss : public LinearSolver
+{
+    void *handle; // DSS handle
+
+    void print(int cod);
+
+public:
+    Dss();
+    ~Dss();
+
+#ifndef SWIG
+    virtual void analyze(Eigen::SparseMatrix<double> const &A) override;
+    virtual void factorize(Eigen::SparseMatrix<double> const &A) override;
+    virtual void solve(Eigen::Map<Eigen::VectorXd> const &b, Eigen::Map<Eigen::VectorXd> &x) override;
+    virtual void compute(Eigen::SparseMatrix<double> const &A, Eigen::Map<Eigen::VectorXd> const &b, Eigen::Map<Eigen::VectorXd> &x) override;
+
+    virtual double getError() override { return std::numeric_limits<double>::epsilon(); }
+    virtual int getIterations() override { return 1; }
+
+    virtual void write(std::ostream &out) const override;
+#endif
+};
+
+} // namespace tbox
+
+#endif // USE_MKL
+
+#endif // WDSS_H
\ No newline at end of file