From cacbc5679dac87a9b389b211d32b6e452770ee5b Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 9 Jun 2015 15:23:08 +0200 Subject: formatting --- unsupported/Eigen/src/IterativeSolvers/GMRES.h | 314 +++++++++++++------------ 1 file changed, 158 insertions(+), 156 deletions(-) (limited to 'unsupported/Eigen/src/IterativeSolvers/GMRES.h') diff --git a/unsupported/Eigen/src/IterativeSolvers/GMRES.h b/unsupported/Eigen/src/IterativeSolvers/GMRES.h index e70a864e1..3cbc2f525 100644 --- a/unsupported/Eigen/src/IterativeSolvers/GMRES.h +++ b/unsupported/Eigen/src/IterativeSolvers/GMRES.h @@ -16,189 +16,191 @@ namespace Eigen { namespace internal { /** - * Generalized Minimal Residual Algorithm based on the - * Arnoldi algorithm implemented with Householder reflections. - * - * Parameters: - * \param mat matrix of linear system of equations - * \param Rhs right hand side vector of linear system of equations - * \param x on input: initial guess, on output: solution - * \param precond preconditioner used - * \param iters on input: maximum number of iterations to perform - * on output: number of iterations performed - * \param restart number of iterations for a restart - * \param tol_error on input: relative residual tolerance - * on output: residuum achieved - * - * \sa IterativeMethods::bicgstab() - * - * - * For references, please see: - * - * Saad, Y. and Schultz, M. H. - * GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems. - * SIAM J.Sci.Stat.Comp. 7, 1986, pp. 856 - 869. - * - * Saad, Y. - * Iterative Methods for Sparse Linear Systems. - * Society for Industrial and Applied Mathematics, Philadelphia, 2003. - * - * Walker, H. F. - * Implementations of the GMRES method. - * Comput.Phys.Comm. 53, 1989, pp. 311 - 320. - * - * Walker, H. F. - * Implementation of the GMRES Method using Householder Transformations. - * SIAM J.Sci.Stat.Comp. 9, 1988, pp. 152 - 163. - * - */ +* Generalized Minimal Residual Algorithm based on the +* Arnoldi algorithm implemented with Householder reflections. +* +* Parameters: +* \param mat matrix of linear system of equations +* \param Rhs right hand side vector of linear system of equations +* \param x on input: initial guess, on output: solution +* \param precond preconditioner used +* \param iters on input: maximum number of iterations to perform +* on output: number of iterations performed +* \param restart number of iterations for a restart +* \param tol_error on input: relative residual tolerance +* on output: residuum achieved +* +* \sa IterativeMethods::bicgstab() +* +* +* For references, please see: +* +* Saad, Y. and Schultz, M. H. +* GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems. +* SIAM J.Sci.Stat.Comp. 7, 1986, pp. 856 - 869. +* +* Saad, Y. +* Iterative Methods for Sparse Linear Systems. +* Society for Industrial and Applied Mathematics, Philadelphia, 2003. +* +* Walker, H. F. +* Implementations of the GMRES method. +* Comput.Phys.Comm. 53, 1989, pp. 311 - 320. +* +* Walker, H. F. +* Implementation of the GMRES Method using Householder Transformations. +* SIAM J.Sci.Stat.Comp. 9, 1988, pp. 152 - 163. +* +*/ template bool gmres(const MatrixType & mat, const Rhs & rhs, Dest & x, const Preconditioner & precond, - Index &iters, const Index &restart, typename Dest::RealScalar & tol_error) { + Index &iters, const Index &restart, typename Dest::RealScalar & tol_error) { - using std::sqrt; - using std::abs; + using std::sqrt; + using std::abs; - typedef typename Dest::RealScalar RealScalar; - typedef typename Dest::Scalar Scalar; - typedef Matrix < Scalar, Dynamic, 1 > VectorType; - typedef Matrix < Scalar, Dynamic, Dynamic > FMatrixType; + typedef typename Dest::RealScalar RealScalar; + typedef typename Dest::Scalar Scalar; + typedef Matrix < Scalar, Dynamic, 1 > VectorType; + typedef Matrix < Scalar, Dynamic, Dynamic > FMatrixType; - RealScalar tol = tol_error; - const Index maxIters = iters; - iters = 0; + RealScalar tol = tol_error; + const Index maxIters = iters; + iters = 0; - const Index m = mat.rows(); + const Index m = mat.rows(); - // residual and preconditioned residual - VectorType p0 = rhs - mat*x; - VectorType r0 = precond.solve(p0); + // residual and preconditioned residual + VectorType p0 = rhs - mat*x; + VectorType r0 = precond.solve(p0); VectorType t(m), v(m), workspace(m); - const RealScalar r0Norm = r0.norm(); + const RealScalar r0Norm = r0.norm(); - // is initial guess already good enough? - if(r0Norm == 0) { - tol_error=0; - return true; - } + // is initial guess already good enough? + if(r0Norm == 0) + { + tol_error = 0; + return true; + } - // storage for Hessenberg matrix and Householder data - FMatrixType H = FMatrixType::Zero(m, restart + 1); - VectorType w = VectorType::Zero(restart + 1); - VectorType tau = VectorType::Zero(restart + 1); + // storage for Hessenberg matrix and Householder data + FMatrixType H = FMatrixType::Zero(m, restart + 1); + VectorType w = VectorType::Zero(restart + 1); + VectorType tau = VectorType::Zero(restart + 1); - // storage for Jacobi rotations - std::vector < JacobiRotation < Scalar > > G(restart); + // storage for Jacobi rotations + std::vector < JacobiRotation < Scalar > > G(restart); - // generate first Householder vector + // generate first Householder vector Ref H0_tail = H.col(0).tail(m - 1); - RealScalar beta; - r0.makeHouseholder(H0_tail, tau.coeffRef(0), beta); - w(0)=(Scalar) beta; - - for (Index k = 1; k <= restart; ++k) { - - ++iters; + RealScalar beta; + r0.makeHouseholder(H0_tail, tau.coeffRef(0), beta); + w(0) = Scalar(beta); - v = VectorType::Unit(m, k - 1); + for (Index k = 1; k <= restart; ++k) + { + ++iters; - // apply Householder reflections H_{1} ... H_{k-1} to v - for (Index i = k - 1; i >= 0; --i) { - v.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data()); - } + v = VectorType::Unit(m, k - 1); - // apply matrix M to v: v = mat * v; - t=mat*v; - v=precond.solve(t); + // apply Householder reflections H_{1} ... H_{k-1} to v + for (Index i = k - 1; i >= 0; --i) { + v.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data()); + } - // apply Householder reflections H_{k-1} ... H_{1} to v - for (Index i = 0; i < k; ++i) { - v.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data()); - } + // apply matrix M to v: v = mat * v; + t=mat*v; + v=precond.solve(t); - if (v.tail(m - k).norm() != 0.0) { - if (k <= restart) { + // apply Householder reflections H_{k-1} ... H_{1} to v + for (Index i = 0; i < k; ++i) { + v.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data()); + } - // generate new Householder vector + if (v.tail(m - k).norm() != 0.0) + { + if (k <= restart) + { + // generate new Householder vector Ref Hk_tail = H.col(k).tail(m - k - 1); - - v.tail(m - k).makeHouseholder(Hk_tail, tau.coeffRef(k), beta); - - // apply Householder reflection H_{k} to v - v.tail(m - k).applyHouseholderOnTheLeft(Hk_tail, tau.coeffRef(k), workspace.data()); - - } - } - - if (k > 1) { - for (Index i = 0; i < k - 1; ++i) { - // apply old Givens rotations to v - v.applyOnTheLeft(i, i + 1, G[i].adjoint()); - } - } - - if (k ().solveInPlace(y); - - // use Horner-like scheme to calculate solution vector - VectorType x_new = y(k - 1) * VectorType::Unit(m, k - 1); - - // apply Householder reflection H_{k} to x_new - x_new.tail(m - k + 1).applyHouseholderOnTheLeft(H.col(k - 1).tail(m - k), tau.coeffRef(k - 1), workspace.data()); - - for (Index i = k - 2; i >= 0; --i) { - x_new += y(i) * VectorType::Unit(m, i); - // apply Householder reflection H_{i} to x_new - x_new.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data()); - } - - x += x_new; - - if (stop) { - return true; - } else { - k=0; - - // reset data for restart - p0 = rhs - mat*x; - r0 = precond.solve(p0); + // apply Householder reflection H_{k} to v + v.tail(m - k).applyHouseholderOnTheLeft(Hk_tail, tau.coeffRef(k), workspace.data()); + } + } - // clear Hessenberg matrix and Householder data - H = FMatrixType::Zero(m, restart + 1); - w = VectorType::Zero(restart + 1); - tau = VectorType::Zero(restart + 1); + if (k > 1) + { + for (Index i = 0; i < k - 1; ++i) + { + // apply old Givens rotations to v + v.applyOnTheLeft(i, i + 1, G[i].adjoint()); + } + } - // generate first Householder vector - r0.makeHouseholder(H0_tail, tau.coeffRef(0), beta); - w(0)=(Scalar) beta; + if (k().solveInPlace(y); + + // use Horner-like scheme to calculate solution vector + VectorType x_new = y(k - 1) * VectorType::Unit(m, k - 1); + + // apply Householder reflection H_{k} to x_new + x_new.tail(m - k + 1).applyHouseholderOnTheLeft(H.col(k - 1).tail(m - k), tau.coeffRef(k - 1), workspace.data()); + + for (Index i = k - 2; i >= 0; --i) + { + x_new += y(i) * VectorType::Unit(m, i); + // apply Householder reflection H_{i} to x_new + x_new.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data()); + } + + x += x_new; + + if (stop) + { + return true; + } + else + { + k=0; + + // reset data for restart + p0 = rhs - mat*x; + r0 = precond.solve(p0); + + // clear Hessenberg matrix and Householder data + H = FMatrixType::Zero(m, restart + 1); + w = VectorType::Zero(restart + 1); + tau = VectorType::Zero(restart + 1); + + // generate first Householder vector + r0.makeHouseholder(H0_tail, tau.coeffRef(0), beta); + w(0)=(Scalar) beta; + } + } + } - return false; + return false; } @@ -314,8 +316,8 @@ public: failed = true; } m_info = failed ? NumericalIssue - : m_error <= Base::m_tolerance ? Success - : NoConvergence; + : m_error <= Base::m_tolerance ? Success + : NoConvergence; m_isInitialized = true; } -- cgit v1.2.3