Update PALM4MSA/hierachical 2020 C++ impl.:

- new argument id (cf. code doc), used in hiearchical calls to palm.
- checking S fro. norm is not zero, otherwise raising an exception before computing lambda (whose the denominator is n²).

src/algorithm/factorization/faust_hierarchical.hpp

@@ -102,7 +102,7 @@ Faust::TransformHelper<FPP,DEVICE>* Faust::hierarchical(const Faust::MatDense<FP
 		}
 		else tmp_sparse = nullptr;
 		Faust::palm4msa2(*tmp_dense, Si_cons, Si_th, lambda_, p.stop_crit_2facts, is_update_way_R2L , use_csr, packing_RL, compute_2norm_on_array,
-				norm2_threshold, norm2_max_iter, constant_step_size, step_size, on_gpu, p.isVerbose);
+				norm2_threshold, norm2_max_iter, constant_step_size, step_size, on_gpu, p.isVerbose, i+1);
 		if(tmp_sparse != nullptr)
 			// the Si factor has been converted into a MatDense in the memory
 			// storage
@@ -138,7 +138,7 @@ Faust::TransformHelper<FPP,DEVICE>* Faust::hierarchical(const Faust::MatDense<FP
 
 		// global optimization
 		Faust::palm4msa2(A, glo_cons, *S, glo_lambda, p.stop_crit_global ,is_update_way_R2L, use_csr, packing_RL, compute_2norm_on_array,
-				norm2_threshold, norm2_max_iter, constant_step_size, step_size, on_gpu, p.isVerbose);
+				norm2_threshold, norm2_max_iter, constant_step_size, step_size, on_gpu, p.isVerbose, i+1);
 	}
 	lambda = glo_lambda;
 	return S;

src/algorithm/factorization/faust_palm4msa2020.h

@@ -64,7 +64,11 @@ namespace Faust
 				const unsigned int norm2_max_iter=FAUST_NORM2_MAX_ITER,
 				const bool constant_step_size=false, const Real<FPP> step_size=FAUST_PRECISION,
 				const bool on_gpu=false,
-				const bool is_verbose=false);
+				const bool is_verbose=false,
+				/* id argument is useful to identify the palm4msa call.
+				 * For example, hierarchical use it to indicate the number of the current factorization.
+				 * */
+				const int id=0);
 
 	/**
 	 * \brief Fill S with nfacts eye() matrices.

src/algorithm/factorization/faust_palm4msa2020_2.hpp

@@ -134,7 +134,7 @@ void Faust::palm4msa2(const Faust::MatDense<FPP,DEVICE>& A,
 		const unsigned int norm2_max_iter,
 		const bool constant_step_size, const Real<FPP> step_size,
 		const bool on_gpu /*=false*/,
-		const bool is_verbose/*=false*/)
+		const bool is_verbose/*=false*/, const int id/*=0*/)
 {
 	std::chrono::time_point<std::chrono::high_resolution_clock> spectral_stop, spectral_start;
 	std::chrono::duration<double> spectral_duration = std::chrono::duration<double>::zero();
@@ -299,8 +299,11 @@ void Faust::palm4msa2(const Faust::MatDense<FPP,DEVICE>& A,
 			}
 			// really update now
 			constraints[f_id]->project<FPP,DEVICE,Real<FPP>>(D);
+
+
 			if(use_csr && dcur_fac != nullptr || !use_csr && scur_fac != nullptr)
 				throw std::runtime_error("Current factor is inconsistent with use_csr.");
+
 			if(use_csr)
 			{
 				spD = D;
@@ -322,8 +325,10 @@ void Faust::palm4msa2(const Faust::MatDense<FPP,DEVICE>& A,
 		//		Faust::MatDense<FPP,DEVICE> A_H_S = A_H_S_.get_product();
 		Real<FPP> trr = std::real(A_H_S.trace());
 		Real<FPP> n = S.normFro();
-		lambda = trr/(n*n);
-		//		std::cout << "debug lambda: " << lambda << std::endl;
+		if(std::numeric_limits<Real<FPP>>::epsilon() >= n)
+			throw std::runtime_error("Faust Frobenius norm is zero, can't compute lambda.");
+		lambda = trr/(n*n); //TODO: raise exception if n == 0
+//		std::cout << "debug lambda: " << lambda << std::endl;
 		i++;
 	}
 	S.update_total_nnz();