libs/numeric/ublas/benchmarks/bench3/bench32.cpp

   1 //
   2 //  Copyright (c) 2000-2002
   3 //  Joerg Walter, Mathias Koch
   4 //
   5 //  Distributed under the Boost Software License, Version 1.0. (See
   6 //  accompanying file LICENSE_1_0.txt or copy at
   7 //  http://www.boost.org/LICENSE_1_0.txt)
   8 //
   9 //  The authors gratefully acknowledge the support of
  10 //  GeNeSys mbH & Co. KG in producing this work.
  11 //
  12
  13 #include "bench3.hpp"
  14
  15 template<class T, int N>
  16 struct bench_c_outer_prod {
  17     typedef T value_type;
  18
  19     void operator () (int runs) const {
  20         try {
  21             static typename c_matrix_traits<T, N, N>::type m;
  22             static typename c_vector_traits<T, N>::type v1, v2;
  23             initialize_c_vector<T, N> () (v1);
  24             initialize_c_vector<T, N> () (v2);
  25             boost::timer t;
  26             for (int i = 0; i < runs; ++ i) {
  27                 for (int j = 0; j < N; ++ j) {
  28                     for (int k = 0; k < N; ++ k) {
  29                         m [j] [k] = - v1 [j] * v2 [k];
  30                     }
  31                 }
  32 //                sink_c_matrix<T, N, N> () (m);
  33             }
  34             footer<value_type> () (N * N, N * N, runs, t.elapsed ());
  35         }
  36         catch (std::exception &e) {
  37             std::cout << e.what () << std::endl;
  38         }
  39     }
  40 };
  41 template<class M, class V, int N>
  42 struct bench_my_outer_prod {
  43     typedef typename M::value_type value_type;
  44
  45     void operator () (int runs, safe_tag) const {
  46         try {
  47             static M m (N, N);
  48             ublas::matrix_range<M> mr (m, ublas::range (0, N), ublas::range (0, N));
  49             static V v1 (N), v2 (N);
  50             ublas::vector_range<V> vr1 (v1, ublas::range (0, N)),
  51                                    vr2 (v2, ublas::range (0, N));
  52             initialize_vector (vr1);
  53             initialize_vector (vr2);
  54             boost::timer t;
  55             for (int i = 0; i < runs; ++ i) {
  56                 mr = - ublas::outer_prod (vr1, vr2);
  57 //                sink_matrix (mr);
  58             }
  59             footer<value_type> () (N * N, N * N, runs, t.elapsed ());
  60         }
  61         catch (std::exception &e) {
  62             std::cout << e.what () << std::endl;
  63         }
  64     }
  65     void operator () (int runs, fast_tag) const {
  66         try {
  67             static M m (N, N);
  68             ublas::matrix_range<M> mr (m, ublas::range (0, N), ublas::range (0, N));
  69             static V v1 (N), v2 (N);
  70             ublas::vector_range<V> vr1 (v1, ublas::range (0, N)),
  71                                    vr2 (v2, ublas::range (0, N));
  72             initialize_vector (vr1);
  73             initialize_vector (vr2);
  74             boost::timer t;
  75             for (int i = 0; i < runs; ++ i) {
  76                 mr.assign (- ublas::outer_prod (vr1, vr2));
  77 //                sink_matrix (mr);
  78             }
  79             footer<value_type> () (N * N, N * N, runs, t.elapsed ());
  80         }
  81         catch (std::exception &e) {
  82             std::cout << e.what () << std::endl;
  83         }
  84     }
  85 };
  86 template<class M, class V, int N>
  87 struct bench_cpp_outer_prod {
  88     typedef typename M::value_type value_type;
  89
  90     void operator () (int runs) const {
  91         try {
  92             static M m (N * N);
  93             static V v1 (N), v2 (N);
  94             initialize_vector (v1);
  95             initialize_vector (v2);
  96             boost::timer t;
  97             for (int i = 0; i < runs; ++ i) {
  98                 for (int j = 0; j < N; ++ j) {
  99                     for (int k = 0; k < N; ++ k) {
 100                         m [N * j + k] = - v1 [j] * v2 [k];
 101                     }
 102                 }
 103 //                sink_vector (m);
 104             }
 105             footer<value_type> () (N * N, N * N, runs, t.elapsed ());
 106         }
 107         catch (std::exception &e) {
 108             std::cout << e.what () << std::endl;
 109         }
 110     }
 111 };
 112
 113 template<class T, int N>
 114 struct bench_c_matrix_vector_prod {
 115     typedef T value_type;
 116
 117     void operator () (int runs) const {
 118         try {
 119             static typename c_matrix_traits<T, N, N>::type m;
 120             static typename c_vector_traits<T, N>::type v1, v2;
 121             initialize_c_matrix<T, N, N> () (m);
 122             initialize_c_vector<T, N> () (v1);
 123             boost::timer t;
 124             for (int i = 0; i < runs; ++ i) {
 125                 for (int j = 0; j < N; ++ j) {
 126                     v2 [j] = 0;
 127                     for (int k = 0; k < N; ++ k) {
 128                         v2 [j] += m [j] [k] * v1 [k];
 129                     }
 130                 }
 131 //                sink_c_vector<T, N> () (v2);
 132             }
 133             footer<value_type> () (N * N, N * (N - 1), runs, t.elapsed ());
 134         }
 135         catch (std::exception &e) {
 136             std::cout << e.what () << std::endl;
 137         }
 138     }
 139 };
 140 template<class M, class V, int N>
 141 struct bench_my_matrix_vector_prod {
 142     typedef typename M::value_type value_type;
 143
 144     void operator () (int runs, safe_tag) const {
 145         try {
 146             static M m (N, N);
 147             ublas::matrix_range<M> mr (m, ublas::range (0, N), ublas::range (0, N));
 148             static V v1 (N), v2 (N);
 149             ublas::vector_range<V> vr1 (v1, ublas::range (0, N)),
 150                                    vr2 (v2, ublas::range (0, N));
 151             initialize_matrix (mr);
 152             initialize_vector (vr1);
 153             boost::timer t;
 154             for (int i = 0; i < runs; ++ i) {
 155                 vr2 = ublas::prod (mr, vr1);
 156 //                sink_vector (vr2);
 157             }
 158             footer<value_type> () (N * N, N * (N - 1), runs, t.elapsed ());
 159         }
 160         catch (std::exception &e) {
 161             std::cout << e.what () << std::endl;
 162         }
 163     }
 164     void operator () (int runs, fast_tag) const {
 165         try {
 166             static M m (N, N);
 167             ublas::matrix_range<M> mr (m, ublas::range (0, N), ublas::range (0, N));
 168             static V v1 (N), v2 (N);
 169             ublas::vector_range<V> vr1 (v1, ublas::range (0, N)),
 170                                    vr2 (v2, ublas::range (0, N));
 171             initialize_matrix (mr);
 172             initialize_vector (vr1);
 173             boost::timer t;
 174             for (int i = 0; i < runs; ++ i) {
 175                 vr2.assign (ublas::prod (mr, vr1));
 176 //                sink_vector (vr2);
 177             }
 178             footer<value_type> () (N * N, N * (N - 1), runs, t.elapsed ());
 179         }
 180         catch (std::exception &e) {
 181             std::cout << e.what () << std::endl;
 182         }
 183     }
 184 };
 185 template<class M, class V, int N>
 186 struct bench_cpp_matrix_vector_prod {
 187     typedef typename M::value_type value_type;
 188
 189     void operator () (int runs) const {
 190         try {
 191             static M m (N * N);
 192             static V v1 (N), v2 (N);
 193             initialize_vector (m);
 194             initialize_vector (v1);
 195             boost::timer t;
 196             for (int i = 0; i < runs; ++ i) {
 197                 for (int j = 0; j < N; ++ j) {
 198                     std::valarray<value_type> row (m [std::slice (N * j, N, 1)]);
 199                     v2 [j] = (row * v1).sum ();
 200                 }
 201 //                sink_vector (v2);
 202             }
 203             footer<value_type> () (N * N, N * (N - 1), runs, t.elapsed ());
 204         }
 205         catch (std::exception &e) {
 206             std::cout << e.what () << std::endl;
 207         }
 208     }
 209 };
 210
 211 template<class T, int N>
 212 struct bench_c_matrix_add {
 213     typedef T value_type;
 214
 215     void operator () (int runs) const {
 216         try {
 217             static typename c_matrix_traits<T, N, N>::type m1, m2, m3;
 218             initialize_c_matrix<T, N, N> () (m1);
 219             initialize_c_matrix<T, N, N> () (m2);
 220             boost::timer t;
 221             for (int i = 0; i < runs; ++ i) {
 222                 for (int j = 0; j < N; ++ j) {
 223                     for (int k = 0; k < N; ++ k) {
 224                         m3 [j] [k] = - (m1 [j] [k] + m2 [j] [k]);
 225                     }
 226                 }
 227 //                sink_c_matrix<T, N, N> () (m3);
 228             }
 229             footer<value_type> () (0, 2 * N * N, runs, t.elapsed ());
 230         }
 231         catch (std::exception &e) {
 232             std::cout << e.what () << std::endl;
 233         }
 234     }
 235 };
 236 template<class M, int N>
 237 struct bench_my_matrix_add {
 238     typedef typename M::value_type value_type;
 239
 240     void operator () (int runs, safe_tag) const {
 241         try {
 242             static M m1 (N, N), m2 (N, N), m3 (N, N);
 243             initialize_matrix (m1);
 244             initialize_matrix (m2);
 245             boost::timer t;
 246             for (int i = 0; i < runs; ++ i) {
 247                 m3 = - (m1 + m2);
 248 //                sink_matrix (m3);
 249             }
 250             footer<value_type> () (0, 2 * N * N, runs, t.elapsed ());
 251         }
 252         catch (std::exception &e) {
 253             std::cout << e.what () << std::endl;
 254         }
 255     }
 256     void operator () (int runs, fast_tag) const {
 257         try {
 258             static M m1 (N, N), m2 (N, N), m3 (N, N);
 259             initialize_matrix (m1);
 260             initialize_matrix (m2);
 261             boost::timer t;
 262             for (int i = 0; i < runs; ++ i) {
 263                 m3.assign (- (m1 + m2));
 264 //                sink_matrix (m3);
 265             }
 266             footer<value_type> () (0, 2 * N * N, runs, t.elapsed ());
 267         }
 268         catch (std::exception &e) {
 269             std::cout << e.what () << std::endl;
 270         }
 271     }
 272 };
 273 template<class M, int N>
 274 struct bench_cpp_matrix_add {
 275     typedef typename M::value_type value_type;
 276
 277     void operator () (int runs) const {
 278         try {
 279             static M m1 (N * N), m2 (N * N), m3 (N * N);
 280             initialize_vector (m1);
 281             initialize_vector (m2);
 282             boost::timer t;
 283             for (int i = 0; i < runs; ++ i) {
 284                 m3 = - (m1 + m2);
 285 //                sink_vector (m3);
 286             }
 287             footer<value_type> () (0, 2 * N * N, runs, t.elapsed ());
 288         }
 289         catch (std::exception &e) {
 290             std::cout << e.what () << std::endl;
 291         }
 292     }
 293 };
 294
 295 // Benchmark O (n ^ 2)
 296 template<class T, int N>
 297 void bench_2<T, N>::operator () (int runs) {
 298     header ("bench_2");
 299
 300     header ("outer_prod");
 301
 302     header ("C array");
 303     bench_c_outer_prod<T, N> () (runs);
 304
 305 #ifdef USE_C_ARRAY
 306     header ("c_matrix, c_vector safe");
 307     bench_my_outer_prod<ublas::c_matrix<T, N, N>,
 308                         ublas::c_vector<T, N>, N> () (runs, safe_tag ());
 309
 310     header ("c_matrix, c_vector fast");
 311     bench_my_outer_prod<ublas::c_matrix<T, N, N>,
 312                         ublas::c_vector<T, N>, N> () (runs, fast_tag ());
 313 #endif
 314
 315 #ifdef USE_BOUNDED_ARRAY
 316     header ("matrix<bounded_array>, vector<bounded_array> safe");
 317     bench_my_outer_prod<ublas::matrix<T, ublas::row_major, ublas::bounded_array<T, N * N> >,
 318                         ublas::vector<T, ublas::bounded_array<T, N> >, N> () (runs, safe_tag ());
 319
 320     header ("matrix<bounded_array>, vector<bounded_array> fast");
 321     bench_my_outer_prod<ublas::matrix<T, ublas::row_major, ublas::bounded_array<T, N * N> >,
 322                         ublas::vector<T, ublas::bounded_array<T, N> >, N> () (runs, fast_tag ());
 323 #endif
 324
 325 #ifdef USE_UNBOUNDED_ARRAY
 326     header ("matrix<unbounded_array>, vector<unbounded_array> safe");
 327     bench_my_outer_prod<ublas::matrix<T, ublas::row_major, ublas::unbounded_array<T> >,
 328                         ublas::vector<T, ublas::unbounded_array<T> >, N> () (runs, safe_tag ());
 329
 330     header ("matrix<unbounded_array>, vector<unbounded_array> fast");
 331     bench_my_outer_prod<ublas::matrix<T, ublas::row_major, ublas::unbounded_array<T> >,
 332                         ublas::vector<T, ublas::unbounded_array<T> >, N> () (runs, fast_tag ());
 333 #endif
 334
 335 #ifdef USE_STD_VALARRAY
 336     header ("matrix<std::valarray>, vector<std::valarray> safe");
 337     bench_my_outer_prod<ublas::matrix<T, ublas::row_major, std::valarray<T> >,
 338                         ublas::vector<T, std::valarray<T> >, N> () (runs, safe_tag ());
 339
 340     header ("matrix<std::valarray>, vector<std::valarray> fast");
 341     bench_my_outer_prod<ublas::matrix<T, ublas::row_major, std::valarray<T> >,
 342                         ublas::vector<T, std::valarray<T> >, N> () (runs, fast_tag ());
 343 #endif
 344
 345 #ifdef USE_STD_VECTOR
 346     header ("matrix<std::vector>, vector<std::vector> safe");
 347     bench_my_outer_prod<ublas::matrix<T, ublas::row_major, std::vector<T> >,
 348                         ublas::vector<T, std::vector<T> >, N> () (runs, safe_tag ());
 349
 350     header ("matrix<std::vector>, vector<std::vector> fast");
 351     bench_my_outer_prod<ublas::matrix<T, ublas::row_major, std::vector<T> >,
 352                         ublas::vector<T, std::vector<T> >, N> () (runs, fast_tag ());
 353 #endif
 354
 355 #ifdef USE_STD_VALARRAY
 356     header ("std::valarray");
 357     bench_cpp_outer_prod<std::valarray<T>, std::valarray<T>, N> () (runs);
 358 #endif
 359
 360     header ("prod (matrix, vector)");
 361
 362     header ("C array");
 363     bench_c_matrix_vector_prod<T, N> () (runs);
 364
 365 #ifdef USE_C_ARRAY
 366     header ("c_matrix, c_vector safe");
 367     bench_my_matrix_vector_prod<ublas::c_matrix<T, N, N>,
 368                                 ublas::c_vector<T, N>, N> () (runs, safe_tag ());
 369
 370     header ("c_matrix, c_vector fast");
 371     bench_my_matrix_vector_prod<ublas::c_matrix<T, N, N>,
 372                                 ublas::c_vector<T, N>, N> () (runs, fast_tag ());
 373 #endif
 374
 375 #ifdef USE_BOUNDED_ARRAY
 376     header ("matrix<bounded_array>, vector<bounded_array> safe");
 377     bench_my_matrix_vector_prod<ublas::matrix<T, ublas::row_major, ublas::bounded_array<T, N * N> >,
 378                                 ublas::vector<T, ublas::bounded_array<T, N> >, N> () (runs, safe_tag ());
 379
 380     header ("matrix<bounded_array>, vector<bounded_array> fast");
 381     bench_my_matrix_vector_prod<ublas::matrix<T, ublas::row_major, ublas::bounded_array<T, N * N> >,
 382                                 ublas::vector<T, ublas::bounded_array<T, N> >, N> () (runs, fast_tag ());
 383 #endif
 384
 385 #ifdef USE_UNBOUNDED_ARRAY
 386     header ("matrix<unbounded_array>, vector<unbounded_array> safe");
 387     bench_my_matrix_vector_prod<ublas::matrix<T, ublas::row_major, ublas::unbounded_array<T> >,
 388                                 ublas::vector<T, ublas::unbounded_array<T> >, N> () (runs, safe_tag ());
 389
 390     header ("matrix<unbounded_array>, vector<unbounded_array> fast");
 391     bench_my_matrix_vector_prod<ublas::matrix<T, ublas::row_major, ublas::unbounded_array<T> >,
 392                                 ublas::vector<T, ublas::unbounded_array<T> >, N> () (runs, fast_tag ());
 393 #endif
 394
 395 #ifdef USE_STD_VALARRAY
 396     header ("matrix<std::valarray>, vector<std::valarray> safe");
 397     bench_my_matrix_vector_prod<ublas::matrix<T, ublas::row_major, std::valarray<T> >,
 398                                 ublas::vector<T, std::valarray<T> >, N> () (runs, safe_tag ());
 399
 400     header ("matrix<std::valarray>, vector<std::valarray> fast");
 401     bench_my_matrix_vector_prod<ublas::matrix<T, ublas::row_major, std::valarray<T> >,
 402                                 ublas::vector<T, std::valarray<T> >, N> () (runs, fast_tag ());
 403 #endif
 404
 405 #ifdef USE_STD_VECTOR
 406     header ("matrix<std::vector>, vector<std::vector> safe");
 407     bench_my_matrix_vector_prod<ublas::matrix<T, ublas::row_major, std::vector<T> >,
 408                                 ublas::vector<T, std::vector<T> >, N> () (runs, safe_tag ());
 409
 410     header ("matrix<std::vector>, vector<std::vector> fast");
 411     bench_my_matrix_vector_prod<ublas::matrix<T, ublas::row_major, std::vector<T> >,
 412                                 ublas::vector<T, std::vector<T> >, N> () (runs, fast_tag ());
 413 #endif
 414
 415 #ifdef USE_STD_VALARRAY
 416     header ("std::valarray");
 417     bench_cpp_matrix_vector_prod<std::valarray<T>, std::valarray<T>, N> () (runs);
 418 #endif
 419
 420     header ("matrix + matrix");
 421
 422     header ("C array");
 423     bench_c_matrix_add<T, N> () (runs);
 424
 425 #ifdef USE_C_ARRAY
 426     header ("c_matrix safe");
 427     bench_my_matrix_add<ublas::c_matrix<T, N, N>, N> () (runs, safe_tag ());
 428
 429     header ("c_matrix fast");
 430     bench_my_matrix_add<ublas::c_matrix<T, N, N>, N> () (runs, fast_tag ());
 431 #endif
 432
 433 #ifdef USE_BOUNDED_ARRAY
 434     header ("matrix<bounded_array> safe");
 435     bench_my_matrix_add<ublas::matrix<T, ublas::row_major, ublas::bounded_array<T, N * N> >, N> () (runs, safe_tag ());
 436
 437     header ("matrix<bounded_array> fast");
 438     bench_my_matrix_add<ublas::matrix<T, ublas::row_major, ublas::bounded_array<T, N * N> >, N> () (runs, fast_tag ());
 439 #endif
 440
 441 #ifdef USE_UNBOUNDED_ARRAY
 442     header ("matrix<unbounded_array> safe");
 443     bench_my_matrix_add<ublas::matrix<T, ublas::row_major, ublas::unbounded_array<T> >, N> () (runs, safe_tag ());
 444
 445     header ("matrix<unbounded_array> fast");
 446     bench_my_matrix_add<ublas::matrix<T, ublas::row_major, ublas::unbounded_array<T> >, N> () (runs, fast_tag ());
 447 #endif
 448
 449 #ifdef USE_STD_VALARRAY
 450     header ("matrix<std::valarray> safe");
 451     bench_my_matrix_add<ublas::matrix<T, ublas::row_major, std::valarray<T> >, N> () (runs, safe_tag ());
 452
 453     header ("matrix<std::valarray> fast");
 454     bench_my_matrix_add<ublas::matrix<T, ublas::row_major, std::valarray<T> >, N> () (runs, fast_tag ());
 455 #endif
 456
 457 #ifdef USE_STD_VECTOR
 458     header ("matrix<std::vector> safe");
 459     bench_my_matrix_add<ublas::matrix<T, ublas::row_major, std::vector<T> >, N> () (runs, safe_tag ());
 460
 461     header ("matrix<std::vector> fast");
 462     bench_my_matrix_add<ublas::matrix<T, ublas::row_major, std::vector<T> >, N> () (runs, fast_tag ());
 463 #endif
 464
 465 #ifdef USE_STD_VALARRAY
 466     header ("std::valarray");
 467     bench_cpp_matrix_add<std::valarray<T>, N> () (runs);
 468 #endif
 469 }
 470
 471 #ifdef USE_FLOAT
 472 template struct bench_2<float, 3>;
 473 template struct bench_2<float, 10>;
 474 template struct bench_2<float, 30>;
 475 template struct bench_2<float, 100>;
 476 #endif
 477
 478 #ifdef USE_DOUBLE
 479 template struct bench_2<double, 3>;
 480 template struct bench_2<double, 10>;
 481 template struct bench_2<double, 30>;
 482 template struct bench_2<double, 100>;
 483 #endif
 484
 485 #ifdef USE_STD_COMPLEX
 486 #ifdef USE_FLOAT
 487 template struct bench_2<std::complex<float>, 3>;
 488 template struct bench_2<std::complex<float>, 10>;
 489 template struct bench_2<std::complex<float>, 30>;
 490 template struct bench_2<std::complex<float>, 100>;
 491 #endif
 492
 493 #ifdef USE_DOUBLE
 494 template struct bench_2<std::complex<double>, 3>;
 495 template struct bench_2<std::complex<double>, 10>;
 496 template struct bench_2<std::complex<double>, 30>;
 497 template struct bench_2<std::complex<double>, 100>;
 498 #endif
 499 #endif