linear_algebra_svd.hpp

Go to the documentation of this file.

/*
 * Copyright(c):
 * Signal Image and Communications (SIC) Department 
 * http://www.sic.sp2mi.univ-poitiers.fr/ 
 *  - University of Poitiers, France http://www.univ-poitiers.fr 
 *  - XLIM  Institute UMR CNRS 7252  http://www.xlim.fr/ 
 * 
 * and 
 * 
 * D2 Fluid, Thermic and Combustion 
 *  - University of Poitiers, France http://www.univ-poitiers.fr 
 *  - PPRIME Institute -  UPR CNRS 3346  http://www.pprime.fr
 *  - ISAE-ENSMA http://www.ensma.fr
 *
 * Contributor(s):
 * The SLIP team,
 * Benoit Tremblais <tremblais_AT_sic.univ-poitiers.fr>,
 * Laurent David <laurent.david_AT_lea.univ-poitiers.fr>, 
 * Ludovic Chatellier <ludovic.chatellier_AT_univ-poitiers.fr>, 
 * Lionel Thomas <lionel.thomas_AT_univ-poitiers.fr>, 
 * Denis Arrivault <arrivault_AT_sic.univ-poitiers.fr>, 
 * Julien Dombre <julien.dombre_AT_univ-poitiers.fr>.
 *
 * Description:
 * The Simple Library of Image Processing (SLIP) is a new image processing 
 * library. It is written in the C++ language following as much as possible 
 * the ISO/ANSI C++ standard. It is consequently compatible with any system  
 * satisfying the ANSI C++ complience. It works on different Unix , Linux ,  
 * Mircrosoft Windows and Mac OS X plateforms. SLIP is a research library that  
 * was created by the Signal, Image and Communications (SIC) departement of  
 * the XLIM, UMR 7252 CNRS Institute in collaboration with the Fluids, Thermic  
 * and Combustion departement of the P', UPR 3346 CNRS Institute of the  
 * University of Poitiers.
 * 
 * The SLIP Library source code has been registered to the APP (French Agency 
 * for the Protection of Programs) by the University of Poitiers and CNRS, 
 * under  registration number IDDN.FR.001.300034.000.S.P.2010.000.21000.

 * http://www.sic.sp2mi.univ-poitiers.fr/slip/
 *
 * This software is governed by the CeCILL-C license under French law and
 * abiding by the rules of distribution of free software.  You can  use, 
 * modify and/ or redistribute the software under the terms of the CeCILL-C 
 * license as circulated by CEA, CNRS and INRIA at the following URL 
 * http://www.cecill.info.
 * As a counterpart to the access to the source code and  rights to copy,
 * modify and redistribute granted by the license, users are provided only
 * with a limited warranty  and the software's author,  the holder of the
 * economic rights,  and the successive licensors  have only  limited
 * liability.
 * 
 * In this respect, the user's attention is drawn to the risks associated
 * with loading,  using,  modifying and/or developing or reproducing the
 * software by the user in light of its specific status of free software,
 * that may mean  that it is complicated to manipulate,  and  that  also
 * therefore means  that it is reserved for developers  and  experienced
 * professionals having in-depth computer knowledge. Users are therefore
 * encouraged to load and test the software's suitability as regards their
 * requirements in conditions enabling the security of their systems and/or
 * data to be ensured and,  more generally, to use and operate it in the 
 * same conditions as regards security. 
 *
 * The fact that you are presently reading this means that you have had
 * knowledge of the CeCILL-C license and that you accept its terms.
 */


#ifndef SLIP_LINEAR_ALGEBRA_SVD_HPP
#define SLIP_LINEAR_ALGEBRA_SVD_HPP
#include <cmath>
#include <algorithm>
#include <limits>
#include "Vector.hpp"
#include "Array2d.hpp"
#include "linear_algebra.hpp"
#include "macros.hpp"
#include "norms.hpp"
#include "linear_algebra_traits.hpp"
#include "complex_cast.hpp"

namespace slip
{

  /* @{ */  
template<typename Matrix1, 
                 typename Vector,
                 typename Matrix2,
                 typename Matrix3>

void BiReduce( Matrix1 &X,
               Vector &d,
               Vector &e,
               Matrix2 &U,
               Matrix3 &V)

 {
   assert(X.rows() >= X.rows());
   assert(U.rows() == X.rows());
   assert(U.cols() == X.cols());
   assert(V.rows() == V.cols());
   assert(V.rows() == X.cols());
   assert(d.size() == X.cols());
   assert((e.size() + 1) == X.cols());
   
   typedef typename Matrix1::value_type value_type;
   typedef typename slip::lin_alg_traits<typename Matrix1::value_type>::value_type norm_type;
   std::size_t X_rows = X.rows();
   std::size_t X_cols = X.cols();

   std::fill(U.begin(),U.end(),value_type(0));
   std::fill(V.begin(),V.end(),value_type(0));

   // temporary variables used for householder methode
   slip::Vector<value_type > a(X_rows,value_type(0));//a
   slip::Vector<value_type > b(X_cols-1,value_type(0));//b
   
   slip::Box2d<int> box;
   slip::Box2d<int> box2;
   
   for(std::size_t k=0;k<X_cols;++k)
     {
       if((k<X_cols-1)||(X_rows>X_cols))
                 { 
                   slip::housegen(X.col_begin(k)+k,X.col_end(k),a.begin()+(k),a.end(),d[k]);
                   //X[k:n,k]=a
                   std::copy(a.begin()+(k),a.end(),X.col_begin(k)+k);
                   box.set_coord(int(k),int(k+1),int(X_rows-1),int(X_cols-1));
                   //X[k:n,k+1:p]
                   slip::left_householder_update(a.begin()+k,a.end(),
                                                                                 X.upper_left(box),X.bottom_right(box));
                                  
                 }

       if(k<X_cols-2)
                 {
                   //housgen(X[k,k+1:p],b^T,e[k])
                   slip::housegen(X.row_begin(k)+(k+1),X.row_end(k),b.begin()+(k),b.end(),e[k]);
                   //X[k,k+1:p]=b^T
                   std::copy(b.begin()+(k),b.end(),X.row_begin(k)+(k+1));
                  box2.set_coord(int(k+1),int(k+1),int(X_rows-1),int(X_cols-1));
                   //M=X[k+1:n,k+1:p]
                   slip::right_householder_update(b.begin()+(k),b.end(),
                                                                                  X.upper_left(box2),X.bottom_right(box2));
                                   
                
                 }

     }

   e[X_cols-2]=X[X_cols-2][X_cols-1];

   if((X_rows)==(X_cols))
     {
       d[X_cols-1]=X[X_cols-1][X_cols-1];
     }

   //Accumulate U
   // U=(   I_p  )
   //    ( 0_n-p,p)
   U.fill(value_type(0));
   std::size_t U_cols = X.cols();
   for (std::size_t i = 0; i < U_cols; ++i)
     {
       U[i][i] = value_type(1);
     }
     
   std::size_t min_np=std::min(X_cols-1,X_rows-2);
   for(int k=min_np;k>=0;k--)
     {
       std::copy(X.col_begin(k)+k,X.col_end(k),a.begin()+k);
       box.set_coord(int(k),int(k),int(X.rows()-1),int(X.cols()-1));//U[k:n,k:p]
       slip::left_householder_update(a.begin()+k,a.end(),
                                                                     U.upper_left(box),U.bottom_right(box));
     }
   //Accumulate V
   //V=I_p()

   V.fill(value_type(0.0));
   std::size_t V_cols = X.cols();
   for (std::size_t i = 0; i < V_cols; ++i)
     {
       V[i][i] = value_type(1.0);
     }

   for(int k=int((X_cols-3));k>=0;k--)
     {
       std::copy(X.row_begin(k)+(k+1),X.row_end(k),b.begin()+k);
       box2.set_coord(int(k+1),int(k+1),int(X_cols-1),int(X_cols-1));
       slip::left_householder_update(b.begin()+(k),b.end(),
                                                     V.upper_left(box2),V.bottom_right(box2));
    
     }

}

template<typename Matrix,typename Vector>
void Bi_complex_to_real(Vector &d,Vector &e,
                        Matrix &U,
                        Matrix &V)
{
  assert(d.size() == U.cols());
  assert((d.size()-1) == e.size());
  assert(V.rows() == V.cols());
  assert(U.cols() == V.rows());
  
   typedef typename Matrix::value_type value_type;
   std::size_t d_size = d.size();//d and X have the same size : d.size()=X.cols();
   value_type mu;
   double a;
  for(std::size_t k=0;k<d_size;k++)
  {
    a=std::abs(d[k]);
    if(a!=0.0)
    {
      mu=slip::conj(d[k])/a;
      d[k]=a;
      if(k!=d_size-1)
       {
        e[k]=mu*e[k];
       }
       slip::vector_scalar_multiplies(U.col_begin(k),U.col_end(k),slip::conj(mu),U.col_begin(k));//U[:,k]=conj(mu)*U[:,k];
    } 
  if(k==(d_size-1)) break;
 
   a=std::abs(e[k]);
  if(a!=value_type(0))
    {
      mu=slip::conj(e[k])/a;
      e[k]=a;
      d[k+1]=mu*d[k+1];
      slip::vector_scalar_multiplies(V.col_begin(k+1),V.col_end(k+1),mu,V.col_begin(k+1));
   }

  }
}



template<typename Matrix,typename Vector, typename Real,typename SizeType >
void BiQR( Vector &d,
           Vector &e,
           Real &sigma,
           const SizeType& row1,
           const SizeType& row2,
           Matrix &U,
           Matrix &V)
  { 
    assert(d.size() == U.cols());
    assert((d.size()-1) == e.size());
    assert(V.rows() == V.cols());
    assert(U.cols() == V.rows());
    typedef typename Vector::value_type value_type;
    Real scl = std::max(std::abs(d[row1]),std::abs(e[row1]));
    scl = std::max(std::abs(scl),std::abs(sigma));
    Real d1  = d[row1]/scl;
    Real e1  = e[row1]/scl;
    Real sig = sigma/scl;
    Real f   = (d1+sig) * (d1-sig);
    Real g   = d1 * e1;
    SizeType row2_m1 = row2 - 1;
    SizeType V_rows_m1  = V.rows() - 1;
    SizeType U_rows_m1  = U.rows() - 1;
   
    //sinus and cosinus
    Real s = Real(0); 
    Real c = Real(0);
    for (SizeType k = row1; k <= row2_m1; ++k)
      {
                slip::rotgen(f,g,c,s);
                if(k != row1)
                  {
                    e[k-1]=f;
                  }
        f = (c*d[k]) + (s*e[k]);
                e[k] = (c*e[k]) - (s*d[k]);
                g = s * d[k+1];
                d[k+1]=c*d[k+1];
                slip::right_givens(V,k,k+1,s,c,SizeType(0),V_rows_m1);
                slip::rotgen(f,g,c,s);  
        d[k] = f;
                f = (c * e[k]) + (s * d[k+1]);
                d[k+1] = (c * d[k+1]) - (s * e[k]);
                                
                if(k < row2_m1)
                {
                  g = s * e[k+1];
                  e[k+1] = c * e[k+1];
                }
                slip::right_givens(U,k,k+1,s,c,SizeType(0),U_rows_m1);
      }
    e[row2-1] = f;
  }

template <typename Real >
void shift(Real &p,
                   Real &r, 
                   Real &q,
                   Real& sigma)
{
  Real ppq = p + q;
  Real pmq = p - q;
  Real r2 = r * r;
  Real a = (ppq * ppq) + r2;
  Real b = (pmq * pmq) + r2;
  Real sigma_max = (std::sqrt(a) + std::sqrt(b)) * Real(0.5);
  if(sigma_max != Real(0.0))
   {
    sigma = (std::abs(p*q)) / sigma_max;
   }

}
template<typename Vector, typename SizeType>
void BdBacksearch( Vector &d,
                   Vector &e,
                   const typename slip::lin_alg_traits<typename Vector::value_type>::value_type tol,
                   SizeType& l,
                   SizeType& i1,
                   SizeType& i2)
 {
   assert(d.size() == (e.size() + 1));
   assert(l < d.size());
   assert(i1 < d.size());
   assert(i2 < d.size());
   
   typedef typename Vector::value_type value_type;
   value_type tau = value_type(0);
   i1 = i2 = l;
   while(i1>0)
     {
       if(i1 == i2)
                 {
                   tau = d[i1];
                 }
       else
                 {
                   if(e[i1-1] != value_type(0.0))
                     {
                       tau = d[i1] * (tau/(tau + std::abs(e[i1-1])));
                     }
                 }
       
         if(std::abs(e[i1-1]) <= (tol*std::abs(tau)))
           {
                     e[i1-1] = value_type(0.0);
                     
                     if(i1 == i2)
                       {
                                 i2 = i1-1;
                                 i1 = i2;
                       }
                     else
                       {
                                 return;
                       }
           }
                 else
           {
                     i1 = i1 - 1;
           }
     }
 }





/* @} */

