Hi there,

I am working on math / matrix library that will be open source (non-computer scientist) engineer friendly and capable of managing really big computations on shared memory machines.

I have created a few template classes describing matrices. The most basic matrix object looks like this:

template <class DATA_TYPE,class INT_TYPE=std::size_t>
class Matrix
{
//--------------------------------------------------------------------------
private:
INT_TYPE Rows, Cols;
DATA_TYPE *Data;
//--------------------------------------------------------------------------
//--------------------------------------------------------------------------
public:

Matrix(); //basic constructor
Matrix(const Matrix<DATA_TYPE,INT_TYPE> &M); //copy cnstructor
~Matrix(); //destructor
//--------------------------------------------------------------------------
inline void Alloc(INT_TYPE _Rows,INT_TYPE _Cols);
inline void Free();
inline void Resize(INT_TYPE _Rows,INT_TYPE _Cols);
inline void ResizeCopyNoFill(INT_TYPE _Rows,INT_TYPE _Cols);
inline void ResizeNoCopyNoFill(INT_TYPE _Rows,INT_TYPE _Cols);
inline void Fill(DATA_TYPE Value=DATA_TYPE());
inline void Copy(const Matrix<DATA_TYPE,INT_TYPE> &M);
//--------------------------------------------------------------------------
inline bool IsEmpty() const;
inline bool IsValid() const;
inline bool IsScalar() const;
inline bool IsVector() const;
inline bool IsRowVector() const;
inline bool IsColumnVector() const;
//--------------------------------------------------------------------------
inline Matrix<DATA_TYPE,INT_TYPE>& operator =(const Matrix<DATA_TYPE,INT_TYPE>& M); //A=B; Direct assignment operator
inline Matrix<DATA_TYPE,INT_TYPE>& operator +=(const Matrix<DATA_TYPE,INT_TYPE>& M); //A+=B; Addition operator
inline Matrix<DATA_TYPE,INT_TYPE>& operator -=(const Matrix<DATA_TYPE,INT_TYPE>& M); //A-=B; Substraction operator
inline Matrix<DATA_TYPE,INT_TYPE>& operator *=(const Matrix<DATA_TYPE,INT_TYPE>& M); //A*=B; Matrix multiplication operator
inline Matrix<DATA_TYPE,INT_TYPE>& operator /=(const Matrix<DATA_TYPE,INT_TYPE>& M); //A/=B; <=> A*=Inverse(B);
inline Matrix<DATA_TYPE,INT_TYPE>& operator ^=(const Matrix<DATA_TYPE,INT_TYPE>& M); //A^=B; Matrix elements to the power A^B
//--------------------------------------------------------------------------
inline Matrix<DATA_TYPE,INT_TYPE>& operator =(DATA_TYPE s); //M=s; Fills the matrix with s entries
inline Matrix<DATA_TYPE,INT_TYPE>& operator +=(DATA_TYPE s); //M+=s; Adds s to every element of M
inline Matrix<DATA_TYPE,INT_TYPE>& operator -=(DATA_TYPE s); //M-=s; Substracts s from every M element
inline Matrix<DATA_TYPE,INT_TYPE>& operator *=(DATA_TYPE s); //M*=s; Multiplies matrix elements by s
inline Matrix<DATA_TYPE,INT_TYPE>& operator /=(DATA_TYPE s); //M/=s; Divides matrix elements by s
inline Matrix<DATA_TYPE,INT_TYPE>& operator ^=(DATA_TYPE s); //A^=s; Matrix elements to the power A^s
//--------------------------------------------------------------------------
inline DATA_TYPE operator ()(INT_TYPE _Row,INT_TYPE _Col) const; // a=M(i,j);
inline DATA_TYPE& operator ()(INT_TYPE _Row,INT_TYPE _Col); // M(i,j)=a;
inline DATA_TYPE operator ()(INT_TYPE ElementIndex) const; // a=M(i);
inline DATA_TYPE& operator ()(INT_TYPE ElementIndex); // M(i)=a;
//--------------------------------------------------------------------------
inline DATA_TYPE* GetDataPtr() const {return Data;}
inline void SetDataPtr(DATA_TYPE* DataPtr) { Data=DataPtr;} //use with extreme caution
//--------------------------------------------------------------------------
inline DATA_TYPE GetElement(INT_TYPE _Row,INT_TYPE _Col) const; //returns element, but index checking is done
inline void SetElement(INT_TYPE _Row,INT_TYPE _Col,DATA_TYPE Value); //sets element, but index checking is done
inline DATA_TYPE GetElement(INT_TYPE ElementIndex) const; //returns element, but index checking is done
inline void SetElement(INT_TYPE ElementIndex,DATA_TYPE Value); //sets element, but no index checking is done
//--------------------------------------------------------------------------
inline INT_TYPE NumOfElements()const {return Rows*Cols;}
inline INT_TYPE NumOfRows() const {return Rows;}
inline INT_TYPE NumOfCols() const {return Cols;}
inline INT_TYPE SetNumOfRows(INT_TYPE _Rows) const { Rows=_Rows;}
inline INT_TYPE SetNumOfCols(INT_TYPE _Cols) const { Cols=_Cols;}
//--------------------------------------------------------------------------
inline DATA_TYPE Sum() const;
inline DATA_TYPE Product() const;
inline DATA_TYPE Max() const;
inline DATA_TYPE Min() const;
inline std::tuple<DATA_TYPE,DATA_TYPE> MinMax() const;
//--------------------------------------------------------------------------
inline void Transpose();
inline void Inverse();
void Print();
//--------------------------------------------------------------------------
};


