Skip to content
Snippets Groups Projects
Commit 80a9182c authored by Christoph Grüninger's avatar Christoph Grüninger
Browse files

Merge branch '40-multirhstest-fails-to-compile-when-arpack-is-installed' into 'master'

Resolve "multirhstest fails to compile when arpack++ is installed"

Closes #40

See merge request core/dune-istl!161

(cherry picked from commit 9d012149)

3fcec40b Fix condition estimate for exotic number types
76b9f4e9 Fix condition estimate for exotic number types, clang edition
parent 961e9d9a
No related branches found
No related tags found
1 merge request!165Merge branch '40-multirhstest-fails-to-compile-when-arpack-is-installed' into 'master'
Pipeline #
......@@ -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
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment