// ============================================================================ // == == // == AMDiS - Adaptive multidimensional simulations == // == == // ============================================================================ // == == // == TU Dresden == // == == // == Institut für Wissenschaftliches Rechnen == // == Zellescher Weg 12-14 == // == 01069 Dresden == // == germany == // == == // ============================================================================ // == == // == https://gforge.zih.tu-dresden.de/projects/amdis/ == // == == // ============================================================================ #ifndef MTL4SOLVER_H #define MTL4SOLVER_H #include "OEMSolver.h" #include "MatrixStreams.h" #include namespace AMDiS { template< typename MTLMatrix, typename MTLVector, typename Worker > class MTL4Solver : public OEMSolver { MTLMatrix matrix; protected: Worker worker; template< typename Matrix, typename Vector, typename Mapper > int solve(const Matrix& A, Vector& x, Vector& b, Mapper& mapper) { if(num_rows(matrix) == 0 || !getMultipleRhs() ) { matrix.change_dim(mapper.nRow(), mapper.nCol()); set_to_zero(matrix); MatMap< const Matrix, Mapper > matMap(A,mapper); matrix << matMap; worker.init(matrix); } MTLVector mtl_x(mapper.nRow()); set_to_zero(mtl_x); VecMap< Vector, Mapper > xVecMap(x, mapper); mtl_x << xVecMap; MTLVector mtl_b(mapper.nRow()); set_to_zero(mtl_b); VecMap< Vector, Mapper> bVecMap(b, mapper); mtl_b << bVecMap; error = worker.solve(matrix ,mtl_x, mtl_b); MTLVector r(mtl_b); r -= matrix * mtl_x; double residual = two_norm(r); MSG("MTL4Solver: ||b-Ax||= %e\n", residual); mtl_x >> xVecMap; return error; } public: MTL4Solver(std::string n): OEMSolver(n), matrix(0,0), worker(this) {} virtual int solveSystem(const SolverMatrix >& A, SystemVector& x, SystemVector& b, VectorialMapper& m) { return solve(A,x,b,m); } }; } #endif