inputs/user/EigenArpack.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/>.
%%%%%%%%%%%%%%%%%%%

\subsection{Arpack++}

{\lib Arpack} is a well-known collection of Fortran77 subroutines designed to solve large scale eigenvalue problems. {\lib Arpack} implements a variant of the Arnoldi process for finding eigenvalues called implicit restarted Arnoldi method (IRAM). In most cases only a compressed matrix or a matrix-vector product $y \leftarrow Ax$ must be supplied by the user.\\
{\lib Arpack++} is a collection of classes that offers c++ programmers an interface to {\lib Arpack}. Because the main features of {\lib Arpack} are preserved by this C++ version, an interface between {\lib Xlife++} and {\lib Arpack++} is designed to enable the use of classes of {\lib Arpack++} to solve eigen problems.\\
In the following section, some instructions are given on how to install Arpack++. The next section discusses briefly how to solve an eigen problem with a basic interface {\lib Xlife++ - Arpack++}. The last one describes some advanced techniques to make use of all features of {\lib Arpack++}

\subsubsection{Basic Usage}

Two goal motivated the development of the interface between {\lib Xlife++} and {\lib Arpack++}: simplicity and homogeneity. The intention of simplicity is to ease the use of {\lib Arpack++} from inside {\lib Xlife++}. To this end, all the complicated parametrization of {\lib Arpack++} are hidden, and solver is invoked to calculate a eigen problem with only the essential parameters. However, user can make use of other parameters of eigen solver totally, which is described in the next section. The concept of homogeneity drives the development of interface to allow users to call {\lib Arpack++} solvers like other iterative solvers. The eigensolver is called with the operator () with similar prototype to other linear solvers. From now on, we call this interface eigensolver Arpack.\\
Like iterative solvers, the eigensolver Arpack is initialized with default value
\vspace{.1cm}
\begin{lstlisting}[]{}
// To solve standard eigenproblem A*x=l*x
EigenSolverStdArpackpp eigenSolverStd;
// To solve gneralized eigenproblem A*x=l*B*x
EigenSolverGenArpackpp eigenSolverGen;
\end{lstlisting}
\vspace{.1cm}
or with user-defined parameters
\vspace{.1cm}
\begin{lstlisting}[]{}
// Corresponding parameters
// EigenSolverStdArpackpp(const String& name, const Number maxOfIt, const Real epsilon)
// name: Name of Solver; maxOfIt: maximum of iteration number; epsilon: Tolerance
EigenSolverStdArpackpp eigenSolverStd("Standard",1000, 1e-8);
EigenSolverGenArpackpp eigenSolverGen("Gen", 1000, 1e-10);
\end{lstlisting}
\vspace{.1cm}
The following examples describe how to use eigensolver Arpack object to solve two  specfic eigen problem in different calculation modes.\\
Supposed that matrix $A$  and matrix $B$ are {\class LargeMatrix} with size $20\times20$. These two matrices define eigen problem. And we also have vector eigenvalue and vector of eigenvector containing computation results.
\vspace{.1cm}
\begin{lstlisting}[]{}
// Eigen problem size
const Number size = 20;
// Number of eigenvalue to retrieve
const Number nev = 5;
// Real sigma which is used in shift and invert mode
Real sigma = 10.0;
// Complex sigma which is used in shift and invert mode
Complex csigma(-700, 20.0);
// Real eigenvalue vector
Vector<Real> rEigValue(size);
// Complex eigenvalue vector
Vector<Complex> cEigValue(size);
// Real eigenvector vector
std::vector<Vector<Real> > rEigVector(nev, rEigValue);
// Complex eigenvector vector
std::vector<Vector<Complex> > cEigVector(nev, cEigValue);
// Standard Eigensolver
EigenSolverStdArpackpp eigenSolverStd("Standard",1000, 1e-8);
// Generalized Eigensolver
EigenSolverGenArpackpp eigenSolverGen("Gen", 1000, 1e-10);
\end{lstlisting}
\vspace{.1cm}

\paragraph*{The standard eigen problem $AX=\lambda X$}

Real symmetric case: Matrix $A$ is real and symmetric\\
Eigenvalues and eigenvectors are always real
\vspace{.1cm}
\begin{lstlisting}[]{}
//Regular mode
eigenSolverStd(A, rEigValue, rEigVector, nev, _regular, sigma);
out << rEigValue << std::endl;
out << rEigVector << std::endl;
//Shift and invert mode
eigenSolverStd(A, rEigValue, rEigVector, nev, _shiftInv, sigma);
out << rEigValue << std::endl;
out << rEigVector << std::endl;
\end{lstlisting}
\vspace{.1cm}
Real non-symmetric case: Matrix $A$ is real and non-symmetric\\
Eigenvalue and eigenvector can be complex, that's why we must use complex vector to get the results.
\vspace{.1cm}
\begin{lstlisting}[]{}
//Regular mode
eigenSolverStd(A, cEigValue, cEigVector, nev, _regular, sigma);
out << cEigValue << std::endl;
out << cEigVector << std::endl;
//Shift and invert mode
eigenSolverStd(A, cEigValue, cEigVector, nev, _shiftInv, sigma);
out << cEigValue << std::endl;
out <<cEigVector << std::endl;
\end{lstlisting}
\vspace{.1cm}
Real non-symmetric case: Matrix $A$ is complex\\
Eigenvalue and eigenvector are always complex, that's why we must use complex vector to get the results.
\vspace{.1cm}
\begin{lstlisting}[]{}
//Regular mode
eigenSolverStd(A, cEigValue, cEigVector, nev, _regular, sigma);
out << cEigValue << std::endl;
out << cEigVector << std::endl;
//Shift and invert mode
eigenSolverStd(A, cEigValue, cEigVector, nev, _shiftInv, sigma);
out << cEigValue << std::endl;
out << cEigVector << std::endl;
\end{lstlisting}
\vspace{.1cm}

\paragraph*{The generalized eigen problem $AX=\lambda BX$}

Real symmetric case: Matrix $A$ is real and symmetric; matrix $B$ is real positive definite.\\
Eigenvalues and eigenvectors are always real
\vspace{.1cm}
\begin{lstlisting}[]{}
//Regular mode
eigenSolverGen(A, B, rEigValue, rEigVector, nev, _regular, sigma);
out << rEigValue << std::endl;
out << rEigVector << std::endl;
//Shift and invert mode
eigenSolverGen(A, B, rEigValue, rEigVector, nev, _shiftInv, sigma);
out << rEigValue << std::endl;
out << rEigVector << std::endl;
//Buckling mode
eigenSolverGen(A, B, rEigValue, rEigVector, nev, _buckling, sigma);
out << rEigValue << std::endl;
out << rEigVector << std::endl;
//Cayley mode
eigenSolverGen(A, B, rEigValue, rEigVector, nev, _cayley, sigma);
out << rEigValue << std::endl;
out << rEigVector << std::endl;
\end{lstlisting}
\vspace{.1cm}
Real non-symmetric case: Matrix $A$ is real and non-symmetric; matrix $B$ is real positive definite.\\
Eigenvalue and eigenvector can be complex, that's why we must use complex vector to get the results.
\vspace{.1cm}
\begin{lstlisting}[]{}
//Regular mode
eigenSolverGen(A, B, cEigValue, cEigVector, nev, _regular, sigma);
out << cEigValue << std::endl;
out << cEigVector << std::endl;
//Shift and invert mode using the real part of shift value
eigenSolverGen(A, B, cEigValue, cEigVector, nev, _cplxShiftInv, 'R', Complex(-700.0, 10));
out << cEigValue << std::endl;
out << cEigVector << std::endl;
//Shift and invert mode using the imaginary part of shift value
eigenSolverGen(A, B, cEigValue, cEigVector, nev, _cplxShiftInv, 'I', Complex(-700.0, 10));
out << cEigValue << std::endl;
out << cEigVector << std::endl;
\end{lstlisting}
\vspace{.1cm}
Real non-symmetric case: Matrix $A$ is complex; matrix $B$ is complex positive definite.\\
Eigenvalue and eigenvector are always complex, that's why we must use complex vector to get the results.
\vspace{.1cm}
\begin{lstlisting}[]{}
//Regular mode
eigenSolverGen(A, B, cEigValue, cEigVector, nev, _regular, csigma);
out << cEigValue << std::endl;
out << cEigVector << std::endl;
//Shift and invert mode
eigenSolverGen(A, B, cEigValue, cEigVector, nev, _shiftInv, csigma);
out << cEigValue << std::endl;
out << cEigVector << std::endl;
\end{lstlisting}
\vspace{.1cm}

\subsubsection{Advanced Usage}

Besides using the available eigensolver classes: {\class EigenSolverStdArpackpp} and {\class EigenSolverGenArpackpp}, users can also take advantage of other functions of {\lib Arpack++} by using {\lib Arpack++} object directly.

\paragraph*{Interface classes between {\lib Xlife++} and {\lib Arpack++}}
To use of {\lib Arpack++} directly from {\lib Xlife++}, users must invoke one of these following class, which play the role of interfaces between these two libraries.
 \begin{itemize}
\item {\class LargeMatrixArpackppStdReg} for standard eigen problem in regular mode
\item {\class LargeMatrixArpackppStdShf} for standard eigen problem in shift and invert mode
\item {\class LargeMatrixArpackppGenReg} for generalized eigen problem in regular mode
\item {\class LargeMatrixArpackppGenShf} for generalized eigen problem in shift and invert, buckling and cayley mode
\end{itemize}
Given an {\lib Arpack++} class and a computational mode, users must
\begin{itemize}
\item create a {\class LargeMatrixArpackppXXX} object
\item create an {\lib Arpack++} object defining an eigen problem and using the previous object to specify the matrix-vector needed.
\end{itemize}
Because {\class LargeMatrixArpackppXXX} is templated, users can instantiate it with their own class that specifies matrix-vector product. However, we highly recommend instantiate this class by taking advantage of the available class {\class LargeMatrix} of {\lib Xlife++}!

\paragraph*{Using {\lib Arpack++} object}

In order to use {\lib Arpack++} object, depending on specific problem, users must include at least one of the {\lib Arpack++} headers file
\begin{lstlisting}[]{}
#include "arssym.h"     // for ARSymStdEig
#include "argsym.h"     // for ARSymGenEig

#include "arsnsym.h"    // for ARNonSymStdEig
#include "argnsym.h"    // for ARNonSymGenEig

#include "arscomp.h"    // for ARCompStdEig
#include "argcomp.h"    // for ARCompGenEig
\end{lstlisting}
After that, they need to define interface object, corresponding to eigen problem and calculation mode then {\lib Arpack++} object. The following example describes how to compute eigenvalue, eigenvector with {\lib Arpack++} object.\\
Supposed we'd like to solve a real symmetric standard problem in regular mode (corresponding {\lib Arpack++} class {\class ARSymStdEig}. Real symmetric matrix $rMatSkylineSymSym$ describes the problem is a {\class LargeMatrix}, the class interface is {\class LargeMatrixArpackppStdReg}
\vspace{.1cm}
\begin{lstlisting}[]{}
  //1. Create interface object
  LargeMatrixArpackppStdReg<Real, LargeMatrix<Real> > intMat(rMatSkylineSymSym, size);

  //2. Create Arpack++ problem object
  ARSymStdEig<Real, LargeMatrixArpackppStdReg<Real, LargeMatrix<Real> > >  arProb(size, nev, &intMat, &LargeMatrixArpackppStdReg<Real, LargeMatrix<Real> >::multOPx, (char *)"LM");

  //3. Compute eigenvalues
  Int nconv = arProb.FindEigenvalues();

  for (Int i = 0; i < nconv; ++i) {
      out << arProb.Eigenvalue(i) << "  ";
  }
  out << "\n";
  //4. Change some parameters
  arProb.ChangeNev(nev-2);
  arProb.ChangeTol(1.e-6);

  nconv = arProb.FindEigenvectors();

  for (Int i = 0; i < nconv; ++i) {
    out << arProb.Eigenvalue(i) << "  ";
  }
  out << "\n";

  //4. Retrieve eigenVector
  Real* eigenVec_p = 0;
  for (Int i = 0; i < nconv; ++i) {
      eigenVec_p = arProb.RawEigenvector(i);
      out << "Eigenvector [" << i << "] :" ;
      for (Int j = 0; j < size; j++) {
          out <<   *eigenVec_p << "  ";
          eigenVec_p++;
      }
      out << "\n";
  }

  out << "\n";
\end{lstlisting}
\vspace{.1cm}