inputs/dev/EigenSparse.tex

%%%%%%%%%%%%%%%%%%%
% XLiFE++ is an extended library of finite elements written in C++
%     Copyright (C) 2014  Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%

\section {{\libtitle eigenSparse} sub-library}

The {\lib eigenSparse}  provides tools that are useful for solving a wide variety
of eigenvalue problems. The solvers currently implemented within {\lib eigenSparse} compute a
partial eigen-decomposition for the generalized eigenvalue problem \\
\[ Ax = Bx \lambda \]
We assume that the matrices A and B are large, possibly sparse, and that only their application to a block of vectors is required. \\
This section outlines the eigenSparse. Three subsections describe the eigenSparse operator/multivector interface, the eigensolver, and the various implementations provided by the eigenSparse. \\
The lib {\lib eigenSparse} is highly independent of other packages of \xlifepp except for {\lib eigenCore} which provides essential functions in performing the dense arithmetic for Rayleigh-Ritz methods.

\subsection {The Operator/Multivector Interface}

{\lib eigenSparse} utilizes traits classes to define interfaces for
the scalar field, multivectors, and matrix operators. This allows generic programming techniques to be used when developing numerical algorithms in the {\lib eigenSparse}. {\lib eigenSparse}’s eigensolver is comprised of abstract
numerical interfaces that are all implemented using templates and the functionality
of the template arguments is provided through their corresponding trait classes.
Most classes in {\lib eigenSparse} accept three template parameters:
\begin{itemize}
\item a scalar type, describing the field over which the vectors and operators are defined;
\item a multivector type over the given scalar field, providing a data structure that
denotes a collection of vectors; and
\item an operator type over the given scalar field, providing linear operators used to
define eigenproblems and preconditioners.
\end{itemize}
We note that the {\lib eigenSparse} was designed to support block methods,
defined as those that apply A or B
to a collection of vectors (a multivector).
Algorithms in {\lib eigenSparse} are developed at a high-level, where the underlying
linear algebra objects are opaque. The choice in linear algebra is made through templating,
and access to the functionality of the underlying objects is provided via the traits classes
{\class MultiVecTraits} and {\class OperatorTraits}. \\
These classes define opaque interfaces, specifying the operations that multivector and
operator classes must support in order to be used in {\lib eigenSparse} without exposing low-level
details of the underlying data structures.\\
The {\class MultiVecTraits} and {\class OperatorTraits} classes specify the operations that
the multivector and operator type must support in order for them to be used
by {\lib eigenSparse}. \\
The {\class OperatorTraits} class only needs to provide one method, described in the Table, that
applies an operator to a multivector. This interface defines the only interaction
required from an operator, even though the underlying operator may be a matrix,
spectral transformation, or preconditioner. \\

\begin{center}
\begin{tabular}{|l|p{11.5cm}|}
\hline
\multicolumn{2}{|c|} {\textbf{OperatorTraits<ST,MV,OP>}} \\
\hline
\textit{Method name} & \textit{Description}\\
\hline
{\cmd apply}(A,X,Y) & Applies the operator A to the multivector X, placing the result in the multivector Y \\
\hline
\end{tabular}
\end{center}

The methods defined by the {\class MultiVecTraits} class, listed in following table, are the nesseary
creational and arithmetic ones.
The creational methods generate empty or populated multivectors from a previously
created multivector. The populated multivectors can be a deep copy, where the
object contains the storage for the multivector entries, or a shallow copy, where the
object has a view of another multivector’s storage. A shallow copy is useful when
only a subset of the columns of a multivector is required for computation, which
is a situation that commonly occurs during the generation of a subspace.
\begin{center}
\begin{tabular}{|l|p{11.5cm}|}
\hline
\multicolumn{2}{|c|} {\textbf{MultiVecTraits<ST,MV>}} \\
\hline
\textit{Method name} & \textit{Description}\\
\hline
{\cmd clone(X,numvecs)}  & Creates a new multivector from X with {\itshape numvecs} vectors \\
{\cmd cloneCopy(X,index)}  & Creates a new multivector with a copy of the contents of a subset of the multivector X (deep copy) \\
{\cmd cloneView(X,index)}  & Creates a new multivector that shares the selected contents of a subset of the multivector X (shallow copy)  \\
\hline
{\cmd getVecLength(X)} & Returns the vector length of the multivector X  \\
{\cmd getNumberVecs(X)} & Returns the number of vectors in the multivector X \\
\hline
{\cmd mvTimesMatAddMv}(alpha,X,  & Applies a dense matrix D to multivector X and \\
\noindent D,beta,Y)                                 &  accumulates the result into multivector Y: \\
						  	        & $Y \leftarrow \alpha XD+\beta Y$ \\
{\cmd mvAddMv}(alpha,X,beta,Y) & Performs multivector: $Y \leftarrow \alpha X+ \beta Y$ \\
{\cmd mvTransMv}(alpha,X,Y,D) & Computes the dense matrix $D \leftarrow \alpha X^{H}Y$ \\
{\cmd mvDot}(X,Y,d) & Computes the corresponding dot products: $d[i] \leftarrow \bar a _{i} y _{i}$ \\
{\cmd mvScale}(X,d) & Scales the i-th column of a multivector X by $ d[i] $ \\
{\cmd mvNorm}(X,d) & Computes the 2-norm of each vector of X : $d[i] \leftarrow ||x _{i}|| _{2}$\\
\hline
{\cmd setBlock}(X,Y,index)& Copies the vectors in X to a subset of vectors in Y \\
{\cmd mvInit}(X,alpha)& Replaces each entry in the multivector X with a scalar $\alpha$ \\
{\cmd mvRandom}(X)& Replaces the entries in the multivector X by random scalars \\
\hline
{\cmd mvPrint}(X) & Print the multivector X \\
\hline
\end{tabular}
\end{center}
The arithmetic methods defined by the {\class multiVecTraits} are essential to the computations required by the Rayleigh-Ritz method and the general eigen-iteration. The {\cmd mvTimesMatAddMv} and {\cmd mvAddMv} methods, for instance, are necessary for updating the approximate eigenpairs and their residuals. {\cmd The mvDot} and {\cmd mvTransMv} methods are required by the orthogonalization procedure of the eigen-iteration. The {\cmd mvScale} and {\cmd mvNorm} methods are  necessary, at the very least, for the computation of approximate eigenpairs and for some termination criteria of the eigen-iteration. Deflation and locking of converged
eigenvectors necessitates the {\cmd setBlock} method in many cases. Initialization of the bases for the eigen-iteration requires methods such as {\cmd mvRandom} and {\cmd mvInit}. The ability to perform error checking and debugging is supported by methods that give dimensional attributes ({\cmd getVecLength}, {\cmd getNumberVecs}) and allow the users to print out a multivector ({\cmd mvPrint}).\\
Calling methods of {\class MultiVecTraits} and {\class OperatorTraits} requires that specializations of these traits classes have been implemented for given template arguments. \xlifepp provides the following classes that can be used with these traits classes:
\begin{itemize}
\item {\class MultiVectorAdapter},  an extended version of {\class Vector}, allowing the process on block of vector.
\item {\class LargeMatrixAdapter}, a wrapper of {\class LargeMatrix} in order to take adavantage of multiplication matrix-vector of sparse matrix
\end{itemize}

\subsubsection{The {\classtitle MultiVectorAdapter} class}

This class provides the implementations of the interfaces along with abstract {\class MultiVec} class. By inheriting directly from {\class MultiVec} class, the implementaion works via run-time polymorphism. This option isn't optimal from performance point of view because of the dispatching cost during run-time. In the near future, we can improve it by implementing a specialization of {\class  MultiVecTraits} class, which allows compile-time polymorphism. \\
By using the pointer to each vector,  {\class MultiVectorAdapter} class allows to create different objects, which share a same view of memory (shallow copy). This feature is essential when we want only to work on specific column vectors of a multivector.
\begin{lstlisting}[deletekeywords={[3] length}]
template <class ScalarType>
class MultiVecAdapter : public MultiVec<ScalarType>
{
  public:
    typedef std::vector<ScalarType>* VectorPointer;

  //! Constructor for a \c numberVecs vectors of length \c length.
  MultiVecAdapter(const Number length, const Dimen numberVecs) :
    length_(length),
    numberVecs_(numberVecs)
  {
    check();

    data_.resize(numberVecs);
    ownership_.resize(numberVecs);

    // Allocates memory to store the vectors.
    for (Dimen v = 0 ; v < numberVecs_ ; ++v)
    {
      data_[v] = new std::vector<ScalarType>(length);
      ownership_[v] = true;
    }

    // Initializes all elements to zero.
    mvInit(0.0);
  }
 ...
\end{lstlisting}

\subsubsection{The {\classtitle LargeMatrixAdapter} class}

This class is simply a wrapper of {\class LargeMatrix} class in order to take advantage of consuming-time matrix-vector multiplication provided by \xlifepp. The method {\cmd apply} implemnt multiplication of a LargeMatrix with a block of vector. Noting that it's only a way to implement application of an operator to vector. Depending on specific eigenproblem, we can create a new class with a different approach to apply operator; for example, operator in case of spectral transformation to find the not-well separated eigenvalues.
\begin{lstlisting}[deletekeywords={[3] x, y}]
  template<typename ScalarType>
  class LargeMatrixAdapter : public Operator<ScalarType>
  {
  private:
      typedef LargeMatrix<ScalarType> LMatrix;
  public:
      LargeMatrixAdapter() { lMatrix_p = NULL; }
      LargeMatrixAdapter(LMatrix* p) { lMatrix_p = p; }

      LargeMatrixAdapter& operator=(LMatrix matrix) { lMatrix_p = &matrix; return *this; }


      //! Destructor.
      virtual ~LargeMatrixAdapter() {};

      void apply(const MultiVec<ScalarType>& x, MultiVec<ScalarType>& y) const
      {
          for (int i = 0; i< x.getNumberVecs(); ++i) {
              multMatrixVector(*lMatrix_p, *(x[i]), *(y[i]));
          }
      }
  private:
      LMatrix*  lMatrix_p;
  };
\end{lstlisting}

\subsubsection{The {\classtitle LargeMatrixAdapterInverse} class}

Similar to {\class LargeMatrixAdapter}, this class is simply a wrapper of {\class LargeMatrix} in order to take advantage of consuming-time matrix-vector multiplication provided by \xlifepp. The class operates on an already factorized matrix to carry out the multiplication of its inverse (if invertible) with a vector. Instead of invoking the function {\cmd multMatrixVector}, the method {\cmd apply} calls the routine {\cmd multInverMatrixVector} of class {\class LargeMatrix} with the type of factorization.

\begin{lstlisting}[deletekeywords={[3] x, y}]
  template<typename MatrixType, typename ScalarType>
  class LargeMatrixAdapterInverse : public Operator<ScalarType>
  {
  private:
      typedef const MatrixType LMatrix;
  public:
//      LargeMatrixAdapterInverse() : lMatrix_p(NULL), fac_(_noFactorization) {}
      LargeMatrixAdapterInverse(LMatrix* p, FactorizationType fac) : lMatrix_p(p), fac_(fac) {}

      LargeMatrixAdapterInverse& operator=(LMatrix matrix) { lMatrix_p = &matrix; fac_ = _noFactorization; return (*this); }


      //! Destructor.
      virtual ~LargeMatrixAdapterInverse() {};

      void apply(const MultiVec<ScalarType>& x, MultiVec<ScalarType>& y) const
      {
          for (int i = 0; i< x.getNumberVecs(); ++i) {
              multInverMatrixVector(*lMatrix_p, *(x[i]), *(y[i]), fac_);
          }
      }

  private:
      LMatrix*  lMatrix_p;
      FactorizationType fac_;
  };
\end{lstlisting}

\subsubsection{The {\classtitle LargeMatrixAdapterInverseGeneralized} class}

On some specific eigen problems, it is better to do spectral transformation, especially, to find not well-seperated eigenvalues. This class serves for this purpose. As indicated by its name, this class works mainly for the spectral transformation of generalized eigen problem. The first operation is multiplication of the one matrix with a vector. The vector result is multiplied by the inverse of the second matrix which is provided in a form of factorized one. Both operations are done in the routine {\cmd apply}.
\begin{lstlisting}[deletekeywords={[3] x, y}]
  template<typename MatrixType, typename ScalarType>
  class LargeMatrixAdapterInverseGeneralized : public Operator<ScalarType>
  {
  private:
      typedef const MatrixType LMatrix;
  public:
      LargeMatrixAdapterInverseGeneralized(LMatrix* pA, LMatrix* pB, FactorizationType fac) : lMatrixA_p(pA), lMatrixB_p(pB), fac_(fac) {}

//      LargeMatrixAdapterInverseGeneralized& operator=(LMatrix matrix) { lMatrix_p = &matrix; fac_ = _noFactorization; return (*this); }


      //! Destructor.
      virtual ~LargeMatrixAdapterInverseGeneralized() {};

      void apply(const MultiVec<ScalarType>& x, MultiVec<ScalarType>& y) const
      {
          Vector<ScalarType> vecTemp(x.getVecLength());
          for (int i = 0; i< x.getNumberVecs(); ++i) {
              multMatrixVector(*lMatrixB_p, *(x[i]), vecTemp);
              multInverMatrixVector(*lMatrixA_p, vecTemp, *(y[i]), fac_);
          }
      }

  private:
      LMatrix*  lMatrixA_p;
      LMatrix*  lMatrixB_p;
      FactorizationType fac_;
  };
\end{lstlisting}

\subsection{{\libtitle eigenSparse} framework}

{\lib eigenSparse} is organized with a multi-tiered access for eigensolver algorithms. Users have the choice of interfacing at one of two levels: either working at a high-level with a eigensolver manager or working at a low-level directly with an eigensolver.\\
Consider as an example the block Davidson iteration. \\
\begin{tabular} {l}
\hline
Block Davidson Algorithm \\
\hline
Require: Set an initial guess \textbf{V} and \textbf{H} = [] \\
1: \textbf{for} $iter=1$  to $iter _{max}$ \textbf{do}\\
2: ~ \textbf{repeat} \\
3: ~  ~ ~ ~Compute \textbf{M}-orthonormal basis \textbf{V} for [\textbf{V},\textbf{H}]\\
4: ~  ~ ~ ~Project \textbf{A} onto \textbf{V}: $\hat A = V^{H}AV$ \\
5: ~  ~ ~ ~Compute selected eigenpairs (\textbf{Q},\textbf{$\Gamma$}) of $\hat A$: $\hat A Q = Q\Gamma $ \\
6: ~  ~ ~ ~Compute Ritz vectors: $X = VQ$ \\
7: ~  ~ ~ ~Compute residuals: $R = AX - MX\Gamma$\\
8: ~  ~ ~ ~Precondition the residuals: $H = NR$ \\
9: ~  ~ ~ ~Check convergence and possibly terminate iteration  \\
10:~ \textbf{until} the matrix \textbf{V} is not expandable \\
11: ~Restart \textbf{V} \\
12: \textbf{end for} \\
\hline
\end{tabular}

\medskip

\begin{figure} [H]
\centering
\includeTikZPict[width=13cm]{solverManager.png}
\caption{\textbf{SolverManager} class collaboration graph}
\end{figure}

The linear operators \textbf{A} and \textbf{M} define the eigenproblem to be solved. The linear operator \textbf{N}
is a preconditioner for the problem, and its application to a block of vectors is required. Example preconditioners \textbf{N}
include the inverse of the diagonal of \textbf{A} (Jacobi preconditioner), an algebraic multigrid preconditioner for \textbf{A}, or
a preconditioned iteration that approximates the solution of a linear set of equations with \textbf{A}, e.g., the preconditioned conjugate gradient algorithm. \\
The \textit{multivector} \textbf{V}, \textbf{X}, and \textbf{R} serve for some purpose. The column-vectors for matrix \textbf{V}
form the basis for the Rayleigh-Ritz approximation conducted at Step 5-6. We make the specific choice that these vectors are orthonormal with respect to the inner product induced by the Hermitian positive-definite matrix \textbf{M}, but this is not a requirement. \\
The block Davidson eigensolver as described above is useful for examining some
of the functionality provided by {\lib eigenSparse}. The algorithm above highlights three levels of
functionality. The first level is given by Steps 1-12 that constitute the eigensolver
strategy to solve problem (2). The second level is given by Steps 2-10 that form
the eigen-iteration. The third level consists of computational steps that can be
implemented in a variety of manners, so encouraging modularization. Step 3 requires an orthonormalization method. The decisions involved in Step 5 require a
determination of the eigenvalues and invariant subspace of interest, as well as a
definition of accuracy. To check convergence in Step 9, several criteria are possible.
For instance, a norm induced by a matrix other than \textbf{M} may be employed. The
restarting of this eigensolver (Step 11) can be performed in a variety of ways and
therefore need not be tightly coupled to the eigensolver. Each of these mechanisms
provide opportunity for decoupling functionality that need not be implemented in
a specific manner. \\
Looking deeper into the algorithm, the iteration repeats until the basis V is full
(in which case it is time to restart) or some stopping criterion has been satisfied. Many valid
stopping criteria exist, as well as many different methods for restarting the basis. Both of
these, however, are distinct from the essential iteration as described above. A user wanting
to perform block Davidson iterations could ask the solver to perform these iterations until
a user-specified stopping criterion was satisfied or the basis was full, at which time the user
would perform a restart. This allows the user complete control over the stopping criteria and the restarting mechanism, and leaves the eigensolver responsible for a relatively simple bit of state and behavior. \\
The eigensolvers (encapsulating an iteration and the associated state) are derived classes of the abstract base class {\class Eigensolver.}
The goals of this class are three-fold:
\begin{itemize}
\item to define an interface used for checking the status of a solver by a status test;
\item to contain the essential iteration associated with a particular eigensolver iteration;
\item to contain the state associated with that iteration.
\end{itemize}
The status tests, assembled to describe a specific stopping criterion and queried by the
eigensolver during the iteration, are represented as subclasses of {\class StatusTest}. The
communication between status test and eigensolver occurs inside of the \textit{iterate()} method
provided by each {\class Eigensolver}. This code generally takes the form:
\begin{lstlisting}[]{}
SomeEigensolver::iterate() {
while ( statustest.checkStatus(this) != _passed ) {
//
// perform eigensolver iterations
//
}
r
\end{lstlisting}
Each {\class StatusTest} provides a method, \textit{checkStatus()}, which through queries to
the methods provided by {\class Eigensolver}, determines whether the solver meets the criteria defined by that particular status test. After a solver returns from \textit{iterate()}, the caller has the option of accessing the state associated with the solver and re-initializing the solver with a new state. \\
While this method of interfacing with the solver is powerful, it can be tedious. This method requires that user construct a number of support classes, in addition to managing call to \textit{Eigensolver::iterate()}. The {\class SolverManager} class was developed to
address this complaint. A solver manager is a class that wraps around an eigensolver, providing additional functionality while also handling lower-level interaction with the eigensolver that a user may not wish to handle. Solver managers are intended to be easy to
use, while still providing the features and flexibility needed to solve real-world eigenvalue problems. For example, the {\class BlockDavidsonSolMgr} takes only two arguments in its constructor: an {\class Eigenproblem} specifying the eigenvalue problem to be solved and a {\class Parameters} of options specific to this solver manager. The solver manager instantiates an eigensolver, along with the status tests and other support classes needed by the eigensolver. To solve the eigenvalue problem, the user simply calls the \textit{solve()} method of the solver manager. The solver manager performs repeated calls to the eigensolver \textit{iterate()}
method, performs restarts and locking, and places the final solution into the eigenproblem.

\subsection{{\libtitle eigenSparse} classes}

The following lists and briefly describes the current classes of {\lib eigenSparse}.

\subsubsection{The {\classtitle EigenProblem} class}

{\class EigenProblem} is a container for the components of an eigenvalue problem, as well as the solutions. By requiring eigen problems to derive from {\class EigenProblem}, {\lib eigenSparse} defines a minimum interface that can be expected of all eigenvalue problems by the classes
that will work with the problems (e.g., eigensolvers and status testers).\\
\begin{lstlisting}[]{}
/*!
 * \class xlifepp::EigenProblem
  \brief This provides a basic implementation for defining standard or
  generalized eigenvalue problems.
*/
  template<class ScalarType, class MV, class OP>
  class EigenProblem
  {
  public:

    //! Empty constructor - allows xlifepp::EigenProblem to be described at a later time through "Set Methods".
    EigenProblem();

    //! Standard Eigenvalue Problem Constructor.
    EigenProblem(const SmartPtr<const OP>& Op, const SmartPtr<MV>& initVec);

    //! Generalized Eigenvalue Problem Constructor.
    EigenProblem(const SmartPtr<const OP>& Op, const SmartPtr<const OP>& B, const SmartPtr<MV>& initVec);

    //! Copy Constructor.
    EigenProblem(const Eigenproblem<ScalarType, MV, OP>& Problem);

    //! Destructor.
    ~EigenProblem() {};
    ...
\end{lstlisting}
Both the eigen problem and the eigensolver in {\lib eigenSparse} are templated on the scalar type,
the multivector type and the operator type. Before declaring an eigen problem, users must
choose classes to represent these entities. Having done so, they can begin to specify the
parameters of the eigen problem. The {\class EigenProblem} class defines set methods for
the parameters of the eigen problem. These methods are:
\begin{itemize}
\item {\cmd setOperator} - set the operator for which the eigenvalues will be computed
\item {\cmd setA} - set the A operator for the eigenvalue problem $Ax = Mx\lambda$
\item {\cmd setM} - set the M operator for the eigenvalue problem $Ax = Mx\lambda$
\item {\cmd setPrec} - set the preconditioner for the eigenvalue problem
\item {\cmd setInitVec} - set the initial iterate
\item {\cmd setAuxVecs} - set the auxiliary vectors, a subspace used to constrain the search space for the solution
\item {\cmd setNEV} - set the number of eigenvalues to be computed
\item {\cmd setHermitian} - specify whether the problem is Hermitian
\end{itemize}
In addition to these {\itshape set} methods, {\class EigenProblem} defines a method {\cmd setProblem()} that gives the class the opportunity to perform any initialization that may be necessary before the problem is handed off to an eigensolver, in addition to verifying that the problem has been adequately defined. \\
For each of the {\itshape set} methods listed above, there is a corresponding {\itshape get} function. These are the functions used by eigensolvers and solver managers to get the necessary information from the eigenvalue problem. In addition, there are two methods for storing and retrieving the results of the eigenvalue computation:
\begin{lstlisting}[]{}
const EigenSolverSolution & getSolution();
void setSolution(const EigenSolverSolution & sol);
\end{lstlisting}
In order to verify if all the parameters are set, {\class EigenProblem} provide a method {\cmd isProblemSet()}

\subsubsection{The {\classtitle EigenSolverSolution} structure}

The {\class EigenSolution} structure was developed in order to facilitate setting and retrieving of solution data.
\begin{lstlisting}[]{}
  //!  Struct for storing an eigenproblem solution.
  template <class ScalarType, class MV>
  struct EigenSolverSolution {
    //! The computed eigenvectors
    SmartPtr<MV> evecs;
    //! An orthonormal basis for the computed eigenspace
    SmartPtr<MV> espace;
    //! The computed eigenvalues
    std::vector<ValueEigenSolver<ScalarType> >  evals;
    /*! \brief An index into evecs to allow compressed storage of eigenvectors for real, non-Hermitian problems.
     *
     *  index has length numVecs, where each entry is 0, +1, or -1. These have the following interpretation:
     *     - index[i] == 0: signifies that the corresponding eigenvector is stored as the i column of evecs. This will usually be the
     *       case when ScalarType is complex, an eigenproblem is Hermitian, or a real, non-Hermitian eigen problem has a real eigenvector.
     *     - index[i] == +1: signifies that the corresponding eigenvector is stored in two vectors: the real part in the i column of evecs and the <i><b>positive</b></i> imaginary part in the i+1 column of evecs.
     *     - index[i] == -1: signifies that the corresponding eigenvector is stored in two vectors: the real part in the i-1 column of evecs and the <i><b>negative</b></i> imaginary part in the i column of evecs
     */
    std::vector<int>         index;
    //! The number of computed eigenpairs
    int numVecs;

    EigenSolverSolution() : evecs(),espace(),evals(0),index(0),numVecs(0) {}
  };
\end{lstlisting}

All {\lib eigenSparse} solver managers place the results of the computation in the {\class EigenProblem} class using an {\class EigenSolverSolution} structure. The number of eigen solutions computed is given by field numVecs. The eigenvalues are always stored as two real values, even
when templated on a complex data type or when the eigenvalues are real. Similarly, the eigen space can always be represented by a multivector of width numVecs, even for non-symmetric eigen problems over the real field. The storage scheme for eigenvectors requires
more finesse.\\
When solving real symmetric eigen problems, the eigenvectors can always be chosen to
be real, and therefore can be stored in a single column of a real multivector. When solving
eigenproblems over a complex field, whether Hermitian or non-Hermitian, the eigenvectors may be complex, but the multivector is defined over the complex field, so that this poses no problem. However, real non-symmetric problems can have complex eigenvectors, which prohibits a one-for-one storage scheme using a real multivector. Fortunately, the eigenvectors
in this scenario occur as complex conjugate pairs, so the pair can be stored in two real
vectors. This permits a compressed storage scheme, which uses an index vector stored in the
{\class EigenSolverSolution}, allowing conjugate pair eigenvectors to be easily retrieved from evecs.
The integers in {\cmd EigenSolverSolution::index} take one of three values: $\{0, +1, -1\}$. These values allow the eigenvectors to be retrieved as follows:
\begin{itemize}
\item $index[i]=0$: the i-th eigenvector is stored uncompressed in column $i$ of evecs
\item $index[i] = +1$: the i-th eigenvector is stored compressed, with the real component in column $i$ of evecs and the positive complex component stored in column $i$ + 1 of evecs
\item $index[i] = -1$: the i-th eigenvector is stored compressed, with the real component in
column $i$ -1 of evecs and the negative complex component stored in column $i$ of evecs
\end{itemize}
Because this storage scheme is only required for non-symmetric problems over the real field,
all other eigenproblems will result in an index vector composed entirely of zeroes. For the
real non-symmetric case, the +1 index will always immediately precede the corresponding -1
index. \\

\subsubsection{The {\classtitle EigenSolver} class}

The {\class EigenSolver} class defines the basic interface that must be met by any eigen solver class in {\lib eigenSparse}. The specific eigensolvers are implemented as derived classes of {\class EigenSolver}. There are two eigen solvers currently implemented in {\lib eigenSparse}
\begin{itemize}
\item {\class BlockDavidson} : A block Davidson solver for Hermitian eigenvalue problems
\item {\class BlockKrylovSchur} : A block Krylov Schur solver for Hermitian or non-Hermitian
eigenvalue problems.
\end{itemize}
The class {\class EigenSolver}, like {\class EigenProblem}, is templated on the scalar type, multivector type and operator type. The options for the eigen solver are passed through the constructor, defined by {\class EigenSolver} to have the following form:
\begin{lstlisting}[]{}
template<class ScalarType, class MV, class OP>
class EigenSolver {
  public:

  //! @name Constructors/Destructor
  //@{

  //! Default Constructor.
  EigenSolver() {};

  //! Basic Constructor.
  /*! This constructor, implemented by all xlifepp eigensolvers, takes an xlifepp::Eigenproblem,
    xlifepp::SortManager, xlifepp::OutputManager, and Parameters as input.  These
    four arguments are sufficient enough for constructing any xlifepp::EigenSolver object.
  */
  EigenSolver( const SmartPtr<EigenProblem<ScalarType,MV,OP> >& problem,
               const SmartPtr<SortManager<ScalarType> >&        sorter,
               const SmartPtr<OutputManager<ScalarType> >&      printer,
               const SmartPtr<StatusTest<ScalarType,MV,OP> >&   tester,
               const SmartPtr<OrthoManager<ScalarType,MV> >&    ortho,
               Parameters& params );
...
\end{lstlisting}
These classes are used as follows:
\begin{itemize}
\item \textbf{problem} - the eigenproblem to be solved; the solver will get the problem operators from here.
\item \textbf{sorter} - the sort manager selects the significant eigenvalues
\item \textbf{printer} - the output manager dictates verbosity level in addition to processing output streams
\item \textbf{tester} - the status tester dictates when the solver should quit iterate() and return to the caller
\item \textbf{ortho} - the orthogonalization manager defines the inner product and other concepts
related to orthogonality, in addition to performing these computations for the solver
\item \textbf{params} - the parameter list specifies eigensolver-specific options; see the documentation
for a list of options support by individual solvers
\end{itemize}
In addition to specifying an iteration, an eigensolver also specifies a concept of state, i.e.  the current data associated with the iteration. After declaring an Eigensolver object, the state of the solver is in an uninitialized state. For most solver, to be initialized mean to be in
a valid state, containing all of the information necessary for performing eigen solver iterations. \\
Eigensolver provides two methods concerning initialization: \textit{isInitialized()} indicates whether the solver is initialized or not, and \textit{initialize()} (with no arguments) instructs the solver to initialize itself using random data or the initial vectors stored in the eigenvalue problem. \\
To ensure that solvers can be used as efficiently as possible, the user needs access to the
state of the solver. To this end, each eigensolver provides low-level methods for getting and
setting the state of the solver:
\begin{itemize}
\item {\cmd getState()} - returns a solver-specific structure with read-only pointers to the current
state of the solver.
\item {\cmd initialize(...)} - accepts a solver-specific structure enabling the user to initialize the
solver with a particular state.
\end{itemize}
The combination of these two methods, along with the flexibility provided by status tests,
allows the user almost total control over eigen solver iterations. \\
Moreover, each eigen solver provides two significant types of methods: status methods and solver-
specific state methods. The status methods are defined by the Eigensolver abstract
base class and represent the information that a generic status test can request from any
eigensolver. A list of these methods is given in following table
\begin{center}
\begin{tabular}{|l|p{11.5cm}|}
  \hline
  \textit{Method name} & \textit{Description}\\
  \hline
  {\cmd getNumIters} & Get the current number of iterations. \\
  {\cmd getRitzValues} & Get the most recent Ritz values. \\
  {\cmd getRitzVectors} & Get the most recent Ritz vectors. \\
  {\cmd getRitzIndex} & Get the Ritz index needed for indexing compressed Ritz vectors. \\
  {\cmd getResNorms} & Get the most recent residual norms, with respect to the OrthoManager. \\
  {\cmd getRes2Norms} & Get the most recent residual 2-norms. \\
  {\cmd getRitzRes2Norms} & Get the most recent Ritz residual 2-norms. \\
  {\cmd getCurSubspaceDim} & Get the current subspace dimension. \\
  {\cmd getMaxSubspaceDim} & Get the maximum subspace dimension. \\
  {\cmd getBlockSize} & Get the block size. \\
  \hline
\end{tabular}
\end{center}

\subsubsection{The {\classtitle SolverManager} class}

Using {\lib eigenSparse} by interfacing directly with eigensolvers is extremely powerful, but can be
tedious. Solver managers provide a way for users to encapsulate specific solving strategies
inside of an easy-to-use class. Novice users may prefer to use existing solver managers, while
advanced user may prefer to write custom solver managers.
{\class SolverManager} defines only two methods: a constructor accepting an {\class EigenProblem} and a parameter list of solver manager-specific options; and a {\cmd solve()} method, taking no arguments and returning either {\tt \_converged} or {\tt \_unconverged}. Consider the following example code:
\begin{lstlisting}[]{}
// Create an eigen problem
SmartPtr<EigenProblem<ST,MV,OP> > problem =
// Create parameter list to pass into the solver manager
  Parameters pl;
  pl.set("Verbosity", (int)verbosity);
  pl.set("Which", which);
  pl.set("Orthogonalization", ortho);
  pl.set("Block Size", blockSize);
  pl.set("Num Blocks", numBlocks);
  pl.set("Maximum Restarts", maxRestarts);
// Create a solver manager
  BlockDavidsonSolMgr<ST,MV,OP> bdSolverMan(problem, pl);
// Solve the eigenvalue problem
  ComputationInfo returnCode = bdSolverMan.solve();
// Get the solution from the problem
  EigenSolverSolution<ST,MV> sol = problem->getSolution();
\end{lstlisting}
As has been stated before, the goal of the solver manager is to create an eigensolver object,
along with the support objects needed by the eigensolver. Another purpose of many solver
managers is to manage and initiate the repeated calls to the underlying solver’s {\cmd iterate()}
method. For solvers that build a Krylov subspace to some maximum dimension (e.g., {\class BlockKrylovSchur} and {\class BlockDavidson}), the solver manager will also assume the task of restarting the solver when the subspace is full. This is something for which multiple approaches are possible. Also, there may be substantial flexibility in creating the support classes (e.g., sort
manager, status tests) for the solver. An aggressive solver manager could even go so far as to
construct a preconditioner for the eigenvalue problem. \\
These examples are meant to illustrate the flexibility that specific solver managers may
have in implementing the {\cmd solve()} routine. Some of these options might best be incorporated
into a single solver manager, which takes orders from the user via the parameter list given in
the constructor. Some of these options may better be contained in multiple solver managers,
for the sake of code simplicity. It is even possible to write solver managers that contain other
solvers managers; motivation for something like this would be to select the optimal solver
manager at runtime based on some expert knowledge, or to create a hybrid method which
uses the output from one solver manager to initialize another one. \\

\subsubsection{The {\classtitle StatusTest} class}

the purpose of the {\class StatusTest} should be clear: to give the user or solver manager flexibility in stopping the eigen solver iterations in order to interact directly with the solver. \\
Many reasons exist for why a user would want to stop the solver from iterating:
\begin{itemize}
\item  some convergence criterion has been satisfied and it is time to quit;
\item some part of the current solution has reached a sufficient accuracy to removed from the
iteration (“locking”);
\item the solver has performed a sufficient or excessive number of iterations.
\end{itemize}
These are just some commonly seen reasons for ceasing the iteration, and each of these can be so varied in implementation/parametrization as to require some abstract mechanism controlling the iteration. \\
Below is a list of status tests:
\begin{itemize}
\item {\class StatusTestCombo} - this status test allows for the boolean combination of other
status tests, creating near unlimited potential for complex status tests.
\item {\class StatusTestOutput} - this status test acts as a wrapper around another status
test, allowing for printing of status information on a call to checkStatus()
\item {\class StatusTestMaxIters} - this status test monitors the number of iterations performed by the solver; it can be used to halt the solver at some maximum number of
iterations or even to require some minimum number of iterations.
\item {\class StatusTestResNorm} - this status test monitors the residual norms of the current iterate.
\item {\class StatusTestWithOrdering} - this status test also monitors the residual norms
of the current iterate, but only considers the residuals along with a set of auxiliary eigenvalues.
\end{itemize}

\subsubsection{The {\classtitle SortManager} class}

The purpose of a sort manager is to separate the eigensolver classes from the sorting functionality required by those classes. This satisfies the flexibility principle sought by {\lib eigenSparse}, by giving users the opportunity to perform the sorting in whatever manner is deemed to be most appropriate. {\lib eigenSparse} defines an abstract class {\class SortManager} with two methods,
one for sorting real values and one for sorting complex values:
\begin{lstlisting}[]{}
// Sorting with real value
virtual void sort(std::vector<MagnitudeType>& evals, SmartPtr<std::vector<int> > perm = _smPtrNull, int n = -1) const = 0;
// Sorting with complex value
virtual void sort(std::vector<MagnitudeType>& rEvals,
                          std::vector<MagnitudeType>& iEvals,
                          SmartPtr<std::vector<int> > perm = _smPtrNull,
                           int n = -1) const = 0;
\end{lstlisting}

\subsubsection{The {\classtitle BasicSort} class}

The concrete derived class implements sorting scheme for sorting real values and one for sorting complex values.
\begin{lstlisting}[]{}
template<class MagnitudeType>
class BasicSort : public SortManager<MagnitudeType>
{
public:
    /*! \brief Parameter list driven constructor

        This constructor accepts a paramter list with the following options:
       - \c "Sort Strategy" - a \c string specifying the desired sorting strategy. See setSortType() for valid options.
    */
    BasicSort(Parameters& pl);

    /*! \brief String driven constructor

        Directly pass the string specifying sort strategy. See setSortType() for valid options.
    */
    BasicSort(const String& which = "LM");

    //! Destructor
    virtual ~BasicSort();
    ...
\end{lstlisting}
By default, sort strategy Largest Magnitude (LM) is used. So as to change the strategy, user can easily specify it with {\cmd setSortType} method.

\subsubsection{The {\classtitle OrthoManager} class}

Orthogonalization and orthonormalization are commonly performed computations in iterative eigensolvers; in fact, for some eigensolvers, they represent the dominant cost. Different scenarios may require different approaches (e.g, Euclidean inner product versus $M$ inner
product, orthogonal projections versus oblique projections). Combined with the plethora of
available methods for performing these computations, Anasazi has left as much leeway to the
users as possible. \\
Orthogonalization of multivectors is performed by derived classes of the abstract class {\class OrthoManager} that provides five methods:
\begin{itemize}
\item {\cmd innerProd(X,Y,Z)} - performs the inner product defined by the manager.
\item {\cmd norm(X)} - computes the norm induced by {\cmd innerProd()}.
\item {\cmd project(X,C,Q)} - given an orthonormal basis Q, projects X onto to the space perpindicular to colspan(Q), optionally returning the coefficients of X in Q.
\item {\cmd normalize(X,B)} - returns an orthonormal basis for colspan(X), optionally returning
the coefficients of X in the computed basis.
\item {\cmd projectAndNormalize(X,C,B,Q)} - computes an orthonormal basis for subspace
colspan(X) - colspan(Q), optionally returning the coefficients of X in Q and the new
basis.
\end{itemize}
It should be noted that a call to {\cmd projectAndNormalize()} is not necessarily equivalent
to a call to {\cmd project()} followed by {\cmd normalize()}. This follows from the fact that, for some
orthogonalization managers, a call to {\cmd normalize()} may augment the column span of a rank-
deficient multivector in order to create an orthonormal basis with the same number of columns
as the input multivector. In this case, the code
\begin{lstlisting}[]{}
orthoman.project(X,C,Q);
orthoman.normalize(X,B);
\end{lstlisting}
could result in an orthonormal basis $X$ that is not orthogonal to the basis in $Q$. \\
\xlifepp provides two orthogonalization managers:
\begin{itemize}
\item {\class BasicOrthoManager} - performs orthogonalization using multiple steps of classical Gram-Schmidt.
\item {\class SVQBOrthoManager} - performs orthogonalization using the SVQB orthogonalization technique described by Stathapoulos and Wu
\end{itemize}
The core function of these orthogonalization managers is {\cmd findBasis}, which find an operator-orthonormal basis for a span, with the option of extending the subspace.

\subsubsection{The {\classtitle OutputManager} class}

The output manager exists to provide flexibility with regard to the verbosity of the eigen solver. Each output manager has two primary concerns: what output is printed and where the output is printed to. When working with the output manager, output is classified
into one of the following message types: \\
\begin{lstlisting}[]{}
enum MsgEigenType
{
  _errorsEigen = 0,               /*!< Errors [ always printed ] */
  _warningsEigen = 0x1,             /*!< Internal warnings */
  _iterationDetailsEigen = 0x2,     /*!< Approximate eigenvalues, errors */
  _orthoDetailsEigen = 0x4,         /*!< Orthogonalization/orthonormalization details */
  _finalSummaryEigen = 0x8,        /*!< Final computational summary */
  _timingDetailsEigen = 0x10,       /*!< Timing details */
  _statusTestDetailsEigen = 0x20,   /*!< Status test details */
  _debugEigen = 0x40               /*!< Debugging information */
};
 \end{lstlisting}
Output manager are subclasses of the abstract base class {\class OutputManager}. This class provides the following output-related methods:
\begin{lstlisting}[deletekeywords={[3] type}]
template <class ScalarType>
class OutputManager {
  public:
   ...
  //! Find out whether we need to print out information for this message type.
  /*! This method is used by the solver to determine whether computations are
      necessary for this message type.
  */
  virtual bool isVerbosity(MsgEigenType type) const = 0;

  //! Send output to the output manager.
  virtual void print(MsgEigenType type, const String output) = 0;

  //! Create a stream for outputting to.
  virtual std::ostream& stream(MsgEigenType type) = 0;
...
 \end{lstlisting}
\xlifepp provides a single output manager, {\class BasicOutputManager}. This class accepts an output stream from the user. The output corresponding to the verbosity level of the manager is sent to this stream.

\subsubsection{The {\classtitle SolverUtils} class}

Most of the large-scale eigenvalue algorithms exploit projection processes so as to extract approximate eigenvectors from a given subspace. Inside this subspace, a small matrix eigenvalue problem is obtained. \xlifepp provides {\class SolverUtils} class, a collection of utilities necessary for the solvers.  These utilities include sorting, orthogonalization, projecting/solving local eigen systems,
\begin{lstlisting}[]{}
  template<class ScalarType, class MV, class OP>
  class SolverUtils
  {
  public:
    typedef typename NumTraits<ScalarType>::RealScalar MagnitudeType;
    typedef typename NumTraits<ScalarType>::Scalar  SCT;
    typedef MatrixEigenDense<SCT> MatrixDense;

    //! @name Constructor/Destructor
    //@{

    //! Constructor.
    SolverUtils();

    //! Destructor.
    virtual ~SolverUtils() {};
    ...
 \end{lstlisting}
In fact, all utilities along with this class make use of functions of {\lib eigenCore}. The most important function of this class is {\cmd directSolver}, a routine for computing eigen pairs of Hermitian pencil. This method is heavily used during the projection of vectors into a subspace and normalization of these vector with a specified norm. \\

\subsection{Examples}

The following gives sample code for solving a symmetric eigenvalue problem using the Block Davidson solver manager.\\
The first step in solving an eigenvalue problem is to define the eigenvalue problem. However, before that, we need to make sure define the operator and the multivector.
\begin{lstlisting}[]{}
// Load LargeMatrix from file
LargeMatrix<Real> rMatSym(dataPathTo(rMatrixSym), _dense, 782, _cs, _symmetric);
// Define the operator (in this case, it's LargeMatrix)
SmartPtr< const LargeMatrixAdapter<ST> > K(new LargeMatrixAdapter<ST>(&rMatSym));
// Sometimes, we need a preconditioner to be able to solve the eigenvalue problem
SmartPtr< const LargeMatrixAdapter<ST> > Prec(new LargeMatrixAdapter<ST>(&matPrec));
// Define multivector
SmartPtr<MV> ivec = SmartPtr<MultiVecAdapter<ST> >(new MultiVecAdapter<ST>(rMatSym.nbCols, blockSize));
// Initialize the multivector  with random values
ivec->mvRandom();
 \end{lstlisting}
Everything is fine then eigenvalue problem is defined as follow:
\begin{lstlisting}
SmartPtr<EigenProblem<ST,MV,OP> > problem = SmartPtr<EigenProblem<ST,MV,OP> >(new EigenProblem<ST,MV,OP>(K,ivec));
// Inform the eigenProblem that the operator cK is symmetric
problem->setHermitian(true);
// Set the number of eigenvalues requested
problem->setNEV(nev);
\end{lstlisting}
The second step is to specify parameters for the solver manager, depending on algorithm (Block Davidson or Block Krylov Schur), we can have a different parameters. Nevertheless, there are some essential ones to be specified
\begin{lstlisting}[]{}
const int nev = 4; //6
int blockSize = 5;//6;
int numBlocks = 5;//14;
int maxRestarts = 25;//8;
Real tol = 1.0e-6;
String which("SM"), ortho("SVQB");
bool locking = true;
bool insitu = false;
MsgEigenType verbosity = _errorsEigen;

Parameters pl;
pl.set("Verbosity", (int)verbosity);
pl.set("Which", which);
pl.set("Orthogonalization", ortho);
pl.set("Block Size", blockSize);
pl.set("Num Blocks", numBlocks);
pl.set("Maximum Restarts", maxRestarts);
pl.set("Convergence Tolerance", tol);
pl.set("Use Locking", locking);
pl.set("Locking Tolerance", tol/10);
pl.set("In Situ Restarting", insitu);
\end{lstlisting}
Till now, we have all the information needed to declare the solver manager and solve the problem:
\begin{lstlisting}
BlockDavidsonSolMgr<ST,MV,OP> bdSolverMan(problem, pl);
ComputationInfo returnCode = bdSolverMan.solve();
 \end{lstlisting}
The return value of the solver indicate whether the algorithm succeeds or not.(i.e: whether the requested number of eigenpairs were found to a sufficient accuracy (as defined by the solver manager). Output from {\cmd solve()} routine in this example might look as follows:
\begin{lstlisting}[language=]
BlockDavidson Eigen Solver Test : Standard symmetric matrix test with Smallest Magnitude
--------------------------------------------------------------------------------
  Direct residual norms computed
          Eigenvalue         Residual(M)
----------------------------------------
        -1.29891e-06         8.63079e-08
        -1.29973e-06         6.27415e-08
        -1.30048e-06         7.99246e-08
        -1.33992e-06         7.09023e-07

 \end{lstlisting}