Merge pull request #148 from TomFischer/CleanUpSolverStructure

GeoLib/AnalyticalGeometry.cpp

@@ -87,7 +87,7 @@ bool lineSegmentIntersect(const GeoLib::Point& a, const GeoLib::Point& b, const
 	rhs[0] = c[0] - a[0];
 	rhs[1] = c[1] - a[1];
 
-	MathLib::GaussAlgorithm lu_solver (mat);
+	MathLib::GaussAlgorithm<MathLib::DenseMatrix<double>, double*> lu_solver (mat);
 	lu_solver.execute (rhs);
 	if (0 <= rhs[0] && rhs[0] <= 1.0 && 0 <= rhs[1] && rhs[1] <= 1.0) {
 		s[0] = a[0] + rhs[0] * (b[0] - a[0]);
@@ -145,7 +145,7 @@ bool isPointInTriangle(const double p[3], const double a[3], const double b[3],
 	mat(1, 1) = c[1] - b[1];
 	double rhs[2] = { p[0] - b[0], p[1] - b[1] };
 
-	MathLib::GaussAlgorithm gauss(mat);
+	MathLib::GaussAlgorithm<MathLib::DenseMatrix<double>, double*> gauss(mat);
 	gauss.execute(rhs);
 
 	if (0 <= rhs[0] && rhs[0] <= 1 && 0 <= rhs[1] && rhs[1] <= 1 && rhs[0] + rhs[1] <= 1)

GeoLib/Triangle.cpp

@@ -98,7 +98,7 @@ bool Triangle::containsPoint (const double *pnt) const
 			mat(1,1) = c[2] - a[2];
 			double y[2] = {pnt[1]-a[1], pnt[2]-a[2]};
 
-			MathLib::GaussAlgorithm gauss (mat);
+			MathLib::GaussAlgorithm<MathLib::DenseMatrix<double>, double*> gauss (mat);
 			gauss.execute (y);
 
 			if (-delta <= y[0] && y[0] <= upper && -delta <= y[1] && y[1] <= upper
@@ -125,7 +125,7 @@ bool Triangle::containsPoint (const double *pnt) const
 			mat(1,1) = c[2] - a[2];
 			double y[2] = {pnt[0]-a[0], pnt[2]-a[2]};
 
-			MathLib::GaussAlgorithm gauss (mat);
+			MathLib::GaussAlgorithm<MathLib::DenseMatrix<double>, double*> gauss (mat);
 			gauss.execute (y);
 
 			if (-delta <= y[0] && y[0] <= upper && -delta <= y[1] && y[1] <= upper && y[0] + y[1] <= upper) {
@@ -146,7 +146,7 @@ bool Triangle::containsPoint (const double *pnt) const
 	mat(1,1) = c[1] - a[1];
 	double y[2] = {pnt[0]-a[0], pnt[1]-a[1]};
 
-	MathLib::GaussAlgorithm gauss (mat);
+	MathLib::GaussAlgorithm<MathLib::DenseMatrix<double>, double*> gauss (mat);
 	gauss.execute (y);
 
 	// check if the solution fulfills the third equation
@@ -175,7 +175,7 @@ bool Triangle::containsPoint2D (const double *pnt) const
 	mat(1,1) = c[1] - a[1];
 	double y[2] = {pnt[0]-a[0], pnt[1]-a[1]};
 
-	MathLib::GaussAlgorithm gauss (mat);
+	MathLib::GaussAlgorithm<MathLib::DenseMatrix<double>, double*> gauss (mat);
 	gauss.execute (y);
 
 	const double delta (std::numeric_limits<double>::epsilon());
@@ -207,7 +207,7 @@ void getPlaneCoefficients(Triangle const& tri, double c[3])
 	c[1] = p1[2];
 	c[2] = p2[2];
 
-	MathLib::GaussAlgorithm gauss (mat);
+	MathLib::GaussAlgorithm<MathLib::DenseMatrix<double>, double*> gauss (mat);
 	gauss.execute (c);
 }
 

MathLib/LinAlg/Dense/DenseMatrix.h

@@ -30,6 +30,10 @@ namespace MathLib {
  */
 template <typename FP_TYPE, typename IDX_TYPE = std::size_t> class DenseMatrix
 {
+public:
+	typedef FP_TYPE FP_T;
+	typedef IDX_TYPE IDX_T;
+
 public:
    DenseMatrix (IDX_TYPE rows, IDX_TYPE cols);
    DenseMatrix (IDX_TYPE rows, IDX_TYPE cols, const FP_TYPE& val);

MathLib/LinAlg/Dense/DenseVector.h

@@ -29,6 +29,9 @@ namespace MathLib
 template <typename T>
 class DenseVector : public std::valarray<T>
 {
+public:
+	typedef T FP_T;
+
 public:
 	using std::valarray<T>::operator=;
 	using std::valarray<T>::operator[];

MathLib/LinAlg/Dense/GlobalDenseMatrix.h

@@ -26,6 +26,9 @@ namespace MathLib
 template<typename FP_TYPE, typename IDX_TYPE = std::size_t>
 class GlobalDenseMatrix: public DenseMatrix<FP_TYPE, IDX_TYPE>
 {
+public:
+	typedef FP_TYPE FP_T;
+
 public:
 	GlobalDenseMatrix (IDX_TYPE rows, IDX_TYPE cols);
 	GlobalDenseMatrix (IDX_TYPE rows, IDX_TYPE cols, const FP_TYPE& val);

MathLib/LinAlg/Solvers/DenseDirectLinearSolver.h

-/**
- * \file
- * \author Thomas Fischer
- * \date   2011-01-07
- * \brief  Definition of the DenseDirectLinearSolver 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 DENSEDIRECTLINEARSOLVER_H_
-#define DENSEDIRECTLINEARSOLVER_H_
-
-#include "DirectLinearSolver.h"
-
-namespace MathLib {
-
-class DenseDirectLinearSolver: public MathLib::DirectLinearSolver {
-public:
-	DenseDirectLinearSolver() {};
-	virtual ~DenseDirectLinearSolver() {};
-};
-
-}
-
-#endif /* DENSEDIRECTLINEARSOLVER_H_ */

MathLib/LinAlg/Solvers/DirectLinearSolver.h

-/**
- * \file
- * \author Thomas Fischer
- * \date   2011-01-07
- * \brief  Definition of the DirectLinearSolver 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 DIRECTLINEARSOLVER_H_
-#define DIRECTLINEARSOLVER_H_
-
-#include "LinearSolver.h"
-
-namespace MathLib {
-
-class DirectLinearSolver : public MathLib::LinearSolver
-{
-public:
-	DirectLinearSolver() {};
-	virtual ~DirectLinearSolver() {};
-};
-
-}
-
-#endif /* DIRECTLINEARSOLVER_H_ */

MathLib/LinAlg/Solvers/GaussAlgorithm.h

@@ -17,11 +17,13 @@
 
 #include <cstddef>
 #include "../Dense/DenseMatrix.h"
-#include "DenseDirectLinearSolver.h"
 #include "TriangularSolve.h"
 
 namespace MathLib {
 
+template <typename MAT_T, typename VEC_T>
+class GaussAlgorithm;
+
 /**
  * This is a class for the direct solution of (dense) systems of
  * linear equations, \f$A x = b\f$. During the construction of
@@ -30,7 +32,13 @@ namespace MathLib {
  * the entries of A change! The solution for a specific
  * right hand side is computed by the method execute().
  */
-class GaussAlgorithm : public MathLib::DenseDirectLinearSolver {
+template <typename MAT_T>
+class GaussAlgorithm <MAT_T, typename MAT_T::FP_T*>
+{
+public:
+	typedef typename MAT_T::FP_T FP_T;
+	typedef typename MAT_T::IDX_T IDX_T;
+
 public:
 	/**
 	 * A direct solver for the (dense) linear system \f$A x = b\f$.
@@ -41,7 +49,7 @@ public:
 	 * Attention: the given matrix will be destroyed!
 	 * @return a object of type GaussAlgorithm
 	 */
-	GaussAlgorithm(DenseMatrix<double> &A);
+	GaussAlgorithm(MAT_T &A);
 	/**
 	 * destructor, deletes the permutation
 	 */
@@ -52,7 +60,7 @@ public:
 	 * using forward solve and backward solve
 	 * @param b at the beginning the right hand side, at the end the solution
 	 */
-	void execute (double *b) const;
+	void execute (FP_T* b) const;
 
 private:
 	/**
@@ -60,22 +68,24 @@ private:
 	 * row permutations of the LU factorization
 	 * @param b the entries of the vector b are permuted
 	 */
-	void permuteRHS (double* b) const;
+	void permuteRHS (FP_T* b) const;
 
 	/**
 	 * a reference to the matrix
 	 */
-	DenseMatrix<double>& _mat;
+	MAT_T& _mat;
 	/**
 	 * the size of the matrix
 	 */
-	std::size_t _n;
+	IDX_T _n;
 	/**
 	 * the permutation of the rows
 	 */
-	std::size_t* _perm;
+	IDX_T* _perm;
 };
 
 } // end namespace MathLib
 
+#include "GaussAlgorithm.tpp"
+
 #endif /* GAUSSALGORITHM_H_ */

MathLib/LinAlg/Solvers/GaussAlgorithm.cpp → MathLib/LinAlg/Solvers/GaussAlgorithm.tpp

@@ -17,18 +17,19 @@
 
 namespace MathLib {
 
-GaussAlgorithm::GaussAlgorithm (DenseMatrix <double> &A) :
-	_mat (A), _n(_mat.getNRows()), _perm (new size_t [_n])
+template <typename MAT_T>
+GaussAlgorithm<MAT_T, typename MAT_T::FP_T*>::GaussAlgorithm (MAT_T &A) :
+	_mat (A), _n(_mat.getNRows()), _perm (new IDX_T [_n])
 {
-	size_t k, i, j, nr (_mat.getNRows()), nc(_mat.getNCols());
-	double l;
+	IDX_T k, i, j, nr (_mat.getNRows()), nc(_mat.getNCols());
+	FP_T l;
 
 	for (k=0; k<nc; k++) {
 		// search pivot
-		double t = fabs(_mat(k, k));
+		FP_T t = std::abs(_mat(k, k));
 		_perm[k] = k;
 		for (i=k+1; i<nr; i++) {
-			if (fabs(_mat(i,k)) > t) {
+			if (std::abs(_mat(i,k)) > t) {
 				t = _mat(i,k);
 				_perm[k] = i;
 			}
@@ -36,7 +37,8 @@ GaussAlgorithm::GaussAlgorithm (DenseMatrix <double> &A) :
 
 		// exchange rows
 		if (_perm[k] != k) {
-			for (j=0; j<nc; j++) std::swap (_mat(_perm[k],j), _mat(k,j));
+			for (j=0; j<nc; j++)
+				std::swap (_mat(_perm[k],j), _mat(k,j));
 		}
 
 		// eliminate
@@ -50,21 +52,24 @@ GaussAlgorithm::GaussAlgorithm (DenseMatrix <double> &A) :
 	}
 }
 
-GaussAlgorithm::~GaussAlgorithm()
+template <typename MAT_T>
+GaussAlgorithm<MAT_T, typename MAT_T::FP_T*>::~GaussAlgorithm()
 {
 	delete [] _perm;
 }
 
-void GaussAlgorithm::execute (double *b) const
+template <typename MAT_T>
+void GaussAlgorithm<MAT_T, typename MAT_T::FP_T*>::execute (typename MAT_T::FP_T *b) const
 {
 	permuteRHS (b);
 	forwardSolve (_mat, b); // L z = b, b will be overwritten by z
 	backwardSolve (_mat, b); // U x = z, b (z) will be overwritten by x
 }
 
-void GaussAlgorithm::permuteRHS (double* b) const
+template <typename MAT_T>
+void GaussAlgorithm<MAT_T, typename MAT_T::FP_T*>::permuteRHS (typename MAT_T::FP_T* b) const
 {
-	for (size_t i=0; i<_n; i++) {
+	for (IDX_T i=0; i<_n; i++) {
 		if (_perm[i] != i) std::swap(b[i], b[_perm[i]]);
 	}
 }

MathLib/LinAlg/Solvers/LinearSolver.h

-/**
- * \file
- * \author Thomas Fischer
- * \date   2011-01-07
- * \brief  Definition of the LinearSolver 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 LINEARSOLVER_H_
-#define LINEARSOLVER_H_
-
-namespace MathLib {
-
-/**
- * Base class for all linear solver classes.
- */
-class LinearSolver {
-public:
-	LinearSolver() {};
-	virtual ~LinearSolver() {};
-};
-
-}
-
-#endif /* LINEARSOLVER_H_ */

MathLib/LinAlg/Solvers/TriangularSolve.h

@@ -25,7 +25,8 @@ namespace MathLib {
  * @param L the lower triangular matrix
  * @param b at beginning the right hand side vector, at the end the solution vector
  */
-void forwardSolve (const DenseMatrix <double> &L, double* b);
+template <typename FP_T>
+void forwardSolve (const DenseMatrix <FP_T> &L, FP_T* b);
 
 /**
  * solves the \f$n \times n\f$ triangular linear system \f$U \cdot x=y\f$,
@@ -33,7 +34,8 @@ void forwardSolve (const DenseMatrix <double> &L, double* b);
  * @param U upper triangular matrix
  * @param y at beginning the right hand side, at the end the solution
  */
-void backwardSolve (const DenseMatrix <double> &U, double* y);
+template <typename FP_T>
+void backwardSolve (const DenseMatrix <FP_T> &U, FP_T* y);
 
 // backwardSolve mat * x = y, mat ... upper triangular matrix
 /**
@@ -43,8 +45,11 @@ void backwardSolve (const DenseMatrix <double> &U, double* y);
  * @param x the solution of the system of linear equations
  * @param b the right hand side
  */
-void backwardSolve ( DenseMatrix<double> const& mat, double* x, double* b);
+template <typename FP_T>
+void backwardSolve ( DenseMatrix<FP_T> const& mat, FP_T* x, FP_T* b);
 
 } // end namespace MathLib
 
+#include "TriangularSolve.tpp"
+
 #endif /* TRIANGULARSOLVE_H_ */

MathLib/LinAlg/Solvers/TriangularSolve.cpp → MathLib/LinAlg/Solvers/TriangularSolve.tpp

@@ -12,44 +12,48 @@
  *
  */
 
-#include "TriangularSolve.h"
-
 namespace MathLib {
 
-void forwardSolve (const DenseMatrix <double> &L, double* b)
+template <typename FP_T>
+void forwardSolve (const DenseMatrix <FP_T> &L, FP_T* b)
 {
-	size_t m (L.getNRows());
-	double t;
+	typedef typename DenseMatrix<FP_T>::IDX_T IDX_T;
+	IDX_T m (L.getNRows());
+	FP_T t;
 
-	for (size_t r=0; r<m; r++) {
+	for (IDX_T r=0; r<m; r++) {
 		t = 0.0;
-		for (size_t c=0; c<r; c++) {
+		for (IDX_T c=0; c<r; c++) {
 			t += L(r,c)*b[c];
 		}
 		b[r] = b[r]-t;
 	}
 }
 
-void backwardSolve (const DenseMatrix <double> &mat, double* b)
+template <typename FP_T>
+void backwardSolve (const DenseMatrix <FP_T> &mat, FP_T* b)
 {
-	double t;
-	size_t m (mat.getNRows()), n(mat.getNCols());
+	FP_T t;
+	typedef typename DenseMatrix<FP_T>::IDX_T IDX_T;
+	IDX_T m (mat.getNRows()), n(mat.getNCols());
 	for (int r=m-1; r>=0; r--) {
 		t = 0.0;
-		for (size_t c=r+1; c<n; c++) {
+		for (IDX_T c=r+1; c<n; c++) {
 			t += mat(r,c)*b[c];
 		}
 		b[r] = (b[r]-t) / mat(r,r);
 	}
 }
 
-void backwardSolve ( DenseMatrix<double> const& mat, double* x, double* b)
+template <typename FP_T>
+void backwardSolve ( DenseMatrix<FP_T> const& mat, FP_T* x, FP_T* b)
 {
-	size_t n_cols (mat.getNCols());
+	typedef typename DenseMatrix<FP_T>::IDX_T IDX_T;
+	IDX_T n_cols (mat.getNCols());
 	for (int r = (n_cols - 1); r >= 0; r--) {
-		double t = 0.0;
+		FP_T t = 0.0;
 
-		for (size_t c = r+1; c < n_cols; c++) {
+		for (IDX_T c = r+1; c < n_cols; c++) {
 			t += mat(r,c) * b[c];
 		}
 		x[r] = (b[r] - t) / mat(r, r);

MathLib/LinAlg/Sparse/CRSMatrix.h

@@ -32,6 +32,9 @@ namespace MathLib {
 template<typename FP_TYPE, typename IDX_TYPE>
 class CRSMatrix: public SparseMatrixBase<FP_TYPE, IDX_TYPE>
 {
+public:
+	typedef FP_TYPE FP_T;
+
 public:
 	explicit CRSMatrix(std::string const &fname) :
 		SparseMatrixBase<FP_TYPE, IDX_TYPE>(),
@@ -164,7 +167,7 @@ public:
      * @param col the column number
      * @return The corresponding matrix entry or 0.0.
      */
-	double getValue(IDX_TYPE row, IDX_TYPE col)
+	FP_TYPE getValue(IDX_TYPE row, IDX_TYPE col)
 	{
 		assert(0 <= row && row < this->_n_rows);
 

MathLib/LinAlg/Sparse/CRSSymMatrix.h

@@ -19,6 +19,9 @@
 
 template<class T> class CRSSymMatrix : public CRSMatrix<T>
 {
+public:
+	typedef T FP_T;
+
 public:
 	CRSSymMatrix(std::string const &fname)
 	: CRSMatrix<T> (fname)

Tests/MathLib/TestDenseGaussAlgorithm.cpp

+/**
+ * @file TestDenseGaussAlgorithm.cpp
+ * @author Thomas Fischer
+ * @date Jun 17, 2013
+ * @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/LICENSE.txt
+ */
+
+#include <gtest/gtest.h>
+
+#include <cstdlib>
+#include <ctime>
+#include <limits>
+#include <algorithm>
+
+#include "LinAlg/Dense/DenseMatrix.h"
+#include "LinAlg/Solvers/GaussAlgorithm.h"
+
+TEST(MathLib, DenseGaussAlgorithm)
+{
+	std::size_t n_rows(50);
+	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(NULL)) );
+	// ** 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
+	double *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);
+
+	MathLib::GaussAlgorithm<MathLib::DenseMatrix<double, std::size_t>, double*> gauss(mat);
+
+	// solve with b0 as right hand side
+	gauss.execute(b0);
+	for (std::size_t i(0); i<n_rows; i++) {
+		ASSERT_NEAR(b0[i], 0.0, 1e5 * std::numeric_limits<double>::epsilon());
+	}
+
+	// solve with b1 as right hand side
+	gauss.execute(b1);
+	for (std::size_t i(0); i<n_rows; i++) {
+		ASSERT_NEAR(b1[i], 1.0, 1e5 * std::numeric_limits<double>::epsilon());
+	}
+
+	// solve with b2 as right hand side
+	gauss.execute(b2);
+	for (std::size_t i(0); i<n_rows; i++) {
+		ASSERT_NEAR(fabs(b2[i]-x[i])/fabs(x[i]), 0.0, 1e5*std::numeric_limits<double>::epsilon());
+	}
+
+	delete [] x;
+	delete [] b0;
+	delete [] b1;
+	delete [] b2;
+}