1 /* Copyright (C) Givaro Team 1999 2 * Copyright (C) LinBox 3 * Written by JG Dumas 4 * 5 * 6 * ========LICENCE======== 7 * This file is part of the library LinBox. 8 * 9 * LinBox is free software: you can redistribute it and/or modify 10 * it under the terms of the GNU Lesser General Public 11 * License as published by the Free Software Foundation; either 12 * version 2.1 of the License, or (at your option) any later version. 13 * 14 * This library is distributed in the hope that it will be useful, 15 * but WITHOUT ANY WARRANTY; without even the implied warranty of 16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 17 * Lesser General Public License for more details. 18 * 19 * You should have received a copy of the GNU Lesser General Public 20 * License along with this library; if not, write to the Free Software 21 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 22 * ========LICENCE======== 23 */ 24 25 26 #ifndef __LINBOX_pp_gauss_H 27 #define __LINBOX_pp_gauss_H 28 29 #include <map> 30 #include <givaro/givconfig.h> // for Signed_Trait 31 #include "linbox/solutions/smith-form.h" 32 #include "linbox/algorithms/gauss.h" 33 34 #ifdef DEBUG 35 # ifndef LINBOX_pp_gauss_intermediate_OUT 36 # define LINBOX_pp_gauss_intermediate_OUT 37 # endif 38 #endif 39 40 // LINBOX_pp_gauss_intermediate_OUT outputs intermediate matrices 41 #ifdef LINBOX_pp_gauss_intermediate_OUT 42 # ifndef LINBOX_pp_gauss_steps_OUT 43 # define LINBOX_pp_gauss_steps_OUT 44 # endif 45 #endif 46 47 // LINBOX_pp_gauss_steps_OUT outputs elimination steps 48 #ifdef LINBOX_pp_gauss_steps_OUT 49 // LINBOX_PRANK_OUT outputs intermediate ranks 50 # ifndef LINBOX_PRANK_OUT 51 # define LINBOX_PRANK_OUT 52 # endif 53 #endif 54 55 56 namespace LinBox 57 { 58 59 template <bool Boolean> struct Boolean_Trait; 60 template <> struct Boolean_Trait<true> { 61 typedef int BooleanType; // int does not matter, only that it differs from float 62 }; 63 template <> struct Boolean_Trait<false> { 64 typedef float BooleanType;// float does not matter, only that it differs from int 65 }; 66 67 enum { 68 // Combine these in binary for use in StaticParameters 69 PRIVILEGIATE_NO_COLUMN_PIVOTING = 1, 70 PRIVILEGIATE_REDUCING_FILLIN = 2, 71 PRESERVE_UPPER_MATRIX = 4 72 }; 73 74 /** \brief Repository of functions for rank modulo 75 * a prime power by elimination on sparse matrices. 76 */ 77 template <class _Field> 78 class PowerGaussDomain : public GaussDomain<_Field> { 79 typedef GaussDomain<_Field> Father_t; 80 typedef _Field Field; 81 typedef typename Field::Element Element; 82 public: 83 84 /** \brief The field parameter is the domain 85 * over which to perform computations 86 */ 87 PowerGaussDomain (const Field &F) : 88 Father_t(F) 89 {} 90 91 //Copy constructor 92 /// 93 PowerGaussDomain (const PowerGaussDomain &M) : 94 Father_t(M) 95 {} 96 97 98 99 // -------------------------------------------- 100 // Modulo operators 101 template<class Modulu> 102 bool isNZero(const Modulu& a ) const { return (bool)a ;} 103 template<class Modulu> 104 bool isZero(const Modulu& a ) const { return a == 0U;} 105 106 template<class Modulo, class Modulo2, class Modulo3> 107 Modulo& MY_Zpz_inv_classic (Modulo& u1, const Modulo2 a, const Modulo3 _p) const 108 { 109 u1 = Modulo(1U); 110 Modulo r0((Modulo)_p), r1((Modulo)a); //! clang complains for examples/smith.C and examples/smithvalence.C 111 Modulo q(r0/r1); 112 113 r0 -= q * r1; 114 if ( isZero(r0) ) return u1; 115 Modulo u0 = q; 116 117 q = r1/r0; 118 r1 -= q * r0; 119 120 while ( isNZero(r1) ) { 121 u1 += q * u0; 122 123 q = r0/r1; 124 r0 -= q * r1; 125 if ( isZero(r0) ) return u1; 126 u0 += q * u1; 127 128 q = r1/r0; 129 r1 -= q * r0; 130 131 } 132 133 return u1=(Modulo)_p-u0; 134 } 135 136 // [On Newton-Raphson iteration for multiplicative 137 // inverses modulo prime powers. J-G. Dumas. 138 // IEEE Trans. on Computers, 63(8), pp 2106-2109, 2014] 139 // http://doi.org/10.1109/TC.2013.94 140 template<class Modulo, class Modulo2, class Modulo3> 141 Modulo& MY_Zpz_inv (Modulo& v1, const Modulo2& a, const Modulo3& _p, const Modulo3& _q, uint64_t exponent) const 142 { 143 Modulo2 ta(a%_p); 144 MY_Zpz_inv_classic(v1, ta, _p); // b=1/a % p 145 Modulo3 ttep2(_q); ttep2 += 2U; // q+2 146 Modulo2 atimeb(a); atimeb*= v1; atimeb %= _q; // ab 147 ttep2 -= atimeb; // 2-ab 148 v1 *= ttep2; v1 %= _q; // b(2-ab) 149 Modulo abmone(atimeb); --abmone; // ab-1 150 for(size_t i=2; i<exponent; i<<=1) { 151 Modulo tmp(abmone); abmone *= tmp; // (ab-1)^(2^i) 152 abmone %= _q; 153 v1 *= ++abmone; --abmone; // in place 154 v1 %= _q; 155 } 156 return v1; 157 } 158 159 template<class Ring1, class Ring2> 160 bool MY_divides(Ring1 a, Ring2 b) const 161 { 162 return (!(b%a)); 163 } 164 165 // ------------------------------------------------ 166 // Pivot Searchers and column strategy 167 // ------------------------------------------------ 168 template<class Modulo, class Vecteur, class D> 169 void SameColumnPivoting(Modulo PRIME, const Vecteur& lignepivot, size_t& indcol, long& indpermut, D& columns, Boolean_Trait<false>::BooleanType ) {} 170 171 172 template<class Modulo, class Vecteur, class D> 173 void SameColumnPivoting(Modulo PRIME, const Vecteur& lignepivot, size_t& indcol, long& indpermut, D& columns, Boolean_Trait<true>::BooleanType ) { 174 // Try first in the same column 175 size_t nj = (size_t) lignepivot.size() ; 176 if (nj && (indcol == lignepivot[0].first) && (! this->MY_divides(PRIME,lignepivot[0].second) ) ) { 177 indpermut = (long) indcol; 178 for(size_t j=nj;j--;) 179 --columns[ lignepivot[(size_t)j].first ]; 180 ++indcol; 181 } 182 } 183 184 template<class Modulo, class BB, class Mmap, class D> 185 bool SameColumnPivotingTrait(Modulo PRIME, size_t& p, const BB& LigneA, const Mmap& psizes, size_t& indcol, long& indpermut, D& columns, Boolean_Trait<false>::BooleanType ) { 186 // Do not try first in the same column 187 return false; 188 } 189 190 template<class Modulo, class BB, class Mmap, class D> 191 bool SameColumnPivotingTrait(Modulo PRIME, size_t& p, const BB& LigneA, const Mmap& psizes, size_t& indcol, long& c, D& columns, Boolean_Trait<true>::BooleanType truetrait) { 192 c=-2; 193 for( typename Mmap::const_iterator iter = psizes.begin(); iter != psizes.end(); ++iter) { 194 p = (size_t)(*iter).second; 195 SameColumnPivoting(PRIME, LigneA[(size_t)p], indcol, c, columns, truetrait ) ; 196 if (c > -2 ) break; 197 } 198 if (c > -2) 199 return true; 200 else 201 return false; 202 203 } 204 205 template<class Modulo, class Vecteur, class D> 206 void CherchePivot(Modulo PRIME, Vecteur& lignepivot, size_t& indcol , long& indpermut, D& columns ) { 207 long nj =(long) lignepivot.size() ; 208 if (nj) { 209 indpermut = (long)lignepivot[0].first; 210 long pp=0; 211 for(;pp<nj;++pp) 212 if (! this->MY_divides(PRIME,lignepivot[(size_t)pp].second) ) break; 213 214 if (pp < nj) { 215 216 long ds = (long)columns[ lignepivot[(size_t)pp].first ],p=pp,j=pp; 217 for(++j;j<nj;++j){ 218 long dl; 219 if ( ( (dl=(long)columns[(size_t)lignepivot[(size_t)j].first] ) < ds ) && (! MY_divides(PRIME,lignepivot[(size_t)j].second) ) ) { 220 ds = dl; 221 p = j; 222 } 223 } 224 if (p != 0) { 225 if (indpermut == (long)indcol) { 226 auto ttm = lignepivot[(size_t)p].second; 227 indpermut = (long)lignepivot[(size_t)p].first; 228 lignepivot[(size_t)p].second = (lignepivot[0].second); 229 lignepivot[0].second = (ttm); 230 } 231 else { 232 auto ttm = lignepivot[(size_t)p]; 233 indpermut = (long)ttm.first; 234 for(long m=p;m;--m) 235 lignepivot[(size_t)m] = lignepivot[(size_t)m-1]; 236 lignepivot[0] = ttm; 237 } 238 } 239 for(j=nj;j--;) 240 --columns[ lignepivot[(size_t)j].first ]; 241 if (indpermut != (long)indcol) { 242 #ifdef LINBOX_pp_gauss_steps_OUT 243 std::cerr << "------CP---- permuting cols " << indcol << " and " << indpermut << " ---" << std::endl; 244 #endif 245 246 lignepivot[0].first = (indcol); 247 } 248 indcol++ ; 249 } 250 else 251 indpermut = -2; 252 } 253 else 254 indpermut = -1; 255 } 256 257 template<class Vecteur> 258 void PermuteColumn(Vecteur& lignecourante, 259 const size_t& nj, 260 const size_t& k, 261 const long& indpermut) { 262 263 REQUIRE( nj > 0 ); 264 REQUIRE( indpermut > (long)k ); 265 266 // Find first non-zero element whose index is 267 // greater than required permutation, if it exists 268 size_t j_head(0); 269 for(; j_head<nj; ++j_head) 270 if (long(lignecourante[(size_t)j_head].first) >= indpermut) break; 271 if ((j_head<nj) && (long(lignecourante[(size_t)j_head].first) == indpermut)) { 272 size_t l(0); 273 for(; l<j_head; ++l) 274 if (lignecourante[(size_t)l].first >= k) break; 275 // ------------------------------------------- 276 // Permutation 277 if (l<j_head && lignecourante[l].first == k) { 278 // non zero <--> non zero 279 auto tmp = lignecourante[l].second ; 280 lignecourante[l].second = (lignecourante[(size_t)j_head].second ); 281 lignecourante[(size_t)j_head].second = (tmp); 282 } else { 283 // zero <--> non zero 284 auto tmp = lignecourante[(size_t)j_head]; 285 tmp.first = (k); 286 for(size_t ll=j_head; ll>l; ll--) 287 lignecourante[ll] = lignecourante[ll-1]; 288 lignecourante[l] = tmp; 289 } 290 } else { 291 // ------------------------------------------- 292 // Permutation 293 size_t l(0); 294 for(; l<nj; ++l) 295 if (lignecourante[(size_t)l].first >= k) break; 296 if ((l<nj) && (lignecourante[(size_t)l].first == k)) { 297 // non zero <--> zero 298 auto tmp = lignecourante[(size_t)l]; 299 tmp.first = (size_t) (indpermut); 300 const size_t bjh(j_head-1); // Position before j_head 301 for(;l<bjh;l++) 302 lignecourante[(size_t)l] = lignecourante[(size_t)l+1]; 303 lignecourante[bjh] = tmp; 304 } // else 305 // zero <--> zero 306 } 307 } 308 309 310 template<class SpMat> 311 void PermuteSubMatrix(SpMat& LigneA, 312 const size_t & start, const size_t & stop, 313 const size_t & currentrank, const long & c) { 314 for(size_t l=start; l < stop; ++l) 315 if ( LigneA[(size_t)l].size() ) 316 PermuteColumn(LigneA[(size_t)l], LigneA[(size_t)l].size(), currentrank, c); 317 } 318 319 template<class SpMat> 320 void PermuteUpperMatrix(SpMat& LigneA, const size_t & k, const size_t & currentrank, const long & c, Boolean_Trait<true>::BooleanType ) { 321 PermuteSubMatrix(LigneA, 0, k, currentrank, c); 322 } 323 324 template<class SpMat> 325 void PermuteUpperMatrix(SpMat&, const size_t &, const size_t &, const long &, Boolean_Trait<false>::BooleanType ) {} 326 327 328 template<class Vecteur> 329 void PreserveUpperMatrixRow(Vecteur& ligne, Boolean_Trait<true>::BooleanType ) {} 330 331 template<class Vecteur> 332 void PreserveUpperMatrixRow(Vecteur& ligne, Boolean_Trait<false>::BooleanType ) { 333 ligne = Vecteur(0); 334 } 335 336 337 template<class Modulo, class Vecteur, class De> 338 void FaireElimination( Modulo MOD, 339 Vecteur& lignecourante, 340 const Vecteur& lignepivot, 341 const typename Signed_Trait<Modulo>::unsigned_type& invpiv, 342 const size_t& k, 343 const long& indpermut, 344 De& columns) { 345 346 // typedef typename Vecteur::coefficientSpace F; 347 // typedef typename Vecteur::value_types E; 348 typedef typename Field::Element F; 349 typedef typename Vecteur::value_type E; 350 351 typedef typename Signed_Trait<Modulo>::unsigned_type UModulo; 352 353 size_t nj = (size_t) lignecourante.size() ; 354 if (nj) { 355 if (lignecourante[0].first == k) { 356 // ------------------------------------------- 357 // Elimination 358 size_t npiv = (size_t) lignepivot.size(); 359 Vecteur construit(nj + npiv); 360 // construit : <-- ci 361 // courante : <-- m 362 // pivot : <-- l 363 typedef typename Vecteur::iterator Viter; 364 Viter ci = construit.begin(); 365 size_t m=1; 366 size_t l(0); 367 368 F headcoeff = MOD-(lignecourante[0].second); 369 headcoeff *= invpiv; 370 headcoeff %= (UModulo)MOD ; 371 lignecourante[0].second = headcoeff; 372 --columns[ lignecourante[0].first ]; 373 374 for(;l<npiv;++l) 375 if (lignepivot[(size_t)l].first > k) break; 376 // for all j such that (j>k) and A[(size_t)k,j]!=0 377 for(;l<npiv;++l) { 378 size_t j_piv; 379 j_piv = (size_t) lignepivot[(size_t)l].first; 380 // if A[(size_t)k,j]=0, then A[(size_t)i,j] <-- A[(size_t)i,j] 381 for (;(m<nj) && (lignecourante[(size_t)m].first < j_piv);) 382 *ci++ = lignecourante[(size_t)m++]; 383 // if A[(size_t)i,j]!=0, then A[(size_t)i,j] <-- A[(size_t)i,j] - A[(size_t)i,k]*A[(size_t)k,j] 384 if ((m<nj) && (lignecourante[(size_t)m].first == j_piv)) { 385 //lignecourante[(size_t)m].second = ( ((UModulo)( headcoeff * lignepivot[(size_t)l].second + lignecourante[(size_t)m].second ) ) % (UModulo)MOD ); 386 lignecourante[(size_t)m].second += ( headcoeff * lignepivot[(size_t)l].second ); 387 lignecourante[(size_t)m].second %= (UModulo)MOD; 388 if (isNZero(lignecourante[(size_t)m].second)) 389 *ci++ = lignecourante[(size_t)m++]; 390 else 391 --columns[ lignecourante[(size_t)m++].first ]; 392 // m++; 393 } 394 else { 395 F tmp(headcoeff); 396 tmp *= lignepivot[(size_t)l].second; 397 tmp %= (UModulo)MOD; 398 if (isNZero(tmp)) { 399 ++columns[(size_t)j_piv]; 400 *ci++ = E(j_piv, tmp ); 401 } 402 } 403 } 404 // if A[(size_t)k,j]=0, then A[(size_t)i,j] <-- A[(size_t)i,j] 405 for (;m<nj;) 406 *ci++ = lignecourante[(size_t)m++]; 407 408 construit.erase(ci,construit.end()); 409 lignecourante = construit; 410 } 411 } 412 } 413 414 // ------------------------------------------------------ 415 // Rank calculators, defining row strategy 416 // ------------------------------------------------------ 417 418 template<class Modulo, class BB, class D, class Container, class Perm, bool PrivilegiateNoColumnPivoting, bool PreserveUpperMatrix> 419 void gauss_rankin(Modulo FMOD, Modulo PRIME, Container& ranks, BB& LigneA, Perm& Q, const size_t Ni, const size_t Nj, const D& density_trait) 420 { 421 linbox_check( Q.coldim() == Q.rowdim() ); 422 linbox_check( Q.coldim() == Nj ); 423 424 commentator().start ("Gaussian elimination with reordering modulo a prime power", 425 "PRGE", Ni); 426 427 ranks.resize(0); 428 429 typedef typename BB::Row Vecteur; 430 typedef typename Signed_Trait<Modulo>::unsigned_type UModulo; 431 432 Modulo MOD = FMOD; 433 uint64_t exponent(1); 434 Modulo tq(FMOD); 435 while(tq > PRIME) { 436 tq /= PRIME; 437 ++exponent; 438 } 439 #ifdef LINBOX_PRANK_OUT 440 std::cerr << "Elimination mod " << MOD << " (=" << PRIME << '^' << exponent << ',' << PrivilegiateNoColumnPivoting << ',' << PreserveUpperMatrix << ')' << std::endl; 441 #endif 442 443 D col_density(Nj); 444 445 // assignment of LigneA with the domain object 446 size_t jj; 447 for(jj=0; jj<Ni; ++jj) { 448 Vecteur tmp = LigneA[(size_t)jj]; 449 Vecteur toto(tmp.size()); 450 size_t k=0,rs=0; 451 for(; k<tmp.size(); ++k) { 452 Modulo r = tmp[(size_t)k].second; 453 if ((r <0) || (r >= MOD)) r %= MOD ; 454 if (r <0) r += MOD ; 455 if (isNZero(r)) { 456 ++col_density[ tmp[(size_t)k].first ]; 457 toto[rs] =tmp[(size_t)k]; 458 toto[rs].second = ( r ); 459 ++rs; 460 } 461 } 462 toto.resize(rs); 463 LigneA[(size_t)jj] = toto; 464 // LigneA[(size_t)jj].reactualsize(Nj); 465 466 } 467 468 size_t last = Ni-1; 469 long c(0); 470 size_t indcol(0); 471 size_t ind_pow = 1; 472 size_t maxout = Ni/100; maxout = (maxout<10 ? 10 : (maxout>1000 ? 1000 : maxout) ); 473 size_t thres = Ni/maxout; thres = (thres >0 ? thres : 1); 474 475 476 for (size_t k=0; k<last;++k) { 477 if ( ! (k % maxout) ) commentator().progress ((long)k); 478 479 480 size_t p=k; 481 for(;;) { 482 483 484 std::multimap< long, long > psizes; 485 for(p=k; p<Ni; ++p) 486 psizes.insert( psizes.end(), std::pair<long,long>( (long) LigneA[(size_t)p].size(), (long) p) ); 487 488 #ifdef LINBOX_pp_gauss_steps_OUT 489 std::cerr << "------------ ordered rows " << k << " -----------" << std::endl; 490 for(auto const& iter: psizes) 491 std::cerr << iter.second << " : #" << iter.first << std::endl; 492 std::cerr << "---------------------------------------" << std::endl; 493 #endif 494 495 if ( SameColumnPivotingTrait(PRIME, p, LigneA, psizes, indcol, c, col_density, typename Boolean_Trait<PrivilegiateNoColumnPivoting>::BooleanType() ) ) 496 break; 497 498 for( typename std::multimap< long, long >::const_iterator iter = psizes.begin(); iter != psizes.end(); ++iter) { 499 p = (size_t)(*iter).second; 500 501 CherchePivot( PRIME, LigneA[(size_t)p], indcol, c , col_density) ; 502 if (c > -2 ) break; 503 } 504 505 if (c > -2) break; 506 for(size_t ii=k;ii<Ni;++ii) 507 for(size_t jjj=LigneA[(size_t)ii].size();jjj--;) 508 LigneA[(size_t)ii][(size_t)jjj].second /= PRIME; 509 MOD /= PRIME; 510 --exponent; 511 512 #ifdef LINBOX_PRANK_OUT 513 std::cerr << "Rank mod " << PRIME << "^" << ind_pow << " : " << indcol << std::endl; 514 if (MOD == 1) std::cerr << "wattadayada inhere ?" << std::endl; 515 #endif 516 ranks.push_back( indcol ); 517 ++ind_pow; 518 519 } 520 if (p != k) { 521 #ifdef LINBOX_pp_gauss_steps_OUT 522 std::cerr << "------------ permuting rows " << p << " and " << k << " ---" << std::endl; 523 #endif 524 Vecteur vtm = LigneA[(size_t)k]; 525 LigneA[(size_t)k] = LigneA[(size_t)p]; 526 LigneA[(size_t)p] = vtm; 527 } 528 if (c != -1) { 529 REQUIRE( indcol > 0); 530 const size_t currentrank(indcol-1); 531 if (c != (long)currentrank) { 532 #ifdef LINBOX_pp_gauss_steps_OUT 533 std::cerr << "------------ permuting cols " << (indcol-1) << " and " << c << " ---" << std::endl; 534 #endif 535 Q.permute(indcol-1,c); 536 std::swap(col_density[currentrank],col_density[c]); 537 PermuteUpperMatrix(LigneA, k, currentrank, c, typename Boolean_Trait<PreserveUpperMatrix>::BooleanType()); 538 PermuteSubMatrix(LigneA, k+1, Ni, currentrank, c); 539 } 540 541 UModulo invpiv; 542 MY_Zpz_inv(invpiv, LigneA[(size_t)k][0].second, PRIME, MOD, exponent); 543 544 for(size_t l=k + 1; (l < Ni) && (col_density[currentrank]); ++l) 545 FaireElimination(MOD, LigneA[(size_t)l], LigneA[(size_t)k], invpiv, currentrank, c, col_density); 546 } 547 548 549 #ifdef LINBOX_pp_gauss_steps_OUT 550 LigneA.write(std::cerr << "step[" << k << "], pivot: " << c << std::endl, Tag::FileFormat::Maple) << std::endl; 551 #endif 552 553 PreserveUpperMatrixRow(LigneA[(size_t)k], typename Boolean_Trait<PreserveUpperMatrix>::BooleanType()); 554 } 555 556 c = -2; 557 SameColumnPivoting(PRIME, LigneA[(size_t)last], indcol, c, col_density, typename Boolean_Trait<PrivilegiateNoColumnPivoting>::BooleanType() ); 558 if (c == -2) CherchePivot( PRIME, LigneA[(size_t)last], indcol, c, col_density ); 559 while( c == -2) { 560 ranks.push_back( indcol ); 561 for(long jjj=(long)LigneA[(size_t)last].size();jjj--;) 562 LigneA[(size_t)last][(size_t)jjj].second /= PRIME; 563 MOD /= PRIME; 564 CherchePivot( PRIME, LigneA[(size_t)last], indcol, c, col_density ); 565 } 566 if (c != -1) { 567 const size_t currentrank(indcol-1); 568 if (c != (long)currentrank) { 569 #ifdef LINBOX_pp_gauss_steps_OUT 570 std::cerr << "------------ permuting cols " << (indcol-1) << " and " << c << " ---" << std::endl; 571 #endif 572 Q.permute(indcol-1,c); 573 PermuteUpperMatrix(LigneA, last, currentrank, c, typename Boolean_Trait<PreserveUpperMatrix>::BooleanType()); 574 } 575 } 576 while( MOD > 1) { 577 MOD /= PRIME; 578 ranks.push_back( indcol ); 579 } 580 581 #ifdef LINBOX_pp_gauss_steps_OUT 582 LigneA.write(std::cerr << "step[" << Ni-1 << "], pivot: " << c << std::endl, Tag::FileFormat::Maple) << std::endl; 583 #endif 584 #ifdef LINBOX_PRANK_OUT 585 std::cerr << "Rank mod " << FMOD << " : " << indcol << std::endl; 586 #endif 587 commentator().stop ("done", 0, "PRGE"); 588 589 } 590 591 template<class Modulo, class BB, class D, class Container, class Perm> 592 void prime_power_rankin (Modulo FMOD, Modulo PRIME, Container& ranks, BB& SLA, Perm& Q, const size_t Ni, const size_t Nj, const D& density_trait, int StaticParameters=PRIVILEGIATE_NO_COLUMN_PIVOTING) 593 { 594 if (PRIVILEGIATE_NO_COLUMN_PIVOTING & StaticParameters) { 595 if (PRESERVE_UPPER_MATRIX & StaticParameters) { 596 gauss_rankin<Modulo,BB,D,Container,Perm,true,true>(FMOD,PRIME,ranks, SLA, Q, Ni, Nj, density_trait); 597 } else { 598 gauss_rankin<Modulo,BB,D,Container,Perm,true,false>(FMOD,PRIME,ranks, SLA, Q, Ni, Nj, density_trait); 599 } 600 } else { 601 if (PRESERVE_UPPER_MATRIX & StaticParameters) { 602 gauss_rankin<Modulo,BB,D,Container,Perm,false,true>(FMOD,PRIME,ranks, SLA, Q, Ni, Nj, density_trait); 603 } else { 604 gauss_rankin<Modulo,BB,D,Container,Perm,false,false>(FMOD,PRIME,ranks, SLA, Q, Ni, Nj, density_trait); 605 } 606 } 607 } 608 609 610 template<class Modulo, class Matrix, class Perm, template<class, class> class Container, template<class> class Alloc> 611 Container<std::pair<Modulo,size_t>, Alloc<std::pair<Modulo,size_t> > >& operator()(Container<std::pair<Modulo,size_t>, Alloc<std::pair<Modulo,size_t> > >& L, Matrix& A, Perm& Q, Modulo FMOD, Modulo PRIME, int StaticParameters=PRIVILEGIATE_NO_COLUMN_PIVOTING) 612 { 613 Container<size_t, Alloc<size_t> > ranks; 614 prime_power_rankin( FMOD, PRIME, ranks, A, Q, A.rowdim(), A.coldim(), std::vector<size_t>(), StaticParameters); 615 L.resize( 0 ) ; 616 Modulo MOD = 1; 617 size_t num = 0; 618 for( typename Container<size_t, Alloc<size_t> >::const_iterator it = ranks.begin(); it != ranks.end(); ++it) { 619 size_t diff(*it-num); 620 if (diff > 0) 621 L.emplace_back(MOD,diff); 622 MOD *= PRIME; 623 num = *it; 624 } 625 return L; 626 } 627 628 }; 629 630 631 } // end of LinBox namespace 632 633 #endif //__LINBOX_pp_gauss_H 634 635 // Local Variables: 636 // mode: C++ 637 // tab-width: 4 638 // indent-tabs-mode: nil 639 // c-basic-offset: 4 640 // End: 641 // vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s 642