Compile-time matrix dimension checking and template friend question

[Ben Ingram]

Hi all,

I am writing a template matrix class in which the template parameters are
the number of rows and number of columns. There are a number of reasons
why this is an appropriate tradeoff for my particular application. One of the
advantages is that the _compiler_ can force inner matrix dimensions used in
multiplication to agree. A _complie-time_ error will be triggered if you
write A * B and the number of coluns in A does not equal the number of
rows in B. Here's simplified code that illustrates the concept:

template<int nRows, int nCols> class Matrix {
public:
double data[nRows][nCols];

// operator*: return this * A
template<int nNewCols> Matrix<nRows,nNewCols>
operator*(const Matrix<nCols,nNewCols> &A) const {
Matrix<nRows,nNewCols> ret;
for(int iRow = 0; iRow < nRows; iRow++) {
for(int iCol = 0; iCol < nNewCols; iCol++) {
double innerProd = 0.0;
for(int iInner = 0; iInner < nCols; iInner++) {
innerProd += data[iRow][iInner] * A.data[iInner][iCol];
}
ret.data[iRow][iCol] = innerProd;
}
}
return ret;
}
};

int main(int argc, char **argv) {
Matrix<4,3> A;
Matrix<3,1> x;
Matrix<2,1> y;

A*x; // compiler creates
// Matrix<4,1> Matrix<4,3>::operator*(const Matrix<3,1>&)
//A*y; // compile-time error thrown since compiler can't create
// Matrix<4,1> Matrix<4,3>::operator*(const Matrix<2,1>&)
// since inner dimensions don't agree

return 0;
}

Note that while this code works under all g++ versions I tested, from
g++-2.7 through g++-3.2, it causes an internal compiler error in MSVC++
6.

The problem with this code is that the data field must be public. Otherwise,
Matrix<4,3>::operator* will not have access to Matrix<4,1>'s or
Matrix<3,1>'s data. Thus, I would like to make Matrix<4,1> and
Matrix<3,1> friends of Matrix<4,3>. Is this possible? If so, what's the
syntax? I tried the code below and some other variations, but they were
all incorrect. The code below generated "partial specialization
`Matrix<nRows, nNewCols>' declared `friend'" and "partial specialization
`Matrix<nCols, nNewCols>' declared `friend'" errors under g++-3.0.4 and
did not give Matrix<4,3> access to the private data of Matrix<4,1> or
of Matrix<3,1>.

template<int nRows, int nCols> class Matrix {
private:
double data[nRows][nCols];

public:
// operator*: return this * A
template<int nNewCols> Matrix<nRows,nNewCols>
operator*(const Matrix<nCols,nNewCols> &A) const {
Matrix<nRows,nNewCols> ret;
for(int iRow = 0; iRow < nRows; iRow++) {
for(int iCol = 0; iCol < nNewCols; iCol++) {
double innerProd = 0.0;
for(int iInner = 0; iInner < nCols; iInner++) {
innerProd += data[iRow][iInner] * A.data[iInner][iCol];
}
ret.data[iRow][iCol] = innerProd;
}
}
return ret;
}

template<int nNewCols> friend class Matrix<nRows,nNewCols>;
template<int nNewCols> friend class Matrix<nCols,nNewCols>;
};

int main(int argc, char **argv) {
Matrix<4,3> A;
Matrix<3,1> x;
Matrix<2,1> y;

A*x; // compiler creates
// Matrix<4,1> Matrix<4,3>::operator*(const Matrix<3,1>&)
//A*y; // compile-time error thrown since compiler can't create
// Matrix<4,1> Matrix<4,3>::operator*(const Matrix<2,1>&)
// since inner dimensions don't agree

return 0;
}

Can this be done and if so, how? Thanks for any help,

Ben :-)

[Chris Theis]

"Ben Ingram" <bt************@catgufu.mit.edu> wrote in message
news:a9******************************@news.sonicne ws.com...

Hi all,

I am writing a template matrix class in which the template parameters are
the number of rows and number of columns. There are a number of reasons
why this is an appropriate tradeoff for my particular application. One of the advantages is that the _compiler_ can force inner matrix dimensions used in multiplication to agree. A _complie-time_ error will be triggered if you
write A * B and the number of coluns in A does not equal the number of
rows in B. Here's simplified code that illustrates the concept:

template<int nRows, int nCols> class Matrix {
public:
double data[nRows][nCols];

// operator*: return this * A
template<int nNewCols> Matrix<nRows,nNewCols>
operator*(const Matrix<nCols,nNewCols> &A) const {
Matrix<nRows,nNewCols> ret;
for(int iRow = 0; iRow < nRows; iRow++) {
for(int iCol = 0; iCol < nNewCols; iCol++) {
double innerProd = 0.0;
for(int iInner = 0; iInner < nCols; iInner++) {
innerProd += data[iRow][iInner] * A.data[iInner][iCol];
}
ret.data[iRow][iCol] = innerProd;
}
}
return ret;
}
};

int main(int argc, char **argv) {
Matrix<4,3> A;
Matrix<3,1> x;
Matrix<2,1> y;

A*x; // compiler creates
// Matrix<4,1> Matrix<4,3>::operator*(const Matrix<3,1>&)
//A*y; // compile-time error thrown since compiler can't create
// Matrix<4,1> Matrix<4,3>::operator*(const Matrix<2,1>&)
// since inner dimensions don't agree

return 0;
}

Note that while this code works under all g++ versions I tested, from
g++-2.7 through g++-3.2, it causes an internal compiler error in MSVC++
6.

The problem with this code is that the data field must be public. Otherwise,
[SNIP]

Having the data field public is not really a good idea and it is not
necessary. I'd recommend to implement the data access via the operator()
like the following for example:

inline double& operator()( unsigned int Row, unsigned int Col ) {
assert( Row < MaxRow && Col < MaxCol ); // you
can store the dimensions at the time of construction!
return data[Row][Col];
}

inline double operator()( int Row, int Col ) const {
assert( Row < MaxRow && Col < MaxCol );
return data[Row][Col];
}

Therefore you can simply write A(3, 5) to access the respective matrix
element. Furthermore the physical way of data management is hidden, because
you might consider to use a 1D array for speed purposes or whatever reason.
However, the problem you're facing with VC++ 6.0 seems to be related to its
shortcomings regarding templates & partial specialization. I think that
using a later version should resolve your problem.

For an easy sample matrix implementation look at:

http://www.parashift.com/c++-faq-lit...html#faq-16.17

Can this be done and if so, how? Thanks for any help,

Ben :-)

HTH
Chris

[Ben Ingram]

On Thu, 29 Jan 2004 01:35:22 -0800, Chris Theis wrote:

"Ben Ingram" <bt************@catgufu.mit.edu> wrote in message
news:a9******************************@news.sonicne ws.com...

Hi all,

I am writing a template matrix class in which the template parameters
are the number of rows and number of columns. There are a number of
reasons why this is an appropriate tradeoff for my particular
application. One of

the

advantages is that the _compiler_ can force inner matrix dimensions
used

in

multiplication to agree. A _complie-time_ error will be triggered if
you write A * B and the number of coluns in A does not equal the number
of rows in B. Here's simplified code that illustrates the concept:

template<int nRows, int nCols> class Matrix { public:
double data[nRows][nCols];

// operator*: return this * A
template<int nNewCols> Matrix<nRows,nNewCols>
operator*(const Matrix<nCols,nNewCols> &A) const {
Matrix<nRows,nNewCols> ret;
for(int iRow = 0; iRow < nRows; iRow++) {
for(int iCol = 0; iCol < nNewCols; iCol++) {
double innerProd = 0.0;
for(int iInner = 0; iInner < nCols; iInner++) {
innerProd += data[iRow][iInner] * A.data[iInner][iCol];
}
ret.data[iRow][iCol] = innerProd;
}
}
return ret;
}
};

int main(int argc, char **argv) {
Matrix<4,3> A;
Matrix<3,1> x;
Matrix<2,1> y;

A*x; // compiler creates
// Matrix<4,1> Matrix<4,3>::operator*(const Matrix<3,1>&)
//A*y; // compile-time error thrown since compiler can't create
// Matrix<4,1> Matrix<4,3>::operator*(const Matrix<2,1>&) //
since inner dimensions don't agree

return 0;
}
}
Note that while this code works under all g++ versions I tested, from
g++-2.7 through g++-3.2, it causes an internal compiler error in MSVC++
6.

The problem with this code is that the data field must be public.

Otherwise,
[SNIP]

