Mentions légales du service

Skip to content
Snippets Groups Projects
Commit 09c887d5 authored by hhakim's avatar hhakim
Browse files

Make GivensFGFTComplex and GivensFGFTComplexParallel available in matfaust...

Make GivensFGFTComplex and GivensFGFTComplexParallel available in matfaust wrapper (fgft_givens/eigtj).
parent fb2e3f27
No related branches found
No related tags found
No related merge requests found
......@@ -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
......@@ -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)
{
}
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment