Update TransformHelper::pack_factors() by using cleaner methods: insert,...

Update TransformHelper::pack_factors() by using cleaner methods: insert, erase, resize, pop_front, pop_back, push_first (defined in Faust::Transform). New signature for get_product().

Other minor changes: debug functions Transform(Helper)::disable_dtor(), update related test on pack_factors...

misc/test/src/C++/unit/test_TH_pack_factors.cpp.in

@@ -88,7 +88,7 @@ int main(int argc, char* argv[])
 		facts[i] = mat;
 	}
 
-	TransformHelper<FPP,Cpu> th0(facts, 1.0, false, false, false); // a Faust that uses facts without copy (I want to be sure factors are not altered on the end)
+	TransformHelper<FPP,Cpu> th0(facts, 1.0, false, true, false); 
 
 	/************** Test packing to the right */
 	cout << "* testing packing factors on the right (from id 3 included)" << endl;
@@ -125,7 +125,7 @@ int main(int argc, char* argv[])
 		assert((*(th2.begin()+i))->norm() == (*(facts.begin()+i))->norm());
 	assert(th2.size() == facts.size());
 
-	// pack factors on the right of factor index 2
+	// pack factors on the left of factor index 4
 	th2.pack_factors(3, PackDir::PACK_LEFT);
 
 	//test it's done properly
@@ -153,7 +153,11 @@ int main(int argc, char* argv[])
 	th3.pack_factors(0,0);
 	th3.pack_factors(4,4);
 	for(int i=0;i<th3.size();i++)
-		assert((*(th3.begin()+i))->norm() == (*(facts.begin()+i))->norm());
+    {
+//        cout << "norm new fact:" << (*(th3.begin()+i))->norm() << " norm old fact:" << (*(facts.begin()+i))->norm() << endl;
+        Real<FPP> dn = Faust::fabs((*(th3.begin()+i))->norm() - (*(facts.begin()+i))->norm());
+		assert(dn <= .0001);
+    }
 	assert(th3.size() == facts.size());
 	cout << "OK" << endl;
 	/************************* test packing all */

src/algorithm/factorization/faust_Params.h

@@ -158,7 +158,7 @@ namespace Faust
         public:
         // Required members
         
-	// data is now independant from the params class,
+	// data is now independent from the params class,
         //Faust::MatDense<FPP,DEVICE> data; 
 	faust_unsigned_int m_nbRow; // number of row of the matrix to be factorized
 	faust_unsigned_int m_nbCol; // number of columns of the matrix to be factorized

src/faust_linear_operator/CPU/faust_Transform.h

@@ -134,7 +134,8 @@ namespace Faust
 				//template<Device DEVICE> class SpBlasHandle;
 				/** \brief Perform the product of all factorized matrix. */
 				Faust::MatDense<FPP,Cpu> get_product(const char opThis='N', const bool isConj=false)const;
-				Faust::MatDense<FPP,Cpu> get_product(Faust::BlasHandle<Cpu> blas_handle,Faust::SpBlasHandle<Cpu> spblas_handle)const;
+				void get_product(Faust::MatDense<FPP,Cpu> &, const char opThis='N', const bool isConj=false)const;
+                Faust::MatDense<FPP,Cpu> get_product(Faust::BlasHandle<Cpu> blas_handle,Faust::SpBlasHandle<Cpu> spblas_handle)const;
 
 
 				/** \brief return a copy of the factor of index id
@@ -208,7 +209,12 @@ namespace Faust
 				 if True, the copy is optimized, the dynamic type of the factor can changed
 				 if False, the dynamic type stay the same
 				 (default value false)*/
-				void push_first(const Faust::MatGeneric<FPP,Cpu>* M, const bool optimizedCopy=false);
+				void push_first(const Faust::MatGeneric<FPP,Cpu>* M, const bool optimizedCopy=false, const bool conjugate=false, const bool copying=true);
+                void insert(faust_unsigned_int i, Faust::MatGeneric<FPP,Cpu>* M);
+                void pop_back();
+                void pop_front();
+                void erase(faust_unsigned_int i);
+                void resize(faust_unsigned_int size);
 				//void pop_back(Faust::MatGeneric<FPP,Cpu>* M);
 				//void pop_first(Faust::MatGeneric<FPP,Cpu>* M);
 				//void pop_first(Faust::MatGeneric<FPP,Cpu>* M) const;
