Changeset 11


Ignore:
Timestamp:
Jun 3, 2003, 6:21:31 PM (20 years ago)
Author:
daniel
Message:

Hmm...

Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/Makefile.am

    r4 r11  
    55lib_LIBRARIES = libc++_tools.a
    66
    7 libc___tools_a_SOURCES = SVD.cc
     7libc___tools_a_SOURCES = matrix.cc vector.cc SVD.cc PCA.cc random_singleton.cc histogram.cc
    88
    99INCLUDES = -I/usr/local_bio/include -I/usr/local_bio

  • trunk/src/SVD.cc

    r7 r11  
    22
    33#include "SVD.h"
    4 #include <iostream>
    5 #include <cmath>   
    6 #include <cstdlib> //To include srand() and rand()
    7 #include <ctime>
    84using namespace thep_cpp_tools;
    95
    … …  
    139
    1410//Constructor initializing the svd object with matrix A to be decomposed
    15 SVD::SVD( const gsl::matrix& Ain ) : A_( Ain )
     11SVD::SVD( const thep_gsl_api::matrix& Ain ) : A_( Ain )
    1612{
    1713  //Get dimensions
    1814  ////////////////
    19   size_t n = Ain.get_rows();
     15  size_t n = Ain.cols();
    2016
    2117  //Instantiate matrix with correct dimenstions
    2218  //and assign all elements the value zero
    2319  /////////////////////////////////////////////
    24   V_.set_dimensions( n, n );
    25   V_.set_zero();
    26 
    27 
     20  V_ = thep_gsl_api::matrix( n, n, true );
     21 
    2822  //Assign vector s_ with correct dimenstion
    2923  /////////////////////////////////////////////////
    30   s_ = gsl::vector( n );
     24  s_ = thep_gsl_api::vector( n );
    3125 
    3226  loaded_ = true;
    … …  
    3933//Ax = b
    4034//The x vector will be reached using the get_x() function.
    41 SVD::SVD( const gsl::matrix& Ain, const gsl::vector& b ) : A_( Ain ), b_( b )
     35SVD::SVD( const thep_gsl_api::matrix& Ain, const thep_gsl_api::vector& b ) : A_( Ain ), b_( b )
    4236{
    4337  //Get dimensions
    4438  ////////////////
    45   size_t n = Ain.get_rows();
     39  size_t n = Ain.cols();
    4640 
    4741  //Instantiate matrix with correct dimenstions
    4842  //and assign all elements the value zero
    4943  /////////////////////////////////////////////
    50   V_.set_dimensions( n, n );
    51   V_.set_zero();
     44    V_ = thep_gsl_api::matrix( n, n, true );
    5245
    5346
    5447  //Assign vector s_ with correct dimenstion
    5548  /////////////////////////////////////////////////
    56   s_ = gsl::vector( n );
     49  s_ = thep_gsl_api::vector( n );
    5750 
    5851 
    5952  //An extra vector holding the answer to equation Ax = b is required
    6053  ///////////////////////////////////////////////////////////////////
    61   x_ = gsl::vector( n );
     54  x_ = thep_gsl_api::vector( n );
    6255
    6356
    … …  
    120113  //required here
    121114  ///////////////////////////////
    122   gsl::vector w( A_.get_rows() );
    123   if( gsl_linalg_SV_decomp( A_.gslobj(), V_.gslobj(),
    124           s_.gslobj(), w.gslobj()
     115  thep_gsl_api::vector w( A_.cols() );
     116  if( gsl_linalg_SV_decomp( A_.get_gsl_matrix(), V_.get_gsl_matrix(),
     117          s_.get_gsl_vector(), w.get_gsl_vector()
    125118          ) != 0 )
    126119    {
    … …  
    141134  //required here
    142135  ///////////////////////////////
    143   size_t n = A_.get_rows();
    144   gsl::vector w( n );
    145   gsl::matrix X;
    146 
    147   X.set_dimensions( n, n );
    148 
    149   if( gsl_linalg_SV_decomp_mod( A_.gslobj(), X.gslobj(),
    150         V_.gslobj(), s_.gslobj(),
    151         w.gslobj()
     136  size_t n = A_.rows();
     137  thep_gsl_api::vector w( n );
     138  thep_gsl_api::matrix X( n, n, true );
     139
     140  if( gsl_linalg_SV_decomp_mod( A_.get_gsl_matrix(), X.get_gsl_matrix(),
     141        V_.get_gsl_matrix(), s_.get_gsl_vector(),
     142        w.get_gsl_vector()
    152143        ) != 0 )
    153144    {
    … …  
    167158  //Perform SVD using GSL library
    168159  ///////////////////////////////
    169   if( gsl_linalg_SV_decomp_jacobi( A_.gslobj(), V_.gslobj(),
    170            s_.gslobj()
     160  if( gsl_linalg_SV_decomp_jacobi( A_.get_gsl_matrix(), V_.get_gsl_matrix(),
     161           s_.get_gsl_vector()
    171162           ) != 0 )
    172163    {
    … …  
    188179 
    189180  //Solve ...
    190   if( gsl_linalg_SV_solve( A_.gslobj(), V_.gslobj(),
    191          s_.gslobj(), b_.gslobj(),
    192          x_.gslobj()
     181  if( gsl_linalg_SV_solve( A_.get_gsl_matrix(), V_.get_gsl_matrix(),
     182         s_.get_gsl_vector(), b_.get_gsl_vector(),
     183         x_.get_gsl_vector()
    193184         ) != 0 )
    194185    {
    … …  
    204195
    205196//This method will create a diagonal matrix of vector s_
    206 gsl::matrix SVD::get_s() const
    207 {
    208   size_t n = A_.get_rows();
     197thep_gsl_api::matrix SVD::get_s() const
     198{
     199  size_t n = A_.rows();
    209200 
    210   gsl::matrix S;
    211   S.set_dimensions( n, n );
    212   S.set_zero();
     201  thep_gsl_api::matrix S( n, n, true );
    213202
    214203  for( size_t i = 0; i < n; ++i )
    215204    {
    216      S.set_element( i, i, s_ ( i ) );
     205     S.set( i, i, s_.get( i ) );
    217206    }
    218207
    … …  
    236225  ///////////////////////////////
    237226    const size_t DIM = 3;
    238     gsl::matrix A( DIM, DIM );
     227    thep_gsl_api::matrix A( DIM, DIM );
    239228    for( size_t i = 0; i < DIM; ++i )
    240229      {
    241230       for( size_t j = 0; j < DIM; ++j )
    242231   {
    243     A( i, j ) =  10 * (double) rand() / RAND_MAX; //Use random numbers between 0 and 10
     232    A.set( i, j, 10 * (double) rand() / RAND_MAX ); //Use random numbers between 0 and 10
    244233   }
    245234      }
    … …  
    249238    //////////
    250239    A_ = A;
    251     size_t n = A_.get_rows();
    252     size_t m = A_.get_cols();
     240    size_t n = A_.rows();
     241    size_t m = A_.cols();
    253242    n = m = DIM;
    254243    loaded_ = true;
    … …  
    257246    //Set work matrices
    258247    ///////////////////
    259     V_.set_dimensions( n, n );
    260     V_.set_zero();
    261 
     248    V_ = thep_gsl_api::matrix( n, n, true );
    262249
    263250    //Assign vector s_ with correct dimensions
    264251    /////////////////////////////////////////////////
    265     s_ = gsl::vector( n );
     252    s_ = thep_gsl_api::vector( n );
    266253
    267254    double error = 1;
    268     gsl::matrix E, S;
     255    thep_gsl_api::matrix E, S;
    269256
    270257
    … …  
    281268    this->process();
    282269    S = this->get_s();
     270
    283271    E = A - A_*S*V_.transpose();
     272   
     273
    284274    error = sqrt( E.norm( 2 ) );
    285275    cout << "error(Unmodified method) = " << error << endl;

  • trunk/src/SVD.h

    r7 r11  
    44#define _THEP_CPPTOOLS_SVD_
    55
    6 //GSLwrap includes
    7 //////////////////
    8 #include <gslwrap/matrix_double.h>
    9 #include <gslwrap/vector_double.h>
     6// C++ tools include
     7/////////////////////
     8#include "matrix.h"
     9#include "vector.h"
     10#include <gsl/gsl_vector.h>
     11
     12
     13// Standard includes
     14////////////////////
     15#include <iostream>
     16#include <cmath>   
     17#include <cstdlib> //To include srand() and rand()
     18#include <ctime>
    1019
    1120
    … …  
    2635  class SVD
    2736  {   
    28     static const double MAXTOL = 1E-10;
     37    static const double MAXTOL = 1E-6;
    2938
    3039    enum SVDtypes
    … …  
    3847  public:
    3948    /**
     49       The default constructor should only be used for testing!!!
     50    */
     51    SVD();
     52
     53   
     54    /**
    4055       Constructs an SVD object using the matrix A as only
    4156       input. The input matrix is copied for further use in the
    4257       object.
    4358    */
    44     SVD( const gsl::matrix& );
     59    SVD( const thep_gsl_api::matrix& );
    4560   
    4661    /**
    4762       Constructor initializing the SVD object for Solver. 
    48        Solver requires an extra-vector paramter \a b: $Ax = b$
    49        The $x$ vector will be reached using the get_x()
     63       Solver requires an extra-vector paramter \a b: \f$Ax = b\f$
     64       The \f$x\f$ vector will be reached using the get_x()
    5065       function.
    5166       @see get_x() function.
    5267    */
    5368
    54     SVD( const gsl::matrix&, const gsl::vector& b);
     69    SVD( const thep_gsl_api::matrix&, const thep_gsl_api::vector& b);
    5570    ~SVD();
    5671
    … …  
    6176       Modified: This method is faster when M>>N. (see GSL documentation
    6277       about gsl_linalg_SV_decomp_mod function)\n
    63        Jakoby: Computes singular values to higer accuracy than default method
     78       Jacobi: Computes singular values to higer accuracy than default method
    6479       (see GSL documentation about gsl_linalg_SV_decomp_jacobi.\n
    6580       Solver: Using Unmodified to perform SVD and then solves the system \f$Ax=b\f$
    … …  
    7186       executed.
    7287    */
    73     gsl::matrix get_u() const { return A_; }
     88    thep_gsl_api::matrix get_u() const { return A_; }
    7489
    7590    /**
    … …  
    7792       executed.
    7893    */
    79     gsl::matrix get_v() const { return V_; }
     94    thep_gsl_api::matrix get_v() const { return V_; }
    8095
    8196
    … …  
    87102    */
    88103
    89     gsl::vector get_x() const { return x_; }
     104    thep_gsl_api::vector get_x() const { return x_; }
    90105
    91106    /**
    … …  
    93108       executed.
    94109    */
    95     gsl::matrix get_s() const;
     110    thep_gsl_api::matrix get_s() const;
     111
    96112
    97113    /**
    98        This method will execute Default, Unmodified and Jakoby methods one at a time
     114       Function returning diagonal matrix S in vector form.
     115    */
     116    thep_gsl_api::vector get_s_vec() const { return s_; }
     117
     118
     119    /**
     120       This method will execute Default, Unmodified and Jacobi methods one at a time
    99121       and calculate the error, \f$ e = \left\Vert A - USV^{T} \right\Vert_{2} \f$.
    100        The method will return "true" if \f$e < 10^{-10}\f$ else false. A test program
     122       The method will return "true" if \f$e < 10^{-6}\f$ else false. A test program
    101123       is available that executes this method, c++_tools/test/svd_test. No test is
    102124       performed for solver but the aim is to add one in the future.
    … …  
    105127   
    106128  private:
    107     SVD();
    108129
    109     gsl::matrix A_, V_;   
    110     gsl::vector s_, b_, x_;
     130    thep_gsl_api::matrix A_, V_;   
     131    thep_gsl_api::vector s_, b_, x_;
    111132    bool loaded_, solver_;
    112133

  • trunk/test/Makefile.am

    r5 r11  
    33## $Id$
    44
    5 check_PROGRAMS = test_svd
     5check_PROGRAMS = test_rnd
     6#test_pca
     7#test_svd
    68
    7 test_svd_SOURCES = test_svd.cc
    8 test_svd_LDADD = -L@top_srcdir@/$(CPP_TOOLS_LIB_LOCATION) $(CPP_TOOLS_LIB) \
     9test_rnd_SOURCES = test_rnd.cc
     10test_rnd_LDADD = -L@top_srcdir@/$(CPP_TOOLS_LIB_LOCATION) $(CPP_TOOLS_LIB) \
    911    -L/usr/local_bio/lib -L/usr/local_bio/lib/Linux_P4SSE2 \
    1012    $(GSLWRAP_LIB) $(GSL_LIB) $(CBLAS_LIB) $(MATH_LIB)
    1113
     14#test_pca_SOURCES = test_pca.cc
     15#test_pca_LDADD = -L@top_srcdir@/$(CPP_TOOLS_LIB_LOCATION) $(CPP_TOOLS_LIB) \
     16#   -L/usr/local_bio/lib -L/usr/local_bio/lib/Linux_P4SSE2 \
     17#   $(GSLWRAP_LIB) $(GSL_LIB) $(CBLAS_LIB) $(MATH_LIB)
     18
     19#test_svd_SOURCES = test_svd.cc
     20#test_svd_LDADD = -L@top_srcdir@/$(CPP_TOOLS_LIB_LOCATION) $(CPP_TOOLS_LIB) \
     21#   -L/usr/local_bio/lib -L/usr/local_bio/lib/Linux_P4SSE2 \
     22#   $(GSLWRAP_LIB) $(GSL_LIB) $(CBLAS_LIB) $(MATH_LIB)
     23
    1224INCLUDES = -I@top_srcdir@/$(CPP_TOOLS_HEADER_LOCATION) \
    1325     -I/usr/local_bio/include
Note: See TracChangeset for help on using the changeset viewer.