Skip to content
Snippets Groups Projects
Commit 9d012149 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 !161
parents 0a5be313 76b9f4e9
No related branches found
No related tags found
1 merge request!161Resolve "multirhstest fails to compile when arpack++ is installed"
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.
Please register or to comment