source: trunk/test/svd_test.cc @ 865

// $Id: svd_test.cc 865 2007-09-10 19:41:04Z peter $

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

  This file is part of the yat library, http://trac.thep.lu.se/trac/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 "yat/random/random.h"
#include "yat/utility/matrix.h"
#include "yat/utility/SVD.h"
#include "yat/utility/vector.h"

using namespace theplu::yat;

double this_norm(const utility::matrix& A)
{
  double sum=0.0;
  for (size_t i=0; i<A.rows(); ++i)
    for (size_t j=0; j<A.columns(); ++j)
      sum += A(i,j)*A(i,j);
  return sum;
}



bool test(size_t m, size_t n, utility::SVD::SVDalgorithm algo)
{
  double MAXTOL=1e-13; // accepted error

  // initialise a random test-matrix
  theplu::yat::random::ContinuousUniform rnd;
  utility::matrix A(m,n);
  for (size_t i=0; i<m; ++i)
    for(size_t j=0; j<n; ++j)
      A(i,j)=1000*rnd();

  utility::SVD svd(A);
  svd.decompose(algo);
  theplu::yat::utility::vector s(svd.s());
  utility::matrix S(s.size(),s.size());
  for (size_t i=0; i<s.size(); ++i)
    S(i,i)=s[i];
  utility::matrix Vtranspose=svd.V();
  Vtranspose.transpose();
  // Reconstructing A = U*S*Vtranspose
  utility::matrix Areconstruct(svd.U());
  Areconstruct*=S;
  Areconstruct*=Vtranspose;
  Areconstruct-=A;  // Expect null matrix
  double error = this_norm(Areconstruct);
  bool testerror=false;
  if (error>MAXTOL) {
    std::cerr << "test_svd: FAILED, algorithm " << algo 
              << " recontruction error ("
              << error << ") > tolerance (" << MAXTOL << "), matrix dimension ("
              << m << ',' << n << ')' << std::endl;
    testerror=true;
  }

  Vtranspose*=svd.V();  // Expect unity matrix
  error=this_norm(Vtranspose)-n;
  if (error>MAXTOL) {
    std::cerr << "test_svd: FAILED, algorithm " << algo 
              << " V orthogonality error ("
              << error << ") > tolerance (" << MAXTOL << ')' << std::endl;
    testerror=true;
  }

  utility::matrix Utranspose(svd.U());
  Utranspose.transpose();
  Utranspose*=svd.U();  // Expect unity matrix
  error=this_norm(Utranspose)-n;
  if (error>MAXTOL) {
    std::cerr << "test_svd: FAILED, algorithm " << algo 
              << " U orthogonality error ("
              << error << ") > tolerance (" << MAXTOL << ')' << std::endl;
    testerror=true;
  }
  return testerror;
}



int main(const int argc,const char* argv[])
{
  bool testfail=false;

  // The GSL Jacobi, Golub-Reinsch, and modified Golub-Reinsch
  // implementations supports rows>=columns matrix dimensions only
  testfail|=test(12,12,utility::SVD::GolubReinsch);
  testfail|=test(12,4,utility::SVD::GolubReinsch);
  testfail|=test(12,12,utility::SVD::ModifiedGolubReinsch);
  testfail|=test(12,4,utility::SVD::ModifiedGolubReinsch);
  testfail|=test(12,12,utility::SVD::Jacobi);
  testfail|=test(12,4,utility::SVD::Jacobi);

  return (testfail ? -1 : 0);
}