4

I'm trying to implement a matrix class that can handle both real and complex matrices. I'm running into problems when I try to overload the multiplication operator. Specifically, when I try to multiply a double matrix by a complex matrix (in that order). The result should be complex, but the * operator is a member of a double matrix in that case, I don't see how to return a complex matrix (I've tried using a friend operator, this doesn't seem to work either). Here's the relevant code snippet:

template <class V> Matrix<T> operator *(const Matrix<V> &b)
{
    long i, j, k;  T temp;  Matrix<T> c(mRows, b.Cols());

    for (i = 0; i < mRows; i++)
        for (j = 0; j < c.Cols(); j++)
        {
            for (temp = 0, k = 0; k < mCols; k++)
                temp += this->Element(i, k)*b.Element(k, j);

            c(i, j) = temp;
        }

        return c;
}

So, if A is real, and B and C are complex, C = BA works fine, but C = AB doesn't. Specifically, the compiler error flags the "temp +=" line, saying there is no global operator found which takes type complex (or no acceptable conversion). I'd like avoid completely specializing the template for real and complex types, is there a way around this?

Thanks in advance.

DJM123
  • 91
  • 3
  • You want to declare `template Matrix operator *(const Matrix &a, const Matrix &b)` – πάντα ῥεῖ Apr 03 '16 at 14:28
  • See also [Operator Overloading](http://stackoverflow.com/questions/4421706/operator-overloading) please. – πάντα ῥεῖ Apr 03 '16 at 14:29
  • Thanks, but this would have to be declared as a friend operator, no? And this creates a new problem, namely that more than one operator matches the operands. – DJM123 Apr 03 '16 at 14:52
  • _"but this would have to be declared as a friend operator, no?"_ Yes, that's one way of solving this. – πάντα ῥεῖ Apr 03 '16 at 14:56
  • Why not implement a global operator for your "temp +=" problem? – Andreas Apr 03 '16 at 15:00
  • How do I implement a global operator? – DJM123 Apr 03 '16 at 15:06
  • @DJM123 In addition, you should overload `*=` instead, and then implement `*` in terms of `*=`. – PaulMcKenzie Apr 03 '16 at 15:20
  • @PaulMcKenzie: yes, in general. But in linear algebra related applications I doubt that would be appropriate. Consider for example `real matrix *= complex matrix` -- that does not work. Further, you can have a lot of changes in the geometry (a vector becomes a matrix, etc.). – davidhigh Apr 03 '16 at 15:31
  • @PaulMcKenzie: agreed in general, but for matrix multiplication, you have to potentially resize the LHS, no? Unlike matrix addition. (Since the new dimensions after multiplication have to conform to the RHS factor, the columns may differ.) So, you have to create local copies anyway, no? – DJM123 Apr 04 '16 at 01:41

1 Answers1

3

As mentioned in the comments, in this case it's better to use a non-member overload of operator* (which is be-friended with the matrix class).

But that is not the reason why it does not work. The problem here is the declaration of

T temp; Matrix<T> c;

in your in-class operator*. You should replace each T with std::common_type_t<T,U>:

using C = std::common_type_t<T,U>;
C temp; Matrix<C> c;

Otherwise the element-type of the returned matrix will always be the same as that of the calling matrix. That is, when you call C = A*B which is resolved as C = A.operator*(B), the element type of C will be the same as that of A.

Then if A is a Matrix<double> and you multiply it with a Matrix<complex<double> >, you would get a Matrix<double> as return -- which is obviously wrong.

davidhigh
  • 12,239
  • 1
  • 34
  • 64
  • Fantastic, love this solution, thanks! It in fact works with the original, binary formulation of overloading *; is there a reason to eschew this approach, i.e. are friend overloadings preferable? – DJM123 Apr 03 '16 at 15:51
  • @DJM123: you're welcome. wrt the question, see [here](http://stackoverflow.com/a/4421729/2412846). – davidhigh Apr 03 '16 at 15:55
  • @DJM123: ... but my link does not give a motivation for using the non-member overload, therefore better see [here](http://stackoverflow.com/a/4622467/2412846). – davidhigh Apr 03 '16 at 15:59