Skip to content
Snippets Groups Projects
Commit 8cc72a18 authored by Tom Fischer's avatar Tom Fischer
Browse files

[MaL] Remove GaussAlgorithm implementation.

parent 72d28467
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
MathLib/CMakeLists.txt
+ 0
− 1
View file @ 8cc72a18
......@@ -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)
......
MathLib/LinAlg/Solvers/GaussAlgorithm-impl.h deleted 100644 → 0
+ 0
− 113
View file @ 72d28467
/**
* \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
MathLib/LinAlg/Solvers/GaussAlgorithm.h deleted 100644 → 0
+ 0
− 101
View file @ 72d28467
/**
* \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"
MathLib/LinAlg/Solvers/TriangularSolve-impl.h deleted 100644 → 0
+ 0
− 64
View file @ 72d28467
/**
* \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
MathLib/LinAlg/Solvers/TriangularSolve.h deleted 100644 → 0
+ 0
− 52
View file @ 72d28467
/**
* \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"
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment