1. Arpack++ classes
2. Interface classes between Melina++ and Arpack++
3. Link between classes
4. Usage
5. Warning

1. Test program
2. Guidelines to create classes

Standard eigenvalue problem of prototype Mass u = l u
Generalized eigenvalue problem of prototype Stiff u = l Mass u

3. Post-processing

The class TermMatrixAR is an interface class between Melina++ and Arpack++ designed to enable the use of Arpack++'s own classes to solve eigenvalue problems in the context of a Melina++ program. This class is derived from Melina++ class TermMatrix and holds informations related to matrices.
It is used along with its companion class TermVectorAR.

This class is an abstract base class so that the classes to be used in the program are the classes that are derived from this one (see 2. below).

1. Arpack++ classes

The Arpack++ classes intended to be used are those that require user-defined matrix-vector products (see documentation file arpack++/doc/arpackpp.doc.gz in the archive available at http://www.ime.unicamp.br/~chico/arpack++/).
They are : ARSymStdEig, ARNonSymStdEig, ARCompStdEig, ARSymGenEig, ARNonSymGenEig and ARCompGenEig.

Thus, several public functions MultXXX are provided to compute the different matrix-vector products needed by the Arpack++ classes. For that, some Melina++ own functions are used.

2. Interface classes between Melina++ and Arpack++

Several user interface classes have been created, allowing to use Arpack++ classes in a rather straightforward way. They are :
- TermMatrixStdReg and TermMatrixStdShf for standard eigenvalue problems,
- TermMatrixSymGenReg, TermMatrixGenReg, TermMatrixGenShf, TermMatrixGenBuc, TermMatrixGenCay and TermMatrixGenCsh for generalized eigenvalue problems.

The hierarchy of the classes is the following:

 (TermMatrixAR) > 1. TermMatrixStdReg
                > 2. TermMatrixStdShf
                >   (TermMatrixGen)  > 1. TermMatrixSymGenReg
                                     > 2. TermMatrixGenReg
                                     > 3. TermMatrixGenShf  > 4. TermMatrixGenBuc
                                                            > 5. TermMatrixGenCay
                                                            > 6. TermMatrixGenCsh

TermMatrixAR and TermMatrixGen are abstract classes.

Given an Arpack++ class and a computational mode, the user has to:
a. create a TermMatrixXXX object defining the matrices involved,
b. create an Arpack++ object defining the problem and using the previous object to specify the matrix-vector products needed.
The name of the matrix-vector product functions are the generic names used in the Arpack++ documentation so that the user's task reduces nearly to a copy-paste operation. Those names are MultOPx, MultAx and MultBx.

3. Link between classes

Here follows the list of the Arpack++ classes along with the name of the matrix-vector product functions and the interface constructor to be used according to the nature of the problem. The arguments of the interface constructors are given under each of them in a short way ; some of them have two forms.

Standard eigenvalue problem Au = l u

===================================================================================
 Kind of problem  Computational  Matrix-vector        Name of member    Interface
 (Arpack++ class)      mode      products required    functions to use  constructor
===================================================================================
 real symmetric  |
 (ARSymStdEig)   | Regular       y <- A * x             MultOPx      TermMatrixStdReg
                 |                                                   (1. A_u)
 real non        |-----------------------------------------------------------------
 symmetric       |
 (ARNonSymStdEig)| Shift
                 | and invert    y <- inv(A-sigma I)*x  MultOPx      TermMatrixStdShf
 complex         |                                                   (2. A_u_s, u_AsI)
 (ARCompStdEig)  |
===================================================================================

Generalized eigenvalue problem Au = l Bu
The matrix A defines the kind of problem.
If A is real, the matrix B is required to be real symmetric positive semi-definite, except in regular mode where it should be positive definite. In bukling mode, the real symmetric matrix A is required to be positive semi-definite while B is only required to be real symmetric indefinite.
If A is complex, the matrix B is required to be hermitian positive semi-definite, except in regular mode where it should be positive definite. Notice that B may still be real and symmetric.

===================================================================================
 Kind of problem  Computational  Matrix-vector        Name of member    Interface
 (Arpack++ class)      mode      products required    functions to use  constructor
===================================================================================
                 | Regular       y <- inv(B)*A * x ;    MultOPx      TermMatrixSymGenReg
                 |                      x <- A * x                   (1. A_B_u)
                 |               y <- B * x             MultBx
                 |-----------------------------------------------------------------
                 | Shift         y <- inv(A-sigma B)*x  MultOPx      TermMatrixGenShf
 real symmetric  | and invert    y <- B * x             MultBx       (3. A_B_u_s, B_u_AsB)
 (ARSymGenEig)   |-----------------------------------------------------------------
                 | Buckling      y <- inv(A-sigma B)*x  MultOPx      TermMatrixGenBuc     
                 |               y <- A * x             MultBx       (4. A_B_u_s)
                 |-----------------------------------------------------------------
                 |               y <- inv(A-sigma B)*x  MultOPx      TermMatrixGenCay
                 | Cayley        y <- A * x             MultAx       (5. A_B_u_s)
                 |               y <- B * x             MultBx
===================================================================================
                 | Regular       y <- inv(B)*A * x      MultOPx      TermMatrixGenReg
                 |               y <- B * x             MultBx       (2. A_B_u)
                 |-----------------------------------------------------------------
 real non        | Real shift    y <- inv(A-sigma B)*x  MultOPx      TermMatrixGenShf
 symmetric       | and invert    y <- B * x             MultBx       (3. A_B_u_s, B_u_AsB)
 (ARNonSymGenEig)|-----------------------------------------------------------------
                 | Complex shift y <- Inv(A-sigma B)*x  MultOPx      TermMatrixGenCsh
                 | and invert    y <- A * x             MultAx       (6. A_B_u_p_sr_si)
                 |               y <- B * x             MultBx
                 | Nota: Inv(A-sigma B) is real{inv(A-sigma B)} or imag{inv(A-sigma B)}
===================================================================================
                 | Regular       y <- inv(B)*A * x      MultOPx      TermMatrixGenReg
                 |               y <- B * x             MultBx       (2. A_B_u)
 complex         |-----------------------------------------------------------------
 (ARCompGenEig)  | Shift         y <- inv(A-sigma B)*x  MultOPx      TermMatrixGenShf
                 | and invert    y <- B * x             MultBx       (3. A_B_u_s, B_u_AsB)
===================================================================================

4. Usage

In the user program, one should include the file TermMatrixAR.h++, which is done by the mean of another header file present in Melina++ distribution:

using namespace std; // mandatory here for the moment...
                     // (modification of Melina++ header files still to do)
#include "Term_s.h++"

and additionally, at least one of the Arpack++ header files:

#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

For example, if we want to solve a real symmetric generalized problem in buckling mode (Arpack++ class ARSymGenEig), we will have to specify two matrix-vector products. The first product must compute inv(A-sigma B)*x and corresponds to the function called MultOPx. The second product must compute A * x and corresponds to the function called MultBx (this is not an error, in all other cases MultBx is supposed to compute B * x). In the documentation, the description of the constructor in shift and invert and buckling modes is:

   ARSymGenEig(char InvertMode, int n, int nev,
               FOP* objOP, TypeOPx MultOPx,
               FB*  objB , TypeBx MultBx  , FLOAT sigma, ... + optional parameters ...)

Thus, in our program, once matrices A, B and vector u have been defined (see Remark 1), and assuming A is real symmetric positive semi-definite and B is symmetric, we'll have to write something like the following:

//   Define interface objet, choosing Buckling mode:
   real_t sigma = ...;
   TermMatrixGenBuc IntMat(A,B,u,sigma);

//   Define Arpack++ problem:
   int n = IntMat.GetN(); // Dimension of the problem
   int nev = ...;         // Number of eigenvalues to be computed
   ARSymGenEig< real_t, TermMatrixGenBuc, TermMatrixGenBuc >
    ArProb('B', n, nev, &IntMat, &TermMatrixGenBuc::MultOPx,
                        &IntMat, &TermMatrixGenBuc::MultBx, sigma);

//    Compute eigenvalues and eigenvectors
    int nconv = ArProb.FindEigenvectors();

//    Extract ith eigenvalue and eigenvector, 0 <= i < nconv
    for (int i = 0; i < nconv; i++) {
      real_t lambda = ArProb.Eigenvalue(i);
      TermVectorAR v (u,ArProb.RawEigenvector(i));
      // use lambda and v ...
    }

Remark 1:
The vector u is a vector associated to the problem. It is typically the right-hand side vector but it may be any vector term compatible with the problem and defined on the same domain. For the computation of the eigenvalues, it is only used to carry the numbering structure of the unknowns.

Remark 2:
The member function GetN() is better to be used to determine the dimension of the problem. It returns N*dim, where N is the number of unknowns and dim the number of components of each unknown.

Remark 3:
With the original version of Arpack++, one cannot include Arpack++ header files in more than one single translation unit because at link time, duplicated symbols would be detected. Since this may be a problem, our own installation procedure of the Arpack++ library takes care of this issue by adding an implementation file that should be specified at link time (the script mmake is provided for that).

5. Warning

The interface matrix and the Arpack++ problem are two objects which are tightly linked and thus should remain coherent.

Once the Arpack++ problem is fully defined, the parameters describing the problem should not be modified. Thus, in particular, the Arpack++ functions ChangeShift and NoShift MUST NOT BE USED because the matrix-vector product (defined in the interface object) would not reflect the change of the shift. For the same reason, the SetXXXMode functions, which are intended to define the problem starting from an object built by the default constructor, should not be used in this context. The functions ChangeMultBx and ChangeMultOPx should not be used either since the matrix-vector product functions are provided by the interface object.

If the shift has to be changed, then the user MUST create a new interface object (TermMatrixXXX object) and use it to define a new Arpack++ problem.

Note however that the SetXXXMode functions can be used in the definition phase of the problem. After the Arpack++ problem object has been defined, one can modify all the other non critical parameters by the mean of the functions ChangeMaxit, ChangeNcv, ChangeNev, ChangeTol and ChangeWhich.

1. Test program

Consider the following test program, given only to show the classes in a plausible environment:

#include "MELINA.h++"   // this does #include "Term_s.h++"
#include "arssym.h"     // for ARSymStdEig
int main ()
{
   init(); // Melina++ initialization  

   Mesh mesh = Square(1,quadrangle,0.,1.);  // 1 element in the unit square
   GeomDomain& Omega(mesh.domain("Omega")); // called Omega

   // Define Finite Element Space over Omega
   FESpace V_h (Omega,Lagrange,P1,"H_1(Omega)");
   // Define Unknown
   FEUnknown u_h ("u_h",V_h);

   // Define Terms to be computed on Mesh's GeomDomains
   FETermMatrix Mass  (Omega, u_h, u_h);
   NVTermVector uzero (Omega, u_h);

   Omega.computation(0);

/*
    Standard eigenvalue problem Mass * u = l * u  solved in regular mode
    --------------------------------------------------------------------
*/
   int nev = 2;   // Number of eigenvalues to be computed

//  1. create interface object
   TermMatrixStdReg< real_t > IntMat(Mass, uzero);

//  2. create Arpack++ problem object
   ARSymStdEig< real_t, TermMatrixStdReg< real_t > >
    ArProb(IntMat.GetN(), nev, &IntMat, &TermMatrixStdReg< real_t >::MultOPx, (char *)"LM");

//  3. compute eigenvalues
   int nconv = ArProb.FindEigenvalues();

   std::cout << std::endl;
   for (int i = 0; i < nconv; i++) {
      std::cout << ArProb.Eigenvalue(i) << std::endl;
   }
}

This program follows the specifications given in the previous section, in the case of a standard eigenvalue problem. It prints:

0.0833333
0.25

A matrix Mass and a vector uzero related to a finite element space V_h are defined. Their values are computed by the statement Omega.computation(0);.
The vector uzero is defined only to carry the finite element structure enabling correct evaluations of matrix-vector products inside the interface class TermMatrixStdReg. Its value is not significant.

2. Guidelines to create classes

In the following, we just give the statements needed to create the interface and Arpack++ objects in all cases, assuming the matrices have the correct type and properties.
The standard eigenvalue problem is always Mass * u = l * u.
The generalized eigenvalue problem is always Stiff * u = l * Mass * u, provided a finite element matrix Stiff has been defined, for example by writing:

   FETermMatrix Stiff (Omega, grad(u_h), grad(u_h));

The tool vector uzero is assumed to be created as shown above. The user just has to adapt his program by replacing these names by his own ones.

In all the cases, the requested eigenvalues are those of Largest Magnitude (argument (char *)"LM"). The interface object is called IntMat and the Arpack++ object is called ArProb.

Standard eigenvalue problem of prototype Mass u = l u

Real symmetric case

{// Regular mode
   TermMatrixStdReg< real_t >  IntMat(Mass, uzero);
   ARSymStdEig< real_t, TermMatrixStdReg< real_t > >
    ArProb(IntMat.GetN(), nev, &IntMat, &TermMatrixStdReg< real_t >::MultOPx, (char *)"LM");
}
{// Shift and invert mode. 1st constructor
   real_t sigma = 0.2; // for example
   TermMatrixStdShf< real_t >  IntMat(Mass, uzero, sigma);
   ARSymStdEig< real_t, TermMatrixStdShf< real_t > >
    ArProb(IntMat.GetN(), nev, &IntMat, &TermMatrixStdShf< real_t >::MultOPx, sigma, (char *)"LM");
}
{// Shift and invert mode. 2nd constructor
   real_t sigma = 0.2; // for example
   TermMatrix MsI(Mass-sigma*identity(Mass));
   TermMatrixStdShf< real_t >  IntMat(uzero, MsI);
   ARSymStdEig< real_t, TermMatrixStdShf< real_t > >
    ArProb(IntMat.GetN(), nev, &IntMat, &TermMatrixStdShf< real_t >::MultOPx, sigma, (char *)"LM");
}

Real non symmetric case

{// Regular mode
   TermMatrixStdReg< real_t >  IntMat(Mass, uzero);
   ARNonSymStdEig< real_t, TermMatrixStdReg< real_t > >
    ArProb(IntMat.GetN(), nev, &IntMat, &TermMatrixStdReg< real_t >::MultOPx, (char *)"LM");
}
{// Shift and invert mode. 1st constructor
   real_t sigma = 0.2; // for example
   TermMatrixStdShf< real_t >  IntMat(Mass, uzero, sigma);
   ARNonSymStdEig< real_t, TermMatrixStdShf< real_t > >
    ArProb(IntMat.GetN(), nev, &IntMat, &TermMatrixStdShf< real_t >::MultOPx, sigma, (char *)"LM");
}
{// Shift and invert mode. 2nd constructor
   real_t sigma = 0.2; // for example
   TermMatrix MsI(Mass-sigma*identity(Mass));
   TermMatrixStdShf< real_t >  IntMat(uzero, MsI);
   ARNonSymStdEig< real_t, TermMatrixStdShf< real_t > >
    ArProb(IntMat.GetN(), nev, &IntMat, &TermMatrixStdShf< real_t >::MultOPx, sigma, (char *)"LM");
}

Complex case

{// Regular mode
   TermMatrixStdReg< complex_t >  IntMat(Mass, uzero);
   ARCompStdEig< real_t, TermMatrixStdReg< complex_t > >
    ArProb(IntMat.GetN(), nev, &IntMat, &TermMatrixStdReg< complex_t >::MultOPx, (char *)"LM");
}
{// Shift and invert mode. 1st constructor
   complex_t sigmac = 0.2; // for example
   TermMatrixStdShf< complex_t >  IntMat(Mass, uzero, sigmac);
   ARCompStdEig< real_t, TermMatrixStdShf< complex_t > >
    ArProb(IntMat.GetN(), nev, &IntMat, &TermMatrixStdShf< complex_t >::MultOPx, sigmac, (char *)"LM");
}
{// Shift and invert mode. 2nd constructor
   complex_t sigmac = 0.2; // for example
   TermMatrix MsI(Mass-sigmac*identity(Mass));
   TermMatrixStdShf< complex_t >  IntMat(uzero, MsI);
   ARCompStdEig< real_t, TermMatrixStdShf< complex_t > >
    ArProb(IntMat.GetN(), nev, &IntMat, &TermMatrixStdShf< complex_t >::MultOPx, sigmac, (char *)"LM");
}

Generalized eigenvalue problem of prototype Stiff u = l Mass u

Real symmetric case

{// Regular mode
   TermMatrixSymGenReg  IntMat(Stiff, Mass, uzero);
   ARSymGenEig< real_t, TermMatrixSymGenReg, TermMatrixSymGenReg >
    ArProb(IntMat.GetN(), nev, &IntMat, &TermMatrixSymGenReg::MultOPx,
                               &IntMat, &TermMatrixSymGenReg::MultBx,    (char *)"LM");
}
{// Shift and invert mode. 1st constructor
   real_t sigma = 2.2; // for example
   TermMatrixGenShf< real_t >  IntMat(Stiff, Mass, uzero, sigma);
   ARSymGenEig< real_t, TermMatrixGenShf< real_t >, TermMatrixGenShf< real_t > >
    ArProb('S', IntMat.GetN(), nev, &IntMat, &TermMatrixGenShf< real_t >::MultOPx,
                                    &IntMat, &TermMatrixGenShf< real_t >::MultBx, sigma, (char *)"LM");
}
{// Shift and invert mode. 2nd constructor
   real_t sigma = 2.2; // for example
   TermMatrix SsM(Stiff-sigma*Mass);
   TermMatrixGenShf< real_t >  IntMat(Mass, uzero, SsM);
   ARSymGenEig< real_t, TermMatrixGenShf< real_t >, TermMatrixGenShf< real_t > >
    ArProb('S', IntMat.GetN(), nev, &IntMat, &TermMatrixGenShf< real_t >::MultOPx,
                                    &IntMat, &TermMatrixGenShf< real_t >::MultBx, sigma, (char *)"LM");
}
{// Buckling mode.
   real_t sigma = 2.2; // for example
   TermMatrixGenBuc  IntMat(Stiff,Mass, uzero, sigma);
   ARSymGenEig< real_t, TermMatrixGenBuc, TermMatrixGenBuc >
    ArProb('B', IntMat.GetN(), nev, &IntMat, &TermMatrixGenBuc::MultOPx,
                                    &IntMat, &TermMatrixGenBuc::MultBx, sigma, (char *)"LM");
}
{// Cayley mode.
   real_t sigma = 2.2; // for example
   TermMatrixGenCay  IntMat(Stiff,Mass, uzero, sigma);
   ARSymGenEig< real_t, TermMatrixGenCay, TermMatrixGenCay >
    ArProb(IntMat.GetN(), nev, &IntMat, &TermMatrixGenCay::MultOPx,
                               &IntMat, &TermMatrixGenCay::MultAx,
                               &IntMat, &TermMatrixGenCay::MultBx, sigma, (char *)"LM");
}

Real non symmetric case

{// Regular mode
   TermMatrixGenReg< real_t >  IntMat(Stiff, Mass, uzero);
   ARNonSymGenEig< real_t, TermMatrixGenReg< real_t >, TermMatrixGenReg< real_t > >
    ArProb(IntMat.GetN(), nev, &IntMat, &TermMatrixGenReg< real_t >::MultOPx,
                               &IntMat, &TermMatrixGenReg< real_t >::MultBx,    (char *)"LM");
}
{// Real shift and invert mode. 1st constructor
   real_t sigma = 2.2; // for example
   TermMatrixGenShf< real_t >  IntMat(Stiff, Mass, uzero, sigma);
   ARNonSymGenEig< real_t, TermMatrixGenShf< real_t >, TermMatrixGenShf< real_t > >
    ArProb(IntMat.GetN(), nev, &IntMat, &TermMatrixGenShf< real_t >::MultOPx,
                               &IntMat, &TermMatrixGenShf< real_t >::MultBx, sigma, (char *)"LM");
}
{// Real shift and invert mode. 2nd constructor
   real_t sigma = 2.2; // for example
   TermMatrix SsM(Stiff-sigma*Mass);
   TermMatrixGenShf< real_t >  IntMat(Mass, uzero, SsM);
   ARNonSymGenEig< real_t, TermMatrixGenShf< real_t >, TermMatrixGenShf< real_t > >
    ArProb(IntMat.GetN(), nev, &IntMat, &TermMatrixGenShf< real_t >::MultOPx,
                               &IntMat, &TermMatrixGenShf< real_t >::MultBx, sigma, (char *)"LM");
}
{// Complex shift and invert mode. Using real part
   real_t sigma  = 2.2; // for example
   real_t sigmaI = 0.;
   TermMatrixGenCsh  IntMat(Stiff, Mass, uzero, 'R', sigma, sigmaI);
   ARNonSymGenEig< real_t, TermMatrixGenCsh, TermMatrixGenCsh >
    ArProb(IntMat.GetN(), nev, &IntMat, &TermMatrixGenCsh::MultOPx,
                               &IntMat, &TermMatrixGenCsh::MultAx,
                               &IntMat, &TermMatrixGenCsh::MultBx, 'R', sigma, sigmaI, (char *)"LM");
}
{// Complex shift and invert mode. Using imaginary part
   real_t sigma  = 2.2; // for example
   real_t sigmaI = 0.5; // should not be 0 here
   TermMatrixGenCsh  IntMat(Stiff, Mass, uzero, 'I', sigma, sigmaI);
   ARNonSymGenEig< real_t, TermMatrixGenCsh, TermMatrixGenCsh >
    ArProb(IntMat.GetN(), nev, &IntMat, &TermMatrixGenCsh::MultOPx,
                               &IntMat, &TermMatrixGenCsh::MultAx,
                               &IntMat, &TermMatrixGenCsh::MultBx, 'I', sigma, sigmaI, (char *)"LM");
}

Complex case

{// Regular mode
   TermMatrixGenReg< complex_t >  IntMat(Stiff, Mass, uzero);
   ARCompGenEig< real_t, TermMatrixGenReg< complex_t >, TermMatrixGenReg< complex_t > >
    ArProb(IntMat.GetN(), nev, &IntMat, &TermMatrixGenReg< complex_t >::MultOPx,
                               &IntMat, &TermMatrixGenReg< complex_t >::MultBx,    (char *)"LM");
}
{// Shift and invert mode. 1st constructor
   complex_t sigmac = 0.2; // for example
   TermMatrixGenShf< complex_t >  IntMat(Stiff, Mass, uzero, sigmac);
   ARCompGenEig< real_t, TermMatrixGenShf< complex_t >, TermMatrixGenShf< complex_t > >
    ArProb(IntMat.GetN(), nev, &IntMat, &TermMatrixGenShf< complex_t >::MultOPx,
                               &IntMat, &TermMatrixGenShf< complex_t >::MultBx, sigmac, (char *)"LM");
}
{// Shift and invert mode. 2nd constructor
   complex_t sigmac = 0.2; // for example
   TermMatrix SsM(Stiff-sigmac*Mass);
   TermMatrixGenShf< complex_t >  IntMat(Mass, uzero, SsM);
   ARCompGenEig< real_t, TermMatrixGenShf< complex_t >, TermMatrixGenShf< complex_t > >
    ArProb(IntMat.GetN(), nev, &IntMat, &TermMatrixGenShf< complex_t >::MultOPx,
                               &IntMat, &TermMatrixGenShf< complex_t >::MultBx, sigmac, (char *)"LM");
}

3. Post-processing

Once the Arpack++ problem has been defined, the computation is ready to be carried on, using member functions of the ArProb object.

In order to make the computation converge, one may tune beforehand the computational process, which is done by altering the default parameters of the Arpack++ problem. Several parameters can be changed as explained in the Arpack++ documentation (see paragraph "Functions that allow changes in problem parameters"). See also the warning written in the previous section Warning.

• If the eigenvalues only are wanted, then the following code can typically be used.

//    Eventually, change critical parameters to allow convergence
   ArProb.ChangeNcv(min(2*nev+1,IntMat.GetN()-1));
   ArProb.ChangeTol(1.e-6);

//    Compute eigenvalues
   int nconv = ArProb.FindEigenvalues();

//    Extract ith eigenvalue, 0 <= i < nconv
   for (int i = 0; i < nconv; i++) {
     real_t lambda = ArProb.Eigenvalue(i);
     // use lambda ...
   }

• If the eigenvectors are also wanted, then the following code can typically be used:

//    Compute eigenvalues and eigenvectors
   int nconv = ArProb.FindEigenvectors();

//    Extract ith eigenvalue and eigenvector, 0 <= i < nconv
   for (int i = 0; i < nconv; i++) {
     real_t lambda = ArProb.Eigenvalue(i);
     TermVectorAR v (uzero,ArProb.RawEigenvector(i));
     // use lambda and v ...
   }

In this case, one may want to write the eigenvectors for a graphical post-processing using medit. This can be achieved in the following way:

//    Compute eigenvalues and eigenvectors
   int nconv = ArProb.FindEigenvectors();

//    Extract ith eigenvalue and eigenvector, 0 <= i < nconv
   for (int i = 0; i < nconv; i++) {
     real_t lambda = ArProb.Eigenvalue(i);
     stringstream ss;
     ss << i+1 << "_" << lambda;
     string OutputFile ("EigVec" + ss.str()); // Root filename is 'EigVecN_lambda'
     TermVectorAR v (uzero,ArProb.RawEigenvector(i));
     meditFormat(OutputFile, v, u_h, Omega);  // Creation of 2 files : 'EigVecN_lambda.mesh' (geometry)
                                              //                   and 'EigVecN_lambda.bb'   (values)
   }