trilinos_eigen_solver.h

Go to the documentation of this file.

// LIC// ====================================================================
// LIC// This file forms part of oomph-lib, the object-oriented,
// LIC// multi-physics finite-element library, available
// LIC// at http://www.oomph-lib.org.
// LIC//
// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
// LIC//
// LIC// This library is free software; you can redistribute it and/or
// LIC// modify it under the terms of the GNU Lesser General Public
// LIC// License as published by the Free Software Foundation; either
// LIC// version 2.1 of the License, or (at your option) any later version.
// LIC//
// LIC// This library is distributed in the hope that it will be useful,
// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// LIC// Lesser General Public License for more details.
// LIC//
// LIC// You should have received a copy of the GNU Lesser General Public
// LIC// License along with this library; if not, write to the Free Software
// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
// LIC// 02110-1301  USA.
// LIC//
// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
// LIC//
// LIC//====================================================================
#ifndef OOMPH_TRILINOS_EIGEN_SOLVER_HEADER
#define OOMPH_TRILINOS_EIGEN_SOLVER_HEADER
 
// oomph-lib includes
#include "double_multi_vector.h"
#include "problem.h"
#include "linear_solver.h"
#include "eigen_solver.h"
 
// Trilinos includes
#include "AnasaziOutputManager.hpp"
#include "AnasaziBasicOutputManager.hpp"
#include "AnasaziConfigDefs.hpp"
#include "AnasaziOperator.hpp"
#include "AnasaziMultiVec.hpp"
#include "AnasaziBasicEigenproblem.hpp"
#include "AnasaziBlockKrylovSchurSolMgr.hpp"
 
namespace Anasazi
{
  //========================================================================
  //========================================================================
  template<>
  class MultiVecTraits<double, oomph::DoubleMultiVector>
  {
  public:
    static Teuchos::RCP<oomph::DoubleMultiVector> Clone(
      const oomph::DoubleMultiVector& mv, const int numvecs)
    {
      return Teuchos::rcp(new oomph::DoubleMultiVector(numvecs, mv));
    }
 
    static Teuchos::RCP<oomph::DoubleMultiVector> CloneCopy(
      const oomph::DoubleMultiVector& mv)
    {
      return Teuchos::rcp(new oomph::DoubleMultiVector(mv));
    }
 
    static Teuchos::RCP<oomph::DoubleMultiVector> CloneCopy(
      const oomph::DoubleMultiVector& mv, const std::vector<int>& index)
    {
      return Teuchos::rcp(new oomph::DoubleMultiVector(mv, index));
    }
 
    static Teuchos::RCP<oomph::DoubleMultiVector> CloneCopy(
      const oomph::DoubleMultiVector& mv, const Teuchos::Range1D& index)
    {
      return Teuchos::rcp(new oomph::DoubleMultiVector(mv, index));
    }
 
    static Teuchos::RCP<oomph::DoubleMultiVector> CloneViewNonConst(
      oomph::DoubleMultiVector& mv, const std::vector<int>& index)
    {
      return Teuchos::rcp(new oomph::DoubleMultiVector(mv, index, false));
    }
 
    static Teuchos::RCP<oomph::DoubleMultiVector> CloneViewNonConst(
      oomph::DoubleMultiVector& mv, const Teuchos::Range1D& index)
    {
      return Teuchos::rcp(new oomph::DoubleMultiVector(mv, index, false));
    }
 
    static Teuchos::RCP<const oomph::DoubleMultiVector> CloneView(
      const oomph::DoubleMultiVector& mv, const std::vector<int>& index)
    {
      return Teuchos::rcp(new oomph::DoubleMultiVector(mv, index, false));
    }
 
    static Teuchos::RCP<oomph::DoubleMultiVector> CloneView(
      oomph::DoubleMultiVector& mv, const std::vector<int>& index)
    {
      return Teuchos::rcp(new oomph::DoubleMultiVector(mv, index, false));
    }
 
    static Teuchos::RCP<oomph::DoubleMultiVector> CloneView(
      oomph::DoubleMultiVector& mv, const Teuchos::Range1D& index)
    {
      return Teuchos::rcp(new oomph::DoubleMultiVector(mv, index, false));
    }
 
    static int GetVecLength(const oomph::DoubleMultiVector& mv)
    {
      return static_cast<int>(mv.nrow());
    }
 
    static int GetNumberVecs(const oomph::DoubleMultiVector& mv)
    {
      return static_cast<int>(mv.nvector());
    }
 
    static void MvTimesMatAddMv(
      const double alpha,
      const oomph::DoubleMultiVector& A,
      const Teuchos::SerialDenseMatrix<int, double>& B,
      const double beta,
      oomph::DoubleMultiVector& mv)
    {
      // For safety let's (deep) copy A
      oomph::DoubleMultiVector C(A);
      // First pre-multiply by the scalars
      mv *= beta;
      C *= alpha;
 
      // Find the number of vectors in each multivector
      int mv_n_vector = mv.nvector();
      int A_n_vector = A.nvector();
      // Find the number of rows
      int n_row_local = A.nrow_local();
      // Now need to multiply by the matrix and add the result
      // to the appropriate entry in the multivector
      for (int i = 0; i < n_row_local; ++i)
      {
        for (int j = 0; j < mv_n_vector; j++)
        {
          double res = 0.0;
          for (int k = 0; k < A_n_vector; k++)
          {
            res += C(k, i) * B(k, j);
          }
          mv(j, i) += res;
        }
      }
    }
 
    static void MvAddMv(const double alpha,
                        const oomph::DoubleMultiVector& A,
                        const double beta,
                        const oomph::DoubleMultiVector& B,
                        oomph::DoubleMultiVector& mv)
    {
      const unsigned n_vector = A.nvector();
      const unsigned n_row_local = A.nrow_local();
      for (unsigned v = 0; v < n_vector; v++)
      {
        for (unsigned i = 0; i < n_row_local; i++)
        {
          mv(v, i) = alpha * A(v, i) + beta * B(v, i);
        }
      }
    }
 
    static void MvScale(oomph::DoubleMultiVector& mv, const double alpha)
    {
      mv *= alpha;
    }
 
    static void MvScale(oomph::DoubleMultiVector& mv,
                        const std::vector<double>& alpha)
    {
      const unsigned n_vector = mv.nvector();
      const unsigned n_row_local = mv.nrow_local();
      for (unsigned v = 0; v < n_vector; v++)
      {
        for (unsigned i = 0; i < n_row_local; i++)
        {
          mv(v, i) *= alpha[v];
        }
      }
    }
 
    static void MvTransMv(const double alpha,
                          const oomph::DoubleMultiVector& A,
                          const oomph::DoubleMultiVector& mv,
                          Teuchos::SerialDenseMatrix<int, double>& B)
    {
      const unsigned A_nvec = A.nvector();
      const unsigned A_nrow_local = A.nrow_local();
      const unsigned mv_nvec = mv.nvector();
      // const unsigned mv_nrow_local = mv.nrow_local();
 
      for (unsigned v1 = 0; v1 < A_nvec; v1++)
      {
        for (unsigned v2 = 0; v2 < mv_nvec; v2++)
        {
          B(v1, v2) = 0.0;
          for (unsigned i = 0; i < A_nrow_local; i++)
          {
            B(v1, v2) += A(v1, i) * mv(v2, i);
          }
          B(v1, v2) *= alpha;
        }
      }
 
      // This will need to be communicated if we're in parallel and distributed
#ifdef OOMPH_HAS_MPI
      if (A.distributed() &&
          A.distribution_pt()->communicator_pt()->nproc() > 1)
      {
        const int n_total_val = A_nvec * mv_nvec;
        // Pack entries into a vector
        double b[n_total_val];
        double b2[n_total_val];
        unsigned count = 0;
        for (unsigned v1 = 0; v1 < A_nvec; v1++)
        {
          for (unsigned v2 = 0; v2 < mv_nvec; v2++)
          {
            b[count] = B(v1, v2);
            b2[count] = b[count];
            ++count;
          }
        }
 
 
        MPI_Allreduce(&b,
                      &b2,
                      n_total_val,
                      MPI_DOUBLE,
                      MPI_SUM,
                      A.distribution_pt()->communicator_pt()->mpi_comm());
 
        count = 0;
        for (unsigned v1 = 0; v1 < A_nvec; v1++)
        {
          for (unsigned v2 = 0; v2 < mv_nvec; v2++)
          {
            B(v1, v2) = b2[count];
            ++count;
          }
        }
      }
#endif
    }
 
 
    static void MvDot(const oomph::DoubleMultiVector& mv,
                      const oomph::DoubleMultiVector& A,
                      std::vector<double>& b)
    {
      mv.dot(A, b);
    }
 
 
    static void MvNorm(const oomph::DoubleMultiVector& mv,
                       std::vector<double>& normvec)
    {
      mv.norm(normvec);
    }
 
 
 
 
    static void SetBlock(const oomph::DoubleMultiVector& A,
                         const std::vector<int>& index,
                         oomph::DoubleMultiVector& mv)
    {
      // Check some stuff
      const unsigned n_index = index.size();
      if (n_index == 0)
      {
        return;
      }
      const unsigned n_row_local = mv.nrow_local();
      for (unsigned v = 0; v < n_index; v++)
      {
        for (unsigned i = 0; i < n_row_local; i++)
        {
          mv(index[v], i) = A(v, i);
        }
      }
    }
 
