Changeset 4 for include/fflas-ffpack/ffpack_frobenius.inl

Timestamp:

03/13/07 10:33:27 (2 years ago)

Author:

pernet

Message:

Clean up the code for charpoly, + misc updates

Files:

• include/fflas-ffpack/ffpack_frobenius.inl (modified) (22 diffs)

include/fflas-ffpack/ffpack_frobenius.inl

@@ -14 +14 @@
 
 //---------------------------------------------------------------------
-// CharpolyArithProg: Las Vegas algorithm to compute the Charpoly normal form over a large field
+// CharpolyArithProg: Las Vegas algorithm to compute the Charpoly normal
+// form over a large field (Z/pZ, s.t.  p > 2n^2)
 //---------------------------------------------------------------------
-
-
-// 
-template <class Field>
-void FFPACK::CompressRowsQK (Field& F, const size_t M,
-                           typename Field::Element * A, const size_t lda,
-                           typename Field::Element * tmp, const size_t ldtmp,
-                           const size_t * d, const size_t deg,const size_t nb_blocs){
-
-        size_t currtmp = 0;
-        size_t currw = d[0]-1;
-        size_t currr = d[0]-1;
-        for (int i = 0; i< int(nb_blocs)-1; ++i){
-                for (int j = d[i]-1; j<deg-1; ++j, currr++, currtmp++)
-                        fcopy(F, M, tmp + currtmp*ldtmp, 1,  A + currr*lda, 1);
-                //cerr<<"d["<<i<<"] = "<<d[i]<<endl;
-                for (int j=0; j < d[i+1] -1; ++j, currr++, currw++){
-                        fcopy(F, M, A + (currw)*lda, 1, A+(currr)*lda, 1);
-                }
-                //cerr<<"Ok"<<endl;
-        }
-
-        //write_field(F,cerr<<"avant tmp"<<endl,A, currtmp, M, lda);
-        for (int i=0; i < currtmp; ++i, currw++){
-                fcopy (F, M, A + (currw)*lda, 1, tmp + i*ldtmp, 1);
-        }
-}
-
-template <class Field>
-void FFPACK::CompressRows (Field& F, const size_t M,
-                             typename Field::Element * A, const size_t lda,
-                             typename Field::Element * tmp, const size_t ldtmp,
-                             const size_t * d, const size_t nb_blocs){
-
-        size_t currd = d[0]-1;
-        size_t curri = d[0]-1;
-        for (int i = 0; i< int(nb_blocs)-1; ++i){
-                fcopy(F, M, tmp + i*ldtmp, 1,  A + currd*lda, 1);
-                for (int j=0; j < d[i+1] -1; ++j){
-                        fcopy(F, M, A + (curri++)*lda, 1, A+(currd+j+1)*lda, 1);
-                }
-                currd += d[i+1];
-        }
-        for (int i=0; i < int(nb_blocs)-1; ++i){
-                fcopy (F, M, A + (curri++)*lda, 1, tmp + i*ldtmp, 1);
-        }
-}
-template <class Field>
-void FFPACK::DeCompressRowsQK (Field& F, const size_t M, const size_t N,
-                               typename Field::Element * A, const size_t lda,
-                               typename Field::Element * tmp, const size_t ldtmp,
-                               const size_t * d, const size_t deg,const size_t nb_blocs){
-        
-        size_t zeroblockdim = 1; // the last block contributes with 1
-        size_t currtmp = 0;
-        for (int i=0; i<int(nb_blocs)-1; ++i)
-                zeroblockdim += deg - d[i];
-        //cerr<<"zeroblockdim = "<<zeroblockdim<<endl;
-        for (int i=0; i < zeroblockdim - 1; ++i, ++currtmp)
-                fcopy(F, M, tmp + currtmp*ldtmp, 1,  A + (N - zeroblockdim +i)*lda, 1);
-        //write_field(F,cerr<<"tmp= "<<endl,tmp,zeroblockdim - 1, M,ldtmp);
-        currtmp--;
-        size_t w_idx = N - 2;
-        size_t r_idx = N - zeroblockdim - 1;
-
-        for (int i = int(nb_blocs)-2; i>=0; --i){
-                for (size_t j = 0; j < d [i+1] - 1; ++j)
-                        fcopy (F, M, A + (w_idx--)*lda, 1, A + (r_idx--)*lda, 1);
-                for (size_t j = 0; j < deg - d[i]; ++j)
-                        fcopy (F, M, A + (w_idx--)*lda, 1, tmp + (currtmp--)*ldtmp, 1);
-        }
-}
-
-template <class Field>
-void FFPACK::DeCompressRows (Field& F, const size_t M, const size_t N,
-                             typename Field::Element * A, const size_t lda,
-                             typename Field::Element * tmp, const size_t ldtmp,
-                             const size_t * d, const size_t nb_blocs){
-        
-        for (int i=0; i<int(nb_blocs)-1; ++i)
-                fcopy(F, M, tmp + i*ldtmp, 1, A + (N-nb_blocs+i)*lda, 1);
-        
-        size_t w_idx = N - 2;
-        size_t r_idx = N - nb_blocs - 1;
-        for (int i = int(nb_blocs)-2; i>=0; --i){
-                for (size_t j = 0; j<d[i+1]-1; ++j)
-                        fcopy (F, M, A + (w_idx--)*lda, 1, A + (r_idx--)*lda, 1);
-                fcopy (F, M, A + (w_idx--)*lda, 1, tmp + i*ldtmp, 1);
-        }
-}
-
-template <class Field>
-void FFPACK::CompressRowsQA (Field& F, const size_t M,
-                             typename Field::Element * A, const size_t lda,
-                             typename Field::Element * tmp, const size_t ldtmp,
-                             const size_t * d, const size_t nb_blocs){
-
-        size_t currd = 0;
-        size_t curri = 0;
-        for (size_t i = 0; i< nb_blocs; ++i){
-                fcopy(F, M, tmp + i*ldtmp, 1,  A + currd*lda, 1);
-                for (size_t j=0; j < d[i] -1; ++j)
-                        fcopy(F, M, A + (curri++)*lda, 1, A+(currd+j+1)*lda, 1);
-                currd += d[i];
-        }
-        for (size_t i=0; i < nb_blocs; ++i)
-                fcopy (F, M, A + (curri++)*lda, 1, tmp + i*ldtmp, 1);
-}
-
-template <class Field>
-void FFPACK::DeCompressRowsQA (Field& F, const size_t M, const size_t N,
-                               typename Field::Element * A, const size_t lda,
-                               typename Field::Element * tmp, const size_t ldtmp,
-                               const size_t * d, const size_t nb_blocs){
-        
-        for (size_t i=0; i<nb_blocs; ++i)
-                fcopy(F, M, tmp + i*ldtmp, 1, A + (N-nb_blocs+i)*lda, 1);
-
-        size_t w_idx = N - 1;
-        size_t r_idx = N - nb_blocs - 1;
-        for (int i = int(nb_blocs)-1; i>=0; --i){
-                for (size_t j = 0; j<d[i]-1; ++j)
-                        fcopy (F, M, A + (w_idx--)*lda, 1, A + (r_idx--)*lda, 1);
-                fcopy (F, M, A + (w_idx--)*lda, 1, tmp + i*ldtmp, 1);
-        }
-}
-
-
 template <class Field, class Polynomial>
 std::list<Polynomial>&
 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){
+                           const size_t N, typename Field::Element * A, const size_t lda,
+                           const size_t c){
         
         static typename Field::Element one, zero, mone;
@@ -155 +29 @@
 
         size_t * rp = new size_t[2*N];
-
         size_t noc = static_cast<size_t>(ceil(double(N)/double(c)));
 
@@ -162 +35 @@
         typename Field::Element *K2 = new typename Field::Element[N*(noc*c)];
         size_t ldk = N;
-        // Permutations such that
-        // A = QA [ I Ac ] PA
-        //        [ 0 Ac ]
-        //
-        // K = [ I Kc ] PK
-        //     [ 0 Kc ]
-        
+
         size_t *dA = new size_t[N]; //PA
-        //      size_t *dK = new size_t[N];
-
-
-        size_t Nrest = N-(c-1)*noc;
+        size_t *dK = new size_t[noc*c];
+        for (size_t i=0; i<noc; ++i)
+                dK[i]=0;
+
+        // Picking a random noc x N block vector U^T
         typename Field::RandIter g (F);
         for (size_t i = 0; i < noc; ++i)
@@ -179 +47 @@
                         g.random( *(K + i*ldk +j) );
 
+        // Computing the bloc Krylov matrix [U AU .. A^(c-1) U]^T
         for (size_t i = 1; i<c; ++i)
                 fgemm( F, FflasNoTrans, FflasTrans,  noc, N, N, one,
                        K+(i-1)*noc*ldk, ldk, A, lda, zero, K+i*noc*ldk, ldk);
-//      fgemm( F, FflasNoTrans, FflasTrans, Nrest, N, N, one,
-//             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);
-        size_t * Pk = new size_t[N];
-        size_t * Qk = new size_t[c*noc];
 
         // K2 <- K (re-ordering)
@@ -196 +58 @@
                         fcopy(F, N, (K2+(w_idx)*ldk), 1, (K+(i+j*noc)*ldk), 1);
 
+        // Copying K <- K2
         for (size_t i=0; i<noc*c; ++i)
                 fcopy (F, N, (K+i*ldk), 1, K2+i*ldk, 1);
         
+        size_t * Pk = new size_t[N];
+        size_t * Qk = new size_t[c*noc];
+        for (size_t i=0; i<c*noc; ++i)
+                Qk[i] = 0;
+        for (size_t i=0; i<N; ++i)
+                Pk[i] = 0;
         size_t R = LUdivine(F, FflasNonUnit, noc*c, N, K, ldk, Pk, Qk, FfpackLQUP);
 
-        //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 * dK = new size_t[noc];
-        for (size_t i=0; i<noc; ++i)
-                dK[i]=0;
-        
         size_t i=0;
-        row_idx=0;
         size_t dold = c;
         size_t nb_full_blocks = 0;
         size_t Mk = 0;
+        // Determining the degree sequence dK
         for (size_t k = 0; k<noc; ++k){
                 size_t d = 0;
@@ -226 +86 @@
                         throw CharpolyFailed();
                 }
-                        
                 dK[k] = dold = d;
                 Mk++;
                 if (d == c)
                         nb_full_blocks++;
-                i = Qk[row_idx];
-        }
-        //cerr<<"Mk = "<<Mk<<endl;
+                if (row_idx < noc*c)
+                        i = Qk[row_idx];
+        }
+
         // 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;
