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