Mentions légales du service

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

Update basisChebyshev and co. (don't skip the degree-0 polynomial anymore and... 

Update basisChebyshev and co. (don't skip the degree-0 polynomial anymore and modify all related functions in consequence).
parent 3e79d1b2
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/faust_linear_operator/CPU/faust_TransformHelperPoly.hpp
+ 67
− 71
View file @ c58b4085
...@@ -14,7 +14,7 @@ namespace Faust ...@@ -14,7 +14,7 @@ namespace Faust
Vect<FPP,Cpu> TransformHelperPoly<FPP>::multiply(const FPP* x, const bool transpose/*=false*/, const bool conjugate/*=false*/) Vect<FPP,Cpu> TransformHelperPoly<FPP>::multiply(const FPP* x, const bool transpose/*=false*/, const bool conjugate/*=false*/)
{ {
int d = L->getNbRow(); int d = L->getNbRow();
int K = this->size(); int K = this->size()-1;
Vect<FPP, Cpu> v0(d*(K+1)); Vect<FPP, Cpu> v0(d*(K+1));
multiply(x, v0.getData(), transpose, conjugate); multiply(x, v0.getData(), transpose, conjugate);
return std::move(v0); return std::move(v0);
...@@ -38,23 +38,22 @@ namespace Faust ...@@ -38,23 +38,22 @@ namespace Faust
* Fx = [v_{0, K} ; v_{1, K} ; v_{2, K}] * Fx = [v_{0, K} ; v_{1, K} ; v_{2, K}]
*/ */
int d = L->getNbRow(); int d = L->getNbRow();
int K = this->size(); int K = this->size()-1;
// seeing how basisChebyshev method is instantiating this class
// K can't be strictly lower than two
Eigen::Map<Eigen::Matrix<FPP, Eigen::Dynamic, 1>> x_vec(const_cast<FPP*>(x), d);
memcpy(y, x, sizeof(FPP)*d); memcpy(y, x, sizeof(FPP)*d);
if(K == 0)
return;
Eigen::Map<Eigen::Matrix<FPP, Eigen::Dynamic, 1>> x_vec(const_cast<FPP*>(x), d);
Eigen::Map<Eigen::Matrix<FPP, Eigen::Dynamic, 1>>v2(const_cast<FPP*>(y+d), d); Eigen::Map<Eigen::Matrix<FPP, Eigen::Dynamic, 1>>v2(const_cast<FPP*>(y+d), d);
v2 = L->mat*x_vec; v2 = L->mat*x_vec;
if(K == 1) // not necessary but clearer
return;
for(int i=3;i<=K+1;i++) for(int i=3;i<=K+1;i++)
{ {
int offset = d*(i-3); Eigen::Map<Eigen::Matrix<FPP, Eigen::Dynamic, 1>>new_v2_(const_cast<FPP*>(y+d*(i-1)), d);
// gemv(L, v2, tmp, FPP(1), FPP(-1));
// v2 = twoL*v2;
Eigen::Map<Eigen::Matrix<FPP, Eigen::Dynamic, 1>>v2_(const_cast<FPP*>(y+d*(i-2)), d); Eigen::Map<Eigen::Matrix<FPP, Eigen::Dynamic, 1>>v2_(const_cast<FPP*>(y+d*(i-2)), d);
Eigen::Map<Eigen::Matrix<FPP, Eigen::Dynamic, 1>>v2__(const_cast<FPP*>(y+d*(i-1)), d); Eigen::Map<Eigen::Matrix<FPP, Eigen::Dynamic, 1>>v1_(const_cast<FPP*>(y+d*(i-3)), d);
Eigen::Map<Eigen::Matrix<FPP, Eigen::Dynamic, 1>>y_vec(const_cast<FPP*>(y+offset), d); new_v2_ = L->mat*v2_*2-v1_;
v2__ = L->mat*v2_*2-y_vec;
} }
} }
...@@ -99,7 +98,7 @@ namespace Faust ...@@ -99,7 +98,7 @@ namespace Faust
{ {
// std::cout << "TransformHelperPoly<FPP>::multiply(MatDense)" << std::endl; // std::cout << "TransformHelperPoly<FPP>::multiply(MatDense)" << std::endl;
int d = L->getNbRow(); int d = L->getNbRow();
int K = this->size(); int K = this->size()-1;
int n = X.getNbCol(); int n = X.getNbCol();
MatDense<FPP,Cpu> V0_ord(d*(K+1), n); MatDense<FPP,Cpu> V0_ord(d*(K+1), n);
auto scale = (K+1)*d; auto scale = (K+1)*d;
...@@ -132,27 +131,26 @@ namespace Faust ...@@ -132,27 +131,26 @@ namespace Faust
template<typename FPP> template<typename FPP>
void TransformHelperPoly<FPP>::poly(int d, int n, const FPP* basisX, const FPP* coeffs, FPP* out) void TransformHelperPoly<FPP>::poly(int d, int n, const FPP* basisX, const FPP* coeffs, FPP* out)
{ {
int K = this->size(); int K = this->size()-1;
Faust::poly(d, K, n, basisX, coeffs, out); Faust::poly(d, K, n, basisX, coeffs, out);
} }
template<typename FPP> template<typename FPP>
TransformHelper<FPP, Cpu>* TransformHelperPoly<FPP>::next() TransformHelper<FPP, Cpu>* TransformHelperPoly<FPP>::next()
{ {
int K = this->size(); int K = this->size()-1;
return next(K+1); return next(K+1);
} }
template<typename FPP> template<typename FPP>
TransformHelper<FPP, Cpu>* TransformHelperPoly<FPP>::next(int K) TransformHelper<FPP, Cpu>* TransformHelperPoly<FPP>::next(int K)
{ {
auto old_K = this->size(); auto old_K = this->size()-1;
auto d = L->getNbRow(); auto d = L->getNbRow();
MatSparse<FPP,Cpu> Id, zero, R; MatSparse<FPP,Cpu> Id, zero, R;
std::vector<MatGeneric<FPP,Cpu>*> facts(K); std::vector<MatGeneric<FPP,Cpu>*> facts(K+1);
facts.resize(K); for(int i=0;i<old_K+1;i++)
for(int i=0;i<old_K;i++) facts[K-i] = this->get_gen_fact_nonconst(old_K-i);
facts[K-i-1] = this->get_gen_fact_nonconst(old_K-1-i);
#pragma omp parallel for private(Id, zero, R) #pragma omp parallel for private(Id, zero, R)
for(int i=old_K+1;i <= K; i++) for(int i=old_K+1;i <= K; i++)
{ {
...@@ -169,7 +167,6 @@ namespace Faust ...@@ -169,7 +167,6 @@ namespace Faust
/* optimizedCopy */ false, /* optimizedCopy */ false,
/* cloning_fact */ false, /* cloning_fact */ false,
/* internal_call */ true); /* internal_call */ true);
//TODO: use reference manager to avoid copies
basisP->L = L; basisP->L = L;
ref_man.acquire(L); ref_man.acquire(L);
basisP->rR = rR; basisP->rR = rR;
...@@ -186,12 +183,11 @@ namespace Faust ...@@ -186,12 +183,11 @@ namespace Faust
MatSparse<FPP,Cpu> Id, Id1; MatSparse<FPP,Cpu> Id, Id1;
MatSparse<FPP,Cpu> coeffDiags; MatSparse<FPP,Cpu> coeffDiags;
auto d = L->getNbRow(); auto d = L->getNbRow();
int K = this->size(); int K = this->size()-1;
Id.resize(d, d, d); Id.resize(d, d, d);
Id.setEyes(); Id.setEyes();
auto Id0 = Id; auto Id0 = Id;
Id0 *= coeffs[0]; Id0 *= coeffs[0];
// K >= 2
for(int i=1;i < K+1; i++) for(int i=1;i < K+1; i++)
{ {
Id1 = Id; Id1 = Id;
...@@ -200,10 +196,10 @@ namespace Faust ...@@ -200,10 +196,10 @@ namespace Faust
coeffDiags.hstack(Id0, Id1); coeffDiags.hstack(Id0, Id1);
Id0 = coeffDiags; Id0 = coeffDiags;
} }
std::vector<MatGeneric<FPP,Cpu>*> facts(K+1); std::vector<MatGeneric<FPP,Cpu>*> facts(this->size()+1);
facts.resize(K+1); facts.resize(this->size()+1);
facts[0] = new MatSparse<FPP,Cpu>(Id0); facts[0] = new MatSparse<FPP,Cpu>(Id0);
for(int i=1;i<=K;i++) for(int i=1;i<=K+1;i++)
{ {
facts[i] = this->get_gen_fact_nonconst(i-1); facts[i] = this->get_gen_fact_nonconst(i-1);
} }
...@@ -239,58 +235,58 @@ namespace Faust ...@@ -239,58 +235,58 @@ namespace Faust
MatSparse<FPP,Cpu> *T1, *T2; MatSparse<FPP,Cpu> *T1, *T2;
auto d = L->getNbRow(); auto d = L->getNbRow();
MatSparse<FPP,Cpu> Id, twoL, minus_Id, *rR, R, zero; MatSparse<FPP,Cpu> Id, twoL, minus_Id, *rR, R, zero;
std::vector<MatGeneric<FPP,Cpu>*> facts(K); //normally it'd be K+1 but T0 is ignored std::vector<MatGeneric<FPP,Cpu>*> facts(K+1); //normally it'd be K+1 but T0 is ignored except if K==0
if(K == 0)
return TransformHelper<FPP,Cpu>::eyeFaust(d,d);
// build the chebyshev polynomials by factor
// K > 1 ignore the first one (for T0 0-degree), which is the identity
// Identity // Identity
Id.resize(d, d, d); Id.resize(d, d, d);
Id.setEyes(); Id.setEyes();
// T1 facts[K] = new MatSparse<FPP,Cpu>(Id);
T1 = new MatSparse<FPP,Cpu>(); if (K > 0)
T1->vstack(Id, *L);
facts[K-1] = T1;
if(K == 1)
return new TransformHelper<FPP,Cpu>(facts, (FPP)1.0,
/* optimizedCopy */ false,
/* cloning_fact */ false,
/* internal_call */ true);
// rR block
twoL = *L;
twoL *= FPP(2);
minus_Id.resize(d,d,d);
minus_Id.setEyes();
minus_Id *= FPP(-1);
rR = new MatSparse<FPP,Cpu>();
rR->hstack(minus_Id, twoL);
// T2
Id.resize(2*d, 2*d, 2*d);
Id.setEyes();
T2 = new MatSparse<FPP,Cpu>();
T2->vstack(Id, *rR);
facts[K-2] = T2;
// T3 to TK
//#pragma omp parallel for private(Id, zero, R)
for(int i=3;i<K+1;i++)
{ {
auto id = i*d; // build the chebyshev polynomials by factor
Id.resize(id, id, id); // K >= 1 ignore the first one (for T0 0-degree), which is the identity
Id.setEyes(); // T1
zero.conservativeResize(d, (i-2)*d); T1 = new MatSparse<FPP,Cpu>();
if(i == 3) T1->vstack(Id, *L);
{ facts[K-1] = T1;
R.hstack(zero, *rR); // rR block
} twoL = *L;
else twoL *= FPP(2);
minus_Id.resize(d,d,d);
minus_Id.setEyes();
minus_Id *= FPP(-1);
rR = new MatSparse<FPP,Cpu>();
rR->hstack(minus_Id, twoL);
if (K > 1)
{ {
R.conservativeResize(d, rR->getNbCol()+(i-2)*d); // T2
for(int i=0;i<R.getNonZeros();i++) Id.resize(2*d, 2*d, 2*d);
R.getColInd()[i] += d; Id.setEyes();
T2 = new MatSparse<FPP,Cpu>();
T2->vstack(Id, *rR);
facts[K-2] = T2;
// T3 to TK
//#pragma omp parallel for private(Id, zero, R)
for(int i=3;i<K+1;i++)
{
auto id = i*d;
Id.resize(id, id, id);
Id.setEyes();
zero.conservativeResize(d, (i-2)*d);
if(i == 3)
{
R.hstack(zero, *rR);
}
else
{
R.conservativeResize(d, rR->getNbCol()+(i-2)*d);
for(int i=0;i<R.getNonZeros();i++)
R.getColInd()[i] += d;
}
auto Ti = new MatSparse<FPP,Cpu>();
Ti->vstack(Id, R);
facts[K-i] = Ti;
}
} }
auto Ti = new MatSparse<FPP,Cpu>();
Ti->vstack(Id, R);
facts[K-i] = Ti;
} }
auto basisP = new TransformHelperPoly<FPP>(facts, (FPP)1.0, auto basisP = new TransformHelperPoly<FPP>(facts, (FPP)1.0,
/* optimizedCopy */ false, /* optimizedCopy */ false,
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment