lapack_sym_eigs.hpp

/*
 * Copyright (c) 2006-2014, Visualization and Multimedia Lab,
 *                          University of Zurich <http://vmml.ifi.uzh.ch>,
 *                          Eyescale Software GmbH,
 *                          Blue Brain Project, EPFL
 *
 * This file is part of VMMLib <https://github.com/VMML/vmmlib/>
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * Redistributions of source code must retain the above copyright notice, this
 * list of conditions and the following disclaimer.  Redistributions in binary
 * form must reproduce the above copyright notice, this list of conditions and
 * the following disclaimer in the documentation and/or other materials provided
 * with the distribution.  Neither the name of the Visualization and Multimedia
 * Lab, University of Zurich nor the names of its contributors may be used to
 * endorse or promote products derived from this software without specific prior
 * written permission.
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 */
#ifndef __VMML__VMMLIB_LAPACK_SYM_EIGS__HPP__
#define __VMML__VMMLIB_LAPACK_SYM_EIGS__HPP__

#include <vmmlib/matrix.hpp>
#include <vmmlib/vector.hpp>
#include <vmmlib/exception.hpp>

#include <vmmlib/lapack_types.hpp>
#include <vmmlib/lapack_includes.hpp>

#include <string>
#include <vector>

namespace vmml
{

namespace lapack
{

    // XYYZZZ
    // X    = data type: S - float, D - double
    // YY   = matrix type, GE - general, TR - triangular
    // ZZZ  = function name


    template< typename float_t >
    struct eigs_params
    {
        char            jobz;
        char            range;
        char            uplo;
        lapack_int      n;
        float_t*        a;
        lapack_int      lda; //leading dimension of input array
        float_t*        vl;
        float_t*        vu;
        lapack_int      il;
        lapack_int      iu;
        float_t         abstol;
        lapack_int      m; //number of found eigenvalues
        float_t*        w; //first m eigenvalues
        float_t*        z; //first m eigenvectors
        lapack_int      ldz; //leading dimension of z
        float_t*        work;
        lapack_int      lwork;
        lapack_int*     iwork;
        lapack_int*     ifail;
        lapack_int      info;

        friend std::ostream& operator << ( std::ostream& os,
                                          const eigs_params< float_t >& p )
        {
            os
            << " (1)\tjobz "     << p.jobz << std::endl
            << " (2)\trange "    << p.range << std::endl
            << " (3)\tuplo "     << p.uplo << std::endl
            << " (4)\tn "        << p.n << std::endl
            << " (5)\ta "        << *p.a << std::endl
            << " (6)\tlda "      << p.lda << std::endl
            << " (7)\tvl "       << p.vl << std::endl
            << " (8)\tvu "       << p.vu << std::endl
            << " (9)\til "       << p.il << std::endl
            << " (10)\tiu "       << p.iu << std::endl
            << " (11)\tabstol "   << p.abstol << std::endl
            << " (12)\tm "        << p.m << std::endl
            << " (13)\tw "        << p.w << std::endl
            << " (14)\tz "        << p.z << std::endl
            << " (15)\tldz "      << p.ldz << std::endl
            << " (16)\twork "     << *p.work << std::endl
            << " (17)\tlwork "    << p.lwork << std::endl
            << " (18)\tiwork "    << *p.iwork << std::endl
            << " (19)\tifail "    << *p.ifail << std::endl
            << " (20)\tinfo "     << p.info
            << std::endl;
            return os;
        }

    };


#if 0
    /* Subroutine */

    int dsyevx_( char *jobz, char *range, char *uplo, integer *n, doublereal *a, integer *lda, doublereal *vl,
                doublereal *vu, integer *il, integer *iu, doublereal *abstol, integer *m, doublereal *w, doublereal *z,
                integer *ldz, doublereal *work, integer* lwork, integer* iwork, integer* ifail, integer* info );

#endif


    template< typename float_t >
    inline void
    sym_eigs_call( eigs_params< float_t >& )
    {
        VMMLIB_ERROR( "not implemented for this type.", VMMLIB_HERE );
    }


    template<>
    inline void
    sym_eigs_call( eigs_params< float >& p )
    {
        //std::cout << "calling lapack sym x eigs (single precision) " << std::endl;
        ssyevx_(
                &p.jobz,
                &p.range,
                &p.uplo,
                &p.n,
                p.a,
                &p.lda,
                p.vl,
                p.vu,
                &p.il,
                &p.iu,
                &p.abstol,
                &p.m,
                p.w,
                p.z,
                &p.ldz,
                p.work,
                &p.lwork,
                p.iwork,
                p.ifail,
                &p.info
                );

    }


