Re: [Boost-users] apply std::pow to vector or matrix

Subject: Re: [Boost-users] apply std::pow to vector or matrix
From: Max S. Kaznady (max.kaznady_at_[hidden])
Date: 2010-07-11 17:35:41

I wrote some benchmarks for Boost and GSL where I initialize a vector
of random numbers and then raise it to some random exponent, and I
redo the operations a certain number of times and report the average
time.

std::transform and std::for_each take the same time, GSL's filling in
manually and Boosts's filling in manually also take the same time BUT,
manual fil is a bit faster than using std::transform and
std::for_each. Any ideas why?

Runtimes:
Boost vector initialization:
0.000000000000000
GSL vector initialization:
0.000000000000000
Boost transform avg time:
0.002141000000000
Boost for_each avg time:
0.002094000000000
Boost manual fill avg time:
0.001811000000000
GSL avg time:
0.001820000000000

The code is below:
/*
 * main.cpp
 *
 *  Created on: 2010-07-08
 *      Author: max
 */

#include <iostream>
#include <vector>
#include <algorithm>
#include <math.h>
#include <fstream>
#include <ctime>            // std::time
#include <boost/bind.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/timer.hpp> // time the runs
// For random number generation
#include <boost/random/linear_congruential.hpp>
#include <boost/random/uniform_int.hpp>
#include <boost/random/uniform_real.hpp>
#include <boost/random/variate_generator.hpp>
// Sun CC doesn't handle boost::iterator_adaptor yet
#if !defined(__SUNPRO_CC) || (__SUNPRO_CC > 0x530)
#include <boost/generator_iterator.hpp>
#endif

#ifdef BOOST_NO_STDC_NAMESPACE
namespace std {
  using ::time;
}
#endif

// This is a typedef for a random number generator.
// Try boost::mt19937 or boost::ecuyer1988 instead of boost::minstd_rand
typedef boost::minstd_rand base_generator_type;

// GSL includes
#include <stdio.h>
#include <gsl/gsl_vector.h>

void my_pow(double & val, double exponent) {
  val = std::pow(val, exponent);
}