Of course it needs a lot of work, as I am at the very beginning. My problem is:

How to define and implement operators like multiplication, addition etc to avoid copying the matrix. For now the multiplication
operator looks like that:

template <class DATA_TYPE,class INT_TYPE>
inline Matrix<DATA_TYPE,INT_TYPE> operator *(const Matrix<DATA_TYPE,INT_TYPE>& A,const Matrix<DATA_TYPE,INT_TYPE>& B)
{
Matrix<DATA_TYPE,INT_TYPE> C;
Multiply<DATA_TYPE,INT_TYPE>(A,B,C);
return C;
}



How to avoid copyying the whole C matrix when such expression is used:

M=A*B;


If for example I would modify direct assignment operator, to copy only *Data, Rows and Cols, that would do, but then A=B would remind more the reference than copy, as it would refer to the same place in memory. What do You think about allowing to use pointers only, like

Matrix<float,int> *M=new Matrix();


Than it would work, as all operations would be done on pointers to matrices, but I somehow have feeling that I didn't notice some kind of disadvantage of such solution (except that using a variable again would require using delete M; first).

So, please help me with even absolutely general ideas about solving this problem, so it would look elegant and be fast (in the sense of avoiding of copying matrices).

Regards,
Misery

Look into move constructors/assignment operator if you're using a C++11-compiler. They are, to describe it very briefly in your context, a different kind of constructors/assignment operators that operate on temporary values that will be immediately destroyed (such as the result of A*B before being assigned to C), which can be implemented very efficiently by just moving the resources from A*B to C (assign the data pointer from A*B to C) instead of copying it (allocate memory for C, copy the actual data from A*B to C, destroy A*B).

Another option could be copy-on-write. Basically, share the internal matrix-data in some reference counted way and split the data once a write operation is performed on the matrix in which case it makes an internal copy of itself. As long as you're just reading a matrix, its internal data can be shared with other matrices. In this case, "other matrices" includes the temporary values that you return.

I would not advise you to combine allocating matrices dynamically though. You now also have to manage the lifetime of the matrix objects and to ensure that they aren't lost in some assignment. Pointers and assignment-syntax of higher level objects doesn't really play well together if you're trying to simulate value-semantics.

Look into move constructors/assignment operator if you're using a C++11-compiler.

That is exactly what I have been looking for!!!! Thank You sooo much!!!! :]
The only remaining thing now is: should I change anything in my operator * (M1, M2) ??? In fact it just seems to work fine :]

And one moe thing. I have followed this simple tutorial. However one of the comments there says it is wrong. What version should I use?

You need to turn C into an R-value reference before returning it so that the return value is move-constructed instead of copy-constructed. Just change the return statement to:
[source]
return std::move( C );
[/source]
I believe the function is defined in the <utility> header. Other than this, your function should be as it is.


template <class DATA_TYPE,class INT_TYPE>
inline Matrix<DATA_TYPE,INT_TYPE> operator *(const Matrix<DATA_TYPE,INT_TYPE>& A,const Matrix<DATA_TYPE,INT_TYPE>& B)
{
Matrix<DATA_TYPE,INT_TYPE> C;
Multiply<DATA_TYPE,INT_TYPE>(A,B,C);
return C;
}


Using Multiply() like that might make it more difficult for compiler to perform return value optimization. Why not just return Multiply(A, B);

For operators which are commutative, you can easily use the resource from rvalue reference since it isn't const and it will get thrown away anyway (most likely a temporary is passed in):

template <class DATA_TYPE,class INT_TYPE>
inline Matrix<DATA_TYPE,INT_TYPE> operator +(Matrix<DATA_TYPE,INT_TYPE>&& A, const Matrix<DATA_TYPE,INT_TYPE>& B)
{
A += B;
return A;
}
template <class DATA_TYPE,class INT_TYPE>
inline Matrix<DATA_TYPE,INT_TYPE> operator +(const Matrix<DATA_TYPE,INT_TYPE>& A, Matrix<DATA_TYPE,INT_TYPE>&& B)
{
B += A;
return B;
}


You should look into Expression Templates, they've been invented exactly for the linear algebra operations you need:
http://www.altdevblogaday.com/2012/01/23/abusing-c-with-expression-templates/
http://en.wikibooks.org/wiki/More_C%2B%2B_Idioms/Expression-template
http://en.wikipedia.org/wiki/Expression_templates
http://www10.informatik.uni-erlangen.de/~pflaum/pflaum/ProSeminar/exprtmpl.html

Expression templates look like a quite complicated alternative to move constructors/operators. Is there any really good reason to use those instead of solution using move operator and/or move constructor?

After reading those articles I've decided to use such an techniqe in my code. I am very grateful for expressions templates, as well as for r-value operations sugesstions.

But in the topic of expression templates one more quesion arises: what if I have a general matrix class. It can be vector or square, or any other rectangular matrix. So how do I make an expression template that will "wisely" expand such cases? Is it worth even trying? Such cases can be O(N^2), what then?