Changeset 3

Timestamp:

03/11/07 15:11:02 (2 years ago)

Author:

pernet

Message:

New algorithm for charpoly: CharpolyArithProg?

- Complete the unfinished preconditionning phase (using complementary subspace comp. and rec. call)
- Exception management for the Las Vegas failure
- set as default for charpoly

Location:

include/fflas-ffpack

Files:

• ffpack.h (modified) (7 diffs)
• ffpack_charpoly.inl (modified) (16 diffs)
• ffpack_frobenius.inl (modified) (16 diffs)
• ffpack_krylovelim.inl (modified) (1 diff)

include/fflas-ffpack/ffpack.h

@@ -48 +48 @@
                                    FfpackKGFast=4,
                                    FfpackDanilevski=5,
-                                   FfpackKGFastG=6
+                                   FfpackArithProg=6,
+                                   FfpackKGFastG=7
         };
+        
+        class CharpolyFailed{};
         
         enum FFPACK_MINPOLY_TAG { FfpackDense=1,
@@ -451 +454 @@
         CharPoly( const Field& F, std::list<Polynomial>& charp, const size_t N,
                   typename Field::Element * A, const size_t lda,
-                  const enum FFPACK_CHARPOLY_TAG CharpTag= FfpackLUK);
+                  const enum FFPACK_CHARPOLY_TAG CharpTag= FfpackArithProg);
         
         //---------------------------------------------------------------------
@@ -469 +472 @@
 
         // Solve L X = B or X L = B in place
-        // L is M*M if Side == FflasLeft and N*N if Side == FflasRigt, B is M*N.
+        // L is M*M if Side == FflasLeft and N*N if Side == FflasRight, B is M*N.
         // Only the R non trivial column of L are stored in the M*R matrix L
         // Requirement :  so that L could  be expanded in-place
@@ -483 +486 @@
                 F.init(one, 1.0);
                 F.init(zero, 0.0);
+                size_t LM = (Side == FflasRight)?N:M;
                 for (int i=R-1; i>=0; --i){
                         if (  Q[i] > (size_t) i){
                                 //for (size_t j=0; j<=Q[i]; ++j)
                                 //F.init( *(L+Q[i]+j*ldl), 0 );
-                                fcopy( F, M-Q[i]-1, L+Q[i]*(ldl+1)+ldl,ldl, L+(Q[i]+1)*ldl+i, ldl );
-                                for ( size_t j=Q[i]*ldl; j<M*ldl; j+=ldl)
+                                //cerr<<"1 deplacement "<<i<<"<-->"<<Q[i]<<endl;
+                                fcopy( F, LM-Q[i]-1, L+Q[i]*(ldl+1)+ldl,ldl, L+(Q[i]+1)*ldl+i, ldl );
+                                for ( size_t j=Q[i]*ldl; j<LM*ldl; j+=ldl)
                                         F.assign( *(L+i+j), zero );
                         }
                 }
                 ftrsm( F, Side, FflasLower, FflasNoTrans, FflasUnit, M, N, one, L, ldl , B, ldb);
-                
+                //write_field(F,cerr<<"dans solveLB "<<endl,L,N,N,ldl);
                 // Undo the permutation of L
                 for (size_t i=0; i<R; ++i){
@@ -499 +504 @@
                                 //for (size_t j=0; j<=Q[i]; ++j)
                                 //F.init( *(L+Q[i]+j*ldl), 0 );
-                                fcopy( F, M-Q[i]-1, L+(Q[i]+1)*ldl+i, ldl, L+Q[i]*(ldl+1)+ldl,ldl );
-                                for ( size_t j=Q[i]*ldl; j<M*ldl; j+=ldl)
+                                fcopy( F, LM-Q[i]-1, L+(Q[i]+1)*ldl+i, ldl, L+Q[i]*(ldl+1)+ldl,ldl );
+                                for ( size_t j=Q[i]*ldl; j<LM*ldl; j+=ldl)
                                         F.assign( *(L+Q[i]+j), zero );
                         }