    template<>
    inline void
    sym_eigs_call( eigs_params< double >& p )
    {
        //std::cout << "calling lapack sym x eigs (double precision) " << std::endl;
        dsyevx_(
                &p.jobz,
                &p.range,
                &p.uplo,
                &p.n,
                p.a,
                &p.lda,
                p.vl,
                p.vu,
                &p.il,
                &p.iu,
                &p.abstol,
                &p.m,
                p.w,
                p.z,
                &p.ldz,
                p.work,
                &p.lwork,
                p.iwork,
                p.ifail,
                &p.info
                );
    }

} // namespace lapack



template< size_t N, typename float_t >
struct lapack_sym_eigs
{

    typedef matrix< N, N, float_t > m_input_type;
    typedef matrix< N, N, float_t > evectors_type;

    typedef vector< N, float_t > evalues_type;
    typedef vector< N, float_t > evector_type;
    typedef typename evalues_type::iterator evalue_iterator;
    typedef typename evalues_type::const_iterator evalue_const_iterator;

    typedef std::pair< float_t, size_t >  eigv_pair_type;

    lapack_sym_eigs();
    ~lapack_sym_eigs();

    // version of reduced sym. eigenvalue decomposition,
    // computes only the x largest magn. eigenvalues and their corresponding eigenvectors
    template< size_t X>
    bool compute_x(
                     const m_input_type& A,
                     matrix< N, X, float_t >& eigvectors_,
                     vector< X, float_t >& eigvalues_
                     );

    // partial sym. eigenvalue decomposition
    // returns only the largest magn. eigenvalue and the corresponding eigenvector
    bool compute_1st(
                   const m_input_type& A,
                   evector_type& eigvector_,
                   float_t& eigvalue_
                   );


    //computes all eigenvalues and eigenvectors for matrix A
    bool compute_all(
                 const m_input_type& A,
                 evectors_type& eigvectors_,
                 evalues_type& eigvalues_
                 );

    inline bool test_success( lapack::lapack_int info );

    lapack::eigs_params< float_t > p;

    const lapack::eigs_params< float_t >& get_params(){ return p; };