@@ -249,11 +255,14 @@ namespace Faust
 #ifdef FAUST_VERBOSE
 					std::cout << "~Transform()" << std::endl;
 #endif
-					for (int i=0;i<data.size();i++)
-						if(this->dtor_delete_data)
-							delete data[i];
-						else
-							ref_man.release(data[i]);
+					if(! this->dtor_disabled)
+					{
+						for (int i=0;i<data.size();i++)
+							if(this->dtor_delete_data)
+								delete data[i];
+							else
+								ref_man.release(data[i]);
+					}
 				}
 
 				/*!
@@ -305,6 +314,10 @@ namespace Faust
 		void print_timers() const;
 #endif
 
+		// should be used very carefully for testing purpose only
+        // (it allows memory leaks)
+		void disable_dtor() { this->dtor_disabled = true; }
+		void enable_dtor() { this->dtor_disabled = false; }
 			private:
 		void disable_data_deletion() { this->dtor_delete_data = false; }
 		void enable_data_deletion() { this->dtor_delete_data = true; }
@@ -314,6 +327,7 @@ namespace Faust
 		static const char * m_className;
 		std::vector<Faust::MatGeneric<FPP,Cpu>*> data;
 		bool dtor_delete_data;
+		bool dtor_disabled;
 		static RefManager ref_man;
 
 #ifdef __COMPILE_TIMERS__

src/faust_linear_operator/CPU/faust_Transform.hpp

@@ -158,7 +158,7 @@ const char * Faust::Transform<FPP,Cpu>::m_className="Faust::Transform<FPP,Cpu>";
 template<typename FPP>
 Faust::Transform<FPP,Cpu>::Transform() :
 	data(std::vector<Faust::MatGeneric<FPP,Cpu>*>()),
-	totalNonZeros(0), dtor_delete_data(false)
+	totalNonZeros(0), dtor_delete_data(false), dtor_disabled(false)
 {
 
 #ifdef __COMPILE_TIMERS__
@@ -171,7 +171,7 @@ Faust::Transform<FPP,Cpu>::Transform() :
 template<typename FPP>
 Faust::Transform<FPP,Cpu>::Transform(const Faust::Transform<FPP,Cpu> & A) :
 	data(std::vector<Faust::MatGeneric<FPP,Cpu>*>()),
-	totalNonZeros(A.totalNonZeros), dtor_delete_data(false)
+	totalNonZeros(A.totalNonZeros), dtor_delete_data(false), dtor_disabled(false)
 {
 	data.resize(0); // to be sure
 	*this = A; // rely on assignment operator (avoid duplicate)
@@ -184,7 +184,7 @@ Faust::Transform<FPP,Cpu>::Transform(const Faust::Transform<FPP,Cpu> & A) :
 template<typename FPP>
 Faust::Transform<FPP,Cpu>::Transform(const std::vector<Faust::MatGeneric<FPP,Cpu> *> & facts, const FPP lambda_ /*default value = 1.0 */,const bool optimizedCopy /*default value = false*/, const bool cloning_fact /* default to true */) :
 	data(std::vector<Faust::MatGeneric<FPP,Cpu>*>()),