Having the data field public is not really a good idea and it is not
necessary. I'd recommend to implement the data access via the operator()
like the following for example:

inline double& operator()( unsigned int Row, unsigned int Col ) {
assert( Row < MaxRow && Col < MaxCol ); //
you
can store the dimensions at the time of construction!
return data[Row][Col];
}
}
inline double operator()( int Row, int Col ) const {
assert( Row < MaxRow && Col < MaxCol ); return data[Row][Col];
}
}
Therefore you can simply write A(3, 5) to access the respective matrix
element. Furthermore the physical way of data management is hidden,
because you might consider to use a 1D array for speed purposes or
whatever reason. However, the problem you're facing with VC++ 6.0 seems
to be related to its shortcomings regarding templates & partial
specialization. I think that using a later version should resolve your
problem.

For an easy sample matrix implementation look at:

http://www.parashift.com/c++-faq-lit...html#faq-16.17

Can this be done and if so, how? Thanks for any help,

Ben :-)

HTH
Chris

Thank you for your suggestion, Chris, but if I can figure out the friend
issue, the template strategy has several advantages. As I mentioned in my
first post, dimensions can be checked at compile-time, not at run-time.
This is true for all operations - addition, subtraction, inversion, etc. -
not just multiplication. Compile-time dimension checking is a boon for me
because the library will be used in an airborne real-time system where
run-time failure is not an option. The template solution I proposed also
makes it easier for the compiler to optimize out bounds-checks that are
gaurantee to pass (like your A(3,5) example) and to perform the named
return value optimizations. Templates also facilitate a custom new() and
delete() strategy that improves the speed of the filter the library is
used for by a factor of 3 by avoiding the creation and deletion of
temporary objects. I know that one can avoid unnecessary creation and
deletion using expression templates, but it seems to me that expression
templates require O(n^4) operations to evaluate A * B * C instead of the
O(n^3) operations that are necessary.

What I'm most curious about isn't whether the template solution is the
best one for my needs, it's how can I use the template solution and keep
the data private at the same time. The template solution may not be the
best one, but I'm still curious if and how this kind of templated friendship
works.

Ben :-)

[red floyd]

Ben Ingram wrote:

Hi all,

I am writing a template matrix class in which the template parameters are
the number of rows and number of columns. There are a number of reasons
why this is an appropriate tradeoff for my particular application. One of the
advantages is that the _compiler_ can force inner matrix dimensions used in
multiplication to agree. A _complie-time_ error will be triggered if you
write A * B and the number of coluns in A does not equal the number of
rows in B. Here's simplified code that illustrates the concept:

template<int nRows, int nCols> class Matrix {
public:
double data[nRows][nCols];

// operator*: return this * A
template<int nNewCols> Matrix<nRows,nNewCols>
operator*(const Matrix<nCols,nNewCols> &A) const {
Matrix<nRows,nNewCols> ret;
for(int iRow = 0; iRow < nRows; iRow++) {
for(int iCol = 0; iCol < nNewCols; iCol++) {
double innerProd = 0.0;
for(int iInner = 0; iInner < nCols; iInner++) {
innerProd += data[iRow][iInner] * A.data[iInner][iCol];
}
ret.data[iRow][iCol] = innerProd;
}
}
return ret;
}
};

int main(int argc, char **argv) {
Matrix<4,3> A;
Matrix<3,1> x;
Matrix<2,1> y;

A*x; // compiler creates
// Matrix<4,1> Matrix<4,3>::operator*(const Matrix<3,1>&)
//A*y; // compile-time error thrown since compiler can't create
// Matrix<4,1> Matrix<4,3>::operator*(const Matrix<2,1>&)
// since inner dimensions don't agree

return 0;
}

Note that while this code works under all g++ versions I tested, from
g++-2.7 through g++-3.2, it causes an internal compiler error in MSVC++
6.

The problem with this code is that the data field must be public. Otherwise,
Matrix<4,3>::operator* will not have access to Matrix<4,1>'s or
Matrix<3,1>'s data. Thus, I would like to make Matrix<4,1> and
Matrix<3,1> friends of Matrix<4,3>. Is this possible? If so, what's the
syntax? I tried the code below and some other variations, but they were
all incorrect. The code below generated "partial specialization
`Matrix<nRows, nNewCols>' declared `friend'" and "partial specialization
`Matrix<nCols, nNewCols>' declared `friend'" errors under g++-3.0.4 and
did not give Matrix<4,3> access to the private data of Matrix<4,1> or
of Matrix<3,1>.

