diff --git a/benchmarks/Makefile.am b/benchmarks/Makefile.am
index 161626145..9b5b92d32 100755
--- a/benchmarks/Makefile.am
+++ b/benchmarks/Makefile.am
@@ -35,7 +35,7 @@ endif
 
 PERFPUBLISHERFILE=benchmarks-report.xml
 
-FFLA_BENCH =    benchmark-fgemm benchmark-fgemm-rns benchmark-wino benchmark-ftrsm  benchmark-fgesv benchmark-ftrsv benchmark-ftrtri  benchmark-inverse benchmark-fsytrf benchmark-fsyrk  benchmark-lqup benchmark-pluq benchmark-charpoly benchmark-charpoly-mp benchmark-fgemm-mp benchmark-fgemv-mp benchmark-ftrsm-mp benchmark-lqup-mp benchmark-checkers benchmark-fadd-lvl2 benchmark-fdot benchmark-fgemv 
+FFLA_BENCH =    benchmark-fgemm benchmark-fgemm-rns benchmark-wino benchmark-ftrsm  benchmark-fgesv benchmark-ftrsv benchmark-ftrtri  benchmark-inverse benchmark-fsytrf benchmark-fsyrk  benchmark-lqup benchmark-pluq benchmark-charpoly benchmark-charpoly-mp benchmark-fgemm-mp benchmark-fgemv-mp benchmark-ftrsm-mp benchmark-lqup-mp benchmark-checkers benchmark-fadd-lvl2 benchmark-fdot benchmark-fgemv benchmark-fgemv-rns
 
 BLAS_BENCH =    benchmark-sgemm$(EXEEXT) benchmark-dgemm benchmark-dtrsm
 LAPA_BENCH =    benchmark-dtrtri benchmark-dgetri benchmark-dgetrf benchmark-dsytrf
@@ -85,6 +85,7 @@ benchmark_checkers_SOURCES = benchmark-checkers.C
 benchmark_fadd_lvl2_SOURCES = benchmark-fadd-lvl2.C
 benchmark_fdot_SOURCES = benchmark-fdot.C
 benchmark_fgemv_SOURCES = benchmark-fgemv.C
+benchmark_fgemv_rns_SOURCES = benchmark-fgemv-rns.C
 
 benchmark_sgemm_CXXFLAGS = $(AM_CXXFLAGS) -D__SGEMM__
 
