Matrix multiply code

[christopher diggins]

Here is some code I wrote for Matrix multiplication for arbitrary
dimensionality known at compile-time. I am curious how practical it is. For
instance, is it common to know the dimensionality of matricies at
compile-time? Any help would be appreciated. Hopefully this code comes in
useful for someone, let me know if you find it useful, or if you have
suggestions on how to improve it.

// Public Domain by Christopher Diggins, 2005

#include <iostream>
#include <numeric>
#include <valarray>
#include <cassert>

template<class Value_T, unsigned int Size_N, unsigned int Stride_N>
class Slice
{
public:
// constructor
Slice(Value_T* x) : p(x) { }
// typedef
typedef Value_T value_type;
// only a forward iterator is provided,
// others are left as an exercise for the reader
struct Slice_forward_i ter {
Slice_forward_i ter(Value_T* x) : p(x) { }
Value_T& operator*() { return *p; }
const Value_T& operator*() const { return *p; }
Value_T& operator++() { Value_T* tmp = p; p += Stride_N; return *tmp; }
Value_T& operator++(int) { return *(p+=Stride_N); }
bool operator==(cons t Slice_forward_i ter& x) { return p == x.p; }
bool operator!=(cons t Slice_forward_i ter& x) { return p != x.p; }
Value_T* p;
};
typedef Slice_forward_i ter iterator;
typedef const Slice_forward_i ter const_iterator;
// public functions
iterator begin() { return iterator(p); }
iterator end() { return iterator(p + (Size_N * Stride_N)); }
value_type& operator[](size_t n) { return *(p + (n * Stride_N)); }
const value_type& operator[](size_t n) const { return *(p + (n *
Stride_N)); }
static size_t size() { return Size_N; }
static size_t stride() { return Stride_N; }
private:
// prevent default construction
Slice() { };
Value_T* p;
};

template<class Value_T, unsigned int Rows_N, unsigned int Cols_N>
class Matrix
{
public:
typedef Slice<Value_T, Cols_N, 1> row_type;
typedef Slice<Value_T, Rows_N, Cols_N> col_type;
const row_type row(unsigned int n) const {
assert(n < rows);
return row_type(data + (n * Cols_N));
}
const col_type column(unsigned int n) const {
assert(n < cols);
return col_type(data + n);
}
row_type operator[](unsigned int n) { return row(n); }
const row_type operator[](unsigned int n) const { return row(n); }
const static unsigned int rows = Rows_N;
const static unsigned int cols = Cols_N;
typedef Value_T value_type;
private:
mutable Value_T data[Rows_N * Cols_N];
};

template<class Matrix1, class Matrix2>
struct MatrixMultiplic ationType {
typedef Matrix<typename Matrix1::value_ type, Matrix1::rows, Matrix2::cols>
type;
};

template<class Slice1_T, class Slice2_T>
typename Slice1_T::value _type dot_product(Sli ce1_T x, Slice2_T y) {
assert(x.size() == y.size());
return std::inner_prod uct(x.begin(), x.end(), y.begin(), typename
Slice1_T::value _type(0));
}

template<class Matrix1, class Matrix2>
typename MatrixMultiplic ationType<Matri x1, Matrix2>::type
matrix_multiply (const Matrix1& m1, const Matrix2& m2)
{
typename MatrixMultiplic ationType<Matri x1, Matrix2>::type result;
assert(Matrix1: :cols == Matrix2::rows);
for (int i=0; i < Matrix1::rows; ++i)
for (int j=0; j < Matrix2::cols; ++j)
result[i][j] = dot_product(m1. row(i), m2.column(j));
return result;
}

template<typena me Matrix_T>
void output_matrix(c onst Matrix_T& x) {
for (int i=0; i < x.rows; ++i) {
for (int j=0; j < x.cols; ++j) {
std::cout << x[i][j] << ", ";
}
std::cout << std::endl;
}
std::cout << std::endl;
}

void test_matrix() {
Matrix<int, 1, 2> m1;
m1[0][0] = 2;
m1[0][1] = 3;
Matrix<int, 2, 2> m2;
m2[0][0] = 2;
m2[0][1] = 0;
m2[1][0] = 0;
m2[1][1] = 2;
output_matrix(m 2);
Matrix<int, 1, 2> m3 = matrix_multiply (m1, m2);
output_matrix(m 3);
}

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

