householder_transform.h (60)

/*
    $Id: householder_transform.h 60 2009-03-01 18:26:41Z 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 _HOUSEHOLDER_TRANSFORM_H
#define _HOUSEHOLDER_TRANSFORM_H

#include <lsp/utils.h>

#include <algorithm>
#include <limits>

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

using namespace boost::numeric::ublas;

namespace lsp{

template< class T > class householder_transform {
public:
        typedef T                                vector_type; 
        typedef typename vector_type::value_type value_type;  
        typedef typename vector_type::size_type  size_type;   

private:
        size_type  m_l;
        size_type  m_p;
        value_type m_s;
        value_type m_h;
        vector_type m_v;

public:
        householder_transform( size_type l, size_type p, vector_type v ):
                m_l( l ), m_p( p ), m_s( 0 ), m_v(v) {
                const size_type m = v.size();

                assert( p < l );
                assert( p >= 0 );
                assert( p < m );

                value_type w;

                if( l < m )
                        w = std::abs( std::max( v(p), *( std::max_element( v.begin() + l, v.end(), less_abs< value_type >() ) ), less_abs< value_type >() ) );
                else {
                        m_h = 2 * v(p);
                        m_s = -v(p);
                        return;
                }

                if( w != 0 ){
                        m_s += ( v(p)/w )*( v(p)/w );
                        for( size_type i = l; i < m; ++i )
                                m_s += ( v(i)/w )*( v(i)/w );
                        m_s = ( v(p) < 0 ? 1 : -1 ) * w * std::pow( m_s, 0.5 );
                } else {
                        m_s = 0;
                }

                m_h = v(p) - m_s;
        }

        template<class M> void apply ( matrix_row<M> v ) const {
                apply( v, vector_tag() );
        }
        template<class M> void apply ( matrix_column<M> v ) const {
                apply( v, vector_tag() );
        }
        template<class U> void apply ( U& v ) const {
                apply( v, vector_tag() );
        }
        template<class U> void apply ( U v, vector_tag ) const {
                typedef U vector2_type;

                assert( m_v.size() == v.size() );

                const value_type b = m_s * m_h;
                if( b == 0 )
                        return;

                value_type s = v(m_p) * m_h;
                for( size_type i = m_l; i < v.size(); i++ )
                        s += v(i) * m_v(i);
                s = s / b;

                v(m_p) += s * m_h;
                for( size_type i = m_l; i < v.size(); i++ )
                        v(i) += s * m_v(i);
        }
        template<class U> void apply ( U& w, row_major_tag ) const {
                for( size_type i = 0; i < w.size2(); ++i )
                        apply( column( w, i ) );
        }
        template<class U> void apply ( U& w, column_major_tag ) const {
                for( size_type i = 0; i < w.size1(); ++i )
                        apply( row( w, i ) );
        }

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

};

#endif // _HOUSEHOLDER_TRANSFORM_H

---