1 #ifndef AMGCL_SOLVERS_BICGSTAB_HPP 2 #define AMGCL_SOLVERS_BICGSTAB_HPP 3 4 /* 5 The MIT License 6 7 Copyright (c) 2012-2021 Denis Demidov <dennis.demidov@gmail.com> 8 9 Permission is hereby granted, free of charge, to any person obtaining a copy 10 of this software and associated documentation files (the "Software"), to deal 11 in the Software without restriction, including without limitation the rights 12 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell 13 copies of the Software, and to permit persons to whom the Software is 14 furnished to do so, subject to the following conditions: 15 16 The above copyright notice and this permission notice shall be included in 17 all copies or substantial portions of the Software. 18 19 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR 20 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 21 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE 22 AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 23 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, 24 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN 25 THE SOFTWARE. 26 */ 27 28 /** 29 * \file amgcl/solver/bicgstab.hpp 30 * \author Denis Demidov <dennis.demidov@gmail.com> 31 * \brief BiCGStab iterative method. 32 */ 33 34 #include <tuple> 35 #include <iostream> 36 37 #include <amgcl/backend/interface.hpp> 38 #include <amgcl/solver/detail/default_inner_product.hpp> 39 #include <amgcl/solver/precond_side.hpp> 40 #include <amgcl/util.hpp> 41 42 namespace amgcl { 43 namespace solver { 44 45 /** BiConjugate Gradient Stabilized (BiCGSTAB) method. 46 * \rst 47 * The BiConjugate Gradient Stabilized method (Bi-CGSTAB) was developed to 48 * solve nonsymmetric linear systems while avoiding the often irregular 49 * convergence patterns of the Conjugate Gradient [Barr94]_. 50 * \endrst 51 */ 52 template < 53 class Backend, 54 class InnerProduct = detail::default_inner_product 55 > 56 class bicgstab { 57 public: 58 typedef Backend backend_type; 59 60 typedef typename Backend::vector vector; 61 typedef typename Backend::value_type value_type; 62 typedef typename Backend::params backend_params; 63 64 typedef typename math::scalar_of<value_type>::type scalar_type; 65 66 typedef typename math::inner_product_impl< 67 typename math::rhs_of<value_type>::type 68 >::return_type coef_type; 69 70 71 /// Solver parameters. 72 struct params { 73 /// Preconditioning kind (left/right). 74 preconditioner::side::type pside; 75 76 /// Maximum number of iterations. 77 size_t maxiter; 78 79 /// Target relative residual error. 80 scalar_type tol; 81 82 /// Target absolute residual error. 83 scalar_type abstol; 84 85 /// Always do at least one iteration. 86 bool check_after; 87 88 /// Ignore the trivial solution x=0 when rhs is zero. 89 //** Useful for searching for the null-space vectors of the system */ 90 bool ns_search; 91 92 /// Verbose output (show iterations and error) 93 bool verbose; 94 paramsamgcl::solver::bicgstab::params95 params() 96 : pside(preconditioner::side::right), maxiter(100), tol(1e-8), 97 abstol(std::numeric_limits<scalar_type>::min()), 98 check_after(false), ns_search(false), verbose(false) 99 {} 100 101 #ifndef AMGCL_NO_BOOST paramsamgcl::solver::bicgstab::params102 params(const boost::property_tree::ptree &p) 103 : AMGCL_PARAMS_IMPORT_VALUE(p, pside), 104 AMGCL_PARAMS_IMPORT_VALUE(p, maxiter), 105 AMGCL_PARAMS_IMPORT_VALUE(p, tol), 106 AMGCL_PARAMS_IMPORT_VALUE(p, abstol), 107 AMGCL_PARAMS_IMPORT_VALUE(p, check_after), 108 AMGCL_PARAMS_IMPORT_VALUE(p, ns_search), 109 AMGCL_PARAMS_IMPORT_VALUE(p, verbose) 110 { 111 check_params(p, {"pside", "maxiter", "tol", "abstol", 112 "check_after", "ns_search", "verbose"}); 113 } 114 getamgcl::solver::bicgstab::params115 void get(boost::property_tree::ptree &p, const std::string &path) const { 116 AMGCL_PARAMS_EXPORT_VALUE(p, path, pside); 117 AMGCL_PARAMS_EXPORT_VALUE(p, path, maxiter); 118 AMGCL_PARAMS_EXPORT_VALUE(p, path, tol); 119 AMGCL_PARAMS_EXPORT_VALUE(p, path, abstol); 120 AMGCL_PARAMS_EXPORT_VALUE(p, path, check_after); 121 AMGCL_PARAMS_EXPORT_VALUE(p, path, ns_search); 122 AMGCL_PARAMS_EXPORT_VALUE(p, path, verbose); 123 } 124 #endif 125 }; 126 127 /// Preallocates necessary data structures for the system of size \p n. bicgstab(size_t n,const params & prm=params (),const backend_params & backend_prm=backend_params (),const InnerProduct & inner_product=InnerProduct ())128 bicgstab( 129 size_t n, 130 const params &prm = params(), 131 const backend_params &backend_prm = backend_params(), 132 const InnerProduct &inner_product = InnerProduct() 133 ) 134 : prm(prm), n(n), 135 r ( Backend::create_vector(n, backend_prm) ), 136 p ( Backend::create_vector(n, backend_prm) ), 137 v ( Backend::create_vector(n, backend_prm) ), 138 s ( Backend::create_vector(n, backend_prm) ), 139 t ( Backend::create_vector(n, backend_prm) ), 140 rh( Backend::create_vector(n, backend_prm) ), 141 T ( Backend::create_vector(n, backend_prm) ), 142 inner_product(inner_product) 143 { } 144 145 /* Computes the solution for the given system matrix \p A and the 146 * right-hand side \p rhs. Returns the number of iterations made and 147 * the achieved residual as a ``std::tuple``. The solution vector 148 * \p x provides initial approximation in input and holds the computed 149 * solution on output. 150 * 151 * The system matrix may differ from the matrix used during 152 * initialization. This may be used for the solution of non-stationary 153 * problems with slowly changing coefficients. There is a strong chance 154 * that a preconditioner built for a time step will act as a reasonably 155 * good preconditioner for several subsequent time steps [DeSh12]_. 156 */ 157 template <class Matrix, class Precond, class Vec1, class Vec2> operator ()(const Matrix & A,const Precond & P,const Vec1 & rhs,Vec2 && x) const158 std::tuple<size_t, scalar_type> operator()( 159 const Matrix &A, const Precond &P, const Vec1 &rhs, Vec2 &&x) const 160 { 161 namespace side = preconditioner::side; 162 163 static const coef_type one = math::identity<coef_type>(); 164 static const coef_type zero = math::zero<coef_type>(); 165 166 ios_saver ss(std::cout); 167 168 scalar_type norm_rhs = norm(rhs); 169 if (norm_rhs < amgcl::detail::eps<scalar_type>(1)) { 170 if (prm.ns_search) { 171 norm_rhs = math::identity<scalar_type>(); 172 } else { 173 backend::clear(x); 174 return std::make_tuple(0, norm_rhs); 175 } 176 } 177 178 if (prm.pside == side::left) { 179 backend::residual(rhs, A, x, *rh); 180 P.apply(*rh, *r); 181 } else { 182 backend::residual(rhs, A, x, *r); 183 } 184 backend::copy(*r, *rh); 185 186 scalar_type eps = std::max(norm_rhs * prm.tol, prm.abstol); 187 scalar_type res = prm.check_after ? 2 * eps : norm(*r); 188 189 coef_type rho1 = zero; 190 coef_type rho2 = zero; 191 coef_type alpha = zero; 192 coef_type omega = zero; 193 194 size_t iter = 0; 195 for(bool first = true; res > eps && iter < prm.maxiter; ++iter) { 196 197 rho2 = rho1; 198 rho1 = inner_product(*r, *rh); 199 200 if (first) { 201 backend::copy(*r, *p); 202 first = false; 203 } else { 204 precondition(!math::is_zero(rho2), "Zero rho in BiCGStab"); 205 coef_type beta = (rho1 * alpha) / (rho2 * omega); 206 backend::axpbypcz(one, *r, -beta * omega, *v, beta, *p); 207 } 208 209 preconditioner::spmv(prm.pside, P, A, *p, *v, *T); 210 211 alpha = rho1 / inner_product(*rh, *v); 212 213 if (prm.pside == side::left) { 214 backend::axpby(alpha, *p, one, x); 215 } else { 216 backend::axpby(alpha, *T, one, x); 217 } 218 219 backend::axpbypcz(one, *r, -alpha, *v, zero, *s); 220 221 if ((res = norm(*s)) > eps) { 222 preconditioner::spmv(prm.pside, P, A, *s, *t, *T); 223 224 omega = inner_product(*t, *s) / inner_product(*t, *t); 225 226 precondition(!math::is_zero(omega), "Zero omega in BiCGStab"); 227 228 if (prm.pside == side::left) { 229 backend::axpby(omega, *s, one, x); 230 } else { 231 backend::axpby(omega, *T, one, x); 232 } 233 234 backend::axpbypcz(one, *s, -omega, *t, zero, *r); 235 236 res = norm(*r); 237 } 238 239 if (prm.verbose && iter % 5 == 0) 240 std::cout << iter << "\t" << std::scientific << res / norm_rhs << std::endl; 241 } 242 243 return std::make_tuple(iter, res / norm_rhs); 244 } 245 246 /* Computes the solution for the given right-hand side \p rhs. The 247 * system matrix is the same that was used for the setup of the 248 * preconditioner \p P. Returns the number of iterations made and the 249 * achieved residual as a ``std::tuple``. The solution vector \p x 250 * provides initial approximation in input and holds the computed 251 * solution on output. 252 */ 253 template <class Precond, class Vec1, class Vec2> operator ()(const Precond & P,const Vec1 & rhs,Vec2 && x) const254 std::tuple<size_t, scalar_type> operator()( 255 const Precond &P, const Vec1 &rhs, Vec2 &&x) const 256 { 257 return (*this)(P.system_matrix(), P, rhs, x); 258 } 259 bytes() const260 size_t bytes() const { 261 return 262 backend::bytes(*r) + 263 backend::bytes(*p) + 264 backend::bytes(*v) + 265 backend::bytes(*s) + 266 backend::bytes(*t) + 267 backend::bytes(*rh) + 268 backend::bytes(*T); 269 } 270 operator <<(std::ostream & os,const bicgstab & s)271 friend std::ostream& operator<<(std::ostream &os, const bicgstab &s) { 272 return os 273 << "Type: BiCGStab" 274 << "\nUnknowns: " << s.n 275 << "\nMemory footprint: " << human_readable_memory(s.bytes()) 276 << std::endl; 277 } 278 public: 279 params prm; 280 281 private: 282 size_t n; 283 284 std::shared_ptr<vector> r; 285 std::shared_ptr<vector> p; 286 std::shared_ptr<vector> v; 287 std::shared_ptr<vector> s; 288 std::shared_ptr<vector> t; 289 std::shared_ptr<vector> rh; 290 std::shared_ptr<vector> T; 291 292 InnerProduct inner_product; 293 294 template <class Vec> norm(const Vec & x) const295 scalar_type norm(const Vec &x) const { 296 return sqrt(math::norm(inner_product(x, x))); 297 } 298 }; 299 300 } // namespace solver 301 } // namespace amgcl 302 303 304 #endif 305