eigen_solver.cc

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//    Version 1.0; svn revision $LastChangedRevision$
//LIC//
//LIC// $LastChangedDate$
//LIC// 
//LIC// Copyright (C) 2006-2016 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//====================================================================
//Non-inline functions for the eigensolvers

#ifdef OOMPH_HAS_MPI
#include "mpi.h"
#endif

//Include cfortran.h and the header for the FORTRAN ARPACK routines
#include "cfortran.h"
#include "arpack.h"
#include "lapack_qz.h"

//Oomph-lib headers
#include "eigen_solver.h"
#include "linear_solver.h"
#include "problem.h"


namespace oomph
{
//===============================================================
/// Constructor, set default values and set the initial 
/// linear solver to be superlu
//===============================================================
 ARPACK::ARPACK() : EigenSolver(), Spectrum(1), NArnoldi(30),
                    Small(true), Compute_eigenvectors(true)
 {Default_linear_solver_pt = Linear_solver_pt = new SuperLUSolver;}

//===============================================================
/// Destructor, delete the default linear solver
//===============================================================
 ARPACK::~ARPACK() {delete Default_linear_solver_pt;}


//==========================================================================
/// Use ARPACK to solve an eigen problem that is assembled by elements in
/// a mesh in a Problem object.
//==========================================================================
void ARPACK::solve_eigenproblem(Problem* const &problem_pt,
                                const int &n_eval,
                                Vector<std::complex<double> > &eigenvalue,
                                Vector<DoubleVector > &eigenvector)
{
 bool Verbose = false;
 
 //Control parameters
 int ido=0; //Reverse communication flag
 char bmat[2]; //Specifies the type of matrix on the RHS
 //Must be 'I' for standard eigenvalue problem
 //        'G' for generalised eigenvalue problem
 int n; //Dimension of the problem
 char which[3]; //Set which eigenvalues are required.
 int nev; //Number of eigenvalues desired
 double tol=0.0; //Stopping criteria
 int ncv; //Number of columns of the matrix V (Dimension of Arnoldi subspace)
 int iparam[11]={0,0,0,0,0,0,0,0,0,0,0}; //Control parameters
 //Pointers to starting locations in the workspace vectors
 int ipntr[14]={0,0,0,0,0,0,0,0,0,0,0,0,0,0}; 
 int info=0; //Setting to zero gives random initial guess for vector to start
 //Arnoldi iteration
 
 //Set up the sizes of the matrix
 n     = problem_pt->ndof(); // Total size of matrix 
 nev   = n_eval; //Number of desired eigenvalues 
 ncv   = NArnoldi < n ? NArnoldi : n; //Number of Arnoldi vectors allowed 
                                      //Maximum possible value is max
                                      //dimension of matrix

 //If we don't have enough Arnoldi vectors to compute the desired number
 //of eigenvalues then complain
 if(nev > ncv) 
  {
   std::ostringstream warning_stream;
   warning_stream << "Number of requested eigenvalues " << nev << "\n"
                  << "is greater than the number of Arnoldi vectors:"
                  << ncv << "\n";
   //Increase the number of Arnoldi vectors
   ncv = nev + 10;
   if(ncv > n) {ncv = n;}

   warning_stream << "Increasing number of Arnoldi vectors to " << ncv
                  << "\n but you may want to increase further using\n"
                  << "ARPACK::narnoldi()\n"
                  << "which will also get rid of this warning.\n";

   OomphLibWarning(warning_stream.str(),
                   OOMPH_CURRENT_FUNCTION,
                   OOMPH_EXCEPTION_LOCATION);
  }

 //Allocate some workspace
 int lworkl  = 3*ncv*ncv + 6*ncv; 
 
 //Array that holds the final set of Arnoldi basis vectors
 double **v = new double*;
 *v = new double[ncv*n];
 //Leading dimension of v (n)
 int ldv = n;
 
 //Residual vector
 double *resid = new double[n];
 //Workspace vector
 double *workd = new double[3*n];
 //Workspace vector
 double *workl = new double[lworkl];
 
 bmat[0]  = 'G';
 //If we require eigenvalues to the left of the shift
 if(Small) {which[0] = 'S';} 
 //Otherwise we require eigenvalues to the right of the shift
 else {which[0] = 'L';}

 //Which part of the eigenvalue interests us
 switch(Spectrum)
  {
   //Imaginary part
  case -1:
   which[1] = 'I';
   break;
   
   //Magnitude
  case 0:
   which[1] = 'M';
   break;

   //Real part
  case 1:
   which[1] = 'R';
   break;

  default:
   std::ostringstream error_stream;
   error_stream 
    << "Spectrum is set to an invalid value. It must be 0, 1 or -1\n"
    << "Not " << Spectrum << std::endl;

   throw OomphLibError(error_stream.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }

//     %---------------------------------------------------%
//     | This program uses exact shifts with respect to    |
//     | the current Hessenberg matrix (IPARAM(1) = 1).    |
//     | IPARAM(3) specifies the maximum number of Arnoldi |
//     | iterations allowed.  Mode 3 of DNAUPD is used     |
//     | (IPARAM(7) = 3). All these options can be changed |
//     | by the user. For details see the documentation in |
//     | DNAUPD.                                           |
//     %---------------------------------------------------%

 int ishifts = 1;
 int maxitr = 300;
 int mode   = 3; //M is symetric and semi-definite

 iparam[0] = ishifts; //Exact shifts wrt Hessenberg matrix H
 iparam[2] = maxitr;  //Maximum number of allowed iterations
 iparam[3] = 1;       //Bloacksize to be used in the recurrence
 iparam[6] = mode;    //Mode is shift and invert in real arithmetic

 //Real and imaginary parts of the shift
 double sigmar = Sigma_real, sigmai = 0.0;

 // TEMPORARY
 // only use non distributed matrice and vectors
 LinearAlgebraDistribution dist(problem_pt->communicator_pt(),n,false);
 this->build_distribution(dist);

 //Before we get into the Arnoldi loop solve the shifted matrix problem
 //Allocated Row compressed matrices for the mass matrix and shifted main 
 //matrix
 CRDoubleMatrix M(this->distribution_pt()), AsigmaM(this->distribution_pt());

 //Assemble the matrices
 problem_pt->get_eigenproblem_matrices(M,AsigmaM,sigmar);

 //Allocate storage for the vectors to be used in matrix vector products
 DoubleVector rhs(this->distribution_pt(),0.0);
 DoubleVector x(this->distribution_pt(),0.0);
 

 bool LOOP_FLAG=true;
 bool First=true;

 //Enable resolves for the linear solver
 Linear_solver_pt->enable_resolve();
 
 //Do not report the time taken
 Linear_solver_pt->disable_doc_time();

 do
  {

//        %---------------------------------------------%
//        | Repeatedly call the routine DNAUPD and take |
//        | actions indicated by parameter IDO until    |
//        | either convergence is indicated or maxitr   |
//        | has been exceeded.                          |
//        %---------------------------------------------%
//

   DNAUPD(ido, bmat, n, which, nev, tol, resid, 
          ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, info);

   if(ido == -1)
    {
     
//           %-------------------------------------------%
//           | Perform matrix vector multiplication      |
//           |                y <--- OP*x                |
//           | The user should supply his/her own        |
//           | matrix vector multiplication routine here |
//           | that takes workd(ipntr(1)) as the input   |
//           | vector, and return the matrix vector      |
//           | product to workd(ipntr(2)).               | 
//           %-------------------------------------------%
//

     //Firstly multiply by mx
     for(int i=0;i<n;i++) {x[i] = workd[ipntr[0]-1+i];}
     M.multiply(x,rhs);

     //Now solve the system
     DoubleVector temp(rhs);
     if(First) {Linear_solver_pt->solve(&AsigmaM,temp,rhs); First=false;}
     else {Linear_solver_pt->resolve(temp,rhs);}
     temp.clear();
     //AsigmaM.lubksub(rhs);
     //Put the solution into the workspace
     for(int i=0;i<n;i++)
      {
       workd[ipntr[1]-1+i] = rhs[i];
      }
     
     continue;
    }
   else if(ido == 1)
    {
     //Already done multiplication by Mx
     //Need to load the rhs vector
     for(int i=0;i<n;i++)
      {
       rhs[i] = workd[ipntr[2]-1+i];
      }
     //Now solve the system
     //AsigmaM.lubksub(rhs);
     DoubleVector temp(rhs);
     if(First) {Linear_solver_pt->solve(&AsigmaM,temp,rhs); First=false;}
     else {Linear_solver_pt->resolve(temp,rhs);}
     //Put the solution into the workspace
     for(int i=0;i<n;i++)
      {
       workd[ipntr[1]-1+i] = rhs[i];
      }
     continue;
    }
   else if(ido == 2)
    {
     //Need to multiply by mass matrix
     //Vector<double> x(n);
     for(int i=0;i<n;i++) {x[i] = workd[ipntr[0]-1+i];}
     M.multiply(x,rhs);
     
     for(int i=0;i<n;i++)
      {
       workd[ipntr[1]-1+i] = rhs[i];
      }
     continue;
    }

   if(info < 0)
    {
     oomph_info << "ERROR" << std::endl;
     oomph_info << "Error with _naupd, info = '" <<  info << std::endl;
     if(info==-7) {oomph_info << "Not enough memory " << std::endl;}
    }

//        %-------------------------------------------%
//        | No fatal errors occurred.                 |
//        | Post-Process using DNEUPD.                |
//        |                                           |
//        | Computed eigenvalues may be extracted.    |
//        |                                           |
//        | Eigenvectors may also be computed now if  |
//        | desired.  (indicated by rvec = .true.)    |
//        %-------------------------------------------%
//
   LOOP_FLAG = false;
  } while(LOOP_FLAG);

 //Report any problems
 if (info == 1) 
  {
   oomph_info << "Maximum number of iterations reached." << std::endl;
  }
 else if(info == 3)
  {
   oomph_info
    << "No shifts could be applied during implicit Arnoldi "
    << "update, try increasing NCV. "
    << std::endl;
  }
 
 //Disable resolves for the linear solver
 Linear_solver_pt->disable_resolve();
 
 //Compute Ritz or Schur vectors, if desired
 int rvec = Compute_eigenvectors;
 //Specify the form of the basis computed
 char howmany[2];
 howmany[0] = 'A';
 //Find out the number of converged eigenvalues 
 int nconv  = iparam[4];

 //Note that there is an error (feature) in ARPACK that
 //means in certain cases, if we request eigenvectors,
 //consequent Schur factorization of the matrix spanned
 //by the Arnoldi vectors leads to more converged eigenvalues
 //than reported by DNAPUD. This is a pain because it
 //causes a segmentation fault and in every case that I've found
 //the eigenvalues/vectors are not those desired.
 //At the moment, I'm just going to leave it, but I note 
 //the problem here to remind myself about it.

 //Array to select which Ritz vectors are computed
 int select[ncv];
 //Array that holds the real part of the eigenvalues
 double *eval_real = new double[nconv+1];
 //Array that holds the imaginary part of the eigenvalues
 double *eval_imag = new double[nconv+1];
 
 //Workspace
 double *workev = new double[3*ncv];
 //Array to hold the eigenvectors
 double **z = new double*;
 *z = new double[(nconv+1)*n];
 
 //Error flag
 int ierr;
 
 DNEUPD ( rvec, howmany, select, eval_real, eval_imag, z, ldv, 
          sigmar, sigmai, workev, bmat, n, which, nev, tol, 
          resid, ncv, v, ldv, iparam, ipntr, workd, workl,
          lworkl, ierr );
//
//        %-----------------------------------------------%
//        | The real part of the eigenvalue is returned   |
//        | in the first column of the two dimensional    |
//        | array D, and the imaginary part is returned   |
//        | in the second column of D.  The corresponding |
//        | eigenvectors are returned in the first NEV    |
//        | columns of the two dimensional array V if     |
//        | requested.  Otherwise, an orthogonal basis    |
//        | for the invariant subspace corresponding to   |
//        | the eigenvalues in D is returned in V.        |
//        %-----------------------------------------------%

 if(ierr != 0)
  {
   oomph_info << "Error with _neupd, info = " << ierr << std::endl;
   }
 else
 //Print the output anyway
 {
  //Resize the eigenvalue output vector
  eigenvalue.resize(nconv);
  for(int j=0;j<nconv;j++)
   {
    //Turn the output from ARPACK into a complex number
    std::complex<double> eigen(eval_real[j],eval_imag[j]);
    //Add the eigenvalues to the output vector
    eigenvalue[j] = eigen;
   }

  if(Compute_eigenvectors)
   {
    //Load the eigenvectors
    eigenvector.resize(nconv);
    for(int j=0;j<nconv;j++)
     {
      eigenvector[j].build(this->distribution_pt(),0.0);
      for(int i=0;i<n;i++)
       {
        eigenvector[j][i] = z[0][j*n+i];
       }
     }
   }
  else
   {
    eigenvector.resize(0);
   }

   //Report the information
   if(Verbose)
    {
     oomph_info << " ARPACK " << std::endl;
     oomph_info << " ====== " << std::endl;
     oomph_info << " Size of the matrix is " <<  n << std::endl;
     oomph_info << " The number of Ritz values requested is " <<  nev 
                << std::endl;
     oomph_info << " The number of Arnoldi vectors generated (NCV) is " <<  ncv
                << std::endl;
     oomph_info << " What portion of the spectrum: " << which << std::endl;
     oomph_info << " The number of converged Ritz values is " 
                << nconv << std::endl;
     oomph_info 
      << " The number of Implicit Arnoldi update iterations taken is "
      << iparam[2] << std::endl;
     oomph_info << " The number of OP*x is " << iparam[8] << std::endl;
     oomph_info << " The convergence criterion is " <<  tol << std::endl;
    }
 }

 //Clean up the allocated memory
 delete[] *z; *z=0;
 delete z; z=0;
 delete[] workev; workev=0;
 delete[] eval_imag; eval_imag=0;
 delete[] eval_real; eval_real=0;

 //Clean up the memory
 delete[] workl; workl = 0;
 delete[] workd; workd = 0;
 delete[] resid; resid = 0;
 delete[] *v; *v=0;
 delete v; v=0;
}


//==========================================================================
/// Use LAPACK to solve an eigen problem that is assembled by elements in
/// a mesh in a Problem object.
//==========================================================================
 void LAPACK_QZ::solve_eigenproblem(Problem* const &problem_pt,
                                    const int &n_eval,
                                    Vector<std::complex<double> > &eigenvalue,
                                    Vector<DoubleVector> &eigenvector)
{
 //Some character identifiers for use in the LAPACK routine 
 //Do not calculate the left eigenvectors
 char no_eigvecs[2] = "N"; 
 //Do caculate the eigenvectors
 char eigvecs[2] = "V"; 

 //Get the dimension of the matrix
 int n = problem_pt->ndof(); // Total size of matrix 

 //Initialise a padding integer
 int padding=0;
 //If the dimension of the matrix is even, then pad the arrays to
 //make the size odd. This somehow sorts out a strange run-time behaviour
 //identified by Rich Hewitt.
 if(n%2==0) {padding=1;}

 //Get the real and imaginary parts of the shift
 double sigmar = Sigma_real;// sigmai = 0.0;

 //Actual size of matrix that will be allocated
 int padded_n = n + padding;

 //Storage for the matrices in the packed form required by the LAPACK
 //routine
 double *M = new double[padded_n*padded_n];
 double *A = new double[padded_n*padded_n];

 // TEMPORARY
 // only use non-distributed matrices and vectors
 LinearAlgebraDistribution dist(problem_pt->communicator_pt(),n,false);
 this->build_distribution(dist);

 //Enclose in a separate scope so that memory is cleaned after assembly
 {
  //Allocated Row compressed matrices for the mass matrix and shifted main 
  //matrix
  CRDoubleMatrix temp_M(this->distribution_pt()),
   temp_AsigmaM(this->distribution_pt());

  //Assemble the matrices
  problem_pt->get_eigenproblem_matrices(temp_M,temp_AsigmaM,sigmar);
  
  //Now convert these matrices into the appropriate packed form
  unsigned index=0;
  for(int i=0;i<n;++i)
   {
    for(int j=0;j<n;++j)
     {
      M[index] = temp_M(j,i);
      A[index] = temp_AsigmaM(j,i);
      ++index;
     }
    //If necessary pad the columns with zeros
    if(padding)
     {
      M[index] = 0.0;
      A[index] = 0.0;
      ++index;
     }
   }
  //No need to pad the final row because it is never accessed by the
  //routine.
 }
  
 // Temporary eigenvalue storage 
 double* alpha_r = new double[n]; 
 double* alpha_i = new double[n]; 
 double* beta = new double[n]; 
 // Temporary eigenvector storage 
 double* vec_left = new double[1]; 
 double* vec_right = new double[n*n]; 

 // Workspace for the LAPACK routine 
 std::vector<double> work(1,0.0);
 //Info flag for the LAPACK routine
 int info = 0;
 
 // Call FORTRAN LAPACK to get the required workspace
 // Note the use of the padding flag 
 LAPACK_DGGEV( no_eigvecs, eigvecs, n, &A[0], padded_n, &M[0], padded_n, 
              alpha_r, alpha_i, beta, vec_left, 1, vec_right, n, 
               &work[0], -1, info ); 
 
 //Get the amount of requires workspace
 int required_workspace = ( int ) work[0]; 
 //If we need it
 work.resize( required_workspace ); 

// call FORTRAN LAPACK again with the optimum workspace 
 LAPACK_DGGEV( no_eigvecs, eigvecs, n, &A[0], padded_n, &M[0], padded_n, 
               alpha_r, alpha_i, beta, vec_left, 1, vec_right, n, &work[0], 
               required_workspace, info ); 
 
 //Now resize storage for the eigenvalues and eigenvectors
 //We get them all!
 eigenvalue.resize(n);
 eigenvector.resize(n);

 //Now convert the output into our format
 for (int i=0; i<n; ++i ) 
  { 
   //N.B. This is naive and dangerous according to the documentation
   //beta could be very small giving over or under flow
   //Remember to fix the shift back again
   eigenvalue[i] = 
    std::complex<double>(sigmar + alpha_r[i]/beta[i],alpha_i[i]/beta[i]);

   //Resize the eigenvector  storage
   eigenvector[i].build(this->distribution_pt(),0.0);
   //Load up the eigenvector (assume that it's real)
   for(int k = 0; k < n; ++k ) 
    { 
     eigenvector[i][k] = vec_right[i*n + k ];
    }
  }

 // Delete all allocated storage
 delete[] vec_right; 
 delete[] vec_left; 
 delete[] beta; 
 delete[] alpha_r; 
 delete[] alpha_i;
 delete[] A;
 delete[] M;
}


//==========================================================================
/// Use LAPACK to solve a complex eigen problem specified by the given 
/// matrices.
//==========================================================================
void LAPACK_QZ::find_eigenvalues(const ComplexMatrixBase &A,
                                 const ComplexMatrixBase &M,
                                 Vector<std::complex<double> > &eigenvalue,
                                 Vector<Vector<std::complex<double> > >
                                 &eigenvector)
{
 //Some character identifiers for use in the LAPACK routine 
 //Do not calculate the left eigenvectors
 char no_eigvecs[2] = "N"; 
 //Do caculate the eigenvectors
 char eigvecs[2] = "V"; 

 //Get the dimension of the matrix
 int n = A.nrow(); // Total size of matrix 

 //Storage for the matrices in the packed form required by the LAPACK
 //routine
 double *M_linear = new double[2*n*n];
 double *A_linear = new double[2*n*n];

 //Now convert the matrices into the appropriate packed form
 unsigned index=0;
 for(int i=0;i<n;++i)
  {
   for(int j=0;j<n;++j)
    {
     M_linear[index] = M(j,i).real();
     A_linear[index] = A(j,i).real();
     ++index;
     M_linear[index] = M(j,i).imag();
     A_linear[index] = A(j,i).imag();
     ++index;
    }
  }
  
 // Temporary eigenvalue storage 
 double* alpha = new double[2*n]; 
 double* beta  = new double[2*n]; 
 // Temporary eigenvector storage 
 double* vec_left  = new double[2]; 
 double* vec_right = new double[2*n*n]; 

 // Workspace for the LAPACK routine 
 std::vector<double> work(2,0.0);
 std::vector<double> rwork(8*n,0.0);

 //Info flag for the LAPACK routine
 int info = 0;
 
 // Call FORTRAN LAPACK to get the required workspace
 // Note the use of the padding flag 
 LAPACK_ZGGEV( no_eigvecs, eigvecs, n, &A_linear[0], n, &M_linear[0], n, 
              alpha, beta, vec_left, 1, vec_right, n, 
              &work[0], -1, &rwork[0], info ); 
 
 //Get the amount of requires workspace
 int required_workspace = ( int ) work[0]; 
 //If we need it
 work.resize( 2*required_workspace ); 

// call FORTRAN LAPACK again with the optimum workspace 
 LAPACK_ZGGEV( no_eigvecs, eigvecs, n, &A_linear[0], n, &M_linear[0], n, 
               alpha, beta, vec_left, 1, vec_right, n, &work[0], 
               required_workspace, &rwork[0], info ); 
 
 //Now resize storage for the eigenvalues and eigenvectors
 //We get them all!
 eigenvalue.resize(n);
 eigenvector.resize(n);

 //Now convert the output into our format
 for (int i=0; i<n; ++i ) 
  { 
   //We have temporary complex numbers
   std::complex<double> num(alpha[2*i],alpha[2*i+1]);
   std::complex<double> den(beta[2*i],beta[2*i+1]);

   //N.B. This is naive and dangerous according to the documentation
   //beta could be very small giving over or under flow
   eigenvalue[i] = num/den;

   //Resize the eigenvector storage
   eigenvector[i].resize(n);
   //Load up the eigenvector (assume that it's real)
   for(int k = 0; k < n; ++k ) 
    { 
     eigenvector[i][k] = std::complex<double>(vec_right[2*i*n + 2*k],
                                              vec_right[2*i*n + 2*k+1]);
    }
  }

 // Delete all allocated storage
 delete[] vec_right; 
 delete[] vec_left; 
 delete[] beta; 
 delete[] alpha; 
 delete[] A_linear;
 delete[] M_linear;
}
}