template<int nRows, int nCols> class Matrix {
private:
double data[nRows][nCols];

public:
// operator*: return this * A
template<int nNewCols> Matrix<nRows,nNewCols>
operator*(const Matrix<nCols,nNewCols> &A) const {
Matrix<nRows,nNewCols> ret;
for(int iRow = 0; iRow < nRows; iRow++) {
for(int iCol = 0; iCol < nNewCols; iCol++) {
double innerProd = 0.0;
for(int iInner = 0; iInner < nCols; iInner++) {
innerProd += data[iRow][iInner] * A.data[iInner][iCol];
}
ret.data[iRow][iCol] = innerProd;
}
}
return ret;
}

template<int nNewCols> friend class Matrix<nRows,nNewCols>;
template<int nNewCols> friend class Matrix<nCols,nNewCols>;
};

int main(int argc, char **argv) {
Matrix<4,3> A;
Matrix<3,1> x;
Matrix<2,1> y;

A*x; // compiler creates
// Matrix<4,1> Matrix<4,3>::operator*(const Matrix<3,1>&)
//A*y; // compile-time error thrown since compiler can't create
// Matrix<4,1> Matrix<4,3>::operator*(const Matrix<2,1>&)
// since inner dimensions don't agree

return 0;
}

Can this be done and if so, how? Thanks for any help,

Ben :-)

I'd make the multiply operator a friend, rather than an internal function:

template<int Rows, int Cols>
class Matrix {
friend template <int R, int C, int N>
Matrix<R,N> operator*(const Matrix<R,C>&, const Matrix<C, N>&);

};

template <int R, int C, int N>
Matrix<R,N> operator*(const Matrix<R,C>&, const Matrix<C, N>&);

[Chris Theis]

"Ben Ingram" <bt************@catgufu.mit.edu> wrote in message
news:44******************************@news.sonicne ws.com...
[SNIP]

What I'm most curious about isn't whether the template solution is the
best one for my needs, it's how can I use the template solution and keep
the data private at the same time. The template solution may not be the
best one, but I'm still curious if and how this kind of templated friendship works.

Ben :-)

Well, to keep the data private without having to use friend declarations,
just use the operator() for access. The advantages of this are also covered
in the FAQ. However, to compile your class with VC++ you might have to
resort to version 7.X as there are many unresolved issures regarding partial
specialization with VC++ 6.0

Cheers
Chris

[Ben Ingram]

On Thu, 29 Jan 2004 08:59:42 -0800, red floyd wrote:

Ben Ingram wrote:

Hi all,

I am writing a template matrix class in which the template parameters
are the number of rows and number of columns. There are a number of
reasons why this is an appropriate tradeoff for my particular
application. One of the advantages is that the _compiler_ can force
inner matrix dimensions used in multiplication to agree. A
_complie-time_ error will be triggered if you write A * B and the
number of coluns in A does not equal the number of rows in B. Here's
simplified code that illustrates the concept:

template<int nRows, int nCols> class Matrix { public:
double data[nRows][nCols];

// operator*: return this * A
template<int nNewCols> Matrix<nRows,nNewCols>
operator*(const Matrix<nCols,nNewCols> &A) const {
Matrix<nRows,nNewCols> ret;
for(int iRow = 0; iRow < nRows; iRow++) {
for(int iCol = 0; iCol < nNewCols; iCol++) {
double innerProd = 0.0;
for(int iInner = 0; iInner < nCols; iInner++) {
innerProd += data[iRow][iInner] * A.data[iInner][iCol];
}
ret.data[iRow][iCol] = innerProd;
}
}
return ret;
}
};

int main(int argc, char **argv) {
Matrix<4,3> A;
Matrix<3,1> x;
Matrix<2,1> y;

A*x; // compiler creates
// Matrix<4,1> Matrix<4,3>::operator*(const Matrix<3,1>&)
//A*y; // compile-time error thrown since compiler can't create
// Matrix<4,1> Matrix<4,3>::operator*(const Matrix<2,1>&) //
since inner dimensions don't agree

return 0;
}
}
Note that while this code works under all g++ versions I tested, from
g++-2.7 through g++-3.2, it causes an internal compiler error in MSVC++
6.

