1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra. 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr> 5c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr> 6c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 7c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This Source Code Form is subject to the terms of the Mozilla 8c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Public License v. 2.0. If a copy of the MPL was not distributed 9c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 10c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#ifndef EIGEN_BICGSTAB_H 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define EIGEN_BICGSTAB_H 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace Eigen { 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \internal Low-level bi conjugate gradient stabilized algorithm 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param mat The matrix A 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param rhs The right hand side vector b 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param x On input and initial solution, on output the computed solution. 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param precond A preconditioner being able to efficiently solve for an 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * approximation of Ax=b (regardless of b) 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param iters On input the max number of iteration, on output the number of performed iterations. 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param tol_error On input the tolerance error, on output an estimation of the relative error. 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \return false in the case of numerical issue, for example a break down of BiCGSTAB. 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, typename Rhs, typename Dest, typename Preconditioner> 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathbool bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x, 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Preconditioner& precond, int& iters, 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename Dest::RealScalar& tol_error) 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath using std::sqrt; 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath using std::abs; 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename Dest::RealScalar RealScalar; 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename Dest::Scalar Scalar; 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar,Dynamic,1> VectorType; 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar tol = tol_error; 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int maxIters = iters; 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int n = mat.cols(); 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorType r = rhs - mat * x; 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorType r0 = r; 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar r0_sqnorm = r0.squaredNorm(); 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar rho = 1; 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar alpha = 1; 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar w = 1; 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorType v = VectorType::Zero(n), p = VectorType::Zero(n); 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorType y(n), z(n); 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorType kt(n), ks(n); 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorType s(n), t(n); 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar tol2 = tol*tol; 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int i = 0; 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath while ( r.squaredNorm()/r0_sqnorm > tol2 && i<maxIters ) 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar rho_old = rho; 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath rho = r0.dot(r); 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (rho == Scalar(0)) return false; /* New search directions cannot be found */ 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar beta = (rho/rho_old) * (alpha / w); 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath p = r + beta * (p - w * v); 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath y = precond.solve(p); 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath v.noalias() = mat * y; 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath alpha = rho / r0.dot(v); 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath s = r - alpha * v; 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath z = precond.solve(s); 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath t.noalias() = mat * z; 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath w = t.dot(s) / t.squaredNorm(); 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x += alpha * y + w * z; 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath r = s - w * t; 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++i; 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tol_error = sqrt(r.squaredNorm()/r0_sqnorm); 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath iters = i; 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return true; 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate< typename _MatrixType, 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> > 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass BiCGSTAB; 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate< typename _MatrixType, typename _Preconditioner> 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct traits<BiCGSTAB<_MatrixType,_Preconditioner> > 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _MatrixType MatrixType; 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _Preconditioner Preconditioner; 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \ingroup IterativeLinearSolvers_Module 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief A bi conjugate gradient stabilized solver for sparse square problems 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This class allows to solve for A.x = b sparse linear problems using a bi conjugate gradient 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * stabilized algorithm. The vectors x and b can be either dense or sparse. 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam _MatrixType the type of the sparse matrix A, can be a dense or a sparse matrix. 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam _Preconditioner the type of the preconditioner. Default is DiagonalPreconditioner 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The maximal number of iterations and tolerance value can be controlled via the setMaxIterations() 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * and setTolerance() methods. The defaults are the size of the problem for the maximal number of iterations 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * and NumTraits<Scalar>::epsilon() for the tolerance. 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This class can be used as the direct solver classes. Here is a typical usage example: 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \code 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * int n = 10000; 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * VectorXd x(n), b(n); 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * SparseMatrix<double> A(n,n); 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * // fill A and b 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * BiCGSTAB<SparseMatrix<double> > solver; 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * solver(A); 126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * x = solver.solve(b); 127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * std::cout << "#iterations: " << solver.iterations() << std::endl; 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * std::cout << "estimated error: " << solver.error() << std::endl; 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * // update b, and solve again 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * x = solver.solve(b); 131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \endcode 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * By default the iterations start with x=0 as an initial guess of the solution. 134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * One can control the start using the solveWithGuess() method. Here is a step by 135c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * step execution example starting with a random guess and printing the evolution 136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * of the estimated error: 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * * \code 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * x = VectorXd::Random(n); 139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * solver.setMaxIterations(1); 140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * int i = 0; 141c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * do { 142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * x = solver.solveWithGuess(b,x); 143c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * std::cout << i << " : " << solver.error() << std::endl; 144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * ++i; 145c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * } while (solver.info()!=Success && i<100); 146c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \endcode 147c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Note that such a step by step excution is slightly slower. 148c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa class SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner 150c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 151c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate< typename _MatrixType, typename _Preconditioner> 152c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass BiCGSTAB : public IterativeSolverBase<BiCGSTAB<_MatrixType,_Preconditioner> > 153c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 154c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef IterativeSolverBase<BiCGSTAB> Base; 155c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath using Base::mp_matrix; 156c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath using Base::m_error; 157c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath using Base::m_iterations; 158c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath using Base::m_info; 159c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath using Base::m_isInitialized; 160c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathpublic: 161c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _MatrixType MatrixType; 162c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 163c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 164c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::RealScalar RealScalar; 165c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _Preconditioner Preconditioner; 166c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 167c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathpublic: 168c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 169c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Default constructor. */ 170c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath BiCGSTAB() : Base() {} 171c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 172c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Initialize the solver with matrix \a A for further \c Ax=b solving. 173c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 174c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This constructor is a shortcut for the default constructor followed 175c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * by a call to compute(). 176c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 177c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \warning this class stores a reference to the matrix A as well as some 178c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * precomputed values that depend on it. Therefore, if \a A is changed 179c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * this class becomes invalid. Call compute() to update it with the new 180c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * matrix A, or modify a copy of A. 181c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 182c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath BiCGSTAB(const MatrixType& A) : Base(A) {} 183c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 184c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ~BiCGSTAB() {} 185c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 186c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A 187c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \a x0 as an initial solution. 188c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 189c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa compute() 190c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 191c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename Rhs,typename Guess> 192c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline const internal::solve_retval_with_guess<BiCGSTAB, Rhs, Guess> 193c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath solveWithGuess(const MatrixBase<Rhs>& b, const Guess& x0) const 194c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 195c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "BiCGSTAB is not initialized."); 196c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(Base::rows()==b.rows() 197c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath && "BiCGSTAB::solve(): invalid number of rows of the right hand side matrix b"); 198c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return internal::solve_retval_with_guess 199c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath <BiCGSTAB, Rhs, Guess>(*this, b.derived(), x0); 200c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 201c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 202c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \internal */ 203c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename Rhs,typename Dest> 204c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void _solveWithGuess(const Rhs& b, Dest& x) const 205c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 206c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool failed = false; 207c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int j=0; j<b.cols(); ++j) 208c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 209c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iterations = Base::maxIterations(); 210c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_error = Base::m_tolerance; 211c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 212c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename Dest::ColXpr xj(x,j); 213c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(!internal::bicgstab(*mp_matrix, b.col(j), xj, Base::m_preconditioner, m_iterations, m_error)) 214c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath failed = true; 215c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 216c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_info = failed ? NumericalIssue 217c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : m_error <= Base::m_tolerance ? Success 218c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : NoConvergence; 219c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized = true; 220c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 221c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 222c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \internal */ 223c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename Rhs,typename Dest> 224c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void _solve(const Rhs& b, Dest& x) const 225c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 226c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x.setZero(); 227c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath _solveWithGuess(b,x); 228c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 229c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 230c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathprotected: 231c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 232c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 233c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 234c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 235c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 236c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 237c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename _MatrixType, typename _Preconditioner, typename Rhs> 238c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct solve_retval<BiCGSTAB<_MatrixType, _Preconditioner>, Rhs> 239c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : solve_retval_base<BiCGSTAB<_MatrixType, _Preconditioner>, Rhs> 240c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 241c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef BiCGSTAB<_MatrixType, _Preconditioner> Dec; 242c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) 243c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 244c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename Dest> void evalTo(Dest& dst) const 245c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 246c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dec()._solve(rhs(),dst); 247c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 248c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 249c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 250c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace internal 251c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 252c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 253c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 254c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_BICGSTAB_H 255