-	totalNonZeros(0), dtor_delete_data(false)
+	totalNonZeros(0), dtor_delete_data(false), dtor_disabled(false)
 {
 	data.resize(facts.size());
 	if(data.size() > 0) {
@@ -219,7 +219,7 @@ Faust::Transform<FPP,Cpu>::Transform(const std::vector<Faust::MatGeneric<FPP,Cpu
 
 
 template<typename FPP>
-Faust::Transform<FPP,Cpu>::Transform(Faust::Transform<FPP,Cpu>&& T) : totalNonZeros(T.totalNonZeros), dtor_delete_data(T.dtor_delete_data)
+Faust::Transform<FPP,Cpu>::Transform(Faust::Transform<FPP,Cpu>&& T) : totalNonZeros(T.totalNonZeros), dtor_delete_data(T.dtor_delete_data), dtor_disabled(T.dtor_disabled)
 {
 	data = std::move(T.data);
 	totalNonZeros = T.totalNonZeros;
@@ -231,7 +231,7 @@ Faust::Transform<FPP,Cpu>::Transform(Faust::Transform<FPP,Cpu>&& T) : totalNonZe
 
 template<typename FPP>
 Faust::Transform<FPP,Cpu>::Transform(const std::vector<Faust::MatDense<FPP,Cpu> >&facts, const bool optimizedCopy /*default value = false*/ ):	data(std::vector<Faust::MatGeneric<FPP,Cpu>*>()),
-	totalNonZeros(0), dtor_delete_data(false)
+	totalNonZeros(0), dtor_delete_data(false), dtor_disabled(false)
 {
 	data.resize(facts.size());
 	for (int i=0 ; i<data.size() ; i++)
@@ -245,7 +245,7 @@ Faust::Transform<FPP,Cpu>::Transform(const std::vector<Faust::MatDense<FPP,Cpu>
 
 template<typename FPP>
 Faust::Transform<FPP,Cpu>::Transform(const std::vector<Faust::MatSparse<FPP,Cpu> >& facts, const bool optimizedCopy /*default value = false*/ ):	data(std::vector<Faust::MatGeneric<FPP,Cpu>*>()),
-	totalNonZeros(0), dtor_delete_data(false)
+	totalNonZeros(0), dtor_delete_data(false), dtor_disabled(false)
 {
 	data.resize(facts.size());
 	for (int i=0 ; i<data.size() ; i++)
@@ -258,7 +258,7 @@ Faust::Transform<FPP,Cpu>::Transform(const std::vector<Faust::MatSparse<FPP,Cpu>
 
 template<typename FPP>
 Faust::Transform<FPP,Cpu>::Transform(const Transform<FPP, Cpu>* A, const bool transpose_A, const bool conj_A, const Transform<FPP, Cpu>* B, const bool transpose_B, const bool conj_B):
-	data(std::vector<Faust::MatGeneric<FPP,Cpu>*>()), totalNonZeros(0), dtor_delete_data(false)
+	data(std::vector<Faust::MatGeneric<FPP,Cpu>*>()), totalNonZeros(0), dtor_delete_data(false), dtor_disabled(false)
 {
 	data.resize(A->size()+B->size());
 	int i = transpose_A?A->size()-1:0;
@@ -459,6 +459,14 @@ void Faust::Transform<FPP,Cpu>::updateNonZeros()
 
 template<typename FPP>
 Faust::MatDense<FPP,Cpu> Faust::Transform<FPP,Cpu>::get_product(const char opThis, const bool isConj)const
+{
+    Faust::MatDense<FPP,Cpu> mat;
+    this->get_product(mat, opThis, isConj);
+    return mat;
+}
+
+template<typename FPP>
+void Faust::Transform<FPP,Cpu>::get_product(Faust::MatDense<FPP,Cpu> &mat, const char opThis, const bool isConj)const
 {
 
 
@@ -478,12 +486,9 @@ Faust::MatDense<FPP,Cpu> Faust::Transform<FPP,Cpu>::get_product(const char opThi
 
 	prod.setEyes();
 
-	Faust::MatDense<FPP,Cpu> p = this->multiply(prod, opThis);
-
-	if(isConj && opThis != 'H') p.conjugate();
-
-	return p;
+	mat = this->multiply(prod, opThis);
 
+	if(isConj && opThis != 'H') mat.conjugate();
 
 	/* modif NB v1102 : factor are no longer MatSparse, they are MatGeneric now
 	   Faust::MatDense<FPP,Cpu> prod(data[0].getNbRow());
@@ -656,6 +661,7 @@ Faust::Transform<FPP,Cpu>& Faust::Transform<FPP,Cpu>::operator=(Faust::Transform
 {
 	data = std::move(T.data);
 	totalNonZeros = T.totalNonZeros;
+	dtor_disabled = T.dtor_disabled;
 	dtor_delete_data = T.dtor_delete_data;
 	T.data.resize(0);
 	return *this;
@@ -902,7 +908,7 @@ void Faust::Transform<FPP,Cpu>::push_back(const Faust::MatGeneric<FPP,Cpu>* M, c
 
 
 	template<typename FPP>
-void Faust::Transform<FPP,Cpu>::push_first(const Faust::MatGeneric<FPP,Cpu>* M, const bool optimizedCopy /*default value = false */)
+void Faust::Transform<FPP,Cpu>::push_first(const Faust::MatGeneric<FPP,Cpu>* M, const bool optimizedCopy /*default value = false */, const bool conjugate, const bool copying /* default to true */ )
 {
 	if (size()>0)
 		if(this->getNbRow()!=M->getNbCol() || M->getNbRow()<1)
@@ -910,16 +916,68 @@ void Faust::Transform<FPP,Cpu>::push_first(const Faust::MatGeneric<FPP,Cpu>* M,
 			handleError(m_className,"push_first : incorrect dimensions");
 		}
 
-	Faust::MatGeneric<FPP,Cpu>* M_copy = M->Clone(optimizedCopy);
-	data.insert(data.begin(),M_copy);
-	if(!dtor_delete_data) ref_man.acquire(M_copy);
-	totalNonZeros += M_copy->getNonZeros();
+	Faust::MatGeneric<FPP,Cpu>* M_;
+	if(copying)
+	{
+		M_ = M->Clone(optimizedCopy);
+		if(conjugate) M_->conjugate();
+	}
+	else
+	{
+		if(conjugate || optimizedCopy) throw runtime_error("copying argument mustn't be true if any of optimizedCopy or conjugate is true.");
+		M_ = const_cast<Faust::MatGeneric<FPP,Cpu>*>(M);
+	}
+	data.insert(data.begin(),M_);
+	if(!dtor_delete_data) ref_man.acquire(M_);
+	totalNonZeros += M_->getNonZeros();
 
 #ifdef __COMPILE_TIMERS__
 	this->t_multiply_vector.insert(this->t_multiply_vector().begin(),Faust::Timer());
 #endif
 }
 
+template<typename FPP>
+void Faust::Transform<FPP,Cpu>::insert(faust_unsigned_int i, Faust::MatGeneric<FPP,Cpu>* M)
+{
+    if(i > size()) throw out_of_range("Faust::Transform<FPP,Cpu>::insert");
+	data.insert(data.begin()+i,M);
+	if(!dtor_delete_data) ref_man.acquire(M);
+	totalNonZeros += M->getNonZeros();
+}
+
+template<typename FPP>
+void Faust::Transform<FPP,Cpu>::pop_back()
+{
+    this->erase(size()-1);
+}
+
+template<typename FPP>
+void Faust::Transform<FPP,Cpu>::pop_front()
+{
+    this->erase(0);
+}
+
+template<typename FPP>
+void Faust::Transform<FPP,Cpu>::erase(faust_unsigned_int i)
+{
+    if(i >= size()) throw out_of_range("Faust::Transform<FPP,Cpu>::erase");
+    totalNonZeros -= (*(begin()+i))->getNonZeros();
+    if(!dtor_delete_data) ref_man.release(*(begin()+i));
+    data.erase(begin()+i);
+}
+
+template<typename FPP>
+void Faust::Transform<FPP,Cpu>::resize(faust_unsigned_int size)
+{
+    if(size < this->size())
+        this->erase(size);
+    else if(this->size() != size)
+    {
+        // size > size()
+        this->data.resize(size);
+    }
+}
+
 /*
    template<typename FPP>
    void Faust::Transform<FPP,Cpu>::pop_back(Faust::MatGeneric<FPP,Cpu>* M)

src/faust_linear_operator/CPU/faust_TransformHelper.h

@@ -116,13 +116,16 @@ namespace Faust {
 
 			void push_back_();
 			void push_back(const MatGeneric<FPP,Cpu>* M, const bool optimizedCopy=false, const bool copying=true);
-
+            void pop_back();
+            void pop_front();
+            void push_first(const Faust::MatGeneric<FPP,Cpu>* M, const bool optimizedCopy=false, const bool copying=true);
 			faust_unsigned_int getNbRow() const;
 			faust_unsigned_int getNbCol() const;
 			faust_unsigned_int getNBytes() const;
 			bool isReal() const;
 			faust_unsigned_int get_total_nnz() const;
 			faust_unsigned_int size() const;
+			void resize(faust_unsigned_int);
 			void display() const;
 			string to_string() const;
 			faust_unsigned_int get_fact_nnz(const faust_unsigned_int id) const;
@@ -160,6 +163,7 @@ namespace Faust {
 					faust_unsigned_int start_col_id, faust_unsigned_int end_col_id);
 			TransformHelper<FPP, Cpu>* fancy_index(faust_unsigned_int* row_ids, faust_unsigned_int num_rows, faust_unsigned_int* col_ids, faust_unsigned_int num_cols);
 			MatDense<FPP,Cpu> get_product() const;
+            void get_product(Faust::MatDense<FPP,Cpu>& prod) const;
 			void save_mat_file(const char* filename) const;
 			double spectralNorm(const int nbr_iter_max, double threshold, int &flag) const;
 			TransformHelper<FPP,Cpu>* transpose();
@@ -196,7 +200,9 @@ namespace Faust {
 			void pack_factors(faust_unsigned_int start_id, faust_unsigned_int end_id);
 			void pack_factors();
 			void pack_factors(const faust_unsigned_int id, const PackDir dir);
-
+            /** for testing purpose only (memory leaks enabled) */
+			void disable_dtor() { this->transform->disable_dtor(); }
+            void enable_dtor() { this->transform->enable_dtor(); }
 			static TransformHelper<FPP,Cpu>* randFaust(RandFaustType t, unsigned int min_num_factors, unsigned int max_num_factors, unsigned int min_dim_size, unsigned int max_dim_size, float density=.1f, bool per_row=true);
 			static TransformHelper<FPP,Cpu>* hadamardFaust(unsigned int n, const bool norma=true);
 			static TransformHelper<FPP,Cpu>* fourierFaust(unsigned int n, const bool norma=true);

src/faust_linear_operator/CPU/faust_TransformHelper.hpp

@@ -539,6 +539,18 @@ namespace Faust {
 			return th;
 		}
 
+	template<typename FPP>
+		void TransformHelper<FPP,Cpu>::pop_back()
+        {
+            transform->pop_back();
+        }
+
+	template<typename FPP>
+		void TransformHelper<FPP,Cpu>::pop_front()
+        {
+            transform->pop_front();
+        }
+
 	template<typename FPP>
 		void TransformHelper<FPP,Cpu>::push_back(const MatGeneric<FPP,Cpu>* M, const bool optimizedCopy /* false by default */, const bool copying /* true to default */)
 		{
@@ -546,6 +558,14 @@ namespace Faust {
 			//this function is here only for python wrapper (TODO: see how to modify that wrapper in order to delete this function after or just use it internally -- not py/matfaust)
 			this->transform->push_back(M, optimizedCopy, is_conjugate, copying); //2nd argument is for opt. (possibly converting dense <-> sparse)
 		}
+    
+	template<typename FPP>
+		void TransformHelper<FPP,Cpu>::push_first(const MatGeneric<FPP,Cpu>* M, const bool optimizedCopy /* false by default */, const bool copying /* true to default */)
+		{
+			//warning: should not be called after initialization of factors (to respect the immutability property)
+			//this function is here only for python wrapper (TODO: see how to modify that wrapper in order to delete this function after or just use it internally -- not py/matfaust)
+			this->transform->push_first(M, optimizedCopy, is_conjugate, copying); //2nd argument is for opt. (possibly converting dense <-> sparse)
+		}
 
 	template<typename FPP>
 		template<typename Head, typename ... Tail>
@@ -627,6 +647,12 @@ namespace Faust {
 			return this->transform->size();
 		}
 
+	template<typename FPP>
+		void TransformHelper<FPP,Cpu>::resize(faust_unsigned_int sz)
+		{
+			return this->transform->resize(sz);
+		}
+
 	template<typename FPP>
 		void TransformHelper<FPP,Cpu>::display() const
 		{
@@ -929,6 +955,11 @@ namespace Faust {
 			return this->transform->get_product(isTransposed2char(), is_conjugate);
 		}
 
+	template<typename FPP>
+		void TransformHelper<FPP,Cpu>::get_product(Faust::MatDense<FPP,Cpu>& prod) const {
+			this->transform->get_product(prod, isTransposed2char(), is_conjugate);
+		}
+
 	template<typename FPP>
 		void TransformHelper<FPP,Cpu>::save_mat_file(const char* filename) const
 		{
@@ -1458,26 +1489,38 @@ namespace Faust {
 				throw out_of_range("start_id is out of range.");
 			if(end_id < 0 || end_id >= size())
 				throw out_of_range("end_id is out of range.");
+            Faust::MatDense<FPP,Cpu> * packed_fac = nullptr;
 			if(end_id == start_id)
-				//nothing to do
-				return;
-			// we have to multiply factors from start_id to end_id into one matrix
-			// simple way to do, 1) create a overhead-free TransformHelper with these factors
-			// 2) call get_product() to override the start_id factors with the result on the end
-			// 3) release through ref_man the factors to pack, from start_id to end_id(that's the dirty part to maybe enhance because normally this is the Transform object's responsibility). The packed factor must be acquired too.
-			// 4) erase factors from start_id to end_id
-			// 1)
-			std::vector<Faust::MatGeneric<FPP,Cpu>*> topack_factors(begin()+start_id, begin()+end_id+1);
-			Faust::TransformHelper<FPP,Cpu> t(topack_factors, 1.0, false, false, true);
-			// 2)
-			Faust::MatDense<FPP,Cpu> * packed_fac = new MatDense<FPP,Cpu>(t.get_product());
+            {
+                //nothing to do except converting to MatDense if start_id
+                //factor is a MatSparse
+                packed_fac = dynamic_cast<Faust::MatDense<FPP,Cpu>*>(*(begin()+start_id));
+                if(packed_fac == nullptr) // factor start_id is not at MatDense, convert it
+                    packed_fac = new MatDense<FPP,Cpu>(*dynamic_cast<Faust::MatSparse<FPP,Cpu>*>(*(begin()+start_id)));
+                else
+                    return; //no change
+            }
+            else
+            {
+                // we have to multiply factors from start_id to end_id into one matrix
+                // simple way to do, 1) create a overhead-free TransformHelper with these factors
+                // 2) call get_product() to override the start_id factors with the result on the end
+                // 3) erase factors from start_id to end_id and insert packed factor too to replace them (that's Transform object responsibility).
+                // 1)
+                std::vector<Faust::MatGeneric<FPP,Cpu>*> topack_factors(begin()+start_id, begin()+end_id+1);
+                Faust::TransformHelper<FPP,Cpu> t(topack_factors, 1.0, false, false, true);
+                // 2)
+                packed_fac = new MatDense<FPP,Cpu>(t.get_product());
+            }
 			// 3)
-			for(auto f = begin()+start_id; f != begin()+end_id+1; f++)
-				Transform<FPP,Cpu>::ref_man.release(*f);
-			// 4)
-			this->transform->data.erase(this->transform->begin()+start_id, this->transform->begin()+end_id+1);
-			this->transform->data.insert(begin()+start_id,  packed_fac);
-			Transform<FPP,Cpu>::ref_man.acquire(packed_fac);
+            faust_unsigned_int i = end_id;
+            while(i>=start_id)
+            {
+                this->transform->erase(i);
+                if(i == 0) break;
+                i--;
+            }
+            this->transform->insert(start_id, packed_fac);
 		}
 
 	template <typename FPP> void TransformHelper<FPP,Cpu>::pack_factors(const faust_unsigned_int id, const PackDir dir)

src/faust_linear_operator/CPU/faust_linear_algebra.hpp

@@ -307,7 +307,7 @@ void Faust::gemv(const Faust::MatDense<FPP,Cpu> & A,const Faust::Vect<FPP,Cpu> &
 	{
 		if(A.is_id() && alpha == FPP(1.0))
 		{
-			std::cout << "gemv identity opt." << std::endl;
+//			std::cout << "gemv identity opt." << std::endl;
 			memcpy(y.getData(), px->getData(), sizeof(FPP)*nbRowOpA);
 		}
 		else if (typeA == 'N')