int main() {
  // Size of assigned vector
  unsigned N = 10000;
  // number of times to repeat the simulation
  // unsigned numrepeat = 10000;
  // exponent to use
  double exponent; // = 12.12345678901234;
  // value to exponentiate
  double fill_val; // = 1.12345678901234;

  // Define a random number generator and initialize it with a reproducible
  // seed.
  // (The seed is unsigned, otherwise the wrong overload may be selected
  // when using mt19937 as the base_generator_type.)
  base_generator_type generator(42u);
  // Define a uniform random number distribution which produces "double"
  // values between 0 and 1 (0 inclusive, 1 exclusive).
  boost::uniform_real<> uni_dist(0, 1);
  boost::variate_generator<base_generator_type&, boost::uniform_real<> > uni(
      generator, uni_dist);
  /*
   * Change seed to something else.
   *
   * Caveat: std::time(0) is not a very good truly-random seed.  When
   * called in rapid succession, it could return the same values, and
   * thus the same random number sequences could ensue.  If not the same
   * values are returned, the values differ only slightly in the
   * lowest bits.  A linear congruential generator with a small factor
   * wrapped in a uniform_smallint (see experiment) will produce the same
   * values for the first few iterations.   This is because uniform_smallint
   * takes only the highest bits of the generator, and the generator itself
   * needs a few iterations to spread the initial entropy from the lowest bits
   * to the whole state.
   */
  // Set the generator later for each piece of code
  //generator.seed(static_cast<unsigned int> (std::time(0)));

  // Vector of times
  boost::numeric::ublas::vector<double> times(N);
  // Boost timer
  boost::timer t;
  std::cout.precision(15);

  t.restart();
  boost::numeric::ublas::vector<double> boost_v(N);
  std::cout << "Boost vector initialization: " << std::endl;
  std::cout << std::fixed << t.elapsed() << std::endl;
  t.restart();
  gsl_vector * gsl_v = gsl_vector_alloc(N);
  std::cout << "GSL vector initialization: " << std::endl;
  std::cout << std::fixed << t.elapsed() << std::endl;

  // Filling Boost using std::transform
  generator.seed(static_cast<unsigned int> (std::time(0)));
  for (unsigned i = 0; i < N; ++i) {
    fill_val = uni();
    exponent = 2.0 + uni();
    t.restart();
    std::fill(boost_v.begin(), boost_v.end(), fill_val);
    std::transform(boost_v.begin(), boost_v.end(), boost_v.begin(),
        boost::bind(pow, _1, exponent));
    times(i) = t.elapsed();
  }
  std::cout << "Boost transform avg time: " << std::endl;
  std::cout << std::fixed << boost::numeric::ublas::sum(times) / (double) N
      << std::endl;

  // Filling Boost using std::for_each
  generator.seed(static_cast<unsigned int> (std::time(0)));
  for (unsigned i = 0; i < N; ++i) {
    fill_val = uni();
    exponent = 2.0 + uni();
    t.restart();
    std::fill(boost_v.begin(), boost_v.end(), fill_val);
    std::for_each(boost_v.begin(), boost_v.end(), boost::bind(my_pow, _1,
        exponent));
    times(i) = t.elapsed();
  }
  std::cout << "Boost for_each avg time: " << std::endl;
  std::cout << std::fixed << boost::numeric::ublas::sum(times) / (double) N
      << std::endl;

  // Filling Boost manually
  generator.seed(static_cast<unsigned int> (std::time(0)));
  for (unsigned i = 0; i < N; ++i) {
    fill_val = uni();
    exponent = 2.0 + uni();
    t.restart();
    for (unsigned j = 0; j<N; ++j) {
      boost_v(j) = std::pow(fill_val, exponent);
    }
    times(i) = t.elapsed();
  }
  std::cout << "Boost manual fill avg time: " << std::endl;
  std::cout << std::fixed << boost::numeric::ublas::sum(times) / (double) N
      << std::endl;

  // Filling GSL by hand
  generator.seed(static_cast<unsigned int> (std::time(0)));
  for (unsigned i = 0; i < N; ++i) {
    fill_val = uni();
    exponent = 2.0 + uni();
    t.restart();
    for (unsigned j = 0; j<N; ++j) {
      gsl_vector_set(gsl_v, j, std::pow(fill_val, exponent));
    }
    times(i) = t.elapsed();
  }
  std::cout << "GSL avg time: " << std::endl;
  std::cout << std::fixed << boost::numeric::ublas::sum(times) / (double) N
      << std::endl;

  //de-allocate GSL
  gsl_vector_free(gsl_v);
  return 0;
}

Cheers,
Max

On Wed, Jul 7, 2010 at 12:13 PM, Robert Jones <robertgbjones_at_[hidden]> wrote:
> On Wed, Jul 7, 2010 at 4:45 PM, Max S. Kaznady <max.kaznady_at_[hidden]>
> wrote:
>>
>> Great example, thanks! Before your post, I wrote something like:
>>
>> namespace my {
>> void my_vec_pow(boost::numeric::ublas::vector<double> &vec, double expon)
>> {
>>  boost::numeric::ublas::vector<double>::iterator it;
>>  for (it = vec.begin(); it != vec.end(); ++it)
>>    *it = std::pow(*it, expon);
>> }
>> }
>>
>> I'm just curious, but from performance perspective, how much more
>> efficient is it to use iterators and functors? I would assume that
>> most compilers can exploit this to improve memory management and speed
>> up the computation... I'll write some benchmarks later on this week,
>> efficiency is crucial for what I'm doing.
>>
>
> By and large, in my experience, it's barely measurable, at least for
> simple cases. When you write
>
> for(i=v.begin(); i!=v.end();....)
>
> v.end() is evaluated on each iteration, which may have some effect,
> depending
> on type of v, but for simple vectors it's probably inlined anyway. (there's
> folks
> on here that know much more about it than I do).
>
> With the new Boost.Range algorithms there's also some syntactic
> simplification
> as you can write
>
> for(v,fn);
>
> which is a real plus if v is complicated expression, since you don't have to
> evaluate
> it twice or even put it in an explicit variable.
>
> Cheers
>
> - Rob.
>
>
> _______________________________________________
> Boost-users mailing list
> Boost-users_at_[hidden]
> http://lists.boost.org/mailman/listinfo.cgi/boost-users
>