Make GivensFGFTComplex and GivensFGFTComplexParallel available in matfaust...

Make GivensFGFTComplex and GivensFGFTComplexParallel available in matfaust wrapper (fgft_givens/eigtj).

wrapper/matlab/+matfaust/+fact/fgft_givens.m

@@ -5,7 +5,7 @@
 %> @b Usage
 %>      @b See fact.eigtj
 %>
-%> @param Lap the Laplacian matrix as a numpy array. Must be real and symmetric.
+%> @param Lap the Laplacian matrix (which is symmetric).
 %> @param maxiter see fact.eigtj
 %> @param 'nGivens_per_fac', integer see fact.eigtj
 %> @param 'tol', number see fact.eigtj
@@ -32,9 +32,9 @@ function [FGFT,D] = fgft_givens(Lap, maxiter, varargin)
 	import matfaust.Faust
 	nGivens_per_fac = 1; % default value
 	verbosity = 0; % default value
-	if(~ ismatrix(Lap) || ~ isreal(Lap))
-		error('Lap must be a real matrix.')
-	end
+%	if(~ ismatrix(Lap) || ~ isreal(Lap))
+%		error('Lap must be a real matrix.')
+%	end
 	if(size(Lap,1) ~= size(Lap,2))
 		error('Lap must be square')
 	end
@@ -98,5 +98,5 @@ function [FGFT,D] = fgft_givens(Lap, maxiter, varargin)
 	end
 	[core_obj, D] = mexfgftgivensReal(Lap, maxiter, nGivens_per_fac, verbosity, tol, relerr, order);
 	D = sparse(diag(D));
-	FGFT = Faust(core_obj, true);
+	FGFT = Faust(core_obj, isreal(Lap));
 end

wrapper/matlab/src/mexfgftgivens.cpp.in

@@ -35,6 +35,8 @@
 #include "mex.h"
 #include "faust_GivensFGFT.h"
 #include "faust_GivensFGFTParallel.h"
+#include "faust_GivensFGFTComplex.h"
+#include "faust_GivensFGFTParallelComplex.h"
 #include "faust_TransformHelper.h"
 #include "class_handle.hpp"
 #include "faust_Vect.h"
@@ -46,11 +48,16 @@
 #include <stdexcept>
 
 //normally Givens is only for double, float scalars
+// update: complex matrices are available with GivensFGFTComplex
 typedef @FAUST_SCALAR@ SCALAR;
 typedef @FACT_FPP@ FPP2;
 
 using namespace Faust;
 
+void fgft_givens(const mxArray* matlab_matrix, int J, int t, double tol, unsigned int verbosity, bool relErr, int order, mxArray **plhs);
+
+void fgft_givens_cplx(const mxArray* matlab_matrix, int J, int t, double tol, unsigned int verbosity, bool relErr, int order, mxArray **plhs);
+
 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
 {
 
@@ -78,6 +85,15 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
 
 	const mxArray* matlab_matrix = prhs[0]; // Laplacian
 
+	if(mxIsComplex(matlab_matrix))
+		fgft_givens_cplx(matlab_matrix, J, t, tol, verbosity, relErr, order, plhs);
+	else
+		fgft_givens(matlab_matrix, J, t, tol, verbosity, relErr, order, plhs);
+
+}
+
+void fgft_givens(const mxArray* matlab_matrix, int J, int t, double tol, unsigned int verbosity, bool relErr, int order, mxArray **plhs)
+{
 	// initialization of the matrix that will be factorized
 	Faust::MatGeneric<SCALAR,Cpu>* Lap;
 	Faust::MatDense<SCALAR,Cpu> dLap;
@@ -87,7 +103,6 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
 
 	Faust::GivensFGFT<SCALAR, Cpu, FPP2>* algo;
 
-
 	try{
 		if (mxIsSparse(matlab_matrix))
 		{
@@ -131,6 +146,60 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
 		mexErrMsgTxt(e.what());
 	}
 
-
 }
 
+void fgft_givens_cplx(const mxArray* matlab_matrix, int J, int t, double tol, unsigned int verbosity, bool relErr, int order, mxArray **plhs)
+{
+	// initialization of the matrix that will be factorized
+	Faust::MatGeneric<complex<SCALAR>,Cpu>* Lap;
+	Faust::MatDense<complex<SCALAR>,Cpu> dLap;
+	Faust::MatSparse<complex<SCALAR>,Cpu> sLap;
+	Faust::BlasHandle<Cpu> blas_handle;
+	Faust::SpBlasHandle<Cpu> spblas_handle;
+
+	Faust::GivensFGFTComplex<complex<SCALAR>, Cpu, FPP2>* algo;
+
+	try{
+/*		if (mxIsSparse(matlab_matrix))
+		{
+			mxArray2FaustspMat(matlab_matrix,sLap);
+			Lap = &sLap;
+			//TODO: blas handles to pass later ? (if you want to do the calculation with blas instead of eigen)
+			if(t <= 1)
+				algo = new GivensFGFTComplex<complex<SCALAR>,Cpu,FPP2>(sLap, J, verbosity, tol, relErr);
+			else
+				algo = new GivensFGFTParallel<complex<SCALAR>,Cpu,FPP2>(sLap, J, t, verbosity, tol, relErr);
+
+		}else
+		{*/
+			mxArray2FaustMat(matlab_matrix,dLap);
+			Lap = &dLap;
+			if(t <= 1)
+				algo = new GivensFGFTComplex<complex<SCALAR>,Cpu,FPP2>(dLap, J, verbosity, tol, relErr);
+			else
+				algo = new GivensFGFTParallelComplex<complex<SCALAR>,Cpu,FPP2>(dLap, J, t, verbosity, tol, relErr);
+		//}
+
+
+
+		algo->compute_facts();
+
+		Faust::Vect<complex<SCALAR>,Cpu> D(Lap->getNbRow());
+		// true below is for ascendant order of eigenvalues
+		algo->get_D(const_cast<complex<SCALAR>*>(D.getData()), order); // don't respect constness for optimization (saving a vector copy)
+
+		plhs[1] = FaustVec2mxArray(D);
+
+		Faust::Transform<complex<SCALAR>,Cpu> trans = std::move(algo->get_transform(order));
+		TransformHelper<complex<SCALAR>,Cpu> *th = new TransformHelper<complex<SCALAR>,Cpu>(trans, true); // true is for moving and not copying the Transform object into TransformHelper (optimization possible cause we know the original object won't be used later)
+
+		plhs[0] = convertPtr2Mat<Faust::TransformHelper<complex<SCALAR>, Cpu>>(th);
+
+		delete algo;
+	}
+	catch (const std::exception& e)
+	{
+
+	} 
+
+}