Changeset 11 for trunk/src/SVD.cc


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

Hmm...

File:
1 edited

Legend:

Unmodified
Added
Removed
  • 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;
Note: See TracChangeset for help on using the changeset viewer.