diff --git a/CMakeLists.txt b/CMakeLists.txt
index afc39bb68171c04e64c70010878655b13e48aca1..0d04b975529508cfa6eedc5112ed09744402faae 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -266,7 +266,6 @@ if( OGS_BUILD_TESTS AND NOT IS_SUBPROJECT )
     if(OGS_USE_MPI)
         add_subdirectory( SimpleTests/MeshTests/MPI )
     else()
-        add_subdirectory( SimpleTests/MatrixTests )
         add_subdirectory( SimpleTests/MeshTests )
     endif()
 endif() # OGS_BUILD_TESTS
diff --git a/GeoLib/AnalyticalGeometry.cpp b/GeoLib/AnalyticalGeometry.cpp
index eea3e423d7563efb4c94899b4c234bd82b1a0b90..b15c1d374cdba993be231fd2e93931a7ce4e55ef 100644
--- a/GeoLib/AnalyticalGeometry.cpp
+++ b/GeoLib/AnalyticalGeometry.cpp
@@ -20,12 +20,13 @@
 
 #include <logog/include/logog.hpp>
 
+#include <Eigen/Dense>
+
 #include "BaseLib/StringTools.h"
 
 #include "Polyline.h"
 #include "PointVec.h"
 
-#include "MathLib/LinAlg/Solvers/GaussAlgorithm.h"
 #include "MathLib/GeometricBasics.h"
 
 extern double orient2d(double *, double *, double *);
