solver_cg.h

// ---------------------------------------------------------------------
//
// 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;



template <typename VectorType = Vector<double> >
class SolverCG : public Solver<VectorType>
{
public:
  typedef types::global_dof_index size_type;

  struct AdditionalData
  {
    bool log_coefficients;

    bool compute_condition_number;

    bool compute_all_condition_numbers;

    bool compute_eigenvalues;

    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;

    AdditionalData();
  };

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

  SolverCG (SolverControl        &cn,
            const AdditionalData &data=AdditionalData());

  virtual ~SolverCG ();

  template <typename MatrixType, typename PreconditionerType>
  void
  solve (const MatrixType         &A,
         VectorType               &x,
         const VectorType         &b,
         const PreconditionerType &precondition);

  boost::signals2::connection
  connect_coefficients_slot(
    const std_cxx11::function<void (double,double)> &slot);

  boost::signals2::connection
  connect_condition_number_slot(const std_cxx11::function<void (double)> &slot,
                                const bool every_iteration=false);

  boost::signals2::connection
  connect_eigenvalues_slot(
    const std_cxx11::function<void (const std::vector<double> &)> &slot,
    const bool every_iteration=false);

protected:
  virtual double criterion();

  virtual void print_vectors(const unsigned int step,
                             const VectorType   &x,
                             const VectorType   &r,
                             const VectorType   &d) const;

  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);

  VectorType *Vr;
  VectorType *Vp;
  VectorType *Vz;

  double res2;

  AdditionalData additional_data;

  boost::signals2::signal<void (double,double)> coefficients_signal;

  boost::signals2::signal<void (double)> condition_number_signal;

  boost::signals2::signal<void (double)> all_condition_numbers_signal;

  boost::signals2::signal<void (const std::vector<double> &)> eigenvalues_signal;

  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,
                                    const VectorType &,
                                    const VectorType &,
                                    const VectorType &) const
{}



template <typename VectorType>
inline void
SolverCG<VectorType>::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 <typename VectorType>
template <typename MatrixType, typename PreconditionerType>
void
SolverCG<VectorType>::solve (const MatrixType         &A,
                             VectorType               &x,
                             const VectorType         &b,
                             const PreconditionerType &precondition)
{
  SolverControl::State conv=SolverControl::iterate;

  deallog.push("cg");

  // Memory allocation
  Vr = this->memory.alloc();
  Vz = this->memory.alloc();
  Vp = this->memory.alloc();
  // Should we build the matrix for
  // eigenvalue computations?
  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
  // computations
  std::vector<double> diagonal;
  std::vector<double> offdiagonal;

  int  it=0;
  double res = -std::numeric_limits<double>::max();

  try
    {
      // define some aliases for simpler access
      VectorType &g = *Vr;
      VectorType &d = *Vz;
      VectorType &h = *Vp;
      // resize the vectors, but do not set
      // the values since they'd be overwritten
      // soon anyway.
      g.reinit(x, true);
      d.reinit(x, true);
      h.reinit(x, true);

      double gh,alpha,beta;

      // compute residual. if vector is
      // zero, then short-circuit the
      // full computation
      if (!x.all_zero())
        {
          A.vmult(g,x);
          g.add(-1.,b);
        }
      else
        g.equ(-1.,b);
      res = g.l2_norm();

      conv = this->iteration_status(0, res, x);
      if (conv != SolverControl::iterate)
        {
          cleanup();
          return;
        }

      if (types_are_equal<PreconditionerType,PreconditionIdentity>::value == false)
        {
          precondition.vmult(h,g);

          d.equ(-1.,h);

          gh = g*h;
        }
      else
        {
          d.equ(-1.,g);
          gh = res*res;
        }

      while (conv == SolverControl::iterate)
        {
          it++;
          A.vmult(h,d);

          alpha = d*h;
          Assert(alpha != 0., ExcDivideByZero());
          alpha = gh/alpha;

          x.add(alpha,d);
          res = std::sqrt(g.add_and_dot(alpha, h, g));

          print_vectors(it, x, g, d);

          conv = this->iteration_status(it, res, x);
          if (conv != SolverControl::iterate)
            break;

          if (types_are_equal<PreconditionerType,PreconditionIdentity>::value
              == false)
            {
              precondition.vmult(h,g);

              beta = gh;
              Assert(beta != 0., ExcDivideByZero());
              gh   = g*h;
              beta = gh/beta;
              d.sadd(beta,-1.,h);
            }
          else
            {
              beta = gh;
              gh = res*res;
              beta = gh/beta;
              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
          // containing the diagonal
          // and the off diagonal of
          // the projected matrix.
          if (do_eigenvalues)
            {
              diagonal.push_back(1./alpha + eigen_beta_alpha);
              eigen_beta_alpha = beta/alpha;
              offdiagonal.push_back(std::sqrt(beta)/alpha);
            }
          compute_eigs_and_cond(diagonal,offdiagonal,all_eigenvalues_signal,
                                all_condition_numbers_signal,false,
                                additional_data.compute_all_condition_numbers);
        }
    }
  catch (...)
    {
      cleanup();
      throw;
    }
  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();
  // in case of failure: throw exception
  if (conv != SolverControl::success)
    AssertThrow(false, SolverControl::NoConvergence (it, res));
  // otherwise exit as normal
}



template<typename VectorType>
boost::signals2::connection
SolverCG<VectorType>::connect_coefficients_slot
(const std_cxx11::function<void(double,double)> &slot)
{
  return coefficients_signal.connect(slot);
}



template<typename VectorType>
boost::signals2::connection
SolverCG<VectorType>::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<typename VectorType>
boost::signals2::connection
SolverCG<VectorType>::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

#endif