diff --git a/MathLib/LinAlg/Lis/LisMatrix.cpp b/MathLib/LinAlg/Lis/LisMatrix.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..2e88de19d1606682c20c43c148f12e91ca482405
--- /dev/null
+++ b/MathLib/LinAlg/Lis/LisMatrix.cpp
@@ -0,0 +1,70 @@
+/**
+ * \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 "LisVector.h"
+#include "LisCheck.h"
+
+namespace MathLib
+{
+
+LisMatrix::LisMatrix(std::size_t dimension, LisOption::MatrixType mat_type)
+: _n_rows(dimension), _max_diag_coeff(.0), _mat_type(mat_type)
+{
+ int ierr = 0;
+ ierr = lis_matrix_create(0, &_AA); checkLisError(ierr);
+ ierr = lis_matrix_set_size(_AA, 0, dimension); checkLisError(ierr);
+ lis_matrix_get_range(_AA, &_is, &_ie);
+}
+
+LisMatrix::~LisMatrix()
+{
+ int ierr = 0;
+ ierr = lis_matrix_destroy(_AA); checkLisError(ierr);
+}
+
+void LisMatrix::setZero()
+{
+ int ierr = 0;
+ // A matrix has to be destroyed and created again because Lis doesn't provide a
+ // function to set matrix entries to zero
+ 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);
+
+ _max_diag_coeff = .0;
+}
+
+void LisMatrix::write(const std::string &filename) const
+{
+ lis_output_matrix(_AA, LIS_FMT_MM, const_cast<char*>(filename.c_str()));
+}
+
+void LisMatrix::matvec (const LisVector &x, LisVector &y) const
+{
+ int ierr = lis_matvec(_AA, const_cast<LisVector*>(&x)->getRawVector(), y.getRawVector()); checkLisError(ierr);
+}
+
+bool finalizeMatrixAssembly(LisMatrix &mat)
+{
+ if (!lis_matrix_is_assembled(mat.getRawMatrix())) {
+ int ierr = lis_matrix_set_type(mat.getRawMatrix(), static_cast<int>(mat.getMatrixType())); checkLisError(ierr);
+ ierr = lis_matrix_assemble(mat.getRawMatrix()); checkLisError(ierr);
+ }
+ return true;
+};
+
+
+} //MathLib
diff --git a/MathLib/LinAlg/Lis/LisMatrix.h b/MathLib/LinAlg/Lis/LisMatrix.h
new file mode 100644
index 0000000000000000000000000000000000000000..a0592c13077e911e0f75bb03a5287284e66467aa
--- /dev/null
+++ b/MathLib/LinAlg/Lis/LisMatrix.h
@@ -0,0 +1,139 @@
+/**
+ * \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 Lis matrix wrapper class
+ *
+ * Lis matrix only supports a regular matrix, i.e. the number of
+ * rows should equal to the number of columns.
+ */
+class LisMatrix
+{
+public:
+ /**
+ * constructor
+ * @param length
+ * @param mat_type default 1 CRS
+ */
+ LisMatrix(std::size_t length, 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)
+ {
+ 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);
+ return 0;
+ }
+
+ /// add value
+ int addValue(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);
+ return 0;
+ }
+
+ /// printout this equation for debugging
+ void write(const std::string &filename) const;
+
+ /// get a maximum value in diagonal entries
+ double getMaxDiagCoeff() const { return _max_diag_coeff; };
+
+ /// 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; };
+
+private:
+ std::size_t _n_rows;
+ double _max_diag_coeff;
+ LisOption::MatrixType _mat_type;
+ LIS_MATRIX _AA;
+ int _is;
+ int _ie;
+};
+
+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_
+