    // comparison functor
    struct eigenvalue_compare
    {
        inline bool operator()( const eigv_pair_type& a, const eigv_pair_type& b )
        {
            return fabs( a.first ) > fabs( b.first );
        }
    };

}; // struct lapack_sym_eigs


template< size_t N, typename float_t >
lapack_sym_eigs< N, float_t >::lapack_sym_eigs()
{
    p.jobz      = 'V'; // Compute eigenvalues and eigenvectors.
    p.range     = 'A'; // all eigenvalues will be found.
    p.uplo      = 'U'; // Upper triangle of A is stored; or Lower triangle of A is stored.
    p.n         = N;
    p.a         = 0; //If UPLO = 'U', the leading N-by-N upper triangular part of A contains the upper triangular part of the matrix A.
    p.lda       = N;
    p.vl        = 0; //Not referenced if RANGE = 'A' or 'I'.
    p.vu        = 0; //Not referenced if RANGE = 'A' or 'I'.
    p.il        = 0; //Not referenced if RANGE = 'A' or 'V'.
    p.iu        = 0; //Not referenced if RANGE = 'A' or 'V'.
    p.abstol    = 0.0001; //lie in an interval [a,b] of width less than or equal to ABSTOL + EPS *   max( |a|,|b| )
    p.m         = N; //The total number of eigenvalues found.  0 <= M <= N.
    p.w         = 0; //first m eigenvalues
    p.z         = 0; //first m eigenvectors
    p.ldz       = N; // The leading dimension of the array Z.  LDZ >= 1, and if JOBZ = 'V', LDZ >= max(1,N).
    p.work      = new float_t;
    //FIXME: check if correct datatype
    p.iwork     = new lapack::lapack_int[5*N]; //[5*N]; // INTEGER array, dimension (5*N)
    p.ifail     = new lapack::lapack_int[N]; //[N];
    p.lwork     = -1; //8N

    // workspace query
    lapack::sym_eigs_call( p );

    p.lwork = static_cast< lapack::lapack_int >( p.work[0] );
    delete p.work;

    p.work = new float_t[ p.lwork ];

}



template< size_t N, typename float_t >
lapack_sym_eigs< N, float_t >::~lapack_sym_eigs()
{
    delete[] p.work;
    delete[] p.iwork;
    delete[] p.ifail;
}

template< size_t N, typename float_t >
bool
lapack_sym_eigs< N, float_t >::compute_all(
                                        const m_input_type& A,
                                        evectors_type& eigvectors_,
                                        evalues_type& eigvalues_
                                        )
{
    // lapack destroys the contents of the input matrix
    m_input_type* AA = new m_input_type( A );

    p.range     = 'A'; // all eigenvalues will be found.
    p.a         = AA->array;
    p.ldz       = N;
    p.w         = eigvalues_.array;
    p.z         = eigvectors_.array;

    //debug std::cout << p << std::endl;

    lapack::sym_eigs_call< float_t >( p );

    delete AA;

    return p.info == 0;
}


template< size_t N, typename float_t >
template< size_t X >
bool
lapack_sym_eigs< N, float_t >::compute_x(
                                        const m_input_type& A,
                                        matrix< N, X, float_t >& eigvectors_,
                                        vector< X, float_t >& eigvalues_
                                        )
{
    //(1) get all eigenvalues and eigenvectors
    evectors_type* all_eigvectors = new evectors_type;
    evalues_type* all_eigvalues = new evalues_type;
    compute_all( A, *all_eigvectors, *all_eigvalues );

    //(2) sort the eigenvalues
    //std::pair< data, original_index >;
    std::vector< eigv_pair_type >* eig_permutations = new std::vector< eigv_pair_type >;

    evalue_const_iterator it = all_eigvalues->begin(), it_end = all_eigvalues->end();
    size_t counter = 0;
    for( ; it != it_end; ++it, ++counter )
    {
        eig_permutations->push_back( eigv_pair_type( *it, counter ) );
    }

    std::sort(
                eig_permutations->begin(),
                eig_permutations->end(),
                eigenvalue_compare()
              );

    delete all_eigvalues;

    //sort the eigenvectors according to eigenvalue permutations
    evectors_type* sorted_eigvectors = new evectors_type;
    evalues_type* sorted_eigvalues = new evalues_type;
    typename std::vector< eigv_pair_type >::const_iterator it2 = eig_permutations->begin(), it2_end = eig_permutations->end();
    evalue_iterator evalues_it = sorted_eigvalues->begin();

    for( counter = 0; it2 != it2_end; ++it2, ++evalues_it, ++counter )
    {
        *evalues_it = it2->first;
        sorted_eigvectors->set_column( counter, all_eigvectors->get_column( it2->second ));
    }
    delete all_eigvectors;
    delete eig_permutations;

    //(3) select the largest magnitude eigenvalues and the corresponding eigenvectors
    typename vector< X, float_t >::iterator it3 = eigvalues_.begin(), it3_end = eigvalues_.end();
    evalues_it = sorted_eigvalues->begin();
    for( ; it3 != it3_end; ++it3, ++evalues_it )
    {
        *it3 = *evalues_it;
    }

    sorted_eigvectors->get_sub_matrix( eigvectors_ );

    delete sorted_eigvectors;
    delete sorted_eigvalues;

    return p.info == 0;
}

template< size_t N, typename float_t >
bool
lapack_sym_eigs< N, float_t >::compute_1st(
                                         const m_input_type& A,
                                         evector_type& eigvector_,
                                         float_t& eigvalue_
                                         )
{
    //(1) get all eigenvalues and eigenvectors
    evectors_type* all_eigvectors = new evectors_type;
    evalues_type* all_eigvalues = new evalues_type;
    compute_all( A, *all_eigvectors, *all_eigvalues );

    //(2) sort the eigenvalues
    //std::pair< data, original_index >;
    std::vector< eigv_pair_type >* eig_permutations = new std::vector< eigv_pair_type >;

    evalue_const_iterator it = all_eigvalues->begin(), it_end = all_eigvalues->end();
    size_t counter = 0;
    for( ; it != it_end; ++it, ++counter )
    {
        eig_permutations->push_back( eigv_pair_type( *it, counter ) );
    }

    std::sort(
              eig_permutations->begin(),
              eig_permutations->end(),
              eigenvalue_compare()
              );

    //(2) select the largest magnitude eigenvalue and the corresponding eigenvector
    typename std::vector< eigv_pair_type >::const_iterator it2 = eig_permutations->begin();

    eigvalue_ = it2->first;
    all_eigvectors->get_column( it2->second, eigvector_ );

    delete eig_permutations;
    delete all_eigvalues;
    delete all_eigvectors;

    return p.info == 0;
}




} // namespace vmml

#endif