/****************************************************************************** * * AMDiS - Adaptive multidimensional simulations * * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved. * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis * * Authors: * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al. * * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. * * * This file is part of AMDiS * * See also license.opensource.txt in the distribution. * ******************************************************************************/ // Written by Thomas Witkowski #ifndef AMDIS_ITL_MINRES_INCLUDE #define AMDIS_ITL_MINRES_INCLUDE #include namespace itl { /// Minimal Residual method template < typename Matrix, typename Vector, typename LeftPreconditioner, typename RightPreconditioner, typename Iteration > int minres(const Matrix &A, Vector &x, const Vector &b, const LeftPreconditioner &L, const RightPreconditioner &R, Iteration& iter) { using std::abs; using math::reciprocal; typedef typename mtl::Collection::value_type Scalar; if (size(b) == 0) throw mtl::logic_error("empty rhs vector"); Scalar zero= math::zero(b[0]), one= math::one(b[0]); Vector v0(size(x), zero), v1(b - A * x), v2(v1), z1(solve(L, v1)), z2(size(x), zero); Vector w0(size(x), zero), w1(size(x), zero), w2(size(x), zero); Scalar s0(zero), s1(zero), c0(one), c1(one), gamma0(one); Scalar gamma1(sqrt(dot(z1, v1))), gamma2(zero), eta(gamma1); Scalar sigma1(one), alpha0(zero), alpha1(zero), alpha2(zero), alpha3(zero); while (!iter.finished(abs(eta))) { z1 *= reciprocal(gamma1); v2 = A * z1; sigma1 = dot(v2, z1); v2 += -(sigma1 / gamma1) * v1 - (gamma1 / gamma0) * v0; z2 = solve(L, v2); gamma2 = sqrt(dot(z2, v2)); alpha0 = c1 * sigma1 - c0 * s1 * gamma1; alpha1 = sqrt(alpha0 * alpha0 + gamma2 * gamma2); alpha2 = s1 * sigma1 + c0 * c1 * gamma1; alpha3 = s0 * gamma1; c0 = c1; c1 = alpha0 / alpha1; s0 = s1; s1 = gamma2 / alpha1; w2 = z1 - alpha3 * w0 - alpha2 * w1; w2 *= reciprocal(alpha1); x += c1 * eta * w2; eta *= -s1; w0 = w1; w1 = w2; v0 = v1; v1 = v2; z1 = z2; gamma0 = gamma1; gamma1 = gamma2; ++iter; } return iter; } } // namespace itl; #endif // AMDIS_ITL_MINRES_INCLUDE