  /* @{ */
template<typename Matrix1,
                 typename Matrix2,
                 typename RandomAccessIterator,
                 typename Matrix3>
void svd(const Matrix1&M, 
                 Matrix2 &U, 
                 RandomAccessIterator S_first, 
                 RandomAccessIterator S_last,
                 Matrix3 &V, 
                 const int max_it = 75)
{
  
  assert(M.rows() >= M.cols());
  assert(M.rows() == U.rows());
  assert(M.cols() == U.cols());
  assert(int(M.cols()) == int((S_last-S_first)));
  assert(V.rows() == M.cols());
  
  typedef typename Matrix1::value_type value_type;
  //typedef typename Matrix2::value_type real;
  typedef typename slip::lin_alg_traits<typename Matrix1::value_type>::value_type real;                         
  slip::Array2d<value_type> X(M.rows(),M.cols());
  std::copy(M.begin(),M.end(),X.begin());             
  std::size_t X_cols = X.cols();
  std::size_t X_cols_m1 = X_cols - 1;
  
  // Reduce X to bidiagonal form, storing the diagonal elements
  // in d and the super-diagonal elements in e. 
  slip::Vector<value_type> d1(X.cols());
  slip::Vector<value_type> e1(X.cols()-1);
  slip::Array2d<value_type> B(X.cols(),X.cols());
  slip::BiReduce(X,d1,e1,U,V);
 
  //if X is complex then apply Bi_complex_to_real
  if(slip::is_complex(value_type()))
    {
      Bi_complex_to_real(d1,e1,U,V);
    }
     
  //copy d1 to d: to work on real data, because the type of d1 
  //by defenition is complex idem for e1
  slip::Vector<real> d(X_cols);
  slip::Vector<real> e(X_cols_m1);
  std::size_t e1_size = e1.size();
  for(std::size_t i = 0; i < e1_size; ++i)
    {
      d[i] = std::real(d1[i]);
      e[i] = std::real(e1[i]);
    }
  //d[d1.size()-1] = std::real(d1[d1.size()-1]);
  d[e1_size] = std::real(d1[e1_size]);
  
  //Main iteration loop for the singular values
  int iter = 0;
  real tol=std::numeric_limits<real>::epsilon();
  std::size_t i1 = 0;
  std::size_t i2 = X_cols_m1;
  //std::size_t l = i2;
  real sigma = real(0);
  std::size_t oldi2 = 0;
  bool error = false;
                                
  while(true)
    {
      iter++;
      if(iter > max_it)
                {
                  error = true;
                  //std::cout<<"error: no convergence"<<std::endl<<std::endl;
                  break;
                }
      oldi2 = i2;
      slip::BdBacksearch(d,e,tol,i2,i1,i2);

      if (i2 == 0) break;
      if(i2 != oldi2) 
                {
                  iter = 0;
                }
      slip::shift(d[i2-1],e[i2-1],d[i2],sigma);
      slip::BiQR(d,e,sigma,i1,i2,U,V);

    } 

 
  // Convergence.
  // Make the singular values positive.
  //if(i1==i2)
  for(std::size_t k = 0; k < X_cols; ++k)
    {
      if (d[k] <= 0.0)
                {
                  d[k] = -d[k];
                  slip::vector_scalar_multiplies(V.col_begin(k),V.col_end(k),
                                                                                 real(-1),V.col_begin(k));
                }
    }
  // sort the singular values
  for(std::size_t k = 0; k < X_cols_m1; ++k)
    {

      if(d[k] <= d[k+1])
                {
                  real t = d[k];
                  d[k] = d[k+1];
                  d[k+1] = t;
                  std::swap_ranges(V.col_begin(k),V.col_end(k),V.col_begin(k+1));
                  std::swap_ranges(U.col_begin(k),U.col_end(k),U.col_begin(k+1));
                }
}

  std::copy(d.begin(),d.end(),S_first);
 
}

template<typename Matrix,typename Matrix2>
void svd(const Matrix&X, 
                 Matrix &U, 
                 Matrix2 &S, 
                 Matrix &V, 
                 const int max_it = 75)
{
  assert(X.rows() >= X.cols());
  assert(X.rows() == U.rows());
  assert(X.cols() == U.cols());
  assert(S.rows() == S.cols());
  assert(X.cols() == S.cols());
  assert(V.rows() == X.cols());
  
  typedef typename slip::lin_alg_traits<typename Matrix::value_type>::value_type real;                          
  

  slip::Vector<real> Svect(X.cols());
  slip::svd(X,U,Svect.begin(),Svect.end(),V,max_it);
  std::fill(S.begin(),S.end(),0);
  slip::set_diagonal(Svect.begin(),Svect.end(),S);
}

/* @} */

/* @{ */
template<typename Matrix1,typename RandomAccessIterator, typename Matrix2, typename Matrix3>
inline
void svd_approx(const Matrix1 &U, 
                                RandomAccessIterator S_first, 
                                RandomAccessIterator S_last,
                                const Matrix2 &V, 
                                const std::size_t K,
                                Matrix3& X)
{
  assert(X.rows() >= X.cols());
  assert(V.rows() == V.cols());
  assert(X.rows() == U.rows());
  assert(X.cols() == U.cols());
  assert(int(S_last-S_first) == int(U.cols()));
  assert(int(S_last-S_first) == int(V.rows()));
  assert(K <= U.cols());
  
  std::size_t U_rows = U.rows();
  //std::size_t U_cols = U.cols();
  std::size_t V_rows = V.rows();
  for(std::size_t i = 0; i < U_rows; ++i)
    {
      for(std::size_t j = 0; j < V_rows; ++j)
                {
                  X[i][j] = 0;
                  for(std::size_t k = 0; k < K; ++k)
                    {
                      X[i][j] += S_first[k] * U[i][k] * slip::conj(V[j][k]);
                    }
                }
    }
  }



