source: trunk/src/SVD.cc @ 42

Last change on this file since 42 was 42, checked in by Jari Häkkinen, 19 years ago

Made a major revision of matrix and vector classes. Everything compiles
but the binaries have not been tested.

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