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 errorinfo type identifies error codes for eigenvector calculations:
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.
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.