<commit_before>// ---------------------------------------------------------------------
//
// Copyright (C) 1998 - 2016 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.
//
// ---------------------------------------------------------------------

#ifndef dealii__solver_cg_h
#define dealii__solver_cg_h


#include <deal.II/base/config.h>
#include <deal.II/lac/tridiagonal_matrix.h>
#include <deal.II/lac/solver.h>
#include <deal.II/lac/solver_control.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/subscriptor.h>
#include <cmath>

DEAL_II_NAMESPACE_OPEN

// forward declaration
class PreconditionIdentity;


/*!@addtogroup Solvers */
/*@{*/

/**
 * Preconditioned cg method for symmetric positive definite matrices. This
 * class is used first in step-3 and step-4, but is used in many other
 * tutorial programs as well. Like all other solver classes, it can work on
 * any kind of vector and matrix as long as they satisfy certain requirements
 * (for the requirements on matrices and vectors in order to work with this
 * class, see the documentation of the Solver base class). The type of the
 * solution vector must be passed as template argument, and defaults to
 * dealii::Vector<double>.
 *
 * Like all other solver classes, this class has a local structure called @p
 * AdditionalData which is used to pass additional parameters to the solver.
 * For this class, there is (among other things) a switch allowing for
 * additional output for the computation of eigenvalues of the matrix.
 *
 * @note This version of CG is taken from <NAME>'s book "Finite Elements".
 * It requires a symmetric preconditioner (i.e., for example, SOR is not a
 * possible choice).
 *
 *
 * <h3>Eigenvalue computation</h3>
 *
 * The cg-method performs an orthogonal projection of the original
 * preconditioned linear system to another system of smaller dimension.
 * Furthermore, the projected matrix @p T is tri-diagonal. Since the
 * projection is orthogonal, the eigenvalues of @p T approximate those of the
 * original preconditioned matrix @p PA. In fact, after @p n steps, where @p n
 * is the dimension of the original system, the eigenvalues of both matrices
 * 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.
 *
 * 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}}.
 * The eigenvalues of this matrix can be computed by postprocessing.
 *
 * @see <NAME>: "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>
 *
 * The solve() function of this class uses the mechanism described in the
 * Solver base class to determine convergence. This mechanism can also be used
 * to observe the progress of the iteration.
 *
 *
 * @author <NAME>, <NAME>, <NAME> and <NAME>
 */
template <typename VectorType = Vector<double> >
class SolverCG : public Solver<VectorType>
{
public:
  /**
   * Declare type for container size.
   */
  typedef types::global_dof_index size_type;

  /**
   * Standardized data struct to pipe additional data to the solver.
   */
  struct AdditionalData
  {
    /**
     * Write coefficients alpha and beta to the log file for later use in
     * eigenvalue estimates.
     */
    bool log_coefficients;

    /**
     * Compute the condition number of the projected matrix.
     *
     * @note Requires LAPACK support.
     */
    bool compute_condition_number;

    /**
     * Compute the condition number of the projected matrix in each step.
     *
     * @note Requires LAPACK support.
     */
    bool compute_all_condition_numbers;

    /**
     * Compute all eigenvalues of the projected matrix.
     *
     * @note Requires LAPACK support.
     */
    bool compute_eigenvalues;

    /**
     * Constructor. Initialize data fields.  Confer the description of those.
     * @deprecated Instead use: connect_coefficients_slot,
     * connect_condition_number_slot, and connect_eigenvalues_slot.
     */
    explicit
    AdditionalData (const bool log_coefficients,
                    const bool compute_condition_number = false,
                    const bool compute_all_condition_numbers = false,
                    const bool compute_eigenvalues = false) DEAL_II_DEPRECATED;

    /**
     * Constructor. Initializes all data fields to false.
     */
    AdditionalData();
  };

  /**
   * Constructor.
   */
  SolverCG (SolverControl            &cn,
            VectorMemory<VectorType> &mem,
            const AdditionalData     &data = AdditionalData());

  /**
   * Constructor. Use an object of type GrowingVectorMemory as a default to
   * allocate memory.
   */
  SolverCG (SolverControl        &cn,
            const AdditionalData &data=AdditionalData());

  /**
   * Virtual destructor.
   */
  virtual ~SolverCG ();

