Reputation: 386
I've written a C++ interface to LAPACK, but I'm running into some memory issues that have made me reconsider some of operator overloading.
Right now, I have overloaded the operator* outside of the the class definition (but as a friend to the Matrix class) that takes two Matrix objects, allocates a third with the proper dimensions, uses D(GE/SY)MM to compute the product (storing into the internal storage of the newly allocated matrix) and then returns the pointer to that new matrix. I.E.
class Matrix {
...
friend Matrix* operator*(const Matrix&, const Matrix&);
...
}
Matrix* operator*(const Matrix& m1, const Matrix& m2) {
Matrix *prod = new Matrix(m1.rows_, m2.cols_);
if(m1.cols_!=m2.rows_) {
throw 3008;
} else {
double alpha = 1.0;
double beta = 0.0;
if(m1.symm_=='G' && m2.symm_=='G'){
dgemm_(&m1.trans_,&m2.trans_,&m1.rows_,&m2.cols_,&m1.cols_,&alpha,m1.data_,
&m1.rows_,m2.data_,&m1.cols_,&beta,prod->data_,&m2.cols_);
} else if(m1.symm_=='S'){
char SIDE = 'L';
char UPLO = 'L';
dsymm_(&SIDE,&UPLO,&m1.rows_,&m2.cols_,&alpha,m1.data_,&m1.rows_,m2.data_,
&m2.cols_,&beta,prod->data_,&m2.cols_);
} else if(m2.symm_=='S'){
char SIDE = 'R';
char UPLO = 'L';
dsymm_(&SIDE,&UPLO,&m2.rows_,&m1.cols_,&alpha,m2.data_,&m2.rows_,m1.data_,
&m1.cols_,&beta,prod->data_,&m1.cols_);
};
}
return prod;
};
Then I utilize
Matrix *A, *B, *C;
// def of A and B
C = (*A)*(*B);
And this works just fine. The problem I'm having is that I have to allocate a new matrix every time I do this. What I'd like to be able to do is allocate the C
matrix once and place the product of A
and B
into the internal storage of C
(C->data_
). From what I've been able to find on operator overloading, I can't find a nice way to do this. I'm aware I can use a member function to do this, (i.e. C->mult(A,B)
) but I'd like to avoid that if at all possible (I'm coding this for ease of development for non-CSE types). Any ideas would be greatly appreciated.
Upvotes: 1
Views: 3758
Reputation: 2172
class Matrix
{
struct Product
{
const Matrix* a;
const Matrix* b;
};
Matrix& operator = (const Product& p)
{
// if this matrix dims differ from required by product of p.a and p.b
// reallocate it first and set dims
// {
// rows = ....; cols = ....;
// delete [] data;
// data = new [rows*cols];
// }
// then calculate product
// data[0] = ...;
// ...
return *this;
}
Product operator * (const Matrix& op) const
{
Product p;
p.a = this;
p.b = &op;
return p;
}
int rows,cols;
double* data;
/// your Matrix stuff
// ...
};
void test()
{
Matrix a(4,2),b(2,4),c;
c = a * b; // (note a*b returns Product without calculating its result
// result is calculated inside = operator
c = a * b; // note that this time c is initialized to correct size so no
// additional reallocations will occur
}
Upvotes: 2