Template Class Specialization

[woessner@gmail.com]

Hi all,

I whipped up a quick class to represent a matrix for use with LAPACK.
It's a template class so it can support the 4 data types supported by
LAPACK (single/double x complex/real). I added a conversion operator
to automatically convert the object to a pointer of the appropriate
type. This makes using LAPACK in C++ a lot easier.

Unfortunately, it does not work well for the complex data types. The
reason is that LAPACK (or, more specifically, f2c) defines its own
complex types, complex and doublecomplex. I could just use LAPACK's
types, but they don't provide any operators for them, making them a
real pain. Instead, I would much prefer to use std::complex.

This is where the problem is. The conversion operator in my matrix
class returns a T* where T is the underlying data type. For complex
data types, it needs to return a pointer to an LAPACK type,
specifically a complex* or doublecomplex*. I thought template
specialization would help me out, but it's turned out to be a real
pain.

Here's the code in question:

template<typena me T>
class LAPACK_Matrix
{
public:
operator T*() { return &m_Data[0]; }
/* other operators */

private:
std::vector<Tm_ Data;
};

So I thought I would just provide a specialized version for the complex
data types:

template<>
class LAPACK_Matrix<s td::complex<dou ble
{
public:
operator doublecomplex*( ) { return &m_Data[0]; }
};

