diff --git a/misc/test/CMakeLists.txt b/misc/test/CMakeLists.txt
index 29f447d44a93b3c45b81deda9a59a9e064510f35..995499548a5bbcc75de0f16412b0a3eb2794a0e2 100755
--- a/misc/test/CMakeLists.txt
+++ b/misc/test/CMakeLists.txt
@@ -156,7 +156,7 @@ endif()
if(NOT NOCPPTESTS)
foreach(TEST_FPP float double)
- foreach(FILE faust_mult faust_mult_cplx test_Vect_min test_MatDense_get_row test_MatDense_lower_upper_tri test_MatDense_nonzeros_indices test_Transform_move test_MatDense_min faust_transform_omp_mul faust_pruneout faust_transform_optimize_storage faust_transform_optimize)
+ foreach(FILE faust_mult faust_mult_cplx test_Vect_min test_MatDense_get_row test_MatDense_lower_upper_tri test_MatDense_nonzeros_indices test_Transform_move test_MatDense_min faust_transform_omp_mul faust_pruneout faust_transform_optimize_storage faust_transform_optimize faust_prox_blockdiag)
set(TEST_BIN_FILE ${FILE}_${TEST_FPP})
set(TEST_FILE_CPP ${TEST_BIN_FILE}.cpp)
message(STATUS ${TEST_FILE_CPP})
diff --git a/misc/test/src/C++/unit/faust_prox_blockdiag.cpp.in b/misc/test/src/C++/unit/faust_prox_blockdiag.cpp.in
new file mode 100644
index 0000000000000000000000000000000000000000..ecbceb0a40f3d400095ed0495378bf3c46e9dd3e
--- /dev/null
+++ b/misc/test/src/C++/unit/faust_prox_blockdiag.cpp.in
@@ -0,0 +1,107 @@
+/****************************************************************************/
+/* Description: */
+/* unitary test for testing multiplication by faust with real scalar */
+/* */
+/* For more information on the FAuST Project, please visit the website */
+/* of the project : <http://faust.inria.fr> */
+/* */
+/* License: */
+/* Copyright (2016): 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 the GNU Affero */
+/* General Public License. */
+/* This program is free software: you can redistribute it and/or modify */
+/* it under the terms of the GNU Affero General Public License as */
+/* published by the Free Software Foundation. */
+/* */
+/* 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. */
+/* See the GNU Affero General Public License for more details. */
+/* */
+/* You should have received a copy of the GNU Affero General Public */
+/* License along with this program. */
+/* If not, see <http://www.gnu.org/licenses/>. */
+/* */
+/* 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 "faust_MatSparse.h"
+#include "faust_HierarchicalFact.h"
+#include "faust_Timer.h"
+#include "faust_Transform.h"
+#include <string>
+#include <sstream>
+#include "faust_BlasHandle.h"
+#include "faust_SpBlasHandle.h"
+
+
+#include <iostream>
+#include <iomanip>
+
+/** \brief unitary test for testing multiplication by faust
+*/
+
+typedef @TEST_FPP@ FPP;
+
+
+//using namespace Faust;
+using namespace std;
+
+int main(int argc, char* argv[])
+{
+ if (typeid(FPP) == typeid(double))
+ {
+ cout<<"floating point precision == double"<<endl;
+ }
+
+ if (typeid(FPP) == typeid(float))
+ {
+ cout<<"floating point precision == float"<<endl;
+ }
+
+ int dim1 = 38;
+ int dim2 = 42;
+
+ Faust::MatDense<FPP,Cpu> * M = Faust::MatDense<FPP,Cpu>::randMat(dim1, dim2);
+
+ Faust::MatDense<FPP,Cpu> M_copy = *M;
+
+ std::vector<faust_unsigned_int> m_vec = {5u, 12u, 25u, 38u};
+ std::vector<faust_unsigned_int> n_vec = {9u, 18u, 32u, 42u};
+
+ Faust::prox_blockdiag(M_copy, m_vec, n_vec);
+ cout << "====================== M:" << endl;
+ for(faust_unsigned_int i = 0; i < M->getNbRow(); i++)
+ {
+ for(faust_unsigned_int j = 0; j < M->getNbCol(); j++)
+ {
+ printf("%3.2f ", (*M)(i,j));
+ }
+ printf("\n");
+ }
+ cout << "====================== M_:"<< endl;
+ for(faust_unsigned_int i = 0; i < M->getNbRow(); i++)
+ {
+ for(faust_unsigned_int j = 0; j < M->getNbCol(); j++)
+ {
+ printf("%3.2f ", M_copy(i,j));
+ assert(M_copy(i,j) == 0.0 || M_copy(i,j) == (*M)(i,j));
+ }
+ printf("\n");
+ }
+ return 0;
+}
diff --git a/src/faust_linear_operator/CPU/faust_MatDense.h b/src/faust_linear_operator/CPU/faust_MatDense.h
index 33662da647f7e6c606b57d34d47d5055b4a88287..f153294880eeea81c4df00bb68e5e6d10ef77be7 100644
--- a/src/faust_linear_operator/CPU/faust_MatDense.h
+++ b/src/faust_linear_operator/CPU/faust_MatDense.h
@@ -417,6 +417,10 @@ void spgemm(const Faust::MatSparse<FPP,Cpu> & A,const Faust::MatDense<FPP,Cpu> &
Faust::MatDense<FPP,Cpu>* get_rows(faust_unsigned_int start_row_id, faust_unsigned_int num_rows) const;
Faust::MatDense<FPP,Cpu>* get_rows(faust_unsigned_int* row_ids, faust_unsigned_int n) const;
+ Faust::MatDense<FPP,Cpu> get_block(faust_unsigned_int i, faust_unsigned_int j, faust_unsigned_int nrows, faust_unsigned_int ncols);
+
+ void set_block(faust_unsigned_int i, faust_unsigned_int j, Faust::MatDense<FPP,Cpu> & block);
+
FPP min_coeff() const;
Faust::Vect<FPP,Cpu> rowwise_min() const;
Faust::Vect<FPP,Cpu> rowwise_min(int* col_indices) const;
diff --git a/src/faust_linear_operator/CPU/faust_MatDense.hpp b/src/faust_linear_operator/CPU/faust_MatDense.hpp
index cb21c038d5997ca3b7e80b468e2700da84f935e0..dc6d208300e95afed098285ecc1f060fd6d482e5 100644
--- a/src/faust_linear_operator/CPU/faust_MatDense.hpp
+++ b/src/faust_linear_operator/CPU/faust_MatDense.hpp
@@ -1136,6 +1136,23 @@ Faust::MatDense<FPP,Cpu>* Faust::MatDense<FPP,Cpu>::get_rows(faust_unsigned_int*
return rows;
}
+template<typename FPP>
+Faust::MatDense<FPP,Cpu> Faust::MatDense<FPP,Cpu>::get_block(faust_unsigned_int i, faust_unsigned_int j, faust_unsigned_int nrows, faust_unsigned_int ncols)
+{
+ Faust::MatDense<FPP,Cpu> block_mat(nrows, ncols);
+ block_mat.mat.block(0,0, nrows, ncols) = mat.block(i, j, nrows, ncols);
+ return block_mat;
+}
+
+template<typename FPP>
+void Faust::MatDense<FPP,Cpu>::set_block(faust_unsigned_int i, faust_unsigned_int j, Faust::MatDense<FPP,Cpu> & block_mat)
+{
+ faust_unsigned_int nrows, ncols;
+ nrows = block_mat.getNbRow();
+ ncols = block_mat.getNbCol();
+ mat.block(i, j, nrows, ncols) = block_mat.mat.block(0, 0, nrows, ncols);
+}
+
template<typename FPP>
FPP Faust::MatDense<FPP,Cpu>::min_coeff() const
{
diff --git a/src/faust_linear_operator/CPU/faust_prox.h b/src/faust_linear_operator/CPU/faust_prox.h
index ef5c96f4987f2de77165f1d511e131f4914591d7..0e6d79169da21d8717b587b14fac79f59745e3d9 100644
--- a/src/faust_linear_operator/CPU/faust_prox.h
+++ b/src/faust_linear_operator/CPU/faust_prox.h
@@ -79,6 +79,9 @@ namespace Faust {
void prox_normlin(Faust::MatDense<FPP,Cpu> & M,FPP2 s, const bool normalized=false, const bool pos=false);
// template<typename FPP>
// void prox_blkdiag(Faust::MatDense<FPP,Cpu> & M,faust_unsigned_int k);
+ template<typename FPP>
+ void prox_blockdiag(Faust::MatDense<FPP,Cpu> & M, std::vector<faust_unsigned_int>& m_vec, std::vector<faust_unsigned_int>& n_vec, const bool normalized = false, const bool pos = false);
+
template<typename FPP> void pre_prox_pos(MatDense<FPP,Cpu> & M);
template<typename FPP> void prox_hankel(Faust::MatDense<FPP, Cpu> & M);
template<typename FPP> void prox_toeplitz(Faust::MatDense<FPP, Cpu> & M);
diff --git a/src/faust_linear_operator/CPU/faust_prox.hpp b/src/faust_linear_operator/CPU/faust_prox.hpp
index c15cc9a695e2903f74382a228f3dcee92c1cf481..8b7a59e907a00dc990c5b0205ef26b5391ba1707 100644
--- a/src/faust_linear_operator/CPU/faust_prox.hpp
+++ b/src/faust_linear_operator/CPU/faust_prox.hpp
@@ -240,6 +240,36 @@ void Faust::prox_normlin(Faust::MatDense<FPP,Cpu> & M,FPP2 s, const bool normali
M.transpose();
}
+template<typename FPP>
+void Faust::prox_blockdiag(Faust::MatDense<FPP,Cpu> & M, std::vector<faust_unsigned_int> & m_vec, std::vector<faust_unsigned_int>& n_vec, const bool normalized /* default to false */, const bool pos)
+{
+// if(M.shape != self.shape): raise ValueError('The dimension of the '
+// 'projector and matrix must '
+// 'agree.')
+// M_ = np.zeros(M.shape)
+// m_ = 0
+// n_ = 0
+// for i,(m,n) in enumerate(zip(self.m_vec, self.n_vec)):
+// print("i=", i, "m=", m, "n=", n)
+// M_[m_:m,n_:n] = M[m_:m,n_:n]
+// m_ = m
+// n_ = n
+// return M_
+//
+ Faust::MatDense<FPP,Cpu> M_(M.getNbRow(), M.getNbCol());
+ M_.setZeros();
+ faust_unsigned_int m_ = 0, n_ = 0;
+ for(faust_unsigned_int i=0; i< m_vec.size();i++)
+ {
+ Faust::MatDense<FPP,Cpu> block_mat = M.get_block(m_, n_, m_vec[i]-m_, n_vec[i]-n_);
+ M_.set_block(m_, n_, block_mat);
+ m_ = m_vec[i];
+ n_ = n_vec[i];
+ }
+ M = M_;
+}
+
+
template<typename FPP>
void
Faust::pre_prox_pos(MatDense<FPP,Cpu> & M)