From 3b3ca784d15aa5087feba92b671b7016437a79ad Mon Sep 17 00:00:00 2001
From: hhakim <hakim.hadj-djilani@inria.fr>
Date: Fri, 20 Dec 2019 17:24:23 +0100
Subject: [PATCH] Implement a new limitation for eigtj: if the nGivens value
implies a size of transform that doesn't worth it in terms of complexity then
the algorithm is stopped with a warning.
A new eigtj option called enable_large_Faust allows to bypass this limitation.
Updating also the unit tests of both pyfaust and matfaust because they trespassed this limitation.
---
misc/test/src/Matlab/FaustFactoryTest.m | 18 +--
misc/test/src/Python/test_FaustPy.py | 16 ++-
.../factorization/faust_GivensFGFT.h | 4 +-
.../factorization/faust_GivensFGFT.hpp | 4 +-
.../factorization/faust_GivensFGFTComplex.h | 4 +-
.../factorization/faust_GivensFGFTComplex.hpp | 4 +-
.../factorization/faust_GivensFGFTGen.h | 9 +-
.../factorization/faust_GivensFGFTGen.hpp | 20 +++-
.../factorization/faust_GivensFGFTParallel.h | 4 +-
.../faust_GivensFGFTParallel.hpp | 4 +-
.../faust_GivensFGFTParallelComplex.h | 4 +-
.../faust_GivensFGFTParallelComplex.hpp | 8 +-
wrapper/matlab/+matfaust/+fact/eigtj.m | 11 +-
wrapper/matlab/src/mexfgftgivens.cpp.in | 51 +++++----
wrapper/python/pyfaust/fact.py | 10 +-
wrapper/python/src/FaustCoreCy.pxd | 34 +++---
wrapper/python/src/FaustFactGivensFGFT.h | 10 +-
wrapper/python/src/FaustFactGivensFGFT.hpp | 67 ++++++-----
wrapper/python/src/_FaustCorePy.pyx | 104 +++++++++++-------
19 files changed, 229 insertions(+), 157 deletions(-)
diff --git a/misc/test/src/Matlab/FaustFactoryTest.m b/misc/test/src/Matlab/FaustFactoryTest.m
index 80134f293..eb479f24a 100644
--- a/misc/test/src/Matlab/FaustFactoryTest.m
+++ b/misc/test/src/Matlab/FaustFactoryTest.m
@@ -163,7 +163,7 @@ classdef FaustFactoryTest < matlab.unittest.TestCase
import matfaust.*
load([this.faust_paths{1} '../../../misc/data/mat/test_GivensDiag_Lap_U_J.mat'])
% Lap and J available
- [F,D] = matfaust.fact.fgft_givens(Lap, 'nGivens', J);%, 0, 'verbosity', 1);
+ [F,D] = matfaust.fact.fgft_givens(Lap, 'nGivens', J, 'enable_large_Faust', true, 'nGivens_per_fac', 1);%, 0, 'verbosity', 1);
this.verifyEqual(size(F), size(Lap))
%disp('norm F: ')
%norm(F, 'fro')
@@ -173,7 +173,7 @@ classdef FaustFactoryTest < matlab.unittest.TestCase
% misc/test/src/C++/GivensFGFT.cpp.in
this.verifyEqual(err, 0.0326529, 'AbsTol', 0.00001)
% verify it works the same using the eigtj() alias function
- [F2,D2] = matfaust.fact.eigtj(Lap, 'nGivens', J);%, 0, 'verbosity', 2);
+ [F2,D2] = matfaust.fact.eigtj(Lap, 'nGivens', J, 'enable_large_Faust', true, 'nGivens_per_fac', 1);%, 0, 'verbosity', 2);
this.verifyEqual(full(F2),full(F))
this.verifyEqual(D,D2)
end
@@ -183,7 +183,7 @@ classdef FaustFactoryTest < matlab.unittest.TestCase
import matfaust.*
load([this.faust_paths{1} '../../../misc/data/mat/test_GivensDiag_Lap_U_J.mat'])
% Lap and J available
- [F,D] = matfaust.fact.fgft_givens(sparse(Lap), 'nGivens', J);%, 0, 'verbosity', 1);
+ [F,D] = matfaust.fact.fgft_givens(sparse(Lap), 'nGivens', J, 'enable_large_Faust', true, 'nGivens_per_fac', 1);%, 0, 'verbosity', 1);
this.verifyEqual(size(F), size(Lap))
%disp('norm F: ')
%norm(F, 'fro')
@@ -193,18 +193,18 @@ classdef FaustFactoryTest < matlab.unittest.TestCase
% misc/test/src/C++/GivensFGFT.cpp.in
this.verifyEqual(err, 0.0326529, 'AbsTol', 0.00001)
% verify it works the same using the eigtj() alias function
- [F2,D2] = matfaust.fact.eigtj(Lap, 'nGivens', J);%, 0, 'verbosity', 2);
+ [F2,D2] = matfaust.fact.eigtj(Lap, 'nGivens', J, 'nGivens_per_fac', 1, 'enable_large_Faust', true);%, 0, 'verbosity', 2);
this.verifyEqual(full(F2),full(F))
this.verifyEqual(D,D2)
end
function test_fgft_givens_parallel(this)
- disp('Test matfaust.fact.fgft_givens_sparse() -- parallel')
+ disp('Test matfaust.fact.fgft_givens() -- parallel')
import matfaust.*
load([this.faust_paths{1} '../../../misc/data/mat/test_GivensDiag_Lap_U_J.mat'])
% Lap and J available
t = size(Lap,1)/2;
- [F,D] = matfaust.fact.fgft_givens(Lap, 'nGivens', J, 'nGivens_per_fac', t); %, 'verbosity', 2);
+ [F,D] = matfaust.fact.fgft_givens(Lap, 'nGivens', J, 'nGivens_per_fac', t, 'enable_large_Faust', true); %, 'verbosity', 2);
this.verifyEqual(size(F), size(Lap))
%disp('norm F: ')
%norm(F, 'fro')
@@ -214,7 +214,7 @@ classdef FaustFactoryTest < matlab.unittest.TestCase
% misc/test/src/C++/GivensFGFTParallel.cpp.in
this.verifyEqual(err,0.0398154, 'AbsTol', 0.00001)
% verify it works the same using the eigtj() alias function
- [F2,D2] = matfaust.fact.eigtj(Lap, 'nGivens', J, 'nGivens_per_fac', t); %, 'verbosity', 2);
+ [F2,D2] = matfaust.fact.eigtj(Lap, 'nGivens', J, 'nGivens_per_fac', t, 'enable_large_Faust', true); %, 'verbosity', 2);
this.verifyEqual(full(F2),full(F))
this.verifyEqual(D,D2)
end
@@ -225,7 +225,7 @@ classdef FaustFactoryTest < matlab.unittest.TestCase
load([this.faust_paths{1} '../../../misc/data/mat/test_GivensDiag_Lap_U_J.mat'])
% Lap and J available
t = size(Lap,1)/2;
- [F,D] = matfaust.fact.fgft_givens(sparse(Lap), 'nGivens', J, 'nGivens_per_fac', t); %, 'verbosity', 2);
+ [F,D] = matfaust.fact.fgft_givens(sparse(Lap), 'nGivens', J, 'nGivens_per_fac', t, 'enable_large_Faust', true); %, 'verbosity', 2);
this.verifyEqual(size(F), size(Lap))
%disp('norm F: ')
%norm(F, 'fro')
@@ -235,7 +235,7 @@ classdef FaustFactoryTest < matlab.unittest.TestCase
% misc/test/src/C++/GivensFGFTParallel.cpp.in
this.verifyEqual(err, 0.0398154, 'AbsTol', 0.00001)
% verify it works the same using the eigtj() alias function
- [F2,D2] = matfaust.fact.eigtj(Lap, 'nGivens', J, 'nGivens_per_fac', t); %, 'verbosity', 2);
+ [F2,D2] = matfaust.fact.eigtj(Lap, 'nGivens', J, 'nGivens_per_fac', t, 'enable_large_Faust', true); %, 'verbosity', 2);
this.verifyEqual(full(F2),full(F))
this.verifyEqual(D,D2)
end
diff --git a/misc/test/src/Python/test_FaustPy.py b/misc/test/src/Python/test_FaustPy.py
index 3ddb1bfbe..f34c88239 100644
--- a/misc/test/src/Python/test_FaustPy.py
+++ b/misc/test/src/Python/test_FaustPy.py
@@ -925,7 +925,7 @@ class TestFaustFactory(unittest.TestCase):
L = L.astype(np.float64)
J = \
int(loadmat(sys.path[0]+"/../../../misc/data/mat/test_GivensDiag_Lap_U_J.mat")['J'])
- D, F = pyfaust.fact.fgft_givens(L, J, nGivens_per_fac=1, verbosity=0)
+ D, F = pyfaust.fact.fgft_givens(L, J, nGivens_per_fac=1, verbosity=0, enable_large_Faust=True)
D = spdiags(D, [0], L.shape[0], L.shape[0])
print("Lap norm:", norm(L, 'fro'))
err = norm((F*D.todense())*F.T.todense()-L,"fro")/norm(L,"fro")
@@ -946,7 +946,8 @@ class TestFaustFactory(unittest.TestCase):
J = \
int(loadmat(sys.path[0]+"/../../../misc/data/mat/test_GivensDiag_Lap_U_J.mat")['J'])
t = int(L.shape[0]/2)
- D, F = fgft_givens(L, J, nGivens_per_fac=t, verbosity=0)
+ D, F = fgft_givens(L, J, nGivens_per_fac=t, verbosity=0,
+ enable_large_Faust=True)
D = spdiags(D, [0], L.shape[0], L.shape[0])
print("Lap norm:", norm(L, 'fro'))
err = norm((F*D.todense())*F.T.todense()-L,"fro")/norm(L,"fro")
@@ -954,7 +955,8 @@ class TestFaustFactory(unittest.TestCase):
# the error reference is from the C++ test,
# misc/test/src/C++/GivensFGFTParallel.cpp.in (_double version)
self.assertAlmostEqual(err, 0.0398154, places=7)
- D2, F2 = eigtj(L, J, nGivens_per_fac=t, verbosity=0)
+ D2, F2 = eigtj(L, J, nGivens_per_fac=t, verbosity=0,
+ enable_large_Faust=True)
D2 = spdiags(D, [0], L.shape[0], L.shape[0])
print("Lap norm:", norm(L, 'fro'))
err2 = norm((F2*D.todense())*F2.T.todense()-L,"fro")/norm(L,"fro")
@@ -977,7 +979,8 @@ class TestFaustFactory(unittest.TestCase):
L = L.astype(np.float64)
J = \
int(loadmat(sys.path[0]+"/../../../misc/data/mat/test_GivensDiag_Lap_U_J.mat")['J'])
- D, F = pyfaust.fact.fgft_givens(csr_matrix(L), J, nGivens_per_fac=0, verbosity=0)
+ D, F = pyfaust.fact.fgft_givens(csr_matrix(L), J, nGivens_per_fac=0,
+ verbosity=0,enable_large_Faust=True)
D = spdiags(D, [0], L.shape[0], L.shape[0])
print("Lap norm:", norm(L, 'fro'))
err = norm((F*D.todense())*F.T.todense()-L,"fro")/norm(L,"fro")
@@ -999,7 +1002,7 @@ class TestFaustFactory(unittest.TestCase):
J = \
int(loadmat(sys.path[0]+"/../../../misc/data/mat/test_GivensDiag_Lap_U_J.mat")['J'])
t = int(L.shape[0]/2)
- D, F = fgft_givens(csr_matrix(L), J, nGivens_per_fac=t, verbosity=0)
+ D, F = fgft_givens(csr_matrix(L), J, nGivens_per_fac=t, verbosity=0, enable_large_Faust=True)
D = spdiags(D, [0], L.shape[0], L.shape[0])
print("Lap norm:", norm(L, 'fro'))
err = norm((F*D.todense())*F.T.todense()-L,"fro")/norm(L,"fro")
@@ -1007,7 +1010,8 @@ class TestFaustFactory(unittest.TestCase):
# the error reference is from the C++ test,
# misc/test/src/C++/GivensFGFTParallel.cpp.in (_double version)
self.assertAlmostEqual(err, 0.0398154, places=7)
- D2, F2 = eigtj(csr_matrix(L), J, nGivens_per_fac=t, verbosity=0)
+ D2, F2 = eigtj(csr_matrix(L), J, nGivens_per_fac=t, verbosity=0,
+ enable_large_Faust=True)
D2 = spdiags(D, [0], L.shape[0], L.shape[0])
print("Lap norm:", norm(L, 'fro'))
err2 = norm((F2*D.todense())*F2.T.todense()-L,"fro")/norm(L,"fro")
diff --git a/src/algorithm/factorization/faust_GivensFGFT.h b/src/algorithm/factorization/faust_GivensFGFT.h
index 255100f32..08a0c22da 100644
--- a/src/algorithm/factorization/faust_GivensFGFT.h
+++ b/src/algorithm/factorization/faust_GivensFGFT.h
@@ -56,8 +56,8 @@ namespace Faust {
* \param Lap The Laplacian matrix to approximate/diagonalize.
* \param J The number of iterations, Givens rotations factors.
* */
- GivensFGFT(Faust::MatSparse<FPP,DEVICE>& Lap, int J, unsigned int verbosity = 0, const double stoppingError = 0.0, const bool errIsRel = true);
- GivensFGFT(Faust::MatDense<FPP,DEVICE>& Lap, int J, unsigned int verbosity = 0, const double stoppingError = 0.0, const bool errIsRel = true);
+ GivensFGFT(Faust::MatSparse<FPP,DEVICE>& Lap, int J, unsigned int verbosity = 0, const double stoppingError = 0.0, const bool errIsRel = true, const bool enable_large_Faust = false);
+ GivensFGFT(Faust::MatDense<FPP,DEVICE>& Lap, int J, unsigned int verbosity = 0, const double stoppingError = 0.0, const bool errIsRel = true, const bool enable_large_Faust = false);
/** Destructor */
virtual ~GivensFGFT() {/*delete[] q_candidates; delete L;*/};
diff --git a/src/algorithm/factorization/faust_GivensFGFT.hpp b/src/algorithm/factorization/faust_GivensFGFT.hpp
index ec920ad66..19e7a90fe 100644
--- a/src/algorithm/factorization/faust_GivensFGFT.hpp
+++ b/src/algorithm/factorization/faust_GivensFGFT.hpp
@@ -379,7 +379,7 @@ void GivensFGFT<FPP,DEVICE,FPP2>::update_err()
}
template<typename FPP, Device DEVICE, typename FPP2>
-GivensFGFT<FPP,DEVICE,FPP2>::GivensFGFT(Faust::MatSparse<FPP,DEVICE>& Lap, int J, unsigned int verbosity /* deft val == 0 */, const double stoppingError /* default to 0.0 */, const bool errIsRel) : GivensFGFTGen<FPP,DEVICE,FPP2>(Lap, J, verbosity, stoppingError, errIsRel), C(Lap.getNbRow(), Lap.getNbCol()), always_theta2(false)
+GivensFGFT<FPP,DEVICE,FPP2>::GivensFGFT(Faust::MatSparse<FPP,DEVICE>& Lap, int J, unsigned int verbosity /* deft val == 0 */, const double stoppingError /* default to 0.0 */, const bool errIsRel, const bool enable_large_Faust/* deft to false */) : GivensFGFTGen<FPP,DEVICE,FPP2>(Lap, J, verbosity, stoppingError, errIsRel, enable_large_Faust), C(Lap.getNbRow(), Lap.getNbCol()), always_theta2(false)
{
//see parent ctor
C.setOnes();
@@ -387,7 +387,7 @@ GivensFGFT<FPP,DEVICE,FPP2>::GivensFGFT(Faust::MatSparse<FPP,DEVICE>& Lap, int J
}
template<typename FPP, Device DEVICE, typename FPP2>
-GivensFGFT<FPP,DEVICE,FPP2>::GivensFGFT(Faust::MatDense<FPP,DEVICE>& Lap, int J, unsigned int verbosity /* deft val == 0 */, const double stoppingError, const bool errIsRel) : GivensFGFTGen<FPP,DEVICE,FPP2>(Lap, J, verbosity, stoppingError, errIsRel), C(Lap.getNbRow(), Lap.getNbCol()), always_theta2(false)
+GivensFGFT<FPP,DEVICE,FPP2>::GivensFGFT(Faust::MatDense<FPP,DEVICE>& Lap, int J, unsigned int verbosity /* deft val == 0 */, const double stoppingError, const bool errIsRel, const bool enable_large_Faust/* deft to false */) : GivensFGFTGen<FPP,DEVICE,FPP2>(Lap, J, verbosity, stoppingError, errIsRel, enable_large_Faust), C(Lap.getNbRow(), Lap.getNbCol()), always_theta2(false)
{
// see parent ctor
C.setOnes();
diff --git a/src/algorithm/factorization/faust_GivensFGFTComplex.h b/src/algorithm/factorization/faust_GivensFGFTComplex.h
index 3c8078a74..15f91b1ec 100644
--- a/src/algorithm/factorization/faust_GivensFGFTComplex.h
+++ b/src/algorithm/factorization/faust_GivensFGFTComplex.h
@@ -53,8 +53,8 @@ namespace Faust {
public:
const static unsigned int ERROR_CALC_PERIOD = 100;
- GivensFGFTComplex(Faust::MatSparse<FPP,DEVICE>& Lap, int J, unsigned int verbosity = 0, const double stoppingError = 0.0, const bool errIsRel = true);
- GivensFGFTComplex(Faust::MatDense<FPP,DEVICE>& Lap, int J, unsigned int verbosity = 0, const double stoppingError = 0.0, const bool errIsRel = true);
+ GivensFGFTComplex(Faust::MatSparse<FPP,DEVICE>& Lap, int J, unsigned int verbosity = 0, const double stoppingError = 0.0, const bool errIsRel = true, const bool enable_large_Faust = false);
+ GivensFGFTComplex(Faust::MatDense<FPP,DEVICE>& Lap, int J, unsigned int verbosity = 0, const double stoppingError = 0.0, const bool errIsRel = true, const bool enable_large_Faust = false);
virtual void next_step();
diff --git a/src/algorithm/factorization/faust_GivensFGFTComplex.hpp b/src/algorithm/factorization/faust_GivensFGFTComplex.hpp
index b44cb3719..fd513a9c2 100644
--- a/src/algorithm/factorization/faust_GivensFGFTComplex.hpp
+++ b/src/algorithm/factorization/faust_GivensFGFTComplex.hpp
@@ -445,14 +445,14 @@ void GivensFGFTComplex<FPP,DEVICE,FPP2>::update_err()
}
template<typename FPP, Device DEVICE, typename FPP2>
-GivensFGFTComplex<FPP,DEVICE,FPP2>::GivensFGFTComplex(Faust::MatSparse<FPP,DEVICE>& Lap, int J, unsigned int verbosity /* deft val == 0 */, const double stoppingError /* default to 0.0 */, const bool errIsRel) : Faust::GivensFGFTGen<typename FPP::value_type, DEVICE, FPP2, FPP>(Lap, J, verbosity, stoppingError, errIsRel), C(Lap.getNbRow(), Lap.getNbCol())
+GivensFGFTComplex<FPP,DEVICE,FPP2>::GivensFGFTComplex(Faust::MatSparse<FPP,DEVICE>& Lap, int J, unsigned int verbosity /* deft val == 0 */, const double stoppingError /* default to 0.0 */, const bool errIsRel, const bool enable_large_Faust) : Faust::GivensFGFTGen<typename FPP::value_type, DEVICE, FPP2, FPP>(Lap, J, verbosity, stoppingError, errIsRel, enable_large_Faust), C(Lap.getNbRow(), Lap.getNbCol())
{
//see parent ctor
this->C.setZeros();
}
template<typename FPP, Device DEVICE, typename FPP2>
-GivensFGFTComplex<FPP,DEVICE,FPP2>::GivensFGFTComplex(Faust::MatDense<FPP,DEVICE>& Lap, int J, unsigned int verbosity /* deft val == 0 */, const double stoppingError, const bool errIsRel) : Faust::GivensFGFTGen<typename FPP::value_type, DEVICE, FPP2, FPP>(Lap, J, verbosity, stoppingError, errIsRel), C(Lap.getNbRow(), Lap.getNbCol())
+GivensFGFTComplex<FPP,DEVICE,FPP2>::GivensFGFTComplex(Faust::MatDense<FPP,DEVICE>& Lap, int J, unsigned int verbosity /* deft val == 0 */, const double stoppingError, const bool errIsRel, const bool enable_large_Faust) : Faust::GivensFGFTGen<typename FPP::value_type, DEVICE, FPP2, FPP>(Lap, J, verbosity, stoppingError, errIsRel, enable_large_Faust), C(Lap.getNbRow(), Lap.getNbCol())
{
// see parent ctor
this->C.setZeros();
diff --git a/src/algorithm/factorization/faust_GivensFGFTGen.h b/src/algorithm/factorization/faust_GivensFGFTGen.h
index 7aa756da7..903d3f8da 100644
--- a/src/algorithm/factorization/faust_GivensFGFTGen.h
+++ b/src/algorithm/factorization/faust_GivensFGFTGen.h
@@ -107,7 +107,8 @@ namespace Faust {
double stoppingError;
/** \brief true if the stopping error is taken as relative error (absolute otherwise). */
bool errIsRel;
-
+ /** \brief (false to default) true to force the computation of the transform even if it doesn't worth it in term of complexity */
+ bool enable_large_Faust;
public:
@@ -117,10 +118,10 @@ namespace Faust {
* \param J The number of iterations, Givens rotations factors.
* TODO: complete argument list
* */
- GivensFGFTGen(Faust::MatGeneric<FPP4,DEVICE>* Lap, int J, unsigned int verbosity /* deft val == 0 */, const double stoppingError, const bool errIsRel);
+ GivensFGFTGen(Faust::MatGeneric<FPP4,DEVICE>* Lap, int J, unsigned int verbosity /* deft val == 0 */, const double stoppingError, const bool errIsRel, const bool enable_large_Faust = false);
- GivensFGFTGen(Faust::MatSparse<FPP4, DEVICE> & Lap, int J, unsigned int verbosity /* deft val == 0 */, const double stoppingError, const bool errIsRel);
- GivensFGFTGen(Faust::MatDense<FPP4, DEVICE> & Lap, int J, unsigned int verbosity /* deft val == 0 */, const double stoppingError, const bool errIsRel);
+ GivensFGFTGen(Faust::MatSparse<FPP4, DEVICE> & Lap, int J, unsigned int verbosity /* deft val == 0 */, const double stoppingError, const bool errIsRel, const bool enable_large_Faust = false);
+ GivensFGFTGen(Faust::MatDense<FPP4, DEVICE> & Lap, int J, unsigned int verbosity /* deft val == 0 */, const double stoppingError, const bool errIsRel, const bool enable_large_Faust = false);
/** Destructor */
virtual ~GivensFGFTGen() {delete[] q_candidates; delete L;};
diff --git a/src/algorithm/factorization/faust_GivensFGFTGen.hpp b/src/algorithm/factorization/faust_GivensFGFTGen.hpp
index f4550e2cf..f8f75935e 100644
--- a/src/algorithm/factorization/faust_GivensFGFTGen.hpp
+++ b/src/algorithm/factorization/faust_GivensFGFTGen.hpp
@@ -446,12 +446,18 @@ void GivensFGFTGen<FPP,DEVICE,FPP2,FPP4>::compute_facts()
is_D_ordered = false; // facts (re)computed then D must be reordered
ite = 0;
bool stopping = false;
- while(J == 0 || ite < facts.size()) // when J == 0 the stopping criterion is the error against Lap
+ if(stopping = !enable_large_Faust && ! stoppingCritIsError && dim_size*dim_size <= J*4)
+ {
+ cerr << "WARNING: the eigtj algorithm stopped because the transform to be computed doesn't worth it according to its complexity (in space and time) relatively to the size of the matrix to decompose. Still, if you want to force the computation, please use the enable_large_Faust flag." << endl;
+ cerr << "nGivens: " << J << " matrix size: " << dim_size << endl;
+ facts.resize(0);
+ }
+ while(! stopping && (J == 0 || ite < facts.size())) // when J == 0 the stopping criterion is the error against Lap
{
next_step();
ite++;
if(stopping = (ite > 1 && stoppingCritIsError && errs.size() > 2 && errs[ite-1]-errs[ite-2] > FLT_EPSILON))
- if(verbosity>0) cerr << "warning: the algorithm stopped because the last error is greater than the previous one." << endl;
+ /*if(verbosity>0) */cerr << "WARNING: the eigtj algorithm stopped because the last error is greater than the previous one." << endl;
if(stopping || stoppingCritIsError && errs.size() > 0 && (*(errs.end()-1) - stoppingError ) < FLT_EPSILON)
{
facts.resize(ite);
@@ -529,8 +535,8 @@ void GivensFGFTGen<FPP,DEVICE,FPP2,FPP4>::compute_facts()
//
template<typename FPP, Device DEVICE, typename FPP2, typename FPP4>
-GivensFGFTGen<FPP,DEVICE,FPP2,FPP4>::GivensFGFTGen(MatGeneric<FPP4,DEVICE>* Lap, int J, unsigned int verbosity /* deft val == 0 */, const double stoppingError, const bool errIsRel) :
-Lap(*Lap), D(Lap->getNbRow()), errs(0), coord_choices(0), q_candidates(new int[Lap->getNbRow()]), is_D_ordered(false), verbosity(verbosity), stoppingCritIsError(stoppingError != 0.0), stoppingError(stoppingError), errIsRel(errIsRel), Lap_squared_fro_norm(0), facts(J>0?J:1), last_fact_permuted(false), J(J), dim_size(Lap->getNbRow())
+GivensFGFTGen<FPP,DEVICE,FPP2,FPP4>::GivensFGFTGen(MatGeneric<FPP4,DEVICE>* Lap, int J, unsigned int verbosity /* deft val == 0 */, const double stoppingError, const bool errIsRel, const bool enable_large_Faust /* deft to false */) :
+Lap(*Lap), D(Lap->getNbRow()), errs(0), coord_choices(0), q_candidates(new int[Lap->getNbRow()]), is_D_ordered(false), verbosity(verbosity), stoppingCritIsError(stoppingError != 0.0), stoppingError(stoppingError), errIsRel(errIsRel), Lap_squared_fro_norm(0), facts(J>0?(J*4<Lap->getNbRow()*Lap->getNbRow()||enable_large_Faust?J:0):1 /* don't allocate if the complexity doesn't worth it and enable_large_Faust is false*/), last_fact_permuted(false), J(J), dim_size(Lap->getNbRow()), enable_large_Faust(enable_large_Faust)
{
if(Lap->getNbCol() != Lap->getNbRow())
handleError("Faust::GivensFGFTComplex", "Laplacian must be a square matrix.");
@@ -549,13 +555,13 @@ Lap(*Lap), D(Lap->getNbRow()), errs(0), coord_choices(0), q_candidates(new int[
}
template<typename FPP, Device DEVICE, typename FPP2, typename FPP4>
-GivensFGFTGen<FPP,DEVICE,FPP2,FPP4>::GivensFGFTGen(Faust::MatSparse<FPP4, DEVICE> & Lap, int J, unsigned int verbosity /* deft val == 0 */, const double stoppingError, const bool errIsRel) : GivensFGFTGen<FPP,DEVICE,FPP2,FPP4>(&Lap, J, verbosity, stoppingError, errIsRel)
+GivensFGFTGen<FPP,DEVICE,FPP2,FPP4>::GivensFGFTGen(Faust::MatSparse<FPP4, DEVICE> & Lap, int J, unsigned int verbosity /* deft val == 0 */, const double stoppingError, const bool errIsRel, const bool enable_large_Faust/* deft to false */) : GivensFGFTGen<FPP,DEVICE,FPP2,FPP4>(&Lap, J, verbosity, stoppingError, errIsRel, enable_large_Faust)
{
L = new MatSparse<FPP4,DEVICE>(Lap);
}
template<typename FPP, Device DEVICE, typename FPP2, typename FPP4>
-GivensFGFTGen<FPP,DEVICE,FPP2,FPP4>::GivensFGFTGen(Faust::MatDense<FPP4, DEVICE> & Lap, int J, unsigned int verbosity /* deft val == 0 */, const double stoppingError, const bool errIsRel) : GivensFGFTGen<FPP,DEVICE,FPP2,FPP4>(&Lap, J, verbosity, stoppingError, errIsRel)
+GivensFGFTGen<FPP,DEVICE,FPP2,FPP4>::GivensFGFTGen(Faust::MatDense<FPP4, DEVICE> & Lap, int J, unsigned int verbosity /* deft val == 0 */, const double stoppingError, const bool errIsRel, const bool enable_large_Faust/* deft to false */) : GivensFGFTGen<FPP,DEVICE,FPP2,FPP4>(&Lap, J, verbosity, stoppingError, errIsRel, enable_large_Faust)
{
L = new MatDense<FPP4,DEVICE>(Lap);
}
@@ -695,6 +701,8 @@ template<typename FPP, Device DEVICE, typename FPP2, typename FPP4>
Faust::Transform<FPP4,DEVICE> GivensFGFTGen<FPP,DEVICE,FPP2,FPP4>::get_transform(int ord)
{
//TODO: an optimization is possible by changing type of facts to vector<MatGeneric*> it would avoid copying facts into Transform and rather use them directly. It will need a destructor that deletes them eventually if they weren't transfered to a Transform object before.
+ if(facts.size() == 0)
+ throw out_of_range("The transform is empty. The algorithm stopped before computing any factor.");
MatSparse<FPP4,DEVICE> & last_fact = *(facts.end()-1);
MatSparse<FPP4,DEVICE> P(last_fact.getNbCol(), last_fact.getNbCol()); //last_fact permutation matrix
// (to reorder eigenvector transform according to ordered D)
diff --git a/src/algorithm/factorization/faust_GivensFGFTParallel.h b/src/algorithm/factorization/faust_GivensFGFTParallel.h
index 8e51c7f42..cc65eceb1 100644
--- a/src/algorithm/factorization/faust_GivensFGFTParallel.h
+++ b/src/algorithm/factorization/faust_GivensFGFTParallel.h
@@ -84,8 +84,8 @@ namespace Faust {
* \param stoppingError defines a stopping criterion based on error (absolute relative error).
*
*/
- GivensFGFTParallel(Faust::MatDense<FPP,DEVICE>& Lap, int J, int t, unsigned int verbosity = 0, const double stoppingCritIsError = 0.0, const bool errIsRel = true);
- GivensFGFTParallel(Faust::MatSparse<FPP,DEVICE>& Lap, int J, int t, unsigned int verbosity = 0, const double stoppingCritIsError = 0.0, const bool errIsRel = true);
+ GivensFGFTParallel(Faust::MatDense<FPP,DEVICE>& Lap, int J, int t, unsigned int verbosity = 0, const double stoppingCritIsError = 0.0, const bool errIsRel = true, const bool enable_large_Faust = false);
+ GivensFGFTParallel(Faust::MatSparse<FPP,DEVICE>& Lap, int J, int t, unsigned int verbosity = 0, const double stoppingCritIsError = 0.0, const bool errIsRel = true, const bool enable_large_Faust = false);
};
diff --git a/src/algorithm/factorization/faust_GivensFGFTParallel.hpp b/src/algorithm/factorization/faust_GivensFGFTParallel.hpp
index 55b745461..dcf12a2d9 100644
--- a/src/algorithm/factorization/faust_GivensFGFTParallel.hpp
+++ b/src/algorithm/factorization/faust_GivensFGFTParallel.hpp
@@ -6,7 +6,7 @@ using namespace Faust;
#endif
template<typename FPP, Device DEVICE, typename FPP2>
-GivensFGFTParallel<FPP,DEVICE,FPP2>::GivensFGFTParallel(Faust::MatDense<FPP,DEVICE>& Lap, int J, int t, unsigned int verbosity, const double stoppingError, const bool errIsRel) : GivensFGFT<FPP,DEVICE,FPP2>(Lap, J, verbosity, stoppingError, errIsRel), t(t), fact_nrots(0)
+GivensFGFTParallel<FPP,DEVICE,FPP2>::GivensFGFTParallel(Faust::MatDense<FPP,DEVICE>& Lap, int J, int t, unsigned int verbosity, const double stoppingError, const bool errIsRel, const bool enable_large_Faust) : GivensFGFT<FPP,DEVICE,FPP2>(Lap, J, verbosity, stoppingError, errIsRel, enable_large_Faust), t(t), fact_nrots(0)
{
if(J > 0) this->facts.resize(round(J/(float)t));
this->always_theta2 = true;
@@ -14,7 +14,7 @@ GivensFGFTParallel<FPP,DEVICE,FPP2>::GivensFGFTParallel(Faust::MatDense<FPP,DEVI
}
template<typename FPP, Device DEVICE, typename FPP2>
-GivensFGFTParallel<FPP,DEVICE,FPP2>::GivensFGFTParallel(Faust::MatSparse<FPP,DEVICE>& Lap, int J, int t, unsigned int verbosity, const double stoppingError, const bool errIsRel) : GivensFGFT<FPP,DEVICE,FPP2>(Lap, J, verbosity, stoppingError, errIsRel), t(t), fact_nrots(0)
+GivensFGFTParallel<FPP,DEVICE,FPP2>::GivensFGFTParallel(Faust::MatSparse<FPP,DEVICE>& Lap, int J, int t, unsigned int verbosity, const double stoppingError, const bool errIsRel, const bool enable_large_Faust) : GivensFGFT<FPP,DEVICE,FPP2>(Lap, J, verbosity, stoppingError, errIsRel, enable_large_Faust), t(t), fact_nrots(0)
{
if(J > 0) this->facts.resize(round(J/(float)t));
this->always_theta2 = true;
diff --git a/src/algorithm/factorization/faust_GivensFGFTParallelComplex.h b/src/algorithm/factorization/faust_GivensFGFTParallelComplex.h
index 67c121d02..4bb116b56 100644
--- a/src/algorithm/factorization/faust_GivensFGFTParallelComplex.h
+++ b/src/algorithm/factorization/faust_GivensFGFTParallelComplex.h
@@ -86,8 +86,8 @@ namespace Faust {
* \param stoppingError defines a stopping criterion based on error (absolute relative error).
*
*/
- GivensFGFTParallelComplex(Faust::MatDense<FPP,DEVICE>& Lap, int J, int t, unsigned int verbosity = 0, const double stoppingCritIsError = 0.0, const bool errIsRel = true);
- GivensFGFTParallelComplex(Faust::MatSparse<FPP,DEVICE>& Lap, int J, int t, unsigned int verbosity = 0, const double stoppingCritIsError = 0.0, const bool errIsRel = true);
+ GivensFGFTParallelComplex(Faust::MatDense<FPP,DEVICE>& Lap, int J, int t, unsigned int verbosity = 0, const double stoppingCritIsError = 0.0, const bool errIsRel = true, const bool enable_large_Faust = false);
+ GivensFGFTParallelComplex(Faust::MatSparse<FPP,DEVICE>& Lap, int J, int t, unsigned int verbosity = 0, const double stoppingCritIsError = 0.0, const bool errIsRel = true, const bool enable_large_Faust = false);
};
diff --git a/src/algorithm/factorization/faust_GivensFGFTParallelComplex.hpp b/src/algorithm/factorization/faust_GivensFGFTParallelComplex.hpp
index 7c71d4cf3..a460341ce 100644
--- a/src/algorithm/factorization/faust_GivensFGFTParallelComplex.hpp
+++ b/src/algorithm/factorization/faust_GivensFGFTParallelComplex.hpp
@@ -6,16 +6,16 @@ using namespace Faust;
#endif
template<typename FPP, Device DEVICE, typename FPP2>
-GivensFGFTParallelComplex<FPP,DEVICE,FPP2>::GivensFGFTParallelComplex(Faust::MatDense<FPP,DEVICE>& Lap, int J, int t, unsigned int verbosity, const double stoppingError, const bool errIsRel) : GivensFGFTComplex<FPP,DEVICE,FPP2>(Lap, J, verbosity, stoppingError, errIsRel), t(t), fact_nrots(0)
+GivensFGFTParallelComplex<FPP,DEVICE,FPP2>::GivensFGFTParallelComplex(Faust::MatDense<FPP,DEVICE>& Lap, int J, int t, unsigned int verbosity, const double stoppingError, const bool errIsRel, const bool enable_large_Faust) : GivensFGFTComplex<FPP,DEVICE,FPP2>(Lap, J, verbosity, stoppingError, errIsRel, enable_large_Faust), t(t), fact_nrots(0)
{
- this->facts.resize(round(J/(float)t));
+ if(J > 0) this->facts.resize(round(J/(float)t));
this->coord_choices.resize(0);
}
template<typename FPP, Device DEVICE, typename FPP2>
-GivensFGFTParallelComplex<FPP,DEVICE,FPP2>::GivensFGFTParallelComplex(Faust::MatSparse<FPP,DEVICE>& Lap, int J, int t, unsigned int verbosity, const double stoppingError, const bool errIsRel) : GivensFGFTComplex<FPP,DEVICE,FPP2>(Lap, J, verbosity, stoppingError, errIsRel), t(t), fact_nrots(0)
+GivensFGFTParallelComplex<FPP,DEVICE,FPP2>::GivensFGFTParallelComplex(Faust::MatSparse<FPP,DEVICE>& Lap, int J, int t, unsigned int verbosity, const double stoppingError, const bool errIsRel, const bool enable_large_Faust) : GivensFGFTComplex<FPP,DEVICE,FPP2>(Lap, J, verbosity, stoppingError, errIsRel, enable_large_Faust), t(t), fact_nrots(0)
{
- this->facts.resize(round(J/(float)t));
+ if(J > 0) this->facts.resize(round(J/(float)t));
this->coord_choices.resize(0);
}
diff --git a/wrapper/matlab/+matfaust/+fact/eigtj.m b/wrapper/matlab/+matfaust/+fact/eigtj.m
index 655708662..0b7a95703 100644
--- a/wrapper/matlab/+matfaust/+fact/eigtj.m
+++ b/wrapper/matlab/+matfaust/+fact/eigtj.m
@@ -75,7 +75,7 @@
%==========================================================================================
function [V,D] = eigtj(M, varargin)
import matfaust.Faust
- nGivens_per_fac = 1; % default value
+ nGivens_per_fac = floor(size(M,1)/2); % default value
verbosity = 0; % default value
% if(~ ismatrix(M) || ~ isreal(M))
% error('M must be a real matrix.')
@@ -88,11 +88,18 @@ function [V,D] = eigtj(M, varargin)
tol = 0;
relerr = true;
verbosity = 0;
+ enable_large_Faust = false;
argc = length(varargin);
order = 1; % ascending order
if(argc > 0)
for i=1:argc
switch(varargin{i})
+ case 'enable_large_Faust'
+ if(argc == i || ~ islogical(varargin{i+1}))
+ error('enable_large_Faust keyword argument is not followed by a logical')
+ else
+ enable_large_Faust = varargin{i+1};
+ end
case 'nGivens'
if(argc == i || ~ isnumeric(varargin{i+1}) || varargin{i+1}-floor(varargin{i+1}) > 0 || varargin{i+1} <= 0 )
error('nGivens keyword arg. is not followed by an integer')
@@ -149,7 +156,7 @@ function [V,D] = eigtj(M, varargin)
if(nGivens > 0)
nGivens_per_fac = min(nGivens_per_fac, nGivens);
end
- [core_obj, D] = mexfgftgivensReal(M, nGivens, nGivens_per_fac, verbosity, tol, relerr, order);
+ [core_obj, D] = mexfgftgivensReal(M, nGivens, nGivens_per_fac, verbosity, tol, relerr, order, enable_large_Faust);
D = sparse(diag(real(D)));
V = Faust(core_obj, isreal(M));
end
diff --git a/wrapper/matlab/src/mexfgftgivens.cpp.in b/wrapper/matlab/src/mexfgftgivens.cpp.in
index d2ac822e2..898dab510 100644
--- a/wrapper/matlab/src/mexfgftgivens.cpp.in
+++ b/wrapper/matlab/src/mexfgftgivens.cpp.in
@@ -55,9 +55,9 @@ 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(const mxArray* matlab_matrix, int J, int t, double tol, unsigned int verbosity, bool rel_err, int order, const bool enable_large_Faust, 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 fgft_givens_cplx(const mxArray* matlab_matrix, int J, int t, double tol, unsigned int verbosity, bool rel_err, int order, const bool enable_large_Faust, mxArray **plhs);
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
@@ -66,11 +66,12 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
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;
+ bool rel_err = true;
+ bool enable_large_Faust = false;
int order;
- if(nrhs < 2 || nrhs > 7)
- mexErrMsgTxt("Bad Number of inputs arguments");
+ if(nrhs < 2 || nrhs > 8)
+ mexErrMsgTxt("Bad number of input arguments");
J = (int) mxGetScalar(prhs[1]);
if(nrhs >= 3)
@@ -80,22 +81,24 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
if(nrhs >= 5)
tol = (double) mxGetScalar(prhs[4]);
if(nrhs >= 6)
- relErr = (bool) mxGetScalar(prhs[5]);
+ rel_err = (bool) mxGetScalar(prhs[5]);
if(nrhs >= 7)
order = (int) mxGetScalar(prhs[6]); //eigenvalues order
+ if(nrhs >= 8)
+ enable_large_Faust = (bool) mxGetScalar(prhs[7]);
tol *= tol; // C++ backend works with squared norm error
const mxArray* matlab_matrix = prhs[0]; // Laplacian
if(mxIsComplex(matlab_matrix))
- fgft_givens_cplx(matlab_matrix, J, t, tol, verbosity, relErr, order, plhs);
+ fgft_givens_cplx(matlab_matrix, J, t, tol, verbosity, rel_err, order, enable_large_Faust, plhs);
else
- fgft_givens(matlab_matrix, J, t, tol, verbosity, relErr, order, plhs);
+ fgft_givens(matlab_matrix, J, t, tol, verbosity, rel_err, order, enable_large_Faust, plhs);
}
-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(const mxArray* matlab_matrix, int J, int t, double tol, unsigned int verbosity, bool rel_err, int order, const bool enable_large_Faust, mxArray **plhs)
{
// initialization of the matrix that will be factorized
Faust::MatGeneric<SCALAR,Cpu>* Lap;
@@ -113,18 +116,18 @@ void fgft_givens(const mxArray* matlab_matrix, int J, int t, double tol, unsigne
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 GivensFGFT<SCALAR,Cpu,FPP2>(sLap, J, verbosity, tol, relErr);
+ algo = new GivensFGFT<SCALAR,Cpu,FPP2>(sLap, J, verbosity, tol, rel_err, enable_large_Faust);
else
- algo = new GivensFGFTParallel<SCALAR,Cpu,FPP2>(sLap, J, t, verbosity, tol, relErr);
+ algo = new GivensFGFTParallel<SCALAR,Cpu,FPP2>(sLap, J, t, verbosity, tol, rel_err, enable_large_Faust);
}else
{
mxArray2FaustMat(matlab_matrix,dLap);
Lap = &dLap;
if(t <= 1)
- algo = new GivensFGFT<SCALAR,Cpu,FPP2>(dLap, J, verbosity, tol, relErr);
+ algo = new GivensFGFT<SCALAR,Cpu,FPP2>(dLap, J, verbosity, tol, rel_err, enable_large_Faust);
else
- algo = new GivensFGFTParallel<SCALAR,Cpu,FPP2>(dLap, J, t, verbosity, tol, relErr);
+ algo = new GivensFGFTParallel<SCALAR,Cpu,FPP2>(dLap, J, t, verbosity, tol, rel_err, enable_large_Faust);
}
@@ -142,16 +145,18 @@ void fgft_givens(const mxArray* matlab_matrix, int J, int t, double tol, unsigne
plhs[0] = convertPtr2Mat<Faust::TransformHelper<SCALAR, Cpu>>(th);
- delete algo;
}
catch (const std::exception& e)
{
mexErrMsgTxt(e.what());
}
+ if(algo != nullptr)
+ delete algo;
+
}
-void fgft_givens_cplx(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 rel_err, int order, const bool enable_large_Faust, mxArray **plhs)
{
// initialization of the matrix that will be factorized
Faust::MatGeneric<complex<SCALAR>,Cpu>* Lap;
@@ -160,7 +165,7 @@ void fgft_givens_cplx(const mxArray* matlab_matrix, int J, int t, double tol, un
Faust::BlasHandle<Cpu> blas_handle;
Faust::SpBlasHandle<Cpu> spblas_handle;
- Faust::GivensFGFTComplex<complex<SCALAR>, Cpu, FPP2>* algo;
+ Faust::GivensFGFTComplex<complex<SCALAR>, Cpu, FPP2>* algo = nullptr;
try{
if (mxIsSparse(matlab_matrix))
@@ -169,22 +174,20 @@ void fgft_givens_cplx(const mxArray* matlab_matrix, int J, int t, double tol, un
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);
+ algo = new GivensFGFTComplex<complex<SCALAR>,Cpu,FPP2>(sLap, J, verbosity, tol, rel_err, enable_large_Faust);
else
- algo = new GivensFGFTParallelComplex<complex<SCALAR>,Cpu,FPP2>(sLap, J, t, verbosity, tol, relErr);
+ algo = new GivensFGFTParallelComplex<complex<SCALAR>,Cpu,FPP2>(sLap, J, t, verbosity, tol, rel_err, enable_large_Faust);
}else
{
mxArray2FaustMat(matlab_matrix,dLap);
Lap = &dLap;
if(t <= 1)
- algo = new GivensFGFTComplex<complex<SCALAR>,Cpu,FPP2>(dLap, J, verbosity, tol, relErr);
+ algo = new GivensFGFTComplex<complex<SCALAR>,Cpu,FPP2>(dLap, J, verbosity, tol, rel_err, enable_large_Faust);
else
- algo = new GivensFGFTParallelComplex<complex<SCALAR>,Cpu,FPP2>(dLap, J, t, verbosity, tol, relErr);
+ algo = new GivensFGFTParallelComplex<complex<SCALAR>,Cpu,FPP2>(dLap, J, t, verbosity, tol, rel_err, enable_large_Faust);
}
-
-
algo->compute_facts();
Faust::Vect<SCALAR,Cpu> D(Lap->getNbRow());
@@ -198,13 +201,15 @@ void fgft_givens_cplx(const mxArray* matlab_matrix, int J, int t, double tol, un
plhs[0] = convertPtr2Mat<Faust::TransformHelper<complex<SCALAR>, Cpu>>(th);
- delete algo;
}
catch (const std::exception& e)
{
+ mexErrMsgTxt(e.what());
}
+ if(algo != nullptr)
+ delete algo;
}
diff --git a/wrapper/python/pyfaust/fact.py b/wrapper/python/pyfaust/fact.py
index 2d2f23e5c..72602d333 100644
--- a/wrapper/python/pyfaust/fact.py
+++ b/wrapper/python/pyfaust/fact.py
@@ -170,7 +170,7 @@ def svdtj(M, nGivens=None, tol=0, order='ascend', relerr=True,
return U, S, V
def eigtj(M, nGivens=None, tol=0, order='ascend', relerr=True,
- nGivens_per_fac=None, verbosity=0):
+ nGivens_per_fac=None, verbosity=0, enable_large_Faust=False):
"""
Performs an approximate eigendecomposition of M and returns the eigenvalues in W along with the corresponding left eigenvectors (as the columns of the Faust object V).
@@ -255,12 +255,13 @@ def eigtj(M, nGivens=None, tol=0, order='ascend', relerr=True,
svdtj
"""
D, core_obj = _FaustCorePy.FaustFact.eigtj(M, nGivens, tol, relerr,
- nGivens_per_fac, verbosity, order)
+ nGivens_per_fac, verbosity, order,
+ enable_large_Faust)
return D, Faust(core_obj=core_obj)
# experimental block start
def fgft_givens(Lap, nGivens=None, tol=0, order='ascend', relerr=True,
- nGivens_per_fac=None, verbosity=0):
+ nGivens_per_fac=None, verbosity=0, enable_large_Faust=False):
"""
Computes the FGFT of the Laplacian matrix Lap (using fact.eigtj).
@@ -292,7 +293,8 @@ def fgft_givens(Lap, nGivens=None, tol=0, order='ascend', relerr=True,
See also:
eigtj, fgft_palm
"""
- return eigtj(Lap, nGivens, tol, order, relerr, nGivens_per_fac, verbosity)
+ return eigtj(Lap, nGivens, tol, order, relerr, nGivens_per_fac, verbosity,
+ enable_large_Faust)
# experimental block end
def _check_fact_mat(funcname, M):
diff --git a/wrapper/python/src/FaustCoreCy.pxd b/wrapper/python/src/FaustCoreCy.pxd
index 92dcc7368..8e8ca85bf 100644
--- a/wrapper/python/src/FaustCoreCy.pxd
+++ b/wrapper/python/src/FaustCoreCy.pxd
@@ -188,20 +188,22 @@ cdef extern from "FaustFact.h":
cdef extern from "FaustFactGivensFGFT.h":
cdef FaustCoreCpp[FPP]* fact_givens_fgft[FPP,FPP2](const FPP* Lap, unsigned int num_rows,
- unsigned int num_cols, unsigned int J,
- unsigned int t, FPP* D, unsigned int verbosity,
+ unsigned int num_cols, unsigned int J,
+ unsigned int t, FPP* D, unsigned int verbosity,
const double stoppingError,
const bool errIsRel,
- const int order)
+ const int order,
+ const bool enable_large_Faust)
cdef FaustCoreCpp[FPP]* fact_givens_fgft_sparse[FPP,FPP2](const FPP* data, int* row_ptr,
- int* id_col, int nnz, unsigned int num_rows,
- unsigned int num_cols, unsigned int J,
+ int* id_col, int nnz, unsigned int num_rows,
+ unsigned int num_cols, unsigned int J,
unsigned int t, FPP* D,
unsigned int verbosity,
const double stoppingError,
const bool errIsRel,
- const int order)
+ const int order,
+ const bool enable_large_Faust)
cdef FaustCoreCpp[FPP]* fact_givens_fgft_cplx[FPP,FPP2](const
FPP* Lap, unsigned int num_rows,
@@ -210,17 +212,19 @@ cdef extern from "FaustFactGivensFGFT.h":
FPP2* D, unsigned int verbosity,
const double stoppingError,
const bool errIsRel,
- const int order)
+ const int order,
+ const bool enable_large_Faust)
cdef FaustCoreCpp[FPP]* fact_givens_fgft_sparse_cplx[FPP, FPP2](const FPP* data, int* row_ptr,
- int* id_col, int nnz, unsigned int num_rows,
- unsigned int num_cols, unsigned int J,
- unsigned int t,
- FPP2* D,
- unsigned int verbosity,
- const double stoppingError,
- const bool errIsRel,
- const int order)
+ int* id_col, int nnz, unsigned int num_rows,
+ unsigned int num_cols, unsigned int J,
+ unsigned int t,
+ FPP2* D,
+ unsigned int verbosity,
+ const double stoppingError,
+ const bool errIsRel,
+ const int order,
+ const bool enable_large_Faust)
cdef void svdtj[FPP, FPP2](FaustCoreCpp[FPP]** U, FaustCoreCpp[FPP] **V, FPP* S,
const FPP* M_data, unsigned int num_rows, unsigned int
diff --git a/wrapper/python/src/FaustFactGivensFGFT.h b/wrapper/python/src/FaustFactGivensFGFT.h
index 701f4a87b..d8174ae51 100644
--- a/wrapper/python/src/FaustFactGivensFGFT.h
+++ b/wrapper/python/src/FaustFactGivensFGFT.h
@@ -7,22 +7,22 @@ FaustCoreCpp<FPP>* fact_givens_fgft(const FPP* Lap, unsigned int num_rows, unsig
template<typename FPP, typename FPP2 = float>
FaustCoreCpp<FPP>* fact_givens_fgft_sparse(const FPP* data, int* row_ptr, int* id_col, int nnz, int nrows, int ncols,
- unsigned int J, unsigned int t /* end of input parameters*/, FPP* D, unsigned int verbosity = 0, const double stoppingError = 0.0, const bool errIsRel = true, const int order = 1 /* ascendant order for eigenvalues */);
+ unsigned int J, unsigned int t /* end of input parameters*/, FPP* D, unsigned int verbosity = 0, const double stoppingError = 0.0, const bool errIsRel = true, const int order = 1 /* ascendant order for eigenvalues */, const bool enable_large_Faust=false);
template<typename FPP, typename FPP2 = float>
-FaustCoreCpp<FPP>* fact_givens_fgft_generic(GivensFGFT<FPP, Cpu, FPP2>* algo, FPP* D, const int order);
+FaustCoreCpp<FPP>* fact_givens_fgft_generic(GivensFGFT<FPP, Cpu, FPP2>* algo, FPP* D, const int order, const bool enable_large_Faust=false);
template<typename FPP, typename FPP2 = float>
-FaustCoreCpp<FPP>* fact_givens_fgft_cplx(const FPP* Lap, unsigned int num_rows, unsigned int num_cols, unsigned int J, unsigned int t, FPP2* D, unsigned int verbosity = 0, const double stoppingError = 0.0, const bool errIsRel = true, const int order = 1 /* ascendant order for eigenvalues */);
+FaustCoreCpp<FPP>* fact_givens_fgft_cplx(const FPP* Lap, unsigned int num_rows, unsigned int num_cols, unsigned int J, unsigned int t, FPP2* D, unsigned int verbosity = 0, const double stoppingError = 0.0, const bool errIsRel = true, const int order = 1 /* ascendant order for eigenvalues */, const bool enable_large_Faust=false);
template<typename FPP, typename FPP2 = float>
FaustCoreCpp<FPP>* fact_givens_fgft_sparse_cplx(const FPP* data, int* row_ptr, int* id_col, int nnz, int nrows, int ncols,
- unsigned int J, unsigned int t, FPP2* D, unsigned int verbosity = 0, const double stoppingError = 0.0, const bool errIsRel = true, const int order = 1 /* ascendant order for eigenvalues */);
+ unsigned int J, unsigned int t, FPP2* D, unsigned int verbosity = 0, const double stoppingError = 0.0, const bool errIsRel = true, const int order = 1 /* ascendant order for eigenvalues */, const bool enable_large_Faust=false);
template<typename FPP, typename FPP2 = float>
-FaustCoreCpp<FPP>* fact_givens_fgft_generic_cplx(GivensFGFTComplex<FPP, Cpu, FPP2>* algo, FPP2* D, const int order);
+FaustCoreCpp<FPP>* fact_givens_fgft_generic_cplx(GivensFGFTComplex<FPP, Cpu, FPP2>* algo, FPP2* D, const int order, const bool enable_large_Faust=false);
template<typename FPP, typename FPP2 = float>
void svdtj( FaustCoreCpp<FPP>** U, FaustCoreCpp<FPP> **V, FPP* S /** end of output parameters */, const FPP* M, unsigned int num_rows, unsigned int num_cols, unsigned int J, unsigned int t, unsigned int verbosity = 0, const double stoppingError = 0.0, const bool errIsRel = true /* end of input parameters*/);
diff --git a/wrapper/python/src/FaustFactGivensFGFT.hpp b/wrapper/python/src/FaustFactGivensFGFT.hpp
index 10dc2eb8f..8a358c680 100644
--- a/wrapper/python/src/FaustFactGivensFGFT.hpp
+++ b/wrapper/python/src/FaustFactGivensFGFT.hpp
@@ -8,23 +8,23 @@ using namespace Faust;
template<typename FPP, typename FPP2>
FaustCoreCpp<FPP>* fact_givens_fgft_sparse(FPP* data, int* row_ptr, int* id_col, int nnz, int nrows, int ncols,
- unsigned int J, unsigned int t /* end of input parameters*/, FPP* D, unsigned int verbosity, const double stoppingError, const bool errIsRel, const int order)
+ unsigned int J, unsigned int t /* end of input parameters*/, FPP* D, unsigned int verbosity, const double stoppingError, const bool errIsRel, const int order, const bool enable_large_Faust)
{
Faust::MatSparse<FPP, Cpu> mat_Lap(nnz, nrows, ncols, data, id_col, row_ptr);
GivensFGFT<FPP, Cpu, FPP2>* algo;
if(t <= 1)
{
- algo = new GivensFGFT<FPP, Cpu, FPP2>(mat_Lap, (int)J, verbosity, stoppingError, errIsRel);
+ algo = new GivensFGFT<FPP, Cpu, FPP2>(mat_Lap, (int)J, verbosity, stoppingError, errIsRel, enable_large_Faust);
}
else
{
- algo = new GivensFGFTParallel<FPP, Cpu, FPP2>(mat_Lap, (int)J, (int) t, verbosity, stoppingError, errIsRel);
+ algo = new GivensFGFTParallel<FPP, Cpu, FPP2>(mat_Lap, (int)J, (int) t, verbosity, stoppingError, errIsRel, enable_large_Faust);
}
return fact_givens_fgft_generic(algo, D, order);
}
template<typename FPP, typename FPP2>
-FaustCoreCpp<FPP>* fact_givens_fgft(const FPP* Lap, unsigned int num_rows, unsigned int num_cols, unsigned int J, unsigned int t /* end of input parameters*/, FPP* D, unsigned int verbosity, const double stoppingError, const bool errIsRel, const int order)
+FaustCoreCpp<FPP>* fact_givens_fgft(const FPP* Lap, unsigned int num_rows, unsigned int num_cols, unsigned int J, unsigned int t /* end of input parameters*/, FPP* D, unsigned int verbosity, const double stoppingError, const bool errIsRel, const int order, const bool enable_large_Faust)
{
//TODO: optimization possible here by avoiding Lap copy in MatDense (by
//just using the data in Lap as underlying pointer of MatDense)
@@ -33,53 +33,60 @@ FaustCoreCpp<FPP>* fact_givens_fgft(const FPP* Lap, unsigned int num_rows, unsig
GivensFGFT<FPP, Cpu, FPP2>* algo;
if(t <= 1)
{
- algo = new GivensFGFT<FPP, Cpu, FPP2>(mat_Lap, (int)J, verbosity, stoppingError, errIsRel);
+ algo = new GivensFGFT<FPP, Cpu, FPP2>(mat_Lap, (int)J, verbosity, stoppingError, errIsRel, enable_large_Faust);
}
else
{
- algo = new GivensFGFTParallel<FPP, Cpu, FPP2>(mat_Lap, (int)J, (int) t, verbosity, stoppingError, errIsRel);
+ algo = new GivensFGFTParallel<FPP, Cpu, FPP2>(mat_Lap, (int)J, (int) t, verbosity, stoppingError, errIsRel, enable_large_Faust);
}
return fact_givens_fgft_generic(algo, D, order);
}
template<typename FPP, typename FPP2>
-FaustCoreCpp<FPP>* fact_givens_fgft_generic(GivensFGFT<FPP, Cpu, FPP2>* algo, FPP* D, const int order)
+FaustCoreCpp<FPP>* fact_givens_fgft_generic(GivensFGFT<FPP, Cpu, FPP2>* algo, FPP* D, const int order, const bool enable_large_Faust)
{
+ FaustCoreCpp<FPP>* fc = nullptr;
algo->compute_facts();
- Faust::Transform<FPP,Cpu> trans = std::move(algo->get_transform(order));
- TransformHelper<FPP,Cpu> *th = new TransformHelper<FPP,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)
+ try {
+ Faust::Transform<FPP,Cpu> trans = std::move(algo->get_transform(order));
+ TransformHelper<FPP,Cpu> *th = new TransformHelper<FPP,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)
- //copy ordered diagonal in buffer D (allocated from the outside)
- algo->get_D(D, order);
+ //copy ordered diagonal in buffer D (allocated from the outside)
+ algo->get_D(D, order);
- delete algo;
+ return new FaustCoreCpp<FPP>(th);
+ }
+ catch(out_of_range e)
+ {
- return new FaustCoreCpp<FPP>(th);
+ }
+ delete algo;
+ return fc;
}
template<typename FPP, typename FPP2>
FaustCoreCpp<FPP>* fact_givens_fgft_sparse_cplx(FPP* data, int* row_ptr, int* id_col, int nnz, int nrows, int ncols,
- unsigned int J, unsigned int t /* end of input parameters*/, FPP2* D, unsigned int verbosity, const double stoppingError, const bool errIsRel, const int order)
+ unsigned int J, unsigned int t /* end of input parameters*/, FPP2* D, unsigned int verbosity, const double stoppingError, const bool errIsRel, const int order, const bool enable_large_Faust)
{
Faust::MatSparse<FPP, Cpu> mat_Lap(nnz, nrows, ncols, data, id_col, row_ptr);
GivensFGFTComplex<FPP, Cpu, FPP2>* algo;
if(t <= 1)
{
- algo = new GivensFGFTComplex<FPP, Cpu, FPP2>(mat_Lap, (int)J, verbosity, stoppingError, errIsRel);
+ algo = new GivensFGFTComplex<FPP, Cpu, FPP2>(mat_Lap, (int)J, verbosity, stoppingError, errIsRel, enable_large_Faust);
}
else
{
- algo = new GivensFGFTParallelComplex<FPP, Cpu, FPP2>(mat_Lap, (int)J, (int) t, verbosity, stoppingError, errIsRel);
+ algo = new GivensFGFTParallelComplex<FPP, Cpu, FPP2>(mat_Lap, (int)J, (int) t, verbosity, stoppingError, errIsRel, enable_large_Faust);
}
return fact_givens_fgft_generic_cplx(algo, D, order);
}
template<typename FPP, typename FPP2>
-FaustCoreCpp<FPP>* fact_givens_fgft_cplx(const FPP* Lap, unsigned int num_rows, unsigned int num_cols, unsigned int J, unsigned int t /* end of input parameters*/, FPP2* D, unsigned int verbosity, const double stoppingError, const bool errIsRel, const int order)
+FaustCoreCpp<FPP>* fact_givens_fgft_cplx(const FPP* Lap, unsigned int num_rows, unsigned int num_cols, unsigned int J, unsigned int t /* end of input parameters*/, FPP2* D, unsigned int verbosity, const double stoppingError, const bool errIsRel, const int order, const bool enable_large_Faust)
{
//TODO: optimization possible here by avoiding Lap copy in MatDense (by
//just using the data in Lap as underlying pointer of MatDense)
@@ -88,32 +95,38 @@ FaustCoreCpp<FPP>* fact_givens_fgft_cplx(const FPP* Lap, unsigned int num_rows,
GivensFGFTComplex<FPP, Cpu, FPP2>* algo;
if(t <= 1)
{
- algo = new GivensFGFTComplex<FPP, Cpu, FPP2>(mat_Lap, (int)J, verbosity, stoppingError, errIsRel);
+ algo = new GivensFGFTComplex<FPP, Cpu, FPP2>(mat_Lap, (int)J, verbosity, stoppingError, errIsRel, enable_large_Faust);
}
else
{
- algo = new GivensFGFTParallelComplex<FPP, Cpu, FPP2>(mat_Lap, (int)J, (int) t, verbosity, stoppingError, errIsRel);
+ algo = new GivensFGFTParallelComplex<FPP, Cpu, FPP2>(mat_Lap, (int)J, (int) t, verbosity, stoppingError, errIsRel, enable_large_Faust);
}
return fact_givens_fgft_generic_cplx(algo, D, order);
}
template<typename FPP, typename FPP2>
-FaustCoreCpp<FPP>* fact_givens_fgft_generic_cplx(GivensFGFTComplex<FPP, Cpu, FPP2>* algo, FPP2* D, const int order)
+FaustCoreCpp<FPP>* fact_givens_fgft_generic_cplx(GivensFGFTComplex<FPP, Cpu, FPP2>* algo, FPP2* D, const int order, const bool enable_large_Faust)
{
-
+ FaustCoreCpp<FPP>* fc = nullptr;
algo->compute_facts();
- Faust::Transform<FPP,Cpu> trans = std::move(algo->get_transform(order));
- TransformHelper<FPP,Cpu> *th = new TransformHelper<FPP,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)
-
+ try
+ {
+ Faust::Transform<FPP,Cpu> trans = std::move(algo->get_transform(order));
+ TransformHelper<FPP,Cpu> *th = new TransformHelper<FPP,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)
- //copy ordered diagonal in buffer D (allocated from the outside)
- algo->get_D(D, order);
+ //copy ordered diagonal in buffer D (allocated from the outside)
+ algo->get_D(D, order);
+ fc =new FaustCoreCpp<FPP>(th);
+ }
+ catch(out_of_range e)
+ {
+ }
delete algo;
- return new FaustCoreCpp<FPP>(th);
+ return fc;
}
template<typename FPP, typename FPP2>
diff --git a/wrapper/python/src/_FaustCorePy.pyx b/wrapper/python/src/_FaustCorePy.pyx
index 4e4e3667d..41112d894 100644
--- a/wrapper/python/src/_FaustCorePy.pyx
+++ b/wrapper/python/src/_FaustCorePy.pyx
@@ -1435,7 +1435,7 @@ cdef class FaustFact:
@staticmethod
def fact_givens_fgft_sparse(Lap, J, t, verbosity=0, stoppingError = 0.0,
- errIsRel=True, order='ascend'):
+ errIsRel=True, order='ascend', enable_large_Faust=False):
from scipy.sparse import spdiags
cdef double [:] data1d #only for csr mat factor
cdef int [:] indices # only for csr mat
@@ -1467,16 +1467,23 @@ cdef class FaustFact:
verbosity,
stoppingError,
errIsRel,
- int(order))
+ int(order),
+ enable_large_Faust)
core._isReal = True
#D_spdiag = spdiags(D, [0], Lap.shape[0], Lap.shape[0])
#return core, D_spdiag
+ if(core._isReal and core.core_faust_dbl == NULL or \
+ not core._isReal and core.core_faust_cplx == NULL):
+ raise Exception("Empty transform (nGivens is too big ? Set"
+ " enable_large_Faust to True to force the computation).")
+
+
return core, D
@staticmethod
def fact_givens_fgft(Lap, J, t, verbosity=0, stoppingError = 0.0,
- errIsRel=True, order=1):
+ errIsRel=True, order=1, enable_large_Faust=False):
isReal = Lap.dtype in [ 'float', 'float128',
'float16', 'float32',
'float64', 'double']
@@ -1506,17 +1513,24 @@ cdef class FaustFact:
&D_view[0], verbosity,
stoppingError,
errIsRel,
- int(order))
+ int(order),
+ enable_large_Faust)
core._isReal = True
#from scipy.sparse import spdiags
#D_spdiag = spdiags(D, [0], Lap.shape[0], Lap.shape[0])
#return core, D_spdiag
+ if(core._isReal and core.core_faust_dbl == NULL or \
+ not core._isReal and core.core_faust_cplx == NULL):
+ raise Exception("Empty transform (nGivens is too big ? Set"
+ " enable_large_Faust to True to force the computation).")
+
+
return core, D
@staticmethod
def fact_givens_fgft_cplx(Lap, J, t, verbosity=0, stoppingError = 0.0,
- errIsRel=True, order=1):
+ errIsRel=True, order=1, enable_large_Faust=False):
isReal = Lap.dtype in [ 'float', 'float128',
'float16', 'float32',
'float64', 'double']
@@ -1541,22 +1555,29 @@ cdef class FaustFact:
core = FaustCore(core=True)
core.core_faust_cplx = FaustCoreCy.fact_givens_fgft_cplx[complex,double](&Lap_view[0,0],
- Lap_num_rows,
- Lap_num_cols, J, t,
- &D_view[0], verbosity,
- stoppingError,
- errIsRel,
- int(order))
+ Lap_num_rows,
+ Lap_num_cols, J, t,
+ &D_view[0], verbosity,
+ stoppingError,
+ errIsRel,
+ int(order),
+ enable_large_Faust)
core._isReal = False
#from scipy.sparse import spdiags
#D_spdiag = spdiags(D, [0], Lap.shape[0], Lap.shape[0])
#return core, D_spdiag
+ if(core._isReal and core.core_faust_dbl == NULL or \
+ not core._isReal and core.core_faust_cplx == NULL):
+ raise Exception("Empty transform (nGivens is too big ? Set"
+ " enable_large_Faust to True to force the computation).")
+
return core, D
@staticmethod
def fact_givens_fgft_sparse_cplx(Lap, J, t, verbosity=0, stoppingError = 0.0,
- errIsRel=True, order='ascend'):
+ errIsRel=True, order='ascend',
+ enable_large_Faust=False):
from scipy.sparse import spdiags
cdef complex[:] data1d #only for csr mat factor
cdef int [:] indices # only for csr mat
@@ -1577,39 +1598,46 @@ cdef class FaustFact:
D_view = D
core = FaustCore(core=True)
- core.core_faust_cplx = FaustCoreCy.fact_givens_fgft_sparse_cplx[complex, double](&data1d[0],
- &indices[0],
- &indptr[0],
- Lap.nnz,
- Lap_num_rows,
- Lap_num_cols, J, t,
- &D_view[0],
- verbosity,
- stoppingError,
- errIsRel,
- int(order))
+ core.core_faust_cplx = \
+ FaustCoreCy.fact_givens_fgft_sparse_cplx[complex, double](&data1d[0],
+ &indices[0],
+ &indptr[0],
+ Lap.nnz,
+ Lap_num_rows,
+ Lap_num_cols, J, t,
+ &D_view[0],
+ verbosity,
+ stoppingError,
+ errIsRel,
+ int(order),
+ enable_large_Faust)
core._isReal = False
#D_spdiag = spdiags(D, [0], Lap.shape[0], Lap.shape[0])
#return core, D_spdiag
+ if(core._isReal and core.core_faust_dbl == NULL or \
+ not core._isReal and core.core_faust_cplx == NULL):
+ raise Exception("Empty transform (nGivens is too big ? Set"
+ " enable_large_Faust to True to force the computation).")
+
return core, D
@staticmethod
- def eigtj(M, maxiter=None, tol=0, relerr=True, nGivens_per_fac=None, verbosity=0,
- order='ascend'):
- if(maxiter == None):
+ def eigtj(M, nGivens=None, tol=0, relerr=True, nGivens_per_fac=None, verbosity=0,
+ order='ascend', enable_large_Faust=False):
+ if(nGivens == None):
if(tol == 0):
- raise Exception("You must specify maxiter or tol argument"
+ raise Exception("You must specify nGivens or tol argument"
" (to define a stopping criterion)")
- maxiter = 0
+ nGivens = 0
if(nGivens_per_fac == None): nGivens_per_fac = int(M.shape[0]/2)
if(isinstance(M, np.ndarray) and \
not np.allclose(np.matrix(M, copy=False).H, M) or M.shape[0] != M.shape[1]):
raise ValueError(" the matrix/array must be symmetric or hermitian.")
- if(not isinstance(maxiter, int)): raise TypeError("maxiter must be a int")
+ if(not isinstance(nGivens, int)): raise TypeError("nGivens must be a int")
if(not isinstance(nGivens_per_fac, int)): raise TypeError("nGivens_per_fac must be a int")
nGivens_per_fac = max(nGivens_per_fac, 1)
- if(maxiter > 0): nGivens_per_fac = min(nGivens_per_fac, maxiter)
+ if(nGivens > 0): nGivens_per_fac = min(nGivens_per_fac, nGivens)
tol *= tol # the C++ impl. works on squared norms to measure errors
M_is_real = M.dtype in [ 'float', 'float128',
@@ -1618,26 +1646,26 @@ cdef class FaustFact:
if(isinstance(M, np.ndarray)):
M = np.asfortranarray(M)
if(M_is_real):
- core_obj,D = FaustFact.fact_givens_fgft(M, maxiter, nGivens_per_fac,
+ core_obj,D = FaustFact.fact_givens_fgft(M, nGivens, nGivens_per_fac,
verbosity, tol,
- relerr, order)
+ relerr, order, enable_large_Faust)
else: #complex
- core_obj,D = FaustFact.fact_givens_fgft_cplx(M, maxiter, nGivens_per_fac,
+ core_obj,D = FaustFact.fact_givens_fgft_cplx(M, nGivens, nGivens_per_fac,
verbosity, tol,
- relerr, order)
+ relerr, order, enable_large_Faust)
elif(isinstance(M, csr_matrix)):
if(M_is_real):
- core_obj,D = FaustFact.fact_givens_fgft_sparse(M, maxiter,
+ core_obj,D = FaustFact.fact_givens_fgft_sparse(M, nGivens,
nGivens_per_fac,
verbosity,
tol,
relerr,
- order)
+ order, enable_large_Faust)
else: #complex
- core_obj,D = FaustFact.fact_givens_fgft_sparse_cplx(M, maxiter, nGivens_per_fac,
+ core_obj,D = FaustFact.fact_givens_fgft_sparse_cplx(M, nGivens, nGivens_per_fac,
verbosity, tol,
- relerr, order)
+ relerr, order, enable_large_Faust)
else:
raise TypeError("The matrix to diagonalize must be a"
" scipy.sparse.csr_matrix or a numpy array.")
--
GitLab