include/ewalena/lac/solver_qr.h

Go to the documentation of this file.

// -------------------------------------------------------------------
// @author Toby D. Young
// @version $Id: solver_qr.h 987 2012-11-26 15:16:27Z oneliefleft $
//
// Copyright 2012 namespace ewalena authors. All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions
// are met:
//
// 1. Redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//
// 2. 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.
//
// THIS SOFTWARE IS PROVIDED BY THE NAMEPSACE EWALENA AUTHORS ``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
// NAMESPACE EWALENA AUTHORS 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.
//
// The views and conclusions contained in the software and
// documentation are those of the authors and should not be
// interpreted as representing official policies, either expressed or
// implied, of the namespace ewalena authors.
// -------------------------------------------------------------------

#ifndef __ewalena_solver_qr_h
#define __ewalena_solver_qr_h

#ifdef EWALENA_HAVE_CONFIG_H
#include <ewalena/base/config.h>
#endif

#include <ewalena/lac/solver_base.h>

#include <algorithm>

namespace ewalena 
{

  template <typename ValueType = double>
    class SolverQR
    : 
    public SolverBase<ValueType>
    {
      
    public:

    SolverQR (SolverControl<ValueType> &solver_control);
    
    virtual void solve (const Matrix<ValueType> &A,
                        Vector<ValueType>       &lambda,
                        VectorBasis<ValueType>  &x);
    
    protected:
    
    void reduce_to_hessenberg_form (Matrix<ValueType> &matrix);
    
    }; /* SolverJacobi */

  /*-------------- Inline and Other Functions -----------------------*/
  
  template <typename ValueType>
    inline
    ValueType
    sign (const ValueType a,
          const ValueType b)
    {
      return (b>ValueType (0)) ? fabs (a) : -fabs (a);
    }
  
  template <typename ValueType>
    inline
    void
    SolverQR<ValueType>::reduce_to_hessenberg_form (Matrix<ValueType> &matrix)
    {
      assert (matrix.n_rows () == matrix.n_cols ());
      assert (matrix.n_cols () != 0);

      /* Trivial Hessenberg form already? */
      if (matrix.n_rows () < 3)
        return;  

      /* m = r+1 in "Numerical Recipes in C++". */
      for (unsigned int m=1; m<(matrix.n_rows ()-1); ++m)
        {
          unsigned int i = m;
          ValueType    x = ValueType (0);

          /* Find a pivot. */
          for (unsigned int j=m; j<matrix.n_rows (); ++j)
            if (std::fabs (matrix(j,m-1)) > std::fabs (x))
              {
                x = matrix(j,m-1); 
                i = j; 
              }

          /* Interchange rows and columns. */
          if (i != m)
            {
              for (unsigned int j=m-1; j<matrix.n_rows (); ++j)
                std::swap (matrix(i,j), matrix(m,j));
              
              for (unsigned int j=0; j<matrix.n_rows (); ++j)
                std::swap (matrix(j,i), matrix(j,m));
            }

          /* Do elimination. */
          if (x)
            {
              for (unsigned int i=m+1; i<matrix.n_rows (); ++i)
                {
                  ValueType y = matrix(i,m-1);
                  if (y != ValueType(0))
                    {
                      y        /= x;
                      matrix(i,m-1)  = y;
                      
                      for (unsigned int j=m; j<matrix.n_rows (); ++j)
                        matrix(i,j) -= y*matrix(m,j);

                      for (unsigned int j=0; j<matrix.n_rows (); ++j)
                        matrix(j,m) += y*matrix(j,i) ;
                    }
                }
            }
          
        }


      /* Explicitly zero out lower elements here (remove round-off
         errors and/or overflow). */
      for (unsigned int j=0; j<matrix.n_cols (); ++j)
        for (unsigned int i=j+2; i<matrix.n_rows (); ++i)
          matrix(i,j) = ValueType(0);
      
    }

} /* namespace ewalena */

#endif /* __ewalena_solver_qr_h */