From 447019ec1ec268754d3e0f2caab4cb1a3434b850 Mon Sep 17 00:00:00 2001
From: Linus Seelinger <S.Linus@gmx.de>
Date: Mon, 26 Jun 2017 12:00:41 +0200
Subject: [PATCH] Add condition number estimate to CG solver
---
dune/istl/eigenvalue/arpackpp.hh | 4 +-
dune/istl/solver.hh | 3 +
dune/istl/solvers.hh | 94 ++++++++++++++++++++++
dune/istl/test/CMakeLists.txt | 2 +
dune/istl/test/cgconditiontest.cc | 129 ++++++++++++++++++++++++++++++
5 files changed, 230 insertions(+), 2 deletions(-)
create mode 100644 dune/istl/test/cgconditiontest.cc
diff --git a/dune/istl/eigenvalue/arpackpp.hh b/dune/istl/eigenvalue/arpackpp.hh
index bee2fed15..7b9e6c9e5 100644
--- a/dune/istl/eigenvalue/arpackpp.hh
+++ b/dune/istl/eigenvalue/arpackpp.hh
@@ -308,7 +308,7 @@ namespace Dune
// allocate memory for variables, set parameters
const int nev = 1; // Number of eigenvalues to compute
- int ncv = std::min(20, nrows) // Number of Arnoldi vectors generated at each iteration (0 == auto)
+ int ncv = std::min(20, nrows); // Number of Arnoldi vectors generated at each iteration (0 == auto)
const Real tol = epsilon; // Stopping tolerance (relative accuracy of Ritz values) (0 == machine precision)
const int maxit = nIterationsMax_*nev; // Maximum number of Arnoldi update iterations allowed (0 == 100*nev)
Real* ev = new Real[nev]; // Computed eigenvalues of A
@@ -849,7 +849,7 @@ namespace Dune
// allocate memory for variables, set parameters
const int nev = 2; // Number of eigenvalues to compute
- const int ncv = std::min(20, nrows); // Number of Arnoldi vectors generated at each iteration (0 == auto)
+ int ncv = std::min(20, nrows); // Number of Arnoldi vectors generated at each iteration (0 == auto)
const Real tol = epsilon; // Stopping tolerance (relative accuracy of Ritz values) (0 == machine precision)
const int maxit = nIterationsMax_*nev; // Maximum number of Arnoldi update iterations allowed (0 == 100*nev)
Real* ev = new Real[nev]; // Computed eigenvalues of A^T*A
diff --git a/dune/istl/solver.hh b/dune/istl/solver.hh
index 672710c6a..1eefa7360 100644
--- a/dune/istl/solver.hh
+++ b/dune/istl/solver.hh
@@ -66,6 +66,9 @@ namespace Dune
/** \brief Convergence rate (average reduction per step) */
double conv_rate;
+ /** \brief Estimate of condition number */
+ double condition_estimate = -1;
+
/** \brief Elapsed time in seconds */
double elapsed;
};
diff --git a/dune/istl/solvers.hh b/dune/istl/solvers.hh
index ebc0c6a04..08b2860f1 100644
--- a/dune/istl/solvers.hh
+++ b/dune/istl/solvers.hh
@@ -26,6 +26,8 @@
#include <dune/common/timer.hh>
#include <dune/common/ftraits.hh>
#include <dune/common/typetraits.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/eigenvalue/arpackpp.hh>
namespace Dune {
/** @addtogroup ISTL_Solvers
@@ -254,6 +256,29 @@ namespace Dune {
// copy base class constructors
using IterativeSolver<X,X>::IterativeSolver;
+ /*!
+ \brief Constructor to initialize a CG solver.
+ \copydetails IterativeSolver::IterativeSolver(LinearOperator<X,Y>&, Preconditioner<X,Y>&, real_type, int, int)
+ \param condition_estimate Whether to calculate an estimate of the condition number.
+ The estimate is given in the InverseOperatorResult returned by apply().
+ */
+ CGSolver (LinearOperator<X,X>& op, Preconditioner<X,X>& prec,
+ real_type reduction, int maxit, int verbose, bool condition_estimate) : IterativeSolver<X,X>(op, prec, reduction, maxit, verbose),
+ condition_estimate_(condition_estimate)
+ {}
+
+ /*!
+ \brief Constructor to initialize a CG solver.
+ \copydetails IterativeSolver::IterativeSolver(LinearOperator<X,Y>&, ScalarProduct<X>&, Preconditioner<X,Y>&, real_type, int, int)
+ \param condition_estimate Whether to calculate an estimate of the condition number.
+ The estimate is given in the InverseOperatorResult returned by apply().
+ */
+ CGSolver (LinearOperator<X,X>& op, ScalarProduct<X>& sp, Preconditioner<X,X>& prec,
+ real_type reduction, int maxit, int verbose, bool condition_estimate) : IterativeSolver<X,X>(op, sp, prec, reduction, maxit, verbose),
+ condition_estimate_(condition_estimate)
+ {}
+
+
/*!
\brief Apply inverse operator.
@@ -311,6 +336,10 @@ namespace Dune {
}
}
+ // Remember lambda and beta values for condition estimate
+ std::vector<real_type> lambdas(0);
+ std::vector<real_type> betas(0);
+
// some local variables
real_type def=def0; // loop variables
field_type rho,rholast,lambda,alpha,beta;
@@ -328,6 +357,8 @@ namespace Dune {
_op->apply(p,q); // q=Ap
alpha = _sp->dot(p,q); // scalar product
lambda = rholast/alpha; // minimization
+ if (condition_estimate_)
+ lambdas.push_back(std::real(lambda));
x.axpy(lambda,p); // update solution
b.axpy(-lambda,q); // update defect
@@ -358,6 +389,8 @@ namespace Dune {
_prec->apply(q,b); // apply preconditioner
rho = _sp->dot(q,b); // orthogonalization
beta = rho/rholast; // scaling factor
+ if (condition_estimate_)
+ betas.push_back(std::real(beta));
p *= beta; // scale old search direction
p += q; // orthogonalization with correction
rholast = rho; // remember rho for recurrence
@@ -382,8 +415,69 @@ namespace Dune {
<< ", TIT=" << res.elapsed/i
<< ", IT=" << i << std::endl;
}
+
+
+ if (condition_estimate_) {
+#if HAVE_ARPACKPP
+
+ // 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);
+
+ 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];
+ }
+
+ 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];
+ }
+ }
+
+ // 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);
+
+ res.condition_estimate = max_eigv / min_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;
+ }
+
+#else
+ std::cerr << "WARNING: Condition estimate was requested. This requires ARPACK, but ARPACK was not found!" << std::endl;
+#endif
+ }
}
+ private:
+ bool condition_estimate_ = false;
+
+ // Matrix and vector types used for condition estimate
+ typedef Dune::BCRSMatrix<Dune::FieldMatrix<real_type,1,1> > COND_MAT;
+ typedef Dune::BlockVector<Dune::FieldVector<real_type,1> > COND_VEC;
+
protected:
using IterativeSolver<X,X>::_op;
using IterativeSolver<X,X>::_prec;
diff --git a/dune/istl/test/CMakeLists.txt b/dune/istl/test/CMakeLists.txt
index a78ee0b83..8dd8970b3 100644
--- a/dune/istl/test/CMakeLists.txt
+++ b/dune/istl/test/CMakeLists.txt
@@ -8,6 +8,8 @@ dune_add_test(SOURCES bcrsassigntest.cc)
dune_add_test(SOURCES bcrsnormtest.cc)
+dune_add_test(SOURCES cgconditiontest.cc)
+
dune_add_test(SOURCES dotproducttest.cc)
dune_add_test(SOURCES complexmatrixtest.cc)
diff --git a/dune/istl/test/cgconditiontest.cc b/dune/istl/test/cgconditiontest.cc
new file mode 100644
index 000000000..9919337c9
--- /dev/null
+++ b/dune/istl/test/cgconditiontest.cc
@@ -0,0 +1,129 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/istl/preconditioners.hh>
+#include <dune/istl/solvers.hh>
+
+typedef Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> > MAT;
+typedef Dune::BlockVector<Dune::FieldVector<double,1> > VEC;
+
+void condition_test (MAT& T) {
+
+ int N = T.N();
+
+ // CG run (noop preconditioner) with condition estimate
+
+ auto prec = std::make_shared<Dune::Richardson<VEC, VEC>> (1.0);
+
+ auto op = std::make_shared<Dune::MatrixAdapter<MAT, VEC, VEC>>(T);
+
+ VEC d (N,N);
+ VEC v (N,N);
+
+ v = 1.0;
+ d = 1.0;
+
+ int verbosity = 2;
+
+ Dune::InverseOperatorResult result;
+ auto solver = std::make_shared<Dune::CGSolver<VEC> >(*op,*prec,1E-6,1000,verbosity,true);
+ solver->apply(v,d,result);
+
+
+
+#if HAVE_ARPACKPP
+ // Actual condition number using ARPACK
+ Dune::ArPackPlusPlus_Algorithms<MAT, VEC> arpack(T);
+
+ double eps = 0.0;
+ VEC eigv;
+ double min_eigv, max_eigv;
+ arpack.computeSymMinMagnitude (eps, eigv, min_eigv);
+ arpack.computeSymMaxMagnitude (eps, eigv, max_eigv);
+
+ std::cout << "Actual condition number from ARPACK: " << max_eigv / min_eigv << std::endl;
+#endif
+}
+
+
+
+int main(int argc, char **argv)
+{
+
+ // Simple stencil
+ {
+ int N = 142;
+
+ MAT T(N, N, 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 < N; ++row) {
+ T[row][row] = 2.0;
+ if (row > 0)
+ T[row][row-1] = -1.0;
+ if (row < N-1)
+ T[row][row+1] = -1.0;
+ }
+
+ condition_test (T);
+ }
+
+
+ // Ill-conditioned stencil
+ {
+ int N = 142;
+
+ MAT T(N, N, 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 < N; ++row) {
+ double factor = 0.8 + 0.2 * ((double)row / N);
+ T[row][row] = 2.0 * factor;
+ if (row > 0)
+ T[row][row-1] = -1.0 * factor;
+ if (row < N-1)
+ T[row][row+1] = -1.0 * factor;
+ }
+
+ condition_test (T);
+ }
+
+
+ // Very ill-conditioned stencil
+ {
+ int N = 40;
+
+ MAT T(N, N, MAT::row_wise);
+
+ for (auto row = T.createbegin(); row != T.createend(); ++row) {
+ for (int i = -10; i <= 10; i++) {
+ if (row.index()+i >= 0 && row.index()+i < T.N())
+ row.insert(row.index()+i);
+ }
+ }
+ for (int row = 0; row < N; ++row) {
+ double factor = 0.8 + 0.2 * ((double)row / N);
+ for (int i = -10; i <= 10; i++) {
+ if (row+i >= 0 && row+i < T.N())
+ T[row][row+i] = -1.0 * factor;
+ }
+ T[row][row] = 20.0 * factor;
+ }
+
+ condition_test (T);
+ }
+
+}
--
GitLab