    static void SetBlock(const oomph::DoubleMultiVector& A,
                         const Teuchos::Range1D& index,
                         oomph::DoubleMultiVector& mv)
    {
      // Check some stuff
      const unsigned n_index = index.size();
      if (n_index == 0)
      {
        return;
      }
      const unsigned n_row_local = mv.nrow_local();
      unsigned range_index = index.lbound();
      for (unsigned v = 0; v < n_index; v++)
      {
        for (unsigned i = 0; i < n_row_local; i++)
        {
          mv(range_index, i) = A(v, i);
        }
        ++range_index;
      }
    }
 
    static void Assign(const oomph::DoubleMultiVector& A,
                       oomph::DoubleMultiVector& mv)
    {
      mv = A;
    }
 
    static void MvRandom(oomph::DoubleMultiVector& mv)
    {
      const unsigned n_vector = mv.nvector();
      const unsigned n_row_local = mv.nrow_local();
      for (unsigned v = 0; v < n_vector; v++)
      {
        for (unsigned i = 0; i < n_row_local; i++)
        {
          mv(v, i) = Teuchos::ScalarTraits<double>::random();
        }
      }
    }
 
    static void MvInit(
      oomph::DoubleMultiVector& mv,
      const double alpha = Teuchos::ScalarTraits<double>::zero())
    {
      mv.initialise(alpha);
    }
 
 
 
 
    static void MvPrint(const oomph::DoubleMultiVector& mv, std::ostream& os)
    {
      mv.output(os);
    }
 
  };
 
 
} // namespace Anasazi
 
namespace oomph
{
  //===================================================================
  //==================================================================
  class DoubleMultiVectorOperator
  {
  public:
    DoubleMultiVectorOperator() {}
 
    virtual ~DoubleMultiVectorOperator() {}
 
    virtual void apply(const DoubleMultiVector& x,
                       DoubleMultiVector& y) const = 0;
  };
 
} // namespace oomph
 
namespace Anasazi
{
  //======================================================================
  //======================================================================
  template<>
  class OperatorTraits<double,
                       oomph::DoubleMultiVector,
                       oomph::DoubleMultiVectorOperator>
  {
  public:
    static void Apply(const oomph::DoubleMultiVectorOperator& Op,
                      const oomph::DoubleMultiVector& x,
                      oomph::DoubleMultiVector& y)
    {
      Op.apply(x, y);
    }
 
  }; // End of the specialised traits class
 
} // namespace Anasazi
 
 
namespace oomph
{
  //====================================================================
  //====================================================================
  class ProblemBasedShiftInvertOperator : public DoubleMultiVectorOperator
  {
  private:
    // Pointer to the problem
    Problem* Problem_pt;
 
    // Pointer to the linear solver used
    LinearSolver* Linear_solver_pt;
 
    // Storage for the shift
    double Sigma;
 
    // Storage for the matrices
    CRDoubleMatrix *M_pt, *AsigmaM_pt;
 
