diff --git a/MathLib/CMakeLists.txt b/MathLib/CMakeLists.txt
index 4864062b78211ff9d83acd126ff5aed0089bde25..43096a270168f4b30a87d12516208117f8f16056 100644
--- a/MathLib/CMakeLists.txt
+++ b/MathLib/CMakeLists.txt
@@ -20,16 +20,10 @@ SET ( SOURCES ${SOURCES} ${SOURCES_LINALG_SOLVERS})
GET_SOURCE_FILES(SOURCES_LINALG_PRECOND LinAlg/Preconditioner)
SET ( SOURCES ${SOURCES} ${SOURCES_LINALG_PRECOND})
-GET_SOURCE_FILES(SOURCES_LINALG_LEQS LinAlg/SystemOfLinearEquations)
-IF (NOT OGS_USE_LIS)
- LIST(REMOVE_ITEM SOURCES_LINALG_LEQS
- LinAlg/SystemOfLinearEquations/LisOption.h
- LinAlg/SystemOfLinearEquations/LisOption.cpp
- LinAlg/SystemOfLinearEquations/LisLinearSystem.h
- LinAlg/SystemOfLinearEquations/LisLinearSystem.cpp
- )
+IF (OGS_USE_LIS)
+ GET_SOURCE_FILES(SOURCES_LINALG_LIS LinAlg/Lis)
+ SET ( SOURCES ${SOURCES} ${SOURCES_LINALG_LIS})
ENDIF ()
-SET ( SOURCES ${SOURCES} ${SOURCES_LINALG_LEQS})
IF (METIS_FOUND)
GET_SOURCE_FILES(SOURCES_LINALG_SPARSE_NESTEDDISSECTION LinAlg/Sparse/NestedDissectionPermutation)
diff --git a/MathLib/LinAlg/Lis/LisCheck.h b/MathLib/LinAlg/Lis/LisCheck.h
new file mode 100644
index 0000000000000000000000000000000000000000..558f4746966c2aa0be43338aa89f1ce691ab3e64
--- /dev/null
+++ b/MathLib/LinAlg/Lis/LisCheck.h
@@ -0,0 +1,43 @@
+/**
+ * \file
+ * \author Norihiro Watanabe
+ * \date 2013-04-16
+ * \brief
+ *
+ * \copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef LISCHECK_H_
+#define LISCHECK_H_
+
+#include <vector>
+#include "lis.h"
+#include "logog/include/logog.hpp"
+
+namespace MathLib
+{
+
+/**
+ * check Lis error codes
+ *
+ * @param err Lis error code
+ * @return success or not
+ */
+inline bool checkLisError(int err)
+{
+ bool ok = (err == LIS_SUCCESS);
+ if (!ok) {
+ ERR("***ERROR: Lis error code = %d", err);
+ }
+ return ok;
+}
+
+} // MathLib
+
+#endif //LISCHECK_H_
+
diff --git a/MathLib/LinAlg/Lis/LisLinearSolver.cpp b/MathLib/LinAlg/Lis/LisLinearSolver.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6f6e44c355468655de6de59286cd4a5cdcdb1752
--- /dev/null
+++ b/MathLib/LinAlg/Lis/LisLinearSolver.cpp
@@ -0,0 +1,133 @@
+/**
+ * \file
+ * \author Norihiro Watanabe
+ * \date 2013-04-16
+ * \brief Implementation of the LisLinearSolver class.
+ *
+ * \copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "LisLinearSolver.h"
+
+#include <string>
+#include <sstream>
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+#include "logog/include/logog.hpp"
+
+#include "LisCheck.h"
+
+namespace MathLib
+{
+
+using boost::property_tree::ptree;
+
+LisLinearSolver::LisLinearSolver(LisMatrix &A, ptree const*const option)
+: _A(A)
+{
+ if (option)
+ setOption(*option);
+}
+
+void LisLinearSolver::setOption(const ptree &option)
+{
+ boost::optional<ptree> ptSolver = option.get_child("LinearSolver");
+ if (!ptSolver)
+ return;
+
+ boost::optional<std::string> solver_type = ptSolver->get_optional<std::string>("solver_type");
+ if (solver_type) {
+ _option.solver_type = _option.getSolverType(*solver_type);
+ }
+ boost::optional<std::string> precon_type = ptSolver->get_optional<std::string>("precon_type");
+ if (precon_type) {
+ _option.precon_type = _option.getPreconType(*precon_type);
+ }
+ boost::optional<std::string> matrix_type = ptSolver->get_optional<std::string>("matrix_type");
+ if (matrix_type) {
+ _option.matrix_type = _option.getMatrixType(*matrix_type);
+ }
+ boost::optional<double> error_tolerance = ptSolver->get_optional<double>("error_tolerance");
+ if (error_tolerance) {
+ _option.error_tolerance = *error_tolerance;
+ }
+ boost::optional<int> max_iteration_step = ptSolver->get_optional<int>("max_iteration_step");
+ if (max_iteration_step) {
+ _option.max_iterations = *max_iteration_step;
+ }
+}
+
+void LisLinearSolver::solve(LisVector &b, LisVector &x)
+{
+ finalizeMatrixAssembly(_A);
+
+ INFO("------------------------------------------------------------------");
+ INFO("*** LIS solver computation");
+#ifdef _OPENMP
+ INFO("-> max number of threads = %d", omp_get_num_procs());
+ INFO("-> number of threads = %d", omp_get_max_threads());
+#endif
+
+ // configure option
+ std::string solver_options;
+ if (_option.solver_precon_arg.empty()) {
+ std::stringstream ss;
+ ss << "-i " << static_cast<int>(_option.solver_type);
+ ss << " -p " << static_cast<int>(_option.precon_type);
+ if (!_option.extra_arg.empty())
+ ss << " " << _option.extra_arg;
+ solver_options = ss.str();
+ } else {
+ solver_options = _option.solver_precon_arg;
+ }
+ std::string tol_option;
+ {
+ std::stringstream ss;
+ ss << "-tol " << _option.error_tolerance;
+ ss << " -maxiter " << _option.max_iterations;
+ ss << " -initx_zeros 0"; //0: use given x as initial guess, 1: x0=0
+#ifdef _OPENMP
+ const int nthreads = omp_get_max_threads();
+ ss << " -omp_num_threads " << nthreads;
+#endif
+ tol_option = ss.str();
+ }
+
+ // Create solver
+ LIS_SOLVER solver;
+ int ierr = lis_solver_create(&solver);
+ checkLisError(ierr);
+ ierr = lis_solver_set_option(const_cast<char*>(solver_options.c_str()), solver);
+ checkLisError(ierr);
+ ierr = lis_solver_set_option(const_cast<char*>(tol_option.c_str()), solver);
+ checkLisError(ierr);
+ ierr = lis_solver_set_option(const_cast<char*>("-print mem"), solver);
+ checkLisError(ierr);
+
+ // solve
+ INFO("-> solve");
+ ierr = lis_solve(_A.getRawMatrix(), b.getRawVector(), x.getRawVector(), solver);
+ checkLisError(ierr);
+
+ int iter = 0;
+ double resid = 0.0;
+ ierr = lis_solver_get_iters(solver, &iter);
+ checkLisError(ierr);
+ ierr = lis_solver_get_residualnorm(solver, &resid);
+ checkLisError(ierr);
+ INFO("\t iteration: %d/%ld\n", iter, _option.max_iterations);
+ INFO("\t residual: %e\n", resid);
+
+ // Clear solver
+ ierr = lis_solver_destroy(solver);
+ checkLisError(ierr);
+ INFO("------------------------------------------------------------------");
+}
+
+} //MathLib
diff --git a/MathLib/LinAlg/Lis/LisLinearSolver.h b/MathLib/LinAlg/Lis/LisLinearSolver.h
new file mode 100644
index 0000000000000000000000000000000000000000..4bcd496b58bb9a8ed48391a5f9dbe8410649a59c
--- /dev/null
+++ b/MathLib/LinAlg/Lis/LisLinearSolver.h
@@ -0,0 +1,83 @@
+/**
+ * \file
+ * \author Norihiro Watanabe
+ * \date 2013-04-16
+ * \brief Definition of the LisLinearSolver class.
+ *
+ * \copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef LISLINEARSOLVER_H_
+#define LISLINEARSOLVER_H_
+
+#include <vector>
+#include <boost/property_tree/ptree.hpp>
+
+#include "lis.h"
+
+#include "LisOption.h"
+#include "LisVector.h"
+#include "LisMatrix.h"
+
+namespace MathLib
+{
+
+/**
+ * \brief Linear solver using Lis (http://www.ssisc.org/lis/)
+ *
+ */
+class LisLinearSolver
+{
+public:
+ /**
+ * Constructor
+ * @param A Coefficient matrix object
+ * @param option A pointer to a linear solver option. In case you omit
+ * this argument, default settings follow those of
+ * LisOption struct.
+ */
+ LisLinearSolver(LisMatrix &A, boost::property_tree::ptree const*const option = nullptr);
+
+ virtual ~LisLinearSolver() {};
+
+ /**
+ * configure linear solvers
+ * @param option
+ */
+ void setOption(const boost::property_tree::ptree &option);
+
+ /**
+ * configure linear solvers
+ * @param option
+ */
+ void setOption(const LisOption &option) { _option = option; }
+
+ /**
+ * get linear solver options
+ * @return
+ */
+ LisOption &getOption() { return _option; }
+
+ /**
+ * solve a given linear equations
+ *
+ * @param b RHS vector
+ * @param x Solution vector
+ */
+ void solve(LisVector &b, LisVector &x);
+
+
+private:
+ LisMatrix& _A;
+ LisOption _option;
+};
+
+} // MathLib
+
+#endif //LISLINEARSOLVER_H_
+
diff --git a/MathLib/LinAlg/Lis/LisMatrix.cpp b/MathLib/LinAlg/Lis/LisMatrix.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..8aae74c5ccbb5dbcaaf72def19b17f5bfacfb18e
--- /dev/null
+++ b/MathLib/LinAlg/Lis/LisMatrix.cpp
@@ -0,0 +1,128 @@
+/**
+ * \file
+ * \author Norihiro Watanabe
+ * \date 2013-04-16
+ * \brief Implementation of the LisMatrix class.
+ *
+ * \copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "LisMatrix.h"
+
+#include <stdexcept>
+
+#include "LisVector.h"
+#include "LisCheck.h"
+
+namespace MathLib
+{
+
+LisMatrix::LisMatrix(std::size_t n_rows, LisOption::MatrixType mat_type)
+ : _n_rows(n_rows), _mat_type(mat_type), _is_assembled(false)
+{
+ int ierr = lis_matrix_create(0, &_AA);
+ checkLisError(ierr);
+ ierr = lis_matrix_set_size(_AA, 0, n_rows);
+ checkLisError(ierr);
+ lis_matrix_get_range(_AA, &_is, &_ie);
+}
+
+LisMatrix::~LisMatrix()
+{
+ int ierr = lis_matrix_destroy(_AA);
+ checkLisError(ierr);
+}
+
+void LisMatrix::setZero()
+{
+ // A matrix has to be destroyed and created again because Lis doesn't provide a
+ // function to set matrix entries to zero
+ int ierr = lis_matrix_destroy(_AA);
+ checkLisError(ierr);
+ ierr = lis_matrix_create(0, &_AA);
+ checkLisError(ierr);
+ ierr = lis_matrix_set_size(_AA, 0, _n_rows);
+ checkLisError(ierr);
+
+ _is_assembled = false;
+}
+
+int LisMatrix::setValue(std::size_t rowId, std::size_t colId, double v)
+{
+ lis_matrix_set_value(LIS_INS_VALUE, rowId, colId, v, _AA);
+ _is_assembled = false;
+ return 0;
+}
+
+int LisMatrix::addValue(std::size_t rowId, std::size_t colId, double v)
+{
+ lis_matrix_set_value(LIS_ADD_VALUE, rowId, colId, v, _AA);
+ _is_assembled = false;
+ return 0;
+}
+
+void LisMatrix::write(const std::string &filename) const
+{
+ if (!_is_assembled)
+ throw std::logic_error("LisMatrix::write(): matrix not assembled.");
+ lis_output_matrix(_AA, LIS_FMT_MM, const_cast<char*>(filename.c_str()));
+}
+
+double LisMatrix::getMaxDiagCoeff()
+{
+ if (!_is_assembled)
+ throw std::logic_error("LisMatrix::getMaxDiagCoeff(): matrix not assembled.");
+ LIS_VECTOR diag;
+ int ierr = lis_vector_create(0, &diag);
+ checkLisError(ierr);
+ ierr = lis_vector_set_size(diag, 0, _n_rows);
+ checkLisError(ierr);
+
+ ierr = lis_matrix_get_diagonal(_AA, diag);
+ checkLisError(ierr);
+
+ double abs_max_entry;
+ ierr = lis_vector_get_value(diag, 0, &abs_max_entry);
+ checkLisError(ierr);
+ abs_max_entry = std::abs(abs_max_entry);
+ for (std::size_t k(1); k<_n_rows; ++k) {
+ double tmp;
+ ierr = lis_vector_get_value(diag, k, &tmp);
+ checkLisError(ierr);
+ if (abs_max_entry < std::abs(tmp)) {
+ abs_max_entry = std::abs(tmp);
+ }
+ }
+
+ return abs_max_entry;
+}
+
+void LisMatrix::matvec (const LisVector &x, LisVector &y) const
+{
+ if (!_is_assembled)
+ throw std::logic_error("LisMatrix::matvec(): matrix not assembled.");
+ int ierr = lis_matvec(_AA, const_cast<LisVector*>(&x)->getRawVector(), y.getRawVector());
+ checkLisError(ierr);
+}
+
+bool finalizeMatrixAssembly(LisMatrix &mat)
+{
+ LIS_MATRIX &A = mat.getRawMatrix();
+
+ if (!mat.isAssembled()) {
+ int ierr = lis_matrix_set_type(A, static_cast<int>(mat.getMatrixType()));
+ checkLisError(ierr);
+ ierr = lis_matrix_assemble(A);
+ checkLisError(ierr);
+ mat._is_assembled = true;
+ }
+ return true;
+}
+
+
+} //MathLib
diff --git a/MathLib/LinAlg/Lis/LisMatrix.h b/MathLib/LinAlg/Lis/LisMatrix.h
new file mode 100644
index 0000000000000000000000000000000000000000..b6a61662ed02f53967f0d3201bf1083947e5d6b4
--- /dev/null
+++ b/MathLib/LinAlg/Lis/LisMatrix.h
@@ -0,0 +1,134 @@
+/**
+ * \file
+ * \author Norihiro Watanabe
+ * \date 2013-04-16
+ * \brief Definition of the LisMatrix class.
+ *
+ * \copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef LISMATRIX_H_
+#define LISMATRIX_H_
+
+#include <iostream>
+#include <cmath>
+#include <vector>
+
+#include "lis.h"
+
+#include "LisOption.h"
+
+namespace MathLib
+{
+class LisVector;
+
+/**
+ * \brief LisMatrix is a wrapper class for matrix types of the
+ * linear iterative solvers library.
+ *
+ * LisMatrix only supports square matrices, i.e. the number of
+ * rows have to be equal to the number of columns.
+ */
+class LisMatrix
+{
+public:
+ /**
+ * constructor
+ * @param n_rows the number of rows (that is equal to the number of columns)
+ * @param mat_type default 1 CRS
+ */
+ LisMatrix(std::size_t n_rows, LisOption::MatrixType mat_type = LisOption::MatrixType::CRS);
+
+ /**
+ *
+ */
+ virtual ~LisMatrix();
+
+ /// return the number of rows
+ std::size_t getNRows() const { return _n_rows; };
+
+ /// return the number of columns
+ std::size_t getNCols() const { return getNRows(); };
+
+ /// return a start index of the active data range
+ std::size_t getRangeBegin() const { return _is; }
+
+ /// return an end index of the active data range
+ std::size_t getRangeEnd() const { return _ie; }
+
+ /// reset this matrix with keeping its original dimension
+ void setZero();
+
+ /// set entry
+ int setValue(std::size_t rowId, std::size_t colId, double v);
+
+ /// add value
+ int addValue(std::size_t rowId, std::size_t colId, double v);
+
+ /// printout this equation for debugging
+ void write(const std::string &filename) const;
+
+ /// get a maximum value in diagonal entries
+ double getMaxDiagCoeff();
+
+ /// return a raw Lis matrix object
+ LIS_MATRIX& getRawMatrix() { return _AA; };
+
+ /// y = mat * x
+ void matvec(const LisVector &x, LisVector &y) const;
+
+ ///
+ template <class T_DENSE_MATRIX>
+ void addSubMatrix(std::vector<std::size_t> const& row_pos,
+ std::vector<std::size_t> const& col_pos, const T_DENSE_MATRIX &sub_matrix,
+ double fkt = 1.0);
+
+ /// get this matrix type
+ LisOption::MatrixType getMatrixType() const { return _mat_type; };
+
+ /// return if this matrix is already assembled or not
+ bool isAssembled() const { return _is_assembled; };
+
+private:
+ std::size_t const _n_rows;
+ LisOption::MatrixType const _mat_type;
+ LIS_MATRIX _AA;
+ bool _is_assembled;
+ int _is; ///< location where the partial matrix _AA starts in global matrix.
+ int _ie; ///< location where the partial matrix _AA ends in global matrix.
+
+ // friend function
+ friend bool finalizeMatrixAssembly(LisMatrix &mat);
+};
+
+template<class T_DENSE_MATRIX>
+void
+LisMatrix::addSubMatrix(std::vector<std::size_t> const& row_pos, std::vector<std::size_t> const& col_pos,
+ const T_DENSE_MATRIX &sub_matrix, double fkt)
+{
+ if (row_pos.size() != sub_matrix.getNRows() || col_pos.size() != sub_matrix.getNCols())
+ return;
+
+ const std::size_t n_rows = row_pos.size();
+ const std::size_t n_cols = col_pos.size();
+ for (std::size_t i = 0; i < n_rows; i++) {
+ const std::size_t row = row_pos[i];
+ for (std::size_t j = 0; j < n_cols; j++) {
+ const std::size_t col = col_pos[j];
+ addValue(row, col, fkt * sub_matrix(i, j));
+ }
+ }
+};
+
+/// finish assembly to make this matrix be ready for use
+bool finalizeMatrixAssembly(LisMatrix &mat);
+
+} // MathLib
+
+#endif //LISMATRIX_H_
+
diff --git a/MathLib/LinAlg/SystemOfLinearEquations/LisOption.cpp b/MathLib/LinAlg/Lis/LisOption.cpp
similarity index 100%
rename from MathLib/LinAlg/SystemOfLinearEquations/LisOption.cpp
rename to MathLib/LinAlg/Lis/LisOption.cpp
diff --git a/MathLib/LinAlg/SystemOfLinearEquations/LisOption.h b/MathLib/LinAlg/Lis/LisOption.h
similarity index 100%
rename from MathLib/LinAlg/SystemOfLinearEquations/LisOption.h
rename to MathLib/LinAlg/Lis/LisOption.h
diff --git a/MathLib/LinAlg/Lis/LisTools.cpp b/MathLib/LinAlg/Lis/LisTools.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c0549b5fb355d6221d2ae70b3a2f8eebcc93297e
--- /dev/null
+++ b/MathLib/LinAlg/Lis/LisTools.cpp
@@ -0,0 +1,47 @@
+/**
+ * \file
+ * \author Norihiro Watanabe
+ * \date 2013-04-16
+ * \brief Implementation of Lis utility functions.
+ *
+ * \copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "LisTools.h"
+
+#include "logog/include/logog.hpp"
+
+#include "LisMatrix.h"
+#include "LisVector.h"
+
+namespace MathLib
+{
+
+void applyKnownSolution(LisMatrix &A, LisVector &b, const std::vector<std::size_t> &vec_knownX_id,
+ const std::vector<double> &vec_knownX_x, double penalty_scaling)
+{
+ //Use penalty parameter
+ const double max_diag_coeff = A.getMaxDiagCoeff();
+ const double penalty = max_diag_coeff * penalty_scaling;
+ INFO("-> max. absolute value of diagonal entries = %e", max_diag_coeff);
+ INFO("-> penalty scaling = %e", penalty_scaling);
+ const std::size_t n_bc = vec_knownX_id.size();
+ for (std::size_t i_bc=0; i_bc<n_bc; i_bc++) {
+ const std::size_t rowId = vec_knownX_id[i_bc];
+ const double x = vec_knownX_x[i_bc];
+ //A(k, k) = penalty
+ A.setValue(rowId, rowId, penalty);
+ //b(k) = x*penalty
+ b.set(rowId, x*penalty);
+ }
+}
+
+} // MathLib
+
+
+
diff --git a/MathLib/LinAlg/Lis/LisTools.h b/MathLib/LinAlg/Lis/LisTools.h
new file mode 100644
index 0000000000000000000000000000000000000000..c5c542f8871c728c0562aad55c549b362391d500
--- /dev/null
+++ b/MathLib/LinAlg/Lis/LisTools.h
@@ -0,0 +1,43 @@
+/**
+ * \file
+ * \author Norihiro Watanabe
+ * \date 2013-04-16
+ * \brief Definition of Lis utility functions.
+ *
+ * \copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef LISTOOLS_H_
+#define LISTOOLS_H_
+
+#include <vector>
+
+namespace MathLib
+{
+class LisMatrix;
+class LisVector;
+
+/**
+ * apply known solutions to a system of linear equations
+ *
+ * This function introduces the constants into the system by the penalty method.
+ *
+ * @param A Coefficient matrix
+ * @param b RHS vector
+ * @param vec_knownX_id a vector of known solution entry IDs
+ * @param vec_knownX_x a vector of known solutions
+ * @param penalty_scaling value for scaling some matrix and right hand side
+ * entries to enforce some conditions
+ */
+void applyKnownSolution(LisMatrix &A, LisVector &b, const std::vector<std::size_t> &_vec_knownX_id,
+ const std::vector<double> &_vec_knownX_x, double penalty_scaling = 1e+10);
+
+} // MathLib
+
+#endif //LISTOOLS_H_
+
diff --git a/MathLib/LinAlg/Lis/LisVector.cpp b/MathLib/LinAlg/Lis/LisVector.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7965f654679b99054e52fa2048d423da016c5acf
--- /dev/null
+++ b/MathLib/LinAlg/Lis/LisVector.cpp
@@ -0,0 +1,68 @@
+/**
+ * \file
+ * \author Norihiro Watanabe
+ * \date 2013-04-16
+ * \brief Implementation of the LisVector class.
+ *
+ * \copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "LisVector.h"
+
+namespace MathLib
+{
+
+LisVector::LisVector(std::size_t length)
+: _length(length)
+{
+ lis_vector_create(0, &_vec);
+ lis_vector_set_size(_vec, 0, length);
+}
+
+LisVector::LisVector(LisVector const &src)
+: _length(src.size())
+{
+ lis_vector_duplicate(src._vec, &_vec);
+ lis_vector_copy(src._vec, _vec);
+}
+
+LisVector::~LisVector()
+{
+ lis_vector_destroy(_vec);
+}
+
+LisVector& LisVector::operator= (const LisVector &src)
+{
+ lis_vector_copy(src._vec, _vec);
+ return *this;
+}
+
+void LisVector::operator+= (const LisVector& v)
+{
+ lis_vector_axpy(1.0, v._vec, _vec);
+}
+
+void LisVector::operator-= (const LisVector& v)
+{
+ lis_vector_axpy(-1.0, v._vec, _vec);
+}
+
+LisVector& LisVector::operator= (double v)
+{
+ lis_vector_set_all(v, _vec);
+ return *this;
+}
+
+void LisVector::write (const std::string &filename) const
+{
+ lis_output_vector(_vec, LIS_FMT_PLAIN, const_cast<char*>(filename.c_str()));
+}
+
+
+} // MathLib
+
diff --git a/MathLib/LinAlg/Lis/LisVector.h b/MathLib/LinAlg/Lis/LisVector.h
new file mode 100644
index 0000000000000000000000000000000000000000..b8a258cebb696c0e2883f7fcdbf755b485cba4b2
--- /dev/null
+++ b/MathLib/LinAlg/Lis/LisVector.h
@@ -0,0 +1,113 @@
+/**
+ * \file
+ * \author Norihiro Watanabe
+ * \date 2013-04-16
+ * \brief Definition of the LisVector class.
+ *
+ * \copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef LISVECTOR_H_
+#define LISVECTOR_H_
+
+#include <iostream>
+#include <vector>
+
+#include "lis.h"
+
+namespace MathLib
+{
+
+/**
+ * \brief Lis vector wrapper class
+ */
+class LisVector
+{
+public:
+ /**
+ * Constructor for initialization of the number of rows
+ * @param length number of rows
+ * @return
+ */
+ explicit LisVector(std::size_t length);
+
+ /// copy constructor
+ LisVector(LisVector const &src);
+
+ /**
+ *
+ */
+ virtual ~LisVector();
+
+ /// return a vector length
+ std::size_t size() const {return _length;};
+
+ /// return a start index of the active data range
+ std::size_t getRangeBegin() const { return 0;}
+
+ /// return an end index of the active data range
+ std::size_t getRangeEnd() const { return _length; }
+
+ /// set all values in this vector
+ LisVector& operator= (double v);
+
+ /// access entry
+ double operator[] (std::size_t rowId) const { return get(rowId); }
+
+ /// get entry
+ double get(std::size_t rowId) const
+ {
+ double v = .0;
+ lis_vector_get_value(_vec, rowId, &v);
+ return v;
+ }
+
+ /// set entry
+ void set(std::size_t rowId, double v)
+ {
+ lis_vector_set_value(LIS_INS_VALUE, rowId, v, _vec);
+ }
+
+ /// add entry
+ void add(std::size_t rowId, double v)
+ {
+ lis_vector_set_value(LIS_ADD_VALUE, rowId, v, _vec);
+ }
+
+ /// printout this equation for debugging
+ void write (const std::string &filename) const;
+
+ /// return a raw Lis vector object
+ LIS_VECTOR& getRawVector() {return _vec; };
+
+ /// vector operation: set data
+ LisVector& operator= (const LisVector &src);
+
+ /// vector operation: add
+ void operator+= (const LisVector& v);
+
+ /// vector operation: subtract
+ void operator-= (const LisVector& v);
+
+ ///
+ template<class T_SUBVEC>
+ void addSubVector(const std::vector<std::size_t> &pos, const T_SUBVEC &sub_vec)
+ {
+ for (std::size_t i=0; i<pos.size(); ++i) {
+ this->add(pos[i], sub_vec[i]);
+ }
+ }
+private:
+ std::size_t const _length;
+ LIS_VECTOR _vec;
+};
+
+} // MathLib
+
+#endif //LISVECTOR_H_
+
diff --git a/MathLib/LinAlg/SystemOfLinearEquations/ISystemOfLinearEquations.h b/MathLib/LinAlg/SystemOfLinearEquations/ISystemOfLinearEquations.h
deleted file mode 100644
index a47d05cb3372f46e7b67d98e565472569f244d9e..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/SystemOfLinearEquations/ISystemOfLinearEquations.h
+++ /dev/null
@@ -1,204 +0,0 @@
-/**
- * \file
- * \author Norihiro Watanabe
- * \date 2012-06-25
- * \brief Definition of the ISystemOfLinearEquations class.
- *
- * \copyright
- * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#ifndef ISYSTEMOFLINEAREQUATIONS_H_
-#define ISYSTEMOFLINEAREQUATIONS_H_
-
-#include <iostream>
-#include <vector>
-#include <boost/property_tree/ptree.hpp>
-#include "MathLib/LinAlg/Sparse/Sparsity.h"
-
-namespace MathLib
-{
-
-/**
- * \brief Interface to any system of linear equations Ax=b
- *
- * This class defines unified interface to a system of linear equations
- * including linear solvers.
- *
- */
-class ISystemOfLinearEquations
-{
-public:
- /// Index representing invalid
- const static std::size_t index_npos = -1;
-
- ///
- virtual ~ISystemOfLinearEquations()
- {
- }
-
- /**
- * set properties
- *
- * @param option
- */
- virtual void setOption(const boost::property_tree::ptree &option) = 0;
-
- /**
- * initialize this system with zero
- *
- * This function initialize the system by setting zero to all entries in
- * a matrix, a RHS vector and solution vector. Dimension of the system is
- * not changed.
- */
- virtual void setZero() = 0;
-
- /// return dimension of this equation
- virtual std::size_t getDimension() const = 0;
-
- /**
- * set an entry in a matrix
- *
- * @param rowId
- * @param colId
- * @param v
- */
- virtual void setMatEntry(std::size_t rowId, std::size_t colId, double v) = 0;
-
- /**
- * add a value into a matrix
- *
- * @param rowId
- * @param colId
- * @param v
- */
- virtual void addMatEntry(std::size_t rowId, std::size_t colId, double v) = 0;
-
- /**
- * add a sub-matrix into a matrix at given row and column positions
- *
- * @param vec_row_pos
- * A vector of global row index. The number of entries of vec_row_pos have
- * to be identical to the number of rows of the local matrix.
- * @param vec_col_pos
- * A vector of global column index. The number of vec_col_pos have to be
- * identical to the number of columns of the local matrix.
- * @param sub_matrix A sub-matrix
- * @param fac scaling parameter
- */
- template <class T_DENSE_MATRIX>
- void addSubMat(const std::vector<std::size_t> &vec_row_pos,
- const std::vector<std::size_t> &vec_col_pos,
- const T_DENSE_MATRIX &sub_matrix, double fac = 1.0);
-
- /**
- * add a sub-matrix into a matrix
- *
- * @param vec_row_col_pos
- * A vector of global row and column index. The number of entries of
- * vec_row_col_pos have to be identical to the number of rows and columns
- * of the local matrix.
- * @param sub_matrix A sub-matrix
- * @param fac scaling parameter
- */
- template <class T_DENSE_MATRIX>
- void addSubMat(const std::vector<std::size_t> &vec_row_col_pos,
- const T_DENSE_MATRIX &sub_matrix, double fac = 1.0);
-
- /**
- * get RHS entry
- *
- * @param rowId
- * @return
- */
- virtual double getRHSVec(std::size_t rowId) const = 0;
-
- /**
- * set RHS entry
- *
- * @param rowId
- * @param v
- */
- virtual void setRHSVec(std::size_t rowId, double v) = 0;
-
- /**
- * add RHS entry
- *
- * @param rowId
- * @param v
- */
- virtual void addRHSVec(std::size_t rowId, double v) = 0;
-
- /**
- * add a sub vector to RHS
- *
- * @param vec_row_pos
- * A vector of global row index. The number of entries of vec_row_pos
- * have to be identical to the number of rows of the sub vector.
- * @param sub_vector
- * Pointer to a sub-vector. Its length must be the same as the length of vec_row_pos.
- * @param fac
- * Scaling parameter. Default value is 1.
- */
- inline virtual void addSubRHS(const std::vector<std::size_t> &vec_row_pos,
- const double* sub_vector, double fac = 1.0);
-
- /**
- * add a sub vector to RHS
- *
- * @param vec_row_pos
- * A vector of global row index. The number of entries of vec_row_pos
- * have to be identical to the number of rows of the sub vector.
- * @param sub_vector A sub-vector
- * @param fac Scaling parameter
- */
- template <class T_DENSE_VECTOR>
- void addSubRHS(const std::vector<std::size_t> &vec_row_pos,
- const T_DENSE_VECTOR &sub_vector, double fac = 1.0);
-
- /**
- * get an entry in a solution vector
- * @param rowId
- * @return
- */
- virtual double getSolVec(std::size_t rowId) = 0;
-
- /**
- * set an entry in a solution vector
- * @param rowId
- * @param v
- */
- virtual void setSolVec(std::size_t rowId, double v) = 0;
-
- /**
- * set prescribed value
- * @param row_id
- * @param x
- */
- virtual void setKnownSolution(std::size_t row_id, double x) = 0;
-
- /**
- * set prescribed value
- * @param vec_id A vector of global row index
- * @param vec_x A vector of prescribed values
- */
- virtual void setKnownSolution(const std::vector<std::size_t> &vec_id,
- const std::vector<double> &vec_x) = 0;
-
- /// solve the linear system
- virtual void solve() = 0;
-
- /// print out the equation for debugging
- virtual void printout(std::ostream &os = std::cout) const = 0;
-
-};
-
-}
-
-#include "ISystemOfLinearEquations.tpp"
-
-#endif //ISYSTEMOFLINEAREQUATIONS_H_
diff --git a/MathLib/LinAlg/SystemOfLinearEquations/ISystemOfLinearEquations.tpp b/MathLib/LinAlg/SystemOfLinearEquations/ISystemOfLinearEquations.tpp
deleted file mode 100644
index aea3ee5961e79f5021abdb7e81675381f4999c72..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/SystemOfLinearEquations/ISystemOfLinearEquations.tpp
+++ /dev/null
@@ -1,65 +0,0 @@
-/**
- * \file
- * \author Norihiro Watanabe
- * \date 2012-06-25
- * \brief Implementation of the ISystemOfLinearEquations class.
- *
- * \copyright
- * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#ifndef ISYSTEMOFLINEAREQUATIONS_TPP_
-#define ISYSTEMOFLINEAREQUATIONS_TPP_
-
-namespace MathLib
-{
-
-template <class T_DENSE_MATRIX>
-void ISystemOfLinearEquations::addSubMat(const std::vector<std::size_t> &vec_row_pos, const std::vector<std::size_t> &vec_col_pos, const T_DENSE_MATRIX &sub_matrix, double fkt)
-{
- const std::size_t n_rows = vec_row_pos.size();
- const std::size_t n_cols = vec_col_pos.size();
- for (std::size_t i=0; i<n_rows; i++) {
- const std::size_t rowId = vec_row_pos[i];
- if (rowId==index_npos) continue;
- for (std::size_t j=0; j<n_cols; j++) {
- const std::size_t colId = vec_col_pos[j];
- if (colId==index_npos) continue;
- addMatEntry(rowId, colId, fkt*sub_matrix(i,j));
- }
- }
-}
-
-template <class T_DENSE_MATRIX>
-void ISystemOfLinearEquations::addSubMat(const std::vector<std::size_t> &vec_pos, const T_DENSE_MATRIX &sub_matrix, double fkt)
-{
- addSubMat(vec_pos, vec_pos, sub_matrix, fkt);
-}
-
-inline void ISystemOfLinearEquations::addSubRHS(const std::vector<std::size_t> &vec_row_pos, const double *sub_vector, double fkt)
-{
- for (std::size_t i=0; i<vec_row_pos.size(); i++) {
- const std::size_t rowId = vec_row_pos[i];
- if (rowId==index_npos) continue;
- addRHSVec(rowId, sub_vector[i]*fkt);
- }
-}
-
-template <class T_DENSE_VECTOR>
-void ISystemOfLinearEquations::addSubRHS(const std::vector<std::size_t> &vec_row_pos, const T_DENSE_VECTOR &sub_vector, double fkt)
-{
- for (std::size_t i=0; i<vec_row_pos.size(); i++) {
- const std::size_t rowId = vec_row_pos[i];
- if (rowId==index_npos) continue;
- addRHSVec(rowId, sub_vector[i]*fkt);
- }
-}
-
-}
-
-#endif //ISYSTEMOFLINEAREQUATIONS_TPP_
-
diff --git a/MathLib/LinAlg/SystemOfLinearEquations/LisLinearSystem.cpp b/MathLib/LinAlg/SystemOfLinearEquations/LisLinearSystem.cpp
deleted file mode 100644
index 97e764abd5561f752d4c0dfd3b8066dac67171e1..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/SystemOfLinearEquations/LisLinearSystem.cpp
+++ /dev/null
@@ -1,191 +0,0 @@
-/**
- * \file
- * \author Norihiro Watanabe
- * \date 2012-06-25
- * \brief Implementation of the LisLinearSystem class.
- *
- * \copyright
- * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#include "LisLinearSystem.h"
-
-#include <string>
-#include <sstream>
-#ifdef _OPENMP
-#include <omp.h>
-#endif
-#include "logog/include/logog.hpp"
-
-namespace MathLib
-{
-
-bool LisLinearSystem::checkError(int err)
-{
- bool ok = (err == LIS_SUCCESS);
- if (!ok) {
- ERR("***ERROR: Lis error code = %d", err);
- }
- return ok;
-}
-
-LisLinearSystem::LisLinearSystem(std::size_t dimension)
-: _dim(dimension), _max_diag_coeff(.0)
-{
- int ierr = 0;
- ierr = lis_matrix_create(0, &_AA); checkError(ierr);
- ierr = lis_matrix_set_size(_AA, 0, dimension); checkError(ierr);
- ierr = lis_vector_duplicate(_AA, &_bb); checkError(ierr);
- ierr = lis_vector_duplicate(_AA, &_xx); checkError(ierr);
-}
-
-LisLinearSystem::~LisLinearSystem()
-{
- int ierr = 0;
- ierr = lis_matrix_destroy(_AA); checkError(ierr);
- ierr = lis_vector_destroy(_bb); checkError(ierr);
- ierr = lis_vector_destroy(_xx); checkError(ierr);
-}
-
-void LisLinearSystem::setOption(const boost::property_tree::ptree &option)
-{
- using boost::property_tree::ptree;
-
- boost::optional<ptree> ptSolver = option.get_child("LinearSolver");
- if (!ptSolver)
- return;
-
- boost::optional<std::string> solver_type = ptSolver->get_optional<std::string>("solver_type");
- if (solver_type) {
- _option.solver_type = _option.getSolverType(*solver_type);
- }
- boost::optional<std::string> precon_type = ptSolver->get_optional<std::string>("precon_type");
- if (precon_type) {
- _option.precon_type = _option.getPreconType(*precon_type);
- }
- boost::optional<std::string> matrix_type = ptSolver->get_optional<std::string>("matrix_type");
- if (matrix_type) {
- _option.matrix_type = _option.getMatrixType(*matrix_type);
- }
- boost::optional<double> error_tolerance = ptSolver->get_optional<double>("error_tolerance");
- if (error_tolerance) {
- _option.error_tolerance = *error_tolerance;
- }
- boost::optional<int> max_iteration_step = ptSolver->get_optional<int>("max_iteration_step");
- if (max_iteration_step) {
- _option.max_iterations = *max_iteration_step;
- }
-}
-
-void LisLinearSystem::setZero()
-{
- int ierr = 0;
- // A matrix has to be removed and created every time because Lis doesn't provide a
- // function to set matrix entries to zero
- ierr = lis_matrix_destroy(_AA); checkError(ierr);
- ierr = lis_matrix_create(0, &_AA); checkError(ierr);
- ierr = lis_matrix_set_size(_AA, 0, getDimension()); checkError(ierr);
- // set zero RHS, x
- ierr = lis_vector_set_all(0.0, _bb); checkError(ierr);
- ierr = lis_vector_set_all(0.0, _xx); checkError(ierr);
-
- _max_diag_coeff = .0;
-}
-
-void LisLinearSystem::applyKnownSolution()
-{
- //Use penalty parameter
- const double penalty_scaling = 1e+10;
- const double penalty = _max_diag_coeff * penalty_scaling;
- INFO("-> max. absolute value of diagonal entries = %e", _max_diag_coeff);
- INFO("-> penalty scaling = %e", penalty_scaling);
- const std::size_t n_bc = _vec_knownX_id.size();
- for (std::size_t i_bc=0; i_bc<n_bc; i_bc++) {
- const std::size_t rowId = _vec_knownX_id[i_bc];
- const double x = _vec_knownX_x[i_bc];
- //A(k, k) = penalty
- setMatEntry(rowId, rowId, penalty);
- //b(k) = x*penalty
- setRHSVec(rowId, x*penalty);
- }
-}
-
-void LisLinearSystem::solve()
-{
- INFO("------------------------------------------------------------------");
- INFO("*** LIS solver computation");
-#ifdef _OPENMP
- INFO("-> max number of threads = %d", omp_get_num_procs());
- INFO("-> number of threads = %d", omp_get_max_threads());
-#endif
-
- applyKnownSolution();
- // assemble a matrix
- int ierr = 0;
- ierr = lis_matrix_set_type(_AA, static_cast<int>(_option.matrix_type)); checkError(ierr);
- ierr = lis_matrix_assemble(_AA); checkError(ierr);
-
- // configure option
- std::string solver_options;
- if (_option.solver_precon_arg.empty()) {
- std::stringstream ss;
- ss << "-i " << static_cast<int>(_option.solver_type);
- ss << " -p " << static_cast<int>(_option.precon_type);
- if (!_option.extra_arg.empty())
- ss << " " << _option.extra_arg;
- solver_options = ss.str();
- } else {
- solver_options = _option.solver_precon_arg;
- }
- std::string tol_option;
- {
- std::stringstream ss;
- ss << "-tol " << _option.error_tolerance;
- ss << " -maxiter " << _option.max_iterations;
- ss << " -initx_zeros 0"; //0: use given x as initial guess, 1: x0=0
-#ifdef _OPENMP
- const int nthreads = omp_get_max_threads();
- ss << " -omp_num_threads " << nthreads;
-#endif
- tol_option = ss.str();
- }
-
- // Create solver
- LIS_SOLVER solver;
- ierr = lis_solver_create(&solver); checkError(ierr);
- ierr = lis_solver_set_option(const_cast<char*>(solver_options.c_str()), solver); checkError(ierr);
- ierr = lis_solver_set_option(const_cast<char*>(tol_option.c_str()), solver); checkError(ierr);
- ierr = lis_solver_set_option(const_cast<char*>("-print mem"), solver); checkError(ierr);
-
- // solve
- INFO("-> solve");
- ierr = lis_solve(_AA, _bb, _xx, solver); checkError(ierr);
- checkError(ierr);
- //lis_output(_AA, _bb, _xx, LIS_FMT_MM, "/home/localadmin/tasks/20120814_ogs6test/matrix1.txt");
-
- int iter = 0;
- double resid = 0.0;
- ierr = lis_solver_get_iters(solver, &iter); checkError(ierr);
- ierr = lis_solver_get_residualnorm(solver, &resid); checkError(ierr);
- INFO("\t iteration: %d/%ld\n", iter, _option.max_iterations);
- INFO("\t residuals: %e\n", resid);
-
- // Clear solver
- ierr = lis_solver_destroy(solver); checkError(ierr);
- INFO("------------------------------------------------------------------");
-}
-
-void LisLinearSystem::printout(std::ostream &os) const
-{
- os << "#A=" << std::endl;
- os << "#x=" << std::endl;
- lis_vector_print(_xx);
- os << "#b=" << std::endl;
- lis_vector_print(_bb);
-}
-
-} //MathLib
diff --git a/MathLib/LinAlg/SystemOfLinearEquations/LisLinearSystem.h b/MathLib/LinAlg/SystemOfLinearEquations/LisLinearSystem.h
deleted file mode 100644
index 8ccf30e0d2b0a7ef339286b1eb80ca52f7ed5677..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/SystemOfLinearEquations/LisLinearSystem.h
+++ /dev/null
@@ -1,185 +0,0 @@
-/**
- * \file
- * \author Norihiro Watanabe
- * \date 2012-06-25
- * \brief Definition of the LisLinearSystem class.
- *
- * \copyright
- * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#ifndef LISLINEARSYSTEM_H_
-#define LISLINEARSYSTEM_H_
-
-#include <iostream>
-#include <string>
-#include <vector>
-#include <boost/property_tree/ptree.hpp>
-
-#include "lis.h"
-
-#include "LisOption.h"
-
-namespace MathLib
-{
-
-/**
- * \brief Linear system using Lis solver (http://www.ssisc.org/lis/)
- *
- * This class utilizes Lis solver.
- */
-class LisLinearSystem
-{
-public:
- //---------------------------------------------------------------
- // realization of ISystemOfLinearEquations
- //---------------------------------------------------------------
-
- /**
- * Create a linear system using Lis
- *
- * @param length System dimension
- */
- LisLinearSystem(std::size_t length);
-
- /**
- *
- */
- virtual ~LisLinearSystem();
-
- /**
- * configure linear solvers
- * @param option
- */
- virtual void setOption(const boost::property_tree::ptree &option);
-
- /// return the system dimension
- virtual std::size_t getDimension() const { return _dim; };
-
- /// reset this equation
- virtual void setZero();
-
- /// set entry in A
- virtual void setMatEntry(std::size_t rowId, std::size_t colId, double v)
- {
- if (rowId==colId)
- _max_diag_coeff = std::max(_max_diag_coeff, std::abs(v));
- lis_matrix_set_value(LIS_INS_VALUE, rowId, colId, v, _AA);
- }
-
- /// add value into A
- virtual void addMatEntry(std::size_t rowId, std::size_t colId, double v)
- {
- if (rowId==colId)
- _max_diag_coeff = std::max(_max_diag_coeff, std::abs(v));
- lis_matrix_set_value(LIS_ADD_VALUE, rowId, colId, v, _AA);
- }
-
- /// get RHS entry
- virtual double getRHSVec(std::size_t rowId) const
- {
- double v;
- lis_vector_get_value(_bb, rowId, &v);
- return v;
- }
-
- /// set RHS entry
- virtual void setRHSVec(std::size_t rowId, double v)
- {
- lis_vector_set_value(LIS_INS_VALUE, rowId, v, _bb);
- }
-
- /// add RHS entry
- virtual void addRHSVec(std::size_t rowId, double v)
- {
- lis_vector_set_value(LIS_ADD_VALUE, rowId, v, _bb);
- }
-
- /// get an entry in a solution vector
- virtual double getSolVec(std::size_t rowId)
- {
- double v;
- lis_vector_get_value(_xx, rowId, &v);
- return v;
- }
-
- /// set a solution vector
- virtual void setSolVec(std::size_t rowId, double v)
- {
- lis_vector_set_value(LIS_INS_VALUE, rowId, v, _xx);
- }
-
- /// set prescribed value
- virtual void setKnownSolution(std::size_t rowId, double x)
- {
- _vec_knownX_id.push_back(rowId);
- _vec_knownX_x.push_back(x);
- }
-
- /// set prescribed values
- virtual void setKnownSolution(const std::vector<std::size_t> &vec_id, const std::vector<double> &vec_x)
- {
- _vec_knownX_id.insert(_vec_knownX_id.end(), vec_id.begin(), vec_id.end());
- _vec_knownX_x.insert(_vec_knownX_x.end(), vec_x.begin(), vec_x.end());
- }
-
- /// solve this equation
- virtual void solve();
-
- /// printout this equation for debugging
- virtual void printout(std::ostream &os=std::cout) const;
-
- //---------------------------------------------------------------
- // specific to this class
- //---------------------------------------------------------------
- /**
- * configure linear solvers
- * @param option
- */
- void setOption(const LisOption &option)
- {
- _option = option;
- }
-
- /**
- * get linear solver options
- * @return
- */
- LisOption &getOption()
- {
- return _option;
- }
-
- /**
- * get a solution vector
- * @param x
- */
- void getSolVec(double* x)
- {
- lis_vector_gather(_xx, x);
- }
-
-private:
- bool checkError(int err);
- void applyKnownSolution();
-
-private:
- const std::size_t _dim;
- double _max_diag_coeff;
- LisOption _option;
- LIS_MATRIX _AA;
- LIS_VECTOR _bb;
- LIS_VECTOR _xx;
- std::vector<std::size_t> _vec_knownX_id;
- std::vector<double> _vec_knownX_x;
-};
-
-
-} // MathLib
-
-#endif //LISLINEARSYSTEM_H_
-
diff --git a/Tests/MathLib/TestLinearSystem.cpp b/Tests/MathLib/TestLinearSystem.cpp
deleted file mode 100644
index a0d61685000cd9e6674ed0e4c26c3cd7a7148ff7..0000000000000000000000000000000000000000
--- a/Tests/MathLib/TestLinearSystem.cpp
+++ /dev/null
@@ -1,129 +0,0 @@
-/**
- * \file
- * \author Norihiro Watanabe
- * \date 2012-08-03
- * \brief Implementation tests.
- *
- * \copyright
- * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#include <gtest/gtest.h>
-#include <boost/property_tree/ptree.hpp>
-
-#include "MathLib/LinAlg/Sparse/Sparsity.h"
-#ifdef USE_LIS
-#include "MathLib/LinAlg/SystemOfLinearEquations/LisLinearSystem.h"
-#endif
-
-namespace
-{
-
-inline void ASSERT_DOUBLE_ARRAY_EQ(const double* Expected, const double* Actual, size_t N, double epsilon=1.0e-8) {
- for (size_t i=0; i<N; i++) \
- ASSERT_NEAR(Expected[i], Actual[i], epsilon);
-}
-
-struct Example1
-{
- std::vector<double> mat;
- std::vector<size_t> list_dirichlet_bc_id;
- std::vector<double> list_dirichlet_bc_value;
- static const size_t dim_eqs = 9;
- std::vector<double> exH;
- MathLib::RowMajorSparsity sparse;
-
- Example1()
- {
- double d_mat[] = {
- 6.66667e-012, -1.66667e-012, 0, -1.66667e-012, -3.33333e-012, 0, 0, 0, 0,
- -1.66667e-012, 1.33333e-011, -1.66667e-012, -3.33333e-012, -3.33333e-012, -3.33333e-012, 0, 0, 0,
- 0, -1.66667e-012, 6.66667e-012, 0, -3.33333e-012, -1.66667e-012, 0, 0, 0,
- -1.66667e-012, -3.33333e-012, 0, 1.33333e-011, -3.33333e-012, 0, -1.66667e-012, -3.33333e-012, 0,
- -3.33333e-012, -3.33333e-012, -3.33333e-012, -3.33333e-012, 2.66667e-011, -3.33333e-012, -3.33333e-012, -3.33333e-012, -3.33333e-012,
- 0, -3.33333e-012, -1.66667e-012, 0, -3.33333e-012, 1.33333e-011, 0, -3.33333e-012, -1.66667e-012,
- 0, 0, 0, -1.66667e-012, -3.33333e-012, 0, 6.66667e-012, -1.66667e-012, 0,
- 0, 0, 0, -3.33333e-012, -3.33333e-012, -3.33333e-012, -1.66667e-012, 1.33333e-011, -1.66667e-012,
- 0, 0, 0, 0, -3.33333e-012, -1.66667e-012, 0, -1.66667e-012, 6.66667e-012
- };
- mat.assign(d_mat, d_mat+dim_eqs*dim_eqs);
- size_t int_dirichlet_bc_id[] = {2,5,8,0,3,6};
- list_dirichlet_bc_id.assign(int_dirichlet_bc_id, int_dirichlet_bc_id+6);
- list_dirichlet_bc_value.resize(6);
- fill(list_dirichlet_bc_value.begin(), list_dirichlet_bc_value.begin()+3, .0);
- fill(list_dirichlet_bc_value.begin()+3, list_dirichlet_bc_value.end(), 1.0);
- exH.resize(9);
- for (size_t i=0; i<9; i++) {
- if (i%3==0) exH[i] = 1.0;
- if (i%3==1) exH[i] = 0.5;
- if (i%3==2) exH[i] = 0.;
- }
- sparse.resize(dim_eqs);
- for (size_t i=0; i<dim_eqs; i++) {
- for (size_t j=0; j<dim_eqs; j++) {
- if (mat[i*dim_eqs+j]!=.0)
- sparse[i].insert(j);
- }
- }
-
- }
-};
-}
-
-#ifdef USE_LIS
-TEST(Math, LinearSystemLis)
-{
- // set a problem
- Example1 ex1;
-
- // create a linear system
- MathLib::LisLinearSystem eqs(ex1.dim_eqs);
-
- // construct
- for (size_t i=0; i<ex1.dim_eqs; i++) {
- for (size_t j=0; j<ex1.dim_eqs; j++) {
- double v = ex1.mat[i*ex1.dim_eqs+j];
- if (v!=.0)
- eqs.addMatEntry(i, j, v);
- }
- }
-
- // apply BC
- eqs.setKnownSolution(ex1.list_dirichlet_bc_id, ex1.list_dirichlet_bc_value);
-
- // set solver options using Boost property tree
- boost::property_tree::ptree t_root;
- boost::property_tree::ptree t_solver;
- t_solver.put("solver_type", "CG");
- t_solver.put("precon_type", "NONE");
- t_solver.put("matrix_type", "CCS");
- t_solver.put("error_tolerance", 1e-15);
- t_solver.put("max_iteration_step", 1000);
- t_root.put_child("LinearSolver", t_solver);
- eqs.setOption(t_root);
-
- // check if the option was correctly parsed
- MathLib::LisOption &lisOption = eqs.getOption();
- ASSERT_EQ(MathLib::LisOption::SolverType::CG, lisOption.solver_type);
- ASSERT_EQ(MathLib::LisOption::PreconType::NONE, lisOption.precon_type);
- ASSERT_EQ(MathLib::LisOption::MatrixType::CCS, lisOption.matrix_type);
- ASSERT_EQ(1e-15, lisOption.error_tolerance);
- ASSERT_EQ(1000, lisOption.max_iterations);
-
- // solve
- eqs.solve();
-
-// eqs.printout();
-
- // check solution
- std::vector<double> vec_x(eqs.getDimension());
- eqs.getSolVec(&vec_x[0]);
- ASSERT_DOUBLE_ARRAY_EQ(&ex1.exH[0], &vec_x[0], ex1.dim_eqs, 1.e-5);
-}
-
-#endif
-