473,322 Members | 1,398 Online
Bytes | Software Development & Data Engineering Community
Post Job

Home Posts Topics Members FAQ

Join Bytes to post your question to a community of 473,322 software developers and data experts.

Diggins PDP #5: KMatrix (fixed row and column matricieis)

// Diggins Public Domain Post #5
// Contributed to Public Domain by Christopher Diggins
// This is an efficient matrix class implementation for when rows and
columns are known at compile-time

#include <algorithm>
#include <iterator>
#include <cassert>

//#include "static_assert.hpp"

template<class Value_T, int Size_N, int Step_N>
class kslice_iter
{
public:
// public typedefs
typedef Value_T value_type;
typedef Value_T& reference;
typedef const Value_T& const_reference;
typedef std::size_t size_type;
typedef std::ptrdiff_t difference_type;
typedef kslice_iter self;
// constructors
explicit kslice_iter(Value_T* x) : m(x) {
}
// operators
self& operator++() {
m += step();
return *this;
}
self operator++(int) {
self tmp = *this;
m += step();
return tmp;
}
Value_T& operator[](size_type n) {
assert(n < size());
return m[n * Step_N];
}
const Value_T& operator[](size_type n) const {
assert(n < N);
return m[n * Step_N];
}
static int size() {
return Size_N;
}
static int step() {
return Step_N;
}
Value_T& operator*() {
return *m;
}
friend difference_type operator-(const self& x, const self& y) {
return (y.m - x.m) / step();
}
friend bool operator==(const self& x, const self& y) {
return x.m == y.m;
}
friend bool operator!=(const self& x, const self& y) {
return x.m != y.m;
}
friend bool operator<(const self& x, const self& y) {
return x.m < y.m;
}
private:
// default constructor hidden from view
kslice_iter() { };
Value_T* m;
};

template<int N, class T, class RowIter_T, class ColIter_T>
T dot_product(RowIter_T x, ColIter_T y)
{
T ret(0);
for (int i=N; i > 0; --i) {
ret += *x++ * *y++;
}
return ret;
}

template<class T, class M1, class M2, class M3>
kmatrix_multiply(const M1& m1, const M2& m2, M3& result)
{
//STATIC_ASSERT(M1::cols == M2::rows);
//STATIC_ASSERT(M3::rows == M1::rows);
//STATIC_ASSERT(M3::cols == M2::cols);
for (int i=0; i < M1::rows; ++i) {
for (int j=0; j < M2::cols; ++j) {
result[i][j] = dot_product<M1::cols, T>(m1.row(i), m2.column(j));
}
}
}

template<class Value_T, int Rows_N, int Cols_N>
class kmatrix
{
public:
// public typedefs
typedef Value_T value_type;
typedef kmatrix self;
typedef Value_T* iterator;
typedef const Value_T* const_iterator;
typedef kslice_iter<Value_T, Cols_N, 1> row_type;
typedef kslice_iter<Value_T, Rows_N, Cols_N> col_type;
// public constants
static const int rows = Rows_N;
static const int cols = Cols_N;
// constructors
kmatrix() {
std::fill(begin(), end(), Value_T(0));
}
kmatrix(Value_T& x) {
std::fill(begin(), end(), x);
}
kmatrix(const self& x) {
operator=(x);
}
// public functions
self& operator=(const self& x) {
std::copy(x.begin(), x.end(), begin());
return *this;
}
self& operator=(Value_T& x) {
std::fill(begin(), end(), x);
return *this;
}
row_type row(int n) const {
assert(n < rows);
return row_type(m + (n * Cols_N));
}
col_type column(int n) const {
assert(n < cols);
return col_type(m + n);
}
row_type operator[](int n) {
return row(n);
}
const row_type operator[](int n) const {
return row(n);
}
iterator begin() {
return m;
}
iterator end() {
return m + size();
}
const_iterator begin() const {
return m;
}
const_iterator end() const {
return m + size();
}
static int size() {
return Rows_N * Cols_N;
}
// operators
self& operator+=(const self& x) {
for (int i=0; i < size(); ++i) {
m[i] += x.m[i];
}
return *this;
}
self& operator-=(const self& x) {
for (int i=0; i < size(); ++i) {
m[i] -= x.m[i];
}
return *this;
}
self& operator+=(value_type x) {
for (int i=0; i < size(); ++i) {
m[i] += x;
}
return *this;
}
self& operator-=(value_type x) {
for (int i=0; i < size(); ++i) {
m[i] -= x;
}
return *this;
}
self& operator*=(value_type x) {
for (int i=0; i < size(); ++i) {
m[i] *= x;
}
return *this;
}
self& operator/=(value_type x) {
for (int i=0; i < size(); ++i) {
m[i] /= x;
}
return *this;
}
self& operator=(value_type x) {
std::fill(begin(), end(), x);
return *this;
}
self operator-() {
self x;
for (int i=0; i < size(); ++i) {
x.m[i] = -m[i];
}
return x;
}
// friends
friend self operator+(const self& x, const self& y) {
return self(x) += y;
}
friend self operator-(const self& x, const self& y) {
return self(x) -= y;
}
friend self operator+(const self& x, value_type y) {
return self(x) += y;
}
friend self operator-(const self& x, value_type y) {
return self(x) -= y;
}
friend self operator*(const self& x, value_type y) {
return self(x) *= y;
}
friend self operator/(const self& x, value_type y) {
return self(x) /= y;
}
template<int Cols2_N>
friend kmatrix<Value_T, Rows_N, Cols2_N>
operator*(const self& x, const kmatrix<Value_T, Cols_N, Cols2_N>& y) {
kmatrix<Value_T, Rows_N, Cols2_N> ret;
kmatrix_multiply(x, y, ret);
return ret;
}
private:
mutable Value_T m[Rows_N * Cols_N];
};
Jul 23 '05 #1
1 1370
// public domain matrix benchmark code