The problem with this code is that the data field must be public.
Otherwise, Matrix<4,3>::operator* will not have access to Matrix<4,1>'s
or Matrix<3,1>'s data. Thus, I would like to make Matrix<4,1> and
Matrix<3,1> friends of Matrix<4,3>. Is this possible? If so, what's the
syntax? I tried the code below and some other variations, but they were
all incorrect. The code below generated "partial specialization
`Matrix<nRows, nNewCols>' declared `friend'" and "partial
specialization `Matrix<nCols, nNewCols>' declared `friend'" errors
under g++-3.0.4 and did not give Matrix<4,3> access to the private data
of Matrix<4,1> or of Matrix<3,1>.

template<int nRows, int nCols> class Matrix { private:
double data[nRows][nCols];

public:
// operator*: return this * A
template<int nNewCols> Matrix<nRows,nNewCols>
operator*(const Matrix<nCols,nNewCols> &A) const {
Matrix<nRows,nNewCols> ret;
for(int iRow = 0; iRow < nRows; iRow++) {
for(int iCol = 0; iCol < nNewCols; iCol++) {
double innerProd = 0.0;
for(int iInner = 0; iInner < nCols; iInner++) {
innerProd += data[iRow][iInner] * A.data[iInner][iCol];
}
ret.data[iRow][iCol] = innerProd;
}
}
return ret;
}
}
template<int nNewCols> friend class Matrix<nRows,nNewCols>;
template<int nNewCols> friend class Matrix<nCols,nNewCols>;
};

int main(int argc, char **argv) {
Matrix<4,3> A;
Matrix<3,1> x;
Matrix<2,1> y;

A*x; // compiler creates
// Matrix<4,1> Matrix<4,3>::operator*(const Matrix<3,1>&)
//A*y; // compile-time error thrown since compiler can't create
// Matrix<4,1> Matrix<4,3>::operator*(const Matrix<2,1>&) //
since inner dimensions don't agree

return 0;
}
}
Can this be done and if so, how? Thanks for any help,

Ben :-)

I'd make the multiply operator a friend, rather than an internal
function:

template<int Rows, int Cols>
class Matrix {
friend template <int R, int C, int N>
Matrix<R,N> operator*(const Matrix<R,C>&, const Matrix<C, N>&);

};

template <int R, int C, int N>
Matrix<R,N> operator*(const Matrix<R,C>&, const Matrix<C, N>&);

Perfect!!! Almost. friend should come after, not before template<...>.
Thanks so much for the help!!! Below is the modified code:

template<int nRows, int nCols> class Matrix {
private:
double data[nRows][nCols];

public:
// operator*: return this * A
template<int multRows, int multInner, int multCols>
friend Matrix<multRows,multCols> operator*(
const Matrix<multRows,multInner> &,
const Matrix<multInner,multCols> &);
};

template <int nRows, int nInner, int nCols>
Matrix<nRows,nCols> operator*(const Matrix<nRows,nInner> &A,
const Matrix<nInner,nCols> &B) {
Matrix<nRows,nCols> ret;
for(int iRow = 0; iRow < nRows; iRow++) {
for(int iCol = 0; iCol < nCols; iCol++) {
double innerProd = 0.0;
for(int iInner = 0; iInner < nInner; iInner++) {
innerProd += A.data[iRow][iInner] * B.data[iInner][iCol];
}
ret.data[iRow][iCol] = innerProd;
}
}
return ret;
}

int main(int argc, char **argv) {
Matrix<4,3> A;
Matrix<3,1> x;
Matrix<2,1> y;

A*x; // compiler creates
// Matrix<4,1> Matrix<4,3>::operator*(const Matrix<3,1>&)
//A*y; // compile-time error thrown since compiler can't create
// Matrix<4,1> Matrix<4,3>::operator*(const Matrix<2,1>&)
// since inner dimensions don't agree

return 0;
}

Ben :-)

[red floyd]

Ben Ingram wrote:

[Redacted for space]

Perfect!!! Almost. friend should come after, not before template<...>.
Thanks so much for the help!!! Below is the modified code:

[Redacted for space]

Glad I could help, Ben! To be honest, I've never tried playing with
friend templates, so I'm amazed that that I got the syntax as close as I
did.

red floyd