diff --git a/dune/istl/solvers.hh b/dune/istl/solvers.hh index d49111eec505a6ca378aeae92bd49aac0a912651..1fabc5a2bba638add83a22392e11078ce47c15e4 100644 --- a/dune/istl/solvers.hh +++ b/dune/istl/solvers.hh @@ -443,52 +443,53 @@ namespace Dune { if (condition_estimate_) { #if HAVE_ARPACKPP + Hybrid::ifElse(enableConditionEstimate_t{}, [&](auto id) { - // Build T matrix which has extreme eigenvalues approximating - // those of the original system - // (see Y. Saad, Iterative methods for sparse linear systems) + // Build T matrix which has extreme eigenvalues approximating + // those of the original system + // (see Y. Saad, Iterative methods for sparse linear systems) - COND_MAT T(i, i, COND_MAT::row_wise); + COND_MAT T(i, i, COND_MAT::row_wise); - for (auto row = T.createbegin(); row != T.createend(); ++row) { - if (row.index() > 0) - row.insert(row.index()-1); - row.insert(row.index()); - if (row.index() < T.N() - 1) - row.insert(row.index()+1); - } - for (int row = 0; row < i; ++row) { - if (row > 0) { - T[row][row-1] = std::sqrt(betas[row-1]) / lambdas[row-1]; + for (auto row = T.createbegin(); row != T.createend(); ++row) { + if (row.index() > 0) + row.insert(row.index()-1); + row.insert(row.index()); + if (row.index() < T.N() - 1) + row.insert(row.index()+1); } - - T[row][row] = 1.0 / lambdas[row]; - if (row > 0) { - T[row][row] += betas[row-1] / lambdas[row-1]; - } - - if (row < i - 1) { - T[row][row+1] = std::sqrt(betas[row]) / lambdas[row]; + for (int row = 0; row < i; ++row) { + if (row > 0) { + T[row][row-1] = std::sqrt(id(betas[row-1])) / lambdas[row-1]; + } + + T[row][row] = 1.0 / id(lambdas[row]); + if (row > 0) { + T[row][row] += betas[row-1] / lambdas[row-1]; + } + + if (row < i - 1) { + T[row][row+1] = std::sqrt(id(betas[row])) / lambdas[row]; + } } - } - - // Compute largest and smallest eigenvalue of T matrix and return as estimate - Dune::ArPackPlusPlus_Algorithms<COND_MAT, COND_VEC> arpack(T); - real_type eps = 0.0; - COND_VEC eigv; - real_type min_eigv, max_eigv; - arpack.computeSymMinMagnitude (eps, eigv, min_eigv); - arpack.computeSymMaxMagnitude (eps, eigv, max_eigv); + // Compute largest and smallest eigenvalue of T matrix and return as estimate + Dune::ArPackPlusPlus_Algorithms<COND_MAT, COND_VEC> arpack(T); - res.condition_estimate = max_eigv / min_eigv; + real_type eps = 0.0; + COND_VEC eigv; + real_type min_eigv, max_eigv; + arpack.computeSymMinMagnitude (eps, eigv, min_eigv); + arpack.computeSymMaxMagnitude (eps, eigv, max_eigv); - if (_verbose > 0) { - std::cout << "Min eigv estimate: " << min_eigv << std::endl; - std::cout << "Max eigv estimate: " << max_eigv << std::endl; - std::cout << "Condition estimate: " << max_eigv / min_eigv << std::endl; - } + res.condition_estimate = id(max_eigv / min_eigv); + if (this->_verbose > 0) { + std::cout << "Min eigv estimate: " << min_eigv << std::endl; + std::cout << "Max eigv estimate: " << max_eigv << std::endl; + std::cout << "Condition estimate: " << max_eigv / min_eigv << std::endl; + } + }); #else std::cerr << "WARNING: Condition estimate was requested. This requires ARPACK, but ARPACK was not found!" << std::endl; #endif