lapack_svd.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_SVD__HPP__
#define __VMML__VMMLIB_LAPACK_SVD__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>

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 svd_params
{
    char            jobu;
    char            jobvt;
    lapack_int      m;
    lapack_int      n;
    float_t*        a;
    lapack_int      lda;
    float_t*        s;
    float_t*        u;
    lapack_int      ldu;
    float_t*        vt;
    lapack_int      ldvt;
    float_t*        work;
    lapack_int      lwork;
    lapack_int      info;

    friend std::ostream& operator << ( std::ostream& os,
        const svd_params< float_t >& p )
    {
        os
            << "jobu "      << p.jobu
            << " jobvt "    << p.jobvt
            << " m "        << p.m
            << " n "        << p.n
            << " lda "      << p.lda
            << " ldu "      << p.ldu
            << " ldvt "     << p.ldvt
            << " lwork "    << p.lwork
            << " info "     << p.info
            << std::endl;
        return os;
    }

};


#if 0
/* Subroutine */ int dgesvd_(char *jobu, char *jobvt, integer *m, integer *n,
    doublereal *a, integer *lda, doublereal *s, doublereal *u, integer *
    ldu, doublereal *vt, integer *ldvt, doublereal *work, integer *lwork,
    integer *info);
#endif


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


template<>
inline void
svd_call( svd_params< float >& p )
{
    //std::cout << "calling lapack svd (single precision) " << std::endl;
    sgesvd_(
        &p.jobu,
        &p.jobvt,
        &p.m,
        &p.n,
        p.a,
        &p.lda,
        p.s,
        p.u,
        &p.ldu,
        p.vt,
        &p.ldvt,
        p.work,
        &p.lwork,
        &p.info
        );

}


template<>
inline void
svd_call( svd_params< double >& p )
{
    //std::cout << "calling lapack svd (double precision) " << std::endl;
    dgesvd_(
        &p.jobu,
        &p.jobvt,
        &p.m,
        &p.n,
        p.a,
        &p.lda,
        p.s,
        p.u,
        &p.ldu,
        p.vt,
        &p.ldvt,
        p.work,
        &p.lwork,
        &p.info
        );
}

} // namespace lapack



template< size_t M, size_t N, typename float_t >
struct lapack_svd
{
    lapack_svd();
    ~lapack_svd();

    // slow version, full SVD, use if all values of U(MXM) and Vt(NXN) are needed
    bool compute_full(
        const matrix< M, N, float_t >& A,
        matrix< M, M, float_t >& U,
        vector< N, float_t >& sigma,
        matrix< N, N, float_t >& Vt
        );

    // version of reduced SVD, computes only most significant left and right singular vectors,
    // i.e., use if U(MXN) and Vt(NXN) are needed
    bool compute(
                 const matrix< M, N, float_t >& A,
                 matrix< M, N, float_t >& U,
                 vector< N, float_t >& sigma,
                 matrix< N, N, float_t >& Vt
                 );

    // overwrites A with the result U,
    bool compute_and_overwrite_input(
        matrix< M, N, float_t >& A_U,
        vector< N, float_t >& sigma
        );

    // fast version, use if only sigma is needed.
    bool compute(
        const matrix< M, N, float_t >& A,
        vector< N, float_t >& sigma
        );

    inline bool test_success( lapack::lapack_int info );

    lapack::svd_params< float_t > p;

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

}; // struct lapack_svd


template< size_t M, size_t N, typename float_t >
lapack_svd< M, N, float_t >::lapack_svd()
{
    p.jobu      = 'N';
    p.jobvt     = 'N';
    p.m         = M;
    p.n         = N;
    p.a         = 0;
    p.lda       = M;
    p.s         = 0;
    p.u         = 0;
    p.ldu       = M;
    p.vt        = 0;
    p.ldvt      = 1;
    p.work      = new float_t;
    p.lwork     = -1;

    // workspace query
    lapack::svd_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 M, size_t N, typename float_t >
lapack_svd< M, N, float_t >::~lapack_svd()
{
    delete[] p.work;
}



template< size_t M, size_t N, typename float_t >
bool
lapack_svd< M, N, float_t >::compute_full(
    const matrix< M, N, float_t >& A,
    matrix< M, M, float_t >& U,
    vector< N, float_t >& S,
    matrix< N, N, float_t >& Vt
    )
{
    // lapack destroys the contents of the input matrix
    typedef matrix< M, N, float_t > m_type;
    m_type* AA = new m_type( A );

    p.jobu      = 'A';
    p.jobvt     = 'A';
    p.a         = AA->array;
    p.u         = U.array;
    p.s         = S.array;
    p.vt        = Vt.array;
    p.ldvt      = N;

    lapack::svd_call< float_t >( p );

    delete AA;

    return p.info == 0;
}

template< size_t M, size_t N, typename float_t >
bool
lapack_svd< M, N, float_t >::compute(
                                     const matrix< M, N, float_t >& A,
                                     matrix< M, N, float_t >& U,
                                     vector< N, float_t >& S,
                                     matrix< N, N, float_t >& Vt
                                     )
{
    // lapack destroys the contents of the input matrix
    typedef matrix< M, N, float_t > m_type;
    m_type* AA = new m_type( A );

    p.jobu      = 'S';
    p.jobvt     = 'S';
    p.a         = AA->array;
    p.u         = U.array;
    p.s         = S.array;
    p.vt        = Vt.array;
    p.ldvt      = N;

    lapack::svd_call< float_t >( p );
    delete AA;
    return p.info == 0;
}


template< size_t M, size_t N, typename float_t >
bool
lapack_svd< M, N, float_t >::compute_and_overwrite_input(
    matrix< M, N, float_t >& A_U,
    vector< N, float_t >& S
    )
{
    p.jobu      = 'O';
    p.jobvt     = 'N';
    p.a         = A_U.array;
    p.s         = S.array;
    p.ldvt      = N;

    lapack::svd_call< float_t >( p );

    return p.info == 0;
}

template< size_t M, size_t N, typename float_t >

bool lapack_svd< M, N, float_t >::compute( const matrix< M, N, float_t >& A,

                                           vector< N, float_t >& S )

{

    // lapack destroys the contents of the input matrix

    typedef matrix< M, N, float_t > m_type;

m_type* AA = new m_type( A );

    p.jobu      = 'N';
    p.jobvt     = 'N';
    p.a         = AA->array;
    p.u         = 0;
    p.s         = S.array;
    p.vt        = 0;

    lapack::svd_call< float_t >( p );

 delete AA;

     return p.info == 0;

}



} // namespace vmml

#endif