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/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"