Changeset 227


Ignore:
Timestamp:
Feb 1, 2005, 1:52:27 AM (18 years ago)
Author:
Jari Häkkinen
Message:

Started reimplementation of the vector class.

Location:
trunk
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/PCA.cc

    r189 r227  
    3535  // Single value decompose the data matrix
    3636  std::auto_ptr<cpptools::SVD> pSVD( new cpptools::SVD( A_center ) );
    37   pSVD->process();
    38   gslapi::matrix U = pSVD->get_u();
    39   gslapi::matrix V = pSVD->get_v();
     37  pSVD->decompose();
     38  gslapi::matrix U = pSVD->U();
     39  gslapi::matrix V = pSVD->V();
    4040
    4141  // Read the eigenvectors and eigenvalues
    4242  eigenvectors_ = U.transpose();
    43   eigenvalues_ = pSVD->get_s_vec();
     43  eigenvalues_ = pSVD->s();
    4444
    4545  // T
    … …  
    8383  // Single value decompose the data matrix
    8484  std::auto_ptr<cpptools::SVD> pSVD( new cpptools::SVD( A_center ) );
    85   pSVD->process();
    86   gslapi::matrix U = pSVD->get_u();
    87   gslapi::matrix V = pSVD->get_v();
     85  pSVD->decompose();
     86  gslapi::matrix U = pSVD->U();
     87  gslapi::matrix V = pSVD->V();
    8888
    8989  // Read the eigenvectors and eigenvalues
    9090  eigenvectors_ = V.transpose();//U.transpose();
    91   eigenvalues_ = pSVD->get_s_vec();
     91  eigenvalues_ = pSVD->s();
    9292
    9393  // Transform back when done with SVD!

  • trunk/src/SVD.cc

    r42 r227  
    11// $Id$
    2 
    3 #include <cstdlib>
    4 #include <cmath>
    5 #include <iostream>
    6 
    7 #include <gsl/gsl_linalg.h>
    82
    93#include "SVD.h"
    … …  
    115
    126namespace theplu {
    13 namespace cpptools { 
    14 
    15 //Constructor that can only be used for test-method
    16 SVD::SVD() : loaded_( false ), solver_( false )
    17 {}
    18 
    19 //Constructor initializing the svd object with matrix A to be decomposed
    20 SVD::SVD( const gslapi::matrix& Ain ) : A_( Ain )
    21 {
    22   //Get dimensions
    23   ////////////////
    24   size_t n = Ain.columns();
    25 
    26   //Instantiate matrix with correct dimenstions
    27   //and assign all elements the value zero
    28   /////////////////////////////////////////////
    29   V_ = gslapi::matrix( n, n, true );
    30  
    31   //Assign vector s_ with correct dimenstion
    32   /////////////////////////////////////////////////
    33   s_ = gslapi::vector( n );
    34  
    35   loaded_ = true;
    36   solver_ = false;
    37 }
     7namespace cpptools {
    388
    399
    40 //Constructor initializing the svd object for Solver
    41 //Solver requires an extra-vector b:
    42 //Ax = b
    43 //The x vector will be reached using the get_x() function.
    44 SVD::SVD( const gslapi::matrix& Ain, const gslapi::vector& b ) : A_( Ain ), b_( b )
    45 {
    46   //Get dimensions
    47   ////////////////
    48   size_t n = Ain.columns();
    49  
    50   //Instantiate matrix with correct dimenstions
    51   //and assign all elements the value zero
    52   /////////////////////////////////////////////
    53     V_ = gslapi::matrix( n, n, true );
    54 
    55 
    56   //Assign vector s_ with correct dimenstion
    57   /////////////////////////////////////////////////
    58   s_ = gslapi::vector( n );
    59  
    60  
    61   //An extra vector holding the answer to equation Ax = b is required
    62   ///////////////////////////////////////////////////////////////////
    63   x_ = gslapi::vector( n );
    64 
    65 
    66   //Also set solver to true
    67   /////////////////////////
    68   solver_ = true;
    69   loaded_ = true;
    70 }
    71 
    72 
    73 //Destructor
    74 SVD::~SVD()
    75 {
    76 }
     10  SVD::SVD(const gslapi::matrix& Ain)
     11    : U_(Ain), V_(Ain.columns(),Ain.columns()), s_(Ain.columns())
     12  {
     13  }
    7714
    7815
    7916
    80 double SVD::norm(gslapi::matrix& A) const
    81 {
    82   double sum=0.0;
    83   for (size_t i=0; i<A.rows(); ++i)
    84     for (size_t j=0; j<A.columns(); ++j)
    85       sum += A(i,j)*A(i,j);
    86   return sqrt(sum);
    87 }
     17  SVD::~SVD(void)
     18  {
     19  }
    8820
    8921
    9022
    91 //Method to process calculation
    92 bool SVD::process( SVDtypes stype )
    93 {
    94   if( !loaded_ )
    95     {
    96      std::cerr << "Use proper constructor, matrix A not loaded!\n";
    97      return false;
     23  int SVD::decompose(SVDalgorithm algo)
     24  {
     25    switch (algo) {
     26        case GolubReinsch:
     27          return golub_reinsch();
     28        case ModifiedGolubReinsch:
     29          return modified_golub_reinsch();
     30        case Jacobi:
     31          return jacobi();
    9832    }
    99 
    100   switch( stype )
    101     {
    102 
    103     case Unmodified:
    104       if( !process_default() ) exit( -1 );
    105       break;
    106 
    107     case Modified:
    108       if( !process_modified() ) exit( -1 );
    109       break;
    110 
    111     case Jacobi:
    112       if( !process_jacobi() ) exit( -1 );
    113       break;
    114 
    115     case Solver:
    116       if( !process_solver() ) exit( -1 );
    117       break;
    118 
    119     default:
    120       std::cerr << "Programmer try to use method not implemented!\n";
    121       return false;
    122       break;
    123 
    124     }
    125   return true;
    126 }
    127 
    128 
    129 //Default-method for SVD
    130 bool SVD::process_default()
    131 {
    132   //Perform SVD using GSL library
    133   //An additional workvector is
    134   //required here
    135   ///////////////////////////////
    136   gslapi::vector w( A_.columns() );
    137   if( gsl_linalg_SV_decomp( A_.gsl_matrix_pointer(), V_.gsl_matrix_pointer(),
    138           s_.gsl_vector_pointer(), w.gsl_vector_pointer()
    139           ) != 0 )
    140     {
    141      std::cerr << "failed to run gsl_linalg_SV_decomp!\n";
    142      return false; 
    143     }
    144 
    145   return true;
    146 }
    147 
    148 
    149 //Modified SVD method
    150 //This method is faster when M>>N
    151 bool SVD::process_modified()
    152 {
    153   //Perform SVD using GSL library
    154   //An additional workmatrix X is
    155   //required here
    156   ///////////////////////////////
    157   size_t n = A_.rows();
    158   gslapi::vector w( n );
    159   gslapi::matrix X( n, n, true );
    160 
    161   if( gsl_linalg_SV_decomp_mod( A_.gsl_matrix_pointer(), X.gsl_matrix_pointer(),
    162         V_.gsl_matrix_pointer(), s_.gsl_vector_pointer(),
    163         w.gsl_vector_pointer()
    164         ) != 0 )
    165     {
    166      std::cerr << "failed to run gsl_linalg_SV_decomp_mod!\n";
    167      return false; 
    168     }
    169 
    170   return true;
    171 }
    172 
    173 
    174 //Default-method for SVD
    175 //Computes singular values to higer
    176 //accuracy
    177 bool SVD::process_jacobi()
    178 {
    179   //Perform SVD using GSL library
    180   ///////////////////////////////
    181   if( gsl_linalg_SV_decomp_jacobi( A_.gsl_matrix_pointer(), V_.gsl_matrix_pointer(),
    182            s_.gsl_vector_pointer()
    183            ) != 0 )
    184     {
    185      std::cerr << "failed to run gsl_linalg_SV_decomp_jacobi!\n";
    186      return false; 
    187     }
    188 
    189   return true;
    190 }
    191 
    192 
    193 //Solver solves the system Ax=b
    194 //using SVD (process_default).
    195 bool SVD::process_solver()
    196 {
    197   //Perform SVD using GSL library
    198   ///////////////////////////////
    199   this->process_default();
    200  
    201   //Solve ...
    202   if( gsl_linalg_SV_solve( A_.gsl_matrix_pointer(), V_.gsl_matrix_pointer(),
    203          s_.gsl_vector_pointer(), b_.gsl_vector_pointer(),
    204          x_.gsl_vector_pointer()
    205          ) != 0 )
    206     {
    207      std::cerr << "failed to run gsl_linalg_SV_solve!\n";
    208      return false; 
    209     }
    210  
    211  
    212   return true;
    213 }
     33    // the program should never end up here, return values should be
     34    // something different from normal GSL return values, or maybe
     35    // throw an exception.
     36    return 0;
     37  }
    21438
    21539
    21640
    217 //This method will create a diagonal matrix of vector s_
    218 gslapi::matrix SVD::get_s() const
    219 {
    220   size_t n = A_.rows();
    221  
    222   gslapi::matrix S( n, n, true );
    223 
    224   for( size_t i = 0; i < n; ++i )
    225      S(i,i)= s_(i);
    226 
    227   return S;
    228 }
    229 
    230 //This function will run the SVD function and print
    231 //matrices in matlab format so that you can verify
    232 //the answer.
    233 bool SVD::test()
    234 {
    235   using namespace std;
     41  int SVD::golub_reinsch(void)
     42  {
     43    gslapi::vector w(U_.columns());
     44    return gsl_linalg_SV_decomp(U_.gsl_matrix_pointer(),
     45                                V_.gsl_matrix_pointer(),
     46                                s_.TEMP_gsl_vector_pointer(),
     47                                w.TEMP_gsl_vector_pointer());
     48  }
    23649
    23750
    238   //Initialize GNU random stuff
    239   ///////////////////////////////////
    240     srand( time( 0 ) );
    24151
    242    
    243   //Initialise random test-matrix
    244   ///////////////////////////////
    245     const size_t DIM = 3;
    246     gslapi::matrix A( DIM, DIM );
    247     for( size_t i = 0; i < DIM; ++i )
    248       {
    249        for( size_t j = 0; j < DIM; ++j )
    250    {
    251      // Jari, change below random numbers to proper RNG
    252     A(i,j)= 10 * (double) rand() / RAND_MAX; //Use random numbers between 0 and 10
    253    }
    254       }
     52  int SVD::modified_golub_reinsch(void)
     53  {
     54    gslapi::vector w(U_.columns());
     55    gslapi::matrix X(U_.columns(),U_.columns());
     56    return gsl_linalg_SV_decomp_mod(U_.gsl_matrix_pointer(),
     57                                    X.gsl_matrix_pointer(),
     58                                    V_.gsl_matrix_pointer(),
     59                                    s_.TEMP_gsl_vector_pointer(),
     60                                    w.TEMP_gsl_vector_pointer());
     61  }
    25562
    25663
    257     //Set data
    258     //////////
    259     A_ = A;
    260     size_t n = A_.rows();
    261     size_t m = A_.columns();
    262     n = m = DIM;
    263     loaded_ = true;
    264 
    265 
    266     //Set work matrices
    267     ///////////////////
    268     V_ = gslapi::matrix( n, n, true );
    269 
    270     //Assign vector s_ with correct dimensions
    271     /////////////////////////////////////////////////
    272     s_ = gslapi::vector( n );
    273 
    274     double error = 1;
    275     gslapi::matrix E, S;
    276 
    277 
    278     //Print matrix to screen
    279     ////////////////////////
    280     cout << "A = " << endl;
    281     cout << A << endl;
    282 
    283     cout.setf( ios::scientific );
    284 
    285 
    286     //Run test 1: Unmodified method
    287     ///////////////////////////////
    288     this->process();
    289     S = this->get_s();
    290 
    291     E = A - A_*S*V_.transpose();
    292    
    293 
    294     error = sqrt(norm(E));
    295     cout << "error(Unmodified method) = " << error << endl;
    296     cout << "to compare with " << MAXTOL << endl;
    297     if( error < MAXTOL ) { cout << "test ok!!!\n"; }
    298     else { cout << "test failed!\n"; return false; }
    299 
    300  
    301     //Run test 2: Modified method
    302     ///////////////////////////////
    303     A_ = A;
    304     this->process( Modified );
    305     S = this->get_s();
    306     E = A - A_*S*V_.transpose();
    307     error = sqrt(norm(E));
    308     cout << "error(Modified method) = " << error << endl;
    309     cout << "to compare with " << MAXTOL << endl;
    310     if( error < MAXTOL ) { cout << "test ok!!!\n"; }
    311     else { cout << "test failed!\n"; return false; }
    312 
    313 
    314     //Run test 2: Modified method
    315     ///////////////////////////////
    316     A_ = A;
    317     this->process( Jacobi );
    318     S = this->get_s();
    319     E = A - A_*S*V_.transpose();
    320     error = sqrt(norm(E));
    321     cout << "error(Jacobi method) = " << error << endl;
    322     cout << "to compare with " << MAXTOL << endl;
    323     if( error < MAXTOL ) { cout << "test ok!!!\n"; }
    324     else { cout << "test failed!\n"; return false; }
    325    
    326     return true; //Test with good outcome!
    327 }
    328 
    32964}} // of namespace cpptools and namespace theplu

  • trunk/src/SVD.h

    r43 r227  
    11// $Id$
    22
    3 #ifndef _theplu_cpptools_svd_ 
    4 #define _theplu_cpptools_svd_ 
     3#ifndef _theplu_cpptools_svd_
     4#define _theplu_cpptools_svd_
    55
    66#include "matrix.h"
    77#include "vector.h"
     8
     9#include <gsl/gsl_linalg.h>
    810
    911namespace theplu {
    … …  
    1315/// namespace.
    1416///
    15 namespace cpptools { 
     17namespace cpptools {
    1618
    17   /**
    18      Class encapsulating methods for Singular Value Decomposition.
    19      \n\n
     19  // Jari check that interface is complete
     20
     21  /**
     22     Class encapsulating GSL methods for singular value decomposition,
     23     SVD.
     24
     25     A = U S V' = (MxN)(NxN)(NxN) = (MxN)\n   
     26
    2027     A = Matrix to be decomposed, size MxN\n
    2128     U = Orthogonal matrix, size MxN\n
    2229     S = Diagonal matrix of singular values, size NxN\n
    2330     V = Orthogonal matrix, size NxN\n
    24      A = U S V' = (MxN)(NxN)(NxN) = (MxN)\n   
    25 
    26      Shouldn't there be a reference to the algortithm here?
    27 
    28      @version 1.0 $Revision$ $Date$
    2931  */
    3032
    3133  class SVD
    3234  {
    33     static const double MAXTOL = 1E-6;
     35  public:
    3436
    35     enum SVDtypes {
    36       Unmodified,
    37       Modified,
    38       Jacobi,
    39       Solver
     37    ///
     38    /// A number of SVD algorithms are implemented in GSL. They have
     39    /// their strengths and weaknesses, check the GSL documentation.
     40    ///
     41    /// There are restrictions on the matrix dimensions depending on
     42    /// which SVD algorithm is used. From the GSL's SVD source code
     43    /// one finds that the Golub-Reinsch algorithm implementation will
     44    /// not work on matrices with fewer rows than columns, the same is
     45    /// also true for the modified Golub-Reinsch algorithm.
     46    ///
     47    /// @see GSL's SVD documentation.
     48    ///
     49    enum SVDalgorithm {
     50      GolubReinsch,
     51      ModifiedGolubReinsch,
     52      Jacobi
    4053    };
    4154
    42   public:
    43     /**
    44        The default constructor should only be used for testing!!!
    45     */
    46     SVD(void);
     55    ///
     56    /// Constructs an SVD object using the matrix A as only input. The
     57    /// input matrix is copied for further use in the object.
     58    ///
     59    SVD(const gslapi::matrix& );
    4760
    48    
    49     /**
    50        Constructs an SVD object using the matrix A as only
    51        input. The input matrix is copied for further use in the
    52        object.
    53     */
    54     SVD( const gslapi::matrix& );
    55    
    56     /**
    57        Constructor initializing the SVD object for Solver. 
    58        Solver requires an extra-vector paramter \a b: \f$Ax = b\f$
    59        The \f$x\f$ vector will be reached using the get_x()
    60        function.
    61        @see get_x() function.
    62     */
     61    ~SVD(void);
    6362
    64     SVD( const gslapi::matrix&, const gslapi::vector& b);
    65     ~SVD();
     63    ///
     64    /// This function will perform SVD with the method specified by \a
     65    /// algo.
     66    ///
     67    /// @return Whatever GSL returns.
     68    ///
     69    int decompose(SVDalgorithm algo=GolubReinsch);
    6670
    67     /**
    68        This function will run the method specified by the SVDtypes. Precently
    69        there are 3 SVD methods implemented and one solver. These are: \n\n
    70        Unmodified: (see GSL documentation about gsl_linalg_SV_decomp function)\n
    71        Modified: This method is faster when M>>N. (see GSL documentation
    72        about gsl_linalg_SV_decomp_mod function)\n
    73        Jacobi: Computes singular values to higer accuracy than default method
    74        (see GSL documentation about gsl_linalg_SV_decomp_jacobi.\n
    75        Solver: Using Unmodified to perform SVD and then solves the system \f$Ax=b\f$
    76     */
    77     bool process( SVDtypes stype = Unmodified );
     71    ///
     72    /// Access to the s vector.
     73    ///
     74    /// @return A copy of the s vector.
     75    ///
     76    /// @note If decompose() has not been run the outcome of the call
     77    /// is undefined.
     78    ///
     79    gslapi::vector s(void) const { return s_; }
    7880
    79     /**
    80        Function will return the matrix U. Require that process() has been
    81        executed.
    82     */
    83     gslapi::matrix get_u() const { return A_; }
     81    ///
     82    /// Solve the system \f$Ax=b\f$ using the decomposition of A.
     83    ///
     84    /// @note If decompose() has not been run the outcome of the call
     85    /// is undefined.
     86    ///
     87    /// @return Whatever GSL returns.
     88    ///
     89    inline int solve(gslapi::vector b, gslapi::vector x)
     90      { return gsl_linalg_SV_solve(U_.gsl_matrix_pointer(),
     91                                   V_.gsl_matrix_pointer(),
     92                                   s_.gsl_vector_pointer(),
     93                                   b.gsl_vector_pointer(),
     94                                   x.TEMP_gsl_vector_pointer());
     95      }
    8496
    85     /**
    86        Function will return the matrix V. Require that process() has been
    87        executed.
    88     */
    89     gslapi::matrix get_v() const { return V_; }
     97    ///
     98    /// Access to the U matrix.
     99    ///
     100    /// @return A copy of the U matrix.
     101    ///
     102    /// @note If decompose() has not been run the outcome of the call
     103    /// is undefined.
     104    ///
     105    gslapi::matrix U(void) const { return U_; }
    90106
     107    ///
     108    /// Access to the V matrix.
     109    ///
     110    /// @return A copy of the V matrix.
     111    ///
     112    /// @note If decompose() has not been run the outcome of the call
     113    /// is undefined.
     114    ///
     115    gslapi::matrix V(void) const { return V_; }
    91116
    92     /**
    93        Function will return the vector x containing the solution
    94        \f$ x=A^{-1}b \f$. Require that process( Solver ) has been
    95        executed. 
    96        @see Read about Solver in the bool process() function.
    97     */
     117  private:
     118    inline int jacobi(void)
     119      { return gsl_linalg_SV_decomp_jacobi(U_.gsl_matrix_pointer(),
     120                                           V_.gsl_matrix_pointer(),
     121                                           s_.TEMP_gsl_vector_pointer());
     122      }
     123    int golub_reinsch(void);
     124    int modified_golub_reinsch(void);
    98125
    99     gslapi::vector get_x() const { return x_; }
    100 
    101     /**
    102        Function returning diagonal matrix S. Require that process() has been
    103        executed.
    104     */
    105     gslapi::matrix get_s() const;
    106 
    107 
    108     /**
    109        Function returning diagonal matrix S in vector form.
    110     */
    111     gslapi::vector get_s_vec() const { return s_; }
    112 
    113 
    114     /**
    115        This method will execute Default, Unmodified and Jacobi methods one at a time
    116        and calculate the error, \f$ e = \left\Vert A - USV^{T} \right\Vert_{2} \f$.
    117        The method will return "true" if \f$e < 10^{-6}\f$ else false. A test program
    118        is available that executes this method, c++_tools/test/svd_test. No test is
    119        performed for solver but the aim is to add one in the future.
    120     */
    121     bool test();
    122    
    123   private:
    124 
    125     gslapi::matrix A_, V_;   
    126     gslapi::vector s_, b_, x_;
    127     bool loaded_, solver_;
    128 
    129     //SVD methods
    130     bool process_default();
    131     bool process_modified();
    132     bool process_jacobi();
    133     bool process_solver();
    134 
    135     double norm(gslapi::matrix& A) const; // move into class gslapi::matrix?
     126    gslapi::matrix U_, V_;
     127    gslapi::vector s_;
    136128  }; 
    137129

  • trunk/src/SVM.cc

    r197 r227  
    119119    }
    120120 
    121     gslapi::vector  E = ( kernel_train*target_train.mul_elements(alpha_train)
     121    gslapi::vector  E = ( kernel_train*target_train.mul(alpha_train)
    122122                          - target_train ) ;
    123123
    … …  
    184184      alpha_train[index.second] = alpha_new2;
    185185   
    186       E = kernel_train*target_train.mul_elements(alpha_train) - target_train;
     186      E = kernel_train*target_train.mul(alpha_train) - target_train;
    187187      gslapi::vector ones(alpha_train.size(),1);
    188188   

  • trunk/src/SVM.h

    r188 r227  
    6262    ///
    6363    inline theplu::gslapi::vector output(void)
    64     {return kernel_.get() * alpha_.mul_elements(target_)+
     64    {return kernel_.get() * alpha_.mul(target_)+
    6565       theplu::gslapi::vector(alpha_.size(),bias_);}
    6666

  • trunk/src/vector.cc

    r132 r227  
    2020
    2121  vector::vector(void)
    22     : v_(NULL)
    23   {
    24   }
    25 
    26 
    27 
    28   vector::vector(const size_t n, const double init_value)
     22    : v_(NULL), view_(NULL)
     23  {
     24  }
     25
     26
     27
     28  vector::vector(const size_t n,const double init_value)
     29    : view_(NULL)
    2930  {
    3031    v_ = gsl_vector_alloc(n);
    … …  
    3435
    3536
    36   vector::vector(const vector& other, const std::vector<size_t>& index)
    37   {
    38     if (index.size()){
    39       v_ = gsl_vector_alloc(index.size());
    40       for (size_t i=0; i<index.size(); i++)
    41         gsl_vector_set(v_, i, other(index[i]));
    42     }
    43     else
    44       v_ = other.gsl_vector_copy();
     37  vector::vector(const vector& other)
     38    : view_(NULL)
     39  {
     40    v_ = other.TEMP_gsl_vector_copy();
    4541  }
    4642
    … …  
    4844
    4945  vector::vector(gsl_vector* vec)
    50     : v_(vec)
     46    : v_(vec), view_(NULL)
    5147  {
    5248  }
    … …  
    5551
    5652  vector::vector(std::istream& is)
     53    : view_(NULL)
    5754  {
    5855    using namespace std;
    … …  
    132129      v_=NULL;
    133130    }
    134   }
    135 
    136 
    137 
    138   gsl_vector* vector::gsl_vector_copy(void) const
     131    delete view_;
     132  }
     133
     134
     135
     136  gsl_vector* vector::TEMP_gsl_vector_copy(void) const
    139137  {
    140138    gsl_vector* vec = gsl_vector_alloc(size());
    … …  
    145143
    146144
     145  std::pair<double,double> vector::minmax(void) const
     146  {
     147    double min, max;
     148    gsl_vector_minmax(v_,&min,&max);
     149    return std::pair<double,double>(min,max);
     150  }
     151
    147152
    148153
    149154  std::pair<size_t,size_t> vector::minmax_index(void) const
    150155  {
    151     size_t min_index=0;
    152     size_t max_index=0;
    153     void gsl_vector_minmax_index (const gsl_vector * v_, size_t *
    154                                   min_index, size_t * max_index);
    155     return std::pair<size_t,size_t>(min_index, max_index);
    156   }
    157 
    158 
    159 
    160 
     156    size_t min_index, max_index;
     157    gsl_vector_minmax_index(v_,&min_index,&max_index);
     158    return std::pair<size_t,size_t>(min_index,max_index);
     159  }
     160
     161
     162
     163/* Jari, remove this section
    161164  std::pair<size_t,size_t> vector::minmax_index(const std::vector<size_t>& subset ) const
    162165  {
    … …  
    171174    return std::pair<size_t,size_t>(min_index, max_index);
    172175  }
    173 
    174 
    175 
    176   vector vector::mul_elements( const vector& other ) const
    177   {
    178     vector res( *this );
    179     gsl_vector_mul( res.v_, other.v_ );
    180     return res;
    181   }
     176*/
    182177
    183178  void vector::shuffle(void) const
    … …  
    196191  {
    197192    double sum = 0;
    198     for (u_int i=0; i<size(); i++)
     193    for (size_t i=0; i<size(); i++)
    199194      sum += gsl_vector_get( v_, i );
    200195    return( sum );
    … …  
    246241      if ( v_ )
    247242        gsl_vector_free( v_ );
    248       v_ = other.gsl_vector_copy();
     243      v_ = other.TEMP_gsl_vector_copy();
    249244    }
    250245    return *this;

  • trunk/src/vector.h

    r142 r227  
    2828
    2929  ///
    30   /// This is the C++ tools interface to GSL vector.
    31   ///
    32 
     30  /// This is the C++ tools interface to GSL vector. 'double' is the
     31  /// only type supported, maybe we should add a 'complex' type as
     32  /// well in the future.
     33  ///
     34  /// \par[File streams] Reading and writing vectors to file streams
     35  /// are of course supported. These are implemented without using GSL
     36  /// functionality, and thus binary read and write to streams are not
     37  /// supported.
     38  ///
     39  /// \par[Vector views] GSL vector views are supported and these are
     40  /// disguised as ordinary gslapi::vectors. A support function is
     41  /// added, gslapi::vector::is_view(), that can be used to check if
     42  /// the vector object is a view. Note that view vectors do not own
     43  /// the undelying data, and is not valid if the original vector is
     44  /// deallocated.
     45  ///
    3346  class vector
    3447  {
    … …  
    4760
    4861    ///
    49     /// The copy constructor, creating a new sub-vector defined by \a
    50     /// index (default is to copy the whole vector).
    51     ///
    52     vector(const vector&,
    53            const std::vector<size_t>& index = std::vector<size_t>());
     62    /// The copy constructor.
     63    ///
     64    /// @note If the object to be copied is a vector view, the values
     65    /// of the view will be copied, i.e. the view is not copied.
     66    ///
     67    vector(const vector&);
    5468
    5569    ///
    … …  
    7387
    7488    ///
    75     /// Create a new copy of the internal GSL vector.
    76     ///
    77     /// Necessary memory for the new GSL vector is allocated and the
    78     /// caller is responsible for freeing the allocated memory.
    79     ///
    80     /// @return A pointer to a copy of the internal GSL vector.
    81     ///
    82     gsl_vector* gsl_vector_copy(void) const;
     89    /// Vector addition, \f$this_i = this_i + other_i \; \forall i\f$.
     90    ///
     91    /// @return .GSL_SUCCESS on normal exit.
     92    ///
     93    // Jari, doxygen group as Vector operators
     94    inline int add(const vector& other) { return gsl_vector_add(v_,other.v_); }
     95
     96    ///
     97    /// Add a constant to a vector, \f$this_i = this_i + term \;
     98    /// \forall i\f$.
     99    ///
     100    /// @return .GSL_SUCCESS on normal exit.
     101    ///
     102    // Jari, doxygen group as Vector operators
     103    inline int
     104    add_constant(double term) { return gsl_vector_add_constant(v_,term); }
     105
     106    ///
     107    /// This function performs element-wise division, \f$this_i =
     108    /// this_i/other_i \; \forall i\f$.
     109    ///
     110    /// @return .GSL_SUCCESS on normal exit.
     111    ///
     112    // Jari, doxygen group as Vector operators
     113    inline int div(const vector& other) { return gsl_vector_div(v_,other.v_); }
     114
     115    ///
     116    /// @return True if all elements in the vector is zero, false
     117    /// othwerwise;
     118    ///
     119    inline bool isnull(void) const { return gsl_vector_isnull(v_); }
     120
     121    ///
     122    /// Check if the vector object is a view (sub vector) to another
     123    /// vector.
     124    ///
     125    /// @return True if the object is a view, false othwerwise;
     126    ///
     127    inline bool isview(void) const { return view_; }
    83128
    84129    ///
    … …  
    90135    /// @return A pointer to the internal GSL vector,
    91136    ///
    92     inline gsl_vector* gsl_vector_pointer(void) { return v_; }
    93 
    94     ///This function returns the indices of the minimum and maximum values in
    95     ///the vector, storing them in imin and imax. When
    96     ///there are several equal minimum or maximum elements then the lowest
    97     ///indices are returned. @return Index corresponding to the smallest
    98     /// and largest value.
    99     ///
     137    inline gsl_vector* TEMP_gsl_vector_pointer(void) { return v_; }
     138
     139    ///
     140    /// @return The maximum value of the vector.
     141    ///
     142    // Jari, doxygen group as Finding maximum and minimum elements
     143    inline double max(void) const { return gsl_vector_max(v_); }
     144
     145    ///
     146    /// @return The element index to the maximum value of the
     147    /// vector. The lowest index has precedence.
     148    ///
     149    // Jari, doxygen group as Finding maximum and minimum elements
     150    inline size_t max_index(void) const { return gsl_vector_max_index(v_); }
     151
     152    ///
     153    /// @return The minimum value of the vector.
     154    ///
     155    // Jari, doxygen group as Finding maximum and minimum elements
     156    inline double min(void) const { return gsl_vector_min(v_); }
     157
     158    ///
     159    /// @return The element index to the minimum value of the
     160    /// vector. The lowest index has precedence.
     161    ///
     162    // Jari, doxygen group as Finding maximum and minimum elements
     163    inline size_t min_index(void) const { return gsl_vector_min_index(v_); }
     164
     165    ///
     166    /// @return The minimum and maximum values of the vector, as the
     167    /// \a first and \a second member of the returned \a pair,
     168    /// respectively.
     169    ///
     170    // Jari, doxygen group as Finding maximum and minimum elements
     171    std::pair<double,double> vector::minmax(void) const;
     172
     173    ///
     174    /// @return The indecies to the minimum and maximum values of the
     175    /// vector, as the \a first and \a second member of the returned
     176    /// \a pair, respectively. The lowest index has precedence.
     177    ///
     178    // Jari, doxygen group as Finding maximum and minimum elements
    100179    std::pair<size_t,size_t> vector::minmax_index(void) const;
    101180
    … …  
    110189    ///
    111190    std::pair<size_t,size_t>
    112     vector::minmax_index(const std::vector<size_t>& subset ) const;
    113 
    114     ///
    115     /// This function multiplies all elements between two vectors and
    116     /// returns a third vector containing the result, \f$a_i = b_i *
    117     /// c_i \; \forall i\f$.
    118     ///
    119     /// @return The result vector.
    120     ///
    121     vector vector::mul_elements( const vector& other ) const;
     191    vector::TEMP_minmax_index(const std::vector<size_t>& subset ) const;
     192
     193    ///
     194    /// This function performs element-wise multiplication, \f$this_i =
     195    /// this_i * other_i \; \forall i\f$.
     196    ///
     197    /// @return .GSL_SUCCESS on normal exit.
     198    ///
     199    // Jari, doxygen group as Vector operators
     200    inline int mul(const vector& other) { return gsl_vector_mul(v_,other.v_); }
     201
     202    ///
     203    /// Reverse the order of elements in the vector.
     204    ///
     205    /// @return .GSL_SUCCESS on normal exit.
     206    ///
     207    // Jari, doxygen group as Exchanging elements
     208    inline int reverse(void) { return gsl_vector_reverse(v_);}
     209
     210    ///
     211    /// Rescale vector, \f$this_i = this_i * factor \; \forall i\f$.
     212    ///
     213    /// @return .GSL_SUCCESS on normal exit.
     214    ///
     215    // Jari, doxygen group as Vector operators
     216    inline int scale(double factor) { return gsl_vector_scale(v_,factor); }
    122217
    123218    ///
    124219    /// Set all elements to \a value.
    125220    ///
     221    // Jari, doxygen group as Initializing vector elements
    126222    inline void set_all(const double& value) { gsl_vector_set_all(v_,value); }
    127223
    … …  
    131227    /// one.
    132228    ///
     229    // Jari, doxygen group as Initializing vector elements
    133230    inline void set_basis(const size_t i) { gsl_vector_set_basis(v_,i); }
    134231   
    … …  
    136233    /// Set all elements to zero.
    137234    ///
     235    // Jari, doxygen group as Initializing vector elements
    138236    inline void set_zero(void) { gsl_vector_set_zero(v_); }
    139237
    … …  
    150248
    151249    ///
     250    /// Vector subtraction, \f$this_i = this_i - other_i \; \forall i\f$.
     251    ///
     252    /// @return .GSL_SUCCESS on normal exit.
     253    ///
     254    // Jari, doxygen group as Vector operators
     255    inline int sub(const vector& other) { return gsl_vector_sub(v_,other.v_); }
     256
     257    ///
    152258    /// Calculate the sum of all vector elements.
    153259    ///
    154260    /// @return The sum.
    155261    ///
     262    // Jari, doxygen group as "Extends GSL".
    156263    double vector::sum(void) const;
    157264
    158265    ///
     266    /// Swap vector elements by copying. The two vectors must have the
     267    /// same length.
     268    ///
     269    /// @return .GSL_SUCCESS on normal exit.
     270    ///
     271    inline int swap(vector& other) { return gsl_vector_swap(v_,other.v_); }
     272
     273    ///
     274    /// Exchange elements \a i and \a j.
     275    ///
     276    /// @return .GSL_SUCCESS on normal exit.
     277    ///
     278    // Jari, doxygen group as Exchanging elements
     279    inline int
     280    swap_elements(size_t i,size_t j) { return gsl_vector_swap_elements(v_,i,j);}
     281
     282    ///
    159283    /// Element access operator.
    160284    ///
    161285    /// @return Reference to element \a i.
    162286    ///
     287    // Jari, doxygen group as Accessing vector elements
    163288    inline double& operator()(size_t i) { return *gsl_vector_ptr(v_,i); }
    164289
    … …  
    168293    /// @return The value of element \a i.
    169294    ///
     295    // Jari, doxygen group as Accessing vector elements
    170296    inline const double& operator()(size_t i) const
    171297      { return *gsl_vector_const_ptr(v_,i); }
    … …  
    176302    /// @return Reference to element \a i.
    177303    ///
     304    // Jari, doxygen group as Accessing vector elements
    178305    inline double& operator[](size_t i) { return *gsl_vector_ptr(v_,i); }
    179306
    … …  
    183310    /// @return The value of element \a i.
    184311    ///
     312    // Jari, doxygen group as Accessing vector elements
    185313    inline const double& operator[](size_t i) const
    186314      { return *gsl_vector_const_ptr(v_,i); }
    … …  
    217345
    218346    ///
    219     /// Multiply and assign operator.
     347    /// Multiply with scalar and assign operator.
    220348    ///
    221349    vector& operator*=(const double d) { gsl_vector_scale(v_,d); return *this; }
    222350
    223351    ///
    224     /// Multiply operator.
    225     ///
    226     vector& operator*(const double d) { gsl_vector_scale(v_,d); return *this; }
     352    /// Vector with scalar multiplication.
     353    ///
     354    /// \note This operator is not implemented!
     355    ///
     356    vector operator*(const double d) const;
    227357
    228358
    229359  private:
    230360
     361    ///
     362    /// Create a new copy of the internal GSL vector.
     363    ///
     364    /// Necessary memory for the new GSL vector is allocated and the
     365    /// caller is responsible for freeing the allocated memory.
     366    ///
     367    /// @return A pointer to a copy of the internal GSL vector.
     368    ///
     369    gsl_vector* TEMP_gsl_vector_copy(void) const;
     370
    231371    gsl_vector* v_;
    232 
     372    gsl_vector_view* view_;
    233373  };
    234374

  • trunk/test/test_svd.cc

    r42 r227  
    11// $Id$
     2//
     3// Testing for a succesful reconstruction of a decomposed random
     4// matrix, aswell that U and V are orthogonal.
    25
    3 //Test program for SVD class.
    46
    57#include "SVD.h"
     8#include "matrix.h"
     9#include "random_singleton.h"
     10#include "vector.h"
    611
    7 int main()
     12using namespace std;
     13using namespace theplu::cpptools;
     14using namespace theplu::gslapi;
     15
     16
     17
     18double this_norm(const matrix& A)
    819{
    9   theplu::cpptools::SVD* mysvd = new theplu::cpptools::SVD;
    10   if( mysvd->test() ) { delete mysvd; return 0; }
    11   else { delete mysvd; return -1; }
     20  double sum=0.0;
     21  for (size_t i=0; i<A.rows(); ++i)
     22    for (size_t j=0; j<A.columns(); ++j)
     23      sum += A(i,j)*A(i,j);
     24  return sum;
    1225}
     26
     27
     28
     29bool test(size_t m, size_t n, SVD::SVDalgorithm algo) {
     30
     31  // accepted error, Jari: should be picked up from GSL
     32  double MAXTOL=1e-13;
     33  random_singleton* rng=random_singleton::get_instance(-1);
     34
     35  // initialise a random test-matrix
     36  matrix A(m,n);
     37  for (size_t i=0; i<m; ++i)
     38    for(size_t j=0; j<n; ++j)
     39      A(i,j)=1000*rng->get_uniform_double();
     40
     41  SVD svd(A);
     42  svd.decompose(algo);
     43  theplu::gslapi::vector s(svd.s());
     44  matrix S(s.size(),s.size());
     45  for (size_t i=0; i<s.size(); ++i)
     46    S(i,i)=s[i];
     47  matrix V=svd.V();
     48  matrix U=svd.U();
     49  double error = this_norm(A-U*S*V.transpose());
     50  bool testerror=false;
     51  if (error>MAXTOL) {
     52    cerr << "test_svd: FAILED, algorithm " << algo << " recontruction error ("
     53         << error << ") > tolerance (" << MAXTOL << "), matrix dimension ("
     54         << m << ',' << n << ')' << endl;
     55    testerror=true;
     56  }
     57
     58  error=this_norm(V.transpose()*V)-n;
     59  if (error>MAXTOL) {
     60    cerr << "test_svd: FAILED, algorithm " << algo << " V orthogonality error ("
     61         << error << ") > tolerance (" << MAXTOL << ')' << endl;
     62    testerror=true;
     63  }
     64
     65  error=this_norm(U.transpose()*U)-n;
     66  if (error>MAXTOL) {
     67    cerr << "test_svd: FAILED, algorithm " << algo << " U orthogonality error ("
     68         << error << ") > tolerance (" << MAXTOL << ')' << endl;
     69    testerror=true;
     70  }
     71  return testerror;
     72}
     73
     74
     75
     76int main(const int argc,const char* argv[])
     77{
     78  bool testfail=false;
     79
     80  // The GSL Jacobi, Golub-Reinsch, and modified Golub-Reinsch
     81  // implementations supports rows>=columns matrix dimensions only
     82  testfail|=test(12,12,SVD::GolubReinsch);
     83  testfail|=test(12,4,SVD::GolubReinsch);
     84  testfail|=test(12,12,SVD::ModifiedGolubReinsch);
     85  testfail|=test(12,4,SVD::ModifiedGolubReinsch);
     86  testfail|=test(12,12,SVD::Jacobi);
     87  testfail|=test(12,4,SVD::Jacobi);
     88
     89  return (testfail ? -1 : 0);
     90}
Note: See TracChangeset for help on using the changeset viewer.