Commit 6454c5ce authored by Nils Rummler's avatar Nils Rummler
Browse files

moved Bandmatrix solver to own file in Math namespace

parent cfcc03d0
/*
* ----------------------------------------------------------------
*
* ITA core libs
* (c) Copyright Institute of Technical Acoustics (ITA)
* RWTH Aachen University, Germany, 2015-2020
*
* ----------------------------------------------------------------
* ____ __________ _______
* // / //__ ___/ // _ |
* // / // / // /_| |
* // / // / // ___ |
* //__/ //__/ //__/ |__|
*
* ----------------------------------------------------------------
*
*/
#ifndef IW_ITA_BASE_LINEAR_ALGEBRA
#define IW_ITA_BASE_LINEAR_ALGEBRA
#include <ITABaseDefinitions.h>
#include <vector>
namespace ITABase
{
namespace Math
{
//! Solves a Tridiagonal Bandmatrix with Thomas algorithm. Expects the 3 diagonals (lower,middel and upper) and the excitation vector.
ITA_BASE_API std::vector<float> BandmatrixSolver(const std::vector<float> vdLower, const std::vector<float> vdMiddle, const std::vector<float> vdUpper, const std::vector<float> vdExcitation);
}
}
#endif // IW_ITA_BASE_LINEAR_ALGEBRA
# $Id:$
set( RelativeDir "include/ITABase/Math" )
set( RelativeSourceGroup "Header Files\\ITABase\\Math" )
set( SubDirs )
set( DirFiles
LinearAlgebra.h
)
set( DirFiles_SourceGroup "${RelativeSourceGroup}" )
set( LocalSourceGroupFiles )
foreach( File ${DirFiles} )
list( APPEND LocalSourceGroupFiles "${RelativeDir}/${File}" )
list( APPEND ProjectSources "${RelativeDir}/${File}" )
endforeach()
source_group( ${DirFiles_SourceGroup} FILES ${LocalSourceGroupFiles} )
set( SubDirFiles "" )
foreach( Dir ${SubDirs} )
list( APPEND SubDirFiles "${RelativeDir}/${Dir}/_SourceFiles.cmake" )
endforeach()
foreach( SubDirFile ${SubDirFiles} )
include( ${SubDirFile} )
endforeach()
......@@ -39,7 +39,7 @@ namespace ITABase
class ITA_BASE_API CPiecewisePolynomial
{
public:
CPiecewisePolynomial( const std::vector<float>& vdBreakPoints, const std::vector<float>& vdCoeffs, const int iOrder );
CPiecewisePolynomial(const std::vector<float>& vdBreakPoints, const std::vector<float>& vdCoeffs, const int iOrder);
//! Returns the polynomial order
inline int Order() const { return iOrder; };
......@@ -88,8 +88,6 @@ namespace ITABase
//! Uses cubic splines to create a piecewise polynomial (order 3) for the given data pairs. vdSupportingPoints need to be sorted in ascending order
ITA_BASE_API CPiecewisePolynomial CubicSpline(const std::vector<float>& vdSupportingPoints, const std::vector<float>& vdDataPoints) { return Spline(vdSupportingPoints, vdDataPoints, 3); };
//! Solves a Tridiagonal Bandmatrix with Thomas algorithm. Expects the 3 diagonals (lower,middel and upper) and the excitation vector.
std::vector<float> BandmatrixSolver(const std::vector<float> vdLower,const std::vector<float> vdMiddle,const std::vector<float> vdUpper,const std::vector<float> vdExcitation);
}
#endif // IW_ITA_BASE_SPLINE
......@@ -2,7 +2,7 @@
set( RelativeDir "include/ITABase" )
set( RelativeSourceGroup "Header Files\\ITABase" )
set( SubDirs )
set( SubDirs Math )
set( DirFiles
ITAProgress.h
......
#include <ITAException.h>
#include <ITABase/Math/LinearAlgebra.h>
#include <string>
using namespace ITABase;
std::vector<float> Math::BandmatrixSolver(const std::vector<float> vdLower,const std::vector<float> vdMiddle, const std::vector<float> vdUpper, const std::vector<float> vdExcitation)
{
//solving bandmatrix with "numerical recipees in C" page 51
if (vdMiddle[0] == 0)
ITA_EXCEPT_INVALID_PARAMETER("Matrix needs to be semi positiv definit.");
if (vdMiddle.size() != vdExcitation.size() || vdLower.size() != vdUpper.size() || vdExcitation.size() != vdLower.size() + 1)
ITA_EXCEPT_INVALID_PARAMETER("Dimensions of arguments missmatch.");
int n = vdMiddle.size();
std::vector<float> vdGamma(n);
float fBeta = vdMiddle[0];
std::vector<float> vdB(n); //solution of bandmatrix
vdB[0] = vdExcitation[0] / vdMiddle[0];
for (int j = 1; j < n; j++) {
vdGamma[j] = vdUpper[j - 1] / fBeta;
fBeta = vdMiddle[j] - vdLower[j - 1] * vdGamma[j];
vdB[j] = (vdExcitation[j] - vdLower[j - 1] * vdB[j-1]) / fBeta;
}
for (int j = n-2; j >= 0; j--) {//backsubstitution
vdB[j] -= vdGamma[j+1] * vdB[j + 1];
}
return vdB;
}
# $Id:$
set( RelativeDir "src/ITABase/Math" )
set( RelativeSourceGroup "Source Files\\ITABase\\Math" )
set( SubDirs )
set( DirFiles
LinearAlgebra.cpp
)
set( DirFiles_SourceGroup "${RelativeSourceGroup}" )
set( LocalSourceGroupFiles )
foreach( File ${DirFiles} )
list( APPEND LocalSourceGroupFiles "${RelativeDir}/${File}" )
list( APPEND ProjectSources "${RelativeDir}/${File}" )
endforeach()
source_group( ${DirFiles_SourceGroup} FILES ${LocalSourceGroupFiles} )
set( SubDirFiles "" )
foreach( Dir ${SubDirs} )
list( APPEND SubDirFiles "${RelativeDir}/${Dir}/_SourceFiles.cmake" )
endforeach()
foreach( SubDirFile ${SubDirFiles} )
include( ${SubDirFile} )
endforeach()
......@@ -4,6 +4,7 @@
#include <string>
#include <algorithm>
#include <cmath>
#include <ITABase/Math/LinearAlgebra.h>
using namespace ITABase;
......@@ -131,7 +132,7 @@ CPiecewisePolynomial ITABase::Spline(const std::vector<float>& vdSupportingPoint
vdR[i] = 3 *( (vdDataPoints[i + 2] - vdDataPoints[i + 1]) / vdH[i + 1] - (vdDataPoints[i + 1] - vdDataPoints[i]) / vdH[i]);
//solve bandmatrix and get the other coefficients
std::vector<float> vdB = BandmatrixSolver(vdLower,vdMiddle,vdUpper,vdR);
std::vector<float> vdB = Math::BandmatrixSolver(vdLower,vdMiddle,vdUpper,vdR);
vdB.insert(vdB.begin(), 0.0);
std::vector<float> vdA(iPoints - 1);
......@@ -161,30 +162,4 @@ CPiecewisePolynomial ITABase::Spline(const std::vector<float>& vdSupportingPoint
return CPiecewisePolynomial(vdSupportingPoints, vdCoefficients, iPolynomialOrder);
}
std::vector<float> ITABase::BandmatrixSolver(const std::vector<float> vdLower,const std::vector<float> vdMiddle, const std::vector<float> vdUpper, const std::vector<float> vdExcitation)
{
//solving bandmatrix with "numerical recipees in C" page 51
if (vdMiddle[0] == 0)
ITA_EXCEPT_INVALID_PARAMETER("Matrix needs to be semi positiv definit.");
if (vdMiddle.size() != vdExcitation.size() || vdLower.size() != vdUpper.size() || vdExcitation.size() != vdLower.size() + 1)
ITA_EXCEPT_INVALID_PARAMETER("Dimensions of arguments missmatch.");
int n = vdMiddle.size();
std::vector<float> vdGamma(n);
float fBeta = vdMiddle[0];
std::vector<float> vdB(n); //solution of bandmatrix
vdB[0] = vdExcitation[0] / vdMiddle[0];
for (int j = 1; j < n; j++) {
vdGamma[j] = vdUpper[j - 1] / fBeta;
fBeta = vdMiddle[j] - vdLower[j - 1] * vdGamma[j];
vdB[j] = (vdExcitation[j] - vdLower[j - 1] * vdB[j-1]) / fBeta;
}
for (int j = n-2; j >= 0; j--) {//backsubstitution
vdB[j] -= vdGamma[j+1] * vdB[j + 1];
}
return vdB;
}
\ No newline at end of file
......@@ -2,7 +2,7 @@
set( RelativeDir "src/ITABase" )
set( RelativeSourceGroup "Source Files\\ITABase" )
set( SubDirs )
set( SubDirs Math )
set( DirFiles
ITAProgress.cpp
......
Markdown is supported
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