-        }
-
+                bk_idx += c;
+        }
         typename Field::Element* K2b = K2+Mk*ldk;
+
         // K <- K A^T 
-
-        //write_field (F,cerr<<"Avant fgemm K2 = "<<endl, K2, noc*c, N, ldk);
         fgemm( F, FflasNoTrans, FflasTrans, Mk, N, N, one,  K2, ldk, A, lda, zero, K2b, ldk);
-
-        //write_field (F,cerr<<"Apres fgemm K2b = "<<endl, K2b, Mk, N, ldk);
 
         // K <- K P^T
         applyP (F, FflasRight, FflasTrans, Mk, 0, R, K2b, ldk, Pk);
 
-        //write_field (F,cerr<<"Apres applyP K2b = "<<endl, K2b, Mk, N, ldk);
         // K <- K U^-1
         ftrsm (F, FflasRight, FflasUpper, FflasNoTrans, FflasNonUnit, Mk, R, one, K, ldk, K2b, ldk);
 
-        //write_field (F,cerr<<"Apres ftrsm K2b = "<<endl, K2b, Mk, N, ldk);
-        // K <- K Q^T
-        //applyP (F, FflasRight, FflasNoTrans, Mk, 0, R, K2b, ldk, Qk);
-
-
-//      row_idx = 0;
-//      for (size_t i=0; i<R; ++i){
-//              for (size_t k=row_idx; k<Qk[i]; ++k) // iterate on the linearly dependent rows
-//                      for (size_t j=0; j<i; ++j)
-//                              F.assign( *(K+j+k*ldk),zero);
-//              row_idx = Qk[i]+1;
-//      }
-        //cerr<<"c, noc, row_idx, N = "<<c<<" "<<noc<<" "<<row_idx<<" "<<N<<endl;
-        //row_idx = noc*c -1;
-//      while (row_idx < noc*c){
-//              for (size_t j=0; j<MIN(row_idx,N); ++j)
-//                      F.assign( *(K + j + (row_idx)*ldk),zero);
-//              row_idx++;
-//      }
-        //write_field (F,cerr<<"Avant solvelb K2b = "<<endl, K2b, Mk, N, ldk);
+        // L <-  Q^T L
+        applyP(F, FflasLeft, FflasNoTrans, N, 0, R, K, ldk, Qk);
+
         // K <- K L^-1