--
Christopher Diggins
Object Oriented Template Library (OOTL)

http://www.ootl.org

[Donovan Rebbechi]

Re: Matrix multiply code

On 2005-05-15, christopher diggins <cdiggins@video tron.ca> wrote:[color=blue]
> Here is some code I wrote for Matrix multiplication for arbitrary
> dimensionality known at compile-time. I am curious how practical it is. For
> instance, is it common to know the dimensionality of matricies at
> compile-time?[/color]

There are some matrices where this is known. It depends on what the matrix
models. If the matrix models some transform, e.g. a 3d coordinate transform,
or a linear transform between two known lists of variables, then the
dimensions may well be known at compile time.

But if the matrix is a data set of some sort (example: some sort of collection
of samples or measurements), the dimensions probably will not be known at
compile time.

Cheers,
--
Donovan Rebbechi

http://pegasus.rutgers.edu/~elflord/

[E. Robert Tisdale]

Re: Matrix multiply code

christopher diggins wrote:
[color=blue]
> Here is some code I wrote for Matrix multiplication for arbitrary
> dimensionality known at compile-time. I am curious how practical it is. For
> instance, is it common to know the dimensionality of matricies at
> compile-time? Any help would be appreciated. Hopefully this code comes in
> useful for someone, let me know if you find it useful, or if you have
> suggestions on how to improve it.
>
> // Public Domain by Christopher Diggins, 2005[/color]

Unless you claim the copyright,
You have no right to give your code away
or even to use it yourself.
For more details, see:

http://www.gnu.org/licenses/licenses.html

You will probably want to use a copyleft:

