source: trunk/src/SVD.cc @ 7

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

* empty log message *

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