-        //compressing L
-        //      write_field (F,cerr<<"Avant compression L = "<<endl, K, noc*c, N, ldk);
-        applyP(F, FflasLeft, FflasNoTrans, N, 0, R, K, ldk, Qk);
-        //write_field (F,cerr<<"Apres compression L = "<<endl, K, noc*c, N, ldk);
-
         ftrsm (F, FflasRight, FflasLower, FflasNoTrans, FflasUnit, Mk, R, one, K, ldk, K2b, ldk);
-        //undoing permutation on K
+
+        //undoing permutation on L
         applyP(F, FflasLeft, FflasTrans, N, 0, R, K, ldk, Qk);
-        //      solveLB2 (F, FflasRight, Mk, noc*c, R, K, ldk, Qk, K2b, ldk);
-        //§!!!!!
-        //applyP (F, FflasRight, FflasTrans, Mk, 0, R, K2b, ldk, Qk);
-        //write_field (F,cerr<<"Apres solvelb K2b = "<<endl, K2b, Mk, N, ldk);
-        
+
         // Recovery of the completed invariant factors
-        //cerr<<"// Recovery of the completed invariant factors"<<endl;
         size_t Ma = Mk;
         size_t Ncurr = R;
         size_t offset = Ncurr-1;
-        size_t oldNcurr = Ncurr;
-        Polynomial *P;
-        //write_field (F,cerr<<"Recovery of the completed blocks"<<endl, K2b, Mk, N, ldk);
         for (size_t i=Mk-1; i>=nb_full_blocks+1;  --i)
                 if (dK[i] >= 1){ 
-                        Polynomial P (dK [i]);
+                        Polynomial P (dK [i]+1);
+                        F.assign(P[dK[i]], one);
                         for (size_t j=0; j < dK [i]; ++j){
                                 F.neg (P [dK [i]-j-1], *(K2b + i*ldk + (offset-j)));
-                                //cerr<<" "<<(*(K2b + i*ldk + (offset-j)));
                         }
-                        //cerr<<endl;
-                        // for (size_t j=0; j<dK[i]; ++j)
-//                              cerr<<" "<<(P)[j];
                         frobeniusForm.push_front(P);
-                        //cerr<<"Storing polynomial "<<i<<endl;
                         offset -= dK [i];
                         Ncurr -= dK [i];
-                        //block_idx--;
                         Ma--;
                 }
         Mk = Ma;
-        //cerr<<"taille frob = "<<frobeniusForm.size()<<endl;
-        //write_field(F,cerr,K,oldNcurr, Mk, ldk);
-        // cerr<<"nb_full_blocks, Mk, Ncurr, oldNcurr, offset  = "
-//                  <<nb_full_blocks<<" "<<Mk<<" "<<Ncurr<<" "<<oldNcurr<<" "<<offset<<endl;
 
 //      for (size_t i= offset+1; i<oldNcurr; ++i)
@@ -338 +149 @@
         if (R<N){
                 
-                std::cerr<<"Preconditionning failed; missing rank = "<<N-R<<" completing the Krylov matrix"
-                         <<std::endl;
-                // exit(-1);
-                
-                Nrest = N-R;
-                //cerr<<"Nrest = "<<Nrest<<endl;
+//              std::cerr<<"Preconditionning failed; missing rank = "<<N-R
+//                       <<" completing the Krylov matrix"
+//                       <<std::endl;
+                size_t Nrest = N-R;
                 typename Field::Element * K21 = K + R*ldk;
                 typename Field::Element * K22 = K21 + R;
                 typename Field::Element * Ki, *Ai;
+
                 //  Compute the n-k last rows of A' = P A^T P^T in K2_
                 // A = A . P^t
                 applyP( F, FflasRight, FflasTrans, N, 0, R, A, lda, Pk);
+
                 // Copy K2_ = (A'_2)^t
                 for (Ki = K21, Ai = A+R; Ki != K21 + Nrest*ldk; Ai++, Ki+=ldk-N)
@@ -355 +166 @@
                                 *(Ki++) = *(Ai+j);
 
-                //write_field (F, cerr<<"A = "<<endl, A, N, N, lda);
-                
-                //write_field (F, cerr<<"K2_ = "<<endl, K21, Nrest, N, ldk);
                 // A = A . P : Undo the permutation on A
                 applyP( F, FflasRight, FflasNoTrans, N, 0, R, A, lda, Pk);
+
                 // K2_ = K2_ . P^t (=  ( P A^t P^t )2_ ) 
                 applyP( F, FflasRight, FflasTrans, Nrest, 0, R, K21, ldk, Pk);
-                //write_field (F, cerr<<"K2_ ( = P A^T P^T) = "<<endl, K21, Nrest, N, ldk);
+                
                 // K21 = K21 . S1^-1
-                //write_field (F, cerr<<"S1  = "<<endl, K,R,R, ldk);
-                ftrsm(F, FflasRight, FflasUpper, FflasNoTrans, FflasNonUnit, Nrest, R,
-                      one, K, ldk, K21, ldk);  
-                //write_field (F, cerr<<"% S1^-1 = "<<endl, K21, Nrest, N, ldk);
+                ftrsm (F, FflasRight, FflasUpper, FflasNoTrans, FflasNonUnit, Nrest, R,
+                       one, K, ldk, K21, ldk);  
+
                 typename Field::Element * Arec = new typename Field::Element[Nrest*Nrest];
                 size_t ldarec = Nrest;
@@ -382 +190 @@
                 std::list<Polynomial> polyList;
                 polyList.clear();
-                //write_field (F, cerr<<"Recursive call with A = "<<endl, Arec, Nrest, Nrest, Nrest);
-
+
+                // Recursive call on the complementary subspace
                 CharpolyArithProg (F, polyList, Nrest, Arec, ldarec, 2);
-                //cerr<<"taille polyList = "<<polyList.size()<<endl;
-                //cerr<<"Apres recursive call"<<endl;
                 delete[] Arec;
                 frobeniusForm.merge(polyList);
-                //cerr<<"taille frob = "<<frobeniusForm.size()<<endl;
-                //              typename list <Polynomial>::iterator pl_it = polyList.begin();
-//              while (pl_it != polyList.end())
-//                      frobeniusForm.push_front(*(pl_it++));
-                
-        }
-        
-        
+        }
+
         delete[] Pk;
         delete[] Qk;
@@ -402 +202 @@
         for (size_t i=0; i<Mk; ++i)
                 dA[i] = dK[i];
-        
         bk_idx = 0;
         
-        //      write_field(F,cerr<<"Apres preconditionnement, K = "<<endl, K2+noc*ldk, noc,N,ldk);
         typename Field::Element *Arp = new typename Field::Element[Ncurr*Ma];
         typename Field::Element *Ac = new typename Field::Element[Ncurr*Ma];
@@ -414 +212 @@
                 for (size_t j=0; j<Ma; ++j)
                         *(K+i*ldk+j) = *(Ac + i*Ma +j) = *(K2b + i + (j)*ldk);
-        //      for (size_t i=0; i < noc; ++i){
-//              for (size_t k=0; k < degini[i]; ++k){
-//                      for (size_t j=0; j < noc; ++j)
-//                              *(Ac + (w_idx)*N + j) = *(K2 + (noc+j)*ldk + i + bk_idx);
-//                      bk_idx += degK[k];
-//                      w_idx++;
-//                      }
-//              bk_idx = 0;
-//      }
-
-        //write_field(F,cerr<<"Apres preconditionnement, Ac = "<<endl, Ac, Ncurr, Ma, Ma);
-        
-        size_t block_idx, it_idx, rp_val;
-
-        //cerr<<"Avant a boucle : Ma, Ncurr = "<<Ma<<" "<<Ncurr<<endl;
+
+        size_t block_idx, it_idx, rp_val;
+
+        // Main loop of the arithmetic progession
         while ((nb_full_blocks >= 1) && (Mk > 1)) {
                 delete[] K;
@@ -436 +223 @@
                 ldk = Ma;
  
-                //cerr<<deg<<endl;
-                //write_field(F, cerr<<"Ac = "<<endl, Ac, Ncurr, Ma, N);
                 // Computation of the rank profile
                 for (size_t i=0; i < Ncurr; ++i)
@@ -448 +233 @@
                         std::cerr<<"FAIL R<Ncurr"<<std::endl;
                         throw CharpolyFailed();
-                        //exit(-1);
-                }
+                }
+
                 // Computation of the degree vector dK