  public:
    ProblemBasedShiftInvertOperator(Problem* const& problem_pt,
                                    LinearSolver* const& linear_solver_pt,
                                    const double& sigma = 0.0)
      : Problem_pt(problem_pt), Linear_solver_pt(linear_solver_pt), Sigma(sigma)
    {
      // Before we get into the Arnoldi loop solve the shifted matrix problem
      // Allocated Row compressed matrices for the mass matrix and shifted main
      // matrix
      M_pt = new CRDoubleMatrix(problem_pt->dof_distribution_pt());
      AsigmaM_pt = new CRDoubleMatrix(problem_pt->dof_distribution_pt());
 
      // Assemble the matrices
      problem_pt->get_eigenproblem_matrices(*M_pt, *AsigmaM_pt, Sigma);
 
      // Do not report the time taken
      Linear_solver_pt->disable_doc_time();
    }
 
 
    // Now specify how to apply the operator
    void apply(const DoubleMultiVector& x, DoubleMultiVector& y) const
    {
      const unsigned n_vec = x.nvector();
      const unsigned n_row_local = x.nrow_local();
      if (n_vec > 1)
      {
        Linear_solver_pt->enable_resolve();
      }
      // Solve the first vector
 
      DoubleVector X(x.distribution_pt());
      // Premultiply by mass matrix
      M_pt->multiply(x.doublevector(0), X);
      // For some reason I need to create a new vector and copy here
      // This is odd and not expected, examine carefully
      DoubleVector Y(x.distribution_pt());
      Linear_solver_pt->solve(AsigmaM_pt, X, Y);
      // Need to synchronise
      //#ifdef OOMPH_HAS_MPI
      //   Problem_pt->synchronise_all_dofs();
      //#endif
      for (unsigned i = 0; i < n_row_local; i++)
      {
        y(0, i) = Y[i];
      }
 
      // Now loop over the others and resolve
      for (unsigned v = 1; v < n_vec; ++v)
      {
        M_pt->multiply(x.doublevector(v), X);
        Linear_solver_pt->resolve(X, Y);
        //#ifdef OOMPH_HAS_MPI
        //     Problem_pt->synchronise_all_dofs();
        //#endif
        for (unsigned i = 0; i < n_row_local; i++)
        {
          y(v, i) = Y[i];
        }
      }
    }
  };
 
 
  //=====================================================================
  //=====================================================================
  class ANASAZI : public EigenSolver
  {
  private:
    typedef double ST;
    typedef Teuchos::ScalarTraits<ST> SCT;
    typedef SCT::magnitudeType MT;
    typedef Anasazi::MultiVec<ST> MV;
    typedef Anasazi::Operator<ST> OP;
    typedef Anasazi::MultiVecTraits<ST, MV> MVT;
    typedef Anasazi::OperatorTraits<ST, MV, OP> OPT;
 
    Anasazi::OutputManager<ST>* Output_manager_pt;
 
    LinearSolver* Linear_solver_pt;
 
    LinearSolver* Default_linear_solver_pt;
 
    int Spectrum;
 
    int NArnoldi;
 
    double Sigma;
 
    bool Small;
 
    bool Compute_eigenvectors;
 
 
  public:
    ANASAZI()
      : Linear_solver_pt(0),
        Default_linear_solver_pt(0),
        Spectrum(0),
        NArnoldi(10),
        Sigma(0.0),
        Compute_eigenvectors(true)
    {
      Output_manager_pt = new Anasazi::BasicOutputManager<ST>();
      // Set verbosity level
      int verbosity =
        Anasazi::Warnings + Anasazi::FinalSummary + Anasazi::TimingDetails;
      Output_manager_pt->setVerbosity(verbosity);
 
      // print greeting
      Output_manager_pt->stream(Anasazi::Warnings)
        << Anasazi::Anasazi_Version() << std::endl
        << std::endl;
    }
 
    ANASAZI(const ANASAZI&) : Sigma(0.0) {}
 
    virtual ~ANASAZI() {}
 
    int& narnoldi()
    {
      return NArnoldi;
    }
 
    const int& narnoldi() const
    {
      return NArnoldi;
    }
 
    void enable_compute_eigenvectors()
    {
      Compute_eigenvectors = true;
    }
 
    void disable_compute_eigenvectors()
    {
      Compute_eigenvectors = false;
    }
 