  template <class Matrix1,class Matrix2,class Matrix3,class Matrix4>
  inline
  void svd_approx(const Matrix1 & U,
                                  const Matrix2 & W,
                                  const Matrix3 & V,
                                  const std::size_t K,
                                  Matrix4 & M)
  {
    slip::Array<typename Matrix1::value_type> Wvect(W.rows());
    slip::get_diagonal(W,Wvect.begin(),Wvect.end());
    slip::svd_approx(U,Wvect.begin(),Wvect.end(),V,K,M);
  }  


/* @} */

  /* @{ */

template <class Matrix1, 
                  typename RandomAccessIterator1, 
                  class Matrix2, 
                  typename RandomAccessIterator2,
                  typename RandomAccessIterator3>
  inline
  void svd_solve(const Matrix1 & U,
                                 RandomAccessIterator1 W_first,
                                 RandomAccessIterator1 W_last,
                                 const Matrix2 & V,
                                 RandomAccessIterator2 B_first,
                                 RandomAccessIterator2 B_last,
                                 RandomAccessIterator3 X_first,
                                 RandomAccessIterator3 X_last)
  {
    assert(U.rows() >= U.cols());
    assert(int(X_last - X_first) == int(U.cols()));
    assert(int(B_last - B_first) == int(U.rows()));
    assert(int(W_last - W_first) == int(U.cols()));
    assert(V.rows() == U.cols());
    typedef typename Matrix1::value_type value_type;
    typedef typename slip::lin_alg_traits<typename Matrix1::value_type>::value_type real;
    //std::size_t U_rows = U.rows();
    std::size_t U_cols = U.cols();
    
    slip::Vector<value_type> tmp(U_cols);
    // Calculate U^HB.
    for (std::size_t j = 0; j < U_cols; ++j) 
      { 
                
                value_type s = value_type(0.0);
                
                // Nonzero result only if wj is nonzero.
                if (W_first[j] != real(0.0)) 
                  { 
                    s = slip::hermitian_inner_product(U.col_begin(j),U.col_end(j),
                                                                                      B_first,value_type(0.0));
                    s /= W_first[j];
                  }
                tmp[j] = s;
      }
    // Matrix multiply by V to get answer.
    slip::matrix_vector_multiplies(V.upper_left(),V.bottom_right(),
                                                                   tmp.begin(),tmp.end(),
                                                                   X_first,X_last);
  }

