Commit 9fbe1a7e authored by Alvarez, Gonzalo's avatar Alvarez, Gonzalo
Browse files

Matrix conversions from closures

parent f3f5d996

drivers/closuresTest.cpp

+5 −0
Original line number Diff line number Diff line
@@ -52,7 +52,12 @@ void testMatrix() 
	m3 = m1 + 1.3*m2;

	std::cout<<m3;

	m3 = 1.3*m2 + m1;

	std::cout<<"--------------\n";

	std::cout<<m3;

	std::cout<<m2;

	m3 = 0.5*(m3-m2);

	std::cout<<m3;

	std::cout<<"--------------\n";



	// minus tests omitted here

src/Matrix.h

+134 −42
Original line number Diff line number Diff line
@@ -77,6 +77,32 @@ public: 
				data_[i+j*nrow_] = m(i,j);

	}



	// ctor closures

	Matrix(const std::ClosureOperator<Matrix,Matrix,std::ClosureOperations::OP_MULT>& c)

	{

		operator=(*this,c.v1_*c.v2_);

	}



	template<typename T1>

	Matrix(const std::ClosureOperator<

	       std::ClosureOperator<T1,Matrix,std::ClosureOperations::OP_MULT>,

	       Matrix,

	       std::ClosureOperations::OP_MULT>& c,

	       typename EnableIf<Loki::TypeTraits<T1>::isArith,int>::Type = 0)

	{

		*this = c.v1_.v1_*c.v1_.v2_*c.v2_;

	}



	template<typename T1>

	Matrix(const std::ClosureOperator<T1,

	       std::ClosureOperator<Matrix,Matrix,std::ClosureOperations::OP_PLUS>,

	       std::ClosureOperations::OP_MULT>& c,

	       typename EnableIf<Loki::TypeTraits<T1>::isArith,int>::Type = 0)

	{

		*this = c.v1_*(c.v2_.v1_+c.v2_.v2_);

	}

	// end all ctors



	template<class Archive>

	void serialize(Archive& ar, const unsigned int)

	{
@@ -221,27 +247,115 @@ public: 
	// start closure memebers



	template<typename T1>

	void operator=(const std::ClosureOperator<T1,Matrix<T>,std::ClosureOperations::OP_MULT>& c)

	Matrix& operator+=(const std::ClosureOperator<T1,Matrix,std::ClosureOperations::OP_MULT>& c)

	{

		nrow_ = c.v2_.nrow_;

		ncol_ = c.v2_.ncol_;

		this->data_ += c.v1_*c.v2_.data_;

		return *this;

	}



	template<typename T1>

	Matrix& operator+=(const std::ClosureOperator<Matrix,T1,std::ClosureOperations::OP_MULT>& c)

	{

		nrow_ = c.v1_.nrow_;

		ncol_ = c.v1_.ncol_;

		this->data_ += c.v2_*c.v1_.data_;

		return *this;

	}



	template<typename T1>

	Matrix& operator=(const std::ClosureOperator<T1,Matrix,std::ClosureOperations::OP_MULT>& c)

	{

		nrow_ = c.v2_.nrow_;

		ncol_ = c.v2_.ncol_;

		this->data_ <= c.v1_*c.v2_.data_;

		return *this;

	}



	template<typename T1>

	Matrix& operator=(const std::ClosureOperator<Matrix,T1,std::ClosureOperations::OP_MULT>& c)

	{

		nrow_ = c.v1_.nrow_;

		ncol_ = c.v1_.ncol_;

		this->data_ <= c.v2_*c.v1_.data_;

		return *this;

	}



	template<typename T1>

	Matrix& operator=(const std::ClosureOperator<

	                  std::ClosureOperator<T1,Matrix,std::ClosureOperations::OP_MULT>,

	                  Matrix,

	                  std::ClosureOperations::OP_MULT>& c)

	{

		const Matrix<T>& a = c.v1_.v2_;

		const Matrix<T>& b = c.v2_;

		nrow_ = a.n_row();

		ncol_ = b.n_col();

		data_.resize(nrow_*ncol_);

		assert(a.n_col()==b.n_row());

		for (SizeType i=0;i<a.n_row();i++) {

			for (SizeType j=0;j<b.n_col();j++) {

				T sum = 0.0;

				for (SizeType k=0;k<a.n_col();k++) {

					sum += a(i,k) * b(k,j);

				}



	void operator=

	(const std::ClosureOperator<Matrix<T>,Matrix<T>,std::ClosureOperations::OP_PLUS>& c)

				this->operator()(i,j) = sum*c.v1_.v1_;

			}

		}



		return *this;

	}



	template<typename T1>

	Matrix& operator=(const std::ClosureOperator<T1,

	                  std::ClosureOperator<Matrix,Matrix,std::ClosureOperations::OP_PLUS>,

	                  std::ClosureOperations::OP_MULT>& c)

	{

		this->nrow_ = c.v2_.v1_.nrow_;

		this->ncol_ = c.v2_.v1_.ncol_;

		this->data_ <= c.v1_*(c.v2_.v1_.data_ + c.v2_.v2_.data_);

		return *this;

	}



	template<typename T1>

	Matrix& operator=(const std::ClosureOperator<T1,

	                  std::ClosureOperator<Matrix,Matrix,std::ClosureOperations::OP_MINUS>,

	                  std::ClosureOperations::OP_MULT>& c)

	{

		this->nrow_ = c.v2_.v1_.nrow_;

		this->ncol_ = c.v2_.v1_.ncol_;

		this->data_ <= c.v1_*(c.v2_.v1_.data_ - c.v2_.v2_.data_);

		return *this;

	}



	Matrix& operator=

	(const std::ClosureOperator<Matrix,Matrix,std::ClosureOperations::OP_PLUS>& c)

	{

		nrow_ = c.v1_.nrow_;

		ncol_ = c.v1_.ncol_;

		assert(nrow_ == c.v2_.nrow_);

		assert(ncol_ == c.v2_.ncol_);

		this->data_ <= c.v1_.data_ + c.v2_.data_;

		return *this;

	}



	Matrix& operator=

	(const std::ClosureOperator<Matrix,Matrix,std::ClosureOperations::OP_MINUS>& c)

	{

		nrow_ = c.v1_.nrow_;

		ncol_ = c.v1_.ncol_;

		assert(nrow_ == c.v2_.nrow_);

		assert(ncol_ == c.v2_.ncol_);

		this->data_ <= c.v1_.data_ - c.v2_.data_;

		return *this;

	}



	template<typename T1>

	void operator=

	(const std::ClosureOperator<Matrix<T>,

	 std::ClosureOperator<T1,Matrix<T>,std::ClosureOperations::OP_MULT>,

	Matrix& operator=

	(const std::ClosureOperator<Matrix,

	 std::ClosureOperator<T1,Matrix,std::ClosureOperations::OP_MULT>,

	 std::ClosureOperations::OP_PLUS>& c)

	{

		nrow_ = c.v1_.nrow_;
@@ -249,13 +363,14 @@ public: 
		assert(nrow_ == c.v2_.v2_.nrow_);

		assert(ncol_ == c.v2_.v2_.ncol_);

		this->data_ <= c.v1_.data_ + c.v2_.v1_*c.v2_.v2_.data_;

		return *this;

	}



	template<typename T1>

	void operator=

	Matrix& operator=

	(const std::ClosureOperator<

	 std::ClosureOperator<T1,Matrix<T>,std::ClosureOperations::OP_MULT>,

	 Matrix<T>,

	 std::ClosureOperator<T1,Matrix,std::ClosureOperations::OP_MULT>,

	 Matrix,

	 std::ClosureOperations::OP_PLUS>& c)

	{

		nrow_ = c.v2_.nrow_;
@@ -263,10 +378,11 @@ public: 
		assert(nrow_ == c.v1_.v2_.nrow_);

		assert(ncol_ == c.v1_.v2_.ncol_);

		this->data_ <= c.v2_.data_ + c.v1_.v1_*c.v1_.v2_.data_;

		return *this;

	}



	void operator=

	(const std::ClosureOperator<Matrix<T>,Matrix<T>,std::ClosureOperations::OP_MULTSAME>& c)

	Matrix& operator=

	(const std::ClosureOperator<Matrix<T>,Matrix,std::ClosureOperations::OP_MULT>& c)

	{

		const Matrix<T>& a = c.v1_;

		const Matrix<T>& b = c.v2_;
@@ -283,6 +399,8 @@ public: 
				this->operator()(i,j) = sum;

			}

		}



		return *this;

	}



	// end closure members
@@ -514,17 +632,12 @@ bool isZero(const Matrix<T>& m) 
// closures start



template<typename T1,typename T2>

std::ClosureOperator<T1,Matrix<T2>,std::ClosureOperations::OP_MULT> operator*(const T1& val,

                                                                              const Matrix<T2>& a)

{

	return std::ClosureOperator<T1,Matrix<T2>,std::ClosureOperations::OP_MULT>(val,a);

}



template<typename T1,typename T2>

std::ClosureOperator<T1,Matrix<T2>,std::ClosureOperations::OP_MULT>

operator*(const Matrix<T2>& a, const T1& val)

typename EnableIf<(IsMatrixLike<T1>::True || IsMatrixLike<T2>::True)

&& !std::IsClosureLike<T1>::True && !std::IsClosureLike<T2>::True,

std::ClosureOperator<T1,T2,std::ClosureOperations::OP_MULT> >::Type operator*(const T1& a,

                                                                              const T2& b)

{

	return val*a;

	return std::ClosureOperator<T1,T2,std::ClosureOperations::OP_MULT>(a,b);

}



template<typename T>
@@ -541,20 +654,6 @@ operator-(const Matrix<T>& a,const Matrix<T>& b) 
	return std::ClosureOperator<Matrix<T>,Matrix<T>,std::ClosureOperations::OP_MINUS>(a,b);

}



