Re: [ublas] ublas atlas bindings

From: Marco (brbromo_at_[hidden])
Date: 2006-09-13 11:20:09

Just another tip... in function

gesvd(Jobu,Jobvt,A,s,U,VT)

described in
http://boost.cvs.sourceforge.net/*checkout*/boost-sandbox/boost-sandbox/libs/numeric/bindings/lapack/doc/index.html#svd-s

it seems that lapack is returning the transpose of matrices U and VT.

I'm attaching some code.

Thanks,
Marco

{CODE}
#include <iostream>

#include <boost/numeric/bindings/atlas/clapack.hpp>
#include <boost/numeric/bindings/traits/ublas_matrix.hpp>
#include <boost/numeric/bindings/traits/std_vector.hpp>
#include <boost/numeric/bindings/lapack/gesvd.hpp>

#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/vector_proxy.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/triangular.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/io.hpp>

namespace ublas  = boost::numeric::ublas;
namespace atlas  = boost::numeric::bindings::atlas;
namespace lapack = boost::numeric::bindings::lapack;

void print(ublas::matrix<double> const& m) {
  std::cout << m(0,0) << "\t"  << m(0,1) << "\t"  << m(0,2) << std::endl;
  std::cout << m(1,0) << "\t"  << m(1,1) << "\t"  << m(1,2) << std::endl;
  std::cout << m(2,0) << "\t"  << m(2,1) << "\t"  << m(2,2) << std::endl;
}

void print(std::vector<double> const& s) {
  std::cout << s[0] << std::endl;
  std::cout << s[1] << std::endl;
  std::cout << s[2] << std::endl;
}

int main() {
  ublas::matrix<double> H (3, 3);

  H(0,0) = 1.00000; H(0,1) = 0.50000; H(0,2) = 0.33333;
  H(1,0) = 0.50000; H(1,1) = 0.33333; H(1,2) = 0.25000;
  H(2,0) = 0.33333; H(2,1) = 0.25000; H(2,2) = 0.20000;

  ublas::matrix<double> U (3, 3);
  std::vector<double> s (3);
  ublas::matrix<double> VT (3, 3);

  lapack::gesvd('A', 'A', H, s, U, VT);

  print(U);
  std::cout << std::endl;
  print(s);
  std::cout << std::endl;
  print(ublas::trans(VT));
  std::cout << std::endl;

  return 0;
 }

{OUTPUT}
-0.827046	-0.459864	-0.323297
0.547445	-0.528278	-0.64902
0.12767	-0.713756	0.68866

1.40832
0.122329
0.00268506

-0.827046	-0.459864	-0.323297
0.547445	-0.528278	-0.64902
0.12767	-0.713756	0.68866

{EXPECTED OUTPUT RESULTS}
octave:2> [U s V] = svd(hilb(3))
U =

  -0.82704   0.54745   0.12766
  -0.45986  -0.52829  -0.71375
  -0.32330  -0.64901   0.68867

s =

   1.40832   0.00000   0.00000
   0.00000   0.12233   0.00000
   0.00000   0.00000   0.00269

V =

  -0.82704   0.54745   0.12766
  -0.45986  -0.52829  -0.71375
  -0.32330  -0.64901   0.68867