 template <class Matrix1, 
                  typename RandomAccessIterator1, 
                  typename RandomAccessIterator2>
  inline
  void svd_solve(const Matrix1 & A,
                                 RandomAccessIterator1 X_first,
                                 RandomAccessIterator1 X_last,
                                 RandomAccessIterator2 B_first,
                                 RandomAccessIterator2 B_last)
  {
    assert(A.rows() >= A.cols());
    assert(int(X_last - X_first) == int(A.cols()));
    assert(int(B_last - B_first) == int(A.rows()));
    
    typedef typename Matrix1::value_type value_type;
    typedef typename slip::lin_alg_traits<typename Matrix1::value_type>::value_type real;                       
  

    slip::Vector<real>  S(A.cols());
    slip::Array2d<value_type> U(A.rows(),A.cols());
    slip::Array2d<value_type> V(A.cols(),A.cols());

    //-----------------------------------------------
    // svd decomposition.
    //-----------------------------------------------
    slip::svd(A,U,S.begin(),S.end(),V);
    //------------------------------------------------------
    // svd truncature (small singular values are put to 0)
    //------------------------------------------------------
    real Smax = S[0];//*std::max_element(S.begin(),S.end());
    real Smin = Smax * slip::epsilon<value_type>();
    std::replace_if(S.begin(),S.end(),std::bind2nd(std::less<real>(), Smin),real(0));
    //-----------------------------------------------
    // svd solve
    //-----------------------------------------------
    slip::svd_solve(U,
                                    S.begin(),S.end(),
                                    V,
                                    B_first,B_last,
                                    X_first,X_last);
  }

