source: trunk/test/svd_test.cc @ 680

// $Id: svd_test.cc 680 2006-10-11 17:49:03Z jari $

/*
  Copyright (C) The authors contributing to this file.

  This file is part of the yat library, http://lev.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)
{
  // accepted error, Jari: should be picked up from GSL
  double MAXTOL=1e-13;

  // 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);
}