Commit 68c8f240 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

Merge branch 'feature/linear-solver-residuum' into 'master'

Add residuum output for linear solvers

See merge request !267
parents b0f55468 e1cc8d4a
...@@ -9,6 +9,7 @@ ...@@ -9,6 +9,7 @@
#include <amdis/linearalgebra/mtl/Traits.hpp> #include <amdis/linearalgebra/mtl/Traits.hpp>
#include <amdis/linearalgebra/mtl/MatrixBackend.hpp> #include <amdis/linearalgebra/mtl/MatrixBackend.hpp>
#include <amdis/linearalgebra/mtl/VectorBackend.hpp> #include <amdis/linearalgebra/mtl/VectorBackend.hpp>
#include <amdis/linearalgebra/mtl/Operations.hpp>
#elif AMDIS_BACKEND == AMDIS_BACKEND_EIGEN #elif AMDIS_BACKEND == AMDIS_BACKEND_EIGEN
...@@ -18,6 +19,7 @@ ...@@ -18,6 +19,7 @@
#include <amdis/linearalgebra/eigen/Traits.hpp> #include <amdis/linearalgebra/eigen/Traits.hpp>
#include <amdis/linearalgebra/eigen/MatrixBackend.hpp> #include <amdis/linearalgebra/eigen/MatrixBackend.hpp>
#include <amdis/linearalgebra/eigen/VectorBackend.hpp> #include <amdis/linearalgebra/eigen/VectorBackend.hpp>
#include <amdis/linearalgebra/eigen/Operations.hpp>
#elif AMDIS_BACKEND == AMDIS_BACKEND_PETSC #elif AMDIS_BACKEND == AMDIS_BACKEND_PETSC
...@@ -26,6 +28,7 @@ ...@@ -26,6 +28,7 @@
#include <amdis/linearalgebra/petsc/Traits.hpp> #include <amdis/linearalgebra/petsc/Traits.hpp>
#include <amdis/linearalgebra/petsc/MatrixBackend.hpp> #include <amdis/linearalgebra/petsc/MatrixBackend.hpp>
#include <amdis/linearalgebra/petsc/VectorBackend.hpp> #include <amdis/linearalgebra/petsc/VectorBackend.hpp>
#include <amdis/linearalgebra/petsc/Operations.hpp>
#elif AMDIS_BACKEND == AMDIS_BACKEND_ISTL #elif AMDIS_BACKEND == AMDIS_BACKEND_ISTL
...@@ -36,6 +39,7 @@ ...@@ -36,6 +39,7 @@
#include <amdis/linearalgebra/istl/Traits.hpp> #include <amdis/linearalgebra/istl/Traits.hpp>
#include <amdis/linearalgebra/istl/MatrixBackend.hpp> #include <amdis/linearalgebra/istl/MatrixBackend.hpp>
#include <amdis/linearalgebra/istl/VectorBackend.hpp> #include <amdis/linearalgebra/istl/VectorBackend.hpp>
#include <amdis/linearalgebra/istl/Operations.hpp>
#endif #endif
......
...@@ -359,8 +359,14 @@ solve(AdaptInfo& /*adaptInfo*/, bool createMatrixData, bool storeMatrixData) ...@@ -359,8 +359,14 @@ solve(AdaptInfo& /*adaptInfo*/, bool createMatrixData, bool storeMatrixData)
linearSolver_->finish(); linearSolver_->finish();
info(2, "solution of discrete system needed {} seconds", stat.elapsed); info(2, "solution of discrete system needed {} seconds", stat.elapsed);
info(1, "Residual norm: ||b-Ax|| = {}",
residuum(systemMatrix_->impl(),solution_->impl(), rhs_->impl()));
if (stat.reduction >= 0.0) if (stat.reduction >= 0.0)
info(1, "Residual norm: ||b-Ax||/||b-Ax^0|| = {}", stat.reduction); info(2, "Relative residual norm: ||b-Ax||/||b|| = {}", stat.reduction);
else
info(2, "Relative residual norm: ||b-Ax||/||b|| = {}",
relResiduum(systemMatrix_->impl(),solution_->impl(), rhs_->impl()));
bool ignoreConverged = false; bool ignoreConverged = false;
Parameters::get(name_ + "->solver->ignore converged", ignoreConverged); Parameters::get(name_ + "->solver->ignore converged", ignoreConverged);
......
...@@ -5,6 +5,7 @@ install(FILES ...@@ -5,6 +5,7 @@ install(FILES
LinearSolver.hpp LinearSolver.hpp
MatrixBackend.hpp MatrixBackend.hpp
MatrixSize.hpp MatrixSize.hpp
Operations.hpp
PreconConfig.hpp PreconConfig.hpp
SolverConfig.hpp SolverConfig.hpp
Solvers.hpp Solvers.hpp
......
#pragma once
#include <amdis/linearalgebra/eigen/MatrixBackend.hpp>
#include <amdis/linearalgebra/eigen/VectorBackend.hpp>
namespace AMDiS {
// ||b - A*x||
template <class T1, int O1, class T2, class T3>
auto residuum(EigenSparseMatrix<T1,O1> const& A, EigenVector<T2> const& x, EigenVector<T3> const& b)
{
auto r = b.vector();
r.noalias() -= A.matrix() * x.vector();
return r.norm();
}
// ||b - A*x|| / ||b||
template <class T1, int O1, class T2, class T3>
auto relResiduum(EigenSparseMatrix<T1,O1> const& A, EigenVector<T2> const& x, EigenVector<T3> const& b)
{
return residuum(A,x,b) / b.vector().norm();
}
} // end namespace AMDiS
...@@ -12,6 +12,7 @@ install(FILES ...@@ -12,6 +12,7 @@ install(FILES
ISTLSolverCreator.hpp ISTLSolverCreator.hpp
LinearSolver.hpp LinearSolver.hpp
MatrixBackend.hpp MatrixBackend.hpp
Operations.hpp
Preconditioners.hpp Preconditioners.hpp
PreconWrapper.hpp PreconWrapper.hpp
Solvers.hpp Solvers.hpp
......
#pragma once
#include <amdis/linearalgebra/istl/MatrixBackend.hpp>
#include <amdis/linearalgebra/istl/VectorBackend.hpp>
namespace AMDiS {
// ||b - A*x||
template <class T1, class C1, class T2, class T3>
auto residuum(ISTLBCRSMatrix<T1,C1> const& A, ISTLBlockVector<T2> const& x, ISTLBlockVector<T3> const& b)
{
auto r = b.vector();
A.matrix().mmv(x.vector(), r);
return r.two_norm();
}
// ||b - A*x|| / ||b||
template <class T1, class C1, class T2, class T3>
auto relResiduum(ISTLBCRSMatrix<T1,C1> const& A, ISTLBlockVector<T2> const& x, ISTLBlockVector<T3> const& b)
{
return residuum(A,x,b) / b.vector().two_norm();
}
} // end namespace AMDiS
...@@ -4,6 +4,7 @@ install(FILES ...@@ -4,6 +4,7 @@ install(FILES
KrylovRunner.hpp KrylovRunner.hpp
LinearSolver.hpp LinearSolver.hpp
MatrixBackend.hpp MatrixBackend.hpp
Operations.hpp
Preconditioner.hpp Preconditioner.hpp
PreconditionerInterface.hpp PreconditionerInterface.hpp
Preconditioners.hpp Preconditioners.hpp
......
#pragma once
#include <boost/numeric/mtl/operation/two_norm.hpp>
#include <amdis/linearalgebra/mtl/MatrixBackend.hpp>
#include <amdis/linearalgebra/mtl/VectorBackend.hpp>
namespace AMDiS {
// ||b - A*x||
template <class T1, class T2, class T3>
auto residuum(MTLSparseMatrix<T1> const& A, MTLVector<T2> const& x, MTLVector<T3> const& b)
{
auto r = b.vector();
r -= A.matrix() * x.vector();
return two_norm(r);
}
// ||b - A*x|| / ||b||
template <class T1, class T2, class T3>
auto relResiduum(MTLSparseMatrix<T1> const& A, MTLVector<T2> const& x, MTLVector<T3> const& b)
{
return residuum(A,x,b) / two_norm(b.vector());
}
} // end namespace AMDiS
...@@ -79,6 +79,7 @@ namespace AMDiS ...@@ -79,6 +79,7 @@ namespace AMDiS
umfpack_error_exit("solve", e.code); umfpack_error_exit("solve", e.code);
} }
stat.iterations = 1;
stat.converged = (code == UMFPACK_OK); stat.converged = (code == UMFPACK_OK);
stat.elapsed = t.elapsed(); stat.elapsed = t.elapsed();
} }
......
...@@ -6,6 +6,7 @@ install(FILES ...@@ -6,6 +6,7 @@ install(FILES
MatrixBackend.hpp MatrixBackend.hpp
MatrixNnzStructure.hpp MatrixNnzStructure.hpp
MatrixNnzStructure.inc.hpp MatrixNnzStructure.inc.hpp
Operations.hpp
Traits.hpp Traits.hpp
VectorBackend.hpp VectorBackend.hpp
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/amdis/linearalgebra/petsc) DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/amdis/linearalgebra/petsc)
#pragma once
#include <amdis/linearalgebra/petsc/MatrixBackend.hpp>
#include <amdis/linearalgebra/petsc/VectorBackend.hpp>
namespace AMDiS {
// ||b - A*x||
template <class DofMap>
auto residuum(PetscMatrix<DofMap> const& A, PetscVector<DofMap> const& x, PetscVector<DofMap> const& b)
{
Vec r;
VecDuplicate(b.vector(), &r);
MatMult(A.matrix(), x.vector(), r);
VecAXPY(r, -1.0, b.vector());
PetscReal res;
VecNorm(r, NORM_2, &res);
return res;
}
// ||b - A*x|| / ||b||
template <class DofMap>
auto relResiduum(PetscMatrix<DofMap> const& A, PetscVector<DofMap> const& x, PetscVector<DofMap> const& b)
{
PetscReal bNrm;
VecNorm(b.vector(), NORM_2, &bNrm);
return residuum(A,x,b) / bNrm;
}
} // end namespace AMDiS
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