| // Boost.Units - A C++ library for zero-overhead dimensional analysis and |
| // unit/quantity manipulation and conversion |
| // |
| // Copyright (C) 2003-2008 Matthias Christian Schabel |
| // Copyright (C) 2008 Steven Watanabe |
| // |
| // Distributed under the Boost Software License, Version 1.0. (See |
| // accompanying file LICENSE_1_0.txt or copy at |
| // http://www.boost.org/LICENSE_1_0.txt) |
| |
| /** |
| \file |
| |
| \brief performance.cpp |
| |
| \details |
| Test runtime performance. |
| |
| Output: |
| @verbatim |
| |
| multiplying ublas::matrix<double>(1000, 1000) : 25.03 seconds |
| multiplying ublas::matrix<quantity>(1000, 1000) : 24.49 seconds |
| tiled_matrix_multiply<double>(1000, 1000) : 1.12 seconds |
| tiled_matrix_multiply<quantity>(1000, 1000) : 1.16 seconds |
| solving y' = 1 - x + 4 * y with double: 1.97 seconds |
| solving y' = 1 - x + 4 * y with quantity: 1.84 seconds |
| |
| @endverbatim |
| **/ |
| |
| #define _SCL_SECURE_NO_WARNINGS |
| |
| #include <cstdlib> |
| #include <ctime> |
| #include <algorithm> |
| #include <iostream> |
| #include <iomanip> |
| |
| #include <boost/config.hpp> |
| #include <boost/timer.hpp> |
| #include <boost/utility/result_of.hpp> |
| |
| #ifdef BOOST_MSVC |
| #pragma warning(push) |
| #pragma warning(disable:4267; disable:4127; disable:4244; disable:4100) |
| #endif |
| |
| #include <boost/numeric/ublas/matrix.hpp> |
| |
| #ifdef BOOST_MSVC |
| #pragma warning(pop) |
| #endif |
| |
| #include <boost/units/quantity.hpp> |
| #include <boost/units/systems/si.hpp> |
| #include <boost/units/cmath.hpp> |
| |
| enum { |
| tile_block_size = 24 |
| }; |
| |
| template<class T0, class T1, class Out> |
| void tiled_multiply_carray_inner(T0* first, |
| T1* second, |
| Out* out, |
| int totalwidth, |
| int width2, |
| int height1, |
| int common) { |
| for(int j = 0; j < height1; ++j) { |
| for(int i = 0; i < width2; ++i) { |
| Out value = out[j * totalwidth + i]; |
| for(int k = 0; k < common; ++k) { |
| value += first[k + totalwidth * j] * second[k * totalwidth + i]; |
| } |
| out[j * totalwidth + i] = value; |
| } |
| } |
| } |
| |
| template<class T0, class T1, class Out> |
| void tiled_multiply_carray_outer(T0* first, |
| T1* second, |
| Out* out, |
| int width2, |
| int height1, |
| int common) { |
| std::fill_n(out, width2 * height1, Out()); |
| int j = 0; |
| for(; j < height1 - tile_block_size; j += tile_block_size) { |
| int i = 0; |
| for(; i < width2 - tile_block_size; i += tile_block_size) { |
| int k = 0; |
| for(; k < common - tile_block_size; k += tile_block_size) { |
| tiled_multiply_carray_inner( |
| &first[k + width2 * j], |
| &second[k * width2 + i], |
| &out[j * width2 + i], |
| width2, |
| tile_block_size, |
| tile_block_size, |
| tile_block_size); |
| } |
| tiled_multiply_carray_inner( |
| &first[k + width2 * j], |
| &second[k * width2 + i], |
| &out[j * width2 + i], |
| width2, |
| tile_block_size, |
| tile_block_size, |
| common - k); |
| } |
| int k = 0; |
| for(; k < common - tile_block_size; k += tile_block_size) { |
| tiled_multiply_carray_inner( |
| &first[k + width2 * j], |
| &second[k * width2 + i], |
| &out[j * width2 + i], |
| width2, width2 - i, |
| tile_block_size, |
| tile_block_size); |
| } |
| tiled_multiply_carray_inner( |
| &first[k + width2 * j], |
| &second[k * width2 + i], |
| &out[j * width2 + i], |
| width2, width2 - i, |
| tile_block_size, |
| common - k); |
| } |
| int i = 0; |
| for(; i < width2 - tile_block_size; i += tile_block_size) { |
| int k = 0; |
| for(; k < common - tile_block_size; k += tile_block_size) { |
| tiled_multiply_carray_inner( |
| &first[k + width2 * j], |
| &second[k * width2 + i], |
| &out[j * width2 + i], |
| width2, |
| tile_block_size, |
| height1 - j, |
| tile_block_size); |
| } |
| tiled_multiply_carray_inner( |
| &first[k + width2 * j], |
| &second[k * width2 + i], |
| &out[j * width2 + i], |
| width2, |
| tile_block_size, |
| height1 - j, |
| common - k); |
| } |
| int k = 0; |
| for(; k < common - tile_block_size; k += tile_block_size) { |
| tiled_multiply_carray_inner( |
| &first[k + width2 * j], |
| &second[k * width2 + i], |
| &out[j * width2 + i], |
| width2, |
| width2 - i, |
| height1 - j, |
| tile_block_size); |
| } |
| tiled_multiply_carray_inner( |
| &first[k + width2 * j], |
| &second[k * width2 + i], |
| &out[j * width2 + i], |
| width2, |
| width2 - i, |
| height1 - j, |
| common - k); |
| } |
| |
| enum { max_value = 1000}; |
| |
| template<class F, class T, class N, class R> |
| R solve_differential_equation(F f, T lower, T upper, N steps, R start) { |
| typedef typename F::template result<T, R>::type f_result; |
| T h = (upper - lower) / (1.0*steps); |
| for(N i = N(); i < steps; ++i) { |
| R y = start; |
| T x = lower + h * (1.0*i); |
| f_result k1 = f(x, y); |
| f_result k2 = f(x + h / 2.0, y + h * k1 / 2.0); |
| f_result k3 = f(x + h / 2.0, y + h * k2 / 2.0); |
| f_result k4 = f(x + h, y + h * k3); |
| start = y + h * (k1 + 2.0 * k2 + 2.0 * k3 + k4) / 6.0; |
| } |
| return(start); |
| } |
| |
| using namespace boost::units; |
| |
| //y' = 1 - x + 4 * y |
| struct f { |
| template<class Arg1, class Arg2> struct result; |
| |
| double operator()(const double& x, const double& y) const { |
| return(1.0 - x + 4.0 * y); |
| } |
| |
| boost::units::quantity<boost::units::si::velocity> |
| operator()(const quantity<si::time>& x, |
| const quantity<si::length>& y) const { |
| using namespace boost::units; |
| using namespace si; |
| return(1.0 * meters / second - |
| x * meters / pow<2>(seconds) + |
| 4.0 * y / seconds ); |
| } |
| }; |
| |
| template<> |
| struct f::result<double,double> { |
| typedef double type; |
| }; |
| |
| template<> |
| struct f::result<quantity<si::time>, quantity<si::length> > { |
| typedef quantity<si::velocity> type; |
| }; |
| |
| |
| |
| //y' = 1 - x + 4 * y |
| //y' - 4 * y = 1 - x |
| //e^(-4 * x) * (dy - 4 * y * dx) = e^(-4 * x) * (1 - x) * dx |
| //d/dx(y * e ^ (-4 * x)) = e ^ (-4 * x) (1 - x) * dx |
| |
| //d/dx(y * e ^ (-4 * x)) = e ^ (-4 * x) * dx - x * e ^ (-4 * x) * dx |
| //d/dx(y * e ^ (-4 * x)) = d/dx((-3/16 + 1/4 * x) * e ^ (-4 * x)) |
| //y * e ^ (-4 * x) = (-3/16 + 1/4 * x) * e ^ (-4 * x) + C |
| //y = (-3/16 + 1/4 * x) + C/e ^ (-4 * x) |
| //y = 1/4 * x - 3/16 + C * e ^ (4 * x) |
| |
| //y(0) = 1 |
| //1 = - 3/16 + C |
| //C = 19/16 |
| //y(x) = 1/4 * x - 3/16 + 19/16 * e ^ (4 * x) |
| |
| |
| |
| int main() { |
| boost::numeric::ublas::matrix<double> ublas_result; |
| { |
| boost::numeric::ublas::matrix<double> m1(max_value, max_value); |
| boost::numeric::ublas::matrix<double> m2(max_value, max_value); |
| std::srand(1492); |
| for(int i = 0; i < max_value; ++i) { |
| for(int j = 0; j < max_value; ++j) { |
| m1(i,j) = std::rand(); |
| m2(i,j) = std::rand(); |
| } |
| } |
| std::cout << "multiplying ublas::matrix<double>(" |
| << max_value << ", " << max_value << ") : "; |
| boost::timer timer; |
| ublas_result = (prod(m1, m2)); |
| std::cout << timer.elapsed() << " seconds" << std::endl; |
| } |
| typedef boost::numeric::ublas::matrix< |
| boost::units::quantity<boost::units::si::dimensionless> |
| > matrix_type; |
| matrix_type ublas_resultq; |
| { |
| matrix_type m1(max_value, max_value); |
| matrix_type m2(max_value, max_value); |
| std::srand(1492); |
| for(int i = 0; i < max_value; ++i) { |
| for(int j = 0; j < max_value; ++j) { |
| m1(i,j) = std::rand(); |
| m2(i,j) = std::rand(); |
| } |
| } |
| std::cout << "multiplying ublas::matrix<quantity>(" |
| << max_value << ", " << max_value << ") : "; |
| boost::timer timer; |
| ublas_resultq = (prod(m1, m2)); |
| std::cout << timer.elapsed() << " seconds" << std::endl; |
| } |
| std::vector<double> cresult(max_value * max_value); |
| { |
| std::vector<double> m1(max_value * max_value); |
| std::vector<double> m2(max_value * max_value); |
| std::srand(1492); |
| for(int i = 0; i < max_value * max_value; ++i) { |
| m1[i] = std::rand(); |
| m2[i] = std::rand(); |
| } |
| std::cout << "tiled_matrix_multiply<double>(" |
| << max_value << ", " << max_value << ") : "; |
| boost::timer timer; |
| tiled_multiply_carray_outer( |
| &m1[0], |
| &m2[0], |
| &cresult[0], |
| max_value, |
| max_value, |
| max_value); |
| std::cout << timer.elapsed() << " seconds" << std::endl; |
| } |
| std::vector< |
| boost::units::quantity<boost::units::si::energy> |
| > cresultq(max_value * max_value); |
| { |
| std::vector< |
| boost::units::quantity<boost::units::si::force> |
| > m1(max_value * max_value); |
| std::vector< |
| boost::units::quantity<boost::units::si::length> |
| > m2(max_value * max_value); |
| std::srand(1492); |
| for(int i = 0; i < max_value * max_value; ++i) { |
| m1[i] = std::rand() * boost::units::si::newtons; |
| m2[i] = std::rand() * boost::units::si::meters; |
| } |
| std::cout << "tiled_matrix_multiply<quantity>(" |
| << max_value << ", " << max_value << ") : "; |
| boost::timer timer; |
| tiled_multiply_carray_outer( |
| &m1[0], |
| &m2[0], |
| &cresultq[0], |
| max_value, |
| max_value, |
| max_value); |
| std::cout << timer.elapsed() << " seconds" << std::endl; |
| } |
| for(int i = 0; i < max_value; ++i) { |
| for(int j = 0; j < max_value; ++j) { |
| double diff = |
| std::abs(ublas_result(i,j) - cresult[i * max_value + j]); |
| if(diff > ublas_result(i,j) /1e14) { |
| std::cout << std::setprecision(15) << "Uh Oh. ublas_result(" |
| << i << "," << j << ") = " << ublas_result(i,j) |
| << std::endl |
| << "cresult[" << i << " * " << max_value << " + " |
| << j << "] = " << cresult[i * max_value + j] |
| << std::endl; |
| return(EXIT_FAILURE); |
| } |
| } |
| } |
| { |
| std::vector<double> values(1000); |
| std::cout << "solving y' = 1 - x + 4 * y with double: "; |
| boost::timer timer; |
| for(int i = 0; i < 1000; ++i) { |
| double x = .1 * i; |
| values[i] = solve_differential_equation(f(), 0.0, x, i * 100, 1.0); |
| } |
| std::cout << timer.elapsed() << " seconds" << std::endl; |
| for(int i = 0; i < 1000; ++i) { |
| double x = .1 * i; |
| double value = 1.0/4.0 * x - 3.0/16.0 + 19.0/16.0 * std::exp(4 * x); |
| if(std::abs(values[i] - value) > value / 1e9) { |
| std::cout << std::setprecision(15) << "i = : " << i |
| << ", value = " << value << " approx = " << values[i] |
| << std::endl; |
| return(EXIT_FAILURE); |
| } |
| } |
| } |
| { |
| using namespace boost::units; |
| using namespace si; |
| std::vector<quantity<length> > values(1000); |
| std::cout << "solving y' = 1 - x + 4 * y with quantity: "; |
| boost::timer timer; |
| for(int i = 0; i < 1000; ++i) { |
| quantity<si::time> x = .1 * i * seconds; |
| values[i] = solve_differential_equation( |
| f(), |
| 0.0 * seconds, |
| x, |
| i * 100, |
| 1.0 * meters); |
| } |
| std::cout << timer.elapsed() << " seconds" << std::endl; |
| for(int i = 0; i < 1000; ++i) { |
| double x = .1 * i; |
| quantity<si::length> value = |
| (1.0/4.0 * x - 3.0/16.0 + 19.0/16.0 * std::exp(4 * x)) * meters; |
| if(abs(values[i] - value) > value / 1e9) { |
| std::cout << std::setprecision(15) << "i = : " << i |
| << ", value = " << value << " approx = " |
| << values[i] << std::endl; |
| return(EXIT_FAILURE); |
| } |
| } |
| } |
| } |
