IETL Documentation - lanczos.h

# IETL library: lanczos.h

This header implements the Lanczos algorithm without reorthogonalization, as invented by Cullum and Willoughby [3,4]. This algorithm is quite complex and we recommend the user to study the book by Cullum and Willoughby  on whose Fortran implementation this generic C++ implementation is based.

## Documentation

### Synopsis

`namespace ietl {   template <class magnitude_type=double>  class Info {  public:    enum errorinfo {ok = 0, no_eigenvalue, not_calculated};              Info() {}            int m1(int i) const;    int m2(int i) const;    int ma(int i) const;    int size() const;    magnitude_type eigenvalue(int i) const;    magnitude_type residual(int i) const;    errorinfo error_info(int i) const;      };    template <class VS>  class Tmatrix {  public:    typedef typename vectorspace_traits<VS>::scalar_type scalar_type;     typedef typename vectorspace_traits<VS>::vector_type vector_type;    typedef typename vectorspace_traits<VS>::magnitude_type magnitude_type;    typedef typename vectorspace_traits<VS>::size_type size_type;        Tmatrix() {}            const std::vector<magnitude_type>& eigenvalues(bool discard_ghosts=true) const;    const std::vector<magnitude_type>& errors(bool discard_ghosts=true) const;     const std::vector<int>& multiplicities(bool discard_ghosts=true) const;  };  template <class MATRIX, class VS>  class lanczos : public Tmatrix<VS> {   public:    typedef typename vectorspace_traits<VS>::vector_type vector_type;    typedef typename vectorspace_traits<VS>::scalar_type scalar_type;    typedef typename vectorspace_traits<VS>::magnitude_type magnitude_type;        lanczos(const MATRIX& matrix, const VS& vec);        template <class IT, class GEN>    void calculate_eigenvalues(IT& iter, GEN gen);        template <class IT>    void more_eigenvalues(IT& iter);        template <class IN, class OUT, class GEN>    void eigenvectors(IN in_eigvals_start, IN in_eigvals_end , OUT eig_vectors,                       Info<magnitude_type>& info, GEN gen, int maxiter=0);      };}`

### The Lanczos information class Info

The Lanczos information class Info is used to return information on eigenvector calculations. It is templated on the magnitude_type of the vector space, which defaults to double.
`  enum errorinfo {ok = 0, no_eigenvalue, not_calculated};      `

The errorinfo type identifies error codes for eigenvector calculations:
• ok indicates that the eigenvector was calculated without any problems.
• no_eigenvalue indicates that the value for which an eigenvector was requested is not an eigenvalue of the matrix.
• not_calculated indicates that the eigenvector could not be calculated since the maximum number of iterations would be exceeded.
`  int m1(int i) const;`
is the T-matrix size at which the i-th eigenvalue has converged for the first time.
`  int m2(int i) const;`

is the T-matrix size at which the i-th eigenvalue has converged for the second time.

`  int ma(int i) const;`

is the T-matrix size used for the calculation of the i-th eigenvector.

`  int size() const;`

is the maximum size of T-matrix calculated by the Lanczos iterations.

`    magnitude_type eigenvalue(int i) const;`

is an improved estimate of the i-th eigenvalue

`    magnitude_type residual(int i) const;`

is the residual of the i-th eigenvector

`    errorinfo error_info(int i) const;  `
returns the error code for the calculation of the i-th eigenvector.

### The T-matrix class Tmatrix

The Lanczos algorithms converts a matrix to a tridiagonal matrix, called the T-matrix, which is stored in the Tmatrix class, templated on the magnitude_type of the vector space.
`    typedef typename vectorspace_traits<VS>::scalar_type scalar_type;     typedef typename vectorspace_traits<VS>::vector_type vector_type;    typedef typename vectorspace_traits<VS>::magnitude_type magnitude_type;    typedef typename vectorspace_traits<VS>::size_type size_type;`

are typedefs as shortcuts.

