NumLinAlg.cxx
Go to the documentation of this file.
1 
14 #include "NumLinAlg.h"
15 
16 #include <fstream>
17 #include <iostream>
18 #include <iomanip>
19 #include <cmath>
20 #include <cfloat>
21 #include <cassert>
22 #include <cstdlib>
23 
24 using std::ofstream;
25 using std::ifstream;
26 using std::setw;
27 using std::setprecision;
28 using std::cout;
29 using std::endl;
30 using std::vector;
31 using std::abs;
32 
33 namespace hippodraw {
34  namespace Numeric {
35 
36 std::vector< double >
37 operator + ( const std::vector< double >& x,
38  const std::vector< double >& y )
39 {
40  int nrows = x.size();
41  vector< double > z;
42 
43  z.resize( nrows, 0.0 );
44 
45  for ( int i = 0; i < nrows; i++)
46  z[i] = x[i] + y[i];
47 
48  return z;
49 }
50 
51 std::vector< double >
52 operator - ( const std::vector< double >& x,
53  const std::vector< double >& y )
54 {
55  int nrows = x.size();
56  vector< double > z;
57 
58  z.resize( nrows, 0.0 );
59 
60  for ( int i = 0; i < nrows; i++)
61  z[i] = x[i] - y[i];
62 
63  return z;
64 }
65 
66 std::vector< double >
67 operator - ( const std::vector< double >& y )
68 {
69  int nrows = y.size();
70  vector< double > z;
71 
72  z.resize( nrows, 0.0 );
73 
74  for ( int i = 0; i < nrows; i++ )
75  z[i] = - y[i];
76 
77  return z;
78 }
79 
80 std::vector< vector< double > >
81 operator + ( const std::vector< std::vector< double > >&A,
82  const std::vector< std::vector< double > >&B )
83 {
84  int nrows = A.size();
85  int ncols = A[0].size();
86 
87  vector< vector< double > > C( nrows );
88  for( int i = 0; i < nrows; i++ )
89  C[i].resize( ncols, 0.0 );
90 
91  for (int i = 0; i < nrows; i++)
92  for (int j = 0; j < ncols; j++)
93  C[i][j] = A[i][j] + B[i][j];
94 
95  return C;
96 }
97 
98 std::vector< vector< double > >
99 operator - ( const std::vector< std::vector< double > >&A,
100  const std::vector< std::vector< double > >&B )
101 {
102  int nrows = A.size();
103  int ncols = A[0].size();
104 
105  vector< vector< double > > C( nrows );
106  for( int i = 0; i < nrows; i++ )
107  C[i].resize( ncols, 0.0 );
108 
109  for (int i = 0; i < nrows; i++)
110  for (int j = 0; j < ncols; j++)
111  C[i][j] = A[i][j] - B[i][j];
112 
113  return C;
114 }
115 
116 
117 std::vector< double >
118 operator * ( double a, const std::vector< double >& x )
119 {
120  int nrows = x.size();
121  vector< double > y;
122 
123  y.resize( nrows, 0.0 );
124 
125  for ( int i = 0; i < nrows; i++)
126  y[i] = a * x[i];
127 
128  return y;
129 }
130 
131 std::vector< double >
132 operator / ( const std::vector< double >& x, double a )
133 {
134  int nrows = x.size();
135  vector< double > y;
136 
137  assert( abs( a ) > DBL_EPSILON );
138 
139  y.resize( nrows, 0.0 );
140 
141  for ( int i = 0; i < nrows; i++)
142  y[i] = x[i] / a;
143 
144  return y;
145 }
146 
147 std::vector< std::vector< double > >
148 operator * ( double a, const std::vector< std::vector< double > >&A )
149 {
150  int nrows = A.size();
151  int ncols = A[0].size();
152 
153  vector< vector< double > > B( nrows );
154  for( int i = 0; i < nrows; i++ )
155  B[i].resize( ncols, 0.0 );
156 
157  for (int i = 0; i < nrows; i++)
158  for (int j = 0; j < ncols; j++)
159  B[i][j] = a * A[i][j];
160 
161  return B;
162 }
163 
164 std::vector< std::vector< double > >
165 operator / ( const std::vector< std::vector< double > >& A, double a )
166 {
167  int nrows = A.size();
168  int ncols = A[0].size();
169 
170  assert( abs( a ) > DBL_EPSILON );
171 
172  vector< vector< double > > B( nrows );
173  for( int i = 0; i < nrows; i++ )
174  B[i].resize( ncols, 0.0 );
175 
176  for (int i = 0; i < nrows; i++)
177  for (int j = 0; j < ncols; j++)
178  B[i][j] = A[i][j]/a;
179 
180  return B;
181 }
182 
183 std::vector< double >
184 operator * ( const std::vector< std::vector< double > >& A,
185  const std::vector< double >& x )
186 {
187  int nrows = A.size();
188  int ncols = A[0].size();
189  vector< double > y;
190 
191  y.resize( nrows, 0.0 );
192 
193  for (int i = 0; i < nrows; i++)
194  for (int j = 0; j < ncols; j++)
195  y[i] += A[i][j] * x[j];
196 
197  return y;
198 }
199 
200 std::vector< double >
201 operator * ( const std::vector< double >& x,
202  const std::vector< std::vector< double > >& A )
203 {
204  int nrows = A.size();
205  int ncols = A[0].size();
206  vector< double > y;
207 
208  y.resize( ncols, 0.0 );
209 
210  for (int j = 0; j < ncols; j++)
211  for (int i = 0; i < nrows; i++)
212  y[j] += A[i][j] * x[i];
213 
214  return y;
215 }
216 
217 std::vector< vector< double > >
218 operator * ( const std::vector< std::vector< double > >&A,
219  const std::vector< std::vector< double > >&B )
220 {
221  int m = A.size();
222  int r = A[0].size();
223  int n = B[0].size();
224 
225  vector< vector< double > > C( m );
226  for( int i = 0; i < m; i++ )
227  C[i].resize( n, 0.0 );
228 
229  for (int i = 0; i < m; i++)
230  for (int j = 0; j < n; j++)
231  for (int k = 0; k < r; k++)
232  C[i][j] += A[i][k] * B[k][j];
233 
234  return C;
235 }
236 
237 double
238 innerProduct( const std::vector< double > & a,
239  const std::vector< double > & b )
240 {
241  double prod = 0.0;
242 
243  for ( unsigned int i = 0; i < a.size(); i++ )
244  prod += a[i] * b[i];
245 
246  return prod;
247 }
248 
249 
250 vector< vector< double > >
251 outerProduct ( const std::vector< double >& a,
252  const std::vector< double >& b )
253 {
254  vector< vector< double> > C( a.size() );
255  for( unsigned int i = 0; i < a.size(); i++ )
256  C[i].resize( b.size() );
257 
258  for( unsigned int i = 0; i < a.size(); i++ )
259  for( unsigned int j = 0; j < b.size(); j++ )
260  C[i][j] = a[i] * b[j];
261 
262  return C;
263 }
264 
265 
266 double
267 quadraticProduct( const std::vector< std::vector< double > > & A,
268  const std::vector< double > x )
269 {
270  double prod = 0.0;
271  unsigned int n = A.size();
272 
273  assert ( x.size() == n );
274 
275  for ( unsigned int i = 0; i < n; i++ )
276  for ( unsigned int j = 0; j < n; j++ )
277  prod += x[i] * A[i][j] * x[j];
278 
279  return prod;
280 }
281 
282 
283 double
284 norm ( const std::vector< double > & a )
285 {
286  return sqrt( innerProduct( a, a ) );
287 }
288 
289 
290 int
291 cholFactor ( std::vector < std::vector< double > > & A )
292 {
293  unsigned int n = A.size();
294 
295  for ( unsigned int i = 0; i < n ; ++i )
296  for ( unsigned int j = 0; j <= i ; ++j)
297  {
298  double sum = A[i][j];
299 
300  A[j][i] = 0;
301 
302  for ( unsigned int k = 0; k < j; ++k )
303  sum -= A[i][k] * A[j][k];
304 
305  if (i == j)
306  {
307  if (sum < 0) return EXIT_FAILURE;
308 
309  sum = sqrt (sum);
310 
311  if ( fabs(sum) < DBL_EPSILON ) return EXIT_FAILURE;
312 
313  A[i][j] = sum;
314  }
315  else
316  A[i][j] = sum / A[j][j];
317  }
318 
319  return EXIT_SUCCESS;
320 }
321 
322 
323 int
324 cholBackSolve ( const std::vector < std::vector< double > > & L,
325  std::vector< double > & x,
326  const std::vector< double > & b )
327 {
328  unsigned int n = L.size();
329 
330  double sum;
331 
332  // Solving L y = b
333  for ( unsigned int i = 0; i < n ; i++)
334  {
335  sum = b [i];
336  for ( unsigned int j = 0; j < i ; ++j)
337  sum -= x[j] * L[i][j];
338  x[i] = sum / L[i][i];
339  }
340 
341  // Solving L' x = y
342  for ( int i = n - 1; i >= 0; i-- )
343  {
344  sum = x[i];
345  for ( unsigned int j = i + 1; j < n ; j++)
346  sum -= x[j] * L[j][i];
347  x[i] = sum / L[i][i];
348  }
349 
350  return EXIT_SUCCESS;
351 }
352 
353 
354 int
355 invertMatrix ( const std::vector < std::vector< double > > & A,
356  std::vector < std::vector < double > > & Ainv )
357 {
358  unsigned int n = A.size();
359  vector< double > x, b;
360  vector< vector< double > > L;
361 
362  // Set appropriate sizes for the vectors and matrices
363  x.resize( n, 0.0 );
364  b.resize( n, 0.0 );
365 
366  L.resize( n );
367  Ainv.resize( n );
368 
369  for ( unsigned int i = 0; i < n; i++ )
370  {
371  L[i].clear ();
372  L[i].resize ( n, 0.0 );
373 
374  Ainv[i].clear ();
375  Ainv[i].resize ( n, 0.0 );
376  }
377 
378  for ( unsigned int i = 0; i < n; i++ )
379  for ( unsigned int j = 0; j < n; j++ )
380  L[i][j] = A[i][j];
381 
382  // Do a cholesky factorization writing the lower triangular factor
383  // in the lower triangular part of L and setting the upper part as 0
384  cholFactor( L );
385 
386  for ( unsigned int j = 0; j < n; j++ )
387  {
388  // Set up right hand side i.e. set b = ei
389  for ( unsigned int i = 0; i < n; i++ )
390  b[i] = ( i == j ) ? 1 : 0;
391 
392  // LL'x= b is being solved here
393  cholBackSolve( L, x, b );
394 
395  // This x constitutes the j th coloumn of the inverse
396  for ( unsigned int i = 0; i < n; i++ )
397  Ainv[i][j] = x[i] ;
398  }
399 
400  return EXIT_SUCCESS;
401 }
402 
403 int
404 eye ( std::vector < std::vector < double > >& I, int n )
405 {
406  I.resize( n );
407  for( int i = 0; i < n; i++ )
408  {
409  I[i].clear();
410  I[i].resize( n, 0.0 );
411  }
412 
413  for( int i = 0; i < n; i++ )
414  I[i][i] = 1.0;
415 
416  return EXIT_SUCCESS;
417 }
418 
419 int
420 write ( const std::vector < double > & a )
421 {
422  unsigned int n = a.size();
423 
424  for ( unsigned int i = 0; i < n; ++i )
425  cout << setprecision( 15 ) << a[i] << endl;
426  cout << endl;
427 
428  return EXIT_SUCCESS;
429 }
430 
431 
432 int
433 write ( const std::vector < std::vector < double > > & A )
434 {
435  unsigned int n = A.size();
436 
437  for ( unsigned int i = 0; i < n; ++i )
438  {
439  for ( unsigned int j = 0; j < n; ++j )
440  cout << setw( 8 ) << setprecision( 4 ) << A[i][j];
441  cout << endl;
442  }
443 
444  cout << endl;
445 
446  return EXIT_SUCCESS;
447 }
448 
449 
450 int
451 allocateMatrix ( std::vector < std::vector < double > > & A,
452  int m, int n )
453 {
454  A.resize( m );
455  for( int i = 0; i < m; i++ )
456  {
457  A[i].clear();
458  A[i].resize( n, 0.0 );
459  }
460 
461  return EXIT_SUCCESS;
462 }
463 
464 
465 int
466 allocateVector ( std::vector < double > & x, int n )
467 {
468  x.clear();
469  x.resize( n, 0.0 );
470 
471  return EXIT_SUCCESS;
472 }
473 
474 
475 /* The driver main subroutine to check a few of above function.
476  Test with the following matlab code:
477  n = 4; L = tril( randn(n,n), -1 ) + eye(n,n); A = L' * L; B = randn(n , n); x = randn(n , 1); y = randn(n , 1); save test.txt -ascii A B x y;
478  Then perform the operations as stated in the following program.
479 */
480 /*int main()
481 {
482  int n = 4;
483  vector< double > x, y, z;
484  vector< vector < double > > A, B, Ainv;
485  ifstream fin( "test.txt" );
486 
487  if( !fin )
488  {
489  cout << " Could not open the file for reading " << endl;
490  return EXIT_SUCCESS;
491  }
492 
493  allocateMatrix( A, n, n );
494  allocateMatrix( B, n, n );
495  allocateVector( x, n );
496  allocateVector( y, n );
497 
498  for( int i = 0; i < n; i++ )
499  for( int j = 0; j < n; j++ )
500  fin >> A[i][j];
501 
502  for( int i = 0; i < n; i++ )
503  for( int j = 0; j < n; j++ )
504  fin >> B[i][j];
505 
506  for( int i = 0; i < n; i++ )
507  fin >> x[i];
508 
509  for( int i = 0; i < n; i++ )
510  fin >> y[i];
511 
512  fin.close();
513 
514  //cout << "x + 3.14 * x + y /1.4141 - 2.8 * A * x - y + (y' * A)'= " << endl;
515  //write( x + 3.141 * x + y / 1.4141 - 2.8 * A * x - y + y * A );
516 
517  //cout << "A + 3.14 * B + A /1.4141 - A * B + 65.0 * x * y' - B = " << endl;
518  //write( A + 3.14 * B + A /1.4141 - A * B + 65.0 * x * y ) - B );
519 
520  //cout << "inv(A) = " << endl;
521  //invertMatrix( A, Ainv );
522  //write( Ainv );
523 
524 
525  double temp = ( 1 + innerProduct( y, B * y ) / ys ) / ys;
526 
527  t1 = ( outerProduct(s, y) * B)/ys + ( m_M * outerProduct(y , s) ) / ys;
528  t2 = temp * outerProduct( s, s );
529  B = B - t1 + t2;
530 
531 
532  return EXIT_SUCCESS;
533 }
534 */
535 
536 } // namespace Numeric
537 } // namespace hippodraw

Generated for HippoDraw Class Library by doxygen