Add MatDiag::multiply(MatDense), MatDiag::multiply(Vect),...

Add MatDiag::multiply(MatDense), MatDiag::multiply(Vect), MatDiag::transpose(),conjugate(), ctor for non-square diag. matrix.

misc/test/src/C++/test_MatDiag.cpp.in

@@ -10,12 +10,37 @@ typedef @TEST_FPP@ FPP;
 typedef @TEST_FPP2@ FPP2;
 
 
+
+void display_all_elts(MatDense<FPP,Cpu> matdense)
+{
+	for(int i=0;i<matdense.getNbRow();i++)
+	{
+		for(int j=0;j<matdense.getNbCol();j++)
+			cout << matdense(i,j) << " ";
+		cout << endl;
+	}
+}
+
+void display_all_elts(Faust::MatSparse<FPP,Cpu> matsp)
+{
+	display_all_elts(MatDense<FPP,Cpu>(matsp));
+}
+
+void display_all_elts(MatDiag<FPP> matdiag)
+{
+	display_all_elts(Faust::MatSparse<FPP,Cpu>(matdiag));
+}
+
 int main()
 {
 	cout << "start of test_MatDiag" << endl;
 	FPP data[10];
 	for(int i=0; i < 10; i++)
+#if(@TEST_IS_COMPLEX@==1)
+		data[i] = FPP(i,i);
+#else
 		data[i] = FPP(i);
+#endif
 	data[4] = FPP(0);
 	MatDiag<FPP> diagMat (10, data);
 	cout << diagMat.to_string(false,true) << endl;
@@ -25,12 +50,25 @@ int main()
 	cout << "try to convert to a MatSparse" << endl;
 	Faust::MatSparse<FPP,Cpu> matsp(diagMat);
 	cout << matsp.to_string(false,true) << endl;
-	for(int i=0;i<10;i++)
-	{
-		for(int j=0;j<10;j++)
-			cout << MatDense<FPP,Cpu>(matsp)(i,j) << " ";
-		cout << endl;
-	}
+	display_all_elts(diagMat);
+	cout << "conjugate()" << endl;
+	diagMat.conjugate();
+	display_all_elts(diagMat);
+	cout << "vec mul.:" << endl;
+	FPP vec_data[diagMat.getNbCol()];
+	for(int i=0;i<diagMat.getNbCol();i++)
+		vec_data[i] = FPP(i);
+	const Faust::Vect<FPP,Cpu> v(diagMat.getNbCol(), vec_data);
+	Faust::Vect<FPP,Cpu> vr = diagMat.multiply(v);
+	vr.Display();
+	cout << "mat. mul.:" << endl;
+	MatDense<FPP,Cpu>* M = MatDense<FPP,Cpu>::randMat(diagMat.getNbCol(), 4);
+	cout << "matrix to mul." << endl;
+	*M *= FPP(20);
+	display_all_elts(*M);
+	diagMat.multiply(*M, 'N');
+	cout << "matrix multiplied:" << endl;
+	display_all_elts(*M);
 	cout << "end of test_MatDiag" << endl;
 	return EXIT_SUCCESS;
 }

src/faust_linear_operator/CPU/faust_MatDense.h

