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);

};

}

enum errorinfo {ok = 0, no_eigenvalue, not_calculated};

The

`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

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

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 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.

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.