source: trunk/src/SVD.cc @ 11

Last change on this file since 11 was 11, checked in by daniel, 19 years ago

Hmm...

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 6.8 KB
Line 
1// $Id: SVD.cc 11 2003-06-03 16:21:31Z daniel $
2
3#include "SVD.h"
4using namespace thep_cpp_tools;
5
6//Constructor that can only be used for test-method
7SVD::SVD() : loaded_( false ), solver_( false )
8{}
9
10//Constructor initializing the svd object with matrix A to be decomposed
11SVD::SVD( const thep_gsl_api::matrix& Ain ) : A_( Ain ) 
12{
13  //Get dimensions
14  ////////////////
15  size_t n = Ain.cols();
16
17  //Instantiate matrix with correct dimenstions
18  //and assign all elements the value zero
19  /////////////////////////////////////////////
20  V_ = thep_gsl_api::matrix( n, n, true );
21 
22  //Assign vector s_ with correct dimenstion
23  /////////////////////////////////////////////////
24  s_ = thep_gsl_api::vector( n );
25 
26  loaded_ = true;
27  solver_ = false;
28}
29
30
31//Constructor initializing the svd object for Solver
32//Solver requires an extra-vector b:
33//Ax = b
34//The x vector will be reached using the get_x() function.
35SVD::SVD( const thep_gsl_api::matrix& Ain, const thep_gsl_api::vector& b ) : A_( Ain ), b_( b ) 
36{
37  //Get dimensions
38  ////////////////
39  size_t n = Ain.cols();
40 
41  //Instantiate matrix with correct dimenstions
42  //and assign all elements the value zero
43  /////////////////////////////////////////////
44    V_ = thep_gsl_api::matrix( n, n, true );
45
46
47  //Assign vector s_ with correct dimenstion
48  /////////////////////////////////////////////////
49  s_ = thep_gsl_api::vector( n );
50 
51 
52  //An extra vector holding the answer to equation Ax = b is required
53  ///////////////////////////////////////////////////////////////////
54  x_ = thep_gsl_api::vector( n );
55
56
57  //Also set solver to true
58  /////////////////////////
59  solver_ = true;
60  loaded_ = true;
61}
62
63
64//Destructor
65SVD::~SVD()
66{
67}
68
69
70//Method to process calculation
71bool SVD::process( SVDtypes stype )
72{
73  if( !loaded_ ) 
74    { 
75     std::cerr << "Use proper constructor, matrix A not loaded!\n"; 
76     return false; 
77    }
78
79  switch( stype )
80    {
81
82    case Unmodified:
83      if( !process_default() ) exit( -1 );
84      break;
85
86    case Modified:
87      if( !process_modified() ) exit( -1 );
88      break;
89
90    case Jacobi:
91      if( !process_jacobi() ) exit( -1 );
92      break;
93
94    case Solver:
95      if( !process_solver() ) exit( -1 );
96      break;
97
98    default:
99      std::cerr << "Programmer try to use method not implemented!\n";
100      return false;
101      break;
102
103    }
104  return true;
105}
106
107
108//Default-method for SVD
109bool SVD::process_default()
110{
111  //Perform SVD using GSL library
112  //An additional workvector is
113  //required here
114  ///////////////////////////////
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() 
118          ) != 0 ) 
119    {
120     std::cerr << "failed to run gsl_linalg_SV_decomp!\n";
121     return false; 
122    }
123
124  return true;
125}
126
127
128//Modified SVD method
129//This method is faster when M>>N
130bool SVD::process_modified()
131{
132  //Perform SVD using GSL library
133  //An additional workmatrix X is
134  //required here
135  ///////////////////////////////
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() 
143        ) != 0 ) 
144    {
145     std::cerr << "failed to run gsl_linalg_SV_decomp_mod!\n";
146     return false; 
147    }
148
149  return true;
150}
151
152
153//Default-method for SVD
154//Computes singular values to higer
155//accuracy
156bool SVD::process_jacobi()
157{
158  //Perform SVD using GSL library
159  ///////////////////////////////
160  if( gsl_linalg_SV_decomp_jacobi( A_.get_gsl_matrix(), V_.get_gsl_matrix(), 
161           s_.get_gsl_vector() 
162           ) != 0 ) 
163    {
164     std::cerr << "failed to run gsl_linalg_SV_decomp_jacobi!\n";
165     return false; 
166    }
167
168  return true;
169}
170
171
172//Solver solves the system Ax=b
173//using SVD (process_default).
174bool SVD::process_solver()
175{
176  //Perform SVD using GSL library
177  ///////////////////////////////
178  this->process_default();
179 
180  //Solve ...
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() 
184         ) != 0 ) 
185    {
186     std::cerr << "failed to run gsl_linalg_SV_solve!\n";
187     return false; 
188    }
189 
190 
191  return true;
192}
193
194
195
196//This method will create a diagonal matrix of vector s_
197thep_gsl_api::matrix SVD::get_s() const 
198{
199  size_t n = A_.rows();
200 
201  thep_gsl_api::matrix S( n, n, true );
202
203  for( size_t i = 0; i < n; ++i )
204    {
205     S.set( i, i, s_.get( i ) );
206    }
207
208  return S;
209}
210
211//This function will run the SVD function and print
212//matrices in matlab format so that you can verify
213//the answer.
214bool SVD::test()
215{
216  using namespace std;
217
218
219  //Initialize GNU random stuff
220  ///////////////////////////////////
221    srand( time( 0 ) );
222
223   
224  //Initialise random test-matrix
225  ///////////////////////////////
226    const size_t DIM = 3;
227    thep_gsl_api::matrix A( DIM, DIM );
228    for( size_t i = 0; i < DIM; ++i )
229      {
230       for( size_t j = 0; j < DIM; ++j )
231   {
232    A.set( i, j, 10 * (double) rand() / RAND_MAX ); //Use random numbers between 0 and 10
233   }
234      }
235
236
237    //Set data
238    //////////
239    A_ = A;
240    size_t n = A_.rows();
241    size_t m = A_.cols();
242    n = m = DIM;
243    loaded_ = true;
244
245
246    //Set work matrices
247    ///////////////////
248    V_ = thep_gsl_api::matrix( n, n, true );
249
250    //Assign vector s_ with correct dimensions
251    /////////////////////////////////////////////////
252    s_ = thep_gsl_api::vector( n );
253
254    double error = 1;
255    thep_gsl_api::matrix E, S;
256
257
258    //Print matrix to screen
259    ////////////////////////
260    cout << "A = " << endl;
261    cout << A << endl;
262
263    cout.setf( ios::scientific );
264
265
266    //Run test 1: Unmodified method
267    ///////////////////////////////
268    this->process();
269    S = this->get_s();
270
271    E = A - A_*S*V_.transpose();
272   
273
274    error = sqrt( E.norm( 2 ) );
275    cout << "error(Unmodified method) = " << error << endl;
276    cout << "to compare with " << MAXTOL << endl;
277    if( error < MAXTOL ) { cout << "test ok!!!\n"; }
278    else { cout << "test failed!\n"; return false; }
279
280 
281    //Run test 2: Modified method
282    ///////////////////////////////
283    A_ = A;
284    this->process( Modified );
285    S = this->get_s();
286    E = A - A_*S*V_.transpose();
287    error = sqrt( E.norm( 2 ) );
288    cout << "error(Modified method) = " << error << endl;
289    cout << "to compare with " << MAXTOL << endl;
290    if( error < MAXTOL ) { cout << "test ok!!!\n"; }
291    else { cout << "test failed!\n"; return false; }
292
293
294    //Run test 2: Modified method
295    ///////////////////////////////
296    A_ = A;
297    this->process( Jacobi );
298    S = this->get_s();
299    E = A - A_*S*V_.transpose();
300    error = sqrt( E.norm( 2 ) );
301    cout << "error(Jacobi method) = " << error << endl;
302    cout << "to compare with " << MAXTOL << endl;
303    if( error < MAXTOL ) { cout << "test ok!!!\n"; }
304    else { cout << "test failed!\n"; return false; }
305   
306    return true; //Test with good outcome!
307}
308
309
Note: See TracBrowser for help on using the repository browser.