SolverCG and SolverGMRES compute quantities such as eigenvalues, condition numbers
and coefficients during the solution process. Add boost::signals2 to be able to
communicate these estimates back to the caller.
* are equal. But, even for small numbers of iteration steps, the condition
* number of @p T is a good estimate for the one of @p PA.
*
- * With the coefficients @p alpha and @p beta written to the log file if
- * <tt>AdditionalData::log_coefficients = true</tt>, the matrix @p T_m after
- * @p m steps is the tri-diagonal matrix with diagonal elements
+ * After @p m steps the matrix T_m can be written in terms of the coefficients
+ * @p alpha and @p beta as the tri-diagonal matrix with diagonal elements
* <tt>1/alpha_0</tt>, <tt>1/alpha_1 + beta_0/alpha_0</tt>, ...,
* <tt>1/alpha_{m-1</tt>+beta_{m-2}/alpha_{m-2}} and off-diagonal elements
* <tt>sqrt(beta_0)/alpha_0</tt>, ..., <tt>sqrt(beta_{m-2</tt>)/alpha_{m-2}}.
* @see Y. Saad: "Iterative methods for Sparse Linear Systems", section 6.7.3
* for details.
*
+ * The coefficients, eigenvalues and condition number (computed as the ratio of
+ * the largest over smallest eigenvalue) can be obtained by connecting a
+ * function as a slot to the solver using one of the functions
+ * @p connect_coefficients_slot, @p connect_eigenvalues_slot and
+ * @p connect_condition_number_slot. These slots will then be called from the
+ * solver with the estimates as argument.
+ *
+ * @deprecated Alternatively these estimates can be written to deallog by
+ * setting flags in @p AdditionalData.
*
* <h3>Observing the progress of linear solver iterations</h3>
*
/**
* Constructor. Initialize data fields. Confer the description of those.
+ * @deprecated Instead use: connect_coefficients_slot,
+ * connect_condition_number_slot, and connect_eigenvalues_slot.
*/
- AdditionalData (const bool log_coefficients = false,
+ AdditionalData (const bool log_coefficients,
const bool compute_condition_number = false,
const bool compute_all_condition_numbers = false,
- const bool compute_eigenvalues = false);
+ const bool compute_eigenvalues = false) DEAL_II_DEPRECATED;
+
+ /**
+ * Constructor. Initializes all data fields to false.
+ */
+ AdditionalData();
};
/**
const VECTOR &b,
const PRECONDITIONER &precondition);
+ /**
+ * Connect a slot to retrieve the CG coefficients. The slot will be called
+ * with alpha as the first argument and with beta as the second argument,
+ * where alpha and beta follow the notation in
+ * Y. Saad: "Iterative methods for Sparse Linear Systems", section 6.7.
+ * Called once per iteration
+ */
+ boost::signals2::connection
+ connect_coefficients_slot(
+ const std_cxx11::function<void (double,double)> &slot);
+
+ /**
+ * Connect a slot to retrieve the estimated condition number.
+ * Called on each iteration if every_iteration=true, otherwise called once
+ * when iterations are ended (i.e., either because convergence has been
+ * achieved, or because divergence has been detected).
+ */
+ boost::signals2::connection
+ connect_condition_number_slot(const std_cxx11::function<void (double)> &slot,
+ const bool every_iteration=false);
+
+ /**
+ * Connect a slot to retrieve the estimated eigenvalues.
+ * Called on each iteration if every_iteration=true, otherwise called once
+ * when iterations are ended (i.e., either because convergence has been
+ * achieved, or because divergence has been detected).
+ */
+ boost::signals2::connection
+ connect_eigenvalues_slot(
+ const std_cxx11::function<void (const std::vector<double> &)> &slot,
+ const bool every_iteration=false);
+
protected:
/**
* Implementation of the computation of the norm of the residual. This can
const VECTOR &r,
const VECTOR &d) const;
+ /**
+ * Estimates the eigenvalues from diagonal and offdiagonal. Uses
+ * these estimate to compute the condition number. Calls the signals
+ * eigenvalues_signal and cond_signal with these estimates as arguments.
+ * Outputs the eigenvalues/condition-number to deallog if
+ * log_eigenvalues/log_cond is true.
+ */
+ static void
+ compute_eigs_and_cond(
+ const std::vector<double> &diagonal,
+ const std::vector<double> &offdiagonal,
+ const boost::signals2::signal<void (const std::vector<double> &)> &eigenvalues_signal,
+ const boost::signals2::signal<void (double)> &cond_signal,
+ const bool log_eigenvalues,
+ const bool log_cond);
+
/**
* Temporary vectors, allocated through the @p VectorMemory object at the
* start of the actual solution process and deallocated at the end.
*/
AdditionalData additional_data;
+ /**
+ * Signal used to retrieve the CG coefficients.
+ * Called on each iteration.
+ */
+ boost::signals2::signal<void (double,double)> coefficients_signal;
+
+ /**
+ * Signal used to retrieve the estimated condition number.
+ * Called once when all iterations are ended.
+ */
+ boost::signals2::signal<void (double)> condition_number_signal;
+
+ /**
+ * Signal used to retrieve the estimated condition numbers.
+ * Called on each iteration.
+ */
+ boost::signals2::signal<void (double)> all_condition_numbers_signal;
+
+ /**
+ * Signal used to retrieve the estimated eigenvalues.
+ * Called once when all iterations are ended.
+ */
+ boost::signals2::signal<void (const std::vector<double> &)> eigenvalues_signal;
+
+ /**
+ * Signal used to retrieve the estimated eigenvalues.
+ * Called on each iteration.
+ */
+ boost::signals2::signal<void (const std::vector<double> &)> all_eigenvalues_signal;
+
private:
void cleanup();
};
+template <class VECTOR>
+inline
+SolverCG<VECTOR>::AdditionalData::
+AdditionalData ()
+ :
+ log_coefficients (false),
+ compute_condition_number(false),
+ compute_all_condition_numbers(false),
+ compute_eigenvalues(false)
+{}
+
+
+
template <class VECTOR>
SolverCG<VECTOR>::SolverCG (SolverControl &cn,
VectorMemory<VECTOR> &mem,
+template <class VECTOR>
+inline void
+SolverCG<VECTOR>::compute_eigs_and_cond(
+ const std::vector<double> &diagonal,
+ const std::vector<double> &offdiagonal,
+ const boost::signals2::signal<void (const std::vector<double> &)> &eigenvalues_signal,
+ const boost::signals2::signal<void (double)> &cond_signal,
+ const bool log_eigenvalues,
+ const bool log_cond)
+{
+ //Avoid computing eigenvalues unless they are needed.
+ if (!cond_signal.empty()|| !eigenvalues_signal.empty() || log_cond ||
+ log_eigenvalues)
+ {
+ TridiagonalMatrix<double> T(diagonal.size(), true);
+ for (size_type i=0; i<diagonal.size(); ++i)
+ {
+ T(i,i) = diagonal[i];
+ if (i< diagonal.size()-1)
+ T(i,i+1) = offdiagonal[i];
+ }
+ T.compute_eigenvalues();
+ //Need two eigenvalues to estimate the condition number.
+ if (diagonal.size()>1)
+ {
+ double condition_number=T.eigenvalue(T.n()-1)/T.eigenvalue(0);
+ cond_signal(condition_number);
+ //Send to deallog
+ if (log_cond)
+ {
+ deallog << "Condition number estimate: " <<
+ condition_number << std::endl;
+ }
+ }
+ //Avoid copying the eigenvalues of T to a vector unless a signal is
+ //connected.
+ if (!eigenvalues_signal.empty())
+ {
+ std::vector<double> eigenvalues(T.n());
+ for (unsigned int j = 0; j < T.n(); ++j)
+ {
+ eigenvalues.at(j)=T.eigenvalue(j);
+ }
+ eigenvalues_signal(eigenvalues);
+ }
+ if (log_eigenvalues)
+ {
+ for (size_type i=0; i<T.n(); ++i)
+ deallog << ' ' << T.eigenvalue(i);
+ deallog << std::endl;
+ }
+ }
+
+}
+
+
+
template <class VECTOR>
template <class MATRIX, class PRECONDITIONER>
void
Vp = this->memory.alloc();
// Should we build the matrix for
// eigenvalue computations?
- bool do_eigenvalues = additional_data.compute_condition_number
- | additional_data.compute_all_condition_numbers
- | additional_data.compute_eigenvalues;
+ const bool do_eigenvalues = !condition_number_signal.empty()
+ |!all_condition_numbers_signal.empty()
+ |!eigenvalues_signal.empty()
+ |!all_eigenvalues_signal.empty()
+ | additional_data.compute_condition_number
+ | additional_data.compute_all_condition_numbers
+ | additional_data.compute_eigenvalues;
double eigen_beta_alpha = 0;
// vectors used for eigenvalue
d.sadd(beta,-1.,g);
}
+ this->coefficients_signal(alpha,beta);
if (additional_data.log_coefficients)
deallog << "alpha-beta:" << alpha << '\t' << beta << std::endl;
// set up the vectors
eigen_beta_alpha = beta/alpha;
offdiagonal.push_back(std::sqrt(beta)/alpha);
}
-
- if (additional_data.compute_all_condition_numbers && (diagonal.size()>1))
- {
- TridiagonalMatrix<double> T(diagonal.size(), true);
- for (size_type i=0; i<diagonal.size(); ++i)
- {
- T(i,i) = diagonal[i];
- if (i< diagonal.size()-1)
- T(i,i+1) = offdiagonal[i];
- }
- T.compute_eigenvalues();
- deallog << "Condition number estimate: " <<
- T.eigenvalue(T.n()-1)/T.eigenvalue(0) << std::endl;
- }
+ compute_eigs_and_cond(diagonal,offdiagonal,all_eigenvalues_signal,
+ all_condition_numbers_signal,false,
+ additional_data.compute_all_condition_numbers);
}
}
catch (...)
cleanup();
throw;
}
-
- // Write eigenvalues or condition number
- if (do_eigenvalues)
- {
- TridiagonalMatrix<double> T(diagonal.size(), true);
- for (size_type i=0; i<diagonal.size(); ++i)
- {
- T(i,i) = diagonal[i];
- if (i< diagonal.size()-1)
- T(i,i+1) = offdiagonal[i];
- }
- T.compute_eigenvalues();
- if (additional_data.compute_condition_number
- && ! additional_data.compute_all_condition_numbers
- && (diagonal.size() > 1))
- deallog << "Condition number estimate: " <<
- T.eigenvalue(T.n()-1)/T.eigenvalue(0) << std::endl;
- if (additional_data.compute_eigenvalues)
- {
- for (size_type i=0; i<T.n(); ++i)
- deallog << ' ' << T.eigenvalue(i);
- deallog << std::endl;
- }
- }
+ compute_eigs_and_cond(diagonal,offdiagonal,eigenvalues_signal,
+ condition_number_signal,
+ additional_data.compute_eigenvalues,
+ (additional_data.compute_condition_number &&
+ !additional_data.compute_all_condition_numbers));
// Deallocate Memory
cleanup();
// otherwise exit as normal
}
+
+
+template<class VECTOR>
+boost::signals2::connection
+SolverCG<VECTOR>::connect_coefficients_slot(
+ const std_cxx11::function<void(double,double)> &slot)
+{
+ return coefficients_signal.connect(slot);
+}
+
+
+
+template<class VECTOR>
+boost::signals2::connection
+SolverCG<VECTOR>::connect_condition_number_slot(
+ const std_cxx11::function<void(double)> &slot,
+ const bool every_iteration)
+{
+ if (every_iteration)
+ {
+ return all_condition_numbers_signal.connect(slot);
+ }
+ else
+ {
+ return condition_number_signal.connect(slot);
+ }
+}
+
+
+
+template<class VECTOR>
+boost::signals2::connection
+SolverCG<VECTOR>::connect_eigenvalues_slot(
+ const std_cxx11::function<void (const std::vector<double> &)> &slot,
+ const bool every_iteration)
+{
+ if (every_iteration)
+ {
+ return all_eigenvalues_signal.connect(slot);
+ }
+ else
+ {
+ return eigenvalues_signal.connect(slot);
+ }
+}
+
+
+
#endif // DOXYGEN
DEAL_II_NAMESPACE_CLOSE
#include <vector>
#include <cmath>
+#include <algorithm>
DEAL_II_NAMESPACE_OPEN
* to observe the progress of the iteration.
*
*
+ * <h3>Eigenvalue and condition number estimates</h3>
+ *
+ * This class can estimate eigenvalues and condition number during the solution
+ * process. This is done by creating the Hessenberg matrix during the inner
+ * iterations. The eigenvalues are estimated as the eigenvalues of the
+ * Hessenberg matrix and the condition number is estimated as the ratio of the
+ * largest and smallest singular value of the Hessenberg matrix. The estimates
+ * can be obtained by connecting a function as a slot using
+ * @p connect_condition_number_slot and @p connect_eigenvalues_slot. These slots
+ * will then be called from the solver with the estimates as argument.
+ *
+ *
* @author Wolfgang Bangerth, Guido Kanschat, Ralf Hartmann.
*/
template <class VECTOR = Vector<double> >
AdditionalData (const unsigned int max_n_tmp_vectors = 30,
const bool right_preconditioning = false,
const bool use_default_residual = true,
- const bool force_re_orthogonalization = false,
- const bool compute_eigenvalues = false);
+ const bool force_re_orthogonalization = false);
+
+ /**
+ * Constructor.
+ * @deprecated To obtain the estimated eigenvalues instead use:
+ * connect_eigenvalues_slot
+ */
+ AdditionalData (const unsigned int max_n_tmp_vectors,
+ const bool right_preconditioning,
+ const bool use_default_residual,
+ const bool force_re_orthogonalization,
+ const bool compute_eigenvalues) DEAL_II_DEPRECATED;
/**
* Maximum number of temporary vectors. This parameter controls the size
const VECTOR &b,
const PRECONDITIONER &precondition);
+ /**
+ * Connect a slot to retrieve the estimated condition number.
+ * Called on each outer iteration if every_iteration=true, otherwise called
+ * once when iterations are ended (i.e., either because convergence has been
+ * achieved, or because divergence has been detected).
+ */
+ boost::signals2::connection
+ connect_condition_number_slot(const std_cxx11::function<void (double)> &slot,
+ const bool every_iteration=false);
+
+ /**
+ * Connect a slot to retrieve the estimated eigenvalues.
+ * Called on each outer iteration if every_iteration=true, otherwise called
+ * once when iterations are ended (i.e., either because convergence has been
+ * achieved, or because divergence has been detected).
+ */
+ boost::signals2::connection
+ connect_eigenvalues_slot(
+ const std_cxx11::function<void (const std::vector<std::complex<double> > &)> &slot,
+ const bool every_iteration=false);
+
+
DeclException1 (ExcTooFewTmpVectors,
int,
<< "The number of temporary vectors you gave ("
<< "any results, and much more for reasonable ones.");
protected:
+
/**
* Includes the maximum number of tmp vectors.
*/
AdditionalData additional_data;
+ /**
+ * Signal used to retrieve the estimated condition number.
+ * Called once when all iterations are ended.
+ */
+ boost::signals2::signal<void (double)> condition_number_signal;
+
+ /**
+ * Signal used to retrieve the estimated condition numbers.
+ * Called on each outer iteration.
+ */
+ boost::signals2::signal<void (double)> all_condition_numbers_signal;
+
+ /**
+ * Signal used to retrieve the estimated eigenvalues.
+ * Called once when all iterations are ended.
+ */
+ boost::signals2::signal<void (const std::vector<std::complex<double> > &)> eigenvalues_signal;
+
+ /**
+ * Signal used to retrieve the estimated eigenvalues.
+ * Called on each outer iteration.
+ */
+ boost::signals2::signal<void (const std::vector<std::complex<double> > &)> all_eigenvalues_signal;
+
+
/**
* Implementation of the computation of the norm of the residual.
*/
Vector<double> &h,
bool &re_orthogonalize);
+ /**
+ * Estimates the eigenvalues from the Hessenberg matrix, H_orig, generated
+ * during the inner iterations. Uses these estimate to compute the condition
+ * number. Calls the signals eigenvalues_signal and cond_signal with these
+ * estimates as arguments. Outputs the eigenvalues to deallog if
+ * log_eigenvalues is true.
+ */
+ static void
+ compute_eigs_and_cond(
+ const FullMatrix<double> &H_orig ,
+ const unsigned int dim,
+ const boost::signals2::signal<void (const std::vector<std::complex<double> > &)> &eigenvalues_signal,
+ const boost::signals2::signal<void(double)> &cond_signal,
+ const bool log_eigenvalues);
+
/**
* Projected system matrix
*/
+template <class VECTOR>
+inline
+SolverGMRES<VECTOR>::AdditionalData::
+AdditionalData (const unsigned int max_n_tmp_vectors,
+ const bool right_preconditioning,
+ const bool use_default_residual,
+ const bool force_re_orthogonalization)
+ :
+ max_n_tmp_vectors(max_n_tmp_vectors),
+ right_preconditioning(right_preconditioning),
+ use_default_residual(use_default_residual),
+ force_re_orthogonalization(force_re_orthogonalization),
+ compute_eigenvalues(false)
+{}
+
+
+
template <class VECTOR>
inline
SolverGMRES<VECTOR>::AdditionalData::
right_preconditioning(right_preconditioning),
use_default_residual(use_default_residual),
force_re_orthogonalization(force_re_orthogonalization),
- compute_eigenvalues (compute_eigenvalues)
+ compute_eigenvalues(compute_eigenvalues)
{}
+
template <class VECTOR>
SolverGMRES<VECTOR>::SolverGMRES (SolverControl &cn,
VectorMemory<VECTOR> &mem,
+template<class VECTOR>
+inline void
+SolverGMRES<VECTOR>::compute_eigs_and_cond(
+ const FullMatrix<double> &H_orig,
+ const unsigned int dim,
+ const boost::signals2::signal<void (const std::vector<std::complex<double> > &)> &eigenvalues_signal,
+ const boost::signals2::signal<void (double)> &cond_signal,
+ const bool log_eigenvalues)
+{
+ //Avoid copying the Hessenberg matrix if it isn't needed.
+ if (!eigenvalues_signal.empty() || !cond_signal.empty() || log_eigenvalues )
+ {
+ LAPACKFullMatrix<double> mat(dim,dim);
+ for (unsigned int i=0; i<dim; ++i)
+ for (unsigned int j=0; j<dim; ++j)
+ mat(i,j) = H_orig(i,j);
+ //Avoid computing eigenvalues if they are not needed.
+ if (!eigenvalues_signal.empty() || log_eigenvalues )
+ {
+ //Copy mat so that we can compute svd below. Necessary since
+ //compute_eigenvalues will leave mat in state LAPACKSupport::unusable.
+ LAPACKFullMatrix<double> mat_eig(mat);
+ mat_eig.compute_eigenvalues();
+ std::vector<std::complex<double> > eigenvalues(dim);
+ for (unsigned int i=0; i<mat_eig.n(); ++i)
+ eigenvalues[i] = mat_eig.eigenvalue(i);
+ eigenvalues_signal(eigenvalues);
+ if (log_eigenvalues)
+ {
+ //Sort eigenvalues for nicer output.
+ std::sort(eigenvalues.begin(), eigenvalues.end(),
+ internal::SolverGMRES::complex_less_pred);
+ deallog << "Eigenvalue estimate: ";
+ for (unsigned int i=0; i<mat_eig.n(); ++i)
+ deallog << ' ' << eigenvalues[i];
+ deallog << std::endl;
+ }
+ }
+ //Calculate condition number, avoid calculating the svd if a slot
+ //isn't connected. Need at least a 2-by-2 matrix to do the estimate.
+ if (!cond_signal.empty() && (mat.n()>1))
+ {
+ mat.compute_svd();
+ double condition_number=mat.singular_value(0)/mat.singular_value(mat.n()-1);
+ cond_signal(condition_number);
+ }
+ }
+}
+
+
+
template<class VECTOR>
template<class MATRIX, class PRECONDITIONER>
void
// restart
unsigned int accumulated_iterations = 0;
+ const bool do_eigenvalues=
+ !condition_number_signal.empty()
+ |!all_condition_numbers_signal.empty()
+ |!eigenvalues_signal.empty()
+ |!all_eigenvalues_signal.empty()
+ |additional_data.compute_eigenvalues;
// for eigenvalue computation, need to collect the Hessenberg matrix (before
// applying Givens rotations)
FullMatrix<double> H_orig;
- if (additional_data.compute_eigenvalues)
+ if (do_eigenvalues)
H_orig.reinit(n_tmp_vectors, n_tmp_vectors-1);
// matrix used for the orthogonalization process later
// for eigenvalues, get the resulting coefficients from the
// orthogonalization process
- if (additional_data.compute_eigenvalues)
+ if (do_eigenvalues)
for (unsigned int i=0; i<dim+1; ++i)
H_orig(i,inner_iteration) = h(i);
for (unsigned int j=0; j<dim; ++j)
H1(i,j) = H(i,j);
- // compute eigenvalues from the Hessenberg matrix generated during the
- // inner iterations
- if (additional_data.compute_eigenvalues)
- {
- LAPACKFullMatrix<double> mat(dim,dim);
- for (unsigned int i=0; i<dim; ++i)
- for (unsigned int j=0; j<dim; ++j)
- mat(i,j) = H_orig(i,j);
- mat.compute_eigenvalues();
- std::vector<std::complex<double> > eigenvalues(dim);
- for (unsigned int i=0; i<mat.n(); ++i)
- eigenvalues[i] = mat.eigenvalue(i);
-
- std::sort(eigenvalues.begin(), eigenvalues.end(),
- internal::SolverGMRES::complex_less_pred);
-
- deallog << "Eigenvalue estimate: ";
- for (unsigned int i=0; i<mat.n(); ++i)
- deallog << ' ' << eigenvalues[i];
- deallog << std::endl;
- }
+ compute_eigs_and_cond(H_orig,dim,all_eigenvalues_signal,
+ all_condition_numbers_signal,
+ additional_data.compute_eigenvalues);
H1.backward(h,gamma);
}
while (iteration_state == SolverControl::iterate);
+ compute_eigs_and_cond(H_orig,dim,eigenvalues_signal,condition_number_signal,
+ false);
if (!use_default_residual)
{
this->memory.free(r);
+template<class VECTOR>
+boost::signals2::connection
+SolverGMRES<VECTOR>::connect_condition_number_slot(
+ const std_cxx11::function<void(double)> &slot,
+ const bool every_iteration)
+{
+ if (every_iteration)
+ {
+ return all_condition_numbers_signal.connect(slot);
+ }
+ else
+ {
+ return condition_number_signal.connect(slot);
+ }
+}
+
+
+
+template<class VECTOR>
+boost::signals2::connection
+SolverGMRES<VECTOR>::connect_eigenvalues_slot(
+ const std_cxx11::function<void (const std::vector<std::complex<double> > &)> &slot,
+ const bool every_iteration)
+{
+ if (every_iteration)
+ {
+ return all_eigenvalues_signal.connect(slot);
+ }
+ else
+ {
+ return eigenvalues_signal.connect(slot);
+ }
+}
+
+
+
template<class VECTOR>
double
SolverGMRES<VECTOR>::criterion ()
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 1998 - 2014 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Connects slots to all signals in solver_cg and solver_gmres and writes all
+//output to deallog.
+
+
+#include "../tests.h"
+#include "testmatrix.h"
+#include <cmath>
+#include <complex>
+#include <fstream>
+#include <iomanip>
+#include <string>
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/vector_memory.h>
+#include <deal.II/lac/solver_control.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/solver_gmres.h>
+#include <deal.II/lac/precondition.h>
+
+
+void output_double_number(double input,const std::string &text)
+{
+ deallog<<text<< input<<std::endl;
+}
+
+void output_coefficients(double alpha,double beta)
+{
+ deallog<<"alpha: "<< alpha << " beta: "<< beta << std::endl;
+}
+
+template<class NUMBER>
+void output_eigenvalues(const std::vector<NUMBER> &eigenvalues,const std::string &text)
+{
+ deallog<< text;
+ for (unsigned int j = 0; j < eigenvalues.size(); ++j)
+ {
+ deallog<< ' ' << eigenvalues.at(j);
+ }
+ deallog << std::endl;
+}
+
+template<class SOLVER, class MATRIX, class VECTOR, class PRECONDITION>
+void
+check_solve( SOLVER &solver, const MATRIX &A,
+ VECTOR &u, VECTOR &f, const PRECONDITION &P)
+{
+ u = 0.;
+ f = 1.;
+ try
+ {
+ solver.solve(A,u,f,P);
+ }
+ catch (dealii::SolverControl::NoConvergence &e)
+ {
+ deallog << "Exception: " << e.get_exc_name() << std::endl;
+ }
+}
+
+
+
+int main()
+{
+ std::ofstream logfile("output");
+ deallog << std::setprecision(4);
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ SolverControl solver_control(100, 1.e-3);
+
+ unsigned int size=30;
+ unsigned int dim = (size-1)*(size-1);
+
+ // Make matrix
+ FDMatrix testproblem(size, size);
+ SparsityPattern structure(dim, dim, 5);
+ testproblem.five_point_structure(structure);
+ structure.compress();
+ SparseMatrix<double> A(structure);
+ testproblem.five_point(A);
+
+ Vector<double> f(dim);
+ Vector<double> u(dim);
+ f = 1.;
+
+ try
+ {
+ SolverCG<> solver_cg(solver_control);
+ //Attach all possible slots.
+ solver_cg.connect_coefficients_slot(&output_coefficients);
+ solver_cg.connect_condition_number_slot(
+ std_cxx11::bind(output_double_number,std_cxx11::_1,"Condition number: "),true);
+ solver_cg.connect_condition_number_slot(
+ std_cxx11::bind(output_double_number,std_cxx11::_1,"Final condition number: "));
+ solver_cg.connect_eigenvalues_slot(
+ std_cxx11::bind(output_eigenvalues<double>,std_cxx11::_1,"Eigenvalues: "),true);
+ solver_cg.connect_eigenvalues_slot(
+ std_cxx11::bind(output_eigenvalues<double>,std_cxx11::_1,"Final Eigenvalues: "));
+ solver_cg.solve(A,u,f,PreconditionIdentity());
+
+ u=0;
+ SolverGMRES<> solver_gmres(solver_control);
+ //Attach all possible slots.
+ solver_gmres.connect_condition_number_slot(
+ std_cxx11::bind(output_double_number,std_cxx11::_1,"Condition number: "),true);
+ solver_gmres.connect_condition_number_slot(
+ std_cxx11::bind(output_double_number,std_cxx11::_1,"Final condition number: "));
+ solver_gmres.connect_eigenvalues_slot(
+ std_cxx11::bind(output_eigenvalues<std::complex<double> >,std_cxx11::_1,"Eigenvalues: "),true);
+ solver_gmres.connect_eigenvalues_slot(
+ std_cxx11::bind(output_eigenvalues<std::complex<double> >,std_cxx11::_1,"Final Eigenvalues: "));
+ solver_gmres.solve(A,u,f,PreconditionIdentity());
+ }
+ catch (std::exception &e)
+ {
+ std::cerr << "Exception: " << e.what() << std::endl;
+ }
+
+}
+
--- /dev/null
+
+DEAL:cg::Starting value 29.00
+DEAL:cg::alpha: 7.250 beta: 6.750
+DEAL:cg::Eigenvalues: 0.1379
+DEAL:cg::alpha: 0.8749 beta: 0.8595
+DEAL:cg::Condition number: 29.00
+DEAL:cg::Eigenvalues: 0.07373 2.138
+DEAL:cg::alpha: 0.7611 beta: 1.026
+DEAL:cg::Condition number: 61.17
+DEAL:cg::Eigenvalues: 0.05344 1.186 3.269
+DEAL:cg::alpha: 0.4874 beta: 0.9388
+DEAL:cg::Condition number: 99.58
+DEAL:cg::Eigenvalues: 0.04473 0.8257 2.584 4.454
+DEAL:cg::alpha: 0.4926 beta: 0.7872
+DEAL:cg::Condition number: 154.6
+DEAL:cg::Eigenvalues: 0.03835 0.5853 1.901 3.408 5.931
+DEAL:cg::alpha: 0.4826 beta: 0.9174
+DEAL:cg::Condition number: 190.9
+DEAL:cg::Eigenvalues: 0.03430 0.4500 1.488 2.865 4.152 6.546
+DEAL:cg::alpha: 0.4614 beta: 0.7765
+DEAL:cg::Condition number: 222.3
+DEAL:cg::Eigenvalues: 0.03119 0.3550 1.172 2.327 3.506 5.277 6.935
+DEAL:cg::alpha: 0.4739 beta: 0.8809
+DEAL:cg::Condition number: 248.2
+DEAL:cg::Eigenvalues: 0.02893 0.2914 0.9558 1.935 3.052 4.057 5.895 7.180
+DEAL:cg::alpha: 0.4500 beta: 0.7691
+DEAL:cg::Condition number: 270.6
+DEAL:cg::Eigenvalues: 0.02716 0.2444 0.7908 1.617 2.606 3.576 4.940 6.324 7.352
+DEAL:cg::alpha: 0.4645 beta: 0.8436
+DEAL:cg::Condition number: 289.5
+DEAL:cg::Eigenvalues: 0.02581 0.2101 0.6685 1.374 2.250 3.186 4.018 5.501 6.632 7.474
+DEAL:cg::alpha: 0.4416 beta: 0.7523
+DEAL:cg::Condition number: 305.7
+DEAL:cg::Eigenvalues: 0.02475 0.1837 0.5730 1.179 1.949 2.805 3.628 4.739 5.912 6.863 7.565
+DEAL:cg::alpha: 0.4535 beta: 0.8013
+DEAL:cg::Condition number: 319.1
+DEAL:cg::Eigenvalues: 0.02393 0.1636 0.4988 1.025 1.704 2.483 3.287 4.000 5.239 6.229 7.040 7.634
+DEAL:cg::alpha: 0.4317 beta: 0.7212
+DEAL:cg::Condition number: 330.1
+DEAL:cg::Eigenvalues: 0.02329 0.1480 0.4398 0.8995 1.500 2.203 2.954 3.669 4.608 5.620 6.480 7.180 7.688
+DEAL:cg::alpha: 0.4389 beta: 0.7451
+DEAL:cg::Condition number: 338.9
+DEAL:cg::Eigenvalues: 0.02281 0.1359 0.3930 0.7982 1.333 1.967 2.663 3.368 3.993 5.055 5.929 6.682 7.291 7.731
+DEAL:cg::alpha: 0.4165 beta: 0.6656
+DEAL:cg::Condition number: 345.7
+DEAL:cg::Eigenvalues: 0.02246 0.1267 0.3557 0.7156 1.193 1.766 2.406 3.073 3.704 4.519 5.405 6.181 6.846 7.381 7.766
+DEAL:cg::alpha: 0.4162 beta: 0.6562
+DEAL:cg::Condition number: 350.8
+DEAL:cg::Eigenvalues: 0.02222 0.1198 0.3265 0.6488 1.078 1.597 2.184 2.810 3.436 3.994 4.921 5.699 6.391 6.983 7.456 7.795
+DEAL:cg::alpha: 0.3885 beta: 0.5568
+DEAL:cg::Condition number: 354.3
+DEAL:cg::Eigenvalues: 0.02207 0.1151 0.3044 0.5959 0.9836 1.455 1.993 2.576 3.174 3.738 4.462 5.246 5.946 6.568 7.097 7.518 7.818
+DEAL:cg::alpha: 0.3719 beta: 0.4876
+DEAL:cg::Condition number: 356.6
+DEAL:cg::Eigenvalues: 0.02198 0.1121 0.2894 0.5567 0.9102 1.341 1.835 2.376 2.942 3.504 4.009 4.829 5.526 6.159 6.720 7.194 7.571 7.838
+DEAL:cg::alpha: 0.3309 beta: 0.3590
+DEAL:cg::Condition number: 358.0
+DEAL:cg::Eigenvalues: 0.02194 0.1107 0.2807 0.5314 0.8581 1.254 1.708 2.209 2.739 3.277 3.780 4.443 5.138 5.770 6.345 6.851 7.279 7.616 7.856
+DEAL:cg::alpha: 0.3081 beta: 0.3265
+DEAL:cg::Condition number: 358.9
+DEAL:cg::Eigenvalues: 0.02193 0.1101 0.2769 0.5181 0.8261 1.193 1.612 2.073 2.567 3.078 3.582 4.046 4.774 5.401 5.979 6.503 6.963 7.350 7.655 7.870
+DEAL:cg::alpha: 0.3281 beta: 0.4916
+DEAL:cg::Condition number: 359.5
+DEAL:cg::Eigenvalues: 0.02193 0.1099 0.2751 0.5103 0.8018 1.136 1.507 1.920 2.374 2.854 3.342 3.801 4.403 5.032 5.612 6.147 6.630 7.053 7.407 7.686 7.882
+DEAL:cg::alpha: 0.4942 beta: 1.139
+DEAL:cg::Condition number: 360.0
+DEAL:cg::Eigenvalues: 0.02192 0.1097 0.2731 0.4961 0.7283 0.9649 1.297 1.703 2.151 2.624 3.106 3.579 4.005 4.676 5.257 5.798 6.296 6.744 7.134 7.459 7.714 7.893
+DEAL:cg::alpha: 0.5077 beta: 0.7756
+DEAL:cg::Condition number: 360.7
+DEAL:cg::Eigenvalues: 0.02192 0.1094 0.2667 0.3974 0.5549 0.8536 1.214 1.617 2.047 2.496 2.953 3.406 3.828 4.389 4.962 5.495 5.993 6.452 6.863 7.219 7.515 7.745 7.905
+DEAL:cg::alpha: 0.3174 beta: 0.3195
+DEAL:cg::Condition number: 361.3
+DEAL:cg::Eigenvalues: 0.02191 0.1092 0.2555 0.3298 0.5368 0.8410 1.198 1.589 1.999 2.417 2.838 3.263 3.681 4.082 4.676 5.200 5.695 6.157 6.581 6.962 7.292 7.565 7.775 7.918
+DEAL:cg::alpha: 0.3428 beta: 0.5580
+DEAL:cg::Condition number: 361.8
+DEAL:cg::Eigenvalues: 0.02191 0.1091 0.2468 0.3123 0.5328 0.8364 1.188 1.561 1.909 2.227 2.599 3.015 3.437 3.832 4.346 4.877 5.373 5.840 6.274 6.673 7.030 7.341 7.599 7.797 7.929
+DEAL:cg::alpha: 0.4717 beta: 0.7491
+DEAL:cg::Condition number: 362.3
+DEAL:cg::Eigenvalues: 0.02191 0.1090 0.2350 0.3000 0.5286 0.8253 1.084 1.255 1.619 2.030 2.452 2.872 3.286 3.684 4.057 4.611 5.100 5.561 5.996 6.399 6.766 7.095 7.384 7.627 7.815 7.940
+DEAL:cg::alpha: 0.3388 beta: 0.3950
+DEAL:cg::Condition number: 362.7
+DEAL:cg::Eigenvalues: 0.02191 0.1090 0.2247 0.2935 0.5238 0.7650 0.8817 1.210 1.599 2.006 2.410 2.796 3.168 3.535 3.880 4.367 4.845 5.300 5.733 6.141 6.519 6.859 7.159 7.422 7.648 7.827 7.948
+DEAL:cg::alpha: 0.3941 beta: 0.7196
+DEAL:cg::Condition number: 363.0
+DEAL:cg::Eigenvalues: 0.02191 0.1089 0.2179 0.2902 0.5179 0.6680 0.8533 1.202 1.587 1.962 2.224 2.520 2.908 3.304 3.685 4.034 4.554 5.012 5.447 5.860 6.248 6.609 6.934 7.216 7.456 7.665 7.836 7.954
+DEAL:cg::alpha: 0.4369 beta: 0.5564
+DEAL:cg::Condition number: 363.2
+DEAL:cg::Eigenvalues: 0.02191 0.1089 0.2098 0.2865 0.4941 0.5708 0.8420 1.185 1.408 1.625 2.023 2.432 2.829 3.206 3.564 3.891 4.348 4.792 5.216 5.621 6.005 6.364 6.700 7.008 7.277 7.497 7.683 7.844 7.959
+DEAL:cg::alpha: 0.2956 beta: 0.4380
+DEAL:cg::Condition number: 363.4
+DEAL:cg::Eigenvalues: 0.02191 0.1089 0.2053 0.2845 0.4628 0.5484 0.8375 1.147 1.267 1.608 2.013 2.416 2.789 3.117 3.433 3.760 4.096 4.563 4.984 5.387 5.774 6.138 6.475 6.783 7.069 7.327 7.536 7.702 7.850 7.962
+DEAL:cg::alpha: 0.4195 beta: 0.5047
+DEAL:cg::Condition number: 363.6
+DEAL:cg::Eigenvalues: 0.02191 0.1089 0.2013 0.2827 0.4285 0.5392 0.8312 1.030 1.217 1.596 1.966 2.127 2.453 2.844 3.219 3.571 3.891 4.326 4.752 5.158 5.548 5.919 6.273 6.603 6.899 7.157 7.393 7.599 7.745 7.863 7.968
+DEAL:cg::alpha: 0.3754 beta: 0.4222
+DEAL:cg::Condition number: 363.9
+DEAL:cg::Eigenvalues: 0.02191 0.1089 0.1988 0.2813 0.4042 0.5349 0.8174 0.9179 1.204 1.555 1.680 2.025 2.429 2.816 3.166 3.479 3.772 4.060 4.476 4.863 5.244 5.616 5.975 6.317 6.642 6.937 7.192 7.419 7.624 7.774 7.875 7.973
+DEAL:cg::alpha: 0.3396 beta: 0.4622
+DEAL:cg::Condition number: 364.0
+DEAL:cg::Eigenvalues: 0.02191 0.1088 0.1974 0.2804 0.3902 0.5326 0.7846 0.8676 1.195 1.432 1.619 2.015 2.410 2.694 2.903 3.243 3.585 3.895 4.315 4.721 5.100 5.446 5.759 6.058 6.362 6.668 6.957 7.207 7.429 7.631 7.784 7.881 7.975
+DEAL:cg::alpha: 0.3734 beta: 0.3948
+DEAL:cg::Condition number: 364.0
+DEAL:cg::Eigenvalues: 0.02191 0.1088 0.1965 0.2797 0.3787 0.5303 0.7245 0.8487 1.175 1.288 1.603 1.978 2.106 2.440 2.828 3.186 3.506 3.791 4.058 4.447 4.825 5.205 5.575 5.928 6.256 6.532 6.750 6.990 7.225 7.437 7.637 7.790 7.885 7.977
+DEAL:cg::alpha: 0.3032 beta: 0.3507
+DEAL:cg::Condition number: 364.0
+DEAL:cg::Eigenvalues: 0.02191 0.1088 0.1961 0.2793 0.3728 0.5288 0.6842 0.8432 1.132 1.234 1.594 1.825 2.032 2.430 2.810 3.115 3.339 3.621 3.910 4.317 4.694 5.015 5.308 5.633 5.971 6.303 6.614 6.873 7.063 7.252 7.446 7.640 7.794 7.887 7.977
+DEAL:cg::alpha: 0.3564 beta: 0.3733
+DEAL:cg::Condition number: 364.1
+DEAL:cg::Eigenvalues: 0.02191 0.1088 0.1959 0.2791 0.3690 0.5273 0.6520 0.8393 1.065 1.214 1.560 1.663 2.018 2.393 2.525 2.840 3.197 3.519 3.798 4.052 4.427 4.803 5.177 5.530 5.828 6.068 6.345 6.647 6.927 7.153 7.317 7.461 7.643 7.796 7.889 7.978
+DEAL:cg::alpha: 0.2812 beta: 0.2570
+DEAL:cg::Condition number: 364.1
+DEAL:cg::Eigenvalues: 0.02191 0.1088 0.1958 0.2790 0.3672 0.5263 0.6336 0.8366 1.012 1.205 1.475 1.618 2.005 2.190 2.442 2.825 3.171 3.450 3.687 3.940 4.335 4.684 4.956 5.250 5.593 5.932 6.235 6.447 6.680 6.950 7.188 7.382 7.493 7.647 7.798 7.889 7.978
+DEAL:cg::alpha: 0.3218 beta: 0.3291
+DEAL:cg::Condition number: 364.1
+DEAL:cg::Eigenvalues: 0.02191 0.1088 0.1958 0.2789 0.3664 0.5257 0.6233 0.8344 0.9748 1.199 1.393 1.607 1.947 2.052 2.430 2.791 2.941 3.216 3.532 3.807 4.050 4.415 4.787 5.148 5.445 5.682 5.979 6.297 6.580 6.762 6.975 7.206 7.413 7.545 7.653 7.799 7.890 7.978
+DEAL:cg::alpha: 0.2821 beta: 0.2217
+DEAL:cg::Condition number: 364.1
+DEAL:cg::Eigenvalues: 0.02191 0.1088 0.1958 0.2789 0.3660 0.5253 0.6175 0.8323 0.9462 1.192 1.322 1.598 1.804 2.024 2.402 2.518 2.834 3.186 3.485 3.728 3.962 4.346 4.672 4.910 5.218 5.562 5.863 6.072 6.333 6.629 6.870 7.024 7.222 7.427 7.590 7.663 7.801 7.891 7.978
+DEAL:cg::Convergence step 40 value 0.0009064
+DEAL:cg::Final condition number: 364.1
+DEAL:cg::Final Eigenvalues: 0.02191 0.1088 0.1958 0.2789 0.3660 0.5253 0.6175 0.8323 0.9462 1.192 1.322 1.598 1.804 2.024 2.402 2.518 2.834 3.186 3.485 3.728 3.962 4.346 4.672 4.910 5.218 5.562 5.863 6.072 6.333 6.629 6.870 7.024 7.222 7.427 7.590 7.663 7.801 7.891 7.978
+DEAL:GMRES::Starting value 29.00
+DEAL:GMRES::Eigenvalues: (7.954,0.000) (7.836,0.000) (7.665,0.000) (7.456,0.000) (7.216,0.000) (6.934,0.000) (6.609,0.000) (6.248,0.000) (5.860,0.000) (5.447,0.000) (5.012,0.000) (4.554,0.000) (0.02191,0.000) (0.1089,0.000) (0.2179,0.000) (0.2902,0.000) (0.5179,0.000) (0.6680,0.000) (0.8533,0.000) (1.202,0.000) (4.034,0.000) (3.685,0.000) (3.304,0.000) (1.587,0.000) (2.908,0.000) (2.520,0.000) (1.962,0.000) (2.224,0.000)
+DEAL:GMRES::Condition number: 363.0
+DEAL:GMRES::Eigenvalues: (7.978,0.000) (7.891,0.000) (7.797,0.000) (7.647,0.000) (7.393,0.000) (7.129,0.000) (6.901,0.000) (6.553,0.000) (6.146,0.000) (5.741,0.000) (5.333,0.000) (4.910,0.000) (4.475,0.000) (4.075,0.000) (3.694,0.000) (3.205,0.000) (2.760,0.000) (0.02765,0.000) (0.1418,0.000) (0.2014,0.000) (0.3484,0.000) (0.5334,0.000) (0.6447,0.000) (1.032,0.000) (1.319,0.000) (2.310,0.000) (1.634,0.000) (1.937,0.000)
+DEAL:GMRES::Condition number: 288.5
+DEAL:GMRES::Convergence step 65 value 0.0009799
+DEAL:GMRES::Eigenvalues: (0.02681,0.000) (0.2333,0.000) (1.178,0.000) (2.161,0.000) (3.554,0.000) (4.930,0.000) (7.726,0.000) (7.170,0.000) (6.307,0.000)
+DEAL:GMRES::Condition number: 288.2
+DEAL:GMRES::Final Eigenvalues: (0.02681,0.000) (0.2333,0.000) (1.178,0.000) (2.161,0.000) (3.554,0.000) (4.930,0.000) (7.726,0.000) (7.170,0.000) (6.307,0.000)
+DEAL:GMRES::Final condition number: 288.2
projects / dealii.git / commitdiff