Lanczos.h

Go to the documentation of this file.

#ifndef REDUKT_LANCZOS_H
#define REDUKT_LANCZOS_H


#include <ietl/interface/ublas.h>
#include <ietl/vectorspace.h>
#include <ietl/lanczos.h>
#include <ietl/iteration.h>

#include <boost/random.hpp>
#include <boost/limits.hpp>

#include <vector>
#include <algorithm>  // for swap()
#include <limits>
#include <cmath>

#include <iostream>
#include <iomanip>

#include <assert.h>


#include <boost/numeric/ublas/io.hpp>


template<class Matrix,class Vector>
class Lanczos
{
public:
        typedef ietl::vectorspace<Vector>  Vectorspace;
        typedef boost::lagged_fibonacci607 RandGen;

        Lanczos( const Matrix& mat, int k )
        {
                assert( mat.size1() == mat.size2() );
                int N = mat.size1();
                
                //
                // Compute eigenvalues
                //
                        
                Vectorspace vec(N);
                RandGen  gen;   
                ietl::lanczos<Matrix,Vectorspace> lanczos( mat, vec );
                
                double rel_tol = 500*std::numeric_limits<double>::epsilon();
                double abs_tol = std::pow(std::numeric_limits<double>::epsilon(),2./3); 
                int max_iter = 2*N;
                
                ietl::lanczos_iteration_nhighest<double> iter( max_iter, k, rel_tol, abs_tol );
                
                std::cout << "Starting Lanczos calculation of " << k << " highest eigenvalues with " << max_iter << " max.iterations" << std::endl;
                try 
                {
                        lanczos.calculate_eigenvalues( iter, gen );
                        eigen = lanczos.eigenvalues();
                        multi = lanczos.multiplicities();
                        err   = lanczos.errors();
                }
                catch( std::runtime_error& e ) 
                {
                        std::cout << e.what() << std::endl;
                        
                        // Printing eigenvalues with error & multiplicities:  
                        std::cout << "#        eigenvalue            error         multiplicity\n";  
                        std::cout.precision(10);
                        for( unsigned int i=0; i < eigen.size(); ++i ) 
                                std::cout << i << "\t" << eigen[i] << "\t" << err[i] << "\t" 
                                                  << multi[i] << "\n";
                        
                        throw std::runtime_error("Problems with ietl lanczos algorithm occured during eigenvalue computation");
                }
                
                std::cout << "Lanczos eigenvalue calculation finished" << std::endl;

                // computed eigenvalues are stored ascending
                std::cout << "Lanczos reordering eigenvalues descending " << std::flush;
                for( unsigned int i=0; i < eigen.size()/2; ++i )
                        std::swap( eigen[i], eigen[eigen.size()-i-1] );
                std::cout << "done" << std::endl;
                
                
                //
                // Compute eigenvectors
                //
                
                std::cout << "Starting Lanczos calculation of eigenvectors" << std::endl;
                try
                {
                        lanczos.eigenvectors( eigen.begin(), eigen.begin() + k,
                                                                  std::back_inserter(eigenvectors),
                                                                  info,
                                                                  gen );
                }
                catch( std::runtime_error& e )
                {
                        std::cout << e.what() << std::endl;                     
                        throw std::runtime_error("Problems with ietl lanczos algorithm occured during eigenvector computation");
                }
                std::cout << "Lanczos eigenvector calculation finished" << std::endl;
        }
        
        void printInfo()
        {
                // Printing eigenvalues with error & multiplicities:  
                std::cout << "#        eigenvalue            error         multiplicity\n";  
                std::cout.precision(10);
                for( unsigned int i=0; i < eigen.size(); ++i ) 
                        std::cout << i << "\t" << eigen[i] << "\t" << err[i] << "\t" 
                                          << multi[i] << "\n";
                
                // Print eigenvectors
                for( unsigned int i=0; i < eigenvectors.size(); ++i )
                        std::cout << "Eigenvector " << i << " = " << eigenvectors[i] << std::endl;
                
                std::cout << std::endl;
                std::cout << " Information about the eigen vectors computations:\n\n";
                for(int i = 0; i < info.size(); i++) 
                {
                        std::cout << " m1(" << i+1 << "): " << info.m1(i) << ", m2(" << i+1 << "): "
                                  << info.m2(i) << ", ma(" << i+1 << "): " << info.ma(i) << " eigenvalue("
                                  << i+1 << "): " << info.eigenvalue(i) << " residual(" << i+1 << "): "
                                  << info.residual(i) << " error_info(" << i+1 << "): "
                                  << info.error_info(i) << std::endl << std::endl;
                }               
        }
        
        const Vector& getEigenvector( unsigned int i ) const
        {
                assert( i < eigenvectors.size() );
                return eigenvectors[i];
        }
        
        double getEigenvalue( unsigned int i ) const
        {
                assert( i < eigenvectors.size() );
                return eigen[i];                
        }

private:
        std::vector<double> eigen;  // eigenvalues
        std::vector<int>    multi;  // multiplicities
        std::vector<double> err;    // error

        std::vector<Vector> eigenvectors; // eigenvectors
        ietl::Info<double>  info;         // (m1, m2, ma, eigenvalue, residual, status).
};

#endif // REDUKT_LANCZOS_H