givens_rotation.h (69)

/*
    $Id: givens_rotation.h 69 2009-08-28 18:32:24Z matwey $
    Copyright (C) 2008  Matwey V. Kornilov <matwey.kornilov@gmail.com>

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

#ifndef _GIVENS_ROTATION_H
#define _GIVENS_ROTATION_H

#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>

using namespace boost::numeric::ublas;

namespace lsp{

template< class T > class givens_rotation{
public:
        typedef T value_type; 

private:
        value_type m_c, m_s;

public:
        givens_rotation( value_type& x, value_type& y ) {
                value_type w,q;

                if( std::abs( x ) <= std::abs( y ) ) {
                        if( y == 0 ) {
                                m_c = value_type( 1 );
                                m_s = value_type( 0 );
                        } else {
                                w = x / y;
                                q = std::sqrt( value_type( 1 ) + w*w );
                                m_s = value_type( 1 ) / q;
                                if( y < 0 )
                                        m_s = -m_s;
                                m_c = w * m_s;
                                x = std::abs( y * q );
                                y = value_type( 0 );
                        }
                } else {
                        w = y / x;
                        q = std::sqrt( value_type( 1 ) + w*w );
                        m_c = value_type( 1 ) / q;
                        if( x < 0 )
                                m_c = -m_c;
                        m_s = w * m_c;
                        x = std::abs( x * q );
                        y = value_type( 0 );
                }
        }

        template<class U> void apply ( U& x, U& y ) const {
                U w ( x * m_c + y * m_s );
                y = x * ( -m_s ) + y * m_c;
                x = w;
        }
        template<class M> void apply ( matrix_row< M > x, matrix_row< M > y ) const {
                typename vector_temporary_traits< matrix_row< M > >::type w ( x * m_c + y * m_s );
                y = x * ( -m_s ) + y * m_c;
                x = w;
        }
        template<class M> void apply ( matrix_column< M > x, matrix_column< M > y ) const {
                typename vector_temporary_traits< matrix_column< M > >::type w ( x * m_c + y * m_s );
                y = x * ( -m_s ) + y * m_c;
                x = w;
        }

        inline const value_type c() const { return m_c; }

        inline const value_type s() const { return m_s; }
};

};

#endif // _GIVENS_ROTATION_H

---