#include "kmatrix.hpp"

#include <cassert>
#include <iostream>
#include <cstdlib>
#include <boost/numeric/ublas/matrix.hpp>

class timer
{
public:
timer() : start(std::clock()) {
}
double elapsed() const {
return (double(std::clock() - start) / CLOCKS_PER_SEC) * 1000;
}
private:
std::clock_t start;
};

void output_time_elapsed(const timer& t, std::ostream& o = std::cout) {
o << static_cast<int>(t.elapsed()) << " msec elapsed" << std::endl;
}

struct reporting_timer {
~reporting_timer() {
output_time_elapsed(m);
}
timer m;
};
template<int Rows_N, int Cols_N, class Matrix_T>
void init_kmatrix(Matrix_T& m)
{
for (int i=0; i < Rows_N; ++i) {
for (int j=0; j < Cols_N; ++j) {
m[i][j] = i * Cols_N + j;
}
}
}

template<int Rows_N, int Cols_N, class Matrix_T>
void init_ublas_matrix(Matrix_T& m)
{
for (int i=0; i < Rows_N; ++i) {
for (int j=0; j < Cols_N; ++j) {
m(i, j) = i * Cols_N + j;
}
}
}

template<class T, int Iters_N, int M, int R, int N>
void kmatrix_performance_test()
{
std::cout << "kmatrix ";
kmatrix<T, M, R> km1;
kmatrix<T, R, N> km2;
kmatrix<T, M, N> km3;
init_kmatrix<M, R>(km1);
init_kmatrix<R, M>(km2);
{
reporting_timer t;
for (int i=0; i < Iters_N; ++i)
{
kmatrix_multiply<T>(km1, km2, km3);
}
}

std::cout << "ublas ";
boost::numeric::ublas::matrix<T> m1(M,R);
boost::numeric::ublas::matrix<T> m2(R,N);
boost::numeric::ublas::matrix<T> m3(M,N);
init_ublas_matrix<M, R>(m1);
init_ublas_matrix<R, M>(m2);
{
reporting_timer t;
for (int i=0; i < Iters_N; ++i)
{
boost::numeric::ublas::noalias(m3) =
boost::numeric::ublas::prod(m1,m2);
}
}

for (int i=0; i<M; ++i) {
for (int j=0; j<N; ++j) {
T x = m3(i,j);
T y = km3[i][j];
ASSURE(x == y);
}
}
}

void benchmarks()
{
std::cout << "comparing integer matricies" << std::endl;
std::cout << "2,2 X 2,2" << std::endl;
kmatrix_performance_test<int, 100000, 2, 2, 2>();
std::cout << "3,3 X 3,3" << std::endl;
kmatrix_performance_test<int, 100000, 3, 3, 3>();
std::cout << "100,100 X 100,100" << std::endl;
kmatrix_performance_test<int, 10, 100, 100, 100>();
std::cout << "100,1 X 1,100" << std::endl;
kmatrix_performance_test<int, 100, 100, 1, 100>();
std::cout << "1,100 X 100, 1" << std::endl;
kmatrix_performance_test<int, 100000, 1, 100, 1>();

std::cout << "comparing double matricies" << std::endl;
std::cout << "2,2 X 2,2" << std::endl;
kmatrix_performance_test<double, 100000, 2, 2, 2>();
std::cout << "3,3 X 3,3" << std::endl;
kmatrix_performance_test<double, 100000, 3, 3, 3>();
std::cout << "100,100 X 100,100" << std::endl;
kmatrix_performance_test<double, 10, 100, 100, 100>();
std::cout << "100,1 X 1,100" << std::endl;
kmatrix_performance_test<double, 100, 100, 1, 100>();
std::cout << "1,100 X 100, 1" << std::endl;
kmatrix_performance_test<double, 100000, 1, 100, 1>();
}

int main()
{
benchmarks();
system("pause");
return 0;
}

