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

#ifndef _theplu_yat_utility_svd_
#define _theplu_yat_utility_svd_

// $Id: SVD.h 1474 2008-09-04 11:43:27Z 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' = (MxM)(MxN)(NxN) = (MxN)\n   

     A = Matrix to be decomposed, size MxN\n
     U = Orthogonal matrix, size MxM\n
     S = Diagonal matrix of singular values, size MxN\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, check the GSL documentation.
    
       There are restrictions on the matrix dimensions depending on
       which SVD algorithm is used. From the GSL's SVD source code one
       finds that the Golub-Reinsch algorithm implementation will not
       work on matrices with fewer rows than columns, the same is also
       true for the modified Golub-Reinsch algorithm.
    
       \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.
    */
    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