-                //cerr<<"// Computation of the degree vector dK"<<endl;
                 it_idx = 0;
                 rp_val = 0;
@@ -458 +242 @@
                 block_idx = 0;
                 nb_full_blocks = 0;
-#if DEBUG
-                std::cerr<<"dK = ";
-#endif
                 while (dtot<Ncurr){
                         do {g++; rp_val++; it_idx++;}
@@ -471 +252 @@
                         dK[block_idx++] = it_idx;
                         dtot += it_idx;
-#if DEBUG
-                        std::cerr<<dK[block_idx-1]<<" ";
-#endif
                         if (it_idx == deg)
                                 nb_full_blocks ++;
@@ -480 +258 @@
                 }
 
-#if DEBUG
-                std::cerr<<std::endl;
-#endif
                 Mk = block_idx;
-                //cerr<<"Mk, nb_full_blocks = "<<Mk<<" "<<nb_full_blocks<<endl;
-
-
                                 
                 // Selection of dense colums of K 
-                //cerr<<"// Selection of dense colums of K "<<endl;
                 for (size_t i=0; i < nb_full_blocks; ++i){
                         fcopy (F, Ncurr, K+i, ldk, Ac+i, ldac);
                 }
-                //write_field(F, cerr<<"K = "<<endl, K, Ncurr, Mk, ldk);
                 
                 // K <- QK K
-                // cerr<<"// K <- QK K "<<endl;
-//              DeCompressRows (F, Mk - nb_full_blocks, Ncurr-nb_full_blocks*deg, K + nb_full_blocks*(deg-1)*ldk, ldk, Arp, N, dK+nb_full_blocks, Mk-nb_full_blocks);
-//              write_field(F, cerr<<"K = "<<endl, K, Ncurr, Mk, ldk);
-
                 size_t pos = nb_full_blocks*(deg-1);
                 for (size_t i = nb_full_blocks; i < Mk; ++i){
-                        //      cerr<<"*(K + "<<i<<") = "<<(*(K+i))<<"ldk = "<<ldk<<endl;
                         for (size_t j=0; j<Ncurr; ++j)
                                 F.assign (*(K + i + j*ldk), zero);
                         F.assign (*(K + i + (pos + dK[i]-1)*ldk), one);
-                        //cerr<<"dA["<<i<<"] = "<<dA[i]<<endl;
                         pos += dA[i];
                 }
-                //write_field(F, cerr<<"K = "<<endl, K, Ncurr, Mk, ldk);
 
                 // Copying K2 <- K