`    const std::vector<magnitude_type>& eigenvalues(bool discard_ghosts=true) const;    const std::vector<magnitude_type>& errors(bool discard_ghosts=true) const;     const std::vector<int>& multiplicities(bool discard_ghosts=true) const;`

return vectors containing the eigenvalues, errors and multiplicities of the T-matrix eigenvalues, sorted from lowest to highest eigenvalue. Instead of  expensive reorthogonalization the Cullum-Willoughby version of the Lanczos algorithm lives with roundoff errors. After the Lanczos iterations it detects incorrect eigenvalues due to roundoff errors, called ghosts. A true argument passed to any of the three functions removes all these spurious eigenvalues.

The errors of the eigenvalues are calculated for isolated single eigenvalues (multiplicity 1) and give an estimate for the error for all unconverged eigenvalues. Once an eigenvalue appears multiple times it is defined to be converged and the errors are set to 0.

The multiplicities of the eigenvalues are not the multiplicities of the original matrix but the multiplicities in the T-matrix, which can be larger because of spurious eigenvalues converging to multiple copies of the real one. Spurious eigenvalues (ghosts) are marked by zero multiplicity if not discarded.

### The Lanczos algorithm class lanczos

The Lanczos algorithm is more complex than the other iterative algorithms and requires two passes of the iterations to calculate the eigenvectors. In the first pass the T-matrix is constructed and eigenvalues of the original matrix are calculated from the T-matrix. The eigenvectors of the T-matrix are in a different basis than the original matrix and a second pass through the iterations is required to transform them back to the original basis. Since information (the T-matrix) needs to be kept from one pass to the other the Lanczos algorithm is implemented as a class, derived from the Tmatrix class. Thus, all member functions of Tmatrix are also accessible to the lanczos class.

The lanczos class is templated on the matrix and vector space types.
`  lanczos(const MATRIX& matrix, const VS& vec);`

The constructor takes a reference to the matrix and the vector space. Since only references are stored inside lanczos, the lifetime of these objects has to be at least the same as that of the lanczos object.

`  template <class IT, class GEN>  void calculate_eigenvalues(IT& iter, GEN gen);`

runs the Lanczos iteration, starting from a random vector generated from the generator gen and runs as long as the Lanczos iteration control object iter tells it to. Note that the Lanczos iteration control object needs to fulfill other concepts than the basic iteration control object since more than one eigenvalue is calculated. Memory requirements are for three vectors. After a call to calculate_eigenvalues the eigenvalues can be accessed through the member functions of the T-matrix base class.

`  template <class IT>  void more_eigenvalues(IT& iter);`

continues the Lanczos iterations from the last two Lanczos vectors stored in the lanczos class.

```  template <class IN, class OUT, class GEN>
void eigenvectors(IN in_eigvals_start, IN in_eigvals_end , OUT eig_vectors,                   Info<magnitude_type>& info, GEN gen, int maxiter=0);    ```

calculates eigenvectors of the eigenvalues passed by a range of iterators. The memory requirements are for 3+in_eigvals_end-in_eigvals_start vectors.

The eigenvector calculation proceeds by running the Lanczos iterations a second time to convert T-matrix eigenvalues back to the original basis. To initialize this second iteration, the same generator object gen used for the eigenvalue calculation needs to be passed.

The resulting eigenvectors are copied to the output iterator for vectors eig_vectors, and information about T-matrix sizes or convergence problems is passed back in the info object.

Eigenvector calculation may need a larger size T-matrix than eigenvalue. The implementation estimates the needed size of T-matrix and attempts eigenvector calculation. If no sufficient accuracy is obtained after increasing the T-matrix by about 50%, the algorithm stops it attempts, does not caculate the eigenvector and sets the error_info for that eigenvalue not not_calculated. If that happens, a larger number of iterations can be allowed by passing the allowable number of additional iterations in the optional maxiter argument.

copyright 2002-2004 by Matthias Troyer and Prakash Dayal