http://www.gnu.org/licenses/licenses...WhatIsCopyleft[color=blue]
>
> #include <iostream>
> #include <numeric>
> #include <valarray>
> #include <cassert>
>
> template<class Value_T, unsigned int Size_N, unsigned int Stride_N>
> class Slice
> {
> public:
> // constructor
> Slice(Value_T* x) : p(x) { }
> // typedef
> typedef Value_T value_type;
> // only a forward iterator is provided,
> // others are left as an exercise for the reader
> struct Slice_forward_i ter {
> Slice_forward_i ter(Value_T* x) : p(x) { }
> Value_T& operator*() { return *p; }
> const Value_T& operator*() const { return *p; }
> Value_T& operator++() { Value_T* tmp = p; p += Stride_N; return *tmp; }
> Value_T& operator++(int) { return *(p+=Stride_N); }
> bool operator==(cons t Slice_forward_i ter& x) { return p == x.p; }
> bool operator!=(cons t Slice_forward_i ter& x) { return p != x.p; }
> Value_T* p;
> };
> typedef Slice_forward_i ter iterator;
> typedef const Slice_forward_i ter const_iterator;
> // public functions
> iterator begin() { return iterator(p); }
> iterator end() { return iterator(p + (Size_N * Stride_N)); }
> value_type& operator[](size_t n) { return *(p + (n * Stride_N)); }
> const value_type& operator[](size_t n) const { return *(p + (n *
> Stride_N)); }
> static size_t size() { return Size_N; }
> static size_t stride() { return Stride_N; }
> private:
> // prevent default construction
> Slice() { };
> Value_T* p;
> };
>
> template<class Value_T, unsigned int Rows_N, unsigned int Cols_N>
> class Matrix
> {
> public:
> typedef Slice<Value_T, Cols_N, 1> row_type;
> typedef Slice<Value_T, Rows_N, Cols_N> col_type;
> const row_type row(unsigned int n) const {
> assert(n < rows);
> return row_type(data + (n * Cols_N));
> }
> const col_type column(unsigned int n) const {
> assert(n < cols);
> return col_type(data + n);
> }
> row_type operator[](unsigned int n) { return row(n); }
> const row_type operator[](unsigned int n) const { return row(n); }
> const static unsigned int rows = Rows_N;
> const static unsigned int cols = Cols_N;
> typedef Value_T value_type;
> private:
> mutable Value_T data[Rows_N * Cols_N];
> };
>
> template<class Matrix1, class Matrix2>
> struct MatrixMultiplic ationType {
> typedef Matrix<typename Matrix1::value_ type, Matrix1::rows, Matrix2::cols>
> type;
> };
>
> template<class Slice1_T, class Slice2_T>
> typename Slice1_T::value _type dot_product(Sli ce1_T x, Slice2_T y) {
> assert(x.size() == y.size());
> return std::inner_prod uct(x.begin(), x.end(), y.begin(), typename
> Slice1_T::value _type(0));
> }
>
> template<class Matrix1, class Matrix2>
> typename MatrixMultiplic ationType<Matri x1, Matrix2>::type
> matrix_multiply (const Matrix1& m1, const Matrix2& m2)
> {
> typename MatrixMultiplic ationType<Matri x1, Matrix2>::type result;
> assert(Matrix1: :cols == Matrix2::rows);
> for (int i=0; i < Matrix1::rows; ++i)
> for (int j=0; j < Matrix2::cols; ++j)
> result[i][j] = dot_product(m1. row(i), m2.column(j));
> return result;
> }
>
> template<typena me Matrix_T>
> void output_matrix(c onst Matrix_T& x) {
> for (int i=0; i < x.rows; ++i) {
> for (int j=0; j < x.cols; ++j) {
> std::cout << x[i][j] << ", ";
> }
> std::cout << std::endl;
> }
> std::cout << std::endl;
> }
>
> void test_matrix() {
> Matrix<int, 1, 2> m1;
> m1[0][0] = 2;
> m1[0][1] = 3;
> Matrix<int, 2, 2> m2;
> m2[0][0] = 2;
> m2[0][1] = 0;
> m2[1][0] = 0;
> m2[1][1] = 2;
> output_matrix(m 2);
> Matrix<int, 1, 2> m3 = matrix_multiply (m1, m2);
> output_matrix(m 3);
> }
>
> int main() {
> test_matrix();
> system("pause") ;
> return 0;
> }
>[/color]
[color=blue]
> g++ -Wall -ansi -pedantic -o main main.cpp[/color]
main.cpp: In function `void output_matrix(c onst Matrix_T&) [with
Matrix_T = Matrix<int, 2u, 2u>]':
main.cpp:112: instantiated from here
main.cpp:94: warning: comparison between signed and unsigned integer
expressionsmain .cpp:112: instantiated from here
main.cpp:95: warning: comparison between signed and unsigned integer
expressionsmain .cpp: In function `typename
MatrixMultiplic ationType<Matri x1, Matrix2>::type matrix_multiply (const
Matrix1&, const Matrix2&) [with Matrix1 = Matrix<int, 1u, 2u>, Matrix2 =
Matrix<int, 2u, 2u>]':
main.cpp:113: instantiated from here
main.cpp:86: warning: comparison between signed and unsigned integer
expressionsmain .cpp:113: instantiated from here
main.cpp:87: warning: comparison between signed and unsigned integer
expressionsmain .cpp: In function `void output_matrix(c onst Matrix_T&)
[with Matrix_T = Matrix<int, 1u, 2u>]':
main.cpp:114: instantiated from here
main.cpp:94: warning: comparison between signed and unsigned integer
expressionsmain .cpp:114: instantiated from here
main.cpp:95: warning: comparison between signed and unsigned integer
expressions

Evidently, your code needs lots of work.
Also, the design is still very poor.
See The Object-Oriented Numerics Page

http://www.oonumerics.org/oon/

to get an idea about how the experts do this.

Good Luck.

[christopher diggins]

Re: Matrix multiply code

"E. Robert Tisdale" <E.Robert.Tisda le@jpl.nasa.gov > wrote in message
news:d66atc$35p $1@nntp1.jpl.na sa.gov...[color=blue]
> christopher diggins wrote:
>[color=green]
>> Here is some code I wrote for Matrix multiplication for arbitrary
>> dimensionality known at compile-time. I am curious how practical it is.
>> For instance, is it common to know the dimensionality of matricies at
>> compile-time? Any help would be appreciated. Hopefully this code comes in
>> useful for someone, let me know if you find it useful, or if you have
>> suggestions on how to improve it.
>>
>> // Public Domain by Christopher Diggins, 2005[/color]
>
> Unless you claim the copyright,[/color]

In many countries (including Canada and the US), copyright doesn't have to
be claimed to be attributed.
see http://en.wikipedia.org/wiki/Copyrig...yright_notices and

http://www.templetons.com/brad/copymyths.html

[color=blue]
> You have no right to give your code away
> or even to use it yourself.
> For more details, see:
>
> http://www.gnu.org/licenses/licenses.html
>
> You will probably want to use a copyleft:
>
> http://www.gnu.org/licenses/licenses...WhatIsCopyleft[/color]

No thank you, I would prefer contributing something to the public domain
rather than restricting usage.
[color=blue][color=green]
> > g++ -Wall -ansi -pedantic -o main main.cpp[/color]
> main.cpp: In function `void output_matrix(c onst Matrix_T&) [with Matrix_T
> = Matrix<int, 2u, 2u>]':
> main.cpp:112: instantiated from here
> main.cpp:94: warning: comparison between signed and unsigned integer
> expressionsmain .cpp:112: instantiated from here
> main.cpp:95: warning: comparison between signed and unsigned integer
> expressionsmain .cpp: In function `typename
> MatrixMultiplic ationType<Matri x1, Matrix2>::type matrix_multiply (const
> Matrix1&, const Matrix2&) [with Matrix1 = Matrix<int, 1u, 2u>, Matrix2 =
> Matrix<int, 2u, 2u>]':
> main.cpp:113: instantiated from here
> main.cpp:86: warning: comparison between signed and unsigned integer
> expressionsmain .cpp:113: instantiated from here
> main.cpp:87: warning: comparison between signed and unsigned integer
> expressionsmain .cpp: In function `void output_matrix(c onst Matrix_T&)
> [with Matrix_T = Matrix<int, 1u, 2u>]':
> main.cpp:114: instantiated from here
> main.cpp:94: warning: comparison between signed and unsigned integer
> expressionsmain .cpp:114: instantiated from here
> main.cpp:95: warning: comparison between signed and unsigned integer
> expressions
>
> Evidently, your code needs lots of work.[/color]

