eigen.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2000 - 2017 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_eigen_h
#define dealii_eigen_h


#include <deal.II/base/config.h>
#include <deal.II/lac/linear_operator.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/solver_control.h>
#include <deal.II/lac/solver_gmres.h>
#include <deal.II/lac/solver.h>
#include <deal.II/lac/solver_minres.h>
#include <deal.II/lac/vector_memory.h>

#include <cmath>

DEAL_II_NAMESPACE_OPEN



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

  struct AdditionalData
  {
    double shift;
    AdditionalData (const double shift = 0.)
      :
      shift(shift)
    {}

  };

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

  virtual ~EigenPower ();

  template <typename MatrixType>
  void
  solve (double           &value,
         const MatrixType &A,
         VectorType       &x);

protected:
  AdditionalData additional_data;
};

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

  struct AdditionalData
  {
    double relaxation;

    unsigned int start_adaption;
    bool use_residual;
    AdditionalData (double relaxation = 1.,
                    unsigned int start_adaption = 6,
                    bool use_residual = true)
      :
      relaxation(relaxation),
      start_adaption(start_adaption),
      use_residual(use_residual)
    {}

  };

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


  virtual ~EigenInverse ();

  template <typename MatrixType>
  void
  solve (double           &value,
         const MatrixType &A,
         VectorType       &x);

protected:
  AdditionalData additional_data;
};

//---------------------------------------------------------------------------


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



template <class VectorType>
EigenPower<VectorType>::~EigenPower ()
{}



template <class VectorType>
template <typename MatrixType>
void
EigenPower<VectorType>::solve (double           &value,
                               const MatrixType &A,
                               VectorType       &x)
{
  SolverControl::State conv=SolverControl::iterate;

  LogStream::Prefix prefix("Power method");

  typename VectorMemory<VectorType>::Pointer Vy (this->memory);
  VectorType &y = *Vy;
  y.reinit (x);
  typename VectorMemory<VectorType>::Pointer Vr (this->memory);
  VectorType &r = *Vr;
  r.reinit (x);

  double length = x.l2_norm ();
  double old_length = 0.;
  x *= 1./length;

  A.vmult (y,x);

  // Main loop
  int iter=0;
  for (; conv==SolverControl::iterate; iter++)
    {
      y.add(additional_data.shift, x);

      // Compute absolute value of eigenvalue
      old_length = length;
      length = y.l2_norm ();

      // do a little trick to compute the sign
      // with not too much effect of round-off errors.
      double entry = 0.;
      size_type i = 0;
      double thresh = length/x.size();
      do
        {
          Assert (i<x.size(), ExcInternalError());
          entry = y (i++);
        }
      while (std::fabs(entry) < thresh);

      --i;

      // Compute unshifted eigenvalue
      value = (entry * x (i) < 0.) ? -length : length;
      value -= additional_data.shift;

      // Update normalized eigenvector
      x.equ (1/length, y);

      // Compute residual
      A.vmult (y,x);

      // Check the change of the eigenvalue
      // Brrr, this is not really a good criterion
      conv = this->iteration_status (iter, std::fabs(1./length-1./old_length), x);
    }

  // in case of failure: throw exception
  AssertThrow(conv == SolverControl::success, SolverControl::NoConvergence (iter,
              std::fabs(1./length-1./old_length)));

  // otherwise exit as normal
}

//---------------------------------------------------------------------------

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



template <class VectorType>
EigenInverse<VectorType>::~EigenInverse ()
{}



template <class VectorType>
template <typename MatrixType>
void
EigenInverse<VectorType>::solve (double           &value,
                                 const MatrixType &A,
                                 VectorType       &x)
{
  LogStream::Prefix prefix("Wielandt");

  SolverControl::State conv=SolverControl::iterate;

  // Prepare matrix for solver
  auto A_op = linear_operator(A);
  double current_shift = -value;
  auto A_s = A_op + current_shift * identity_operator(A_op);

  // Define solver
  ReductionControl inner_control (5000, 1.e-16, 1.e-5, false, false);
  PreconditionIdentity prec;
  SolverGMRES<VectorType>
  solver(inner_control, this->memory);

  // Next step for recomputing the shift
  unsigned int goal = additional_data.start_adaption;

  // Auxiliary vector
  typename VectorMemory<VectorType>::Pointer Vy (this->memory);
  VectorType &y = *Vy;
  y.reinit (x);
  typename VectorMemory<VectorType>::Pointer Vr (this->memory);
  VectorType &r = *Vr;
  r.reinit (x);

  double length = x.l2_norm ();
  double old_value = value;

  x *= 1./length;

  // Main loop
  double res = -std::numeric_limits<double>::max();
  size_type iter=0;
  for (; conv==SolverControl::iterate; iter++)
    {
      solver.solve (A_s, y, x, prec);

      // Compute absolute value of eigenvalue
      length = y.l2_norm ();

      // do a little trick to compute the sign
      // with not too much effect of round-off errors.
      double entry = 0.;
      size_type i = 0;
      double thresh = length/x.size();
      do
        {
          Assert (i<x.size(), ExcInternalError());
          entry = y (i++);
        }
      while (std::fabs(entry) < thresh);

      --i;

      // Compute unshifted eigenvalue
      value = (entry * x(i) < 0. ? -1. : 1.) / length - current_shift;

      if (iter==goal)
        {
          const auto &relaxation = additional_data.relaxation;
          const double new_shift =
            relaxation * (-value) + (1. - relaxation) * current_shift;

          A_s = A_op + new_shift * identity_operator(A_op);
          current_shift = new_shift;

          ++goal;
        }

      // Update normalized eigenvector
      x.equ (1./length, y);
      // Compute residual
      if (additional_data.use_residual)
        {
          y.equ (value, x);
          A.vmult(r,x);
          r.sadd(-1., value, x);
          res = r.l2_norm();
          // Check the residual
          conv = this->iteration_status (iter, res, x);
        }
      else
        {
          res = std::fabs(1./value-1./old_value);
          conv = this->iteration_status (iter, res, x);
        }
      old_value = value;
    }

  // in case of failure: throw
  // exception
  AssertThrow (conv == SolverControl::success,
               SolverControl::NoConvergence (iter,
                                             res));
  // otherwise exit as normal
}

DEAL_II_NAMESPACE_CLOSE

#endif