Fix issue #280 on MatDiag multiply (adjoint case) and add a unit test for it.

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

@@ -31,9 +31,33 @@ void display_all_elts(MatDiag<FPP> matdiag)
 	display_all_elts(Faust::MatSparse<FPP,Cpu>(matdiag));
 }
 
+void test_mul_adjoint()
+{
+	cout  << "test_mul_adjoint" << endl;
+	Vect<FPP, Cpu> v(10);
+	v.setRand();
+	MatDiag<FPP> dmat(10, v.getData());
+
+
+	Vect<FPP, Cpu> v2(10);
+	Vect<FPP, Cpu> test_vec(10);
+	v2.setRand();
+	test_vec = v2;
+	dmat.multiply(test_vec, 'H');
+
+	Vect<FPP, Cpu> ref_vec(10);
+	for(int i = 0; i < 10; i++)
+		ref_vec.getData()[i] = conj(v(i)) * v2(i);
+
+	test_vec -= ref_vec;
+	assert(test_vec.norm() < 1e-6);
+	cout  << "OK" << endl;
+}
+
 int main()
 {
 	cout << "start of test_MatDiag" << endl;
+	test_mul_adjoint();
 	FPP data[10];
 	for(int i=0; i < 10; i++)
 #if(@TEST_IS_COMPLEX@==1)

src/faust_linear_operator/CPU/faust_MatDiag.hpp

@@ -44,7 +44,7 @@ template<typename FPP>
 void  Faust::MatDiag<FPP>::multiply(Faust::Vect<FPP,Cpu> & vec, char opThis) const
 {
 	if(opThis = 'H')
-		vec.vec = mat.conjugate() * vec.vec;
+		vec.vec = mat.diagonal().conjugate().array() * vec.vec.array();
 	else //if (opThis == 'N' || opThis == 'T')
 		vec.vec = mat * vec.vec;
 }
@@ -53,7 +53,7 @@ template<typename FPP>
 void  Faust::MatDiag<FPP>::multiply(MatDense<FPP,Cpu> & M, char opThis) const
 {
 	if(opThis = 'H')
-		M.mat = mat.conjugate() * M.mat;
+		M.mat = mat.diagonal().conjugate().asDiagonal() * M.mat;
 	else //if (opThis == 'N' || opThis == 'T')
 		M.mat = mat * M.mat;
 }