    void solve_eigenproblem(Problem* const& problem_pt,
                            const int& n_eval,
                            Vector<std::complex<double>>& eigenvalue,
                            Vector<DoubleVector>& eigenvector)
    {
      // No access to sigma, so set from sigma real
      Sigma = Sigma_real;
 
      // Initially be dumb here
      Linear_solver_pt = problem_pt->linear_solver_pt();
 
      // Let's make the initial vector
      Teuchos::RCP<DoubleMultiVector> initial = Teuchos::rcp(
        new DoubleMultiVector(1, problem_pt->dof_distribution_pt()));
      Anasazi::MultiVecTraits<double, DoubleMultiVector>::MvRandom(*initial);
 
      // Make the operator
      Teuchos::RCP<DoubleMultiVectorOperator> Op_pt =
        Teuchos::rcp(new ProblemBasedShiftInvertOperator(
          problem_pt, this->linear_solver_pt(), Sigma));
 
      // Create the basic problem
      Teuchos::RCP<Anasazi::BasicEigenproblem<double,
                                              DoubleMultiVector,
                                              DoubleMultiVectorOperator>>
        anasazi_pt = Teuchos::rcp(
          new Anasazi::BasicEigenproblem<double,
                                         DoubleMultiVector,
                                         DoubleMultiVectorOperator>(Op_pt,
                                                                    initial));
 
      // Think I have it?
      anasazi_pt->setHermitian(false);
 
      // set the number of eigenvalues requested
      anasazi_pt->setNEV(n_eval);
 
      // Inform the eigenproblem that you are done passing it information
      if (anasazi_pt->setProblem() != true)
      {
        oomph_info << "Anasazi::BasicEigenproblem::setProblem() had an error."
                   << std::endl
                   << "End Result: TEST FAILED" << std::endl;
      }
 
      // Create the solver manager
      // No need to have ncv specificed, Triliinos has a sensible default
      //  int ncv = 10;
      MT tol = 1.0e-10;
      int verbosity =
        Anasazi::Warnings + Anasazi::FinalSummary + Anasazi::TimingDetails;
 
      Teuchos::ParameterList MyPL;
      MyPL.set("Which", "LM");
      MyPL.set("Block Size", 1);
      //  MyPL.set( "Num Blocks", ncv );
      MyPL.set("Maximum Restarts", 500);
      MyPL.set("Orthogonalization", "DGKS");
      MyPL.set("Verbosity", verbosity);
      MyPL.set("Convergence Tolerance", tol);
      Anasazi::BlockKrylovSchurSolMgr<double,
                                      DoubleMultiVector,
                                      DoubleMultiVectorOperator>
        BKS(anasazi_pt, MyPL);
 
      // Solve the problem to the specified tolerances or length
      Anasazi::ReturnType ret = BKS.solve();
      bool testFailed = false;
      if (ret != Anasazi::Converged)
      {
        testFailed = true;
      }
 
      if (testFailed)
      {
        oomph_info << "Eigensolver not converged\n";
      }
 
      // Get the eigenvalues and eigenvectors from the eigenproblem
      Anasazi::Eigensolution<double, DoubleMultiVector> sol =
        anasazi_pt->getSolution();
      std::vector<Anasazi::Value<double>> evals = sol.Evals;
      Teuchos::RCP<DoubleMultiVector> evecs = sol.Evecs;
 
      eigenvalue.resize(evals.size());
      eigenvector.resize(evals.size());
 
      for (unsigned i = 0; i < evals.size(); i++)
      {
        // Undo shift and invert
        double a = evals[i].realpart;
        double b = evals[i].imagpart;
        double det = a * a + b * b;
        eigenvalue[i] = std::complex<double>(a / det + Sigma, -b / det);
 
        // Now set the eigenvectors, I hope
        eigenvector[i].build(evecs->distribution_pt());
        unsigned nrow_local = evecs->nrow_local();
        // Would be faster with pointers, but I'll sort that out later!
        for (unsigned n = 0; n < nrow_local; n++)
        {
          eigenvector[i][n] = (*evecs)(i, n);
        }
      }
    }
 
    void get_eigenvalues_left_of_shift()
    {
      Small = true;
    }
 
    void get_eigenvalues_right_of_shift()
    {
      Small = false;
    }
 
    void track_eigenvalue_real_part()
    {
      Spectrum = 1;
    }
 
    void track_eigenvalue_imaginary_part()
    {
      Spectrum = -1;
    }
 
    void track_eigenvalue_magnitude()
    {
      Spectrum = 0;
    }
 
    LinearSolver*& linear_solver_pt()
    {
      return Linear_solver_pt;
    }
 
    LinearSolver* const& linear_solver_pt() const
    {
      return Linear_solver_pt;
    }
  };
 
} // namespace oomph
 
 
#endif