diff --git a/Applications/ApplicationsLib/UncoupledProcessesTimeLoop.cpp b/Applications/ApplicationsLib/UncoupledProcessesTimeLoop.cpp
index 27f0f01802be099aef2925af02ce0aab500270bf..ab3c0efedb63b887c116cc628fb6ace0f1869426 100644
--- a/Applications/ApplicationsLib/UncoupledProcessesTimeLoop.cpp
+++ b/Applications/ApplicationsLib/UncoupledProcessesTimeLoop.cpp
@@ -88,7 +88,7 @@ void UncoupledProcessesTimeLoop::setInitialConditions(
auto& x0 = *_process_solutions[pcs_idx];
pcs.setInitialConditions(x0);
- MathLib::BLAS::finalizeAssembly(x0);
+ MathLib::LinAlg::finalizeAssembly(x0);
time_disc.setInitialState(t0, x0); // push IC
diff --git a/MathLib/LinAlg/BLAS.cpp b/MathLib/LinAlg/LinAlg.cpp
similarity index 97%
rename from MathLib/LinAlg/BLAS.cpp
rename to MathLib/LinAlg/LinAlg.cpp
index 84d2376fbec6415e18679dc1c8059487fc409882..08527728ddc1ee537541cc3d899667e4ddd55a75 100644
--- a/MathLib/LinAlg/BLAS.cpp
+++ b/MathLib/LinAlg/LinAlg.cpp
@@ -7,9 +7,9 @@
*
*/
-#include "BLAS.h"
+#include "LinAlg.h"
-// TODO reorder BLAS function signatures?
+// TODO reorder LinAlg function signatures?
// Dense Eigen matrices and vectors ////////////////////////////////////////
#ifdef OGS_USE_EIGEN
@@ -18,7 +18,7 @@
namespace MathLib
{
-namespace BLAS
+namespace LinAlg
{
// Explicit specialization
@@ -65,7 +65,7 @@ double normMax(Eigen::VectorXd const& x)
#include "MathLib/LinAlg/PETSc/PETScVector.h"
#include "MathLib/LinAlg/PETSc/PETScMatrix.h"
-namespace MathLib { namespace BLAS
+namespace MathLib { namespace LinAlg
{
// Vector
@@ -220,7 +220,7 @@ void finalizeAssembly(PETScVector& x)
#include "MathLib/LinAlg/Eigen/EigenVector.h"
#include "MathLib/LinAlg/Eigen/EigenMatrix.h"
-namespace MathLib { namespace BLAS
+namespace MathLib { namespace LinAlg
{
// Vector
@@ -347,7 +347,7 @@ void finalizeAssembly(EigenMatrix& x)
x.getRawMatrix().makeCompressed();
}
-} // namespace BLAS
+} // namespace LinAlg
} // namespace MathLib
diff --git a/MathLib/LinAlg/BLAS.h b/MathLib/LinAlg/LinAlg.h
similarity index 96%
rename from MathLib/LinAlg/BLAS.h
rename to MathLib/LinAlg/LinAlg.h
index 5c261401c3b2d1901ea1d706fb83b30025b568db..8bf87bfbd9791c17899f9602038c1a4711edcab0 100644
--- a/MathLib/LinAlg/BLAS.h
+++ b/MathLib/LinAlg/LinAlg.h
@@ -7,16 +7,16 @@
*
*/
-#ifndef MATHLIB_BLAS_H
-#define MATHLIB_BLAS_H
+#ifndef MATHLIB_LINALG_H
+#define MATHLIB_LINALG_H
-#include<cassert>
+#include <cassert>
namespace MathLib
{
-/*! \namespace MathLib::BLAS
+/*! \namespace MathLib::LinAlg
* Some general linear algebra functionality.
*
* By using the provided functions linear algebra capabilities can be
@@ -26,7 +26,7 @@ namespace MathLib
* For documentation, refer to that of the templated method. All specializations
* or overload must behave in the same way.
*/
-namespace BLAS
+namespace LinAlg
{
// Matrix or Vector
@@ -128,7 +128,7 @@ namespace MathLib {
class PETScMatrix;
class PETScVector;
-namespace BLAS
+namespace LinAlg
{
// Vector
@@ -186,7 +186,7 @@ namespace MathLib {
class EigenMatrix;
class EigenVector;
-namespace BLAS
+namespace LinAlg
{
// Vector
@@ -232,10 +232,10 @@ void matMultAdd(EigenMatrix const& A, EigenVector const& v1,
void finalizeAssembly(EigenMatrix& A);
-} // namespace BLAS
+} // namespace LinAlg
} // namespace MathLib
#endif
-#endif // MATHLIB_BLAS_H
+#endif // MATHLIB_LINALG_H
diff --git a/NumLib/DOF/SimpleMatrixVectorProvider-impl.h b/NumLib/DOF/SimpleMatrixVectorProvider-impl.h
index 8e0672aec4236f1429d0299f4ce2633c33ed7642..30c2a742b9abba34d0c9a9c70a2ff1d2495718fb 100644
--- a/NumLib/DOF/SimpleMatrixVectorProvider-impl.h
+++ b/NumLib/DOF/SimpleMatrixVectorProvider-impl.h
@@ -11,7 +11,7 @@
#include <logog/include/logog.hpp>
#include "BaseLib/Error.h"
-#include "MathLib/LinAlg/BLAS.h"
+#include "MathLib/LinAlg/LinAlg.h"
#include "MathLib/LinAlg/MatrixVectorTraits.h"
#include "SimpleMatrixVectorProvider.h"
@@ -140,7 +140,7 @@ getMatrix(Matrix const& A)
std::size_t id = 0u;
auto const& res = getMatrix_<false>(id, A);
if (!res.second) // no new object has been created
- BLAS::copy(A, *res.first);
+ LinAlg::copy(A, *res.first);
return *res.first;
}
@@ -151,7 +151,7 @@ getMatrix(Matrix const& A, std::size_t& id)
{
auto const& res = getMatrix_<true>(id, A);
if (!res.second) // no new object has been created
- BLAS::copy(A, *res.first);
+ LinAlg::copy(A, *res.first);
return *res.first;
}
@@ -222,7 +222,7 @@ getVector(Vector const& x)
std::size_t id = 0u;
auto const& res = getVector_<false>(id, x);
if (!res.second) // no new object has been created
- BLAS::copy(x, *res.first);
+ LinAlg::copy(x, *res.first);
return *res.first;
}
@@ -233,7 +233,7 @@ getVector(Vector const& x, std::size_t& id)
{
auto const& res = getVector_<true>(id, x);
if (!res.second) // no new object has been created
- BLAS::copy(x, *res.first);
+ LinAlg::copy(x, *res.first);
return *res.first;
}
diff --git a/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator-impl.h b/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator-impl.h
index b95f39cb29127f7ebb503293b209ea1ce25a664a..0b3121af7b01b56849f229c4dae196a4e38d6a9d 100644
--- a/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator-impl.h
+++ b/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator-impl.h
@@ -12,7 +12,7 @@
#include <logog/include/logog.hpp>
#include <Eigen/Core>
-#include "MathLib/LinAlg/BLAS.h"
+#include "MathLib/LinAlg/LinAlg.h"
#include "NumLib/Assembler/SerialExecutor.h"
#include "MathLib/LinAlg/MatrixVectorTraits.h"
#include "NumLib/Function/Interpolation.h"
@@ -40,7 +40,7 @@ extrapolate(LocalAssemblers const& local_assemblers, PropertyTag const property)
NumLib::SerialExecutor::executeMemberDereferenced(
*this, &Self::extrapolateElement, local_assemblers, property, *counts);
- MathLib::BLAS::componentwiseDivide(_nodal_values, _nodal_values, *counts);
+ MathLib::LinAlg::componentwiseDivide(_nodal_values, _nodal_values, *counts);
}
template<typename PropertyTag, typename LocalAssembler>
diff --git a/NumLib/ODESolver/MatrixTranslator.cpp b/NumLib/ODESolver/MatrixTranslator.cpp
index bedfcea39240cfcb2c1ad6ad2f9898c7ac674154..b79d9976a68057c75195aa6c75c970c29954e2a3 100644
--- a/NumLib/ODESolver/MatrixTranslator.cpp
+++ b/NumLib/ODESolver/MatrixTranslator.cpp
@@ -9,7 +9,7 @@
#include "MatrixTranslator.h"
-#include "MathLib/LinAlg/BLAS.h"
+#include "MathLib/LinAlg/LinAlg.h"
namespace NumLib
{
@@ -17,26 +17,26 @@ void MatrixTranslatorGeneral<ODESystemTag::FirstOrderImplicitQuasilinear>::
computeA(GlobalMatrix const& M, GlobalMatrix const& K,
GlobalMatrix& A) const
{
- namespace BLAS = MathLib::BLAS;
+ namespace LinAlg = MathLib::LinAlg;
auto const dxdot_dx = _time_disc.getNewXWeight();
// A = M * dxdot_dx + K
- BLAS::copy(M, A);
- BLAS::aypx(A, dxdot_dx, K);
+ LinAlg::copy(M, A);
+ LinAlg::aypx(A, dxdot_dx, K);
}
void MatrixTranslatorGeneral<ODESystemTag::FirstOrderImplicitQuasilinear>::
computeRhs(const GlobalMatrix& M, const GlobalMatrix& /*K*/,
const GlobalVector& b, GlobalVector& rhs) const
{
- namespace BLAS = MathLib::BLAS;
+ namespace LinAlg = MathLib::LinAlg;
auto& tmp = MathLib::GlobalVectorProvider<GlobalVector>::provider.getVector(_tmp_id);
_time_disc.getWeightedOldX(tmp);
// rhs = M * weighted_old_x + b
- BLAS::matMultAdd(M, tmp, b, rhs);
+ LinAlg::matMultAdd(M, tmp, b, rhs);
MathLib::GlobalVectorProvider<GlobalVector>::provider.releaseVector(tmp);
}
@@ -46,43 +46,43 @@ void MatrixTranslatorGeneral<ODESystemTag::FirstOrderImplicitQuasilinear>::
GlobalVector const& b, GlobalVector const& x_new_timestep,
GlobalVector const& xdot, GlobalVector& res) const
{
- namespace BLAS = MathLib::BLAS;
+ namespace LinAlg = MathLib::LinAlg;
auto const& x_curr = _time_disc.getCurrentX(x_new_timestep);
// res = M * x_dot + K * x_curr - b
- BLAS::matMult(M, xdot, res); // the local vector x_dot seems to be
+ LinAlg::matMult(M, xdot, res); // the local vector x_dot seems to be
// necessary because of this multiplication
- BLAS::matMultAdd(K, x_curr, res, res);
- BLAS::axpy(res, -1.0, b);
+ LinAlg::matMultAdd(K, x_curr, res, res);
+ LinAlg::axpy(res, -1.0, b);
}
void MatrixTranslatorGeneral<ODESystemTag::FirstOrderImplicitQuasilinear>::
computeJacobian(GlobalMatrix const& Jac_in, GlobalMatrix& Jac_out) const
{
- namespace BLAS = MathLib::BLAS;
+ namespace LinAlg = MathLib::LinAlg;
- BLAS::copy(Jac_in, Jac_out);
+ LinAlg::copy(Jac_in, Jac_out);
}
void MatrixTranslatorForwardEuler<ODESystemTag::FirstOrderImplicitQuasilinear>::
computeA(GlobalMatrix const& M, GlobalMatrix const& /*K*/,
GlobalMatrix& A) const
{
- namespace BLAS = MathLib::BLAS;
+ namespace LinAlg = MathLib::LinAlg;
auto const dxdot_dx = _fwd_euler.getNewXWeight();
// A = M * dxdot_dx
- BLAS::copy(M, A);
- BLAS::scale(A, dxdot_dx);
+ LinAlg::copy(M, A);
+ LinAlg::scale(A, dxdot_dx);
}
void MatrixTranslatorForwardEuler<ODESystemTag::FirstOrderImplicitQuasilinear>::
computeRhs(const GlobalMatrix& M, const GlobalMatrix& K,
const GlobalVector& b, GlobalVector& rhs) const
{
- namespace BLAS = MathLib::BLAS;
+ namespace LinAlg = MathLib::LinAlg;
auto& tmp = MathLib::GlobalVectorProvider<GlobalVector>::provider.getVector(_tmp_id);
_fwd_euler.getWeightedOldX(tmp);
@@ -90,9 +90,9 @@ void MatrixTranslatorForwardEuler<ODESystemTag::FirstOrderImplicitQuasilinear>::
auto const& x_old = _fwd_euler.getXOld();
// rhs = b + M * weighted_old_x - K * x_old
- BLAS::matMult(K, x_old, rhs); // rhs = K * x_old
- BLAS::aypx(rhs, -1.0, b); // rhs = b - K * x_old
- BLAS::matMultAdd(M, tmp, rhs, rhs); // rhs += M * weighted_old_x
+ LinAlg::matMult(K, x_old, rhs); // rhs = K * x_old
+ LinAlg::aypx(rhs, -1.0, b); // rhs = b - K * x_old
+ LinAlg::matMultAdd(M, tmp, rhs, rhs); // rhs += M * weighted_old_x
MathLib::GlobalVectorProvider<GlobalVector>::provider.releaseVector(tmp);
}
@@ -102,22 +102,22 @@ void MatrixTranslatorForwardEuler<ODESystemTag::FirstOrderImplicitQuasilinear>::
GlobalVector const& b, GlobalVector const& x_new_timestep,
GlobalVector const& xdot, GlobalVector& res) const
{
- namespace BLAS = MathLib::BLAS;
+ namespace LinAlg = MathLib::LinAlg;
auto const& x_curr = _fwd_euler.getCurrentX(x_new_timestep);
// res = M * x_dot + K * x_curr - b
- BLAS::matMult(M, xdot, res);
- BLAS::matMultAdd(K, x_curr, res, res);
- BLAS::axpy(res, -1.0, b);
+ LinAlg::matMult(M, xdot, res);
+ LinAlg::matMultAdd(K, x_curr, res, res);
+ LinAlg::axpy(res, -1.0, b);
}
void MatrixTranslatorForwardEuler<ODESystemTag::FirstOrderImplicitQuasilinear>::
computeJacobian(GlobalMatrix const& Jac_in, GlobalMatrix& Jac_out) const
{
- namespace BLAS = MathLib::BLAS;
+ namespace LinAlg = MathLib::LinAlg;
- BLAS::copy(Jac_in, Jac_out);
+ LinAlg::copy(Jac_in, Jac_out);
}
void MatrixTranslatorCrankNicolson<
@@ -125,17 +125,17 @@ void MatrixTranslatorCrankNicolson<
computeA(GlobalMatrix const& M, GlobalMatrix const& K,
GlobalMatrix& A) const
{
- namespace BLAS = MathLib::BLAS;
+ namespace LinAlg = MathLib::LinAlg;
auto const dxdot_dx = _crank_nicolson.getNewXWeight();
auto const theta = _crank_nicolson.getTheta();
// A = theta * (M * dxdot_dx + K) + dxdot_dx * _M_bar
- BLAS::copy(M, A);
- BLAS::aypx(A, dxdot_dx, K); // A = M * dxdot_dx + K
+ LinAlg::copy(M, A);
+ LinAlg::aypx(A, dxdot_dx, K); // A = M * dxdot_dx + K
- BLAS::scale(A, theta); // A *= theta
- BLAS::axpy(A, dxdot_dx, _M_bar); // A += dxdot_dx * _M_bar
+ LinAlg::scale(A, theta); // A *= theta
+ LinAlg::axpy(A, dxdot_dx, _M_bar); // A += dxdot_dx * _M_bar
}
void MatrixTranslatorCrankNicolson<
@@ -143,7 +143,7 @@ void MatrixTranslatorCrankNicolson<
computeRhs(const GlobalMatrix& M, const GlobalMatrix& /*K*/,
const GlobalVector& b, GlobalVector& rhs) const
{
- namespace BLAS = MathLib::BLAS;
+ namespace LinAlg = MathLib::LinAlg;
auto& tmp = MathLib::GlobalVectorProvider<GlobalVector>::provider.getVector(_tmp_id);
_crank_nicolson.getWeightedOldX(tmp);
@@ -152,11 +152,11 @@ void MatrixTranslatorCrankNicolson<
// rhs = theta * (b + M * weighted_old_x) + _M_bar * weighted_old_x -
// _b_bar;
- BLAS::matMultAdd(M, tmp, b, rhs); // rhs = b + M * weighted_old_x
+ LinAlg::matMultAdd(M, tmp, b, rhs); // rhs = b + M * weighted_old_x
- BLAS::scale(rhs, theta); // rhs *= theta
- BLAS::matMultAdd(_M_bar, tmp, rhs, rhs); // rhs += _M_bar * weighted_old_x
- BLAS::axpy(rhs, -1.0, _b_bar); // rhs -= b
+ LinAlg::scale(rhs, theta); // rhs *= theta
+ LinAlg::matMultAdd(_M_bar, tmp, rhs, rhs); // rhs += _M_bar * weighted_old_x
+ LinAlg::axpy(rhs, -1.0, _b_bar); // rhs -= b
MathLib::GlobalVectorProvider<GlobalVector>::provider.releaseVector(tmp);
}
@@ -167,33 +167,33 @@ void MatrixTranslatorCrankNicolson<
GlobalVector const& b, GlobalVector const& x_new_timestep,
GlobalVector const& xdot, GlobalVector& res) const
{
- namespace BLAS = MathLib::BLAS;
+ namespace LinAlg = MathLib::LinAlg;
auto const& x_curr = _crank_nicolson.getCurrentX(x_new_timestep);
auto const theta = _crank_nicolson.getTheta();
// res = theta * (M * x_dot + K*x_curr - b) + _M_bar * x_dot + _b_bar
- BLAS::matMult(M, xdot, res); // res = M * x_dot
- BLAS::matMultAdd(K, x_curr, res, res); // res += K * x_curr
- BLAS::axpy(res, -1.0, b); // res = M * x_dot + K * x_curr - b
+ LinAlg::matMult(M, xdot, res); // res = M * x_dot
+ LinAlg::matMultAdd(K, x_curr, res, res); // res += K * x_curr
+ LinAlg::axpy(res, -1.0, b); // res = M * x_dot + K * x_curr - b
- BLAS::aypx(res, theta, _b_bar); // res = res * theta + _b_bar
- BLAS::matMultAdd(_M_bar, xdot, res, res); // rs += _M_bar * x_dot
+ LinAlg::aypx(res, theta, _b_bar); // res = res * theta + _b_bar
+ LinAlg::matMultAdd(_M_bar, xdot, res, res); // rs += _M_bar * x_dot
}
void MatrixTranslatorCrankNicolson<
ODESystemTag::FirstOrderImplicitQuasilinear>::
computeJacobian(GlobalMatrix const& Jac_in, GlobalMatrix& Jac_out) const
{
- namespace BLAS = MathLib::BLAS;
+ namespace LinAlg = MathLib::LinAlg;
auto const dxdot_dx = _crank_nicolson.getNewXWeight();
auto const theta = _crank_nicolson.getTheta();
// J = theta * Jac + dxdot_dx * _M_bar
- BLAS::copy(Jac_in, Jac_out);
- BLAS::scale(Jac_out, theta);
- BLAS::axpy(Jac_out, dxdot_dx, _M_bar);
+ LinAlg::copy(Jac_in, Jac_out);
+ LinAlg::scale(Jac_out, theta);
+ LinAlg::axpy(Jac_out, dxdot_dx, _M_bar);
}
void MatrixTranslatorCrankNicolson<
@@ -201,7 +201,7 @@ void MatrixTranslatorCrankNicolson<
pushMatrices(GlobalMatrix const& M, GlobalMatrix const& K,
GlobalVector const& b)
{
- namespace BLAS = MathLib::BLAS;
+ namespace LinAlg = MathLib::LinAlg;
auto const theta = _crank_nicolson.getTheta();
@@ -213,13 +213,13 @@ void MatrixTranslatorCrankNicolson<
auto const& x_old = _crank_nicolson.getXOld();
// _M_bar = (1.0-theta) * M;
- BLAS::copy(M, _M_bar);
- BLAS::scale(_M_bar, 1.0 - theta);
+ LinAlg::copy(M, _M_bar);
+ LinAlg::scale(_M_bar, 1.0 - theta);
// _b_bar = (1.0-theta) * (K * x_old - b)
- BLAS::matMult(K, x_old, _b_bar);
- BLAS::axpy(_b_bar, -1.0, b);
- BLAS::scale(_b_bar, 1.0 - theta);
+ LinAlg::matMult(K, x_old, _b_bar);
+ LinAlg::axpy(_b_bar, -1.0, b);
+ LinAlg::scale(_b_bar, 1.0 - theta);
}
} // NumLib
diff --git a/NumLib/ODESolver/MatrixTranslator.h b/NumLib/ODESolver/MatrixTranslator.h
index 47111033a1111a1c6fed3cc5802b229f57ebb8e9..576c1b29641204f0acd03957355c6d92e0b9f888 100644
--- a/NumLib/ODESolver/MatrixTranslator.h
+++ b/NumLib/ODESolver/MatrixTranslator.h
@@ -15,7 +15,7 @@
#include "TimeDiscretization.h"
#include "Types.h"
-#include "MathLib/LinAlg/BLAS.h"
+#include "MathLib/LinAlg/LinAlg.h"
namespace NumLib
{
diff --git a/NumLib/ODESolver/NonlinearSolver.cpp b/NumLib/ODESolver/NonlinearSolver.cpp
index f41741504af9deb7c15491dd1dca7b78b89d15d5..87949c8ea211610266fa42fdc07b7cc34970ead0 100644
--- a/NumLib/ODESolver/NonlinearSolver.cpp
+++ b/NumLib/ODESolver/NonlinearSolver.cpp
@@ -16,7 +16,7 @@
#include "BaseLib/ConfigTree.h"
#include "BaseLib/Error.h"
-#include "MathLib/LinAlg/BLAS.h"
+#include "MathLib/LinAlg/LinAlg.h"
#include "MathLib/LinAlg/VectorNorms.h"
#include "NumLib/DOF/GlobalMatrixProviders.h"
@@ -32,7 +32,7 @@ bool NonlinearSolver<NonlinearSolverTag::Picard>::solve(
GlobalVector& x,
std::function<void(unsigned, GlobalVector const&)> const& postIterationCallback)
{
- namespace BLAS = MathLib::BLAS;
+ namespace LinAlg= MathLib::LinAlg;
auto& sys = *_equation_system;
auto& A =
@@ -45,7 +45,7 @@ bool NonlinearSolver<NonlinearSolverTag::Picard>::solve(
bool error_norms_met = false;
- BLAS::copy(x, x_new); // set initial guess, TODO save the copy
+ LinAlg::copy(x, x_new); // set initial guess, TODO save the copy
unsigned iteration = 1;
for (; iteration <= _maxiter; ++iteration)
@@ -86,7 +86,7 @@ bool NonlinearSolver<NonlinearSolverTag::Picard>::solve(
// Copy new solution to x.
// Thereby the failed solution can be used by the caller for
// debugging purposes.
- BLAS::copy(x_new, x);
+ LinAlg::copy(x_new, x);
break;
case IterationResult::REPEAT_ITERATION:
INFO(
@@ -104,17 +104,17 @@ bool NonlinearSolver<NonlinearSolverTag::Picard>::solve(
break;
}
- auto const norm_x = BLAS::norm2(x);
+ auto const norm_x = LinAlg::norm2(x);
// x is used as delta_x in order to compute the error.
- BLAS::aypx(x, -1.0, x_new); // x = _x_new - x
- auto const error_dx = BLAS::norm2(x);
+ LinAlg::aypx(x, -1.0, x_new); // x = _x_new - x
+ auto const error_dx = LinAlg::norm2(x);
INFO(
"Picard: Iteration #%u |dx|=%.4e, |x|=%.4e, |dx|/|x|=%.4e,"
" tolerance(dx)=%.4e",
iteration, error_dx, norm_x, error_dx / norm_x, _tol);
// Update x s.t. in the next iteration we will compute the right delta x
- BLAS::copy(x_new, x);
+ LinAlg::copy(x_new, x);
if (error_dx < _tol)
{
@@ -156,7 +156,7 @@ bool NonlinearSolver<NonlinearSolverTag::Newton>::solve(
GlobalVector& x,
std::function<void(unsigned, GlobalVector const&)> const& postIterationCallback)
{
- namespace BLAS = MathLib::BLAS;
+ namespace LinAlg = MathLib::LinAlg;
auto& sys = *_equation_system;
auto& res = MathLib::GlobalVectorProvider<GlobalVector>::provider.getVector(
@@ -171,7 +171,7 @@ bool NonlinearSolver<NonlinearSolverTag::Newton>::solve(
// TODO be more efficient
// init _minus_delta_x to the right size and 0.0
- BLAS::copy(x, minus_delta_x);
+ LinAlg::copy(x, minus_delta_x);
minus_delta_x.setZero();
unsigned iteration = 1;
@@ -186,7 +186,7 @@ bool NonlinearSolver<NonlinearSolverTag::Newton>::solve(
sys.getJacobian(J);
sys.applyKnownSolutionsNewton(J, res, minus_delta_x);
- auto const error_res = BLAS::norm2(res);
+ auto const error_res = LinAlg::norm2(res);
// std::cout << " J:\n" << Eigen::MatrixXd(J) << std::endl;
@@ -204,7 +204,7 @@ bool NonlinearSolver<NonlinearSolverTag::Newton>::solve(
auto& x_new =
MathLib::GlobalVectorProvider<GlobalVector>::provider.getVector(
x, _x_new_id);
- BLAS::axpy(x_new, -_alpha, minus_delta_x);
+ LinAlg::axpy(x_new, -_alpha, minus_delta_x);
if (postIterationCallback)
postIterationCallback(iteration, x_new);
@@ -231,7 +231,7 @@ bool NonlinearSolver<NonlinearSolverTag::Newton>::solve(
// TODO could be done via swap. Note: that also requires swapping
// the ids. Same for the Picard scheme.
- BLAS::copy(x_new, x); // copy new solution to x
+ LinAlg::copy(x_new, x); // copy new solution to x
MathLib::GlobalVectorProvider<GlobalVector>::provider.releaseVector(
x_new);
}
@@ -243,8 +243,8 @@ bool NonlinearSolver<NonlinearSolverTag::Newton>::solve(
break;
}
- auto const error_dx = BLAS::norm2(minus_delta_x);
- auto const norm_x = BLAS::norm2(x);
+ auto const error_dx = LinAlg::norm2(minus_delta_x);
+ auto const norm_x = LinAlg::norm2(x);
INFO(
"Newton: Iteration #%u |dx|=%.4e, |r|=%.4e, |x|=%.4e, "
"|dx|/|x|=%.4e,"
diff --git a/NumLib/ODESolver/TimeDiscretization.h b/NumLib/ODESolver/TimeDiscretization.h
index 2f3a36ac382153fd1701b50503605181c8753520..eb1f48253816ac87454389205950c4e7e764e759 100644
--- a/NumLib/ODESolver/TimeDiscretization.h
+++ b/NumLib/ODESolver/TimeDiscretization.h
@@ -12,7 +12,7 @@
#include <vector>
-#include "MathLib/LinAlg/BLAS.h"
+#include "MathLib/LinAlg/LinAlg.h"
#include "NumLib/DOF/GlobalMatrixProviders.h"
#include "Types.h"
@@ -148,13 +148,13 @@ public:
//! \f$.
void getXdot(GlobalVector const& x_at_new_timestep, GlobalVector& xdot) const
{
- namespace BLAS = MathLib::BLAS;
+ namespace LinAlg = MathLib::LinAlg;
auto const dxdot_dx = getNewXWeight();
// xdot = dxdot_dx * x_at_new_timestep - x_old
getWeightedOldX(xdot);
- BLAS::axpby(xdot, dxdot_dx, -1.0, x_at_new_timestep);
+ LinAlg::axpby(xdot, dxdot_dx, -1.0, x_at_new_timestep);
}
//! Returns \f$ \alpha = \partial \hat x / \partial x_N \f$.
@@ -220,13 +220,13 @@ public:
void setInitialState(const double t0, GlobalVector const& x0) override
{
_t = t0;
- MathLib::BLAS::copy(x0, _x_old);
+ MathLib::LinAlg::copy(x0, _x_old);
}
void pushState(const double /*t*/, GlobalVector const& x,
InternalMatrixStorage const&) override
{
- MathLib::BLAS::copy(x, _x_old);
+ MathLib::LinAlg::copy(x, _x_old);
}
void nextTimestep(const double t, const double delta_t) override
@@ -239,11 +239,11 @@ public:
double getNewXWeight() const override { return 1.0 / _delta_t; }
void getWeightedOldX(GlobalVector& y) const override
{
- namespace BLAS = MathLib::BLAS;
+ namespace LinAlg = MathLib::LinAlg;
// y = x_old / delta_t
- BLAS::copy(_x_old, y);
- BLAS::scale(y, 1.0 / _delta_t);
+ LinAlg::copy(_x_old, y);
+ LinAlg::scale(y, 1.0 / _delta_t);
}
private:
@@ -270,13 +270,13 @@ public:
{
_t = t0;
_t_old = t0;
- MathLib::BLAS::copy(x0, _x_old);
+ MathLib::LinAlg::copy(x0, _x_old);
}
void pushState(const double /*t*/, GlobalVector const& x,
InternalMatrixStorage const&) override
{
- MathLib::BLAS::copy(x, _x_old);
+ MathLib::LinAlg::copy(x, _x_old);
}
void nextTimestep(const double t, const double delta_t) override
@@ -300,11 +300,11 @@ public:
double getNewXWeight() const override { return 1.0 / _delta_t; }
void getWeightedOldX(GlobalVector& y) const override
{
- namespace BLAS = MathLib::BLAS;
+ namespace LinAlg = MathLib::LinAlg;
// y = x_old / delta_t
- BLAS::copy(_x_old, y);
- BLAS::scale(y, 1.0 / _delta_t);
+ LinAlg::copy(_x_old, y);
+ LinAlg::scale(y, 1.0 / _delta_t);
}
bool isLinearTimeDisc() const override { return true; }
@@ -344,13 +344,13 @@ public:
void setInitialState(const double t0, GlobalVector const& x0) override
{
_t = t0;
- MathLib::BLAS::copy(x0, _x_old);
+ MathLib::LinAlg::copy(x0, _x_old);
}
void pushState(const double, GlobalVector const& x,
InternalMatrixStorage const& strg) override
{
- MathLib::BLAS::copy(x, _x_old);
+ MathLib::LinAlg::copy(x, _x_old);
strg.pushMatrices();
}
@@ -364,11 +364,11 @@ public:
double getNewXWeight() const override { return 1.0 / _delta_t; }
void getWeightedOldX(GlobalVector& y) const override
{
- namespace BLAS = MathLib::BLAS;
+ namespace LinAlg = MathLib::LinAlg;
// y = x_old / delta_t
- BLAS::copy(_x_old, y);
- BLAS::scale(y, 1.0 / _delta_t);
+ LinAlg::copy(_x_old, y);
+ LinAlg::scale(y, 1.0 / _delta_t);
}
bool needsPreload() const override { return true; }
@@ -439,7 +439,7 @@ public:
void pushState(const double, GlobalVector const& x,
InternalMatrixStorage const&) override
{
- namespace BLAS = MathLib::BLAS;
+ namespace LinAlg = MathLib::LinAlg;
// TODO use boost cirular buffer?
// until _xs_old is filled, lower-order BDF formulas are used.
@@ -450,7 +450,7 @@ public:
}
else
{
- BLAS::copy(x, *_xs_old[_offset]);
+ LinAlg::copy(x, *_xs_old[_offset]);
_offset = (_offset + 1) %
_num_steps; // treat _xs_old as a circular buffer
}
@@ -471,22 +471,22 @@ public:
void getWeightedOldX(GlobalVector& y) const override
{
- namespace BLAS = MathLib::BLAS;
+ namespace LinAlg = MathLib::LinAlg;
auto const k = eff_num_steps();
auto const* const BDFk = detail::BDF_Coeffs[k - 1];
// compute linear combination \sum_{i=0}^{k-1} BDFk_{k-i} \cdot x_{n+i}
- BLAS::copy(*_xs_old[_offset], y); // _xs_old[offset] = x_n
- BLAS::scale(y, BDFk[k]);
+ LinAlg::copy(*_xs_old[_offset], y); // _xs_old[offset] = x_n
+ LinAlg::scale(y, BDFk[k]);
for (unsigned i = 1; i < k; ++i)
{
auto const off = (_offset + i) % k;
- BLAS::axpy(y, BDFk[k - i], *_xs_old[off]);
+ LinAlg::axpy(y, BDFk[k - i], *_xs_old[off]);
}
- BLAS::scale(y, 1.0 / _delta_t);
+ LinAlg::scale(y, 1.0 / _delta_t);
}
private:
diff --git a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
index 37bbb469ed6236ebc667283c201810ea7a2cc633..e83405d9f85110e3961c061d62b92ca27a8a0159 100644
--- a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
+++ b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
@@ -68,7 +68,7 @@ void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
NonlinearSolverTag::Newton>::
assembleResidualNewton(const GlobalVector& x_new_timestep)
{
- namespace BLAS = MathLib::BLAS;
+ namespace LinAlg = MathLib::LinAlg;
auto const t = _time_disc.getCurrentTime();
auto const& x_curr = _time_disc.getCurrentX(x_new_timestep);
@@ -79,16 +79,16 @@ void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
_ode.assemble(t, x_curr, *_M, *_K, *_b);
- BLAS::finalizeAssembly(*_M);
- BLAS::finalizeAssembly(*_K);
- BLAS::finalizeAssembly(*_b);
+ LinAlg::finalizeAssembly(*_M);
+ LinAlg::finalizeAssembly(*_K);
+ LinAlg::finalizeAssembly(*_b);
}
void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
NonlinearSolverTag::Newton>::
assembleJacobian(const GlobalVector& x_new_timestep)
{
- namespace BLAS = MathLib::BLAS;
+ namespace LinAlg = MathLib::LinAlg;
auto const t = _time_disc.getCurrentTime();
auto const& x_curr = _time_disc.getCurrentX(x_new_timestep);
@@ -102,7 +102,7 @@ void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
_ode.assembleJacobian(t, x_curr, xdot, dxdot_dx, *_M, dx_dx, *_K, *_Jac);
- MathLib::BLAS::finalizeAssembly(*_Jac);
+ MathLib::LinAlg::finalizeAssembly(*_Jac);
MathLib::GlobalVectorProvider<GlobalVector>::provider.releaseVector(xdot);
}
@@ -190,7 +190,7 @@ void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
NonlinearSolverTag::Picard>::
assembleMatricesPicard(const GlobalVector& x_new_timestep)
{
- namespace BLAS = MathLib::BLAS;
+ namespace LinAlg = MathLib::LinAlg;
auto const t = _time_disc.getCurrentTime();
auto const& x_curr = _time_disc.getCurrentX(x_new_timestep);
@@ -201,9 +201,9 @@ void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
_ode.assemble(t, x_curr, *_M, *_K, *_b);
- BLAS::finalizeAssembly(*_M);
- BLAS::finalizeAssembly(*_K);
- BLAS::finalizeAssembly(*_b);
+ LinAlg::finalizeAssembly(*_M);
+ LinAlg::finalizeAssembly(*_K);
+ LinAlg::finalizeAssembly(*_b);
}
void TimeDiscretizedODESystem<
diff --git a/Tests/MathLib/TestGlobalMatrixInterface.cpp b/Tests/MathLib/TestGlobalMatrixInterface.cpp
index aef7a1d382d863de335e770007ced7ca44ebce67..af0dc0486c8748c0646ecfedc8e0bf9732091326 100644
--- a/Tests/MathLib/TestGlobalMatrixInterface.cpp
+++ b/Tests/MathLib/TestGlobalMatrixInterface.cpp
@@ -15,7 +15,7 @@
#include <gtest/gtest.h>
-#include "MathLib/LinAlg/BLAS.h"
+#include "MathLib/LinAlg/LinAlg.h"
#if defined(USE_PETSC)
#include "MathLib/LinAlg/PETSc/PETScMatrix.h"
@@ -28,7 +28,7 @@
#include "NumLib/NumericsConfig.h"
-using namespace MathLib::BLAS;
+using namespace MathLib::LinAlg;
namespace
{
diff --git a/Tests/MathLib/TestGlobalVectorInterface.cpp b/Tests/MathLib/TestGlobalVectorInterface.cpp
index 1e242179a53673ce5f6253b10e1ba0e90832d276..6535fa19010ccc1a26b7f6a1eeab549dd6adf5cd 100644
--- a/Tests/MathLib/TestGlobalVectorInterface.cpp
+++ b/Tests/MathLib/TestGlobalVectorInterface.cpp
@@ -16,7 +16,7 @@
#include <gtest/gtest.h>
#include "../TestTools.h"
-#include "MathLib/LinAlg/BLAS.h"
+#include "MathLib/LinAlg/LinAlg.h"
#if defined(USE_PETSC)
#include "MathLib/LinAlg/PETSc/PETScVector.h"
@@ -27,7 +27,7 @@
#include "MathLib/LinAlg/FinalizeVectorAssembly.h"
#include "NumLib/NumericsConfig.h"
-using namespace MathLib::BLAS;
+using namespace MathLib::LinAlg;
namespace
{
diff --git a/Tests/MathLib/TestLinearSolver.cpp b/Tests/MathLib/TestLinearSolver.cpp
index 5707cf37638af5b76c5dd10950f0f7834b7e35bc..f462937cc5810509e8990c59de899a07cab7a56e 100644
--- a/Tests/MathLib/TestLinearSolver.cpp
+++ b/Tests/MathLib/TestLinearSolver.cpp
@@ -15,7 +15,7 @@
#include <gtest/gtest.h>
-#include "MathLib/LinAlg/BLAS.h"
+#include "MathLib/LinAlg/LinAlg.h"
#include "MathLib/LinAlg/Dense/DenseMatrix.h"
#include "MathLib/LinAlg/FinalizeMatrixAssembly.h"
#include "MathLib/LinAlg/ApplyKnownSolution.h"
@@ -170,7 +170,7 @@ void checkLinearSolverInterface(T_MATRIX& A, T_VECTOR& b,
double x1[6];
x.getGlobalVector(x0);
- MathLib::BLAS::matMult(A, x, b);
+ MathLib::LinAlg::matMult(A, x, b);
// apply BC
std::vector<int> bc_id; // Type must be int to match Petsc_Int
diff --git a/Tests/NumLib/ODEs.h b/Tests/NumLib/ODEs.h
index 05f8e10bff6e3dcda65312595cea319bf04fea6b..0b01130f63587072d8e8846f9479521a0141b9ab 100644
--- a/Tests/NumLib/ODEs.h
+++ b/Tests/NumLib/ODEs.h
@@ -10,7 +10,7 @@
#ifndef TESTS_NUMLIB_ODES_H
#define TESTS_NUMLIB_ODES_H
-#include "MathLib/LinAlg/BLAS.h"
+#include "MathLib/LinAlg/LinAlg.h"
#include "MathLib/LinAlg/UnifiedMatrixSetters.h"
#include "NumLib/ODESolver/ODESystem.h"
@@ -40,13 +40,13 @@ public:
const GlobalMatrix& M, const double dx_dx,
const GlobalMatrix& K, GlobalMatrix& Jac) override
{
- namespace BLAS = MathLib::BLAS;
+ namespace LinAlg = MathLib::LinAlg;
// compute Jac = M*dxdot_dx + dx_dx*K
- BLAS::copy(M, Jac);
- BLAS::scale(Jac, dxdot_dx);
+ LinAlg::copy(M, Jac);
+ LinAlg::scale(Jac, dxdot_dx);
if (dx_dx != 0.0)
- BLAS::axpy(Jac, dx_dx, K);
+ LinAlg::axpy(Jac, dx_dx, K);
}
MathLib::MatrixSpecifications getMatrixSpecifications() const override
@@ -69,14 +69,14 @@ public:
static void setIC(GlobalVector& x0)
{
MathLib::setVector(x0, { 1.0, 0.0 });
- MathLib::BLAS::finalizeAssembly(x0);
+ MathLib::LinAlg::finalizeAssembly(x0);
}
static GlobalVector solution(const double t)
{
GlobalVector v(2);
MathLib::setVector(v, {cos(t), sin(t)});
- MathLib::BLAS::finalizeAssembly(v);
+ MathLib::LinAlg::finalizeAssembly(v);
return v;
}
@@ -108,18 +108,18 @@ public:
GlobalMatrix const& M, const double dx_dx,
GlobalMatrix const& K, GlobalMatrix& Jac) override
{
- namespace BLAS = MathLib::BLAS;
+ namespace LinAlg = MathLib::LinAlg;
// compute Jac = M*dxdot_dx + dK_dx + dx_dx*K
- BLAS::copy(M, Jac);
- BLAS::scale(Jac, dxdot_dx);
+ LinAlg::copy(M, Jac);
+ LinAlg::scale(Jac, dxdot_dx);
MathLib::addToMatrix(Jac, { x[0] }); // add dK_dx
if (dx_dx != 0.0)
{
- BLAS::finalizeAssembly(Jac);
- BLAS::axpy(Jac, dx_dx, K);
+ LinAlg::finalizeAssembly(Jac);
+ LinAlg::axpy(Jac, dx_dx, K);
}
}
@@ -143,14 +143,14 @@ public:
static void setIC(GlobalVector& x0)
{
MathLib::setVector(x0, { 1.0 });
- MathLib::BLAS::finalizeAssembly(x0);
+ MathLib::LinAlg::finalizeAssembly(x0);
}
static GlobalVector solution(const double t)
{
GlobalVector v(1);
MathLib::setVector(v, { 1.0/t });
- MathLib::BLAS::finalizeAssembly(v);
+ MathLib::LinAlg::finalizeAssembly(v);
return v;
}
@@ -207,12 +207,12 @@ public:
auto const dx = xdot[0];
auto const dz = xdot[2];
- namespace BLAS = MathLib::BLAS;
+ namespace LinAlg = MathLib::LinAlg;
// Compute Jac = M dxdot/dx + dM/dx xdot + K dx/dx + dK/dx x - db/dx
- BLAS::copy(M, Jac);
- BLAS::scale(Jac, dxdot_dx); // Jac = M * dxdot_dx
+ LinAlg::copy(M, Jac);
+ LinAlg::scale(Jac, dxdot_dx); // Jac = M * dxdot_dx
/* dx_dx == 0 holds if and only if the ForwardEuler scheme is used.
*
@@ -233,7 +233,7 @@ public:
0.0, t*dz, 0.0,
omega*t*dx+omega*dz, 0.0, 0.0 });
- BLAS::axpy(Jac, dx_dx, K); // add K \cdot dx_dx
+ LinAlg::axpy(Jac, dx_dx, K); // add K \cdot dx_dx
// add dK/dx \cdot \dot x
MathLib::
@@ -284,7 +284,7 @@ public:
MathLib::setVector(x0, { sin(omega*t0)/omega/t0,
1.0/t0,
cos(omega*t0) });
- MathLib::BLAS::finalizeAssembly(x0);
+ MathLib::LinAlg::finalizeAssembly(x0);
// std::cout << "IC:\n" << Eigen::VectorXd(x0.getRawVector()) << "\n";
}
@@ -297,7 +297,7 @@ public:
MathLib::setVector(v, { sin(omega*t)/omega/t,
1.0/t,
cos(omega*t) });
- MathLib::BLAS::finalizeAssembly(v);
+ MathLib::LinAlg::finalizeAssembly(v);
return v;
}
diff --git a/Tests/NumLib/TestExtrapolation.cpp b/Tests/NumLib/TestExtrapolation.cpp
index c475a34e0c9be942ba8524a054f5c17d3b05ff6a..e58769e378f78953f65b54dcb3abc2044b326d8a 100644
--- a/Tests/NumLib/TestExtrapolation.cpp
+++ b/Tests/NumLib/TestExtrapolation.cpp
@@ -15,7 +15,7 @@
#include "MathLib/LinAlg/UnifiedMatrixSetters.h"
#include "NumLib/Assembler/VectorMatrixAssembler.h"
-#include "MathLib/LinAlg/BLAS.h"
+#include "MathLib/LinAlg/LinAlg.h"
#include "MeshLib/MeshGenerators/MeshGenerator.h"
@@ -264,7 +264,7 @@ void extrapolate(TestProcess const& pcs,
expected_extrapolated_global_nodal_values,
std::size_t const nnodes, std::size_t const nelements)
{
- namespace BLAS = MathLib::BLAS;
+ namespace LinAlg = MathLib::LinAlg;
auto const tolerance_dx = 20.0 * std::numeric_limits<double>::epsilon();
auto const tolerance_res = 5.0 * std::numeric_limits<double>::epsilon();
@@ -276,15 +276,15 @@ void extrapolate(TestProcess const& pcs,
ASSERT_EQ(nnodes, x_extra.size());
ASSERT_EQ(nelements, residual.size());
- auto const res_norm = BLAS::normMax(residual);
+ auto const res_norm = LinAlg::normMax(residual);
DBUG("maximum norm of residual: %g", res_norm);
EXPECT_GT(tolerance_res, res_norm);
auto delta_x = MathLib::MatrixVectorTraits<GlobalVector>::newInstance(
expected_extrapolated_global_nodal_values);
- BLAS::axpy(*delta_x, -1.0, x_extra); // delta_x = x_expected - x_extra
+ LinAlg::axpy(*delta_x, -1.0, x_extra); // delta_x = x_expected - x_extra
- auto const dx_norm = BLAS::normMax(*delta_x);
+ auto const dx_norm = LinAlg::normMax(*delta_x);
DBUG("maximum norm of delta x: %g", dx_norm);
EXPECT_GT(tolerance_dx, dx_norm);
}
@@ -308,7 +308,7 @@ TEST(NumLib, DISABLED_Extrapolation)
for (unsigned integration_order : {2, 3, 4})
{
- namespace BLAS = MathLib::BLAS;
+ namespace LinAlg = MathLib::LinAlg;
const double mesh_length = 1.0;
const double mesh_elements_in_each_direction = 5.0;
@@ -339,7 +339,7 @@ TEST(NumLib, DISABLED_Extrapolation)
// expect 2*x as extraplation result for derived quantity
auto two_x = MathLib::MatrixVectorTraits<GlobalVector>::newInstance(*x);
- BLAS::axpy(*two_x, 1.0, *x); // two_x = x + x
+ LinAlg::axpy(*two_x, 1.0, *x); // two_x = x + x
// test extrapolation of a quantity that is derived from some integration
// point values