#include <iostream>
#include <cstdlib> // rand
#include <ctime> // time
#include <sys/time.h> // gettimeofday
class Timer
{
private:
timeval t1, t2;
public:
Timer() {}
~Timer() {}
void start() { gettimeofday(&t1, NULL); }
void stop() { gettimeofday(&t2, NULL); }
double elapsedTime() { return (t2.tv_sec - t1.tv_sec)*1000.0 + (t2.tv_usec - t1.tv_usec)/1000; }
};
template<typename T>
class Vector
{
private:
T *data;
const size_t num;
public:
Vector(const size_t num) : num(num) { data = new T[num]; }
~Vector() { delete[] data; }
inline T& operator() (const size_t i) { return data[i]; }
inline const T& operator() (const size_t i) const { return data[i]; }
size_t size() const { return num; }
};
template<typename T>
class Matrix
{
private:
T *data;
const size_t nrows, ncols;
public:
Matrix(const size_t nr, const size_t nc) : nrows(nr), ncols(nc) { data = new T[nrows * ncols]; }
~Matrix() { delete[] data; }
inline T& operator() (const size_t r, const size_t c) { return data[c*nrows + r]; }
inline const T& operator() (const size_t r, const size_t c) const { return data[c*nrows + r]; }
size_t size1() const { return nrows; }
size_t size2() const { return ncols; }
};
inline double rand_double(double min=0.0, double max=1.0)
{
return (max - min) * (static_cast<double>(rand()) / RAND_MAX) + min;
}
int main() {
// seed random number generator
srand( static_cast<unsigned int>(time(NULL)) );
// intialize data
const int m = 1000, n = 40000;
Matrix<double> A(m,n);
Vector<double> B(m);
for(size_t i=0; i<A.size1(); i++) {
B(i) = rand_double();
for(size_t j=0; j<A.size2(); j++) {
A(i,j) = rand_double();
}
}
// measure timing
Timer timer;
timer.start();
// in MATLAB: count = sum(bsxfun(@ne, A, B))
Vector<double> count(n);
#pragma omp parallel for
for(int j=0; j<n; ++j) {
count(j) = 0.0;
for(int i=0; i<m; i++) {
count(j) += (A(i,j) != B(i));
}
}
timer.stop();
// elapsed time in milliseconds
std::cout << "Elapsed time is " << timer.elapsedTime() << " milliseconds." << std::endl;
return 0;
}
I2luY2x1ZGUgPGlvc3RyZWFtPgojaW5jbHVkZSA8Y3N0ZGxpYj4gICAgICAgIC8vIHJhbmQKI2luY2x1ZGUgPGN0aW1lPiAgICAgICAgICAvLyB0aW1lCiNpbmNsdWRlIDxzeXMvdGltZS5oPiAgICAgLy8gZ2V0dGltZW9mZGF5CgpjbGFzcyBUaW1lcgp7CnByaXZhdGU6CiAgICB0aW1ldmFsIHQxLCB0MjsKcHVibGljOgogICAgVGltZXIoKSB7fQogICAgflRpbWVyKCkge30KICAgIHZvaWQgc3RhcnQoKSB7IGdldHRpbWVvZmRheSgmdDEsIE5VTEwpOyB9CiAgICB2b2lkIHN0b3AoKSB7IGdldHRpbWVvZmRheSgmdDIsIE5VTEwpOyB9CiAgICBkb3VibGUgZWxhcHNlZFRpbWUoKSB7IHJldHVybiAodDIudHZfc2VjIC0gdDEudHZfc2VjKSoxMDAwLjAgKyAodDIudHZfdXNlYyAtIHQxLnR2X3VzZWMpLzEwMDA7IH0KfTsKCnRlbXBsYXRlPHR5cGVuYW1lIFQ+CmNsYXNzIFZlY3Rvcgp7CnByaXZhdGU6CiAgICBUICpkYXRhOwogICAgY29uc3Qgc2l6ZV90IG51bTsKcHVibGljOgogICAgVmVjdG9yKGNvbnN0IHNpemVfdCBudW0pIDogbnVtKG51bSkgeyBkYXRhID0gbmV3IFRbbnVtXTsgfQogICAgflZlY3RvcigpIHsgZGVsZXRlW10gZGF0YTsgfQogICAgaW5saW5lIFQmIG9wZXJhdG9yKCkgKGNvbnN0IHNpemVfdCBpKSB7IHJldHVybiBkYXRhW2ldOyB9CiAgICBpbmxpbmUgY29uc3QgVCYgb3BlcmF0b3IoKSAoY29uc3Qgc2l6ZV90IGkpIGNvbnN0IHsgcmV0dXJuIGRhdGFbaV07IH0KICAgIHNpemVfdCBzaXplKCkgY29uc3QgeyByZXR1cm4gbnVtOyB9Cn07Cgp0ZW1wbGF0ZTx0eXBlbmFtZSBUPgpjbGFzcyBNYXRyaXgKewpwcml2YXRlOgogICAgVCAqZGF0YTsKICAgIGNvbnN0IHNpemVfdCBucm93cywgbmNvbHM7CnB1YmxpYzoKICAgIE1hdHJpeChjb25zdCBzaXplX3QgbnIsIGNvbnN0IHNpemVfdCBuYykgOiBucm93cyhuciksIG5jb2xzKG5jKSB7IGRhdGEgPSBuZXcgVFtucm93cyAqIG5jb2xzXTsgfQogICAgfk1hdHJpeCgpIHsgZGVsZXRlW10gZGF0YTsgfQogICAgaW5saW5lIFQmIG9wZXJhdG9yKCkgKGNvbnN0IHNpemVfdCByLCBjb25zdCBzaXplX3QgYykgeyByZXR1cm4gZGF0YVtjKm5yb3dzICsgcl07IH0KICAgIGlubGluZSBjb25zdCBUJiBvcGVyYXRvcigpIChjb25zdCBzaXplX3QgciwgY29uc3Qgc2l6ZV90IGMpIGNvbnN0IHsgcmV0dXJuIGRhdGFbYypucm93cyArIHJdOyB9CiAgICBzaXplX3Qgc2l6ZTEoKSBjb25zdCB7IHJldHVybiBucm93czsgfQogICAgc2l6ZV90IHNpemUyKCkgY29uc3QgeyByZXR1cm4gbmNvbHM7IH0KfTsKCmlubGluZSBkb3VibGUgcmFuZF9kb3VibGUoZG91YmxlIG1pbj0wLjAsIGRvdWJsZSBtYXg9MS4wKQp7CiAgICByZXR1cm4gKG1heCAtIG1pbikgKiAoc3RhdGljX2Nhc3Q8ZG91YmxlPihyYW5kKCkpIC8gUkFORF9NQVgpICsgbWluOwp9CgppbnQgbWFpbigpIHsKICAgIC8vIHNlZWQgcmFuZG9tIG51bWJlciBnZW5lcmF0b3IKICAgIHNyYW5kKCBzdGF0aWNfY2FzdDx1bnNpZ25lZCBpbnQ+KHRpbWUoTlVMTCkpICk7CgogICAgLy8gaW50aWFsaXplIGRhdGEKICAgIGNvbnN0IGludCBtID0gMTAwMCwgbiA9IDQwMDAwOwogICAgTWF0cml4PGRvdWJsZT4gQShtLG4pOwogICAgVmVjdG9yPGRvdWJsZT4gQihtKTsKICAgIGZvcihzaXplX3QgaT0wOyBpPEEuc2l6ZTEoKTsgaSsrKSB7CiAgICAgICAgQihpKSA9IHJhbmRfZG91YmxlKCk7CiAgICAgICAgZm9yKHNpemVfdCBqPTA7IGo8QS5zaXplMigpOyBqKyspIHsKICAgICAgICAgICAgQShpLGopID0gcmFuZF9kb3VibGUoKTsKICAgICAgICB9CiAgICB9CgogICAgLy8gbWVhc3VyZSB0aW1pbmcKICAgIFRpbWVyIHRpbWVyOwogICAgdGltZXIuc3RhcnQoKTsKCiAgICAvLyBpbiBNQVRMQUI6IGNvdW50ID0gc3VtKGJzeGZ1bihAbmUsIEEsIEIpKQogICAgVmVjdG9yPGRvdWJsZT4gY291bnQobik7CiAgICAjcHJhZ21hIG9tcCBwYXJhbGxlbCBmb3IKICAgIGZvcihpbnQgaj0wOyBqPG47ICsraikgewogICAgICAgIGNvdW50KGopID0gMC4wOwogICAgICAgIGZvcihpbnQgaT0wOyBpPG07IGkrKykgewogICAgICAgICAgICBjb3VudChqKSArPSAoQShpLGopICE9IEIoaSkpOwogICAgICAgIH0KICAgIH0KCiAgICB0aW1lci5zdG9wKCk7CgogICAgLy8gZWxhcHNlZCB0aW1lIGluIG1pbGxpc2Vjb25kcwogICAgc3RkOjpjb3V0IDw8ICJFbGFwc2VkIHRpbWUgaXMgIiA8PCB0aW1lci5lbGFwc2VkVGltZSgpIDw8ICIgbWlsbGlzZWNvbmRzLiIgPDwgc3RkOjplbmRsOwoKICAgIHJldHVybiAwOwp9Cg==