Simply because I overlooked a few comparison between signed and unsigned
integers warnings? I don't doubt my code has errors and needs work, but the
pedantic warnings from GCC do not indicate where work is actually needed. I
do know that you do have much more to offer me than that ;-)
[color=blue]
> Also, the design is still very poor.[/color]

Could you be more specific? I do appreciate criticism, but I need something
more tangible.
[color=blue]
> See The Object-Oriented Numerics Page
>
> http://www.oonumerics.org/oon/
>
> to get an idea about how the experts do this.[/color]

Thank you, I have looked closely at several libraries, and my original
questions still stand.
[color=blue]
> Good Luck.[/color]

Thank you.

- Christopher Diggins

[E. Robert Tisdale]

Re: Matrix multiply code

christopher diggins wrote:
[color=blue]
> E. Robert Tisdale wrote:
>[color=green]
>>christopher diggins wrote:
>>
>>[color=darkred]
>>>Here is some code I wrote for Matrix multiplication for arbitrary
>>>dimensionali ty known at compile-time. I am curious how practical it is.
>>>For instance, is it common to know the dimensionality of matricies at
>>>compile-time? Any help would be appreciated. Hopefully this code comes in
>>>useful for someone, let me know if you find it useful, or if you have
>>>suggestion s on how to improve it.
>>>
>>>// Public Domain by Christopher Diggins, 2005[/color]
>>
>>Unless you claim the copyright,[/color]
>
>
> In many countries (including Canada and the US), copyright doesn't have to
> be claimed to be attributed.
> see http://en.wikipedia.org/wiki/Copyrig...yright_notices and
> http://www.templetons.com/brad/copymyths.html
>
>[color=green]
>>You have no right to give your code away
>>or even to use it yourself.
>>For more details, see:
>>
>>http://www.gnu.org/licenses/licenses.html
>>
>>You will probably want to use a copyleft:
>>
>>http://www.gnu.org/licenses/licenses...WhatIsCopyleft[/color]
>
>
> No thank you, I would prefer contributing something to the public domain
> rather than restricting usage.
>
>[color=green][color=darkred]
>>>g++ -Wall -ansi -pedantic -o main main.cpp[/color]
>>
>>main.cpp: In function `void output_matrix(c onst Matrix_T&) [with Matrix_T
>>= Matrix<int, 2u, 2u>]':
>>main.cpp:11 2: instantiated from here
>>main.cpp:94 : warning: comparison between signed and unsigned integer
>>expressionsma in.cpp:112: instantiated from here
>>main.cpp:95 : warning: comparison between signed and unsigned integer
>>expressionsma in.cpp: In function `typename
>>MatrixMultipl icationType<Mat rix1, Matrix2>::type matrix_multiply (const
>>Matrix1&, const Matrix2&) [with Matrix1 = Matrix<int, 1u, 2u>, Matrix2 =
>>Matrix<int, 2u, 2u>]':
>>main.cpp:11 3: instantiated from here
>>main.cpp:86 : warning: comparison between signed and unsigned integer
>>expressionsma in.cpp:113: instantiated from here
>>main.cpp:87 : warning: comparison between signed and unsigned integer
>>expressionsma in.cpp: In function `void output_matrix(c onst Matrix_T&)
>>[with Matrix_T = Matrix<int, 1u, 2u>]':
>>main.cpp:11 4: instantiated from here
>>main.cpp:94 : warning: comparison between signed and unsigned integer
>>expressionsma in.cpp:114: instantiated from here
>>main.cpp:95 : warning: comparison between signed and unsigned integer
>>expressions
>>
>>Evidently, your code needs lots of work.[/color]
>
>
> Simply because I overlooked a few comparison between signed and unsigned
> integers warnings? I don't doubt my code has errors and needs work, but the
> pedantic warnings from GCC do not indicate where work is actually needed. I
> do know that you do have much more to offer me than that ;-)
>
>[color=green]
>>Also, the design is still very poor.[/color]
>
>
> Could you be more specific? I do appreciate criticism, but I need something
> more tangible.
>
>[color=green]
>>See The Object-Oriented Numerics Page
>>
>>http://www.oonumerics.org/oon/
>>
>>to get an idea about how the experts do this.[/color]
>
>
> Thank you, I have looked closely at several libraries, and my original
> questions still stand.
>
>[color=green]
>>Good Luck.[/color]
>
>
> Thank you.
>
> - Christopher Diggins
>
>[/color]

Take a look at the blitz++ Tiny Vector Matrix library

http://www.oonumerics.org/blitz/

You can also look at the Array classes in
The C++ Scalar, Vector, Matrix and Tensor class Library

http://www.netwood.net/~edwin/svmtl/

I also have an old recursive template definitions
for Matrix classes that you can get from

ftp.cs.ucla.edu/pub/Tensor.tar.Z

All of these things exist
mainly to help answer the questions that you have asked.

Strictly speaking, your questions are off-topic in this forum
but probably dead on in the object oriented numerics mailing list

http://www.oonumerics.org/mailman/listinfo.cgi/oon-list/

Try there. I don't think that you'll be disappointed.

[christopher diggins]

Re: Matrix multiply code

> Evidently, your code needs lots of work.[color=blue]
> Also, the design is still very poor.
> See The Object-Oriented Numerics Page
>
> http://www.oonumerics.org/oon/
>
> to get an idea about how the experts do this.[/color]

I was curious how an "expert" implementation compares, so I compared my code
to boost::ublas. Preliminary tests show that my implementation runs over
twice as fast, and I haven't done any profiling or optimizing whatsoever.
This makes sense, because the compiler can optimize the indexing using
constants rather than doing variable lookups. This still leaves my original
question: How practical is it? And do people often no matrix dimensions at
compile time?

For those interested here is my first attempt at benchmarking.:

#include <boost/progress.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>

// ... see previous post ...

using namespace boost::numeric: :ublas;

template<class Number, unsigned int M, unsigned int R, unsigned int N,
unsigned int I>
void compare_naive_a nd_ublas()
{
Matrix<Number, M, N> result1;
matrix<Number> result2;
{
std::cout << "about to run Naive test " << std::endl;
Matrix<Number, M, R> m1;
Matrix<Number, R, N> m2;
int n = 0;
for (unsigned i = 0; i < M; ++i)
for (unsigned j = 0; j < R; ++j)
m1[i][j] = ++n;
n = 0;
for (unsigned i = 0; i < R; ++i)
for (unsigned j = 0; j < N; ++j)
m2[i][j] = ++n;
boost::progress _timer t;
for (int i=0; i < I; i++) {
result1 = matrix_multiply (m1, m2);
}
std::cout << "naive time elapsed " << t.elapsed() << std::endl;
}
{
std::cout << "about to run uBlas test " << std::endl;
matrix<Number> m1(M, R);
matrix<Number> m2(R, N);
int n = 0;
for (unsigned i = 0; i < M; ++i)
for (unsigned j = 0; j < R; ++j)
m1(i, j) = ++n;
n = 0;
for (unsigned i = 0; i < R; ++i)
for (unsigned j = 0; j < N; ++j)
m2(i, j) = ++n;
boost::progress _timer t;
for (int i=0; i < I; i++) {
result2 = prod(m1, m2);
}
std::cout << "boost time elapsed " << t.elapsed() << std::endl;
}
for (unsigned i = 0; i < M; ++i)
for (unsigned j = 0; j < N; ++j)
assert(result1[i][j] == result2(i, j));
}

int main() {
compare_naive_a nd_ublas<int,17 ,23,31,1000>();
system("pause") ;
return 0;
}

[christopher diggins]

Re: Matrix multiply code

"christophe r diggins" <cdiggins@video tron.ca> wrote in message
news:bJxhe.1019 16$Tz3.1530622@ weber.videotron .net...[color=blue]
> Here is some code I wrote for Matrix multiplication for arbitrary
> dimensionality known at compile-time. I am curious how practical it is.
> For instance, is it common to know the dimensionality of matricies at
> compile-time? Any help would be appreciated. Hopefully this code comes in
> useful for someone, let me know if you find it useful, or if you have
> suggestions on how to improve it.[/color]

I should be saying rows and columns, not dimensionality.

Christopher Diggins

[christopher diggins]

Re: Matrix multiply code

> Take a look at the blitz++ Tiny Vector Matrix library[color=blue]
>
> http://www.oonumerics.org/blitz/[/color]

This only works for very small matricies, my method works for very
arbitrarily large matricies.
[color=blue]
> You can also look at the Array classes in
> The C++ Scalar, Vector, Matrix and Tensor class Library
>
> http://www.netwood.net/~edwin/svmtl/[/color]

I don't see how this code tracks the rows and columns at compile-time.
[color=blue]
> I also have an old recursive template definitions
> for Matrix classes that you can get from
>
> ftp.cs.ucla.edu/pub/Tensor.tar.Z[/color]

I can't open it.
[color=blue]
> All of these things exist
> mainly to help answer the questions that you have asked.
>
> Strictly speaking, your questions are off-topic in this forum
> but probably dead on in the object oriented numerics mailing list[/color]

Okay let me rephrase it. I have presented a way to do very fast matrix
multiplication in C++ by representing the rows and columns as template
parameters. I use some simple template metaprogramming to compute the type
of the result of the matrix multiplication. The core of the technique is:

template&lt;cla ss Matrix1, class Matrix2>
struct MatrixMultiplic ationType {
typedef Matrix<typename Matrix1::value_ type, Matrix1::rows, Matrix2::cols>
type;
};

template<class Slice1_T, class Slice2_T>
typename Slice1_T::value _type dot_product(Sli ce1_T x, Slice2_T y) {
assert(x.size() == y.size());
return std::inner_prod uct(x.begin(), x.end(), y.begin(), typename
Slice1_T::value _type(0));
}

template<class Matrix1, class Matrix2>
typename MatrixMultiplic ationType<Matri x1, Matrix2>::type
matrix_multiply (const Matrix1& m1, const Matrix2& m2)
{
typename MatrixMultiplic ationType<Matri x1, Matrix2>::type result;
assert(Matrix1: :cols == Matrix2::rows);
for (int i=0; i < Matrix1::rows; ++i)
for (int j=0; j < Matrix2::cols; ++j)
result[i][j] = dot_product(m1. row(i), m2.column(j));
return result;
}

My question is: am I duplicating prior work and does anyone find this
interesting?
[color=blue]
> http://www.oonumerics.org/mailman/li....cgi/oon-list/
>
> Try there. I don't think that you'll be disappointed.[/color]

Thanks.

--
Christopher Diggins
Object Oriented Template Library (OOTL)

http://www.ootl.org

[Mark P]

Re: Matrix multiply code

christopher diggins wrote:[color=blue]
> Here is some code I wrote for Matrix multiplication for arbitrary
> dimensionality known at compile-time. I am curious how practical it is. For
> instance, is it common to know the dimensionality of matricies at
> compile-time? Any help would be appreciated. Hopefully this code comes in
> useful for someone, let me know if you find it useful, or if you have
> suggestions on how to improve it.
>[/color]

You might want to look into Strassen's algorithm for matrix
multiplication which is asymptotically faster than the conventional
(row.column) approach.

[Ufit]

Re: Matrix multiply code


"christophe r diggins" <cdiggins@video tron.ca> wrote in message news:bJxhe.1019 16$Tz3.1530622@ weber.videotron .net...[color=blue]
> Here is some code I wrote for Matrix multiplication for arbitrary
> dimensionality known at compile-time. I am curious how practical it is. For[/color]
---------

BTW if I may interfere - now mostly they do it in assembly language to
boost the performance. I'd like to see some bench comparisons for
normal C++ algorithms because I'd say there's not much difference.

UF

[christopher diggins]

Re: Matrix multiply code (comparing to boost::ublas)

"Ufit" <kot_tmp0SPAMSP AM@NOpoczta.fm> wrote in message
news:d67ifs$qb$ 1@news.onet.pl. ..[color=blue]
>
> "christophe r diggins" <cdiggins@video tron.ca> wrote in message
> news:bJxhe.1019 16$Tz3.1530622@ weber.videotron .net...[color=green]
>> Here is some code I wrote for Matrix multiplication for arbitrary
>> dimensionality known at compile-time. I am curious how practical it is.
>> For[/color]
> ---------
>
> BTW if I may interfere[/color]

Of course! :-)
[color=blue]
> - now mostly they do it in assembly language to
> boost the performance. I'd like to see some bench comparisons for
> normal C++ algorithms because I'd say there's not much difference.
>
> UF[/color]

Here are the results for 1000 multiplications of a 17x23 matrix by a 23x31
matrix of integers by my code compared to ublas. Compiled on Visual C++ 7.1,

about to run Naive test
1.33 s

about to run uBlas test
3.38 s

(YMMV)

There are probably optimizations left to do on the code. However, the code
is mostly pointer arithmetic, so there might not be that much improvement to
be gained by doing it in assembly, though I don't know yet. Any optimization
suggestions would be welcome.
--
Christopher Diggins
Object Oriented Template Library (OOTL)

http://www.ootl.org

[christopher diggins]

Re: Matrix multiply code

"Mark P" <not@my.real.em ail> wrote in message
news:fuEhe.1709 9$J12.3953@news svr14.news.prod igy.com...[color=blue]
> christopher diggins wrote:[color=green]
>> Here is some code I wrote for Matrix multiplication for arbitrary
>> dimensionality known at compile-time. I am curious how practical it is.
>> For instance, is it common to know the dimensionality of matricies at
>> compile-time? Any help would be appreciated. Hopefully this code comes in
>> useful for someone, let me know if you find it useful, or if you have
>> suggestions on how to improve it.
>>[/color]
>
> You might want to look into Strassen's algorithm for matrix multiplication
> which is asymptotically faster than the conventional (row.column)
> approach.[/color]

That is great stuff, thanks for pointing that out. Even though it is only
applicable to n x n matricies, it can still be leveraged by using partial
specializations . What could be done is something like the following (not
tested):

template<class Matrix1, class Matrix2>
struct MatrixMultiplic ationType {
typedef Matrix<typename Matrix1::value_ type, Matrix1::rows, Matrix2::cols>
type;
};

template<class Matrix1, class Matrix2>
typename MatrixMultiplic ationType<Matri x1, Matrix2>::type
matrix_multiply (const Matrix1& m1, const Matrix2& m2)
{
return matrix_multiply _algorithm_choo ser
<
Matrix1,
Matrix2,
(Matrix1::rows == Matrix1::cols)
&& (Matrix1::cols == Matrix2::rows)
&& (Matrix2::rows == Matrix2::cols)[color=blue]
>[/color]
(m1, m2);
}

template<class Matrix1, class Matrix2, bool UseStrassen_B = false>
typename MatrixMultiplic ationType<Matri x1, Matrix2>::type
matrix_multiply _algorithm_choo ser(const Matrix1& m1, const Matrix2& m2)
{
// normal multiplication algorithm
}

template<class Matrix1, class Matrix2>
typename MatrixMultiplic ationType<Matri x1, Matrix2, true>::type
matrix_multiply _algorithm_choo ser(const Matrix1& m1, const Matrix2& m2)
{
// strassen multiplication algorithm
}

This brings to mind other possibilities for other special cases of
matricies.

--
Christopher Diggins
Object Oriented Template Library (OOTL)

http://www.ootl.org

[Jeff Flinn]

Re: Matrix multiply code


"christophe r diggins" <cdiggins@video tron.ca> wrote in message
news:HhChe.1149 39$Tz3.1666576@ weber.videotron .net...

....
[color=blue]
> My question is: am I duplicating prior work and does anyone find this
> interesting?[/color]

IIRC, Dave Abrahams is currently working on Vector/Matrix Linear Algebra
topics over at boost.

Jeff Flinn

[christopher diggins]

Re: Matrix multiply code

"Mark P" <not@my.real.em ail> wrote in message
news:fuEhe.1709 9$J12.3953@news svr14.news.prod igy.com...[color=blue]
> christopher diggins wrote:[color=green]
>> Here is some code I wrote for Matrix multiplication for arbitrary
>> dimensionality known at compile-time. I am curious how practical it is.
>> For instance, is it common to know the dimensionality of matricies at
>> compile-time? Any help would be appreciated. Hopefully this code comes in
>> useful for someone, let me know if you find it useful, or if you have
>> suggestions on how to improve it.
>>[/color]
>
> You might want to look into Strassen's algorithm for matrix multiplication
> which is asymptotically faster than the conventional (row.column)
> approach.[/color]

Upon further research it is important to point out to the interested reader
that Strassen's is only useful for multiplication of square matricies for n[color=blue]
> 654.[/color]

For reference : http://mathworld.wolfram.com/StrassenFormulas.html

Christopher Diggins

[Mark P]

Re: Matrix multiply code

christopher diggins wrote:[color=blue]
> "Mark P" <not@my.real.em ail> wrote in message[color=green]
>>
>>You might want to look into Strassen's algorithm for matrix multiplication
>>which is asymptotically faster than the conventional (row.column)
>>approach.[/color]
>
>
> Upon further research it is important to point out to the interested reader
> that Strassen's is only useful for multiplication of square matricies for n[color=green]
> > 654.[/color]
>[/color]

It's not quite so bad.

That result of 654 is obtained by assuming that Strassen is applied
recursively all the way down to matrices of size 1. This is actually a
pretty horrible idea :) It's much better to apply Strassen down to a
certain minimum size n0, and then use conventional multiplication for
these small matrices.

From my own calculations, the theoretically optimal choice is to apply
Strassen down to a size of about 16. In practice, due to non-algorithmic
overhead (memory access patterns, extra copying, etc.), this number may
be a few times larger. I implemented this once in Java and found n0
around 45, but my implementation may not have been ideal (perhaps a
certainty, given that I used Java :).

The other point is that the algorithm can be extended to deal with
non-square matrices without great difficulty. Assuming M < N, an MxN
matrix can be viewed as floor(N/M) adjacent MxM square matrices and an
adjacent (N%M)xM "non-square" matrix. You can use Strassen to multiply
all combinations of the square matrices and then deal with the
combinations involving non-square matrices separately (either recurse or
use conventional multiply).

Finally, you may want to consider asking about this on a group with
algorithms in the name. Surely you'll find people with better practical
knowledge of this than I can offer.

-Mark