@@ -547 +552 @@
                         int j=R-1;
                         while ( j >=0 ) {
+                                //cerr<<"j="<<j<<endl;
                                 k = ib = Q[j];
                                 while ( (j>=0) &&  (Q[j] == k)  ) {--k;--j;}
                                 Ldim = ib-k;
+                                //cerr<<"Ldim, ib, k, N = "<<Ldim<<" "<<ib<<" "<<k<<" "<<N<<endl;
                                 Lcurr = L + j+1 + (k+1)*ldl;
                                 Bcurr = B + ib;
                                 Rcurr = Lcurr + Ldim*ldl;
-                                fgemm( F, FflasNoTrans, FflasNoTrans, M, N-ib-1, Ldim, Mone, Bcurr, ldb, Rcurr, ldl,  one, Bcurr-Ldim, ldb);
+                                fgemm (F, FflasNoTrans, FflasNoTrans, M,  Ldim, N-ib-1, Mone, Bcurr, ldb, Rcurr, ldl,  one, Bcurr-Ldim, ldb);
                                 //cerr<<"j avant="<<j<<endl;
                                 //cerr<<"k, ib, j, R "<<k<<" "<<ib<<" "<<j<<" "<<R<<endl;
                                 //cerr<<"M,k="<<M<<" "<<k<<endl;
                                 //cerr<<" ftrsm with M, N="<<Ldim<<" "<<N<<endl;
-                                ftrsm( F, Side, FflasLower, FflasNoTrans, FflasUnit, M, Ldim, one, Lcurr, ldl , Bcurr-Ldim, ldb );
+                                ftrsm (F, Side, FflasLower, FflasNoTrans, FflasUnit, M, Ldim, one, Lcurr, ldl , Bcurr-Ldim, ldb );
                                 //cerr<<"M,k="<<M<<" "<<k<<endl;
                                 //cerr<<" fgemm with M, N, K="<<M-k<<" "<<N<<" "<<Ldim<<endl;
@@ -581 +588 @@
         template <class Field, class Polynomial>
         static std::list<Polynomial>&
-        Frobenius (const Field& F, std::list<Polynomial>& frobeniusForm, 
-                   const size_t N, typename Field::Element * A, const size_t lda, const size_t c);
+        CharpolyArithProg (const Field& F, std::list<Polynomial>& frobeniusForm, 
+                           const size_t N, typename Field::Element * A, const size_t lda, const size_t c);
         template <class Field>
         static void CompressRows (Field& F, const size_t M,

include/fflas-ffpack/ffpack_charpoly.inl

@@ -11 +11 @@
 template <class Field, class Polynomial>
 std::list<Polynomial>&
-FFPACK::CharPoly( const Field& F, std::list<Polynomial>& charp, const size_t N,
-                            typename Field::Element * A, const size_t lda,
-                            const enum FFPACK_CHARPOLY_TAG CharpTag ){
-        switch ( CharpTag ) {
+FFPACK::CharPoly (const Field& F, std::list<Polynomial>& charp, const size_t N,
+                  typename Field::Element * A, const size_t lda,
+                  const enum FFPACK_CHARPOLY_TAG CharpTag){
+        switch (CharpTag) {
         case FfpackLUK:{
                 typename Field::Element * X = new typename Field::Element[N*(N+1)];
-                LUKrylov( F, charp, N, A, lda, X, N);
+                LUKrylov (F, charp, N, A, lda, X, N);
                 delete[] X;
                 return charp;
         }
         case FfpackKG:{
-                return KellerGehrig( F, charp, N, A, lda );
+                return KellerGehrig (F, charp, N, A, lda);
                 break;
         }
@@ -31 +31 @@
         case FfpackKGFast:{
                 size_t mc, mb, j;
-                if (KGFast( F, charp, N, A, lda, &mc, &mb, &j )){
+                if (KGFast (F, charp, N, A, lda, &mc, &mb, &j)){
                         std::cerr<<"NON GENERIC MATRIX PROVIDED TO KELLER-GEHRIG-FAST"<<std::endl;
                 }
@@ -43 +43 @@
         case FfpackHybrid:{
                 typename Field::Element * X = new typename Field::Element[N*(N+1)];
-                LUKrylov_KGFast( F, charp, N, A, lda, X, N);
+                LUKrylov_KGFast (F, charp, N, A, lda, X, N);
                 delete[] X;
+                return charp;
+                break;
+        }
+
+        case FfpackArithProg:{
+                size_t attempts=0;
+                bool cont = false;
+                do{
+                        try {
+                                CharpolyArithProg (F, charp, N, A, lda, 30);
+                        }
+                        catch (CharpolyFailed){
+                                if (attempts<2)
+                                        cont = true;
+                                else
+                                        CharPoly(F, charp, N, A, lda, FfpackLUK);
+
+                        }
+                } while (cont);
                 return charp;
                 break;
@@ -50 +69 @@
         default:{
                 typename Field::Element * X = new typename Field::Element[N*(N+1)];
-                LUKrylov( F, charp, N, A, lda, X, N);
+                LUKrylov (F, charp, N, A, lda, X, N);
                 delete[] X;
                 return charp;
@@ -60 +79 @@
 template <class Field, class Polynomial>
 std::list<Polynomial>&
-FFPACK::LUKrylov( const Field& F, std::list<Polynomial>& charp, const size_t N,
-                  typename Field::Element * A, const size_t lda,
-                  typename Field::Element * X, const size_t ldx){
+FFPACK::LUKrylov (const Field& F, std::list<Polynomial>& charp, const size_t N,
+                 typename Field::Element * A, const size_t lda,
+                 typename Field::Element * X, const size_t ldx){
         
         typedef typename Field::Element elt;
@@ -73 +92 @@
         charp.clear();
         int nbfac = 0;
-        while( Ncurr > 0 ){
+        while (Ncurr > 0){
                 size_t P[Ncurr];
-                Polynomial minP;//=new  Polynomial();
-                MinPoly( F, minP, Ncurr, A, lda, X2, ldx, P );
+                Polynomial minP;//=new Polynomial();
+                MinPoly (F, minP, Ncurr, A, lda, X2, ldx, P);
                 int k = minP.size()-1; // degre of minpoly
-                if ( (k==1) && F.isZero( (minP)[0] ) ){ // minpoly is X
+                if ((k==1) && F.isZero ((minP)[0])){ // minpoly is X
                         Ai = A;
                         int j = Ncurr*Ncurr;
-                        while ( j-- && F.isZero(*(Ai++)) );
-                        if ( !j ){ // A is 0, CharPoly=X^n
+                        while (j-- && F.isZero(*(Ai++)));
+                        if (!j){ // A is 0, CharPoly=X^n
                                 minP.resize(Ncurr+1);
                                 (minP)[1] = zero;
@@ -90 +109 @@
                 }
                 nbfac++;
-                charp.push_front( minP );
-                if ( k==Ncurr ){
+                charp.push_front (minP);
+                if (k==Ncurr){
                         return charp;
                 }
@@ -97 +116 @@
                 elt * X21 = X2 + k*ldx;
                 elt * X22 = X21 + k;
-                //  Compute the n-k last rows of A' = PA^tP^t in X2_
+                // Compute the n-k last rows of A' = PA^tP^t in X2_
                 // A = A . P^t
-                applyP( F, FflasRight, FflasTrans, Ncurr, 0, k, A, lda, P);
+                applyP (F, FflasRight, FflasTrans, Ncurr, 0, k, A, lda, P);
                 // Copy X2_ = (A'_2)^t
-                for ( Xi = X21, Ai = A+k; Xi != X21 + Nrest*ldx; Ai++, Xi+=ldx-Ncurr )
-                        for ( size_t jj=0; jj<Ncurr*lda; jj+=lda )
+                for (Xi = X21, Ai = A+k; Xi != X21 + Nrest*ldx; Ai++, Xi+=ldx-Ncurr)
+                        for (size_t jj=0; jj<Ncurr*lda; jj+=lda)
                                 *(Xi++) = *(Ai+jj);
                 // A = A . P : Undo the permutation on A
-                applyP( F, FflasRight, FflasNoTrans, Ncurr, 0, k, A, lda, P);
-                // X2_ = X2_ . P^t (=  ( P A^t P^t )2_ ) 
-                applyP( F, FflasRight, FflasTrans, Nrest, 0, k, X21, ldx, P);
+                applyP (F, FflasRight, FflasNoTrans, Ncurr, 0, k, A, lda, P);
+                // X2_ = X2_ . P^t (=  (P A^t P^t)2_) 
+                applyP (F, FflasRight, FflasTrans, Nrest, 0, k, X21, ldx, P);
                 // X21 = X21 . S1^-1
                 ftrsm(F, FflasRight, FflasUpper, FflasNoTrans, FflasUnit, Nrest, k,
-                      one, X2, ldx, X21, ldx);  
+                      one, X2, ldx, X21, ldx); 
                 // Creation of the matrix A2 for recurise call 
-                for ( Xi = X22,  Ai = A;
-                      Xi != X22 + Nrest*ldx;
-                      Xi += (ldx-Nrest), Ai += (lda-Nrest) )
-                        for ( size_t jj=0; jj<Nrest; ++jj )
+                for (Xi = X22, Ai = A;
+                     Xi != X22 + Nrest*ldx;
+                     Xi += (ldx-Nrest), Ai += (lda-Nrest))
+                        for (size_t jj=0; jj<Nrest; ++jj)
                                 *(Ai++) = *(Xi++);
-                fgemm( F, FflasNoTrans, FflasNoTrans, Nrest, Nrest, k, Mone,
-                       X21, ldx, X2+k, ldx, one, A, lda );
+                fgemm (F, FflasNoTrans, FflasNoTrans, Nrest, Nrest, k, Mone,
+                       X21, ldx, X2+k, ldx, one, A, lda);
                 X2 = X22;
                 Ncurr = Nrest;
@@ -127 +146 @@
 template <class Field, class Polynomial>
 std::list<Polynomial>&
-FFPACK::LUKrylov_KGFast( const Field& F, std::list<Polynomial>& charp, const size_t N,
-                                   typename Field::Element * A, const size_t lda,
-                                   typename Field::Element * X, const size_t ldx){
+FFPACK::LUKrylov_KGFast (const Field& F, std::list<Polynomial>& charp, const size_t N,
+                         typename Field::Element * A, const size_t lda,
+                         typename Field::Element * X, const size_t ldx){
         
         typedef typename Field::Element elt;
@@ -139 +158 @@
         size_t kg_mc, kg_mb, kg_j;
         
-        if (!KGFast( F, charp, N, A, lda, &kg_mc, &kg_mb, &kg_j ))
+        if (!KGFast (F, charp, N, A, lda, &kg_mc, &kg_mb, &kg_j))
                 return charp;
         else{// Matrix A is not generic
@@ -147 +166 @@
                 size_t *P = new size_t[N];
 
-                MinPoly( F, *minP, N, A, lda, X, ldx, P, FfpackKGF, kg_mc, kg_mb, kg_j );
+                MinPoly (F, *minP, N, A, lda, X, ldx, P, FfpackKGF, kg_mc, kg_mb, kg_j);
                 size_t k = minP->size()-1; // degre of minpoly
-                if ( (k==1) && F.isZero( (*minP)[0] ) ){ // minpoly is X
+                if ((k==1) && F.isZero ((*minP)[0])){ // minpoly is X
                         Ai = A;
                         int j = N*N;
-                        while ( j-- && F.isZero(*(Ai++)) );
-                        if ( !j ){ // A is 0, CharPoly=X^n
+                        while (j-- && F.isZero(*(Ai++)));
+                        if (!j){ // A is 0, CharPoly=X^n
                                 minP->resize(N+1);
                                 (*minP)[1] = zero;
@@ -161 +180 @@
                 }
                 
-                if ( k==N ){
+                if (k==N){
                         charp.clear();
                         charp.push_front(*minP); // CharPoly = MinPoly
@@ -176 +195 @@
                 size_t imax = kg_mc+kg_mb;
                 // First Id
-                for ( size_t j = 0; j < lambda; ++j){
+                for (size_t j = 0; j < lambda; ++j){
                         for (size_t i=0; i<imax; ++i)
-                                F.assign( *(A+j+i*lda), zero);
-                        F.assign( *(A+j+imax*lda), one);
+                                F.assign (*(A+j+i*lda), zero);
+                        F.assign (*(A+j+imax*lda), one);
                         for (size_t i=imax+1; i<N; ++i)
-                                F.assign( *(A+j+i*lda), zero);
+                                F.assign (*(A+j+i*lda), zero);
                         ++imax;
                 }
                 // Column block B
                 for (typename Field::Element* Ai=A; Ai<A+N*lda; Ai+=lda)
-                        fcopy( F, kg_mb, Ai+lambda, 1, Ai+N-kg_mc-kg_mb, 1 );
+                        fcopy (F, kg_mb, Ai+lambda, 1, Ai+N-kg_mc-kg_mb, 1);
 
                 // Second Id block
@@ -192 +211 @@
                 for (size_t j = 0; j< kg_j*kg_mc; ++j){
                         for (size_t i = 0; i<imax; ++i)
-                                F.assign( *(A+lambda+kg_mb+j+i*lda), zero );
-                        F.assign( *(A+lambda+kg_mb+j+imax*lda), one );
+                                F.assign (*(A+lambda+kg_mb+j+i*lda), zero);
+                        F.assign (*(A+lambda+kg_mb+j+imax*lda), one);
                         for (size_t i = imax+1; i<N; ++i)
-                                F.assign( *(A+lambda+kg_mb+j+i*lda), zero );
+                                F.assign (*(A+lambda+kg_mb+j+i*lda), zero);
                         ++imax;
                 }
                 
-                //  Compute the n-k last rows of A' = PA^tP^t in X2_
+                // Compute the n-k last rows of A' = PA^tP^t in X2_
                 
                 // A = P . A 
-                applyP( F, FflasLeft, FflasNoTrans, N, 0, k, 
+                applyP (F, FflasLeft, FflasNoTrans, N, 0, k, 
                         const_cast<typename Field::Element* &>(A), lda, P);
                 
                 // Copy X2_ = (A'2_)
-                for ( Xi = X21, Ai = A+k*lda; Xi != X21 + Nrest*ldx; Ai+=lda-N, Xi+=ldx-N ){
-                        for ( size_t jj=0; jj<N; ++jj ){
+                for (Xi = X21, Ai = A+k*lda; Xi != X21 + Nrest*ldx; Ai+=lda-N, Xi+=ldx-N){
+                        for (size_t jj=0; jj<N; ++jj){
                                 *(Xi++) = *(Ai++);
                         }
@@ -213 +232 @@
 
                 // A = P^t . A : Undo the permutation on A
-                applyP( F, FflasLeft, FflasTrans, N, 0, k, 
+                applyP (F, FflasLeft, FflasTrans, N, 0, k, 
                         const_cast<typename Field::Element* &>(A), lda, P);
         
-                // X2_ = X2_ . P^t (=  ( P A P^t )2_ ) 
-                applyP( F, FflasRight, FflasTrans, Nrest, 0, k, X21, ldx, P);
+                // X2_ = X2_ . P^t (=  (P A P^t)2_) 
+                applyP (F, FflasRight, FflasTrans, Nrest, 0, k, X21, ldx, P);
 
                 // X21 = X21 . S1^-1
@@ -226 +245 @@
                 elt * A2 = new elt[Nrest*Nrest];
         
-                for ( Xi = X22,  A2i = A2;
-                      Xi != X22 + Nrest*ldx;
-                      Xi += (ldx-Nrest) ){
-                        for ( size_t jj=0; jj<Nrest; ++jj ){
+                for (Xi = X22, A2i = A2;
+                     Xi != X22 + Nrest*ldx;
+                     Xi += (ldx-Nrest)){
+                        for (size_t jj=0; jj<Nrest; ++jj){
                                 *(A2i++) = *(Xi++);
                         }
                 }
-                fgemm( F, FflasNoTrans, FflasNoTrans, Nrest, Nrest, k, Mone,
-                       X21, ldx, X+k, ldx, one, A2, Nrest );
+                fgemm (F, FflasNoTrans, FflasNoTrans, Nrest, Nrest, k, Mone,
+                       X21, ldx, X+k, ldx, one, A2, Nrest);
         
                 // Recursive call on X22
-                LUKrylov_KGFast( F, charp, Nrest, A2, Nrest, X22, ldx );
-                charp.push_front( *minP );
+                LUKrylov_KGFast (F, charp, Nrest, A2, Nrest, X22, ldx);
+                charp.push_front (*minP);
                 delete[] P;
                 delete[] A2;

include/fflas-ffpack/ffpack_frobenius.inl

@@ -14 +14 @@
 
 //---------------------------------------------------------------------
-// Frobenius: Las Vegas algorithm to compute the Frobenius normal form over a large field
+// CharpolyArithProg: Las Vegas algorithm to compute the Charpoly normal form over a large field
 //---------------------------------------------------------------------
 
@@ -146 +146 @@
 template <class Field, class Polynomial>
 std::list<Polynomial>&
-FFPACK::Frobenius (const Field& F, std::list<Polynomial>& frobeniusForm, 
-                   const size_t N, typename Field::Element * A, const size_t lda, const size_t c){
+FFPACK::CharpolyArithProg (const Field& F, std::list<Polynomial>& frobeniusForm, 
+                           const size_t N, typename Field::Element * A, const size_t lda, const size_t c){
         
         static typename Field::Element one, zero, mone;
@@ -156 +156 @@
         size_t * rp = new size_t[2*N];
 
-        size_t noc = static_cast<size_t>(ceil(N/double(c)));
-
-        // Building the workplace matrix Aw
+        size_t noc = static_cast<size_t>(ceil(double(N)/double(c)));
+
+        // Building the workplace matrix 
         typename Field::Element *K = new typename Field::Element[N*(noc*c)];
         typename Field::Element *K2 = new typename Field::Element[N*(noc*c)];
@@ -170 +170 @@
         
         size_t *dA = new size_t[N]; //PA
-        size_t *dK = new size_t[N];
+        //      size_t *dK = new size_t[N];
 
 
@@ -185 +185 @@
 //             K+(c-2)*noc*ldk, ldk, A, lda, zero, K+(c-1)*noc*ldk, ldk);
 
-//      write_field(F,cerr<<"K = "<<endl,K, c*noc,N,ldk);
-//      write_field(F,cerr<<"A = "<<endl,A, N,N,ldk);
+        //write_field(F,cerr<<"K = "<<endl,K, c*noc,N,ldk);
+        //write_field(F,cerr<<"A = "<<endl,A, N,N,ldk);
         size_t * Pk = new size_t[N];
         size_t * Qk = new size_t[c*noc];
@@ -193 +193 @@
         size_t w_idx = 0;
         for (size_t i=0; i<noc; ++i)
-                for (size_t j=0; j<c; ++j)
-                        fcopy(F, N, (K2+(w_idx++)*ldk), 1, (K+(i+j*noc)*ldk), 1);
+                for (size_t j=0; j<c; ++j, w_idx++)
+                        fcopy(F, N, (K2+(w_idx)*ldk), 1, (K+(i+j*noc)*ldk), 1);
 
         for (size_t i=0; i<noc*c; ++i)
@@ -200 +200 @@
         
         size_t R = LUdivine(F, FflasNonUnit, noc*c, N, K, ldk, Pk, Qk, FfpackLQUP);
-        if (R<N){
-                std::cerr<<"Preconditionning failed; not yet implemented"<<std::endl;
-                exit(-1);
-        }
-        
+
+        //write_field(F,cerr<<"Apres LUDivine K = "<<endl,K, c*noc,N,ldk);
+        // cerr<<"Q = ";
+//      for (size_t i=0; i<noc*c; ++i)
+//              cerr<<Qk[i]<<" ";
+//      cerr<<endl;
         size_t row_idx = 0;
-        size_t * degini = new size_t[noc];
+        size_t * dK = new size_t[noc];
         for (size_t i=0; i<noc; ++i)
-                degini[i]=0;
-        // size_t * degK = new size_t[c];
-//      for (size_t i=0; i<c; ++i)
-//              degK[i]=0;
+                dK[i]=0;
         
         size_t i=0;
         row_idx=0;
+        size_t dold = c;
+        size_t nb_full_blocks = 0;
+        size_t Mk = 0;
         for (size_t k = 0; k<noc; ++k){
                 size_t d = 0;
                 while ( (d<c) && (row_idx<R) && (Qk[row_idx] == i)) {i++; row_idx++; d++;}
-                degini[k] = d;
+                if (d > dold){
+                        std::cerr << "FAIL in preconditionning phase:"
+                             << " degree sequence is not monotonically not increasing"
+                             << std::endl;
+                        //exit (-1);
+                        throw CharpolyFailed();
+                }
+                        
+                dK[k] = dold = d;
+                Mk++;
+                if (d == c)
+                        nb_full_blocks++;
                 i = Qk[row_idx];
         }
-        
-//      for (size_t i = 0; i<N; ++i)
-//              if (Qk[row_idx] == i){
-//                      degini[i % noc] ++;
-//      //              degK[i/noc] ++;
-//                      row_idx ++;
+        //cerr<<"Mk = "<<Mk<<endl;
+        // Selection of the last iterate of each block
+        size_t bk_idx = 0;
+        //write_field (F,cerr<<"Avant fgemm K2 = "<<endl, K2, noc*c, N, ldk);
+
+        for (size_t i = 0; i < Mk; ++i){
+#if DEBUG
+                std::cerr<<"dK ["<<i<<"] = "<<dK[i]<<std::endl;
+#endif
+                fcopy (F, N, (K2+i*ldk), 1, (K2 + (bk_idx + dK[i]-1)*ldk), 1);
+                //cerr<<"(bk_idx+dK[i]-1) = "<<(bk_idx+dK[i]-1)<<endl;
+                bk_idx+= c;
+        }
+
+        typename Field::Element* K2b = K2+Mk*ldk;
+        // K <- K A^T 
+