-                //cerr<<"// Copying K2 <- K"<<endl;
                 for (size_t i=0; i<Mk; ++i)
                         fcopy (F, Ncurr, K2+i, ldk, K+i, ldk);
-                CompressRowsQK (F, Mk, K2 + nb_full_blocks*(deg-1)*ldk, ldk, Arp, ldarp, dK+nb_full_blocks, deg, Mk-nb_full_blocks);
-
-                
+                CompressRowsQK (F, Mk, K2 + nb_full_blocks*(deg-1)*ldk, ldk,
+                                Arp, ldarp, dK+nb_full_blocks, deg, Mk-nb_full_blocks);
+
                 // K <- PA K
-                //cerr<<"// K <- PA K"<<endl;
                 CompressRows (F, nb_full_blocks, K, ldk, Arp, ldarp, dA, Ma);
-                //write_field(F, cerr<<"K = "<<endl, K, Ncurr, Mk, ldk);
                 
                 // A <- newQA^T K (compress)
-                //cerr<<"// K <- newQA^T K (compress)"<<endl;
                 CompressRowsQA (F, Ma, Ac, ldac, Arp, ldarp, dA, Ma);
-                //write_field(F, cerr<<"K = "<<endl, K, Ncurr, Mk, ldk);
-
+                
                 // K <- A K
-                // cerr<<"// K <- A K"<<endl;
-//              cerr<<"C <- AB + C"<<endl;
-                // write_field(F,cerr<<"A = "<<endl, Ac, Ncurr-Ma, Ma, N);
-//              write_field(F,cerr<<"B = "<<endl, K+(Ncurr-Ma)*ldk, Ma,nb_full_blocks, ldk);
-//              write_field(F,cerr<<"C = "<<endl, K, Ncurr-Ma, nb_full_blocks, ldk);
-                
-                fgemm (F, FflasNoTrans, FflasNoTrans, Ncurr-Ma, nb_full_blocks, Ma, one, Ac, ldac, K+(Ncurr-Ma)*ldk, ldk, one, K, ldk);
-                fgemm (F, FflasNoTrans, FflasNoTrans, Ma, nb_full_blocks, Ma, one, Ac+(Ncurr-Ma)*ldac, ldac, K+(Ncurr-Ma)*ldk, ldk, zero, Arp, ldarp);
+                fgemm (F, FflasNoTrans, FflasNoTrans, Ncurr-Ma, nb_full_blocks, Ma, one,
+                       Ac, ldac, K+(Ncurr-Ma)*ldk, ldk, one, K, ldk);
+                fgemm (F, FflasNoTrans, FflasNoTrans, Ma, nb_full_blocks, Ma, one,
+                       Ac+(Ncurr-Ma)*ldac, ldac, K+(Ncurr-Ma)*ldk, ldk, zero, Arp, ldarp);
                 for (size_t i=0; i< Ma; ++i)
                         fcopy(F, nb_full_blocks, K+(Ncurr-Ma+i)*ldk, 1, Arp+i*ldarp, 1);
-                //write_field(F, cerr<<"K = "<<endl, K, Ncurr, Mk, ldk);
                 
                 // Copying the last rows of A times K
-                //cerr<<"// Copying the last rows of k "<<nb_full_blocks<<" "<<Mk<<endl;
-                 offset = (deg-2)*nb_full_blocks;
+                offset = (deg-2)*nb_full_blocks;
                 for (size_t i = nb_full_blocks; i < Mk; ++i) {
-                        //cerr<<"dK["<<i<<"] = "<<dK[i]<<" deg = "<<deg<<endl;
                         for (size_t j=0; j<Ncurr; ++j)
                                 F.assign(*(K+i+j*ldk), zero);
@@ -552 +303 @@
                         else{
                                 F.assign (*(K + i + (offset+dK[i]-1)*ldk), one);
-                                //cerr<<"*(K + "<<i<<" + ("<<offset <<"+dK["<<i<<"]-1)*ldk) := 1 "<<endl;
                         }
                         offset += dA[i]-1;
                 }
-                //write_field(F, cerr<<"K = "<<endl, K, Ncurr, Mk, ldk);