template<typename T>

std::ClosureOperator<Matrix<T>,Matrix<T>,std::ClosureOperations::OP_MULTSAME>

operator*(const Matrix<T>& a,const Matrix<T>& b)

{

	return std::ClosureOperator<Matrix<T>,Matrix<T>,std::ClosureOperations::OP_MULTSAME>(a,b);

}



template<typename T>

std::ClosureOperator<Matrix<T>,std::vector<T>,std::ClosureOperations::OP_MULT>

operator*(const Matrix<T>& a, const std::vector<T>& v)

{

	return std::ClosureOperator<Matrix<T>,std::vector<T>,std::ClosureOperations::OP_MULT>(a,v);

}



template<typename T,typename A>

void operator<=(std::vector<T,A>& v, const std::ClosureOperator<Matrix<T>,

                std::vector<T,A>,
@@ -571,13 +670,6 @@ void operator<=(std::vector<T,A>& v, const std::ClosureOperator<Matrix<T>, 
	}

}



template<typename T>

std::ClosureOperator<std::vector<T>,Matrix<T>,std::ClosureOperations::OP_MULT>

operator*(const std::vector<T>& v, const Matrix<T>& a)

{

	return std::ClosureOperator<std::vector<T>,Matrix<T>,std::ClosureOperations::OP_MULT>(v,a);

}



template<typename T,typename A>

