Changeset 11 for trunk/src/SVD.cc
 Timestamp:
 Jun 3, 2003, 6:21:31 PM (19 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/src/SVD.cc
r7 r11 2 2 3 3 #include "SVD.h" 4 #include <iostream>5 #include <cmath>6 #include <cstdlib> //To include srand() and rand()7 #include <ctime>8 4 using namespace thep_cpp_tools; 9 5 … … 13 9 14 10 //Constructor initializing the svd object with matrix A to be decomposed 15 SVD::SVD( const gsl::matrix& Ain ) : A_( Ain )11 SVD::SVD( const thep_gsl_api::matrix& Ain ) : A_( Ain ) 16 12 { 17 13 //Get dimensions 18 14 //////////////// 19 size_t n = Ain. get_rows();15 size_t n = Ain.cols(); 20 16 21 17 //Instantiate matrix with correct dimenstions 22 18 //and assign all elements the value zero 23 19 ///////////////////////////////////////////// 24 V_.set_dimensions( n, n ); 25 V_.set_zero(); 26 27 20 V_ = thep_gsl_api::matrix( n, n, true ); 21 28 22 //Assign vector s_ with correct dimenstion 29 23 ///////////////////////////////////////////////// 30 s_ = gsl::vector( n );24 s_ = thep_gsl_api::vector( n ); 31 25 32 26 loaded_ = true; … … 39 33 //Ax = b 40 34 //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 )35 SVD::SVD( const thep_gsl_api::matrix& Ain, const thep_gsl_api::vector& b ) : A_( Ain ), b_( b ) 42 36 { 43 37 //Get dimensions 44 38 //////////////// 45 size_t n = Ain. get_rows();39 size_t n = Ain.cols(); 46 40 47 41 //Instantiate matrix with correct dimenstions 48 42 //and assign all elements the value zero 49 43 ///////////////////////////////////////////// 50 V_.set_dimensions( n, n ); 51 V_.set_zero(); 44 V_ = thep_gsl_api::matrix( n, n, true ); 52 45 53 46 54 47 //Assign vector s_ with correct dimenstion 55 48 ///////////////////////////////////////////////// 56 s_ = gsl::vector( n );49 s_ = thep_gsl_api::vector( n ); 57 50 58 51 59 52 //An extra vector holding the answer to equation Ax = b is required 60 53 /////////////////////////////////////////////////////////////////// 61 x_ = gsl::vector( n );54 x_ = thep_gsl_api::vector( n ); 62 55 63 56 … … 120 113 //required here 121 114 /////////////////////////////// 122 gsl::vector w( A_.get_rows() );123 if( gsl_linalg_SV_decomp( A_.g slobj(), V_.gslobj(),124 s_.g slobj(), 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() 125 118 ) != 0 ) 126 119 { … … 141 134 //required here 142 135 /////////////////////////////// 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() 152 143 ) != 0 ) 153 144 { … … 167 158 //Perform SVD using GSL library 168 159 /////////////////////////////// 169 if( gsl_linalg_SV_decomp_jacobi( A_.g slobj(), V_.gslobj(),170 s_.g slobj()160 if( gsl_linalg_SV_decomp_jacobi( A_.get_gsl_matrix(), V_.get_gsl_matrix(), 161 s_.get_gsl_vector() 171 162 ) != 0 ) 172 163 { … … 188 179 189 180 //Solve ... 190 if( gsl_linalg_SV_solve( A_.g slobj(), V_.gslobj(),191 s_.g slobj(), b_.gslobj(),192 x_.g slobj()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() 193 184 ) != 0 ) 194 185 { … … 204 195 205 196 //This method will create a diagonal matrix of vector s_ 206 gsl::matrix SVD::get_s() const207 { 208 size_t n = A_. get_rows();197 thep_gsl_api::matrix SVD::get_s() const 198 { 199 size_t n = A_.rows(); 209 200 210 gsl::matrix S; 211 S.set_dimensions( n, n ); 212 S.set_zero(); 201 thep_gsl_api::matrix S( n, n, true ); 213 202 214 203 for( size_t i = 0; i < n; ++i ) 215 204 { 216 S.set _element( i, i, s_( i ) );205 S.set( i, i, s_.get( i ) ); 217 206 } 218 207 … … 236 225 /////////////////////////////// 237 226 const size_t DIM = 3; 238 gsl::matrix A( DIM, DIM );227 thep_gsl_api::matrix A( DIM, DIM ); 239 228 for( size_t i = 0; i < DIM; ++i ) 240 229 { 241 230 for( size_t j = 0; j < DIM; ++j ) 242 231 { 243 A ( i, j ) = 10 * (double) rand() / RAND_MAX; //Use random numbers between 0 and 10232 A.set( i, j, 10 * (double) rand() / RAND_MAX ); //Use random numbers between 0 and 10 244 233 } 245 234 } … … 249 238 ////////// 250 239 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(); 253 242 n = m = DIM; 254 243 loaded_ = true; … … 257 246 //Set work matrices 258 247 /////////////////// 259 V_.set_dimensions( n, n ); 260 V_.set_zero(); 261 248 V_ = thep_gsl_api::matrix( n, n, true ); 262 249 263 250 //Assign vector s_ with correct dimensions 264 251 ///////////////////////////////////////////////// 265 s_ = gsl::vector( n );252 s_ = thep_gsl_api::vector( n ); 266 253 267 254 double error = 1; 268 gsl::matrix E, S;255 thep_gsl_api::matrix E, S; 269 256 270 257 … … 281 268 this>process(); 282 269 S = this>get_s(); 270 283 271 E = A  A_*S*V_.transpose(); 272 273 284 274 error = sqrt( E.norm( 2 ) ); 285 275 cout << "error(Unmodified method) = " << error << endl;
Note: See TracChangeset
for help on using the changeset viewer.