  /**
   * Solve the linear system $Ax=b$ for x.
   */
  template <typename MatrixType, typename PreconditionerType>
  void
  solve (const MatrixType         &A,
         VectorType               &x,
         const VectorType         &b,
         const PreconditionerType &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
   * be replaced by a more problem oriented functional in a derived class.
   */
  virtual double criterion();

  /**
   * Interface for derived class. This function gets the current iteration
   * vector, the residual and the update vector in each step. It can be used
   * for a graphical output of the convergence history.
   */
  virtual void print_vectors(const unsigned int step,
                             const VectorType   &x,
                             const VectorType   &r,
                             const VectorType   &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.
   */
  VectorType *Vr;
  VectorType *Vp;
  VectorType *Vz;

  /**
   * Within the iteration loop, the square of the residual vector is stored in
   * this variable. The function @p criterion uses this variable to compute
   * the convergence value, which in this class is the norm of the residual
   * vector and thus the square root of the @p res2 value.
   */
  double res2;

  /**
   * Additional parameters.
   */
  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();
};

/*@}*/

/*------------------------- Implementation ----------------------------*/

#ifndef DOXYGEN

template <typename VectorType>
inline
SolverCG<VectorType>::AdditionalData::
AdditionalData (const bool log_coefficients,
                const bool compute_condition_number,
                const bool compute_all_condition_numbers,
                const bool compute_eigenvalues)
  :
  log_coefficients (log_coefficients),
  compute_condition_number(compute_condition_number),
  compute_all_condition_numbers(compute_all_condition_numbers),
  compute_eigenvalues(compute_eigenvalues)
{}



template <typename VectorType>
inline
SolverCG<VectorType>::AdditionalData::
AdditionalData ()
  :
  log_coefficients (false),
  compute_condition_number(false),
  compute_all_condition_numbers(false),
  compute_eigenvalues(false)
{}



template <typename VectorType>
SolverCG<VectorType>::SolverCG (SolverControl        &cn,
                                VectorMemory<VectorType> &mem,
                                const AdditionalData     &data)
  :
  Solver<VectorType>(cn,mem),
  additional_data(data)
{}



template <typename VectorType>
SolverCG<VectorType>::SolverCG (SolverControl        &cn,
                                const AdditionalData &data)
  :
  Solver<VectorType>(cn),
  additional_data(data)
{}



template <typename VectorType>
SolverCG<VectorType>::~SolverCG ()
{}



template <typename VectorType>
double
SolverCG<VectorType>::criterion()
{
  return std::sqrt(res2);
}



template <typename VectorType>
void
SolverCG<VectorType>::cleanup()
{
  this->memory.free(Vr);
  this->memory.free(Vp);
  this->memory.free(Vz);
  deallog.pop();
}



template <typename VectorType>
void
SolverCG<VectorType>::print_vectors(const unsigned int,<commit_msg>Remove SolverCG::res2.<commit_after>// ---------------------------------------------------------------------
//
// Copyright (C) 1998 - 2016 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.
//
// ---------------------------------------------------------------------

#ifndef dealii__solver_cg_h
#define dealii__solver_cg_h


#include <deal.II/base/config.h>
#include <deal.II/lac/tridiagonal_matrix.h>
#include <deal.II/lac/solver.h>
#include <deal.II/lac/solver_control.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/subscriptor.h>
#include <cmath>

DEAL_II_NAMESPACE_OPEN

// forward declaration
class PreconditionIdentity;


/*!@addtogroup Solvers */
/*@{*/

/**
 * Preconditioned cg method for symmetric positive definite matrices. This
 * class is used first in step-3 and step-4, but is used in many other
 * tutorial programs as well. Like all other solver classes, it can work on
 * any kind of vector and matrix as long as they satisfy certain requirements
 * (for the requirements on matrices and vectors in order to work with this
 * class, see the documentation of the Solver base class). The type of the
 * solution vector must be passed as template argument, and defaults to
 * dealii::Vector<double>.
 *
 * Like all other solver classes, this class has a local structure called @p
 * AdditionalData which is used to pass additional parameters to the solver.
 * For this class, there is (among other things) a switch allowing for
 * additional output for the computation of eigenvalues of the matrix.
 *
 * @note This version of CG is taken from <NAME>'s book "Finite Elements".
 * It requires a symmetric preconditioner (i.e., for example, SOR is not a
 * possible choice).
 *
 *
 * <h3>Eigenvalue computation</h3>
 *
 * The cg-method performs an orthogonal projection of the original
 * preconditioned linear system to another system of smaller dimension.
 * Furthermore, the projected matrix @p T is tri-diagonal. Since the
 * projection is orthogonal, the eigenvalues of @p T approximate those of the
 * original preconditioned matrix @p PA. In fact, after @p n steps, where @p n
 * is the dimension of the original system, the eigenvalues of both matrices
 * 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.
 *
 * 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}}.
 * The eigenvalues of this matrix can be computed by postprocessing.
 *
 * @see <NAME>: "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>
 *
 * The solve() function of this class uses the mechanism described in the
 * Solver base class to determine convergence. This mechanism can also be used
 * to observe the progress of the iteration.
 *
 *
 * @author <NAME>, <NAME>, <NAME> and <NAME>
 */
template <typename VectorType = Vector<double> >
class SolverCG : public Solver<VectorType>
{
public:
  /**
   * Declare type for container size.
   */
  typedef types::global_dof_index size_type;

  /**
   * Standardized data struct to pipe additional data to the solver.
   */
  struct AdditionalData
  {
    /**
     * Write coefficients alpha and beta to the log file for later use in
     * eigenvalue estimates.
     */
    bool log_coefficients;

    /**
     * Compute the condition number of the projected matrix.
     *
     * @note Requires LAPACK support.
     */
    bool compute_condition_number;

    /**
     * Compute the condition number of the projected matrix in each step.
     *
     * @note Requires LAPACK support.
     */
    bool compute_all_condition_numbers;

    /**
     * Compute all eigenvalues of the projected matrix.
     *
     * @note Requires LAPACK support.
     */
    bool compute_eigenvalues;

    /**
     * Constructor. Initialize data fields.  Confer the description of those.
     * @deprecated Instead use: connect_coefficients_slot,
     * connect_condition_number_slot, and connect_eigenvalues_slot.
     */
    explicit
    AdditionalData (const bool log_coefficients,
                    const bool compute_condition_number = false,
                    const bool compute_all_condition_numbers = false,
                    const bool compute_eigenvalues = false) DEAL_II_DEPRECATED;

    /**
     * Constructor. Initializes all data fields to false.
     */
    AdditionalData();
  };

  /**
   * Constructor.
   */
  SolverCG (SolverControl            &cn,
            VectorMemory<VectorType> &mem,
            const AdditionalData     &data = AdditionalData());

  /**
   * Constructor. Use an object of type GrowingVectorMemory as a default to
   * allocate memory.
   */
  SolverCG (SolverControl        &cn,
            const AdditionalData &data=AdditionalData());

  /**
   * Virtual destructor.
   */
  virtual ~SolverCG ();

  /**
   * Solve the linear system $Ax=b$ for x.
   */
  template <typename MatrixType, typename PreconditionerType>
  void
  solve (const MatrixType         &A,
         VectorType               &x,
         const VectorType         &b,
         const PreconditionerType &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:
  /**
   * Interface for derived class. This function gets the current iteration
   * vector, the residual and the update vector in each step. It can be used
   * for a graphical output of the convergence history.
   */
  virtual void print_vectors(const unsigned int step,
                             const VectorType   &x,
                             const VectorType   &r,
                             const VectorType   &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.
   */
  VectorType *Vr;
  VectorType *Vp;
  VectorType *Vz;

  /**
   * Additional parameters.
   */
  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();
};

/*@}*/

/*------------------------- Implementation ----------------------------*/

#ifndef DOXYGEN

template <typename VectorType>
inline
SolverCG<VectorType>::AdditionalData::
AdditionalData (const bool log_coefficients,
                const bool compute_condition_number,
                const bool compute_all_condition_numbers,
                const bool compute_eigenvalues)
  :
  log_coefficients (log_coefficients),
  compute_condition_number(compute_condition_number),
  compute_all_condition_numbers(compute_all_condition_numbers),
  compute_eigenvalues(compute_eigenvalues)
{}



template <typename VectorType>
inline
SolverCG<VectorType>::AdditionalData::
AdditionalData ()
  :
  log_coefficients (false),
  compute_condition_number(false),
  compute_all_condition_numbers(false),
  compute_eigenvalues(false)
{}



template <typename VectorType>
SolverCG<VectorType>::SolverCG (SolverControl        &cn,
                                VectorMemory<VectorType> &mem,
                                const AdditionalData     &data)
  :
  Solver<VectorType>(cn,mem),
  additional_data(data)
{}



template <typename VectorType>
SolverCG<VectorType>::SolverCG (SolverControl        &cn,
                                const AdditionalData &data)
  :
  Solver<VectorType>(cn),
  additional_data(data)
{}



template <typename VectorType>
SolverCG<VectorType>::~SolverCG ()
{}



template <typename VectorType>
void
SolverCG<VectorType>::cleanup()
{
  this->memory.free(Vr);
  this->memory.free(Vp);
  this->memory.free(Vz);
  deallog.pop();
}



template <typename VectorType>
void
SolverCG<VectorType>::print_vectors(const unsigned int,
{'fraction_non_alphanumeric': 0.10216679786220326, 'fraction_numerical': 0.004391597961635663, 'mean_word_length': 3.0718623481781377, 'pattern_counts': {'":': 0, '<': 139, '<?xml version=': 0, '>': 145, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4645706', 'n_tokens_mistral': 6854, 'n_tokens_neox': 6502, 'n_words': 2422}