I thought this would be sufficient, but (as I'm sure you all know), the
specialized class does not "inherit" anything from the unspecialized
class. So in order to make this work, I would have to re-implement all
my other operators in the specialized class. That doesn't seem like a
good solution to me.

So I've been coming up with a reasonable way to avoid having to rewrite
a bunch of code and I'm just stumped. I thought maybe inheritance
could solve my problem, but I don't see how. I thought about trying to
combine templates and in heritance, but I'm not sure how that would
work, either.

If anyone has any suggestions, I'd love to hear them.

Thanks in Advance,
Bill

[Victor Bazarov]

Re: Template Class Specialization

woessner@gmail. com wrote:

[...]
Here's the code in question:
>
template<typena me T>
class LAPACK_Matrix
{
public:
operator T*() { return &m_Data[0]; }
/* other operators */
>
private:
std::vector<Tm_ Data;
};
>
So I thought I would just provide a specialized version for the
complex data types:
>
template<>
class LAPACK_Matrix<s td::complex<dou ble
{
public:
operator doublecomplex*( ) { return &m_Data[0]; }
};
>
I thought this would be sufficient, but (as I'm sure you all know),
the specialized class does not "inherit" anything from the
unspecialized class. So in order to make this work, I would have to
re-implement all my other operators in the specialized class. That
doesn't seem like a good solution to me.

You could inherit...

So I've been coming up with a reasonable way to avoid having to
rewrite a bunch of code and I'm just stumped. I thought maybe
inheritance could solve my problem, but I don't see how. I thought
about trying to combine templates and in heritance, but I'm not sure
how that would work, either.
>
If anyone has any suggestions, I'd love to hear them.

Yes, inherit. Lessee...

template<class Tstruct MyT_Base {
// all the stuff
};

template<class Tstruct MyT : MyT_Base<T{
operator T*() { return 0; } // special stuff
};

template<struct MyT<double: MyT_Base<double {
operator float*() { return 42.; } // special stuff
};

V
--
Please remove capital 'A's when replying by e-mail
I do not respond to top-posted replies, please don't ask

[woessner@gmail.com]

Re: Template Class Specialization

Victor Bazarov wrote:

template<class Tstruct MyT_Base {
// all the stuff
};
template<class Tstruct MyT : MyT_Base<T{
operator T*() { return 0; } // special stuff
};
template<struct MyT<double: MyT_Base<double {
operator float*() { return 42.; } // special stuff
};

Victor, this is a great idea. And it works really well. I just have
one lingering question.

The base class has a constructor that takes 2 ints. As far as I can
tell, I have to implement similiar constructors in the derived classes
that do nothing but call the base class's constructor. This isn't a
really big deal (and it's a VAST improvement over the way it was
before), but I'm wondering if there's any way around it. I tried a
using clause. That didn't get me very far. :-)

Thanks,
Bill

[Noah Roberts]

Re: Template Class Specialization


woessner@gmail. com wrote:

Victor Bazarov wrote:

template<class Tstruct MyT_Base {
// all the stuff
};
template<class Tstruct MyT : MyT_Base<T{
operator T*() { return 0; } // special stuff
};
template<struct MyT<double: MyT_Base<double {
operator float*() { return 42.; } // special stuff
};

>
Victor, this is a great idea. And it works really well. I just have
one lingering question.
>
The base class has a constructor that takes 2 ints. As far as I can
tell, I have to implement similiar constructors in the derived classes
that do nothing but call the base class's constructor. This isn't a
really big deal (and it's a VAST improvement over the way it was
before), but I'm wondering if there's any way around it. I tried a
using clause. That didn't get me very far. :-)

You can inherit in the other direction. If you need access to members
of your object in the cast operator you can use the curriously
reoccuring template pattern...look it up in google. Also look up
policy based programming - you would basically be creating a "casting
policy" and inheriting from it to gain that behavior.

[Victor Bazarov]

Re: Template Class Specialization

woessner@gmail. com wrote:

Victor Bazarov wrote:

>template<cla ss Tstruct MyT_Base {
> // all the stuff
>};
>template<cla ss Tstruct MyT : MyT_Base<T{
> operator T*() { return 0; } // special stuff
>};
>template<struc t MyT<double: MyT_Base<double {
> operator float*() { return 42.; } // special stuff
>};

>
Victor, this is a great idea. And it works really well. I just have
one lingering question.
>
The base class has a constructor that takes 2 ints. As far as I can
tell, I have to implement similiar constructors in the derived classes
that do nothing but call the base class's constructor. This isn't a
really big deal (and it's a VAST improvement over the way it was
before), but I'm wondering if there's any way around it. I tried a
using clause. That didn't get me very far. :-)

Yes, you'll have to reimplement. Constructors are not inherited.

V
--
Please remove capital 'A's when replying by e-mail
I do not respond to top-posted replies, please don't ask

[Howard Gardner]

Re: Template Class Specialization

woessner@gmail. com wrote:

The base class has a constructor that takes 2 ints. As far as I can
tell, I have to implement similiar constructors in the derived classes
that do nothing but call the base class's constructor. This isn't a
really big deal (and it's a VAST improvement over the way it was
before), but I'm wondering if there's any way around it. I tried a
using clause. That didn't get me very far. :-)

If you have to write the base class constructor(s), then you have to
write the derived class constructor(s).

It might be possible to avoid having to write the base class
constructor(s), in which case you won't have to write the derived class
constructor(s).

Can you show us the code? If so, it would be interesting to see the
member variables, the constructor(s), the class copy assignment
operator, and the destructor.

[woessner@gmail.com]

Re: Template Class Specialization

Howard Gardner wrote:

Can you show us the code? If so, it would be interesting to see the
member variables, the constructor(s), the class copy assignment
operator, and the destructor.

Sure. Just... be kind. I'm a mathematician, not a programmer. :-p
In all honesty, it's already great the way it is. This is supposed to
be a really light-weight class that just makes it easier to do stuff
with LAPACK.

Here's what I have thus far:

template<typena me T>
class MatrixBase
{
public:
T& operator()(int iRow, int iColumn) { return m_Data[m_iRows *
iColumn + iRow]; }

int GetRows() const { return m_iRows; }

int GetColumns() const { return m_iColumns; }

protected:
MatrixBase(int iRows, int iColumns): m_iRows(iRows),
m_iColumns(iCol umns), m_Data(iRows * iColumns) {}

std::vector<Tm_ Data;

private:
int m_iRows;
int m_iColumns;
};

template<typena me T>
class Matrix: public MatrixBase<T>
{
public:
Matrix(int iRows, int iColumns): MatrixBase<T>(i Rows, iColumns) {}

operator T*() { return &m_Data[0]; }
};

template<>
class Matrix<std::com plex<double: public
MatrixBase<std: :complex<double
{
public:
Matrix(int iRows, int iColumns): MatrixBase<std: :complex<double >

>(iRows, iColumns) {}

operator doublecomplex*( ) { return
reinterpret_cas t<doublecomplex *>(&m_Data[0]); }
};

[Howard Gardner]

Re: Template Class Specialization

There are some potentially troublesome things in your implementation.
Whether they are real problems depends (as always) on what you mean to
do with it. If it is something that you intend to use once and then
forget about, they probably aren't. If it is something at the heart of
an important program (ie, your thesis work), then they may well be.

Before I get into that, though: have you looked at blitz++? It has an
excellent reputation.

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

Here's another implementation.

The main difference between this one and yours is that this introduces
the "non type template parameters" ROWS and COLUMNS. Using this
implementation, C++ will treat a 2x2 matrix and a 3x3 matrix as
different types. Using your implementation, it treated them as the same
type. The benefit is a cleaner implementation of the functionality that
is here, and I think it will also improve the implementation of the
functionality (the matrix operators) that are not here. The cost is that
this version might use more memory than your version.

A second important difference is that I've bounds-checked the operator()
for accessing the elements. If you specify an index out of bounds then
it will throw an exception. There is a performance penalty for this.

I've commented the reinterpret_cas t out because I don't have the library
that you're using. That function is indispensable to your purpose, but
it contains two assumptions: first that doublecomplex is bit compatible
with std::complex<do uble>, and second that std::vector is storing its
data as an array. Your program works, so the assumptions are true for
your current combination of compiler/libraries. It might not be true
with a different combination.

#include<vector >
#include<comple x>
#include<algori thm>
#include<stdexc ept>

#include<ostrea m// for testing

// ----------------------------------
template<typena me T, unsigned ROWS, unsigned COLUMNS>
class MatrixBase
{
public:
MatrixBase();
T& operator()(int iRow, int iColumn);
static const unsigned rows = ROWS;
static const unsigned columns = COLUMNS;
protected:
std::vector<Tm_ Data;
~MatrixBase();
};

// ---
template<typena me T, unsigned ROWS, unsigned COLUMNS>
MatrixBase<T, ROWS, COLUMNS>::Matri xBase()
:
m_Data(rows * columns)
{
}

// ---
template<typena me T, unsigned ROWS, unsigned COLUMNS>
T& MatrixBase<T, ROWS, COLUMNS>::opera tor()(int iRow, int iColumn)
{
if(iRow >= rows) throw std::range_erro r("row");
if(iColumn >= columns) throw std::range_erro r("columns");
return m_Data[rows * iColumn + iRow];
}

// ---
template<typena me T, unsigned ROWS, unsigned COLUMNS>
MatrixBase<T, ROWS, COLUMNS>::~Matr ixBase()
{
}

// ---
template<typena me S, typename T, unsigned ROWS, unsigned COLUMNS>
S& operator<<(S & stream, MatrixBase<T, ROWS, COLUMNS& matrix)
{
for(unsigned row = 0; row < matrix.rows; ++row)
{
for(unsigned column = 0; column < matrix.columns; ++column)
{
stream << matrix( row, column ) << ' ';
}
stream << '\n';
};

return stream;
}

// ----------------------------------
template<typena me T, unsigned ROWS, unsigned COLUMNS>
class Matrix: public MatrixBase<T, ROWS, COLUMNS>
{
public:
operator T*() { return &MatrixBase< T, ROWS, COLUMNS>::m_Dat a[0]; }
};

// ----------------------------------
template< unsigned ROWS, unsigned COLUMNS >
class Matrix<std::com plex<double>, ROWS, COLUMNS>: public
MatrixBase<std: :complex<double >, ROWS, COLUMNS>
{
public:
// operator doublecomplex*( );
};

// ---
//Matrix<std::com plex<double::op erator doublecomplex*( )
//{
// return reinterpret_cas t<doublecomplex *>(&m_Data[0]);
//}


// Test code below this line //
template< typename T, unsigned ROWS, unsigned COLUMNS >
Matrix< T, ROWS, COLUMNS func( Matrix< T, ROWS, COLUMNS & a )
{
return a;
}

int main(void)
{
using std::cout;
using std::endl;

Matrix< int, 2, 2 a;
a( 0, 1 ) = 1;
a( 1, 0 ) = 2;
a( 1, 1 ) = 3;
cout << a << endl;

Matrix< int, 2, 2 b(a);
Matrix< int, 2, 2 c;
c = a;

Matrix< std::complex< double >, 20, 20 d;
d( 5, 5 ) = std::complex< double >( 1, 1 );
cout << d << endl;

Matrix< std::complex< double >, 20, 20 e(d);
cout << e << endl;
}