diff --git a/benchmarks/benchmark-fgemv-mp.C b/benchmarks/benchmark-fgemv-mp.C
index 4cb25ee06..2758990ab 100644
--- a/benchmarks/benchmark-fgemv-mp.C
+++ b/benchmarks/benchmark-fgemv-mp.C
@@ -1,3 +1,4 @@
+
 /*
  * Copyright (C) FFLAS-FFPACK
  * Written by Pascal Giorgi <pascal.giorgi@lirmm.fr>
@@ -28,6 +29,7 @@
 // everywhere in the call stack
 #define __FFLASFFPACK_OPENBLAS_NT_ALREADY_SET 1
 
+
 #if not defined(MG_DEFAULT)
 #define MG_DEFAULT MG_ACTIVE
 #endif
@@ -35,7 +37,6 @@
 #define STD_RECINT_SIZE 8
 #endif
 
-
 #include "fflas-ffpack/fflas-ffpack-config.h"
 #include <iostream>
 #include <typeinfo>
@@ -43,8 +44,9 @@
 #include <string>
 using namespace std;
 
-#include "fflas-ffpack/utils/timer.h"
 #include "fflas-ffpack/fflas/fflas.h"
+#include "fflas-ffpack/utils/fflas_io.h"
+#include "fflas-ffpack/utils/timer.h"
 #include "fflas-ffpack/utils/args-parser.h"
 #include "givaro/modular-integer.h"
 #include "givaro/givcaster.h"
@@ -53,33 +55,16 @@ using namespace std;
 #include "recint/recint.h"
 #endif
 
-
-template<typename T>
-std::ostream& write_matrix(std::ostream& out, Givaro::Integer p, size_t m, size_t n, T* C, size_t ldc){
-
-    size_t www(size_t((double(p.bitsize())*log(2.))/log(10.)));
-    out<<"Matrix("<<m<<','<<n<<",[[";
-    out.width(www+1);
-    out<<std::right<<C[0];
-    for (size_t j=1;j<n;++j){
-        out<<',';
-        out.width(www);
-        out<<std::right<<C[j];
-    }
-    out<<']';
-    for (size_t i=1;i<m;++i){
-        out<<endl<<",[";
-        out.width(www+1);
-        out<<std::right<<C[i*ldc];
-        for (size_t j=1;j<n;++j){
-            out<<',';
-            out.width(www);
-            out<<std::right<<C[i*ldc+j];
-        }
-        out<<']';
-    }
-    return out<<"])";
+#ifdef	BENCH_FLINT
+#define __GMP_BITS_PER_MP_LIMB 64
+extern "C" {
+#include "flint/longlong.h"
+#include "flint/long_extras.h"
+#include "flint/fmpz_mat.h"
+#include "flint/fmpz.h"
+#include "flint/flint.h"
 }
+#endif
 
 
 static size_t iters = 3 ;
@@ -87,6 +72,7 @@ static Givaro::Integer q = -1 ;
 static unsigned long b = 512 ;
 static size_t m = 512 ;
 static size_t k = 512 ;
+static size_t n = 512 ;
 static int nbw = -1 ;
 static size_t seed= time(NULL);
 static Argument as[] = {
@@ -94,6 +80,7 @@ static Argument as[] = {
     { 'b', "-b B", "Set the bitsize of the random characteristic.",         TYPE_INT , &b },
     { 'm', "-m M", "Set the dimension m of the matrix.",                    TYPE_INT , &m },
     { 'k', "-k K", "Set the dimension k of the matrix.",                    TYPE_INT , &k },
+    { 'n', "-n N", "Set the dimension n of the matrix.",                    TYPE_INT , &n },
     { 'w', "-w N", "Set the number of winograd levels (-1 for random).",    TYPE_INT , &nbw },
     { 'i', "-i R", "Set number of repetitions.",                            TYPE_INT , &iters },
     { 's', "-s S", "Sets seed.",                            				TYPE_INT , &seed },
@@ -109,7 +96,10 @@ int tmain(){
     typedef Givaro::Modular<Ints> Field;
     Givaro::Integer p;
     FFLAS::Timer chrono, TimFreivalds;
-    double time=0.;
+    double time=0.,timev=0.;
+#ifdef BENCH_FLINT
+    double timeFlint=0.;
+#endif
     for (size_t loop=0;loop<iters;loop++){
         Givaro::Integer::random_exact_2exp(p, b);
         Givaro::IntPrimeDom IPD;
@@ -120,8 +110,8 @@ int tmain(){
         Field F(ip);
         size_t lda,ldb,ldc;
         lda=k;
-        ldb=1;
-        ldc=1;
+        ldb=n;
+        ldc=n;
 
         typename Field::RandIter Rand(F,seed);
         typename Field::Element_ptr A,B,C;
@@ -140,8 +130,8 @@ int tmain(){
         // 				Rand.random(C[i*ldc+j]);
 
         PAR_BLOCK { FFLAS::pfrand(F,Rand, m,k,A,m/size_t(MAX_THREADS)); }
-        PAR_BLOCK { FFLAS::pfrand(F,Rand, k,1,B,k/MAX_THREADS); }
-        PAR_BLOCK { FFLAS::pfzero(F, m,1,C,m/MAX_THREADS); }
+        PAR_BLOCK { FFLAS::pfrand(F,Rand, k,n,B,k/MAX_THREADS); }
+        PAR_BLOCK { FFLAS::pfzero(F, m,n,C,m/MAX_THREADS); }
 
 
         Ints alpha,beta;
@@ -149,34 +139,95 @@ int tmain(){
         beta=F.zero;
 
 
+#ifdef	BENCH_FLINT
+        // FLINT MUL //
+        fmpz_t modp,tmp;
+        fmpz_init(modp);
+        fmpz_init(tmp);
+        fmpz_set_mpz(modp, *(reinterpret_cast<const mpz_t*>(&p)));
+        fmpz_mat_t AA,BB,CC,DD;
+        fmpz_mat_init (AA, m, k);
+        fmpz_mat_init (BB, k, n);
+        fmpz_mat_init (CC, m, n);
+        fmpz_mat_init (DD, m, n);
+        fmpz_t aalpha, bbeta;
+        fmpz_set_mpz(aalpha,*(reinterpret_cast<const mpz_t*>(&alpha)));
+        fmpz_set_mpz(bbeta,*(reinterpret_cast<const mpz_t*>(&beta)));
+
+        for (size_t i=0;i<m;++i)
+            for (size_t j=0;j<k;++j)
+                fmpz_set_mpz(fmpz_mat_entry(AA,i,j),*(reinterpret_cast<const mpz_t*>(A+i*lda+j)));
+        for (size_t i=0;i<k;++i)
+            for (size_t j=0;j<n;++j)
+                fmpz_set_mpz(fmpz_mat_entry(BB,i,j),*(reinterpret_cast<const mpz_t*>(B+i*ldb+j)));
+        for (size_t i=0;i<m;++i)
+            for (size_t j=0;j<n;++j)
+                fmpz_set_mpz(fmpz_mat_entry(CC,i,j),*(reinterpret_cast<const mpz_t*>(C+i*ldc+j)));
+        chrono.clear();chrono.start();
+        // DD= A.B
+        fmpz_mat_mul(DD,AA,BB);
+        // CC = beta.C
+        fmpz_mat_scalar_mul_fmpz(CC,CC,bbeta);
+        // CC = CC + DD.alpha
+        fmpz_mat_scalar_addmul_fmpz(CC,DD,aalpha);
+        // CC = CC mod p
+        for (size_t i=0;i<m;++i)
+            for (size_t j=0;j<n;++j)
+                fmpz_mod(fmpz_mat_entry(CC,i,j),fmpz_mat_entry(CC,i,j),modp);
+
+        chrono.stop();
+        timeFlint+=chrono.usertime();
+        fmpz_mat_clear(AA);
+        fmpz_mat_clear(BB);
+#endif
+        //END FLINT CODE //
         using  FFLAS::CuttingStrategy::Recursive;
         using  FFLAS::StrategyParameter::TwoDAdaptive;
         // RNS MUL_LA
         chrono.clear();chrono.start();
+        // 		PAR_BLOCK{
+        //             FFLAS::fgemm(F,FFLAS::FflasNoTrans,FFLAS::FflasNoTrans,m,n,k,alpha,A,lda,B,ldb,beta,C,ldc, SPLITTER(NUM_THREADS,Recursive,TwoDAdaptive) );
+        // 		}
         {
-            FFLAS::ParSeqHelper::Sequential seqH;
-            FFLAS::fgemv(F,FFLAS::FflasNoTrans,m,k,alpha,A,lda,B,ldb,beta,C,ldc,seqH);
+            FFLAS::fgemm(F,FFLAS::FflasNoTrans,FFLAS::FflasNoTrans,m,n,k,alpha,A,lda,B,ldb,beta,C,ldc,FFLAS::ParSeqHelper::Sequential());
         }
+
         chrono.stop();
         time+=chrono.realtime();
 
+        TimFreivalds.start();
+        bool pass = FFLAS::freivalds(F, FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, m,n,k, alpha, A, k, B, n, C,n);
+        TimFreivalds.stop();
+        timev+=TimFreivalds.usertime();
+        if (!pass) {
+            std::cout<<"FAILED"<<std::endl;
+            std::cout << "p:=" << p << ';'<<std::endl;
+            FFLAS::WriteMatrix (std::cout<<"A:=",F,m,k,A,lda)<<';'<<std::endl;
+            FFLAS::WriteMatrix(std::cout<<"B:=",F,k,n,B,ldb)<<';'<<std::endl;
+            FFLAS::WriteMatrix(std::cout<<"C:=",F,m,n,C,ldc)<<';'<<std::endl;
+        }
+
         FFLAS::fflas_delete(A);
         FFLAS::fflas_delete(B);
         FFLAS::fflas_delete(C);
 
     }
 
-    double Mflops=((2.*double(m)-1)/1000.*double(k)/1000.0) /time * double(iters);
-    // 	Mflops*=p.bitsize()/16.;
-    cout << "Time: "<< (time/double(iters))  <<" Gfops: "<<Mflops*1.0/1000.0
+    double Gflops=(2.*double(m)/1000.*double(n)/1000.*double(k)/1000.0) / time * double(iters);
+    // 	Gflops*=p.bitsize()/16.;
+    cout  << "Time: "<< (time/double(iters))
+    <<" Gfops: "<<Gflops
     << " (total:" << time <<") "
     <<typeid(Ints).name()
-    <<" perword: "<< (Mflops*double(p.bitsize()))/64. ;
-    FFLAS::writeCommandString(std::cout << " | " << p << " (" << p.bitsize()<<")|", as)  << std::endl;
-    return 0;
-}
+    <<"  | perword: "<< (Gflops*double(p.bitsize()))/64. ;
 
+    FFLAS::writeCommandString(std::cout << '|' << p << " (" << p.bitsize()<<")|", as) << "  | Freivalds: "<< timev/double(iters) << std::endl;
 
+#ifdef BENCH_FLINT
+    cout<<"Time FLINT: "<<timeFlint<<endl;
+#endif
+    return 0;
+}
 
 int main(int argc, char** argv){
 
diff --git a/benchmarks/benchmark-fgemv-rns.C b/benchmarks/benchmark-fgemv-rns.C
new file mode 100644
index 000000000..2f7f916bf
--- /dev/null
+++ b/benchmarks/benchmark-fgemv-rns.C
@@ -0,0 +1,261 @@
+/* Copyright (c) FFLAS-FFPACK
+ * ========LICENCE========
+ * This file is part of the library FFLAS-FFPACK.
+ *
+ * FFLAS-FFPACK is free software: you can redistribute it and/or modify
+ * it under the terms of the  GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * This library 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
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this library; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
+ * ========LICENCE========
+ */
+//#define PROFILE_FGEMM_MP
+//#include "goto-def.h"
+
+#include "fflas-ffpack/fflas-ffpack-config.h"
+#include <iostream>
+#include <givaro/modular-balanced.h>
+
+#include "fflas-ffpack/config-blas.h"
+#include "fflas-ffpack/fflas/fflas.h"
+#include "fflas-ffpack/utils/timer.h"
+#include "fflas-ffpack/utils/args-parser.h"
+
+#include "fflas-ffpack/utils/fflas_io.h"
+#include "fflas-ffpack/utils/test-utils.h"
+
+#include "fflas-ffpack/utils/timer.h"
+#include "givaro/modular-integer.h"
+#include "givaro/givcaster.h"
+
+using namespace FFPACK;
+
+using namespace std;
+using namespace FFLAS;
+
+template <typename Field>
+struct need_field_characteristic { static constexpr bool value = false; };
+template <typename Field>
+struct need_field_characteristic<Givaro::Modular<Field>>{ static constexpr bool value = true; };
+template <typename Field>
+struct need_field_characteristic<Givaro::ModularBalanced<Field>>{ static constexpr bool value = true; };
+
+template <typename Field>
+struct compatible_data_type { static constexpr bool value = true; };
+template <>
+struct compatible_data_type<Givaro::ZRing<float>>{ static constexpr bool value = false; };
+template <>
+struct compatible_data_type<Givaro::ZRing<double>>{ static constexpr bool value = false; };
+
+
+template <class Field, class RandIter, class Matrix, class Vector>
+void fill_value(Field& F, RandIter& Rand,
+		Matrix& A, Vector& X, Vector& Y,
+		size_t m, size_t k, size_t incX, size_t incY, size_t lda, int NBK){
+
+  SYNCH_GROUP(
+	      FORBLOCK1D(iter, m, SPLITTER(NBK, CuttingStrategy::Row, StrategyParameter::Threads),
+			 TASK(MODE(CONSTREFERENCE(F,Rand,A)),
+			      {
+                    frand(F, Rand, iter.end()-iter.begin(), k, A+iter.begin()*lda, lda);
+			      }
+			      );
+			 );
+	      );
+  //FFLAS::pfrand(F,Rand, m,k,A,m/NBK);
+  FFLAS::frand(F,Rand, k,1,X,incX);
+  FFLAS::fzero(F, m,1,Y,incY);
+}
+
+template <class Field, class Matrix, class Vector>
+void genData(Field& F,
+	     Matrix& A, Vector& X, Vector& Y,
+	     size_t m, size_t k, size_t incX, size_t incY, size_t lda, int NBK,
+	     int bitsize, uint64_t seed){
+  typename Field::RandIter Rand(F,seed,bitsize); //Field::RandIter's parameters order has been changed between seed and bitsize
+  fill_value(F, Rand, A, X, Y, m, k, incX, incY, lda, NBK);
+}
+
+template <class Field, class Matrix, class Vector>
+bool check_result(Field& F, size_t m, size_t lda, Matrix& A, Vector& X, size_t incX, Vector& Y, size_t incY){
+  //Naive result checking by comparing result from pfgemv against the one from fgemv
+  typename Field::Element_ptr Y2 = FFLAS::fflas_new(F,m,1);
+  FFLAS::fgemv(F, FFLAS::FflasNoTrans, m, lda, F.one, A, lda, X, incX, F.zero, Y2,  incY);
+
+  for(size_t j=0; j<m; ++j){
+    if(!F.areEqual(Y2[j],Y[j])){
+      FFLAS::fflas_delete(Y2);
+      return false;
+    }
+  }
+  FFLAS::fflas_delete(Y2);
+  return true;
+}
+
+
+template <class Field, class Matrix, class Vector>
+bool benchmark_with_timer(Field& F, int p, Matrix& A, Vector& X, Vector& Y, size_t m, size_t k, size_t incX,
+			  size_t incY, size_t lda, size_t iters, int t, double& time, size_t GrainSize){
+  Timer chrono;
+  bool pass = true;
+  for (size_t i=0;i<=iters;++i){
+
+    chrono.clear();
+
+    if (p){
+
+      //typedef CuttingStrategy::Row row;
+      typedef CuttingStrategy::Recursive rec;
+      typedef StrategyParameter::Threads threads;
+      typedef StrategyParameter::Grain grain;
+
+      if (i) { chrono.start(); }
+
+      switch (p){
+          case 1:{
+	        ParSeqHelper::Parallel<FFLAS::CuttingStrategy::RNSModulus, grain> H(GrainSize);
+	        FFLAS::fgemv(F, FFLAS::FflasNoTrans, m, lda, F.one, A, lda, X, incX, F.zero, Y,  incY, H);
+	        break;
+          }
+          case 2:{
+	        ParSeqHelper::Parallel<FFLAS::CuttingStrategy::RNSModulus, threads> H(t);
+	        FFLAS::fgemv(F, FFLAS::FflasNoTrans, m, lda, F.one, A, lda, X, incX, F.zero, Y,  incY, H);
+	        break;
+          }
+          case 3:{
+            ParSeqHelper::Compose<ParSeqHelper::Parallel<FFLAS::CuttingStrategy::RNSModulus, grain>, ParSeqHelper::Parallel<rec, StrategyParameter::TwoDAdaptive>> H(GrainSize,t);
+
+            FFLAS::fgemv(F, FFLAS::FflasNoTrans, m, lda, F.one, A, lda, X, incX, F.zero, Y,  incY, H);
+            break;
+          }
+          default:{
+	        FFLAS::fgemv(F, FFLAS::FflasNoTrans, m, lda, F.one, A, lda, X, incX, F.zero, Y,  incY);
+	        break;
+          }
+      }
+      if (i) {chrono.stop(); time+=chrono.realtime();}
+    }else{
+      if (i) chrono.start();
+      FFLAS::fgemv(F, FFLAS::FflasNoTrans, m, lda, F.one, A, lda, X, incX, F.zero, Y,  incY);
+      if (i) {chrono.stop(); time+=chrono.realtime();}
+    }
+  }
+  if(!check_result(F, m, lda,  A,  X, incX,  Y, incY)) pass = false;
+  return pass;
+}
+
+template <class Field, class arg>
+void benchmark_disp(Field& F, bool pass, double& time, size_t iters, int p,  size_t m, size_t k, arg& as){
+  if(pass){
+    std::cout << "Time: " << time / double(iters)
+	      << " Gflops: " << (2.*double(m)/1000.*double(k)/1000.0/1000.0) / time * double(iters);
+    writeCommandString(std::cout, as) << std::endl;
+  }else{
+    std::cout<<"FAILED for "<<typeid(Field).name()<<std::endl;
+    std::cout << "p:=" << p << ';'<<std::endl;
+  }
+}
+
+
+template <class Field, class arg>
+void benchmark_in_Field(Field& F, int p,  size_t m, size_t k, int NBK, int bitsize, uint64_t seed, size_t iters,
+			int t, arg& as, size_t GrainSize){
+  double time=0.0;
+  size_t lda,incX,incY;
+  lda=k;
+  incX=1;
+  incY=1;
+  typename Field::Element_ptr A,X,Y;
+
+    A = FFLAS::fflas_new(F,m,lda);
+    X = FFLAS::fflas_new(F,k,incX);
+    Y = FFLAS::fflas_new(F,m,incY);
+
+    genData(F, A, X, Y, m, k, incX, incY, lda, NBK, bitsize, seed);
+
+    bool pass=benchmark_with_timer( F, p, A, X, Y, m, k, incX, incY, lda, iters, t, time, GrainSize);
+
+    benchmark_disp(F, pass, time, iters, p, m, k, as);
+
+    FFLAS::fflas_delete(A);
+    FFLAS::fflas_delete(X);
+    FFLAS::fflas_delete(Y);
+
+}
+
+template <class Field,  class arg >
+void benchmark_with_field(int p,  size_t m, size_t k, int NBK, int bitsize, uint64_t seed, size_t iters,
+			  int t, arg& as, size_t GrainSize){
+  Field F;
+  //static assert to raise compile time error for Non ZRing without providing a characteristic
+  static_assert(!need_field_characteristic<Field>::value,
+		"A field characteristic should be provided for Non ZRing data type !");
+  //static assert to raise compile time error for ZRing with either float or double that could lead to inconsistent result
+  static_assert(compatible_data_type<Field>::value,
+		"The provided data type for ZRing is not compatible for the desired operation and could lead to inconsistent result !");
+
+  benchmark_in_Field(F, p,  m, k, NBK, bitsize, seed, iters, t, as, GrainSize);
+
+}
+
+template <class Field, class arg>
+void benchmark_with_field(const Givaro::Integer& q, int p,  size_t m, size_t k,
+			  int NBK, int bitsize, uint64_t seed, size_t iters, int t,
+			  arg& as, size_t GrainSize){
+    Field  F(q);
+    benchmark_in_Field(F, p,  m, k, NBK, bitsize, seed, iters, t, as, GrainSize);
+}
+
+int main(int argc, char** argv) {
+
+  int p=0;
+
+  size_t iters = 3;
+  Givaro::Integer q = 131071;
+  size_t m = 4000;
+  size_t k = 4000;
+
+  uint64_t seed = getSeed();
+  int t;
+  PAR_BLOCK { t = NUM_THREADS; }
+  int NBK = -1;
+  int b=100;
+  size_t GrainSize = 64;
+
+  Argument as[] = {
+    { 'q', "-q Q", "Set the field characteristic (-1 for random).",                 TYPE_INTEGER , &q },
+    { 'b', "-b B", "Set the bitsize of input.",                                     TYPE_INT , &b },
+    { 'p', "-p P", "0 for sequential, 1 for <RNSModulus,Grain>, 2 for <RNSModulus,Thread>, 3 for Compose<<RNSModulus, grain>, <Recursive, TwoDAdaptive>>.",
+                                                                                    TYPE_INT , &p },
+    { 'm', "-m M", "Set the dimension m of the matrix.",                            TYPE_INT , &m },
+    { 'k', "-k K", "Set the dimension k of the matrix.",                            TYPE_INT , &k },
+    { 't', "-t T", "number of virtual threads to drive the partition.",             TYPE_INT , &t },
+    { 'N', "-n N", "number of numa blocks per dimension for the numa placement",    TYPE_INT , &NBK },
+    { 'i', "-i R", "Set number of repetitions.",                                    TYPE_INT , &iters },
+    { 's', "-s S", "Sets seed.",                            				        TYPE_INT , &seed },
+    { 'g', "-g G", "Sets GrainSize.",                            			        TYPE_INT , &GrainSize },
+    END_OF_ARGUMENTS
+  };
+
+  parseArguments(argc,argv,as);
+
+  if (NBK==-1) NBK = t;
+
+  PAR_BLOCK {
+    //benchmark_with_field<Givaro::Modular<Givaro::Integer>>( p,  m, k, NBK, b, seed, iters, t, as, GrainSize);
+    benchmark_with_field<Givaro::ZRing<Givaro::Integer>>( p,  m, k, NBK, b, seed, iters, t, as, GrainSize);
+  }
+
+
+  return 0;
+}
+/* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
+// vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
diff --git a/benchmarks/benchmark-fgemv.C b/benchmarks/benchmark-fgemv.C
index 5fb04802b..544c5c3c3 100644
--- a/benchmarks/benchmark-fgemv.C
+++ b/benchmarks/benchmark-fgemv.C
@@ -151,12 +151,12 @@ bool benchmark_with_timer(Field& F, int p, Matrix& A, Vector& X, Vector& Y, size
       FFLAS::fgemv(F, FFLAS::FflasNoTrans, m, lda, F.one, A, lda, X, incX, F.zero, Y,  incY);
       if (i) {chrono.stop(); time+=chrono.realtime();}
     }
-
+/*
     if(!check_result(F, m, lda,  A,  X, incX,  Y, incY)){
       pass = false;
       break;
     }
-
+*/
   }
   return pass;
 }
