Mentions légales du service

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

svdtj C++ version (WIP).

parent c04aaeae
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
wrapper/matlab/CMakeLists.txt
+ 2
− 0
View file @ e34b939b
...@@ -102,12 +102,14 @@ foreach(SCALAR_AND_FSUFFIX double:Real std::complex<double>:Cplx) # TODO: float ...@@ -102,12 +102,14 @@ foreach(SCALAR_AND_FSUFFIX double:Real std::complex<double>:Cplx) # TODO: float
# Givens FGFT is only for double and float scalars # Givens FGFT is only for double and float scalars
if(NOT ${FAUST_SCALAR} MATCHES "complex") if(NOT ${FAUST_SCALAR} MATCHES "complex")
configure_file(${FAUST_MATLAB_MEX_SRC_DIR}/mexfgftgivens.cpp.in ${FAUST_MATLAB_MEX_SRC_DIR}/mexfgftgivens${FSUFFIX}.cpp @ONLY) configure_file(${FAUST_MATLAB_MEX_SRC_DIR}/mexfgftgivens.cpp.in ${FAUST_MATLAB_MEX_SRC_DIR}/mexfgftgivens${FSUFFIX}.cpp @ONLY)
configure_file(${FAUST_MATLAB_MEX_SRC_DIR}/mexsvdtj.cpp.in ${FAUST_MATLAB_MEX_SRC_DIR}/mexsvdtj${FSUFFIX}.cpp @ONLY)
endif() endif()
# copy the *.m for factorization now, because we have the FSUFFIX in hands # copy the *.m for factorization now, because we have the FSUFFIX in hands
configure_file(${FAUST_MATLAB_DOC_SRC_DIR}/mexHierarchical_fact.m.in ${FAUST_MATLAB_DOC_SRC_DIR}/mexHierarchical_fact${FSUFFIX}.m COPYONLY) configure_file(${FAUST_MATLAB_DOC_SRC_DIR}/mexHierarchical_fact.m.in ${FAUST_MATLAB_DOC_SRC_DIR}/mexHierarchical_fact${FSUFFIX}.m COPYONLY)
configure_file(${FAUST_MATLAB_DOC_SRC_DIR}/mexPalm4MSA.m.in ${FAUST_MATLAB_DOC_SRC_DIR}/mexPalm4MSA${FSUFFIX}.m COPYONLY) configure_file(${FAUST_MATLAB_DOC_SRC_DIR}/mexPalm4MSA.m.in ${FAUST_MATLAB_DOC_SRC_DIR}/mexPalm4MSA${FSUFFIX}.m COPYONLY)
if(NOT ${FAUST_SCALAR} MATCHES "complex") if(NOT ${FAUST_SCALAR} MATCHES "complex")
configure_file(${FAUST_MATLAB_DOC_SRC_DIR}/mexfgftgivens.m.in ${FAUST_MATLAB_DOC_SRC_DIR}/mexfgftgivens${FSUFFIX}.m COPYONLY) configure_file(${FAUST_MATLAB_DOC_SRC_DIR}/mexfgftgivens.m.in ${FAUST_MATLAB_DOC_SRC_DIR}/mexfgftgivens${FSUFFIX}.m COPYONLY)
configure_file(${FAUST_MATLAB_DOC_SRC_DIR}/mexsvdtj.m.in ${FAUST_MATLAB_DOC_SRC_DIR}/mexsvdtj${FSUFFIX}.m COPYONLY)
endif() endif()
endforeach() endforeach()
......
wrapper/matlab/doc/mexsvdtj.m.in 0 → 100644
+ 0
− 0
View file @ e34b939b
wrapper/matlab/src/mexfgftgivens.cpp.in
+ 4
− 1
View file @ e34b939b
...@@ -38,6 +38,7 @@ ...@@ -38,6 +38,7 @@
#include "faust_GivensFGFTComplex.h" #include "faust_GivensFGFTComplex.h"
#include "faust_GivensFGFTParallelComplex.h" #include "faust_GivensFGFTParallelComplex.h"
#include "faust_TransformHelper.h" #include "faust_TransformHelper.h"
#include "faust_linear_algebra.h"
#include "class_handle.hpp" #include "class_handle.hpp"
#include "faust_Vect.h" #include "faust_Vect.h"
#include <vector> #include <vector>
...@@ -200,6 +201,8 @@ void fgft_givens_cplx(const mxArray* matlab_matrix, int J, int t, double tol, un ...@@ -200,6 +201,8 @@ void fgft_givens_cplx(const mxArray* matlab_matrix, int J, int t, double tol, un
catch (const std::exception& e) catch (const std::exception& e)
{ {
} }
} }
wrapper/matlab/src/mexsvdtj.cpp.in 0 → 100644
+ 179
− 0
View file @ e34b939b
/****************************************************************************/
/* Description: */
/* For more information on the FAuST Project, please visit the website */
/* of the project : <http://faust.inria.fr> */
/* */
/* License: */
/* Copyright (2019): Hakim HADJ-DJILANI, */
/* Nicolas Bellot, Adrien Leman, Thomas Gautrais, */
/* Luc Le Magoarou, Remi Gribonval */
/* INRIA Rennes, FRANCE */
/* http://www.inria.fr/ */
/* */
/* The FAuST Toolbox is distributed under the terms of INRIA. */
/* */
/* This program is distributed in the hope that it will be useful, but */
/* WITHOUT ANY WARRANTY; without even the implied warranty of */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. */
/* */
/* More information about license: http://faust.inria.fr */
/* */
/* Contacts: */
/* Nicolas Bellot : nicolas.bellot@inria.fr */
/* Adrien Leman : adrien.leman@inria.fr */
/* Thomas Gautrais : thomas.gautrais@inria.fr */
/* Luc Le Magoarou : luc.le-magoarou@inria.fr */
/* Remi Gribonval : remi.gribonval@inria.fr */
/* */
/* References: */
/* [1] Le Magoarou L. and Gribonval R., "Flexible multi-layer sparse */
/* approximations of matrices and applications", Journal of Selected */
/* Topics in Signal Processing, 2016. */
/* <https://hal.archives-ouvertes.fr/hal-01167948v1> */
#include "mex.h"
#include "faust_GivensFGFT.h"
#include "faust_GivensFGFTParallel.h"
#include "faust_GivensFGFTComplex.h"
#include "faust_GivensFGFTParallelComplex.h"
#include "faust_TransformHelper.h"
#include "faust_linear_algebra.h"
#include "class_handle.hpp"
#include "faust_Vect.h"
#include <vector>
#include <string>
#include <algorithm>
#include "mx2Faust.h"
#include "faust2Mx.h"
#include <stdexcept>
typedef @FAUST_SCALAR@ SCALAR;
typedef @FACT_FPP@ FPP2;
using namespace Faust;
//void svdtj_cplx(const mxArray* matlab_matrix, int J, int t, double tol, unsigned int verbosity, bool relErr, int order, mxArray **plhs);
void svdtj(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[])
{
int J;
int t = 1; // default value for non-parallel Givens FGFT
unsigned int verbosity = 0; //default verbosity (no info. displayed)
double tol = 0;
bool relErr = true;
int order;
if(nrhs < 2 || nrhs > 7)
mexErrMsgTxt("Bad Number of inputs arguments");
J = (int) mxGetScalar(prhs[1]);
if(nrhs >= 3)
t = (int) mxGetScalar(prhs[2]);
if(nrhs >= 4)
verbosity = (int) mxGetScalar(prhs[3]);
if(nrhs >= 5)
tol = (double) mxGetScalar(prhs[4]);
if(nrhs >= 6)
relErr = (bool) mxGetScalar(prhs[5]);
if(nrhs >= 7)
order = (int) mxGetScalar(prhs[6]); //eigenvalues order
const mxArray* matlab_matrix = prhs[0]; // Laplacian
// if(mxIsComplex(matlab_matrix))
// svdtj_cplx(matlab_matrix, J, t, tol, verbosity, relErr, order, plhs);
// else
svdtj(matlab_matrix, J, t, tol, verbosity, relErr, order, plhs);
}
void svdtj(const mxArray* matlab_matrix, int J, int t, double tol, unsigned int verbosity, bool relErr, int order, mxArray **plhs)
{
Faust::MatGeneric<SCALAR,Cpu>* M;
Faust::MatDense<SCALAR,Cpu> dM;
Faust::MatDense<SCALAR, Cpu> dM_M, dMM_; // M'*M, M*M'
Faust::MatSparse<SCALAR,Cpu> sM;
Faust::MatDense<SCALAR, Cpu> sM_M, sMM_; // M'*M, M*M'
Faust::BlasHandle<Cpu> blas_handle;
Faust::SpBlasHandle<Cpu> spblas_handle;
Faust::GivensFGFT<SCALAR, Cpu, FPP2>* algoW1;
Faust::GivensFGFT<SCALAR, Cpu, FPP2>* algoW2;
try{
if (mxIsSparse(matlab_matrix))
{
mxArray2FaustspMat(matlab_matrix,sM);
dM = sM; //TODO: optimize
spgemm(sM, dM, dM_M, 1.0, 0.0, 'T', 'N');
spgemm(sM, dM, dMM_, 1.0, 0.0, 'N', 'T');
M = &sM;
if(t <= 1)
{
algoW1 = new GivensFGFT<SCALAR,Cpu,FPP2>(dMM_, J, verbosity, tol, relErr);
algoW2 = new GivensFGFT<SCALAR,Cpu,FPP2>(dM_M, J, verbosity, tol, relErr);
}
else
{
algoW1 = new GivensFGFTParallel<SCALAR,Cpu,FPP2>(dMM_, J, t, verbosity, tol, relErr);
algoW2 = new GivensFGFT<SCALAR,Cpu,FPP2>(dM_M, J, verbosity, tol, relErr);
}
}else
{
mxArray2FaustMat(matlab_matrix, dM);
gemm(dM, dM, dM_M, 1.0, 0.0, 'T', 'N');
gemm(dM, dM, dMM_, 1.0, 0.0, 'N', 'T');
M = &dM;
if(t <= 1)
{
algoW1 = new GivensFGFT<SCALAR,Cpu,FPP2>(dMM_, J, verbosity, tol, relErr);
algoW2 = new GivensFGFT<SCALAR,Cpu,FPP2>(dM_M, J, verbosity, tol, relErr);
}
else
{
algoW1 = new GivensFGFTParallel<SCALAR,Cpu,FPP2>(dMM_, J, t, verbosity, tol, relErr);
algoW2 = new GivensFGFT<SCALAR,Cpu,FPP2>(dM_M, J, verbosity, tol, relErr);
}
}
//TODO: parallelize
algoW1->compute_facts();
algoW2->compute_facts();
Faust::Vect<SCALAR,Cpu> S(M->getNbRow());
Faust::Transform<SCALAR,Cpu> transW1 = std::move(algoW1->get_transform(order));
TransformHelper<SCALAR,Cpu> *thW1 = new TransformHelper<SCALAR,Cpu>(transW1, 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<SCALAR, Cpu>>(thW1);
Faust::Transform<SCALAR,Cpu> transW2 = std::move(algoW1->get_transform(order));
TransformHelper<SCALAR,Cpu> *thW2 = new TransformHelper<SCALAR,Cpu>(transW2, 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[2] = convertPtr2Mat<Faust::TransformHelper<SCALAR, Cpu>>(thW2);
// compute S = W1'*M*W2 = W1'*(W2^T*M)^T
thW2->transpose();
Faust::MatDense<SCALAR,Cpu> MW2 = thW2->multiply(dM, /* transpose */ true);
thW1->transpose();
Faust::MatDense<SCALAR,Cpu> W1_MW2 = thW1->multiply(MW2, /* transpose */ true);
// create diagonal vector
for(int i=0;i<S.size();i++)
S.getData()[i] = W1_MW2(i,i);
plhs[1] = FaustVec2mxArray(S);
delete algoW1;
delete algoW2;
}
catch (const std::exception& e)
{
mexErrMsgTxt(e.what());
}
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment