1 #ifndef AMGCL_PRECONDITIONER_SCHUR_PRESSURE_CORRECTION_HPP 2 #define AMGCL_PRECONDITIONER_SCHUR_PRESSURE_CORRECTION_HPP 3 4 /* 5 The MIT License 6 7 Copyright (c) 2012-2021 Denis Demidov <dennis.demidov@gmail.com> 8 Copyright (c) 2016, Riccardo Rossi, CIMNE (International Center for Numerical Methods in Engineering) 9 10 Permission is hereby granted, free of charge, to any person obtaining a copy 11 of this software and associated documentation files (the "Software"), to deal 12 in the Software without restriction, including without limitation the rights 13 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell 14 copies of the Software, and to permit persons to whom the Software is 15 furnished to do so, subject to the following conditions: 16 17 The above copyright notice and this permission notice shall be included in 18 all copies or substantial portions of the Software. 19 20 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR 21 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 22 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE 23 AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 24 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, 25 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN 26 THE SOFTWARE. 27 */ 28 29 /** 30 * \file amgcl/preconditioner/schur_pressure_correction.hpp 31 * \author Denis Demidov <dennis.demidov@gmail.com> 32 * \brief Schur-complement pressure correction preconditioning scheme. 33 * 34 * [1] Elman, Howard, et al. "A taxonomy and comparison of parallel block 35 * multi-level preconditioners for the incompressible Navier–Stokes 36 * equations." Journal of Computational Physics 227.3 (2008): 1790-1808. 37 * [2] Gmeiner, Björn, et al. "A quantitative performance analysis for Stokes 38 * solvers at the extreme scale." arXiv preprint arXiv:1511.02134 (2015). 39 * [3] Vincent, C., and R. Boyer. "A preconditioned conjugate gradient 40 * Uzawa‐type method for the solution of the Stokes problem by mixed Q1–P0 41 * stabilized finite elements." International journal for numerical methods 42 * in fluids 14.3 (1992): 289-298. 43 */ 44 45 #include <vector> 46 47 #include <memory> 48 49 #include <amgcl/backend/builtin.hpp> 50 #include <amgcl/backend/detail/mixing.hpp> 51 #include <amgcl/util.hpp> 52 #include <amgcl/io/mm.hpp> 53 54 namespace amgcl { 55 namespace preconditioner { 56 57 /// Schur-complement pressure correction preconditioner 58 template <class USolver, class PSolver> 59 class schur_pressure_correction { 60 static_assert( 61 backend::backends_compatible< 62 typename USolver::backend_type, 63 typename PSolver::backend_type 64 >::value, 65 "Backends for pressure and flow preconditioners should coincide!" 66 ); 67 public: 68 typedef 69 typename backend::detail::common_scalar_backend< 70 typename USolver::backend_type, 71 typename PSolver::backend_type 72 >::type 73 backend_type; 74 75 typedef typename backend_type::value_type value_type; 76 typedef typename backend_type::matrix matrix; 77 typedef typename backend_type::vector vector; 78 typedef typename backend_type::params backend_params; 79 80 typedef typename backend::builtin<value_type>::matrix build_matrix; 81 82 struct params { 83 typedef typename USolver::params usolver_params; 84 typedef typename PSolver::params psolver_params; 85 86 usolver_params usolver; 87 psolver_params psolver; 88 89 std::vector<char> pmask; 90 91 // Variant of block preconditioner to use in apply() 92 // 1: schur pressure correction: 93 // S p = fp - Kpu Kuu^-1 fu 94 // Kuu u = fu - Kup p 95 // 2: Block triangular: 96 // S p = fp 97 // Kuu u = fu - Kup p 98 int type; 99 100 // Approximate Kuu^-1 with inverted diagonal of Kuu during 101 // construction of matrix-less Schur complement. 102 // When false, USolver is used instead. 103 bool approx_schur; 104 105 // Adjust preconditioner matrix for the Schur complement system. 106 // That is, use 107 // Kpp when adjust_p == 0, 108 // Kpp - dia(Kpu * dia(Kuu)^-1 * Kup) when adjust_p == 1, 109 // Kpp - Kpu * dia(Kuu)^-1 * Kup when adjust_p == 2 110 int adjust_p; 111 112 // Use 1/sum_j(abs(Kuu_{i,j})) instead of dia(Kuu)^-1 113 // as approximation for the Kuu^-1 (as in SIMPLEC algorithm) 114 bool simplec_dia; 115 116 int verbose; 117 paramsamgcl::preconditioner::schur_pressure_correction::params118 params() : type(1), approx_schur(false), adjust_p(1), simplec_dia(true), verbose(0) {} 119 120 #ifndef AMGCL_NO_BOOST paramsamgcl::preconditioner::schur_pressure_correction::params121 params(const boost::property_tree::ptree &p) 122 : AMGCL_PARAMS_IMPORT_CHILD(p, usolver), 123 AMGCL_PARAMS_IMPORT_CHILD(p, psolver), 124 AMGCL_PARAMS_IMPORT_VALUE(p, type), 125 AMGCL_PARAMS_IMPORT_VALUE(p, approx_schur), 126 AMGCL_PARAMS_IMPORT_VALUE(p, adjust_p), 127 AMGCL_PARAMS_IMPORT_VALUE(p, simplec_dia), 128 AMGCL_PARAMS_IMPORT_VALUE(p, verbose) 129 { 130 size_t n = 0; 131 132 n = p.get("pmask_size", n); 133 134 precondition(n > 0, 135 "Error in schur_complement parameters: " 136 "pmask_size is not set"); 137 138 if (p.count("pmask_pattern")) { 139 pmask.resize(n, 0); 140 141 std::string pattern = p.get("pmask_pattern", std::string()); 142 switch (pattern[0]) { 143 case '%': 144 { 145 int start = std::atoi(pattern.substr(1).c_str()); 146 int stride = std::atoi(pattern.substr(3).c_str()); 147 for(size_t i = start; i < n; i += stride) pmask[i] = 1; 148 } 149 break; 150 case '<': 151 { 152 size_t m = std::atoi(pattern.c_str()+1); 153 for(size_t i = 0; i < std::min(m, n); ++i) pmask[i] = 1; 154 } 155 break; 156 case '>': 157 { 158 size_t m = std::atoi(pattern.c_str()+1); 159 for(size_t i = m; i < n; ++i) pmask[i] = 1; 160 } 161 break; 162 default: 163 precondition(false, "Unknown pattern in pmask_pattern"); 164 } 165 } else if (p.count("pmask")) { 166 void *pm = 0; 167 pm = p.get("pmask", pm); 168 pmask.assign(static_cast<char*>(pm), static_cast<char*>(pm) + n); 169 } else { 170 precondition(false, 171 "Error in schur_complement parameters: " 172 "neither pmask_pattern, nor pmask is set" 173 ); 174 } 175 176 check_params(p, {"usolver", "psolver", "type", "approx_schur", "adjust_p", "simplec_dia", "pmask_size", "verbose"}, 177 {"pmask", "pmask_pattern"}); 178 } 179 getamgcl::preconditioner::schur_pressure_correction::params180 void get(boost::property_tree::ptree &p, const std::string &path = "") const 181 { 182 AMGCL_PARAMS_EXPORT_CHILD(p, path, usolver); 183 AMGCL_PARAMS_EXPORT_CHILD(p, path, psolver); 184 AMGCL_PARAMS_EXPORT_VALUE(p, path, type); 185 AMGCL_PARAMS_EXPORT_VALUE(p, path, approx_schur); 186 AMGCL_PARAMS_EXPORT_VALUE(p, path, adjust_p); 187 AMGCL_PARAMS_EXPORT_VALUE(p, path, simplec_dia); 188 AMGCL_PARAMS_EXPORT_VALUE(p, path, verbose); 189 } 190 #endif 191 } prm; 192 193 template <class Matrix> schur_pressure_correction(const Matrix & K,const params & prm=params (),const backend_params & bprm=backend_params ())194 schur_pressure_correction( 195 const Matrix &K, 196 const params &prm = params(), 197 const backend_params &bprm = backend_params() 198 ) 199 : prm(prm), n(backend::rows(K)), np(0), nu(0) 200 { 201 init(std::make_shared<build_matrix>(K), bprm); 202 } 203 schur_pressure_correction(std::shared_ptr<build_matrix> K,const params & prm=params (),const backend_params & bprm=backend_params ())204 schur_pressure_correction( 205 std::shared_ptr<build_matrix> K, 206 const params &prm = params(), 207 const backend_params &bprm = backend_params() 208 ) 209 : prm(prm), n(backend::rows(*K)), np(0), nu(0) 210 { 211 init(K, bprm); 212 } 213 214 template <class Vec1, class Vec2> apply(const Vec1 & rhs,Vec2 && x) const215 void apply(const Vec1 &rhs, Vec2 &&x) const { 216 backend::spmv(1, *x2u, rhs, 0, *rhs_u); 217 backend::spmv(1, *x2p, rhs, 0, *rhs_p); 218 219 if (prm.type == 1) { 220 // Kuu u = rhs_u 221 backend::clear(*u); 222 report("U1", (*U)(*rhs_u, *u)); 223 224 // rhs_p -= Kpu u 225 backend::spmv(-1, *Kpu, *u, 1, *rhs_p); 226 227 // S p = rhs_p 228 backend::clear(*p); 229 report("P1", (*P)(*this, *rhs_p, *p)); 230 231 // rhs_u -= Kup p 232 backend::spmv(-1, *Kup, *p, 1, *rhs_u); 233 234 // Kuu u = rhs_u 235 backend::clear(*u); 236 report("U2", (*U)(*rhs_u, *u)); 237 } else if (prm.type == 2) { 238 // S p = fp 239 backend::clear(*p); 240 report("P", (*P)(*this, *rhs_p, *p)); 241 242 // Kuu u = fu - Kup p 243 backend::spmv(-1, *Kup, *p, 1, *rhs_u); 244 backend::clear(*u); 245 report("U", (*U)(*rhs_u, *u)); 246 } 247 248 backend::spmv(1, *u2x, *u, 0, x); 249 backend::spmv(1, *p2x, *p, 1, x); 250 } 251 252 template <class Alpha, class Vec1, class Beta, class Vec2> spmv(Alpha alpha,const Vec1 & x,Beta beta,Vec2 & y) const253 void spmv(Alpha alpha, const Vec1 &x, Beta beta, Vec2 &y) const { 254 // y = beta y + alpha S x, where S = Kpp - Kpu Kuu^-1 Kup 255 if (prm.adjust_p == 1) { 256 backend::spmv( alpha, P->system_matrix(), x, beta, y); 257 backend::vmul( alpha, *Ld, x, 1, y); 258 } else if (prm.adjust_p == 2) { 259 backend::spmv( alpha, *Lm, x, beta, y); 260 } else { 261 backend::spmv( alpha, P->system_matrix(), x, beta, y); 262 } 263 264 backend::spmv(1, *Kup, x, 0, *tmp); 265 266 if (prm.approx_schur) { 267 backend::vmul(1, *M, *tmp, 0, *u); 268 } else { 269 backend::clear(*u); 270 (*U)(*tmp, *u); 271 } 272 273 backend::spmv(-alpha, *Kpu, *u, 1, y); 274 } 275 system_matrix_ptr() const276 std::shared_ptr<matrix> system_matrix_ptr() const { 277 return K; 278 } 279 system_matrix() const280 const matrix& system_matrix() const { 281 return *K; 282 } 283 bytes() const284 size_t bytes() const { 285 size_t b = 0; 286 287 b += backend::bytes(*K); 288 b += backend::bytes(*Kup); 289 b += backend::bytes(*Kpu); 290 b += backend::bytes(*x2u); 291 b += backend::bytes(*x2p); 292 b += backend::bytes(*u2x); 293 b += backend::bytes(*p2x); 294 b += backend::bytes(*rhs_u); 295 b += backend::bytes(*rhs_p); 296 b += backend::bytes(*u); 297 b += backend::bytes(*p); 298 b += backend::bytes(*tmp); 299 b += backend::bytes(*U); 300 b += backend::bytes(*P); 301 302 if (M) b += backend::bytes(*M); 303 if (Ld) b += backend::bytes(*Ld); 304 if (Lm) b += backend::bytes(*Lm); 305 306 return b; 307 } 308 309 private: 310 size_t n, np, nu; 311 312 std::shared_ptr<matrix> K, Lm, Kup, Kpu, x2u, x2p, u2x, p2x; 313 std::shared_ptr<vector> rhs_u, rhs_p, u, p, tmp; 314 std::shared_ptr<typename backend_type::matrix_diagonal> M; 315 std::shared_ptr<typename backend_type::matrix_diagonal> Ld; 316 317 std::shared_ptr<USolver> U; 318 std::shared_ptr<PSolver> P; 319 init(const std::shared_ptr<build_matrix> & K,const backend_params & bprm)320 void init(const std::shared_ptr<build_matrix> &K, const backend_params &bprm) 321 { 322 this->K = backend_type::copy_matrix(K, bprm); 323 324 // Extract matrix subblocks. 325 auto Kuu = std::make_shared<build_matrix>(); 326 auto Kpu = std::make_shared<build_matrix>(); 327 auto Kup = std::make_shared<build_matrix>(); 328 auto Kpp = std::make_shared<build_matrix>(); 329 330 std::vector<ptrdiff_t> idx(n); 331 332 for(size_t i = 0; i < n; ++i) 333 idx[i] = (prm.pmask[i] ? np++ : nu++); 334 335 Kuu->set_size(nu, nu, true); 336 Kup->set_size(nu, np, true); 337 Kpu->set_size(np, nu, true); 338 Kpp->set_size(np, np, true); 339 340 #pragma omp parallel for 341 for(ptrdiff_t i = 0; i < static_cast<ptrdiff_t>(n); ++i) { 342 ptrdiff_t ci = idx[i]; 343 char pi = prm.pmask[i]; 344 for(auto k = backend::row_begin(*K, i); k; ++k) { 345 char pj = prm.pmask[k.col()]; 346 347 if (pi) { 348 if (pj) { 349 ++Kpp->ptr[ci+1]; 350 } else { 351 ++Kpu->ptr[ci+1]; 352 } 353 } else { 354 if (pj) { 355 ++Kup->ptr[ci+1]; 356 } else { 357 ++Kuu->ptr[ci+1]; 358 } 359 } 360 } 361 } 362 363 Kuu->set_nonzeros(Kuu->scan_row_sizes()); 364 Kup->set_nonzeros(Kup->scan_row_sizes()); 365 Kpu->set_nonzeros(Kpu->scan_row_sizes()); 366 Kpp->set_nonzeros(Kpp->scan_row_sizes()); 367 368 #pragma omp parallel for 369 for(ptrdiff_t i = 0; i < static_cast<ptrdiff_t>(n); ++i) { 370 ptrdiff_t ci = idx[i]; 371 char pi = prm.pmask[i]; 372 373 ptrdiff_t uu_head = 0, up_head = 0, pu_head = 0, pp_head = 0; 374 375 if(pi) { 376 pu_head = Kpu->ptr[ci]; 377 pp_head = Kpp->ptr[ci]; 378 } else { 379 uu_head = Kuu->ptr[ci]; 380 up_head = Kup->ptr[ci]; 381 } 382 383 for(auto k = backend::row_begin(*K, i); k; ++k) { 384 ptrdiff_t j = k.col(); 385 value_type v = k.value(); 386 ptrdiff_t cj = idx[j]; 387 char pj = prm.pmask[j]; 388 389 if (pi) { 390 if (pj) { 391 Kpp->col[pp_head] = cj; 392 Kpp->val[pp_head] = v; 393 ++pp_head; 394 } else { 395 Kpu->col[pu_head] = cj; 396 Kpu->val[pu_head] = v; 397 ++pu_head; 398 } 399 } else { 400 if (pj) { 401 Kup->col[up_head] = cj; 402 Kup->val[up_head] = v; 403 ++up_head; 404 } else { 405 Kuu->col[uu_head] = cj; 406 Kuu->val[uu_head] = v; 407 ++uu_head; 408 } 409 } 410 } 411 } 412 413 if (prm.verbose >= 2) { 414 io::mm_write("Kuu.mtx", *Kuu); 415 io::mm_write("Kpp.mtx", *Kpp); 416 } 417 418 std::shared_ptr<backend::numa_vector<value_type>> Kuu_dia; 419 420 if (prm.simplec_dia) { 421 Kuu_dia = std::make_shared<backend::numa_vector<value_type>>(nu); 422 #pragma omp parallel for 423 for(ptrdiff_t i = 0; i < static_cast<ptrdiff_t>(nu); ++i) { 424 value_type s = math::zero<value_type>(); 425 for(ptrdiff_t j = Kuu->ptr[i], e = Kuu->ptr[i+1]; j < e; ++j) { 426 s += math::norm(Kuu->val[j]); 427 } 428 (*Kuu_dia)[i] = math::inverse(s); 429 } 430 } else { 431 Kuu_dia = diagonal(*Kuu, /*invert = */true); 432 } 433 434 if (prm.adjust_p == 1) { 435 // Use (Kpp - dia(Kpu * dia(Kuu)^-1 * Kup)) 436 // to setup the P preconditioner. 437 auto L = std::make_shared<backend::numa_vector<value_type>>(np, false); 438 439 #pragma omp parallel for 440 for(ptrdiff_t i = 0; i < static_cast<ptrdiff_t>(np); ++i) { 441 value_type s = math::zero<value_type>(); 442 for(ptrdiff_t j = Kpu->ptr[i], e = Kpu->ptr[i+1]; j < e; ++j) { 443 ptrdiff_t k = Kpu->col[j]; 444 value_type v = Kpu->val[j]; 445 for(ptrdiff_t jj = Kup->ptr[k], ee = Kup->ptr[k+1]; jj < ee; ++jj) { 446 if (Kup->col[jj] == i) { 447 s += v * (*Kuu_dia)[k] * Kup->val[jj]; 448 break; 449 } 450 } 451 } 452 453 (*L)[i] = s; 454 for(ptrdiff_t j = Kpp->ptr[i], e = Kpp->ptr[i+1]; j < e; ++j) { 455 if (Kpp->col[j] == i) { 456 Kpp->val[j] -= s; 457 break; 458 } 459 } 460 } 461 Ld = backend_type::copy_vector(L, bprm); 462 } else if (prm.adjust_p == 2) { 463 Lm = backend_type::copy_matrix(Kpp, bprm); 464 465 // Use (Kpp - Kpu * dia(Kuu)^-1 * Kup) 466 // to setup the P preconditioner. 467 backend::numa_vector<value_type> val(Kup->nnz); 468 469 #pragma omp parallel for 470 for(ptrdiff_t i = 0; i < static_cast<ptrdiff_t>(nu); ++i) { 471 value_type d = (*Kuu_dia)[i]; 472 for(ptrdiff_t j = Kup->ptr[i], e = Kup->ptr[i+1]; j < e; ++j) { 473 val[j] = d * Kup->val[j]; 474 } 475 } 476 477 build_matrix Kup_hat; 478 479 Kup_hat.own_data = false; 480 Kup_hat.nrows = nu; 481 Kup_hat.ncols = np; 482 Kup_hat.nnz = Kup->nnz; 483 Kup_hat.ptr = Kup->ptr; 484 Kup_hat.col = Kup->col; 485 Kup_hat.val = val.data(); 486 487 Kpp = backend::sum( 488 math::identity<value_type>(), *Kpp, 489 -math::identity<value_type>(), *backend::product(*Kpu, Kup_hat)); 490 } 491 492 U = std::make_shared<USolver>(*Kuu, prm.usolver, bprm); 493 P = std::make_shared<PSolver>(*Kpp, prm.psolver, bprm); 494 495 this->Kup = backend_type::copy_matrix(Kup, bprm); 496 this->Kpu = backend_type::copy_matrix(Kpu, bprm); 497 498 rhs_u = backend_type::create_vector(nu, bprm); 499 rhs_p = backend_type::create_vector(np, bprm); 500 501 u = backend_type::create_vector(nu, bprm); 502 p = backend_type::create_vector(np, bprm); 503 504 tmp = backend_type::create_vector(nu, bprm); 505 506 if (prm.approx_schur) 507 M = backend_type::copy_vector(Kuu_dia, bprm); 508 509 // Scatter/Gather matrices 510 auto x2u = std::make_shared<build_matrix>(); 511 auto x2p = std::make_shared<build_matrix>(); 512 auto u2x = std::make_shared<build_matrix>(); 513 auto p2x = std::make_shared<build_matrix>(); 514 515 x2u->set_size(nu, n, true); 516 x2p->set_size(np, n, true); 517 u2x->set_size(n, nu, true); 518 p2x->set_size(n, np, true); 519 520 { 521 ptrdiff_t x2u_head = 0, x2u_idx = 0; 522 ptrdiff_t x2p_head = 0, x2p_idx = 0; 523 ptrdiff_t u2x_head = 0, u2x_idx = 0; 524 ptrdiff_t p2x_head = 0, p2x_idx = 0; 525 526 for(size_t i = 0; i < n; ++i) { 527 if (prm.pmask[i]) { 528 x2p->ptr[++x2p_idx] = ++x2p_head; 529 ++p2x_head; 530 } else { 531 x2u->ptr[++x2u_idx] = ++x2u_head; 532 ++u2x_head; 533 } 534 535 p2x->ptr[++p2x_idx] = p2x_head; 536 u2x->ptr[++u2x_idx] = u2x_head; 537 } 538 } 539 540 x2u->set_nonzeros(); 541 x2p->set_nonzeros(); 542 u2x->set_nonzeros(); 543 p2x->set_nonzeros(); 544 545 { 546 ptrdiff_t x2u_head = 0; 547 ptrdiff_t x2p_head = 0; 548 ptrdiff_t u2x_head = 0; 549 ptrdiff_t p2x_head = 0; 550 551 for(size_t i = 0; i < n; ++i) { 552 ptrdiff_t j = idx[i]; 553 554 if (prm.pmask[i]) { 555 x2p->col[x2p_head] = i; 556 x2p->val[x2p_head] = math::identity<value_type>(); 557 ++x2p_head; 558 559 p2x->col[p2x_head] = j; 560 p2x->val[p2x_head] = math::identity<value_type>(); 561 ++p2x_head; 562 } else { 563 x2u->col[x2u_head] = i; 564 x2u->val[x2u_head] = math::identity<value_type>(); 565 ++x2u_head; 566 567 u2x->col[u2x_head] = j; 568 u2x->val[u2x_head] = math::identity<value_type>(); 569 ++u2x_head; 570 } 571 } 572 } 573 574 this->x2u = backend_type::copy_matrix(x2u, bprm); 575 this->x2p = backend_type::copy_matrix(x2p, bprm); 576 this->u2x = backend_type::copy_matrix(u2x, bprm); 577 this->p2x = backend_type::copy_matrix(p2x, bprm); 578 } 579 operator <<(std::ostream & os,const schur_pressure_correction & p)580 friend std::ostream& operator<<(std::ostream &os, const schur_pressure_correction &p) { 581 os << "Schur complement (two-stage preconditioner)" << std::endl; 582 os << " Unknowns: " << p.n << "(" << p.np << ")" << std::endl; 583 os << " Nonzeros: " << backend::nonzeros(p.system_matrix()) << std::endl; 584 os << " Memory: " << human_readable_memory(p.bytes()) << std::endl; 585 os << std::endl; 586 os << "[ U ]\n" << *p.U << std::endl; 587 os << "[ P ]\n" << *p.P << std::endl; 588 589 return os; 590 } 591 592 template <typename I, typename E> report(const std::string & name,const std::tuple<I,E> & c) const593 void report(const std::string &name, const std::tuple<I, E> &c) const { 594 if (prm.verbose >= 1) { 595 std::cout << name << " (" << std::get<0>(c) << ", " << std::get<1>(c) << ")\n"; 596 } 597 } 598 }; 599 600 } // namespace preconditioner 601 602 namespace backend { 603 604 template <class US, class PS, class Alpha, class Beta, class Vec1, class Vec2> 605 struct spmv_impl< Alpha, preconditioner::schur_pressure_correction<US, PS>, Vec1, Beta, Vec2> 606 { applyamgcl::backend::spmv_impl607 static void apply(Alpha alpha, const preconditioner::schur_pressure_correction<US, PS> &A, const Vec1 &x, Beta beta, Vec2 &y) 608 { 609 A.spmv(alpha, x, beta, y); 610 } 611 }; 612 613 template <class US, class PS, class Vec1, class Vec2, class Vec3> 614 struct residual_impl< preconditioner::schur_pressure_correction<US, PS>, Vec1, Vec2, Vec3> 615 { applyamgcl::backend::residual_impl616 static void apply(const Vec1 &rhs, const preconditioner::schur_pressure_correction<US, PS> &A, const Vec2 &x, Vec3 &r) 617 { 618 backend::copy(rhs, r); 619 A.spmv(-1, x, 1, r); 620 } 621 }; 622 623 } // namespace backend 624 } // namespace amgcl 625 626 #endif 627