From 21bbdec20f80a7bde2229cb80a1956df5b09eff4 Mon Sep 17 00:00:00 2001
From: acrovato <39187559+acrovato@users.noreply.github.com>
Date: Tue, 1 Mar 2022 14:58:45 +0100
Subject: [PATCH] Add interface for Intel DSS
---
tbox/_src/tboxw.h | 1 +
tbox/_src/tboxw.i | 2 +
tbox/src/tbox.h | 1 +
tbox/src/wDss.cpp | 100 ++++++++++++++++++++++++++++++++++++++++++++++
tbox/src/wDss.h | 61 ++++++++++++++++++++++++++++
5 files changed, 165 insertions(+)
create mode 100644 tbox/src/wDss.cpp
create mode 100644 tbox/src/wDss.h
diff --git a/tbox/_src/tboxw.h b/tbox/_src/tboxw.h
index b9feacd..3da0640 100644
--- a/tbox/_src/tboxw.h
+++ b/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"
diff --git a/tbox/_src/tboxw.i b/tbox/_src/tboxw.i
index 36e9d5a..94acd7f 100644
--- a/tbox/_src/tboxw.i
+++ b/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);
diff --git a/tbox/src/tbox.h b/tbox/src/tbox.h
index b6814d4..d136502 100644
--- a/tbox/src/tbox.h
+++ b/tbox/src/tbox.h
@@ -99,6 +99,7 @@ class Linesearch;
// Linear solvers
class LinearSolver;
class Pardiso;
+class Dss;
class Mumps;
class SparseLu;
class Gmres;
diff --git a/tbox/src/wDss.cpp b/tbox/src/wDss.cpp
new file mode 100644
index 0000000..c412f3c
--- /dev/null
+++ b/tbox/src/wDss.cpp
@@ -0,0 +1,100 @@
+/*
+ * 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
diff --git a/tbox/src/wDss.h b/tbox/src/wDss.h
new file mode 100644
index 0000000..a7e3d2d
--- /dev/null
+++ b/tbox/src/wDss.h
@@ -0,0 +1,61 @@
+/*
+ * 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
--
GitLab