source: trunk/yat/utility/SVD.h @ 1478

#ifndef _theplu_yat_utility_svd_
#define _theplu_yat_utility_svd_

// $Id: SVD.h 1478 2008-09-06 18:55:46Z peter $

/*
  Copyright (C) 2003 Daniel Dalevi, Jari Häkkinen
  Copyright (C) 2004 Jari Häkkinen
  Copyright (C) 2005 Jari Häkkinen, Peter Johansson
  Copyright (C) 2006 Jari Häkkinen, Markus Ringnér
  Copyright (C) 2007 Jari Häkkinen, Peter Johansson
  Copyright (C) 2008 Peter Johansson

  This file is part of the yat library, http://dev.thep.lu.se/yat

  The yat library 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 2 of the
  License, or (at your option) any later version.

  The yat library 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, write to the Free Software
  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
  02111-1307, USA.
*/

#include "Matrix.h"
#include "Vector.h"

#include <gsl/gsl_linalg.h>

namespace theplu {
namespace yat {
namespace utility {

  class VectorBase;

  /**
     @brief Singular Value Decomposition

     Class encapsulating GSL methods for singular value
     decomposition, SVD.

     A = U S V' = (MxN)(NxN)(NxN) = (MxN)\n   

     A = Matrix to be decomposed, size MxN\n
     U = Orthogonal matrix, size MxN\n
     S = Diagonal matrix of singular values, size NxN\n
     V = Orthogonal matrix, size NxN\n
  */
  class SVD
  {
  public:

    /**
       A number of SVD algorithms are implemented in GSL. They have
       their strengths and weaknesses.
    
       \see GSL's SVD documentation. 
    */
    enum SVDalgorithm {
      GolubReinsch,
      ModifiedGolubReinsch,
      Jacobi
    };

    /**
       \brief Constructs an SVD object using the matrix Ain as only
       input. The input matrix is copied for further use in the
       object.

       \note Number of rows must be equal or larger than number of columns.
    */
    SVD(const utility::Matrix& Ain);

    /**
       \brief The destructor
    */
    ~SVD(void);

    /**
       \brief This function will perform SVD with the method specified
       by \a algo.

       \throw GSL_error if the underlying GSL function fails.
    */
    void decompose(SVDalgorithm algo=GolubReinsch);

    /**
       \brief Access to the s vector.

       \return A copy of the s vector.

       \note If decompose() has not been run the outcome of the call
       is undefined.
    */
    const utility::Vector& s(void) const;

    /**
       \brief Solve the system \f$ Ax=b \f$ using the decomposition of
       A.

       \note If decompose() has not been run the outcome of the call
       is undefined.

       \throw GSL_error if the underlying GSL function fails.
    */
    void solve(const utility::VectorBase& b, utility::Vector& x);

    /**
       \brief Access to the U matrix.

       \return A copy of the U matrix.

       \note If decompose() has not been run the outcome of the call
       is undefined.
    */
    const utility::Matrix& U(void) const;

    /**
       \brief Access to the V matrix.

       \return A copy of the V matrix.

       \note If decompose() has not been run the outcome of the call
       is undefined.
    */
    const utility::Matrix& V(void) const;

  private:
    /**
       \brief Call GSL's Jacobi algorithm for SVD.

       \return gsl_error status.
    */
    int jacobi(void);

    /**
       \brief Call GSL's Golub-Reinsch algorithm for SVD.

       \return gsl_error status.
    */
    int golub_reinsch(void);

    /**
       \brief Call GSL's modified Golub-Reinch algorithm for SVD.

       \return gsl_error status.
    */
    int modified_golub_reinsch(void);

    utility::Matrix U_, V_;
    utility::Vector s_;
  };  

}}} // of namespace utility, yat, and theplu

#endif