Overloading the multiplication operator in c++

Question

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.

Answer 1

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

}