Sample Code: sample69.cpp Source File

00001 
00002 #include <iostream>
00003 
00004 #include <boost/timer.hpp>
00005 #include <boost/numeric/ublas/matrix.hpp>
00006 #include <boost/numeric/ublas/banded.hpp>
00007 #include <boost/numeric/ublas/matrix_sparse.hpp>
00008 #include <boost/numeric/ublas/triangular.hpp>
00009 #include <boost/numeric/ublas/operation.hpp>
00010 #include <boost/numeric/ublas/io.hpp>
00011 #include <boost/numeric/ublas/lu.hpp>
00012 
00013 using std::size_t;
00014 using std::cout;
00015 using std::endl;
00016 
00017 using namespace boost::numeric::ublas;
00018 
00019 template <class Matrix>
00020 void spy(std::ostream& os, const Matrix& A)
00021 {
00022   const double eps = 1e-10;
00023   typename Matrix::size_type i,j,m,n;
00024   m = A.size1();
00025   n = A.size2();
00026   for (i=0; i<m; ++i) {
00027         for (j=0; j<n; ++j) {
00028           if (A(i,j) < -eps) {
00029                 os << "-";
00030           } else if (A(i,j) > eps) {
00031                 os << "+";
00032           } else /* if A(i,j) == 0 */ {
00033                 os << " ";
00034           }
00035         }
00036         os << "\n";
00037   }
00038   os << "\n";
00039 }
00040 
00041 // #define DENSE
00042 // #define BANDED
00043 // #define SPARSE
00044 // #define COMPRESSED
00045 
00046 #ifndef NDEBUG
00047 const size_t default_size = 100;
00048 #else
00049 const size_t default_size = 4000;
00050 #endif
00051 const size_t sub1 = 1;
00052 const size_t sub2 = 10;
00053 const size_t sup1 = 1;
00054 const size_t sup2 = 10;
00055 
00056 void test_column_major(size_t size) {
00057   typedef permutation_matrix<std::size_t> MY_PERMA;
00058 #ifdef DENSE
00059   typedef matrix<double, column_major> MY_STIMA;
00060 #endif
00061 #ifdef BANDED
00062   typedef banded_matrix<double, column_major> MY_STIMA;
00063 #endif
00064 #ifdef SPARSE
00065   typedef sparse_matrix<double, column_major> MY_STIMA;
00066 #endif
00067 #ifdef COMPRESSED
00068   typedef compressed_matrix<double, column_major> MY_STIMA;
00069 #endif
00070   typedef vector<double> MY_VEC;
00071 
00072 #ifdef DENSE
00073   MY_STIMA M(size,size);
00074 #elif BANDED
00075   MY_STIMA M(size,size,sub2, sup2);
00076 #else
00077   MY_STIMA M(size,size,(sup2+1+sub2)*size);
00078 #endif
00079   MY_VEC   x(size);
00080   MY_VEC   y(size);
00081 
00082   boost::timer t;
00083 
00084   t.restart();
00085   for (size_t i=sup2; i<size; ++i)  M(i-sup2,i     ) =-1.0;
00086   for (size_t i=sup1; i<size; ++i)  M(i-sup1,i     ) =-1.0;
00087   for (size_t i=0;    i<size; ++i)  M(i     ,i     ) = 4.0;
00088   for (size_t i=sub1; i<size; ++i)  M(i     ,i-sub1) =-1.0;
00089   for (size_t i=sub2; i<size; ++i)  M(i     ,i-sub2) =-1.0;
00090 #if defined (DENSE) || defined (BANDED)
00091   cout << "assembly " << t.elapsed() << " " << endl;
00092 #else
00093   cout << "assembly " << t.elapsed() << " " << M.non_zeros() << endl;
00094 #endif
00095 
00096   //cout << "input M = " << M << "\n";
00097   if (size < 72) {
00098         // cout << "matrix spy of input matrix\n";
00099         t.restart();
00100         spy(cout, M);
00101 #if defined (DENSE) || defined (BANDED)
00102         cout << "spy input " << t.elapsed() << endl;
00103 #else
00104         cout << "spy input " << t.elapsed() << " " << M.non_zeros() << endl;
00105 #endif
00106   }
00107 
00108   for (size_t i=0; i<x.size(); ++i)  x(i) = 1.0;
00109 
00110   t.restart();
00111   axpy_prod(M, x, y, true);
00112 #if defined (DENSE) || defined (BANDED)
00113   cout << "axpy_prod " << t.elapsed() << endl;
00114 #else
00115   cout << "axpy_prod " << t.elapsed() << " " << M.non_zeros() << endl;
00116 #endif
00117 
00118   t.restart();
00119   MY_PERMA P (size);
00120   size_t rc = lu_factorize(M, P);
00121   assert(!rc);
00122 #if defined (DENSE) || defined (BANDED)
00123   cout << "lu_factorize " << t.elapsed() << endl;
00124 #else
00125   cout << "lu_factorize " << t.elapsed() << " " << M.non_zeros() << endl;
00126 #endif
00127 
00128   if (size < 72) {
00129         // cout << "matrix spy after LU decomposition\n";
00130         t.restart();
00131         spy(cout, M);
00132 #if defined (DENSE) || defined (BANDED)
00133         cout << "spy LU " << t.elapsed() << endl;
00134 #else
00135         cout << "spy LU " << t.elapsed() << " " << M.non_zeros() << endl;
00136 #endif
00137   }
00138   // cout << "output M = " << M << endl;
00139 
00140   // cout << "start solver\n";
00141 
00142   // why can't I use inplace_solve(M, y, unit_low_tag()) ?
00143   t.restart();
00144   lu_substitute(M, P, y);
00145 #if defined (DENSE) || defined (BANDED)
00146   cout << "lu_substitute " << t.elapsed() << endl;
00147 #else
00148   cout << "lu_substitute " << t.elapsed() << " " << M.non_zeros() << endl;
00149 #endif
00150 
00151   t.restart();
00152 #if defined (DENSE) || defined (BANDED)
00153   cout << "norm_2 " << t.elapsed() << endl;
00154 #else
00155   cout << "norm_2 " << t.elapsed() << " " << M.non_zeros() << endl;
00156 #endif
00157   cout << "L2 error : " << norm_2(x-y) << "\n";
00158 }
00159 
00160 void test_row_major(size_t size) {
00161   typedef permutation_matrix<std::size_t> MY_PERMA;
00162 #ifdef DENSE
00163   typedef matrix<double, row_major> MY_STIMA;
00164 #endif
00165 #ifdef BANDED
00166   typedef banded_matrix<double, row_major> MY_STIMA;
00167 #endif
00168 #ifdef SPARSE
00169   typedef sparse_matrix<double, row_major> MY_STIMA;
00170 #endif
00171 #ifdef COMPRESSED
00172   typedef compressed_matrix<double, row_major> MY_STIMA;
00173 #endif
00174   typedef vector<double> MY_VEC;
00175 
00176 #ifdef DENSE
00177   MY_STIMA M(size,size);
00178 #elif BANDED
00179   MY_STIMA M(size,size,sub2, sup2);
00180 #else
00181   MY_STIMA M(size,size,(sup2+1+sub2)*size);
00182 #endif
00183   MY_VEC   x(size);
00184   MY_VEC   y(size);
00185 
00186   boost::timer t;
00187 
00188   t.restart();
00189 // #define ORIGINAL
00190 #ifdef ORIGINAL
00191   for (size_t i=sup2; i<size; ++i)  M(i-sup2,i     ) =-1.0;
00192   for (size_t i=sup1; i<size; ++i)  M(i-sup1,i     ) =-1.0;
00193   for (size_t i=0;    i<size; ++i)  M(i     ,i     ) = 4.0;
00194   for (size_t i=sub1; i<size; ++i)  M(i     ,i-sub1) =-1.0;
00195   for (size_t i=sub2; i<size; ++i)  M(i     ,i-sub2) =-1.0;
00196 #else
00197   for (size_t i=sup2; i<size; ++i)  M(i     ,i-sup2) =-1.0;
00198   for (size_t i=sup1; i<size; ++i)  M(i     ,i-sup1) =-1.0;
00199   for (size_t i=0;    i<size; ++i)  M(i     ,i) = 4.0;
00200   for (size_t i=sub1; i<size; ++i)  M(i-sub1,i) =-1.0;
00201   for (size_t i=sub2; i<size; ++i)  M(i-sub2,i) =-1.0;
00202 #endif
00203 #if defined (DENSE) || defined (BANDED)
00204   cout << "assembly " << t.elapsed() << endl;
00205 #else
00206   cout << "assembly " << t.elapsed() << " " << M.non_zeros() << endl;
00207 #endif
00208 
00209   //cout << "input M = " << M << "\n";
00210   if (size < 72) {
00211         // cout << "matrix spy of input matrix\n";
00212         t.restart();
00213         spy(cout, M);
00214 #if defined (DENSE) || defined (BANDED)
00215         cout << "spy input " << t.elapsed() << endl;
00216 #else
00217         cout << "spy input " << t.elapsed() << " " << M.non_zeros() << endl;
00218 #endif
00219   }
00220 
00221   for (size_t i=0; i<x.size(); ++i)  x(i) = 1.0;
00222 
00223   t.restart();
00224 #ifdef ORIGINAL
00225   axpy_prod(M, x, y, true);
00226 #else
00227   axpy_prod(x, M, y, true);
00228 #endif
00229 #if defined (DENSE) || defined (BANDED)
00230   cout << "axpy_prod " << t.elapsed() << endl;
00231 #else
00232   cout << "axpy_prod " << t.elapsed() << " " << M.non_zeros() << endl;
00233 #endif
00234 
00235   t.restart();
00236   MY_PERMA P (size);
00237   size_t rc = lu_factorize(M, P);
00238   assert(!rc);
00239 #if defined (DENSE) || defined (BANDED)
00240   cout << "lu_factorize " << t.elapsed() << endl;
00241 #else
00242   cout << "lu_factorize " << t.elapsed() << " " << M.non_zeros() << endl;
00243 #endif
00244 
00245   if (size < 72) {
00246         // cout << "matrix spy after LU decomposition\n";
00247         t.restart();
00248         spy(cout, M);
00249 #if defined (DENSE) || defined (BANDED)
00250         cout << "spy LU " << t.elapsed() << endl;
00251 #else
00252         cout << "spy LU " << t.elapsed() << " " << M.non_zeros() << endl;
00253 #endif
00254   }
00255   // cout << "output M = " << M << endl;
00256 
00257   // cout << "start solver\n";
00258 
00259   // why can't I use inplace_solve(M, y, unit_low_tag()) ?
00260   t.restart();
00261 #ifdef ORIGINAL
00262   lu_substitute(M, P, y);
00263 #else
00264   lu_substitute(y, M, P);
00265 #endif
00266 #if defined (DENSE) || defined (BANDED)
00267   cout << "lu_substitute " << t.elapsed() << endl;
00268 #else
00269   cout << "lu_substitute " << t.elapsed() << " " << M.non_zeros() << endl;
00270 #endif
00271 
00272   t.restart();
00273 #if defined (DENSE) || defined (BANDED)
00274   cout << "norm_2 " << t.elapsed() << endl;
00275 #else
00276   cout << "norm_2 " << t.elapsed() << " " << M.non_zeros() << endl;
00277 #endif
00278   cout << "L2 error : " << norm_2(x-y) << "\n";
00279 }
00280 
00281 void test_column_major_axpy(size_t size) {
00282   typedef permutation_matrix<std::size_t> MY_PERMA;
00283 #ifdef DENSE
00284   typedef matrix<double, column_major> MY_STIMA;
00285 #endif
00286 #ifdef BANDED
00287   typedef banded_matrix<double, column_major> MY_STIMA;
00288 #endif
00289 #ifdef SPARSE
00290   typedef sparse_matrix<double, column_major> MY_STIMA;
00291 #endif
00292 #ifdef COMPRESSED
00293   typedef compressed_matrix<double, column_major> MY_STIMA;
00294 #endif
00295   typedef vector<double> MY_VEC;
00296 
00297 #ifdef DENSE
00298   MY_STIMA M(size,size);
00299 #elif BANDED
00300   MY_STIMA M(size,size,sub2, sup2);
00301 #else
00302   MY_STIMA M(size,size,(sup2+1+sub2)*size);
00303 #endif
00304   MY_VEC   x(size);
00305   MY_VEC   y(size);
00306 
00307   boost::timer t;
00308 
00309   t.restart();
00310   for (size_t i=sup2; i<size; ++i)  M(i-sup2,i     ) =-1.0;
00311   for (size_t i=sup1; i<size; ++i)  M(i-sup1,i     ) =-1.0;
00312   for (size_t i=0;    i<size; ++i)  M(i     ,i     ) = 4.0;
00313   for (size_t i=sub1; i<size; ++i)  M(i     ,i-sub1) =-1.0;
00314   for (size_t i=sub2; i<size; ++i)  M(i     ,i-sub2) =-1.0;
00315 #if defined (DENSE) || defined (BANDED)
00316   cout << "assembly " << t.elapsed() << " " << endl;
00317 #else
00318   cout << "assembly " << t.elapsed() << " " << M.non_zeros() << endl;
00319 #endif
00320 
00321   //cout << "input M = " << M << "\n";
00322   if (size < 72) {
00323         // cout << "matrix spy of input matrix\n";
00324         t.restart();
00325         spy(cout, M);
00326 #if defined (DENSE) || defined (BANDED)
00327         cout << "spy input " << t.elapsed() << endl;
00328 #else
00329         cout << "spy input " << t.elapsed() << " " << M.non_zeros() << endl;
00330 #endif
00331   }
00332 
00333   for (size_t i=0; i<x.size(); ++i)  x(i) = 1.0;
00334 
00335   t.restart();
00336   axpy_prod(M, x, y, true);
00337 #if defined (DENSE) || defined (BANDED)
00338   cout << "axpy_prod " << t.elapsed() << endl;
00339 #else
00340   cout << "axpy_prod " << t.elapsed() << " " << M.non_zeros() << endl;
00341 #endif
00342 
00343   t.restart();
00344   MY_PERMA P (size);
00345   size_t rc = axpy_lu_factorize(M, P);
00346   assert(!rc);
00347 #if defined (DENSE) || defined (BANDED)
00348   cout << "axpy_lu_factorize " << t.elapsed() << endl;
00349 #else
00350   cout << "axpy_lu_factorize " << t.elapsed() << " " << M.non_zeros() << endl;
00351 #endif
00352 
00353   if (size < 72) {
00354         // cout << "matrix spy after LU decomposition\n";
00355         t.restart();
00356         spy(cout, M);
00357 #if defined (DENSE) || defined (BANDED)
00358         cout << "spy LU " << t.elapsed() << endl;
00359 #else
00360         cout << "spy LU " << t.elapsed() << " " << M.non_zeros() << endl;
00361 #endif
00362   }
00363   // cout << "output M = " << M << endl;
00364 
00365   // cout << "start solver\n";
00366 
00367   // why can't I use inplace_solve(M, y, unit_low_tag()) ?
00368   t.restart();
00369   lu_substitute(M, P, y);
00370 #if defined (DENSE) || defined (BANDED)
00371   cout << "lu_substitute " << t.elapsed() << endl;
00372 #else
00373   cout << "lu_substitute " << t.elapsed() << " " << M.non_zeros() << endl;
00374 #endif
00375 
00376   t.restart();
00377 #if defined (DENSE) || defined (BANDED)
00378   cout << "norm_2 " << t.elapsed() << endl;
00379 #else
00380   cout << "norm_2 " << t.elapsed() << " " << M.non_zeros() << endl;
00381 #endif
00382   cout << "L2 error : " << norm_2(x-y) << "\n";
00383 }
00384 
00385 void test_row_major_axpy(size_t size) {
00386   typedef permutation_matrix<std::size_t> MY_PERMA;
00387 #ifdef DENSE
00388   typedef matrix<double, row_major> MY_STIMA;
00389 #endif
00390 #ifdef BANDED
00391   typedef banded_matrix<double, row_major> MY_STIMA;
00392 #endif
00393 #ifdef SPARSE
00394   typedef sparse_matrix<double, row_major> MY_STIMA;
00395 #endif
00396 #ifdef COMPRESSED
00397   typedef compressed_matrix<double, row_major> MY_STIMA;
00398 #endif
00399   typedef vector<double> MY_VEC;
00400 
00401 #ifdef DENSE
00402   MY_STIMA M(size,size);
00403 #elif BANDED
00404   MY_STIMA M(size,size,sub2, sup2);
00405 #else
00406   MY_STIMA M(size,size,(sup2+1+sub2)*size);
00407 #endif
00408   MY_VEC   x(size);
00409   MY_VEC   y(size);
00410 
00411   boost::timer t;
00412 
00413   t.restart();
00414 // #define ORIGINAL
00415 #ifdef ORIGINAL
00416   for (size_t i=sup2; i<size; ++i)  M(i-sup2,i     ) =-1.0;
00417   for (size_t i=sup1; i<size; ++i)  M(i-sup1,i     ) =-1.0;
00418   for (size_t i=0;    i<size; ++i)  M(i     ,i     ) = 4.0;
00419   for (size_t i=sub1; i<size; ++i)  M(i     ,i-sub1) =-1.0;
00420   for (size_t i=sub2; i<size; ++i)  M(i     ,i-sub2) =-1.0;
00421 #else
00422   for (size_t i=sup2; i<size; ++i)  M(i     ,i-sup2) =-1.0;
00423   for (size_t i=sup1; i<size; ++i)  M(i     ,i-sup1) =-1.0;
00424   for (size_t i=0;    i<size; ++i)  M(i     ,i) = 4.0;
00425   for (size_t i=sub1; i<size; ++i)  M(i-sub1,i) =-1.0;
00426   for (size_t i=sub2; i<size; ++i)  M(i-sub2,i) =-1.0;
00427 #endif
00428 #if defined (DENSE) || defined (BANDED)
00429   cout << "assembly " << t.elapsed() << endl;
00430 #else
00431   cout << "assembly " << t.elapsed() << " " << M.non_zeros() << endl;
00432 #endif
00433 
00434   //cout << "input M = " << M << "\n";
00435   if (size < 72) {
00436         // cout << "matrix spy of input matrix\n";
00437         t.restart();
00438         spy(cout, M);
00439 #if defined (DENSE) || defined (BANDED)
00440         cout << "spy input " << t.elapsed() << endl;
00441 #else
00442         cout << "spy input " << t.elapsed() << " " << M.non_zeros() << endl;
00443 #endif
00444   }
00445 
00446   for (size_t i=0; i<x.size(); ++i)  x(i) = 1.0;
00447 
00448   t.restart();
00449 #ifdef ORIGINAL
00450   axpy_prod(M, x, y, true);
00451 #else
00452   axpy_prod(x, M, y, true);
00453 #endif
00454 #if defined (DENSE) || defined (BANDED)
00455   cout << "axpy_prod " << t.elapsed() << endl;
00456 #else
00457   cout << "axpy_prod " << t.elapsed() << " " << M.non_zeros() << endl;
00458 #endif
00459 
00460   t.restart();
00461   MY_PERMA P (size);
00462   size_t rc = axpy_lu_factorize(M, P);
00463   assert(!rc);
00464 #if defined (DENSE) || defined (BANDED)
00465   cout << "axpy_lu_factorize " << t.elapsed() << endl;
00466 #else
00467   cout << "axpy_lu_factorize " << t.elapsed() << " " << M.non_zeros() << endl;
00468 #endif
00469 
00470   if (size < 72) {
00471         // cout << "matrix spy after LU decomposition\n";
00472         t.restart();
00473         spy(cout, M);
00474 #if defined (DENSE) || defined (BANDED)
00475         cout << "spy LU " << t.elapsed() << endl;
00476 #else
00477         cout << "spy LU " << t.elapsed() << " " << M.non_zeros() << endl;
00478 #endif
00479   }
00480   // cout << "output M = " << M << endl;
00481 
00482   // cout << "start solver\n";
00483 
00484   // why can't I use inplace_solve(M, y, unit_low_tag()) ?
00485   t.restart();
00486 #ifdef ORIGINAL
00487   lu_substitute(M, P, y);
00488 #else
00489   lu_substitute(y, M, P);
00490 #endif
00491 #if defined (DENSE) || defined (BANDED)
00492   cout << "lu_substitute " << t.elapsed() << endl;
00493 #else
00494   cout << "lu_substitute " << t.elapsed() << " " << M.non_zeros() << endl;
00495 #endif
00496 
00497   t.restart();
00498 #if defined (DENSE) || defined (BANDED)
00499   cout << "norm_2 " << t.elapsed() << endl;
00500 #else
00501   cout << "norm_2 " << t.elapsed() << " " << M.non_zeros() << endl;
00502 #endif
00503   cout << "L2 error : " << norm_2(x-y) << "\n";
00504 }
00505 
00506 int main(size_t argc, char *argv[])
00507 {
00508   size_t size = default_size;
00509   if (argc > 1)
00510       size = ::atoi (argv [1]);
00511   test_column_major(size);
00512   test_row_major(size);
00513   test_column_major_axpy(size);
00514   test_row_major_axpy(size);
00515   return EXIT_SUCCESS;
00516 };