@@ -240,7 +240,7 @@ int main(int argc, char** argv) {
   int t;
   PAR_BLOCK { t = NUM_THREADS; }
   int NBK = -1;
-  uint64_t b=0;
+  uint64_t b=100;
   size_t GrainSize = 64;
 
   Argument as[] = {
@@ -268,15 +268,15 @@ int main(int argc, char** argv) {
     }
   }else{
     PAR_BLOCK {
-      //benchmark_with_field<Givaro::Modular<float>>(q, p,  m, k, NBK, b, seed, iters, t, as);
-      //benchmark_with_field<Givaro::Modular<double>>(q, p,  m, k, NBK, b, seed, iters, t, as);
-      //benchmark_with_field<Givaro::Modular<int32_t>>(q, p,  m, k, NBK, b, seed, iters, t, as);
+      //benchmark_with_field<Givaro::Modular<float>>(q, p,  m, k, NBK, b, seed, iters, t, as, GrainSize);
+      //benchmark_with_field<Givaro::Modular<double>>(q, p,  m, k, NBK, b, seed, iters, t, as, GrainSize);
+      //benchmark_with_field<Givaro::Modular<int32_t>>(q, p,  m, k, NBK, b, seed, iters, t, as, GrainSize);
 
-      //benchmark_with_field<Givaro::Modular<Givaro::Integer>>(q, p,  m, k, NBK, b, seed, iters, t, as);
+      //benchmark_with_field<Givaro::Modular<Givaro::Integer>>(q, p,  m, k, NBK, b, seed, iters, t, as, GrainSize);
 
-      //benchmark_with_field<Givaro::ModularBalanced<float>>(q, p,  m, k, NBK, b, seed, iters, t, as);
+      //benchmark_with_field<Givaro::ModularBalanced<float>>(q, p,  m, k, NBK, b, seed, iters, t, as, GrainSize);
       benchmark_with_field<Givaro::ModularBalanced<double>>(q, p,  m, k, NBK, b, seed, iters, t, as, GrainSize);
-      //benchmark_with_field<Givaro::ModularBalanced<int32_t>>(q, p,  m, k, NBK, b, seed, iters, t, as);
+      //benchmark_with_field<Givaro::ModularBalanced<int32_t>>(q, p,  m, k, NBK, b, seed, iters, t, as, GrainSize);
     }
   }
 
diff --git a/configure.ac b/configure.ac
index 04df5b84f..cd6a8ee62 100644
--- a/configure.ac
+++ b/configure.ac
@@ -158,7 +158,7 @@ AC_PROG_LIBTOOL
 AC_PROG_EGREP
 AC_PROG_SED
 # newer libtool...
-LT_PREREQ([2.4.3])
+LT_PREREQ([2.4.2])
 LT_INIT
 
 
diff --git a/fflas-ffpack/fflas/fflas_bounds.inl b/fflas-ffpack/fflas/fflas_bounds.inl
index 896cd2484..4cea3e8fa 100644
--- a/fflas-ffpack/fflas/fflas_bounds.inl
+++ b/fflas-ffpack/fflas/fflas_bounds.inl
@@ -37,6 +37,10 @@
 #include <givaro/modular.h>
 #include <givaro/modular-balanced.h>
 
+#ifdef PROFILE_FGEMM_MP
+#include "fflas-ffpack/utils/timer.h"
+#endif
+
 namespace FFLAS { namespace Protected {
 
     template <class Field>
@@ -115,11 +119,27 @@ namespace FFLAS {
     inline Givaro::Integer
     InfNorm (const size_t M, const size_t N, const Givaro::Integer* A, const size_t lda){
         Givaro::Integer max = 0;
-        for (size_t i=0; i<M; ++i)
-            for (size_t j=0; j<N; ++j) {
-                const Givaro::Integer & x(A[i*lda+j]);
-                if (Givaro::absCompare(x,max)>0) max = x;
-            }
+        std::vector<Givaro::Integer> vmax(M,0);
+        auto sp=SPLITTER(NUM_THREADS,FFLAS::CuttingStrategy::Row,FFLAS::StrategyParameter::Threads);
+        SYNCH_GROUP({
+              FORBLOCK1D(iter, M, sp,
+                          TASK(MODE(CONSTREFERENCE(A,max,vmax) ),
+                                {
+                                for(auto i=iter.begin(); i!=iter.end(); ++i)
+                                {
+                                    for (size_t j=0; j<N; ++j) {
+                                        const Givaro::Integer & x(A[i*lda+j]);
+                                        if (Givaro::absCompare(x,vmax[i])>0){ vmax[i] = x;}
+                                    }
+
+                                }
+                                })
+                          );
+        });
+        max=vmax[0];
+        for (size_t i=0; i<M; ++i){
+            if (Givaro::absCompare(vmax[i],max)>0){ max = vmax[i];}
+        }
         return abs(max);
     }
 
diff --git a/fflas-ffpack/fflas/fflas_fgemm.inl b/fflas-ffpack/fflas/fflas_fgemm.inl
index a87f8b21e..ac7ac1ed9 100644
--- a/fflas-ffpack/fflas/fflas_fgemm.inl
+++ b/fflas-ffpack/fflas/fflas_fgemm.inl
@@ -410,7 +410,9 @@ namespace FFLAS {
             else if (!std::is_same<Field,Givaro::ModularBalanced<float> >::value){
                 if (F.characteristic() < DOUBLE_TO_FLOAT_CROSSOVER)
                     return Protected::fgemm_convert<Givaro::ModularBalanced<float>,Field>(F,ta,tb,m,n,k,alpha,A,lda,B,ldb,beta,C,ldc,H);
-                else if (!std::is_same<Field,Givaro::ModularBalanced<double> >::value && 16*F.cardinality() < Givaro::ModularBalanced<double>::maxCardinality())
+                else if (!std::is_same<Field,Givaro::ModularBalanced<double> >::value &&
+			 !std::is_same<Field,Givaro::ModularBalanced<double> >::value &&
+			 16*F.cardinality() < Givaro::ModularBalanced<double>::maxCardinality())
                     return Protected::fgemm_convert<Givaro::ModularBalanced<double>,Field>(F,ta,tb,m,n,k,alpha,A,lda,B,ldb,beta,C,ldc,H);
             }
         }
diff --git a/fflas-ffpack/fflas/fflas_fgemm/fgemm_classical_mp.inl b/fflas-ffpack/fflas/fflas_fgemm/fgemm_classical_mp.inl
index 4fef8e4e4..b772783d5 100644
--- a/fflas-ffpack/fflas/fflas_fgemm/fgemm_classical_mp.inl
+++ b/fflas-ffpack/fflas/fflas_fgemm/fgemm_classical_mp.inl
@@ -192,7 +192,7 @@ namespace FFLAS {
     }
 
     // fgemm for RnsInteger: handle the moduli in parallel
-    template<typename RNS, typename ParSeqTrait>
+  template<typename RNS, class param, typename ParSeqTrait>
     inline  typename FFPACK::RNSInteger<RNS>::Element_ptr
     fgemm (const FFPACK::RNSInteger<RNS> &F,
            const FFLAS_TRANSPOSE ta,
@@ -203,28 +203,110 @@ namespace FFLAS {
            typename FFPACK::RNSInteger<RNS>::ConstElement_ptr Bd, const size_t ldb,
            const typename FFPACK::RNSInteger<RNS>::Element beta,
            typename FFPACK::RNSInteger<RNS>::Element_ptr Cd, const size_t ldc,
-           MMHelper<FFPACK::RNSInteger<RNS>, MMHelperAlgo::Classic, ModeCategories::DefaultTag, ParSeqHelper::Compose<ParSeqHelper::Parallel<CuttingStrategy::RNSModulus, StrategyParameter::Threads>, ParSeqTrait> > & H)
+           MMHelper<FFPACK::RNSInteger<RNS>, MMHelperAlgo::Classic, ModeCategories::DefaultTag, ParSeqHelper::Compose<ParSeqHelper::Parallel<CuttingStrategy::RNSModulus, param>, ParSeqTrait> > & H)
     {
 #ifdef PROFILE_FGEMM_MP
         Givaro::Timer t;t.start();
 #endif
         size_t rns_size = F.size();
-        typedef MMHelper<typename RNS::ModField, MMHelperAlgo::Winograd, typename ModeTraits<typename RNS::ModField>::value, ParSeqTrait> SubHelper;
-        FORBLOCK1D(iter, rns_size, H.parseq.first_component(),
+        typedef MMHelper<typename RNS::ModField, MMHelperAlgo::Auto, typename ModeTraits<typename RNS::ModField>::value, ParSeqTrait> SubHelper;
+
+        SYNCH_GROUP({
+              FORBLOCK1D(iter, rns_size, H.parseq.first_component(),
+		            TASK(MODE(CONSTREFERENCE(F,H)),
+                                {
+                                for(auto i=iter.begin(); i!=iter.end(); ++i)
+                                {
+                                SubHelper Hsub (F.rns()._field_rns[i], H.recLevel, H.parseq.second_component());
+                                FFLAS::fgemm(F.rns()._field_rns[i],ta,tb,
+                                             m, n, k, alpha._ptr[i*alpha._stride],
+                                             Ad._ptr+i*Ad._stride, lda,
+                                             Bd._ptr+i*Bd._stride, ldb,
+                                             beta._ptr[i*beta._stride],
+                                             Cd._ptr+i*Cd._stride, ldc
+                                             , Hsub
+                                             );
+                                }
+                                })
+                          );
+
+        });
+
+#ifdef PROFILE_FGEMM_MP
+        t.stop();
+        std::cerr<<"=========================================="<<std::endl
+        <<"Pointwise fgemm : "<<t.realtime()<<" ("<<F.size()<<") moduli "<<std::endl
+        <<"=========================================="<<std::endl;
+#endif
+        return Cd;
+    }
+
+
+    // fgemm for RnsInteger default parallel version
+    template<typename RNS, typename Cut, typename Param>
+    inline  typename FFPACK::RNSInteger<RNS>::Element_ptr
+    fgemm (const FFPACK::RNSInteger<RNS> &F,
+           const FFLAS_TRANSPOSE ta,
+           const FFLAS_TRANSPOSE tb,
+           const size_t m, const size_t n,const size_t k,
+           const typename FFPACK::RNSInteger<RNS>::Element alpha,
+           typename FFPACK::RNSInteger<RNS>::ConstElement_ptr Ad, const size_t lda,
+           typename FFPACK::RNSInteger<RNS>::ConstElement_ptr Bd, const size_t ldb,
+           const typename FFPACK::RNSInteger<RNS>::Element beta,
+           typename FFPACK::RNSInteger<RNS>::Element_ptr Cd, const size_t ldc,
+           MMHelper<FFPACK::RNSInteger<RNS>, MMHelperAlgo::Classic, ModeCategories::DefaultTag, ParSeqHelper::Parallel<Cut,Param> > & H)
+    {
+        // compute each fgemm componentwise
+        size_t rns_size = F.size();
+        size_t nt = H.parseq.numthreads();
+        size_t loop_nt = std::min (rns_size, nt);
+        size_t iter_nt = nt / loop_nt;
+        size_t leftover_nt = nt % loop_nt;
+        ParSeqHelper::Parallel<CuttingStrategy::RNSModulus, StrategyParameter::Threads> Hloop (loop_nt);
+#ifdef PROFILE_FGEMM_MP
+        Givaro::Timer t;t.start();
+#endif
+        typedef MMHelper<typename RNS::ModField,
+                MMHelperAlgo::Winograd,
+                typename ModeTraits<typename RNS::ModField>::value,
+                ParSeqHelper::Parallel<Cut,Param> > SubPar;
+
+        typedef MMHelper<typename RNS::ModField,
+                MMHelperAlgo::Winograd,
+                typename ModeTraits<typename RNS::ModField>::value,
+                ParSeqHelper::Sequential> SubSeq;
+
+        FORBLOCK1D(iter, rns_size, Hloop,
                    TASK(MODE(CONSTREFERENCE(F,H)),
                         {
                         for(auto i=iter.begin(); i!=iter.end(); ++i)
                         {
-                        SubHelper Hsub (F.rns()._field_rns[i], H.recLevel, H.parseq.second_component());
-                        FFLAS::fgemm(F.rns()._field_rns[i],ta,tb,
-                                     m, n, k, alpha._ptr[i*alpha._stride],
-                                     Ad._ptr+i*Ad._stride, lda,
-                                     Bd._ptr+i*Bd._stride, ldb,
-                                     beta._ptr[i*beta._stride],
-                                     Cd._ptr+i*Cd._stride, ldc, Hsub);
+                        size_t fgemm_nt = iter_nt;
+                        if (i < leftover_nt)
+                        fgemm_nt++;
+                        if (fgemm_nt>1) // Running a parallel fgemm
+                        {
+                        SubPar H2(F.rns()._field_rns[i], H.recLevel, ParSeqHelper::Parallel<Cut,Param>(fgemm_nt));
+                        fgemm(F.rns()._field_rns[i], ta, tb, m, n, k,
+                              alpha._ptr[i*alpha._stride], Ad._ptr+i*Ad._stride,
+                              lda, Bd._ptr+i*Bd._stride, ldb,
+                              beta._ptr[i*beta._stride], Cd._ptr+i*Cd._stride,
+                              ldc, H2);
+                        }
+                        else // Running a sequential fgemm
+                        {
+                        SubSeq H2(F.rns()._field_rns[i], H.recLevel, ParSeqHelper::Sequential());
+                        fgemm(F.rns()._field_rns[i], ta, tb, m, n, k,
+                              alpha._ptr[i*alpha._stride], Ad._ptr+i*Ad._stride,
+                              lda, Bd._ptr+i*Bd._stride, ldb,
+                              beta._ptr[i*beta._stride], Cd._ptr+i*Cd._stride,
+                              ldc, H2);
                         }
-                        })
-                  );
+                        }
+                        }); // TASK
+        ); // FLORBLOCK1D
+
+
 #ifdef PROFILE_FGEMM_MP
         t.stop();
         std::cerr<<"=========================================="<<std::endl
@@ -234,7 +316,8 @@ namespace FFLAS {
         return Cd;
     }
 
-    // fgemm for RnsInteger default parallel version
+
+    // Specialization of fgemm for RnsInteger parallel version for CuttingStrategy::RNSModulus
     template<typename RNS, typename Cut, typename Param>
     inline  typename FFPACK::RNSInteger<RNS>::Element_ptr
     fgemm (const FFPACK::RNSInteger<RNS> &F,
@@ -246,7 +329,7 @@ namespace FFLAS {
            typename FFPACK::RNSInteger<RNS>::ConstElement_ptr Bd, const size_t ldb,
            const typename FFPACK::RNSInteger<RNS>::Element beta,
            typename FFPACK::RNSInteger<RNS>::Element_ptr Cd, const size_t ldc,
-           MMHelper<FFPACK::RNSInteger<RNS>, MMHelperAlgo::Classic, ModeCategories::DefaultTag, ParSeqHelper::Parallel<Cut,Param> > & H)
+           MMHelper<FFPACK::RNSInteger<RNS>, MMHelperAlgo::Classic, ModeCategories::DefaultTag, ParSeqHelper::Parallel<CuttingStrategy::RNSModulus,Param> > & H)
     {
         // compute each fgemm componentwise
         size_t rns_size = F.size();
@@ -298,9 +381,9 @@ namespace FFLAS {
                         }); // TASK
         ); // FLORBLOCK1D
 
+
 #ifdef PROFILE_FGEMM_MP
         t.stop();
-
         std::cerr<<"=========================================="<<std::endl
         <<"Pointwise fgemm : "<<t.realtime()<<" ("<<rns_size<<") moduli "<<std::endl
         <<"=========================================="<<std::endl;
@@ -322,17 +405,19 @@ namespace FFLAS {
            Givaro::Integer* C, const size_t ldc,
            MMHelper<Givaro::ZRing<Givaro::Integer>, MMHelperAlgo::Classic, ModeCategories::ConvertTo<ElementCategories::RNSElementTag>, ParSeq >  & H)
     {
-        //std::cerr<<"Entering fgemm<ZRing<Integer>> ParSeq"<<std::endl;
+      //std::cerr<<"Entering fgemm<ZRing<Integer>> ParSeq"<<std::endl;
 #ifdef PROFILE_FGEMM_MP
         Timer chrono;
         chrono.start();
 #endif
-        if (alpha == 0){
+
+	if (alpha == 0){
             fscalin(F,m,n,beta,C,ldc);
             return C;
         }
 
         if (k==0) return C;
+
         // compute bit size of feasible prime for FFLAS
         size_t _k=k,lk=0;
         while ( _k ) {_k>>=1; ++lk;}
@@ -352,6 +437,7 @@ namespace FFLAS {
 
         mC = 2*uint64_t(k)*H.normA*H.normB*abs(alpha); // need to use 2x bound to reach both positive and negative
 
+
         // A or B is zero, no need to modify C
         if (mC == 0) return C;
 
@@ -381,6 +467,7 @@ namespace FFLAS {
         chrono.start();
 #endif
 
+
         // convert the input matrices to RNS representation
         finit_rns(Zrns,Arowd,Acold,(logA/16)+((logA%16)?1:0),A,lda,Ap);
         finit_rns(Zrns,Browd,Bcold,(logB/16)+((logB%16)?1:0),B,ldb,Bp);
@@ -422,6 +509,7 @@ namespace FFLAS {
         std::cout<<"-------------------------------"<<std::endl;
 #endif
 
+
         return C;
     }
 
@@ -595,6 +683,7 @@ namespace FFLAS {
         Timer chrono;
         chrono.start();
 #endif
+
         if (alpha == 0){
             fscalin(F,m,n,beta,C,ldc);
             return C;
@@ -671,7 +760,6 @@ namespace FFLAS {
         chrono.start();
 #endif
 
-
         // convert the RNS output to integer representation (C=beta.C+ RNS^(-1)(Cp) ) modulo mA
         //fconvert_rns(Zrns,m,n,beta,C,ldc,Cp);
         RNS.convert(m,n,beta,C,ldc,Cp._ptr,Cp._stride,mA);
diff --git a/fflas-ffpack/fflas/fflas_fgemv.inl b/fflas-ffpack/fflas/fflas_fgemv.inl
index b7aba54b1..d0f0ee8ff 100644
--- a/fflas-ffpack/fflas/fflas_fgemv.inl
+++ b/fflas-ffpack/fflas/fflas_fgemv.inl
@@ -33,7 +33,7 @@
 #if defined(__FFLASFFPACK_HAVE_SSE4_1_INSTRUCTIONS) and defined(__x86_64__)
 #include "fflas-ffpack/fflas/fflas_igemm/igemm.h"
 #endif
-
+//!Function converting to float field for the provided field
 namespace FFLAS{ namespace Protected {
     template <typename FloatElement, class Field>
     inline typename Field::Element_ptr
@@ -82,7 +82,7 @@ namespace FFLAS{ namespace Protected {
     }
 }// Protected
 }// FFLAS
-
+//!Convert to either float or double according to field's cardinality
 namespace FFLAS {
     template<class Field>
     inline  typename Field::Element_ptr
@@ -114,7 +114,7 @@ namespace FFLAS {
     // Computes  Y <- alpha.op(A).X + beta.Y
     // A is M*N,
     //---------------------------------------------------------------------
-
+    //! Performs Matrix Vector Multiplication with delayed mod reductions. Ensures result is reduced.
     template<class Field>
     inline typename Field::Element_ptr
     fgemv (const Field& F, const FFLAS_TRANSPOSE ta,
@@ -126,7 +126,6 @@ namespace FFLAS {
            typename Field::Element_ptr Y, const size_t incY,
            MMHelper<Field, MMHelperAlgo::Classic, ModeCategories::DelayedTag> & H)
     {
-
         if (!M) {return Y;}
         size_t Ydim = (ta == FflasNoTrans)?M:N;
         size_t Xdim = (ta == FflasNoTrans)?N:M;
@@ -404,6 +403,8 @@ namespace FFLAS{
 #endif
         return Y;
     }
+
+    //specialization for ZRing<double>
     inline Givaro::DoubleDomain::Element_ptr
     fgemv (const Givaro::DoubleDomain& F, const FFLAS_TRANSPOSE ta,
            const size_t M, const size_t N,
@@ -442,6 +443,7 @@ namespace FFLAS{
         return fgemv(F, ta, M, N, alpha, A, lda, X, incX, beta, Y, incY, Hb);
     }
 
+    ////specialization for ZRing<float>
     inline Givaro::FloatDomain::Element_ptr
     fgemv (const Givaro::FloatDomain& F, const FFLAS_TRANSPOSE ta,
            const size_t M, const size_t N,
@@ -463,6 +465,45 @@ namespace FFLAS{
         return Y;
     }
 
+    //specialization for Givaro::ZRing<Givaro::Integer> with ParSeqHelper::Compose
+    template<class... ComposeArgs>
+    Givaro::Integer*
+    fgemv(const Givaro::ZRing<Givaro::Integer>& F,
+           const FFLAS_TRANSPOSE ta,
+           const size_t m,
+           const size_t n,
+           const Givaro::Integer alpha,
+           const Givaro::Integer* A, const size_t lda,
+           const Givaro::Integer* X, const size_t incX,
+           const Givaro::Integer beta,
+           Givaro::Integer* Y, const size_t incY,
+           ParSeqHelper::Compose<ComposeArgs...>& parH){
+        MMHelper<Givaro::ZRing<Givaro::Integer>, MMHelperAlgo::Auto, FFLAS::ModeTraits<Givaro::ZRing<Givaro::Integer>>::value, ParSeqHelper::Compose<ComposeArgs...>> pH (F,m,n,1,parH);
+        fgemv(F, ta, m, n, alpha, A, lda, X, incX, beta, Y, incY, pH);
+        return Y;
+    }
+
+    //specialization for Givaro::Modular<Givaro::Integer> with ParSeqHelper::Compose
+    template<class... ComposeArgs>
+    Givaro::Integer*
+    fgemv(const Givaro::Modular<Givaro::Integer>& F,
+           const FFLAS_TRANSPOSE ta,
+           const size_t m,
+           const size_t n,
+           const Givaro::Integer alpha,
+           const Givaro::Integer* A, const size_t lda,
+           const Givaro::Integer* X, const size_t incX,
+           const Givaro::Integer beta,
+           Givaro::Integer* Y, const size_t incY,
+           ParSeqHelper::Compose<ComposeArgs...>& parH){
+        MMHelper<Givaro::Modular<Givaro::Integer>, MMHelperAlgo::Auto, FFLAS::ModeTraits<Givaro::Modular<Givaro::Integer>>::value, ParSeqHelper::Compose<ComposeArgs...>> pH (F,m,n,1,parH);
+        fgemv(F, ta, m, n, alpha, A, lda, X, incX, beta, Y, incY, pH);
+        return Y;
+    }
+
+
+
+    //Common interface for fgemv with ParSeqHelper::Parallel input parameter in which the corresponding parallel implementation will be called for the given field (ref. pfgemv.inl)
     template<class Field, class Cut, class Param>
     typename Field::Element_ptr
     fgemv(const Field& F,
@@ -475,10 +516,13 @@ namespace FFLAS{
            const typename Field::Element beta,
            typename Field::Element_ptr Y, const size_t incY,
            ParSeqHelper::Parallel<Cut,Param>& parH){
-        MMHelper<Field, MMHelperAlgo::Auto, typename FFLAS::ModeTraits<Field>::value, ParSeqHelper::Parallel<Cut,Param> > pH (F,m,n,1,parH);
-        return fgemv(F, ta, m, n, alpha, A, lda, X, incX, beta, Y, incY, pH);
+        MMHelper<Field, MMHelperAlgo::Classic, typename FFLAS::ModeTraits<Field>::value, ParSeqHelper::Parallel<Cut,Param> > pH (F,m,n,1,parH);
+        fgemv(F, ta, m, n, alpha, A, lda, X, incX, beta, Y, incY, pH);
+        return Y;
     }
 
+
+    //Common interface for fgemv with ParSeqHelper::Sequential input parameter in which the corresponding sequential implementation will be called for the given field type either for common field implementated as above or multiprcesion field ref. fflas_fgemv_mp.inl
     template<class Field>
     typename Field::Element_ptr
     fgemv(const Field& F,
@@ -491,11 +535,13 @@ namespace FFLAS{
            const typename Field::Element beta,
            typename Field::Element_ptr Y, const size_t incY,
            ParSeqHelper::Sequential& seqH ){
-        MMHelper<Field, MMHelperAlgo::Classic, ModeCategories::DefaultTag> pH(F,m,n,1,seqH);
+        MMHelper<Field, MMHelperAlgo::Auto, ModeCategories::DefaultTag> pH(F,m,n,1,seqH);
         return fgemv(F, ta, m, n, alpha, A, lda, X, incX, beta, Y, incY, pH);
     }
+
 }
 
 #endif //  __FFLASFFPACK_fgemv_INL
 /* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
 // vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
+
diff --git a/fflas-ffpack/fflas/fflas_fgemv_mp.inl b/fflas-ffpack/fflas/fflas_fgemv_mp.inl
index 4f7017f40..b04867456 100644
--- a/fflas-ffpack/fflas/fflas_fgemv_mp.inl
+++ b/fflas-ffpack/fflas/fflas_fgemv_mp.inl
@@ -31,7 +31,6 @@
 
 namespace FFLAS {
 
-
     // specialization of the fgemv function for the field RNSInteger<rns_double>
     inline FFPACK::rns_double::Element_ptr
     fgemv (const FFPACK::RNSInteger<FFPACK::rns_double>& F, const FFLAS_TRANSPOSE ta,
@@ -57,7 +56,6 @@ namespace FFLAS {
         return Y;
     }
 
-
     // specialization of the fgemv function for the field RNSIntegerMod<rns_double>
     inline FFPACK::rns_double::Element_ptr
     fgemv (const FFPACK::RNSIntegerMod<FFPACK::rns_double>& F, const FFLAS_TRANSPOSE ta,
@@ -70,7 +68,7 @@ namespace FFLAS {
            MMHelper<FFPACK::RNSIntegerMod<FFPACK::rns_double>, MMHelperAlgo::Classic, ModeCategories::DefaultTag> & H)
     {
         //std::cout<<"HERE 1"<<std::endl;
-        MMHelper<FFPACK::RNSInteger<FFPACK::rns_double>, MMHelperAlgo::Classic, ModeCategories::DefaultTag >  H2;
+        MMHelper<FFPACK::RNSInteger<FFPACK::rns_double>, MMHelperAlgo::Classic, ModeCategories::DefaultTag >  H2(H);
         //std::cout<<"HERE 2"<<std::endl;
         fgemv(F.delayed(),ta,M,N,alpha,A,lda,X,incX, beta,Y,incY,H2);
         //std::cout<<"HERE 3"<<std::endl;
@@ -79,27 +77,56 @@ namespace FFLAS {
         return Y;
     }
 
-
     // BB hack. might not work.
     // Calling fgemm, TODO: really specialize fgemv
     // specialization of the fgemv function for the field Givaro::ZRing<Givaro::Integer>
+    template<class AlgoT>
     inline Givaro::Integer* fgemv (const Givaro::ZRing<Givaro::Integer>& F,
                                    const FFLAS_TRANSPOSE ta,
                                    const size_t m, const size_t n,
                                    const Givaro::Integer alpha,
-                                   Givaro::Integer* A, const size_t lda,
-                                   Givaro::Integer* X, const size_t ldx,
+                                   const Givaro::Integer* A, const size_t lda, // @fixme Why not originally const?
+                                   const Givaro::Integer* X, const size_t ldx,
                                    Givaro::Integer beta,
                                    Givaro::Integer* Y, const size_t ldy,
-                                   MMHelper<Givaro::ZRing<Givaro::Integer>, MMHelperAlgo::Classic, ModeCategories::ConvertTo<ElementCategories::RNSElementTag> > & H)
-    {
-        MMHelper<Givaro::ZRing<Givaro::Integer>, MMHelperAlgo::Classic, ModeCategories::ConvertTo<ElementCategories::RNSElementTag>, ParSeqHelper::Sequential> H2;
-        fgemm(F,ta,FFLAS::FflasNoTrans, (ta==FFLAS::FflasNoTrans)?m:n, 1,(ta==FFLAS::FflasNoTrans)?n:m, alpha,A,lda,X,ldx,beta,Y,ldy,H2);
-        return Y;
+                                   MMHelper<Givaro::ZRing<Givaro::Integer>, AlgoT, ModeCategories::ConvertTo<ElementCategories::RNSElementTag>, ParSeqHelper::Sequential> & H) {
+      fgemm(F,ta,FFLAS::FflasNoTrans, (ta==FFLAS::FflasNoTrans)?m:n, 1,(ta==FFLAS::FflasNoTrans)?n:m, alpha,A,lda,X,ldx,beta,Y,ldy,H);
+      return Y;
+    }
+  template<class AlgoT, class Param>
+    inline Givaro::Integer* fgemv (const Givaro::ZRing<Givaro::Integer>& F,
+                                   const FFLAS_TRANSPOSE ta,
+                                   const size_t m, const size_t n,
+                                   const Givaro::Integer alpha,
+                                   const Givaro::Integer* A, const size_t lda, // @fixme Why not originally const?
+                                   const Givaro::Integer* X, const size_t ldx,
+                                   Givaro::Integer beta,
+                                   Givaro::Integer* Y, const size_t ldy,
+                                   MMHelper<Givaro::ZRing<Givaro::Integer>, AlgoT, ModeCategories::ConvertTo<ElementCategories::RNSElementTag>, ParSeqHelper::Parallel<CuttingStrategy::RNSModulus,Param>> & H){
+    ParSeqHelper::Compose<ParSeqHelper::Parallel<CuttingStrategy::RNSModulus, Param>,ParSeqHelper::Sequential> CompHelper (H.parseq, ParSeqHelper::Sequential());
+    MMHelper<Givaro::ZRing<Givaro::Integer>, AlgoT, ModeCategories::ConvertTo<ElementCategories::RNSElementTag>, ParSeqHelper::Compose<ParSeqHelper::Parallel<CuttingStrategy::RNSModulus, Param>,ParSeqHelper::Sequential>>  Hc(F,m,1,n, CompHelper);
+    fgemv(F,ta, m, n, alpha,A,lda,X,ldx,beta,Y,ldy,Hc);
+    return Y;
+  }
+  template<class AlgoT, class Param, class ... ComposeArgs>
+    inline Givaro::Integer* fgemv (const Givaro::ZRing<Givaro::Integer>& F,
+                                   const FFLAS_TRANSPOSE ta,
+                                   const size_t m, const size_t n,
+                                   const Givaro::Integer alpha,
+                                   const Givaro::Integer* A, const size_t lda, // @fixme Why not originally const?
+                                   const Givaro::Integer* X, const size_t ldx,
+                                   Givaro::Integer beta,
+                                   Givaro::Integer* Y, const size_t ldy,
+                                   MMHelper<Givaro::ZRing<Givaro::Integer>, AlgoT, ModeCategories::ConvertTo<ElementCategories::RNSElementTag>, ParSeqHelper::Compose<ParSeqHelper::Parallel<CuttingStrategy::RNSModulus, Param>,ComposeArgs...>> & H){
+    fgemm(F,ta,FFLAS::FflasNoTrans, (ta==FFLAS::FflasNoTrans)?m:n, 1,(ta==FFLAS::FflasNoTrans)?n:m, alpha,A,lda,X,ldx,beta,Y,ldy,H);
+
+      return Y;
     }
 
+
     // specialization of the fgemv function for the field Givaro::Modular<Givaro::Integer>
     // Calling fgemm, TODO: really specialize fgemv
+    template<class ParSeq>
     inline Givaro::Integer* fgemv (const Givaro::Modular<Givaro::Integer>& F,
                                    const FFLAS_TRANSPOSE ta,
                                    const size_t m, const size_t n,
@@ -108,16 +135,17 @@ namespace FFLAS {
                                    Givaro::Integer* X, const size_t ldx,
                                    Givaro::Integer beta,
                                    Givaro::Integer* Y, const size_t ldy,
-                                   MMHelper<Givaro::Modular<Givaro::Integer>, MMHelperAlgo::Classic, ModeCategories::ConvertTo<ElementCategories::RNSElementTag> > & H)
+                                   MMHelper<Givaro::Modular<Givaro::Integer>, MMHelperAlgo::Classic, ModeCategories::ConvertTo<ElementCategories::RNSElementTag>, ParSeq> & H)
     {
-        MMHelper<Givaro::Modular<Givaro::Integer>, MMHelperAlgo::Classic, ModeCategories::ConvertTo<ElementCategories::RNSElementTag>, ParSeqHelper::Sequential> H2;
+        MMHelper<Givaro::Modular<Givaro::Integer>, MMHelperAlgo::Classic, ModeCategories::ConvertTo<ElementCategories::RNSElementTag>, ParSeq> H2(H);
         fgemm(F,ta,FFLAS::FflasNoTrans,(ta==FFLAS::FflasNoTrans)?m:n,1,(ta==FFLAS::FflasNoTrans)?n:m,alpha,A,lda,X,ldx,beta,Y,ldy,H2);
         return Y;
     }
 
     // specialization of the fgemv function for the field Givaro::Modular<RecInt::ruint<K>>
     // Calling fgemm, TODO: really specialize fgemv
-    template <size_t K1, size_t K2, class ParSeq>
+    //@QuickFix: This is only the sequential implementation and any call to parallel fgemv for the field Givaro::Modular<RecInt::ruint<K>> will refer to the implementation in the pfgemv.inl file
+    template <size_t K1, size_t K2>
     inline RecInt::ruint<K1>*
     fgemv (const Givaro::Modular<RecInt::ruint<K1>,RecInt::ruint<K2> >& F,
            const FFLAS_TRANSPOSE ta,
@@ -130,8 +158,8 @@ namespace FFLAS {
            MMHelper<Givaro::Modular<RecInt::ruint<K1>,RecInt::ruint<K2> >,
            MMHelperAlgo::Classic,
            ModeCategories::ConvertTo<ElementCategories::RNSElementTag>,
-           ParSeq >  & H) {
-        MMHelper<Givaro::Modular<RecInt::ruint<K1>,RecInt::ruint<K2> >, MMHelperAlgo::Classic, ModeCategories::ConvertTo<ElementCategories::RNSElementTag>, ParSeqHelper::Sequential> H2;
+           ParSeqHelper::Sequential >  & H) {
+        MMHelper<Givaro::Modular<RecInt::ruint<K1>,RecInt::ruint<K2> >, MMHelperAlgo::Classic, ModeCategories::ConvertTo<ElementCategories::RNSElementTag>, ParSeqHelper::Sequential> H2(H);
         fgemm (F,ta,FflasNoTrans,(ta==FFLAS::FflasNoTrans)?m:n,1,(ta==FFLAS::FflasNoTrans)?n:m,alpha,A,lda,X,incx,beta,Y,incy,H2);
         return Y;
     }
diff --git a/fflas-ffpack/field/rns-double.h b/fflas-ffpack/field/rns-double.h
index 8286c8768..69bffb28c 100644
--- a/fflas-ffpack/field/rns-double.h
+++ b/fflas-ffpack/field/rns-double.h
@@ -53,7 +53,7 @@
 namespace FFPACK {
 
     /* Structure that handles rns representation given a bound and bitsize for prime moduli
-     * support sign representation (i.e. the bound must be twice larger then ||A||)
+     * support sign representation (i.e. the bound must be twice larger than ||A||)
      */
     struct rns_double {
         typedef Givaro::Integer integer;
@@ -209,7 +209,7 @@ namespace FFPACK {
             // std::cout<<"t2="<<t2<<std::endl;
             // std::cout<<"t3="<<t3<<std::endl;
         }
-
+//@TODO: Find out which init is called and add detailed timing for that
         // Arns must be an array of m*n*_size
         // abs(||A||) <= maxA
         template<typename T>
diff --git a/fflas-ffpack/field/rns-double.inl b/fflas-ffpack/field/rns-double.inl
index 7713129d3..4c3d103ee 100644
--- a/fflas-ffpack/field/rns-double.inl
+++ b/fflas-ffpack/field/rns-double.inl
@@ -29,13 +29,16 @@
 #define __FFLASFFPACK_field_rns_double_INL
 
 #include "fflas-ffpack/fflas/fflas_freduce.h"
-
+#ifdef PROFILE_FGEMM_MP
+#include "fflas-ffpack/utils/timer.h"
+#endif
 namespace FFPACK {
 
     // Arns must be an array of m*n*_size
     // abs(||A||) < 2^(16k)
     inline void rns_double::init(size_t m, size_t n, double* Arns, size_t rda, const integer* A, size_t lda, size_t k, bool RNS_MAJOR) const
     {
+
         if (k>_ldm){
             FFPACK::failure()(__func__,__FILE__,__LINE__,"rns_double [init] -> rns basis is too small to handle integers with 2^(16*k) values ");
             std::cerr<<"with k="<<k<<" _ldm="<<_ldm<<std::endl;
@@ -45,7 +48,7 @@ namespace FFPACK {
         double *A_beta = FFLAS::fflas_new<double >(mn*k);
         const integer* Aiter=A;
         // split A into A_beta according to a Kronecker transform in base 2^16
-        //		auto sp=SPLITTER(MAX_THREADS,FFLAS::CuttingStrategy::Column,FFLAS::StrategyParameter::Threads);
+        auto sp=SPLITTER(NUM_THREADS,FFLAS::CuttingStrategy::Row,FFLAS::StrategyParameter::Threads);
 
         Givaro::Timer tkr; tkr.start();
         // #ifndef __FFLASFFPACK_SEQUENTIAL
@@ -57,9 +60,10 @@ namespace FFPACK {
         //       TASK(MODE(READ(Aiter[0]) READWRITE(A_beta[0])),
         //for(size_t i=0;i<m;i++)
         //PAR_BLOCK{
-        //			FOR1D(i,m,sp,
-        PARFOR1D(i,m,SPLITTER(NUM_THREADS),
 
+        FOR1D(i,m,sp,
+
+//	    PARFOR1D(i,m,SPLITTER(NUM_THREADS),
                  for(size_t j=0;j<n;j++){
                  size_t idx=j+i*n;
                  const mpz_t*    m0     = reinterpret_cast<const mpz_t*>(Aiter+j+i*lda);
@@ -92,7 +96,8 @@ namespace FFPACK {
 
                  // 	   );
                  }
-        );
+        )
+
 
         tkr.stop();
         //if(m>1 && n>1) std::cerr<<"Kronecker : "<<tkr.realtime()<<std::endl;
@@ -117,6 +122,8 @@ namespace FFPACK {
             cblas_dgemm(CblasRowMajor,CblasNoTrans,CblasTrans,(int)mn,(int)_size,(int)k,1.0,A_beta,(int)k,_crt_in.data(),(int)_ldm,0.,Arns,(int)_size);
 #endif
         }
+
+
         Givaro::Timer tred; tred.start();
 
         reduce(mn,Arns,rda,RNS_MAJOR);
@@ -142,7 +149,10 @@ namespace FFPACK {
                 }
         std::cout<<"RNS freduce ... "<<(ok?"OK":"ERROR")<<std::endl;
 #endif
+
+
         }
+
     }
 
     // Arns must be an array of m*n*_size
@@ -527,10 +537,12 @@ namespace FFPACK {
             // #else
             // 			auto sp=SPLITTER(1);
             // #endif
-            PARFOR1D(i,_size,SPLITTER(NUM_THREADS),
-                     //for(size_t i=0;i<_size;i++)
+
+            auto sp=SPLITTER(NUM_THREADS);
+            FOR1D(i,_size,sp,{
                      FFLAS::freduce (_field_rns[i],n,Arns+i*rda,1);
-                    );
+            });
+
         }
 
     }
diff --git a/fflas-ffpack/field/rns-integer.h b/fflas-ffpack/field/rns-integer.h
index 27782c83f..053a6a144 100644
--- a/fflas-ffpack/field/rns-integer.h
+++ b/fflas-ffpack/field/rns-integer.h
@@ -110,9 +110,9 @@ namespace FFPACK {
             return x;
         }
         Element& init(Element& x, const Givaro::Integer& y) const{
-            init(x);
-            size_t k =(y.bitsize())/16+((y.bitsize())%16?1:0);
-            _rns->init(1,1,x._ptr,x._stride, &y,1,k);
+	        init(x);
+	        size_t k =(y.bitsize())/16+((y.bitsize())%16?1:0);
+	        _rns->init(1,1,x._ptr,x._stride, &y,1,k);
             return x;
         }
         Element& reduce (Element& x, const Element& y) const {return assign (x,y);}
diff --git a/fflas-ffpack/paladin/blockcuts.inl b/fflas-ffpack/paladin/blockcuts.inl
index 96b337c22..1ffd02d09 100644
--- a/fflas-ffpack/paladin/blockcuts.inl
+++ b/fflas-ffpack/paladin/blockcuts.inl
@@ -40,7 +40,7 @@ namespace FFLAS {
         struct Column{};
         struct Block{};
         struct Recursive{};
-        typedef Row RNSModulus;
+        struct RNSModulus{};
     }
 
     namespace StrategyParameter{
@@ -96,10 +96,10 @@ namespace FFLAS {
         struct Compose{
 
             Compose() : _comp1 (), _comp2 () {}
-            Compose(const Compose & other) : _comp1 (other.first_component()), _comp2 (other.second_component()) {}
+	  Compose(const Compose & other) : _comp1 (other.first_component()), _comp2 (other.second_component()) {}
             Compose(const Sequential & S) : _comp1 (1), _comp2 (1) {}
             Compose(size_t th1, size_t th2) : _comp1 (th1), _comp2 (th2) {}
-            Compose(const H1 & o1, const H2 & o2) : _comp1 (o1), _comp2 (o2) {}
+	  Compose(const H1 & o1, const H2 & o2) : _comp1 (o1), _comp2 (o2) {}
 
             H1 first_component () const { return _comp1; }
             H2 second_component () const { return _comp2; }
diff --git a/fflas-ffpack/paladin/pfgemv.inl b/fflas-ffpack/paladin/pfgemv.inl
index c6f0234ad..440a1f650 100644
--- a/fflas-ffpack/paladin/pfgemv.inl
+++ b/fflas-ffpack/paladin/pfgemv.inl
@@ -25,7 +25,7 @@
 namespace FFLAS
 {
 
-
+    // specialization of the fgemv function for the MMHelper with CuttingStrategy::Recursive but templated for all possible field type so that the corresponding templated sequential implementation will be invoked in the parallel code region
     template<class Field, class AlgoT, class FieldTrait>
     typename Field::Element_ptr
     fgemv(const Field& F,
@@ -38,7 +38,6 @@ namespace FFLAS
            const typename Field::Element beta,
            typename Field::Element_ptr Y, const size_t incY,
            MMHelper<Field, AlgoT, FieldTrait, ParSeqHelper::Parallel<CuttingStrategy::Recursive, StrategyParameter::Threads> > & H){
-
         if (H.parseq.numthreads()==1 || m <= 1){
             fgemv(F, ta,  m, n,  alpha, A, lda, X, incX, beta, Y, incY);
 
@@ -75,7 +74,7 @@ namespace FFLAS
         return Y;
     }
 
-
+    // specialization of the fgemv function for the MMHelper with CuttingStrategy::Row but templated for all possible field type so that the corresponding templated sequential implementation will be invoked in the parallel code region
     template<class Field, class AlgoT, class FieldTrait, class Cut>
     typename Field::Element_ptr
     fgemv(const Field& F,
@@ -106,7 +105,6 @@ namespace FFLAS
         return Y;
     }
 
-
 } // FFLAS
 
 /* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 4326820ec..a623fe333 100755
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -67,6 +67,7 @@ BASIC_TESTS =               \
 		test-fgesv             \
 		test-simd \
 		test-fgemv \
+		test-fgemv-rns \
 		test-nullspace \
 		regression-check
 
@@ -162,6 +163,7 @@ regression_check_SOURCES = regression-check.C
 test_solve_SOURCES = test-solve.C
 test_simd_SOURCES = test-simd.C
 test_fgemv_SOURCES = test-fgemv.C
+test_fgemv_rns_SOURCES = test-fgemv-rns.C
 #test_pfgemm_DSL_SOURCES = test-pfgemm-DSL.C
 
 
diff --git a/tests/test-fgemv-rns.C b/tests/test-fgemv-rns.C
new file mode 100644
index 000000000..75c357214
--- /dev/null
+++ b/tests/test-fgemv-rns.C
@@ -0,0 +1,419 @@
+/*
+ * Copyright (C) the FFLAS-FFPACK group
+ * Written by Clément Pernet
+ *            Brice Boyer (briceboyer) <boyer.brice@gmail.com>
+ * This file is Free Software and part of FFLAS-FFPACK.
+ *
+ * ========LICENCE========
+ * This file is part of the library FFLAS-FFPACK.
+ *
+ * FFLAS-FFPACK is free software: you can redistribute it and/or modify
+ * it under the terms of the  GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * This library 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
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this library; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
+ * ========LICENCE========
+ *.
+ */
+
+// #ifndef NEWINO
+// #define NEWWINO
+// #endif
+
+// #define WINOTHRESHOLD 100
+// #define OLD_DYNAMIC_PEELING
+
+
+
+#include "fflas-ffpack/fflas-ffpack-config.h"
+
+#include <iomanip>
+#include <iostream>
+
+#include <givaro/modular.h>
+
+#include <recint/rint.h>
+
+#include "fflas-ffpack/utils/timer.h"
+#include "fflas-ffpack/fflas/fflas.h"
+
+#include "fflas-ffpack/utils/args-parser.h"
+#include "fflas-ffpack/utils/test-utils.h"
+
+using namespace std;
+using namespace FFPACK;
+using namespace FFLAS;
+
+using Givaro::Modular;
+using Givaro::ModularBalanced;
+
+template <typename Field>
+struct rsn_compatible_data_type { static constexpr bool value = false; };
+template <>
+struct rsn_compatible_data_type<Givaro::ZRing<Givaro::Integer>>{ static constexpr bool value = true; };
+template <typename Field>
+struct rsn_compatible_data_type<Givaro::Modular<Field>>{ static constexpr bool value = false; };
+template <typename Field>
+struct rsn_compatible_data_type<Givaro::ModularBalanced<Field>>{ static constexpr bool value = false; };
+
+// checks that D = beta . Y + alpha . A ^ta * X
+template<class Field>
+bool check_MV(const Field                   & F,
+              const typename Field::Element_ptr  Cd, // c0
+              enum FFLAS_TRANSPOSE   & ta,
+              const size_t                    m,
+              const size_t                    k,
+              const typename Field::Element & alpha,
+              const typename Field::Element_ptr  A, size_t lda,
+              const typename Field::Element_ptr  X, size_t incX,
+              const typename Field::Element & beta,
+              const typename Field::Element_ptr Y, size_t incY)
+{
+    bool wrong = false;
+    typename Field::Element_ptr D;
+    if (ta == FflasNoTrans){
+        D = fflas_new(F,m);
+        fassign (F, m, Cd, 1, D, 1);
+        for (size_t i=0; i<m; i++){
+            F.mulin (D[i], beta);
+            typename Field::Element tmp;
+            F.init(tmp);
+            for (size_t j=0; j<k; j++){
+                F.axpyin (tmp, A[i*lda+j], X[j*incX]);
+            }
+            F.axpyin(D[i],alpha,tmp);
+        }
+        wrong = !fequal(F, m, D, 1, Y, incY);
+    } else {
+        D = fflas_new(F,k);
+        fassign (F, k, Cd, 1, D, 1);
+        for (size_t i=0; i<k; i++){
+            F.mulin (D[i], beta);
+            typename Field::Element tmp;
+            F.init(tmp);
+            for (size_t j=0; j<m; j++){
+                F.axpyin (tmp, A[i+j*lda], X[j*incX]);
+            }
+            F.axpyin(D[i],alpha,tmp);
+        }
+        wrong = !fequal(F, k, D, 1, Y, incY);
+    }
+    size_t Ydim = (ta==FflasNoTrans)? m : k;
+    if ( wrong ){
+        size_t canprint = 20 ;
+        std::cerr<<"FAIL"<<std::endl;
+        std::cerr << "alpha   :" << alpha<<", beta   : " << beta << std::endl;
+        std::cerr << "m   :" << m   << ", k   : " << k << std::endl;
+        std::cerr << "ldA :" << lda << ", incX : " << incX << ", incY : " << incY << std::endl;
+        for (size_t i=0; i<Ydim && canprint; ++i){
+            if (!F.areEqual( Y[i*incY], D[i] ) ) {
+                std::cerr<<"Error Y["<<i<<"]="<<Y[i*incY]<<" D["<<i<<"]="<<D[i]<<std::endl;
+                canprint--;
+            }
+        }
+        if (Ydim<80) {
+            for (size_t i=0; i<Ydim ; ++i){
+                if (!F.areEqual( Y[i*incY], D[i] ) )
+                    std::cout << 'X' ;
+                else
+                    std::cout << '.' ;
+            }
+            std::cout << std::endl;
+        }
+    }
+
+    fflas_delete (D);
+
+    return !wrong ;
+
+}
+
+
+template<class Field, class RandIter>
+bool launch_MV(const Field & F,
+               const size_t   m,
+               const size_t   k,
+               const typename Field::Element alpha,
+               const typename Field::Element beta,
+               const size_t lda,
+               enum FFLAS_TRANSPOSE    ta,
+               const size_t incX,
+               const size_t incY,
+               size_t iters,
+               bool par,
+               RandIter& G)
+{
+
+    bool ok = true;
+
+    typedef typename Field::Element_ptr Element_ptr;
+    Element_ptr A ;
+    for(size_t i = 0;i<iters;++i){
+        FFLASFFPACK_check(lda >= k);
+        A = fflas_new (F, m, lda);
+        fzero(F,m,lda,A,lda);
+        RandomMatrix(F, m, k, A, lda, G);
+        size_t Xdim = (ta == FflasNoTrans)? k : m;
+        size_t Ydim = (ta == FflasNoTrans)? m : k;
+        Element_ptr X = fflas_new (F, Xdim, incX);
+        Element_ptr Y = fflas_new (F, Ydim, incY);
+        fzero (F, Xdim, incX, X, incX);
+        fzero (F, Ydim, incY, Y, incY);
+        Element_ptr D = fflas_new (F, Ydim);
+
+        RandomMatrix (F, Xdim, 1, X, incX, G);
+        RandomMatrix (F, Ydim, 1, Y, incY, G);
+        fassign (F, Ydim, Y, incY, D, 1);
+
+        //Y will be modified so keep a copy of Y as Y2
+        Element_ptr Y2 =  fflas_new (F, Ydim, incY);
+        fassign (F, Ydim, Y, incY, Y2, incY);
+
+        fassign (F, Ydim, Y2, incY, D, 1);
+        fassign (F, Ydim, Y2, incY, Y, incY);
+        if (par){
+            {
+                ParSeqHelper::Parallel<CuttingStrategy::RNSModulus, StrategyParameter::Grain >  WH(1);
+
+                PAR_BLOCK{
+                    fgemv(F, ta, m,k,alpha, A,lda, X, incX, beta, Y, incY, WH);
+                }
+            }
+        }else{
+            //MMHelper<Field,MMHelperAlgo::Auto,typename ModeTraits<Field>::value> WH(F,nbw,ParSeqHelper::Sequential());
+            fgemv(F, ta, m, k,alpha, A,lda, X, incX, beta, Y, incY);
+        }
+
+        ok = ok && check_MV(F, D, ta, m, k,alpha, A, lda, X, incX, beta, Y, incY);
+
+
+
+        if (!ok){
+            fflas_delete (A, X, Y, Y2, D);
+            break;
+        }
+
+
+        fassign (F, Ydim, Y2, incY, D, 1);
+        fassign (F, Ydim, Y2, incY, Y, incY);
+        if (par){
+            {
+
+                PAR_BLOCK{
+                    ParSeqHelper::Parallel<CuttingStrategy::RNSModulus, StrategyParameter::Threads >  WH(NUM_THREADS);
+                    fgemv(F, ta, m,k,alpha, A,lda, X, incX, beta, Y, incY, WH);
+                }
+            }
+        }else{
+            //MMHelper<Field,MMHelperAlgo::Auto,typename ModeTraits<Field>::value> WH(F,nbw,ParSeqHelper::Sequential());
+            fgemv(F, ta, m, k,alpha, A,lda, X, incX, beta, Y, incY);
+        }
+
+        ok = ok && check_MV(F, D, ta, m, k,alpha, A, lda, X, incX, beta, Y, incY);
+
+
+
+        if (!ok){
+            fflas_delete (A, X, Y, Y2, D);
+            break;
+        }
+
+        fassign (F, Ydim, Y2, incY, D, 1);
+        fassign (F, Ydim, Y2, incY, Y, incY);
+        if (par){
+            {
+
+                PAR_BLOCK{
+                    ParSeqHelper::Compose<ParSeqHelper::Parallel<FFLAS::CuttingStrategy::RNSModulus, StrategyParameter::Grain>, ParSeqHelper::Parallel<CuttingStrategy::Recursive, StrategyParameter::TwoDAdaptive>>  WH(1,NUM_THREADS);
+                    fgemv(F, ta, m,k,alpha, A,lda, X, incX, beta, Y, incY, WH);
+                }
+            }
+        }else{
+            //MMHelper<Field,MMHelperAlgo::Auto,typename ModeTraits<Field>::value> WH(F,nbw,ParSeqHelper::Sequential());
+            fgemv(F, ta, m, k,alpha, A,lda, X, incX, beta, Y, incY);
+        }
+
+        ok = ok && check_MV(F, D, ta, m, k,alpha, A, lda, X, incX, beta, Y, incY);
+
+
+
+        if (!ok){
+            fflas_delete (A, X, Y, Y2, D);
+            break;
+        }
+
+
+    }
+    return ok ;
+}
+
+
+template<class Field, class RandIter>
+bool launch_MV_dispatch(const Field &F,
+                        const int mm,
+                        const int kk,
+                        const typename Field::Element alpha,
+                        const typename Field::Element beta,
+                        const size_t iters,
+                        const bool par,
+                        RandIter& G)
+{
+    bool ok = true;
+    size_t m,k;
+    size_t lda,incX, incY;
+    size_t ld = 13 ;
+    {
+        //FFLAS_TRANSPOSE ta = FflasNoTrans ;
+        //if (! par) {
+        //if (random()%2) ta = FflasTrans ;
+        //}
+
+        if (mm<0)
+            m = 1+(size_t)random() % -mm;
+        else m = mm;
+        if (kk<0)
+            k = 1+(size_t)random() % -kk;
+        else k = kk;
+
+        lda = k+(size_t)random()%ld;
+        incX = 1+(size_t)random()%ld;
+        incY = 1+(size_t)random()%ld;
+
+        ok = ok && launch_MV (F, m, k, alpha,beta, lda, FflasNoTrans, incX, incY, iters, par, G);
+        ok = ok && launch_MV (F, m, k, alpha,beta, lda, FflasTrans, incX, incY, iters, par, G);
+    }
+    return ok ;
+}
+
+template <class Field>
+bool run_with_field (Givaro::Integer q, uint64_t b, int m, int k, size_t iters, bool par, uint64_t seed){
+    bool ok = true ;
+
+    int nbit=(int)iters;
+
+    while (ok &&  nbit){
+        typedef typename Field::Element Element ;
+        // choose Field
+        Field* F= chooseField<Field>(q,b,seed);
+        if (F==nullptr)
+            return true;
+
+        std::ostringstream oss;
+        F->write(oss);
+        std::cout.fill('.');
+        std::cout<<"Checking ";
+        std::cout.width(50);
+        std::cout<<oss.str();
+        std::cout<<" ... ";
+
+#ifdef __FFLASFFPACK_DEBUG
+        F->write(std::cerr) << std::endl;
+#endif
+        typedef typename Field::Element  Element ;
+        typename Field::RandIter R(*F,seed++);
+        typename Field::NonZeroRandIter NZR(R);
+
+        //size_t k = 0 ;
+        //std::cout << k << "/24" << std::endl; ++k;
+        ok = ok && launch_MV_dispatch<Field>(*F,m,k,F->one,F->zero,iters, par, R);
+        //std::cout << k << "/24" << std::endl; ++k;
+        ok = ok && launch_MV_dispatch<Field>(*F,m,k,F->zero,F->zero,iters, par, R);
+        //std::cout << k << "/24" << std::endl; ++k;
+        ok = ok && launch_MV_dispatch<Field>(*F,m,k,F->mOne,F->zero,iters, par, R);
+        //std::cout << k << "/24" << std::endl; ++k;
+        ok = ok && launch_MV_dispatch<Field>(*F,m,k,F->one ,F->one,iters, par, R);
+        //std::cout << k << "/24" << std::endl; ++k;
+        ok = ok && launch_MV_dispatch<Field>(*F,m,k,F->zero,F->one,iters, par, R);
+        //std::cout << k << "/24" << std::endl; ++k;
+        ok = ok && launch_MV_dispatch<Field>(*F,m,k,F->mOne,F->one,iters, par, R);
+        //std::cout << k << "/24" << std::endl; ++k;
+        ok = ok && launch_MV_dispatch<Field>(*F,m,k,F->one ,F->mOne,iters, par, R);
+        //std::cout << k << "/24" << std::endl; ++k;
+        ok = ok && launch_MV_dispatch<Field>(*F,m,k,F->zero,F->mOne,iters, par, R);
+        //std::cout << k << "/24" << std::endl; ++k;
+        ok = ok && launch_MV_dispatch<Field>(*F,m,k,F->mOne,F->mOne,iters, par, R);
+        //std::cout << k << "/24" << std::endl; ++k;
+
+        Element alpha,beta ;
+        NZR.random(alpha);
+        ok = ok && launch_MV_dispatch<Field>(*F,m,k,F->one ,alpha,iters, par, R);
+        //std::cout << k << "/24" << std::endl; ++k;
+        ok = ok && launch_MV_dispatch<Field>(*F,m,k,F->zero,alpha,iters, par, R);
+        //std::cout << k << "/24" << std::endl; ++k;
+        ok = ok && launch_MV_dispatch<Field>(*F,m,k,F->mOne,alpha,iters, par, R);
+        //std::cout << k << "/24" << std::endl; ++k;
+        ok = ok && launch_MV_dispatch<Field>(*F,m,k,alpha,F->one ,iters, par, R);
+        //std::cout << k << "/24" << std::endl; ++k;
+        ok = ok && launch_MV_dispatch<Field>(*F,m,k,alpha,F->zero,iters, par, R);
+        //std::cout << k << "/24" << std::endl; ++k;
+        ok = ok && launch_MV_dispatch<Field>(*F,m,k,alpha,F->mOne,iters, par, R);
+        //std::cout << k << "/24" << std::endl; ++k;
+
+        for (size_t j = 0 ; j < 3 ; ++j) {
+            R.random(alpha);
+            R.random(beta);
+            ok = ok && launch_MV_dispatch<Field>(*F,m,k,alpha,beta,iters, par, R);
+            //std::cout << k << "/24" << std::endl; ++k;
+        }
+        //std::cout<<std::endl;
+        nbit--;
+        if ( !ok )
+            //std::cout << "\033[1;31mFAILED\033[0m "<<std::endl;
+            std::cout << "FAILED "<<std::endl;
+        else
+            //std::cout << "\033[1;32mPASSED\033[0m "<<std::endl;
+            std::cout << "PASSED "<<std::endl;
+        delete F;
+    }
+    return ok;
+}
+int main(int argc, char** argv)
+{
+    std::cout<<setprecision(17);
+    std::cerr<<setprecision(17);
+
+    uint64_t seed = getSeed();
+    size_t iters = 3 ;
+    Givaro::Integer q = -1 ;
+    uint64_t b = 100 ;
+    int m = -50 ;
+    int k = -50 ;
+    int nbw = -1 ;
+    bool loop = false;
+    bool p = false;
+    Argument as[] = {
+        { 'q', "-q Q", "Set the field characteristic (-1 for random).",         TYPE_INTEGER , &q },
+        { 'b', "-b B", "Set the bitsize of the random characteristic.",         TYPE_INT , &b },
+        { 'm', "-m M", "Set the dimension of the matrix (negative values, mean, any random value between 0 and |n|).",      TYPE_INT , &m },
+        { 'k', "-k K", "Set the dimension of the matrix (negative values, mean, any random value between 0 and |k|).",      TYPE_INT , &k },
+        { 'w', "-w N", "Set the number of winograd levels (-1 for random).",    TYPE_INT , &nbw },
+        { 'i', "-i R", "Set number of repetitions.",            TYPE_INT , &iters },
+        { 'l', "-l Y/N", "run the test in an infinte loop.", TYPE_BOOL , &loop },
+        { 'p', "-p Y/N", "run the parallel fgemv.", TYPE_BOOL , &p },
+        { 's', "-s seed", "Set seed for the random generator", TYPE_UINT64, &seed },
+        END_OF_ARGUMENTS
+    };
+
+    parseArguments(argc,argv,as);
+
+    bool ok = true;
+    srand(seed);
+    do{
+        ok = ok && run_with_field<Givaro::ZRing<Givaro::Integer> >(0,(b?b:512_ui64),m,k,iters,p, seed);
+    } while (loop && ok);
+
+
+
+
+    return !ok ;
+}
+/* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
+// vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
diff --git a/tests/test-fgemv.C b/tests/test-fgemv.C
index 7ddf8a14b..6a250f7f6 100644
--- a/tests/test-fgemv.C
+++ b/tests/test-fgemv.C
@@ -55,6 +55,14 @@ using namespace FFLAS;
 using Givaro::Modular;
 using Givaro::ModularBalanced;
 
+template <typename Field>
+struct rsn_compatible_data_type { static constexpr bool value = false; };
+template <>
+struct rsn_compatible_data_type<Givaro::ZRing<Givaro::Integer>>{ static constexpr bool value = true; };
+template <typename Field>
+struct rsn_compatible_data_type<Givaro::Modular<Field>>{ static constexpr bool value = false; };
+template <typename Field>
+struct rsn_compatible_data_type<Givaro::ModularBalanced<Field>>{ static constexpr bool value = false; };
 
 // checks that D = beta . Y + alpha . A ^ta * X
 template<class Field>
@@ -243,8 +251,6 @@ bool launch_MV(const Field & F,
             break;
         }
 
-
-
     }
     return ok ;
 }
@@ -409,7 +415,6 @@ int main(int argc, char** argv)
         ok = ok && run_with_field<ModularBalanced<int32_t> >(q,b,m,k,iters,p, seed);
         ok = ok && run_with_field<Modular<int64_t> >(q,b,m,k,iters, p, seed);
         ok = ok && run_with_field<ModularBalanced<int64_t> >(q,b,m,k,iters, p, seed);
-
         ok = ok && run_with_field<Modular<RecInt::rint<8> > >(q,b?b:127_ui64,m,k,iters, p, seed);
         ok = ok && run_with_field<Modular<RecInt::ruint<7>,RecInt::ruint<8> > >(q,b?b:127_ui64,m,k,iters, p, seed);
         ok = ok && run_with_field<Modular<Givaro::Integer> >(q,(b?b:512_ui64),m,k,iters,p, seed);