diff --git a/MathLib/CMakeLists.txt b/MathLib/CMakeLists.txt
index 763b6aabc5bc1b5f4326651484007fc579e3122f..05603a7f63c0fe792d896c0610357efb1019733b 100644
--- a/MathLib/CMakeLists.txt
+++ b/MathLib/CMakeLists.txt
@@ -5,7 +5,6 @@ append_source_files(SOURCES Curve)
 append_source_files(SOURCES InterpolationAlgorithms)
 append_source_files(SOURCES Integration)
 append_source_files(SOURCES LinAlg)
-append_source_files(SOURCES LinAlg/Dense)
 if(CVODE_FOUND)
     append_source_files(SOURCES ODE)
 endif()
diff --git a/MathLib/LinAlg/Dense/DenseMatrix-impl.h b/MathLib/LinAlg/Dense/DenseMatrix-impl.h
deleted file mode 100644
index b83c89b0b6c5055fea960fed6b749d185d5fae58..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Dense/DenseMatrix-impl.h
+++ /dev/null
@@ -1,380 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer and Haibing Shao
- * \date Jun 10, 2013
- *
- * \copyright
- * Copyright (c) 2012-2020, 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 <cassert>
-
-namespace MathLib
-{
-
-template<typename FP_TYPE, typename IDX_TYPE>
-DenseMatrix<FP_TYPE, IDX_TYPE>::DenseMatrix(IDX_TYPE rows, IDX_TYPE cols) :
-        _n_rows(rows), _n_cols(cols), _data(new FP_TYPE[_n_rows * _n_cols])
-{}
-
-template<typename FP_TYPE, typename IDX_TYPE>
-
-DenseMatrix<FP_TYPE, IDX_TYPE>::DenseMatrix(IDX_TYPE rows, IDX_TYPE cols,
-        FP_TYPE const& initial_value) :
-        _n_rows(rows), _n_cols(cols), _data(new FP_TYPE[_n_rows * _n_cols])
-{
-    const IDX_TYPE n(_n_rows * _n_cols);
-    for (IDX_TYPE k(0); k < n; k++)
-    {
-        _data[k] = initial_value;
-    }
-}
-
-template<typename FP_TYPE, typename IDX_TYPE>
-DenseMatrix<FP_TYPE, IDX_TYPE>::DenseMatrix (const DenseMatrix<FP_TYPE, IDX_TYPE>& src) :
-        _n_rows(src.getNumberOfRows ()), _n_cols(src.getNumberOfColumns ()), _data (new FP_TYPE[_n_rows * _n_cols])
-{
-    std::copy(src._data, src._data+_n_rows*_n_cols, _data);
-}
-
-template<typename FP_TYPE, typename IDX_TYPE>
-DenseMatrix<FP_TYPE, IDX_TYPE>::DenseMatrix (DenseMatrix<FP_TYPE, IDX_TYPE> &&src) :
-       _n_rows(src.getNumberOfRows()), _n_cols(src.getNumberOfColumns())
-{
-    src._n_rows = 0;
-    src._n_cols = 0;
-    _data = src._data;
-    src._data = nullptr;
-}
-
-
-template <typename FP_TYPE, typename IDX_TYPE>
-DenseMatrix<FP_TYPE, IDX_TYPE>::~DenseMatrix ()
-{
-   delete [] _data;
-}
-
-template <typename FP_TYPE, typename IDX_TYPE>
-DenseMatrix<FP_TYPE, IDX_TYPE>&
-DenseMatrix<FP_TYPE, IDX_TYPE>::operator=(DenseMatrix<FP_TYPE, IDX_TYPE> const& rhs)
-{
-    if (this == &rhs)
-    {
-        return *this;
-    }
-
-    if (_n_rows != rhs.getNumberOfRows() || _n_cols != rhs.getNumberOfColumns()) {
-        std::string msg("DenseMatrix::operator=(DenseMatrix const& rhs), Dimension mismatch, ");
-        msg += " left hand side: " + std::to_string(_n_rows) + " x "
-                + std::to_string(_n_cols);
-        msg += " right hand side: " + std::to_string(rhs.getNumberOfRows()) + " x "
-                + std::to_string(rhs.getNumberOfColumns());
-        throw std::range_error(msg);
-        return *this;
-    }
-
-    std::copy(rhs._data, rhs._data + _n_rows * _n_cols, _data);
-
-    return *this;
-}
-
-template<typename FP_TYPE, typename IDX_TYPE>
-DenseMatrix<FP_TYPE, IDX_TYPE>&
-DenseMatrix<FP_TYPE, IDX_TYPE>::operator=(DenseMatrix && rhs)
-{
-    _n_rows = rhs._n_rows;
-    _n_cols = rhs._n_cols;
-    _data = rhs._data;
-
-    rhs._n_rows = 0;
-    rhs._n_cols = 0;
-    rhs._data = nullptr;
-    return *this;
-}
-
-template<typename FP_TYPE, typename IDX_TYPE>
-DenseMatrix<FP_TYPE, IDX_TYPE>&
-DenseMatrix<FP_TYPE, IDX_TYPE>::operator=(FP_TYPE const& v)
-{
-    std::fill(this->_data, this->_data + this->_n_rows * this->_n_cols, v);
-    return *this;
-}
-
-template<typename FP_TYPE, typename IDX_TYPE>
-void DenseMatrix<FP_TYPE, IDX_TYPE>::axpy(FP_TYPE alpha, const FP_TYPE* x, FP_TYPE beta,
-        FP_TYPE* y) const
-{
-    for (IDX_TYPE i(0); i < _n_rows; i++) {
-        y[i] += beta * y[i];
-        for (IDX_TYPE j(0); j < _n_cols; j++) {
-            y[i] += alpha * _data[address(i, j)] * x[j];
-        }
-    }
-}
-
-template<typename FP_TYPE, typename IDX_TYPE>
-FP_TYPE* DenseMatrix<FP_TYPE, IDX_TYPE>::operator* (FP_TYPE* const& x) const
-{
-    return this->operator*(static_cast<FP_TYPE const*>(x));
-}
-
-template<typename FP_TYPE, typename IDX_TYPE>
-FP_TYPE* DenseMatrix<FP_TYPE, IDX_TYPE>::operator* (FP_TYPE const* const& x) const
-{
-    auto* y(new FP_TYPE[_n_rows]);
-    for (IDX_TYPE i(0); i < _n_rows; i++) {
-        y[i] = 0.0;
-        for (IDX_TYPE j(0); j < _n_cols; j++) {
-            y[i] += _data[address(i, j)] * x[j];
-        }
-    }
-
-    return y;
-}
-
-template<typename FP_TYPE, typename IDX_TYPE>
-template <typename V>
-V DenseMatrix<FP_TYPE, IDX_TYPE>::operator* (V const& x) const
-{
-    V y(_n_rows);
-    for (IDX_TYPE i(0); i < _n_rows; i++) {
-        y[i] = 0.0;
-        for (IDX_TYPE j(0); j < _n_cols; j++) {
-            y[i] += _data[address(i, j)] * x[j];
-        }
-    }
-
-    return y;
-}
-
-template<typename FP_TYPE, typename IDX_TYPE>
-MathLib::Vector3
-DenseMatrix<FP_TYPE, IDX_TYPE>::operator*(MathLib::Vector3 const& x) const
-{
-    assert(_n_rows>2);
-
-    MathLib::Vector3 y;
-    for (IDX_TYPE i(0); i < _n_rows; i++) {
-        y[i] = 0.0;
-        for (IDX_TYPE j(0); j < _n_cols; j++) {
-            y[i] += _data[address(i, j)] * x[j];
-        }
-    }
-
-    return y;
-}
-
-template<typename FP_TYPE, typename IDX_TYPE>
-DenseMatrix<FP_TYPE, IDX_TYPE>*
-DenseMatrix<FP_TYPE, IDX_TYPE>::operator+(const DenseMatrix<FP_TYPE, IDX_TYPE>& mat) const
-{
-    // make sure the two matrices have the same dimension.
-    if (_n_rows != mat.getNumberOfRows() || _n_cols != mat.getNumberOfColumns())
-    {
-        throw std::range_error("DenseMatrix::operator+, illegal matrix size!");
-    }
-
-    DenseMatrix<FP_TYPE, IDX_TYPE>* y(new DenseMatrix<FP_TYPE, IDX_TYPE>(_n_rows, _n_cols));
-    for (IDX_TYPE i = 0; i < _n_rows; i++) {
-        for (IDX_TYPE j = 0; j < _n_cols; j++) {
-            (*y)(i, j) = _data[address(i, j)] + mat(i, j);
-        }
-    }
-
-    return y;
-}
-
-template<typename FP_TYPE, typename IDX_TYPE>
-DenseMatrix<FP_TYPE, IDX_TYPE>*
-DenseMatrix<FP_TYPE, IDX_TYPE>::operator-(const DenseMatrix<FP_TYPE, IDX_TYPE>& mat) const
-{
-    // make sure the two matrices have the same dimension.
-    if (_n_rows != mat.getNumberOfRows() || _n_cols != mat.getNumberOfColumns())
-    {
-        throw std::range_error("DenseMatrix::operator-, illegal matrix size!");
-    }
-
-    DenseMatrix<FP_TYPE, IDX_TYPE>* y(new DenseMatrix<FP_TYPE, IDX_TYPE>(_n_rows, _n_cols));
-    for (IDX_TYPE i = 0; i < _n_rows; i++) {
-        for (IDX_TYPE j = 0; j < _n_cols; j++) {
-            (*y)(i, j) = _data[address(i, j)] - mat(i, j);
-        }
-    }
-
-    return y;
-}
-
-template<typename FP_TYPE, typename IDX_TYPE>
-DenseMatrix<FP_TYPE, IDX_TYPE>*
-DenseMatrix<FP_TYPE, IDX_TYPE>::operator*(const DenseMatrix<FP_TYPE, IDX_TYPE>& mat) const
-{
-    // make sure the two matrices have the same dimension.
-    if (_n_cols != mat.getNumberOfRows())
-    {
-        throw std::range_error(
-            "DenseMatrix::operator*, number of rows and cols should be the "
-            "same!");
-    }
-
-    IDX_TYPE y_cols(mat.getNumberOfColumns());
-    DenseMatrix<FP_TYPE, IDX_TYPE>* y(
-            new DenseMatrix<FP_TYPE, IDX_TYPE>(_n_rows, y_cols, FP_TYPE(0)));
-
-    for (IDX_TYPE i = 0; i < _n_rows; i++) {
-        for (IDX_TYPE j = 0; j < y_cols; j++) {
-            for (IDX_TYPE k = 0; k < _n_cols; k++)
-            {
-                (*y)(i, j) += _data[address(i, k)] * mat(k, j);
-            }
-        }
-    }
-
-    return y;
-}
-
-template<typename FP_TYPE, typename IDX_TYPE>
-DenseMatrix<FP_TYPE, IDX_TYPE>*
-DenseMatrix<FP_TYPE, IDX_TYPE>::transpose() const
-{
-    DenseMatrix<FP_TYPE, IDX_TYPE>* y(new DenseMatrix<FP_TYPE, IDX_TYPE>(_n_cols, _n_rows));
-
-    for (IDX_TYPE i = 0; i < _n_rows; i++) {
-        for (IDX_TYPE j = 0; j < _n_cols; j++) {
-            (*y)(j, i) = _data[address(i, j)];
-        }
-    }
-    return y;
-}
-
-template<typename FP_TYPE, typename IDX_TYPE>
-void
-DenseMatrix<FP_TYPE, IDX_TYPE>::transposeInPlace()
-{
-    if (_n_rows==_n_cols) { // square matrix
-        for (IDX_TYPE i = 0; i < _n_rows; i++)
-        {
-            for (IDX_TYPE j = i + 1; j < _n_cols; j++)
-            {
-                std::swap(_data[address(i, j)], _data[address(j, i)]);
-            }
-        }
-    } else { // non-square matrix
-        const DenseMatrix<FP_TYPE, IDX_TYPE> org(*this);
-        std::swap(_n_rows, _n_cols);
-        for (IDX_TYPE i = 0; i < _n_rows; i++) {
-            for (IDX_TYPE j = 0; j < _n_cols; j++) {
-                _data[address(i, j)] = org(j, i);
-            }
-        }
-    }
-
-}
-
-template<typename FP_TYPE, typename IDX_TYPE>
-DenseMatrix<FP_TYPE, IDX_TYPE>*
-DenseMatrix<FP_TYPE, IDX_TYPE>::getSubMatrix(
-        IDX_TYPE b_row, IDX_TYPE b_col,
-        IDX_TYPE e_row, IDX_TYPE e_col) const
-{
-    if (b_row >= e_row | b_col >= e_col)
-    {
-        throw std::range_error(
-            "DenseMatrix::getSubMatrix() illegal sub matrix");
-    }
-    if (e_row > _n_rows | e_col > _n_cols)
-    {
-        throw std::range_error(
-            "DenseMatrix::getSubMatrix() illegal sub matrix");
-    }
-
-    DenseMatrix<FP_TYPE, IDX_TYPE>* y(
-            new DenseMatrix<FP_TYPE, IDX_TYPE>(e_row - b_row, e_col - b_col));
-    for (IDX_TYPE i = b_row; i < e_row; i++) {
-        for (IDX_TYPE j = b_col; j < e_col; j++) {
-            (*y)(i - b_row, j - b_col) = _data[address(i, j)];
-        }
-    }
-    return y;
-}
-
-template<typename FP_TYPE, typename IDX_TYPE>
-void
-DenseMatrix<FP_TYPE, IDX_TYPE>::setSubMatrix(IDX_TYPE b_row, IDX_TYPE b_col,
-        const DenseMatrix<FP_TYPE, IDX_TYPE>& sub_mat)
-{
-    if (b_row + sub_mat.getNumberOfRows() > _n_rows |
-        b_col + sub_mat.getNumberOfColumns() > _n_cols)
-    {
-        throw std::range_error("DenseMatrix::setSubMatrix() sub matrix to big");
-    }
-
-    for (IDX_TYPE i = 0; i < sub_mat.getNumberOfRows(); i++) {
-        for (IDX_TYPE j = 0; j < sub_mat.getNumberOfColumns(); j++) {
-            _data[address(i + b_row, j + b_col)] = sub_mat(i, j);
-        }
-    }
-}
-
-template<typename FP_TYPE, typename IDX_TYPE>
-FP_TYPE&
-DenseMatrix<FP_TYPE, IDX_TYPE>::operator() (IDX_TYPE row, IDX_TYPE col)
-{
-    assert((row < _n_rows) && (col < _n_cols));
-    return _data [address(row,col)];
-}
-
-
-template<typename FP_TYPE, typename IDX_TYPE>
-FP_TYPE const&
-DenseMatrix<FP_TYPE, IDX_TYPE>::operator() (IDX_TYPE row, IDX_TYPE col) const
-{
-    assert((row < _n_rows) && (col < _n_cols));
-    return _data[address(row, col)];
-}
-
-template <typename FP_TYPE, typename IDX_TYPE>
-void
-DenseMatrix<FP_TYPE, IDX_TYPE>::write (std::ostream &out) const
-{
-    out << _n_rows << " " << _n_cols << "\n";
-    for (IDX_TYPE i = 0; i < _n_rows; i++) {
-        for (IDX_TYPE j = 0; j < _n_cols; j++) {
-            out << _data[address(i, j)] << "\t";
-        }
-        out << "\n";
-    }
-}
-
-template <typename FP_TYPE, typename IDX_TYPE>
-void
-DenseMatrix<FP_TYPE, IDX_TYPE>::setIdentity()
-{
-    (*this) = 0.0;
-    const IDX_TYPE n_square_rows = std::min(_n_rows, _n_cols);
-    for (IDX_TYPE i = 0; i < n_square_rows; i++)
-    {
-        _data[address(i, i)] = 1.0;
-    }
-}
-
-template <typename FP_TYPE, typename IDX_TYPE>
-FP_TYPE
-sqrFrobNrm (const DenseMatrix<FP_TYPE, IDX_TYPE> &mat)
-{
-    FP_TYPE nrm(static_cast<FP_TYPE>(0));
-    IDX_TYPE i, j;
-    for (j = 0; j < mat.getNumberOfColumns(); j++)
-    {
-        for (i = 0; i < mat.getNumberOfRows(); i++)
-        {
-            nrm += mat(i, j) * mat(i, j);
-        }
-    }
-
-    return nrm;
-}
-
-} // end namespace MathLib
diff --git a/MathLib/LinAlg/Dense/DenseMatrix.h b/MathLib/LinAlg/Dense/DenseMatrix.h
deleted file mode 100644
index 2d160d01692f9964749ecd58937982a2e3c1d727..0000000000000000000000000000000000000000
--- a/MathLib/LinAlg/Dense/DenseMatrix.h
+++ /dev/null
@@ -1,190 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer and Haibing Shao
- * \date   2011-05-24
- * \brief  Definition of the DenseMatrix class.
- *
- * \copyright
- * Copyright (c) 2012-2020, 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 "MathLib/Vector3.h"
-
-namespace MathLib {
-
-/**
- * Matrix represents a dense matrix for a numeric data type.
- */
-template <typename FP_TYPE, typename IDX_TYPE = std::size_t> class DenseMatrix
-{
-public:
-    using FP_T = FP_TYPE;
-    using IDX_T = IDX_TYPE;
-
-public:
-   /// Dense square matrix constructor.
-   explicit DenseMatrix (IDX_TYPE rows) : DenseMatrix(rows, rows) {}
-
-   /// Dense rectangular matrix constructor.
-   DenseMatrix (IDX_TYPE rows, IDX_TYPE cols);
-   DenseMatrix(IDX_TYPE rows, IDX_TYPE cols, const FP_TYPE& initial_value);
-   DenseMatrix (const DenseMatrix &src);
-   /**
-    * Move constructor.
-    * @param src The original DenseMatrix object. After applying the
-    * move constructor the object src has no rows and columns anymore.
-    *
-    */
-   DenseMatrix (DenseMatrix &&src);
-
-   virtual ~DenseMatrix ();
-
-   /**
-    * Get all entries, which are stored in an array.
-   */
-   const FP_TYPE *getEntries() const
-   {
-       return _data;
-   }
-
-   /**
-    * Get all entries, which are stored in an array.
-   */
-   const FP_TYPE *data() const
-   {
-       return _data;
-   }
-
-   /**
-    * Assignment operator, makes a copy of the internal data of the object.
-    * @param rhs The DenseMatrix object to the right side of the assignment symbol.
-    */
-   DenseMatrix& operator=(DenseMatrix const& rhs);
-
-   /**
-    * This is the move assignment operator.
-    * @param rhs This is the right hand side of a assignment operation.
-    * After applying this operation the object src has no rows and columns anymore.
-    */
-   DenseMatrix& operator=(DenseMatrix && rhs);
-
-   /**
-    * Assignment operator
-    */
-   DenseMatrix& operator=(FP_TYPE const& v);
-
-   /**
-    * \f$ y = \alpha \cdot A x + \beta y\f$
-    */
-   void axpy (FP_TYPE alpha, const FP_TYPE* x, FP_TYPE beta, FP_TYPE* y) const;
-
-   /**
-    * DenseMatrix vector multiplication
-    */
-   FP_TYPE* operator* (FP_TYPE* const& x) const;
-   FP_TYPE* operator* (FP_TYPE const* const& x) const;
-   template <typename V> V operator* (V const& x) const;
-   MathLib::Vector3 operator*(MathLib::Vector3 const& x) const;
-
-   /**
-    * DenseMatrix matrix addition.
-    */
-   DenseMatrix* operator+ (const DenseMatrix& mat) const;
-   /**
-    * DenseMatrix matrix subtraction
-    */
-   DenseMatrix* operator- (const DenseMatrix& mat) const;
-
-   /**
-    * DenseMatrix matrix multiplication \f$ C = A \cdot B\f$
-    * @param mat the matrix \f$ B \f$
-    * @return the matrix \f$ C \f$
-    */
-   DenseMatrix* operator* (const DenseMatrix& mat) const;
-
-   /**
-    * matrix transpose
-    * @return the transpose of the matrix
-    */
-   DenseMatrix* transpose() const; // HB & ZC
-
-   /**
-    * transpose the matrix in place
-    */
-   void transposeInPlace();
-
-   DenseMatrix* getSubMatrix (IDX_TYPE b_row, IDX_TYPE b_col, IDX_TYPE e_row, IDX_TYPE e_col) const;
-
-   /**
-    * overwrites values of the matrix with the given sub matrix
-    * @param b_row the first row
-    * @param b_col the first column
-    * @param sub_mat the sub matrix
-    */
-   void setSubMatrix (IDX_TYPE b_row, IDX_TYPE b_col, const DenseMatrix& sub_mat);
-
-   inline FP_TYPE & operator() (IDX_TYPE row, IDX_TYPE col);
-   inline FP_TYPE const& operator() (IDX_TYPE row, IDX_TYPE col) const;
-
-   /**
-    * writes the matrix entries into the output stream
-    * @param out the output stream
-    */
-   void write (std::ostream& out) const;
-
-    /**
-     * get the number of rows
-     * @return the number of rows
-     */
-    IDX_TYPE getNumberOfRows () const { return _n_rows; }
-    /**
-     * get the number of columns
-     * @return the number of columns
-     */
-    IDX_TYPE getNumberOfColumns () const { return _n_cols; }
-
-    /**
-     * get the number of entries in the matrix
-     */
-    IDX_TYPE size() const { return _n_rows*_n_cols; }
-
-    /**
-     * set the identity matrix
-     */
-    void setIdentity();
-
-protected:
-    /**
-     * the number of rows
-     */
-    IDX_TYPE _n_rows;
-    /**
-     * the number of columns
-     */
-    IDX_TYPE _n_cols;
-
-   // zero based addressing, but Fortran storage layout
-   //inline IDX_TYPE address(IDX_TYPE i, IDX_TYPE j) const { return j*rows+i; }
-   // zero based addressing, C storage layout
-   inline IDX_TYPE address(IDX_TYPE i, IDX_TYPE j) const { return i*_n_cols+j; }
-
-   FP_TYPE *_data;
-};
-
-/** overload the output operator for class DenseMatrix */
-template <typename FP_TYPE, typename IDX_TYPE>
-std::ostream& operator<< (std::ostream &os, const DenseMatrix<FP_TYPE, IDX_TYPE> &mat)
-{
-    mat.write (os);
-    return os;
-}
-
-} // end namespace MathLib
-
-#include "DenseMatrix-impl.h"
diff --git a/Tests/MathLib/TestDenseMatrix.cpp b/Tests/MathLib/TestDenseMatrix.cpp
deleted file mode 100644
index 9d8db7a26fdbb0db89f7bae61ff99019c854b894..0000000000000000000000000000000000000000
--- a/Tests/MathLib/TestDenseMatrix.cpp
+++ /dev/null
@@ -1,79 +0,0 @@
-/**
- * \copyright
- * Copyright (c) 2012-2020, 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 "MathLib/LinAlg/Dense/DenseMatrix.h"
-#include "Tests/TestTools.h"
-
-using namespace MathLib;
-
-TEST(MathLib, DenseMatrixTransposeInPlace)
-{
-    const double eps(std::numeric_limits<double>::epsilon());
-
-    // square matrix
-    DenseMatrix<double> m1(3,3);
-    unsigned cnt = 0;
-    for (unsigned i = 0; i < m1.getNumberOfRows(); i++)
-    {
-        for (unsigned j = 0; j < m1.getNumberOfColumns(); j++)
-        {
-            m1(i, j) = cnt++;
-        }
-    }
-    double expected_m1[] = {0, 1, 2, 3, 4, 5, 6, 7, 8};
-    ASSERT_EQ(3u, m1.getNumberOfRows());
-    ASSERT_EQ(3u, m1.getNumberOfColumns());
-    ASSERT_ARRAY_NEAR(expected_m1, m1.data(), m1.size(), eps);
-    m1.transposeInPlace();
-    ASSERT_EQ(3u, m1.getNumberOfRows());
-    ASSERT_EQ(3u, m1.getNumberOfColumns());
-    double expected_m1t[] = {0, 3, 6, 1, 4, 7, 2, 5, 8};
-    ASSERT_ARRAY_NEAR(expected_m1t, m1.data(), m1.size(), eps);
-
-    // non-square matrix 1
-    DenseMatrix<double> m2(2,3);
-    cnt = 0;
-    for (unsigned i = 0; i < m2.getNumberOfRows(); i++)
-    {
-        for (unsigned j = 0; j < m2.getNumberOfColumns(); j++)
-        {
-            m2(i, j) = cnt++;
-        }
-    }
-    ASSERT_EQ(2u, m2.getNumberOfRows());
-    ASSERT_EQ(3u, m2.getNumberOfColumns());
-    double expected_m2[] = {0, 1, 2, 3, 4, 5};
-    ASSERT_ARRAY_NEAR(expected_m2, m2.data(), m2.size(), eps);
-    m2.transposeInPlace();
-    ASSERT_EQ(3u, m2.getNumberOfRows());
-    ASSERT_EQ(2u, m2.getNumberOfColumns());
-    double expected_m2t[] = {0, 3, 1, 4, 2, 5};
-    ASSERT_ARRAY_NEAR(expected_m2t, m2.data(), m2.size(), eps);
-
-    // non-square matrix 2
-    DenseMatrix<double> m3(3,2);
-    cnt = 0;
-    for (unsigned i = 0; i < m3.getNumberOfRows(); i++)
-    {
-        for (unsigned j = 0; j < m3.getNumberOfColumns(); j++)
-        {
-            m3(i, j) = cnt++;
-        }
-    }
-    ASSERT_EQ(3u, m3.getNumberOfRows());
-    ASSERT_EQ(2u, m3.getNumberOfColumns());
-    double expected_m3[] = {0, 1, 2, 3, 4, 5};
-    ASSERT_ARRAY_NEAR(expected_m3, m3.data(), m3.size(), eps);
-    m3.transposeInPlace();
-    ASSERT_EQ(2u, m3.getNumberOfRows());
-    ASSERT_EQ(3u, m3.getNumberOfColumns());
-    double expected_m3t[] = {0, 2, 4, 1, 3, 5};
-    ASSERT_ARRAY_NEAR(expected_m3t, m3.data(), m3.size(), eps);
-}
diff --git a/Tests/MathLib/TestGlobalMatrixInterface.cpp b/Tests/MathLib/TestGlobalMatrixInterface.cpp
index 161d83e22210041b700f6ef33851c10ed6d740c1..48970c58e24dd1772c2cd2aadc5d21a98eec1bcb 100644
--- a/Tests/MathLib/TestGlobalMatrixInterface.cpp
+++ b/Tests/MathLib/TestGlobalMatrixInterface.cpp
@@ -15,6 +15,8 @@
 
 #include <gtest/gtest.h>
 
+#include <Eigen/Eigen>
+
 #include "MathLib/LinAlg/LinAlg.h"
 
 #if defined(USE_PETSC)
@@ -23,7 +25,6 @@
 #include "MathLib/LinAlg/Eigen/EigenMatrix.h"
 #endif
 
-#include "MathLib/LinAlg/Dense/DenseMatrix.h"
 #include "MathLib/LinAlg/FinalizeMatrixAssembly.h"
 
 #include "NumLib/NumericsConfig.h"
@@ -45,7 +46,8 @@ void checkGlobalMatrixInterface(T_MATRIX &m)
     m.add(0, 0, 1.0);
     m.setZero();
 
-    MathLib::DenseMatrix<double> local_m(2, 2, 1.0);
+    Eigen::Matrix2d local_m;
+    local_m << 1.0, 1.0, 1.0, 1.0;
     std::vector<GlobalIndexType> vec_pos(2);
     vec_pos[0] = 1;
     vec_pos[1] = 3;
@@ -77,7 +79,7 @@ void checkGlobalMatrixInterfaceMPI(T_MATRIX &m, T_VECTOR &v)
     ASSERT_EQ(m.getNumberOfColumns(), gathered_cols);
 
     // Add entries
-    MathLib::DenseMatrix<double> loc_m(2, 2);
+    Eigen::Matrix2d loc_m(2, 2);
     loc_m(0, 0) = 1.;
     loc_m(0, 1) = 2.;
     loc_m(1, 0) = 3.;
@@ -140,7 +142,7 @@ void checkGlobalRectangularMatrixInterfaceMPI(T_MATRIX &m, T_VECTOR &v)
     ASSERT_EQ(m.getNumberOfColumns(), gathered_cols);
 
     // Add entries
-    MathLib::DenseMatrix<double> loc_m(2, 3);
+    Eigen::Matrix<double, 2, 3> loc_m;
     loc_m(0, 0) = 1.;
     loc_m(0, 1) = 2.;
     loc_m(0, 2) = 3.;