@@ -172,16 +173,16 @@ bool lineSegmentIntersect(GeoLib::LineSegment const& s0,
     const double sqr_len_v(v.getSqrLength());
     const double sqr_len_w(w.getSqrLength());
 
-    MathLib::DenseMatrix<double> mat(2,2);
+    Eigen::Matrix2d mat;
     mat(0,0) = sqr_len_v;
     mat(0,1) = -1.0 * MathLib::scalarProduct(v,w);
     mat(1,1) = sqr_len_w;
     mat(1,0) = mat(0,1);
 
-    double rhs[2] = {MathLib::scalarProduct(v,qp), MathLib::scalarProduct(w,pq)};
+    Eigen::Vector2d rhs;
+    rhs << MathLib::scalarProduct(v, qp), MathLib::scalarProduct(w, pq);
 
-    MathLib::GaussAlgorithm<MathLib::DenseMatrix<double>, double*> lu;
-    lu.solve(mat, rhs, true);
+    rhs = mat.partialPivLu().solve(rhs);
 
     // no theory for the following tolerances, determined by testing
     // lower tolerance: little bit smaller than zero
@@ -505,16 +506,15 @@ std::vector<MathLib::Point3d> lineSegmentIntersect2d(
     // at this point it is sure that there is an intersection and the system of
     // linear equations will be invertible
     // solve the two linear equations (b-a, c-d) (t, s)^T = (c-a) simultaneously
-    MathLib::DenseMatrix<double, std::size_t> mat(2,2);
+    Eigen::Matrix2d mat;
     mat(0,0) = b[0]-a[0];
     mat(0,1) = c[0]-d[0];
     mat(1,0) = b[1]-a[1];
     mat(1,1) = c[1]-d[1];
-    std::vector<double> rhs = {{c[0]-a[0], c[1]-a[1]}};
+    Eigen::Vector2d rhs;
+    rhs << c[0] - a[0], c[1] - a[1];
 
-    MathLib::GaussAlgorithm<
-        MathLib::DenseMatrix<double, std::size_t>, std::vector<double>> solver;
-    solver.solve(mat, rhs);
+    rhs = mat.partialPivLu().solve(rhs);
     if (0 <= rhs[1] && rhs[1] <= 1.0) {
         return { MathLib::Point3d{std::array<double,3>{{
                 c[0]+rhs[1]*(d[0]-c[0]), c[1]+rhs[1]*(d[1]-c[1]),
diff --git a/GeoLib/Triangle.cpp b/GeoLib/Triangle.cpp
index 1b6ec36420ee5d865b6c49c49ac89a6939bb6752..87ee3491eefea2bee20628298ddf4ae78baf4616 100644
--- a/GeoLib/Triangle.cpp
+++ b/GeoLib/Triangle.cpp
@@ -10,10 +10,11 @@
 
 #include "Triangle.h"
 
+#include <Eigen/Dense>
+
 #include "Point.h"
 #include "AnalyticalGeometry.h"
 
-#include "MathLib/LinAlg/Solvers/GaussAlgorithm.h"
 #include "MathLib/GeometricBasics.h"
 
 namespace GeoLib {
@@ -49,15 +50,15 @@ bool Triangle::containsPoint2D (Point const& pnt) const
     GeoLib::Point const& c (*(_pnts[_pnt_ids[2]]));
 
     // criterion: p-a = u0 * (b-a) + u1 * (c-a); 0 <= u0, u1 <= 1, u0+u1 <= 1
-    MathLib::DenseMatrix<double> mat (2,2);
+    Eigen::Matrix2d mat;
     mat(0,0) = b[0] - a[0];
     mat(0,1) = c[0] - a[0];
     mat(1,0) = b[1] - a[1];
     mat(1,1) = c[1] - a[1];
-    double y[2] = {pnt[0]-a[0], pnt[1]-a[1]};
+    Eigen::Vector2d y;
+    y << pnt[0]-a[0], pnt[1]-a[1];
 
-    MathLib::GaussAlgorithm<MathLib::DenseMatrix<double>, double*> gauss;
-    gauss.solve(mat, y, true);
+    y = mat.partialPivLu().solve(y);
 
     const double delta (std::numeric_limits<double>::epsilon());
     const double upper (1+delta);
@@ -72,7 +73,7 @@ void getPlaneCoefficients(Triangle const& tri, double c[3])
     GeoLib::Point const& p0 (*(tri.getPoint(0)));
     GeoLib::Point const& p1 (*(tri.getPoint(1)));
     GeoLib::Point const& p2 (*(tri.getPoint(2)));
-    MathLib::DenseMatrix<double> mat (3,3);
+    Eigen::Matrix3d mat;
     mat(0,0) = p0[0];
     mat(0,1) = p0[1];
     mat(0,2) = 1.0;
@@ -82,12 +83,10 @@ void getPlaneCoefficients(Triangle const& tri, double c[3])
     mat(2,0) = p2[0];
     mat(2,1) = p2[1];
     mat(2,2) = 1.0;
-    c[0] = p0[2];
-    c[1] = p1[2];
-    c[2] = p2[2];
+    Eigen::Vector3d y;
+    y << p0[2], p1[2], p2[2];
 
-    MathLib::GaussAlgorithm<MathLib::DenseMatrix<double>, double*> gauss;
-    gauss.solve (mat, c);
+    Eigen::Map<Eigen::Vector3d>(c,3) = mat.partialPivLu().solve(y);
 }
 
 } // end namespace GeoLib
diff --git a/MathLib/CMakeLists.txt b/MathLib/CMakeLists.txt
index 1695d4b435d26bfb7366d3709ee454a8e3223535..b26bbc4b792736749206d3f176c31ca3aa9b5a90 100644
--- a/MathLib/CMakeLists.txt
+++ b/MathLib/CMakeLists.txt
@@ -6,7 +6,6 @@ APPEND_SOURCE_FILES(SOURCES InterpolationAlgorithms)
 APPEND_SOURCE_FILES(SOURCES Integration)
 APPEND_SOURCE_FILES(SOURCES LinAlg)
 APPEND_SOURCE_FILES(SOURCES LinAlg/Dense)
-APPEND_SOURCE_FILES(SOURCES LinAlg/Solvers)
 APPEND_SOURCE_FILES(SOURCES ODE)
 APPEND_SOURCE_FILES(SOURCES Nonlinear)
 
diff --git a/MathLib/GeometricBasics.cpp b/MathLib/GeometricBasics.cpp
index a80e22aaccc87004563c606fae227ed1a81ca5f1..44ce5bd5941e07de9e15d5d85f5a820ba30d17c9 100644
--- a/MathLib/GeometricBasics.cpp
+++ b/MathLib/GeometricBasics.cpp
@@ -10,8 +10,8 @@
 
 #include <logog/include/logog.hpp>
 
-#include "MathLib/LinAlg/Dense/DenseMatrix.h"
-#include "MathLib/LinAlg/Solvers/GaussAlgorithm.h"
+#include <Eigen/Dense>
+
 #include "Point3d.h"
 #include "Vector3.h"
 
@@ -112,17 +112,16 @@ bool gaussPointInTriangle(MathLib::Point3d const& q,
     MathLib::Vector3 const v(a, b);
     MathLib::Vector3 const w(a, c);
 
-    MathLib::DenseMatrix<double> mat(2, 2);
+    Eigen::Matrix2d mat;
     mat(0, 0) = v.getSqrLength();
     mat(0, 1) = v[0] * w[0] + v[1] * w[1] + v[2] * w[2];
     mat(1, 0) = mat(0, 1);
     mat(1, 1) = w.getSqrLength();
-    double y[2] = {
-        v[0] * (q[0] - a[0]) + v[1] * (q[1] - a[1]) + v[2] * (q[2] - a[2]),
-        w[0] * (q[0] - a[0]) + w[1] * (q[1] - a[1]) + w[2] * (q[2] - a[2])};
+    Eigen::Vector2d y;
+    y << v[0] * (q[0] - a[0]) + v[1] * (q[1] - a[1]) + v[2] * (q[2] - a[2]),
+        w[0] * (q[0] - a[0]) + w[1] * (q[1] - a[1]) + w[2] * (q[2] - a[2]);
 
-    MathLib::GaussAlgorithm<MathLib::DenseMatrix<double>, double*> gauss;
-    gauss.solve(mat, y);
+    y = mat.partialPivLu().solve(y);
 
     const double lower(eps_pnt_out_of_tri);
     const double upper(1 + lower);
@@ -171,15 +170,15 @@ bool isPointInTriangleXY(MathLib::Point3d const& p,
                          MathLib::Point3d const& c)
 {
     // criterion: p-a = u0 * (b-a) + u1 * (c-a); 0 <= u0, u1 <= 1, u0+u1 <= 1
-    MathLib::DenseMatrix<double> mat(2, 2);
+    Eigen::Matrix2d mat;
     mat(0, 0) = b[0] - a[0];
     mat(0, 1) = c[0] - a[0];
     mat(1, 0) = b[1] - a[1];
     mat(1, 1) = c[1] - a[1];
-    double y[2] = {p[0] - a[0], p[1] - a[1]};
+    Eigen::Vector2d y;
+    y << p[0] - a[0], p[1] - a[1];
 
-    MathLib::GaussAlgorithm<MathLib::DenseMatrix<double>, double*> gauss;
-    gauss.solve(mat, y, true);
+    y = mat.partialPivLu().solve(y);
 
     // check if u0 and u1 fulfills the condition
     return 0 <= y[0] && y[0] <= 1 && 0 <= y[1] && y[1] <= 1 && y[0] + y[1] <= 1;
diff --git a/MathLib/LinAlg/Solvers/GaussAlgorithm-impl.h b/MathLib/LinAlg/Solvers/GaussAlgorithm-impl.h
deleted file mode 100644
index b32b81fdaa2ae7a3f2d41df051e241f2d26b9c23..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Solvers/GaussAlgorithm-impl.h
+++ /dev/null
@@ -1,113 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date   2011-05-06
- * \brief  Implementation of the GaussAlgorithm class.
- *
- * \copyright
- * Copyright (c) 2012-2017, 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 <cmath>
-#include <algorithm>
-#include "GaussAlgorithm.h"
-
-namespace MathLib
-{
-
-template <typename MAT_T, typename VEC_T>
-void GaussAlgorithm<MAT_T, VEC_T>::performLU(MAT_T& A)
-{
-    IDX_T const nr(A.getNumberOfRows());
-    IDX_T const nc(A.getNumberOfColumns());
-
-    for (IDX_T k=0; k<nc; k++) {
-        // search pivot
-        FP_T t = std::abs(A(k, k));
-        _perm[k] = k;
-        for (IDX_T i=k+1; i<nr; i++) {
-            FP_T const s = std::abs(A(i,k));
-            if (s > t) {
-                t = s;
-                _perm[k] = i;
-            }
-        }
-
-        // exchange rows
-        if (_perm[k] != k) {
-            for (IDX_T j=0; j<nc; j++)
-                std::swap (A(_perm[k],j), A(k,j));
-        }
-
-        // eliminate
-        for (IDX_T i=k+1; i<nr; i++) {
-            FP_T const l = A(i,k)/A(k,k);
-            for (IDX_T j=k; j<nc; j++) {
-                A(i,j) -= A(k,j) * l;
-            }
-            A(i,k) = l;
-        }
-    }
-}
-
-template <typename MAT_T, typename VEC_T>
-template <typename V>
-void GaussAlgorithm<MAT_T, VEC_T>::
-solve (MAT_T& A, V& b, bool decompose)
-{
-    _perm.resize(A.getNumberOfRows());
-
-    if (decompose)
-        performLU(A);
-    permuteRHS (b);
-    forwardSolve (A, b); // L z = b, b will be overwritten by z
-    backwardSolve (A, b); // U x = z, b (z) will be overwritten by x
-}
-
-template <typename MAT_T, typename VEC_T>
-void GaussAlgorithm<MAT_T, VEC_T>::
-solve (MAT_T& A, FP_T* & b, bool decompose)
-{
-    _perm.resize(A.getNumberOfRows());
-
-    if (decompose)
-        performLU(A);
-    permuteRHS (b);
-    forwardSolve (A, b); // L z = b, b will be overwritten by z
-    backwardSolve (A, b); // U x = z, b (z) will be overwritten by x
-}
-
-template <typename MAT_T, typename VEC_T>
-void GaussAlgorithm<MAT_T, VEC_T>::solve (
-        MAT_T& A, VEC_T const& b, VEC_T & x,
-        bool decompose)
-{
-    for (std::size_t k(0); k<A.getNumberOfRows(); k++)
-        x[k] = b[k];
-    solve(A, x, decompose);
-}
-
-template <typename MAT_T, typename VEC_T>
-template <typename V>
-void GaussAlgorithm<MAT_T, VEC_T>::permuteRHS (V & b) const
-{
-    for (IDX_T i=0; i<_perm.size(); i++) {
-        if (_perm[i] != i)
-            std::swap(b[i], b[_perm[i]]);
-    }
-}
-
-template <typename MAT_T, typename VEC_T>
-void GaussAlgorithm<MAT_T, VEC_T>::permuteRHS (VEC_T& b) const
-{
-    for (IDX_T i=0; i<_perm.size(); i++) {
-        if (_perm[i] != i)
-            std::swap(b[i], b[_perm[i]]);
-    }
-}
-
-} // end namespace MathLib
diff --git a/MathLib/LinAlg/Solvers/GaussAlgorithm.h b/MathLib/LinAlg/Solvers/GaussAlgorithm.h
deleted file mode 100644
index 068b7494d7d786adcf8ca8c849ef425d26b0e3c1..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Solvers/GaussAlgorithm.h
+++ /dev/null
@@ -1,101 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date   2011-05-06
- * \brief  Definition of the GaussAlgorithm class.
- *
- * \copyright
- * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- */
-
-#pragma once
-
-#include <vector>
-
-
-#include "BaseLib/ConfigTree.h"
-#include "../Dense/DenseMatrix.h"
-#include "TriangularSolve.h"
-
-namespace MathLib {
-
-/**
- * This is a class for the direct solution of (dense) systems of
- * linear equations, \f$A x = b\f$. During the construction of
- * the object the matrix A is factorized in matrices L and U using
- * Gauss-Elimination with partial pivoting (rows are exchanged). In doing so
- * the entries of A change! The solution for a specific
- * right hand side is computed by the method execute().
- */
-template <typename MAT_T, typename VEC_T = typename MAT_T::FP_T*>
-class GaussAlgorithm
-{
-public:
-    using FP_T = typename MAT_T::FP_T;
-    using IDX_T = typename MAT_T::IDX_T;
-
-public:
-    /**
-     * A direct solver for the (dense) linear system \f$A x = b\f$.
-     * @param solver_name A name used as a prefix for command line options
-     *                    if there are such options available.
-     * @param option For some solvers the user can give parameters to the
-     * algorithm. GaussAlgorithm has to fulfill the common interface
-     * of all solvers of systems of linear equations. For this reason the
-     * second argument was introduced.
-     */
-    GaussAlgorithm(const std::string solver_name = "",
-                   BaseLib::ConfigTree const*const option = nullptr)
-    {
-        (void) solver_name; (void) option; // silence both compiler and doxygen warnings.
-    }
-
-    /**
-     * Method solves the linear system \f$A x = b\f$ (based on the LU factorization)
-     * using forward solve and backward solve.
-     * @param A the coefficient matrix
-     * @param b at the beginning the right hand side, at the end the solution
-     * @param decompose Flag that signals if the LU decomposition should be
-     *        performed or not. If the matrix \f$A\f$ does not change, the LU
-     *        decomposition needs to be performed once only!
-     * @attention The entries of the given matrix will be changed!
-     */
-    template <typename V>
-    void solve (MAT_T& A, V & b, bool decompose = true);
-
-    void solve(MAT_T& A, FP_T* & b, bool decompose = true);
-
-    /**
-     * Method solves the linear system \f$A x = b\f$ (based on the LU factorization)
-     * using forward solve and backward solve.
-     * @param A (input) the coefficient matrix
-     * @param b (input) the right hand side
-     * @param x (output) the solution
-     * @param decompose see documentation of the other solve methods.
-     * @attention The entries of the given matrix will be changed!
-     */
-    void solve(MAT_T& A, VEC_T const& b, VEC_T & x, bool decompose = true);
-
-private:
-    // void solve (MAT_T& A, VEC_T const& b, bool decompose);
-
-    void performLU(MAT_T& A);
-    /**
-     * permute the right hand side vector according to the
-     * row permutations of the LU factorization
-     * @param b the entries of the vector b are permuted
-     */
-    template <typename V> void permuteRHS(V & b) const;
-    void permuteRHS (VEC_T& b) const;
-
-    //! the permutation of the rows
-    std::vector<IDX_T> _perm;
-};
-
-} // end namespace MathLib
-
-#include "GaussAlgorithm-impl.h"
diff --git a/MathLib/LinAlg/Solvers/TriangularSolve-impl.h b/MathLib/LinAlg/Solvers/TriangularSolve-impl.h
deleted file mode 100644
index 0836ce16e65e214ce9dbb86f566359f230795d5d..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Solvers/TriangularSolve-impl.h
+++ /dev/null
@@ -1,64 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date   2011-05-05
- * \brief  Implementation of triangular solver functions.
- *
- * \copyright
- * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- */
-
-namespace MathLib {
-
-template <typename FP_T, typename VEC_T>
-void forwardSolve (const DenseMatrix <FP_T> &L, VEC_T& b)
-{
-    using IDX_T = typename DenseMatrix<double>::IDX_T;
-    IDX_T m (L.getNumberOfRows());
-    FP_T t;
-
-    for (IDX_T r=0; r<m; r++) {
-        t = 0.0;
-        for (IDX_T c=0; c<r; c++) {
-            t += L(r,c)*b[c];
-        }
-        b[r] = b[r]-t;
-    }
-}
-
-template <typename FP_T, typename VEC_T>
-void backwardSolve (const DenseMatrix <FP_T> &mat, VEC_T& b)
-{
-    FP_T t;
-    using IDX_T = typename DenseMatrix<double>::IDX_T;
-    IDX_T m (mat.getNumberOfRows()), n(mat.getNumberOfColumns());
-    for (int r=m-1; r>=0; r--) {
-        t = 0.0;
-        for (IDX_T c=r+1; c<n; c++) {
-            t += mat(r,c)*b[c];
-        }
-        b[r] = (b[r]-t) / mat(r,r);
-    }
-}
-
-template <typename FP_T, typename VEC_T>
-void backwardSolve ( DenseMatrix<FP_T> const& mat, VEC_T& x, VEC_T const& b)
-{
-    using IDX_T = typename DenseMatrix<FP_T>::IDX_T;
-    IDX_T n_cols (mat.getNumberOfColumns());
-    for (int r = (n_cols - 1); r >= 0; r--) {
-        FP_T t = 0.0;
-
-        for (IDX_T c = r+1; c < n_cols; c++) {
-            t += mat(r,c) * b[c];
-        }
-        x[r] = (b[r] - t) / mat(r, r);
-    }
-}
-
-
-} // end namespace MathLib
diff --git a/MathLib/LinAlg/Solvers/TriangularSolve.h b/MathLib/LinAlg/Solvers/TriangularSolve.h
deleted file mode 100644
index 51543ce397fc84a67ebda53f28b5e3ffd944c20b..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Solvers/TriangularSolve.h
+++ /dev/null
@@ -1,52 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date   2011-05-06
- * \brief  Definition of triangular solver functions.
- *
- * \copyright
- * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- */
-
-#pragma once
-
-#include "../Dense/DenseMatrix.h"
-
-namespace MathLib {
-
-/**
- * solves the \f$n \times n\f$ triangular linear system \f$L \cdot y = b\f$,
- * assumes \f$L_{ii} = 1.0\f$, \f$i=1,...,n\f$, \f$b\f$ is destroyed
- * @param L the lower triangular matrix
- * @param b at beginning the right hand side vector, at the end the solution vector
- */
-template <typename FP_T, typename VEC_T>
-void forwardSolve (const DenseMatrix <FP_T> &L, VEC_T& b);
-
-/**
- * solves the \f$n \times n\f$ triangular linear system \f$U \cdot x=y\f$,
- * \f$U\f$, where \f$U\f$ is a upper triangular matrix.
- * @param U upper triangular matrix
- * @param y at beginning the right hand side, at the end the solution
- */
-template <typename FP_T, typename VEC_T>
-void backwardSolve (const DenseMatrix <FP_T> &U, VEC_T& y);
-
-// backwardSolve mat * x = y, mat ... upper triangular matrix
-/**
- * backward solve the system of linear equations \f$ U \cdot x = y\f$,
- * where \f$U\f$ is a upper triangular matrix
- * @param mat the upper triangular matrix
- * @param x the solution of the system of linear equations
- * @param b the right hand side
- */
-template <typename FP_T, typename VEC_T>
-void backwardSolve ( DenseMatrix<FP_T> const& mat, VEC_T& x, VEC_T const& b);
-
-} // end namespace MathLib
-
-#include "TriangularSolve-impl.h"
diff --git a/SimpleTests/MatrixTests/CMakeLists.txt b/SimpleTests/MatrixTests/CMakeLists.txt
deleted file mode 100644
index f5db03f53f0d95ecfcd3fd16312c192f92720a89..0000000000000000000000000000000000000000
--- a/SimpleTests/MatrixTests/CMakeLists.txt
+++ /dev/null
@@ -1,12 +0,0 @@
-# Create the executable
-add_executable(DenseGaussEliminationChecker
-    DenseGaussEliminationChecker.cpp
-    ${SOURCES}
-    ${HEADERS}
-)
-set_target_properties(DenseGaussEliminationChecker PROPERTIES FOLDER SimpleTests)
-target_link_libraries(DenseGaussEliminationChecker
-    logog
-    BaseLib
-    MathLib
-)
diff --git a/SimpleTests/MatrixTests/DenseGaussEliminationChecker.cpp b/SimpleTests/MatrixTests/DenseGaussEliminationChecker.cpp
deleted file mode 100644
index 6586c0affc0d4cfcee90395f43eb5879a78e39ae..0000000000000000000000000000000000000000
--- a/SimpleTests/MatrixTests/DenseGaussEliminationChecker.cpp
+++ /dev/null
@@ -1,83 +0,0 @@
-/**
- * \date   2014-06-11
- * \brief  Implementation of tests.
- *
- * \copyright
- * Copyright (c) 2012-2017, 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 <fstream>
-#include <sstream>
-
-#include <tclap/CmdLine.h>
-#include <logog/include/logog.hpp>
-#include <logog/include/formatter.hpp>
-
-#include "BaseLib/LogogSimpleFormatter.h"
-#include "MathLib/LinAlg/Dense/DenseMatrix.h"
-#include "MathLib/LinAlg/Solvers/GaussAlgorithm.h"
-
-int main(int argc, char *argv[])
-{
-    LOGOG_INITIALIZE();
-    auto* custom_format(new BaseLib::LogogSimpleFormatter);
-    auto* logogCout(new logog::Cout);
-    logogCout->SetFormatter(*custom_format);
-
-    TCLAP::CmdLine cmd("Simple direct matrix solver test.\n\
-        It consists of the following steps:\n\
-        (1) Read a matrix A from ascii format\n\
-        (2) Set all entries of a vector x to one and compute b = A * x\n\
-        (3) Solve the system of linear equations -> result have to be (1,...,1)", ' ', "0.1");
-    TCLAP::ValueArg<std::string> matrix_arg("m", "matrix", "input matrix file (ascii format)", true, "", "string");
-    cmd.add( matrix_arg );
-    cmd.parse( argc, argv );
-
-    // *** reading dense matrix in ascii format from file
-    std::string const fname_mat(matrix_arg.getValue());
-    std::ifstream in(fname_mat.c_str());
-    if (!in) {
-        INFO("error reading matrix from %s", fname_mat.c_str());
-        return -1;
-    }
-    INFO("reading matrix from %s ...", fname_mat.c_str());
-    std::size_t n_rows(0), n_cols(0);
-    in >> n_rows;
-    in >> n_cols;
-    MathLib::DenseMatrix<double, std::size_t> mat(n_rows, n_cols);
-    for (std::size_t i(0); i<mat.getNumberOfRows(); ++i) {
-        for (std::size_t j(0); j<mat.getNumberOfColumns(); ++j) {
-            in >> mat(i,j);
-        }
-    }
-    {
-        std::stringstream stream;
-        stream << mat;
-        INFO("read matrix:\n%s", stream.str().c_str());
-    }
-
-    std::vector<double> x(n_cols,1.0), b;
-    b.resize(n_rows);
-    b = mat * x;
-
-    MathLib::GaussAlgorithm<MathLib::DenseMatrix<double, std::size_t>> gauss;
-    gauss.solve(mat, b, true);
-
-    {
-        std::stringstream stream;
-        std::copy(b.begin(), b.end(), std::ostream_iterator<double>(stream, " "));
-        stream << std::endl;
-        INFO("solution vector:\n%s", stream.str().c_str());
-    }
-
-    delete custom_format;
-    delete logogCout;
-    LOGOG_SHUTDOWN();
-
-    return 0;
-}
-
diff --git a/Tests/MathLib/TestDenseGaussAlgorithm.cpp b/Tests/MathLib/TestDenseGaussAlgorithm.cpp
deleted file mode 100644
index f7a9703501ae4b6526639d66f3218124fe0d32f8..0000000000000000000000000000000000000000
--- a/Tests/MathLib/TestDenseGaussAlgorithm.cpp
+++ /dev/null
@@ -1,97 +0,0 @@
-/**
- * @file TestDenseGaussAlgorithm.cpp
- * @author Thomas Fischer
- * @date Jun 17, 2013
- * @brief
- *
- * @copyright
- * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/LICENSE.txt
- */
-
-#include <cstdlib>
-#include <ctime>
-#include <limits>
-#include <algorithm>
-
-#include <gtest/gtest.h>
-
-#include "MathLib/LinAlg/Solvers/GaussAlgorithm.h"
-
-TEST(MathLib, DenseGaussAlgorithm)
-{
-    std::size_t n_rows(100);
-    std::size_t n_cols(n_rows);
-
-    MathLib::DenseMatrix<double,std::size_t> mat(n_rows, n_cols);
-
-    // *** fill matrix with arbitrary values
-    // ** initialize random seed
-    srand ( static_cast<unsigned>(time(nullptr)) );
-    // ** loop over rows and columns
-    for (std::size_t i(0); i<n_rows; i++) {
-        for (std::size_t j(0); j<n_cols; j++) {
-            mat(i,j) = rand()/static_cast<double>(RAND_MAX);
-        }
-    }
-
-    // *** create solution vector, set all entries to 0.0
-    auto* x(new double[n_cols]);
-    std::fill(x,x+n_cols, 0.0);
-    double *b0(mat * x);
-
-    // *** create other right hand sides,
-    // set all entries of the solution vector to 1.0
-    std::fill(x,x+n_cols, 1.0);
-    double *b1(mat * x);
-
-    std::generate(x,x+n_cols, std::rand);
-    double *b2(mat * x);
-
-    // right hand side and solution vector with random entries
-    double *b3(mat * x);
-    double *b3_copy(mat * x);
-    auto* x3(new double[n_cols]);
-    std::generate(x3,x3+n_cols, std::rand);
-
-    MathLib::GaussAlgorithm<MathLib::DenseMatrix<double, std::size_t>, double*> gauss;
-
-    // solve with b0 as right hand side
-    gauss.solve(mat, b0, true);
-    for (std::size_t i(0); i<n_rows; i++) {
-        ASSERT_NEAR(b0[i], 0.0, std::numeric_limits<float>::epsilon());
-    }
-
-    // solve with b1 as right hand side
-    gauss.solve(mat, b1, false);
-    for (std::size_t i(0); i<n_rows; i++) {
-        ASSERT_NEAR(b1[i], 1.0, std::numeric_limits<float>::epsilon());
-    }
-
-    // solve with b2 as right hand side
-    gauss.solve(mat, b2, false);
-    for (std::size_t i(0); i<n_rows; i++) {
-        ASSERT_NEAR(fabs(b2[i]-x[i])/fabs(x[i]), 0.0, std::numeric_limits<float>::epsilon());
-    }
-
-    // solve with b3 as right hand side and x3 as solution vector
-    gauss.solve(mat, b3, x3, false);
-    for (std::size_t i(0); i<n_rows; i++) {
-        ASSERT_NEAR(fabs(x3[i]-x[i])/fabs(x[i]), 0.0, std::numeric_limits<float>::epsilon());
-    }
-    // assure entries of vector b3 are not changed
-    for (std::size_t i(0); i<n_rows; i++) {
-        ASSERT_NEAR(fabs(b3[i]-b3_copy[i])/fabs(b3[i]), 0.0, std::numeric_limits<float>::epsilon());
-    }
-
-    delete [] x;
-    delete [] b0;
-    delete [] b1;
-    delete [] b2;
-    delete [] b3;
-    delete [] x3;
-    delete [] b3_copy;
-}
-
diff --git a/Tests/MathLib/TestLinearSolver.cpp b/Tests/MathLib/TestLinearSolver.cpp
index 23f69595d56c0da5f0440b41c43aec9943e687a0..85940cbc316b3d3c8b0efab21ca1eb6e01cb1df9 100644
--- a/Tests/MathLib/TestLinearSolver.cpp
+++ b/Tests/MathLib/TestLinearSolver.cpp
@@ -19,7 +19,6 @@
 #include "MathLib/LinAlg/Dense/DenseMatrix.h"
 #include "MathLib/LinAlg/FinalizeMatrixAssembly.h"
 #include "MathLib/LinAlg/ApplyKnownSolution.h"
-#include "MathLib/LinAlg/Solvers/GaussAlgorithm.h"
 
 #ifdef OGS_USE_EIGEN
 #include "MathLib/LinAlg/Eigen/EigenMatrix.h"
diff --git a/Tests/NumLib/TestSerialLinearSolver.cpp b/Tests/NumLib/TestSerialLinearSolver.cpp
index ab5ab0758f4f332389a3b5474d6410d4054f6ef3..1f3da8cef622f46d428e9ead15e7b9b904d695f9 100644
--- a/Tests/NumLib/TestSerialLinearSolver.cpp
+++ b/Tests/NumLib/TestSerialLinearSolver.cpp
@@ -20,7 +20,6 @@
 #include "MathLib/LinAlg/GlobalMatrixVectorTypes.h"
 #include "MathLib/LinAlg/MatrixSpecifications.h"
 #include "MathLib/LinAlg/MatrixVectorTraits.h"
-#include "MathLib/LinAlg/Solvers/GaussAlgorithm.h"
 #include "MathLib/LinAlg/FinalizeMatrixAssembly.h"
 #include "MathLib/MathTools.h"