11 replies to this topic

[Misery]

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[i]^B[i]
//--------------------------------------------------------------------------
		 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[i]^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

Edited by Misery, 28 June 2012 - 04:33 AM.

[Brother Bob]

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.

[Misery]

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?

Edited by Misery, 28 June 2012 - 06:01 AM.

[Brother Bob]

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:

return std::move( C );

I believe the function is defined in the <utility> header. Other than this, your function should be as it is.

Edited by Brother Bob, 28 June 2012 - 05:56 AM.

[japro]

Also look at expression templates. Those are used to avoid creation of temporary objects to begin with and also help the compiler to optimize operations.

[Codarki]

		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;
}

Edited by Codarki, 29 June 2012 - 03:05 AM.

[Matt-D]

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

[Misery]

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?

Edited by Misery, 02 July 2012 - 02:07 AM.

[japro]

They don't to the exact same thing. Simply put they allow to "coalesce" loops. So when you write:

Vector a(1024), b(1024), c(1024);
a = b*3.14+c;

then a naive implementation gets compiled to something like:

Vector tmp(1024);
for(int i = 0;i<1024;++i)
   tmp[i] = b[i] * 3.14;
for(int i = 0;i<1024;++i)
   a = tmp[i] + c[i];

but expression templates can give you:

for(int i = 0;i<1024;++i)
   a[i] = b[i] * 3.14 + c[i];

[Misery]

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?

Edited by Misery, 02 July 2012 - 06:53 AM.

[Matt-D]

I don't think you need or even should do everything for general matrices -- in general, don't do things that don't make sense in linear algebra.
For instance, for multiplication, operator*, you only need to concern yourself with the conformable ones.
For addition, operator+ (and subtraction, operator-), it only makes sense to look at ones of the same size. So, only implement the code for the right size, as in here: http://www.drdobbs.c...lates/184401656
// Notice that there are two ways to check for the right size -- if you encode the size in your type (for fixed-size vectors/arrays -- for instance like it's done in std::array<T,N>) you only accept the right-sized (N elements) arguments and thus you're checking at compile-time; OTOH, if you store the size (for resizable vectors/arrays -- for instance, like it's done in std::vector<T>) you have to assert on entrance to the operator (or just crash if that's your preference). All of these are choices and have various trade-offs involved :-)

You may consider how it's done in the other libraries -- a prominent example is Eigen (very fast):
http://eigen.tuxfami...EigenTypes.html
http://eigen.tuxfami...Evaluation.html
http://eigen.tuxfami...ralProduct.html
http://eigen.tuxfami...ionXprTemplates
It's open-source, so you can look at the code, but, again -- this is a state-of-the-art, highly-optimized library, so it may not be a good idea if you're still learning (it also uses plenty of SIMD extensions and, thus, intrinsics).

Perhaps consider looking at some of these resources:
http://www.issi.uned...isCzarnecki.pdf -- Section 10.3.7.1 provides sample code with discussion (and introduction in preceding sections)
http://dlib.net/matr...ons_ex.cpp.html // expression templates (ET) are thoroughly discussed in the comments -- you can click on the header files to follow the source code
http://tvmet.sourceforge.net/ // "This Tiny Vector and Matrix template library uses Meta Templates (MT) and Expression Templates (ET) to evaluate results at compile time -- which makes it fast for low order (tiny) systems."
http://www.chem.utor...ode/vecmat3v1.h
http://arxiv.org/abs/1104.1729v1 // in particular, "DMatDMatMultExpr"; note, that they propose a variation upon ET, called "smart" ET, which may or may not be what you need -- tellingly, they don't include Eigen in their benchmarks, but Boost.uBLAS (which is much slower than Eigen) and Blitz++ (which hasn't been maintained for years)
http://www10.informa...echRep09-17.pdf // again, this is also from authors "selling" the "smart" ET approach, so expect some criticism of ET -- but may still be useful; in their case they make use of the fact that a matrix-vector product is again a vector and make use of/express this fact as follows:
class DMatDVecMultExpr : public DenseVector < DMatDVecMultExpr <MT,VT> >, private Expression
- similarly, matrix-matrix product is a matrix:
class DMatDMatMultExpr : public DenseMatrix < DMatDMatMultExpr <MT1,MT2> >, private Expression

In general, rely on what you know from linear algebra and don't try to make anything too weird/fancy -- and you should be fine! :-)

Edited by Matt-D, 02 July 2012 - 09:39 AM.

[Misery]

Thanks a lot Gyus! :]