/*
my results:

comparing integer matricies
2,2 X 2,2
kmatrix 78 msec elapsed
ublas 172 msec elapsed
3,3 X 3,3
kmatrix 219 msec elapsed
ublas 437 msec elapsed
100,100 X 100,100
kmatrix 516 msec elapsed
ublas 1156 msec elapsed
100,1 X 1,100
kmatrix 125 msec elapsed
ublas 203 msec elapsed
1,100 X 100, 1
kmatrix 500 msec elapsed
ublas 1172 msec elapsed
comparing double matricies
2,2 X 2,2
kmatrix 79 msec elapsed
ublas 171 msec elapsed
3,3 X 3,3
kmatrix 250 msec elapsed
ublas 422 msec elapsed
100,100 X 100,100
kmatrix 656 msec elapsed
ublas 1344 msec elapsed
100,1 X 1,100
kmatrix 156 msec elapsed
ublas 219 msec elapsed
1,100 X 100, 1
kmatrix 703 msec elapsed
ublas 1172 msec elapsed
Press any key to continue . . .
*/
Jul 23 '05 #2

This thread has been closed and replies have been disabled. Please start a new discussion.

Similar topics

179
by: SoloCDM | last post by:
How do I keep my entire web page at a fixed width? ********************************************************************* Signed, SoloCDM
6
by: Mason A. Clark | last post by:
Masters: On two or three-column layouts, one column often has a list of links. Scrolling the page hides them. I'm aware there's supposed to be the ability to fix the column (frame-like). I...
4
by: christopher diggins | last post by:
Welcome to the first installment of the Diggins Public Domain Posts (PDP). There is a significant dearth of good public domain C++ code. All of it seems to come with one kind of a license or...
1
by: christopher diggins | last post by:
// meta_binder.hpp // The Diggins PDP (Public Domain Post) #2 // Public Domain code by Christopher Diggins, May 22, 2005 // // Description: // A meta-binder binds function objects to function...
0
by: christopher diggins | last post by:
// binary_aritmetic.hpp // Diggins PDP (Public Domain Post) #1.2 // by Christopher Diggins, May 24, 2005 // // Comment: // Provides several algorithms for unsigned binary arithmetic // //...
0
by: christopher diggins | last post by:
// big_int.hpp // The Diggins PDP (Public Domain Post) #3 // Public Domain code by Christopher Diggins, May 22, 2005 // // Description: // A naively implemented unsigned integer type for...
0
by: christopher diggins | last post by:
// The Diggins PDP (Public Domain Post) #4 // public domain code for computing the FFT // contributed by Christopher Diggins, 2005 template<int N> unsigned int bit_reverse(unsigned int x) {...
4
by: Otie | last post by:
Hello, I am using the MSFlexGrd Control in VB5. I have 1 fixed row and one fixed column. I am trying to do a sort when the user clicks a column in the FIXED ROW. But when I capture the row...
4
by: Jeff | last post by:
Hey I'm wondering how the Fixed-Width Text Format is What I know is that the top line in this text format will contain column names. and each row beneath the top line represent for example a...
1
isladogs
by: isladogs | last post by:
The next Access Europe meeting will be on Wednesday 6 Mar 2024 starting at 18:00 UK time (6PM UTC) and finishing at about 19:15 (7.15PM). In this month's session, we are pleased to welcome back...
0
by: Vimpel783 | last post by:
Hello! Guys, I found this code on the Internet, but I need to modify it a little. It works well, the problem is this: Data is sent from only one cell, in this case B5, but it is necessary that data...
0
by: jfyes | last post by:
As a hardware engineer, after seeing that CEIWEI recently released a new tool for Modbus RTU Over TCP/UDP filtering and monitoring, I actively went to its official website to take a look. It turned...
0
by: ArrayDB | last post by:
The error message I've encountered is; ERROR:root:Error generating model response: exception: access violation writing 0x0000000000005140, which seems to be indicative of an access violation...
1
by: CloudSolutions | last post by:
Introduction: For many beginners and individual users, requiring a credit card and email registration may pose a barrier when starting to use cloud servers. However, some cloud server providers now...
1
by: Defcon1945 | last post by:
I'm trying to learn Python using Pycharm but import shutil doesn't work
1
by: Shællîpôpï 09 | last post by:
If u are using a keypad phone, how do u turn on JavaScript, to access features like WhatsApp, Facebook, Instagram....
0
by: Faith0G | last post by:
I am starting a new it consulting business and it's been a while since I setup a new website. Is wordpress still the best web based software for hosting a 5 page website? The webpages will be...
0
isladogs
by: isladogs | last post by:
The next Access Europe User Group meeting will be on Wednesday 3 Apr 2024 starting at 18:00 UK time (6PM UTC+1) and finishing by 19:30 (7.30PM). In this session, we are pleased to welcome former...

By using Bytes.com and it's services, you agree to our Privacy Policy and Terms of Use.

To disable or enable advertisements and analytics tracking please visit the manage ads & tracking page.