@@ -167,12 +167,16 @@ void spgemm(const Faust::MatSparse<FPP,Cpu> & A,const Faust::MatDense<FPP,Cpu> &
     class Transform;
     //template<Device DEVICE> class BlasHandle;
 
+    template<typename FPP>
+    class MatDiag;
+
     template<typename FPP>
     class MatDense<FPP,Cpu> : public Faust::MatGeneric<FPP,Cpu>
     {
 
-	friend class MatSparse<FPP,Cpu>;
-        
+		friend class MatSparse<FPP,Cpu>;
+		friend void  MatDiag<FPP>::multiply(MatDense<FPP,Cpu> & M, char opThis) const;
+
 	/// All derived class template of MatDense are considered as friends
         template<class,Device> friend class MatDense;
 

src/faust_linear_operator/CPU/faust_MatDiag.h

@@ -22,16 +22,25 @@ namespace Faust
 
 			public:
 			MatDiag(faust_unsigned_int n): MatGeneric<FPP,Cpu>(n,n), isZeros(true) {}
+			MatDiag(faust_unsigned_int nrows, faust_unsigned_int ncols): MatGeneric<FPP,Cpu>(nrows,ncols), isZeros(true) {}
 
-			MatDiag(faust_unsigned_int n, FPP* data): MatGeneric<FPP,Cpu>(n,n), isZeros(false)
+			MatDiag(faust_unsigned_int n, const FPP* data): MatGeneric<FPP,Cpu>(n,n), isZeros(false)
 			{
 				Matrix<FPP, Dynamic, 1> v(n);
 				memcpy(v.data(), data, sizeof(FPP)*n);
 				mat = v.asDiagonal();
 //				cout << mat.diagonal().size() << endl;
 //				cout << mat.rows() << " " << mat.cols() << endl;
+				isZeros = getNonZeros() == 0;
 			}
 
+			MatDiag(faust_unsigned_int nrows, faust_unsigned_int ncols, const FPP* data): MatDiag<FPP>(min(nrows,ncols), data)
+			{
+				this->dim1 = nrows;
+				this->dim2 = ncols;
+			}
+
+
 			MatType getType() const { return Diag; }
 
 			MatGeneric<FPP,Cpu>* Clone(const bool isOptimize=false) const;
@@ -39,9 +48,9 @@ namespace Faust
 			void multiply(Vect<FPP,Cpu> & vec, char opThis='N') const;
 
 			void multiply(MatDense<FPP,Cpu> & M, char opThis) const;
-			void transpose() { /* nothing to do */ }
-			void conjugate() {  };//TODO:
-			faust_unsigned_int getNonZeros() const { return mat.diagonal().nonZeros(); } //TODO
+			void transpose() { faust_unsigned_int tmp; tmp = this->dim1; this->dim1 = this->dim2; this->dim2 = tmp; }
+			void conjugate() { mat = mat.diagonal().conjugate().asDiagonal(); }
+			faust_unsigned_int getNonZeros() const { return mat.diagonal().nonZeros(); }
 
 			matvar_t* toMatIOVar(bool transpose, bool conjugate) const;
 			FPP normL1(const bool transpose) const;

src/faust_linear_operator/CPU/faust_MatDiag.hpp

@@ -39,20 +39,31 @@ void MatDiag<FPP>::faust_gemm(const Faust::MatDense<FPP,Cpu> & B, Faust::MatDens
 template<typename FPP>
 Vect<FPP,Cpu> MatDiag<FPP>::multiply(const Vect<FPP,Cpu> &v) const
 {
-	Vect<FPP,Cpu> v_(v);
-	return v;
+	Vect<FPP,Cpu> v_(this->getNbRow());
+	v_.vec = mat * v.vec;
+	return v_;
 }
 
 template<typename FPP>
 void  MatDiag<FPP>::multiply(Vect<FPP,Cpu> & vec, char opThis) const
 {
-	//TODO
+	if(opThis == 'T')
+	{
+		MatDiag<FPP> tmat(this->getNbCol(), this->getNbRow(), this->getData());
+		vec.vec = tmat.mat * vec.vec;
+	}
+	else vec.vec = mat * vec.vec;
 }
 
 template<typename FPP>
 void  MatDiag<FPP>::multiply(MatDense<FPP,Cpu> & M, char opThis) const
 {
-	//TODO
+	if (opThis == 'N')
+		M.mat = this->mat * M.mat;
+	else {
+		MatDiag<FPP> tmat(this->getNbCol(), this->getNbRow(), this->getData());
+		M.mat = tmat.mat * M.mat;
+	}
 }
 
 template<typename FPP>

src/faust_linear_operator/CPU/faust_Vect.h

@@ -81,6 +81,10 @@ namespace Faust
     template<typename FPP,Device DEVICE>
     class MatSparse;
 
+
+template<typename FPP>
+    class MatDiag;
+
    // friend function of faust_linear_algebra.h
    template<typename FPP>
    void gemv(const Faust::MatDense<FPP,Cpu> & A,const Faust::Vect<FPP,Cpu> & x,Faust::Vect<FPP,Cpu> & y,const FPP & alpha, const FPP & beta, char typeA);
@@ -94,6 +98,9 @@ namespace Faust
     {
         template<class,Device> friend class Vect;
 		friend double Faust::Transform<FPP,Cpu>::normL1(const bool transpose) const;
+		friend Vect<FPP,Cpu> Faust::MatDiag<FPP>::multiply(const Vect<FPP,Cpu> & vec) const;
+		friend void  MatDiag<FPP>::multiply(Vect<FPP,Cpu> & vec, char opThis) const;
+
 
         public :
         Vect() : dim(0), vec() {}