1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra. 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 42b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// Copyright (C) 2011-2014 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, 302b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const Preconditioner& precond, Index& 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; 392b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Index maxIters = iters; 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 412b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Index n = mat.cols(); 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorType r = rhs - mat * x; 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorType r0 = r; 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar r0_sqnorm = r0.squaredNorm(); 467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar rhs_sqnorm = rhs.squaredNorm(); 477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(rhs_sqnorm == 0) 487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x.setZero(); 507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return true; 517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar rho = 1; 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar alpha = 1; 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar w = 1; 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorType v = VectorType::Zero(n), p = VectorType::Zero(n); 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorType y(n), z(n); 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorType kt(n), ks(n); 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorType s(n), t(n); 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 622b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang RealScalar tol2 = tol*tol*rhs_sqnorm; 637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar eps2 = NumTraits<Scalar>::epsilon()*NumTraits<Scalar>::epsilon(); 642b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Index i = 0; 652b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Index restarts = 0; 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 672b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang while ( r.squaredNorm() > tol2 && i<maxIters ) 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar rho_old = rho; 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath rho = r0.dot(r); 727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (abs(rho) < eps2*r0_sqnorm) 737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 742b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // The new residual vector became too orthogonal to the arbitrarily chosen direction r0 757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Let's restart with a new r0: 762b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang r = rhs - mat * x; 777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez r0 = r; 787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez rho = r0_sqnorm = r.squaredNorm(); 797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(restarts++ == 0) 807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez i = 0; 817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar beta = (rho/rho_old) * (alpha / w); 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath p = r + beta * (p - w * v); 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath y = precond.solve(p); 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath v.noalias() = mat * y; 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath alpha = rho / r0.dot(v); 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath s = r - alpha * v; 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath z = precond.solve(s); 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath t.noalias() = mat * z; 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar tmp = t.squaredNorm(); 967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(tmp>RealScalar(0)) 977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez w = t.dot(s) / tmp; 987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else 997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez w = Scalar(0); 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x += alpha * y + w * z; 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath r = s - w * t; 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++i; 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez tol_error = sqrt(r.squaredNorm()/rhs_sqnorm); 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath iters = i; 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return true; 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate< typename _MatrixType, 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> > 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass BiCGSTAB; 114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate< typename _MatrixType, typename _Preconditioner> 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct traits<BiCGSTAB<_MatrixType,_Preconditioner> > 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _MatrixType MatrixType; 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _Preconditioner Preconditioner; 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \ingroup IterativeLinearSolvers_Module 127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief A bi conjugate gradient stabilized solver for sparse square problems 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This class allows to solve for A.x = b sparse linear problems using a bi conjugate gradient 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * stabilized algorithm. The vectors x and b can be either dense or sparse. 131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam _MatrixType the type of the sparse matrix A, can be a dense or a sparse matrix. 133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam _Preconditioner the type of the preconditioner. Default is DiagonalPreconditioner 134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 1352b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \implsparsesolverconcept 1362b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The maximal number of iterations and tolerance value can be controlled via the setMaxIterations() 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * and setTolerance() methods. The defaults are the size of the problem for the maximal number of iterations 139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * and NumTraits<Scalar>::epsilon() for the tolerance. 140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 1412b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * The tolerance corresponds to the relative residual error: |Ax-b|/|b| 1422b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 1432b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \b Performance: when using sparse matrices, best performance is achied for a row-major sparse matrix format. 1442b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * Moreover, in this case multi-threading can be exploited if the user code is compiled with OpenMP enabled. 1452b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * See \ref TopicMultiThreading for details. 1462b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 147c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This class can be used as the direct solver classes. Here is a typical usage example: 1482b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \include BiCGSTAB_simple.cpp 149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 150c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * By default the iterations start with x=0 as an initial guess of the solution. 151a829215e078ace896f52702caa0c27608f40e3b0Miao Wang * One can control the start using the solveWithGuess() method. 152c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 1532b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * BiCGSTAB can also be used in a matrix-free context, see the following \link MatrixfreeSolverExample example \endlink. 1542b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 155c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa class SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner 156c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 157c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate< typename _MatrixType, typename _Preconditioner> 158c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass BiCGSTAB : public IterativeSolverBase<BiCGSTAB<_MatrixType,_Preconditioner> > 159c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 160c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef IterativeSolverBase<BiCGSTAB> Base; 1612b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang using Base::matrix; 162c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath using Base::m_error; 163c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath using Base::m_iterations; 164c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath using Base::m_info; 165c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath using Base::m_isInitialized; 166c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathpublic: 167c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _MatrixType MatrixType; 168c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 169c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::RealScalar RealScalar; 170c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _Preconditioner Preconditioner; 171c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 172c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathpublic: 173c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 174c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Default constructor. */ 175c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath BiCGSTAB() : Base() {} 176c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 177c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Initialize the solver with matrix \a A for further \c Ax=b solving. 178c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 179c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This constructor is a shortcut for the default constructor followed 180c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * by a call to compute(). 181c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 182c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \warning this class stores a reference to the matrix A as well as some 183c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * precomputed values that depend on it. Therefore, if \a A is changed 184c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * this class becomes invalid. Call compute() to update it with the new 185c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * matrix A, or modify a copy of A. 186c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1872b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename MatrixDerived> 1882b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang explicit BiCGSTAB(const EigenBase<MatrixDerived>& A) : Base(A.derived()) {} 189c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 190c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ~BiCGSTAB() {} 1912b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 192c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \internal */ 193c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename Rhs,typename Dest> 1942b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang void _solve_with_guess_impl(const Rhs& b, Dest& x) const 195c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 196c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool failed = false; 1972b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang for(Index j=0; j<b.cols(); ++j) 198c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 199c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_iterations = Base::maxIterations(); 200c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_error = Base::m_tolerance; 201c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 202c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename Dest::ColXpr xj(x,j); 2032b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if(!internal::bicgstab(matrix(), b.col(j), xj, Base::m_preconditioner, m_iterations, m_error)) 204c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath failed = true; 205c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 206c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_info = failed ? NumericalIssue 207c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : m_error <= Base::m_tolerance ? Success 208c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : NoConvergence; 209c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized = true; 210c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 211c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 212c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \internal */ 2132b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang using Base::_solve_impl; 214c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename Rhs,typename Dest> 2152b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang void _solve_impl(const MatrixBase<Rhs>& b, Dest& x) const 216c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 2172b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang x.resize(this->rows(),b.cols()); 2182b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang x.setZero(); 2192b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _solve_with_guess_impl(b,x); 220c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 221c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 222c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathprotected: 223c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 224c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 225c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 226c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 227c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 228c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_BICGSTAB_H 229