void operator<=(std::vector<T,A>& v, const std::ClosureOperator<std::vector<T,A>,

                Matrix<T>,

src/Vector.h

+37 −12
Original line number Diff line number Diff line
@@ -31,6 +31,7 @@ Please see full open source license included in file LICENSE. 
#include <stdexcept>

#include "Complex.h"

#include "AllocatorCpu.h"

#include "../loki/TypeTraits.h"



namespace PsimagLite {
@@ -93,7 +94,7 @@ vector<T2,A> operator*(const vector<vector<T1,A>,AA>& v1, 


struct ClosureOperations {



	enum {OP_PLUS,OP_MINUS,OP_MULT,OP_DIVIDE,OP_MULTSAME,OP_CONJ};

	enum {OP_PLUS,OP_MINUS,OP_MULT,OP_DIVIDE,OP_CONJ};

};



template<typename T1, typename T2,int type>
@@ -134,21 +135,17 @@ ClosureOperator<T1,vector<T2,A>,ClosureOperations::OP_MULT> operator*(const T1& 
}



// vector * scalar

template<typename T1,typename T2,typename A>

ClosureOperator<T1,

ClosureOperator<T2,std::vector<T2,A>,ClosureOperations::OP_MULT>,

ClosureOperations::OP_MULT> operator*(const T1& v1,

                                      const ClosureOperator<T2,

                                      std::vector<T2,A>,

                                      ClosureOperations::OP_MULT>& v2)

template<typename T1,typename T2>

typename PsimagLite::EnableIf<IsClosureLike<T1>::True || IsClosureLike<T2>::True,

ClosureOperator<T1,T2,ClosureOperations::OP_MULT> >::Type

operator*(const T1& v1,const T2& v2)

{

	return ClosureOperator<T1,

	        ClosureOperator<T2,std::vector<T2,A>,ClosureOperations::OP_MULT>,

	        ClosureOperations::OP_MULT>(v1,v2);

	return ClosureOperator<T1,T2,ClosureOperations::OP_MULT>(v1,v2);

}



template<typename T1,typename T2,typename A>

void operator<=(vector<T2,A>& v,

typename PsimagLite::EnableIf<Loki::TypeTraits<T1>::isArith,void>::Type

operator<=(vector<T2,A>& v,

                const ClosureOperator<T1,vector<T2,A>,ClosureOperations::OP_MULT>& c)

{

	v = c.v2_;
@@ -218,6 +215,16 @@ void operator<=(vector<T1,A1>& v, 
	for (SizeType i=0;i<v.size();i++) v[i] = c.v1_.v1_*c.v1_.v2_[i] + c.v2_[i];

}



template<typename T1,typename T2,typename A>

void operator<=(vector<T2,A>& v,

                const ClosureOperator<T1,

                ClosureOperator<vector<T2,A>,vector<T2,A>, ClosureOperations::OP_PLUS>,

                ClosureOperations::OP_MULT>& c)

{

	v.resize(c.v2_.v2_.size());

	for (SizeType i=0;i<v.size();i++) v[i] = c.v1_*(c.v2_.v1_[i] + c.v2_.v2_[i]);

}



template<typename T1,typename T2,typename A>

void operator<=(vector<T2,A>& v,

                const ClosureOperator<
@@ -260,6 +267,14 @@ operator-(const T1& v1,const T2& v2) 
	return ClosureOperator<T1,T2,ClosureOperations::OP_MINUS>(v1,v2);

}



template<typename T,typename A>

void operator<=(vector<T,A>& v,

                const ClosureOperator<vector<T,A>,vector<T,A>,ClosureOperations::OP_MINUS>& c)

{

	v.resize(c.v1_.size());

	for (SizeType i=0;i<v.size();i++) v[i] = c.v1_[i] - c.v2_[i];

}



template<typename T1,typename T2,typename A1, typename A2>

void operator<=(vector<T1,A1>& v,

                const ClosureOperator<vector<T1,A1>,
@@ -270,6 +285,16 @@ void operator<=(vector<T1,A1>& v, 
	for (SizeType i=0;i<v.size();i++) v[i] = c.v1_[i] - c.v2_.v1_*c.v2_.v2_[i];

}



template<typename T1,typename T2,typename A>

void operator<=(vector<T2,A>& v,

                const ClosureOperator<T1,

                ClosureOperator<vector<T2,A>,vector<T2,A>, ClosureOperations::OP_MINUS>,

                ClosureOperations::OP_MULT>& c)

{

	v.resize(c.v2_.v2_.size());

	for (SizeType i=0;i<v.size();i++) v[i] = c.v1_*(c.v2_.v1_[i] - c.v2_.v2_[i]);

}



template<typename T1, typename T2, typename A>

void operator<=(vector<T2,A>& v,

                const ClosureOperator<