  template <class Matrix1, 
                    class Vector1, 
                    class Vector2>
  inline
  void svd_solve(const Matrix1 & A,
                                 Vector1 & X,
                                 const Vector2 & B)
  {
    
     slip::svd_solve(A,X.begin(),X.end(),B.begin(),B.end());
  }

/* @} */

  /* @{ */
template <class Matrix1, 
                  typename RandomAccessIterator1, 
                  class Matrix2, 
                  class Matrix3>
  inline
  void svd_inv(const Matrix1 & U,
                       RandomAccessIterator1 W_first,
                       RandomAccessIterator1 W_last,
                       const Matrix2 & V,
                       Matrix3& Ainv)

  {
    typedef typename Matrix1::value_type value_type;
    typedef typename slip::lin_alg_traits<typename Matrix1::value_type>::value_type real;
    std::size_t U_rows = U.rows();
    std::size_t U_cols = U.cols();

    //------------------------------------------------------
    // svd truncature (small singular values are put to 0)
    //------------------------------------------------------
    real Wmax = *W_first;
    real Wmin = Wmax * slip::epsilon<value_type>();
    std::replace_if(W_first,W_last,std::bind2nd(std::less<real>(), Wmin),real(0));
    //------------------------------------------------------
    // find the first 0 value
    //------------------------------------------------------
    RandomAccessIterator1 it = std::find(W_first,W_last,real(0));
    std::size_t K = std::size_t(it - W_first);
     
    
    for(std::size_t i = 0; i < U_cols; ++i)
      {
                for(std::size_t j = 0; j < U_rows; ++j)
                  {
                    Ainv[i][j] = 0;
                    
                    for(std::size_t k = 0; k < K; ++k)
                      {

                                Ainv[i][j] += (real(1)/W_first[k]) * slip::conj(U[j][k]) * V[i][k];

                      }
                  }
      }
  }

template <class Matrix1, 
                  class Matrix2>
inline
void svd_inv(const Matrix1 & A,
                     Matrix2& Ainv)
{
  typedef typename Matrix1::value_type value_type;
    typedef typename slip::lin_alg_traits<typename Matrix1::value_type>::value_type real;                       
  

    slip::Vector<real>  S(A.cols());
    slip::Array2d<value_type> U(A.rows(),A.cols());
    slip::Array2d<value_type> V(A.cols(),A.cols());
    
    //-----------------------------------------------
    // svd decomposition.
    //-----------------------------------------------
    slip::svd(A,U,S.begin(),S.end(),V);
    //-----------------------------------------------
    // svd inverse
    //-----------------------------------------------
    slip::svd_inv(U,S.begin(),S.end(),V,Ainv);
}

template <class Matrix1, 
                  class Matrix2>
inline
void pinv(const Matrix1 & A,
                  Matrix2& Ainv)
{
  slip::svd_inv(A,Ainv);
}
/* @} */
} // end slip

#endif //SLIP_LINEAR_ALGEBRA_SVD_HPP