1 /*
2 * algorithms/smith-form-valence.h
3 * Copyright (c) Linbox
4 * ========LICENCE========
5 * This file is part of the library LinBox.
6 *
7 * LinBox is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU Lesser General Public
9 * License as published by the Free Software Foundation; either
10 * version 2.1 of the License, or (at your option) any later version.
11 *
12 * This library is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 * Lesser General Public License for more details.
16 *
17 * You should have received a copy of the GNU Lesser General Public
18 * License along with this library; if not, write to the Free Software
19 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
20 * ========LICENCE========
21 */
22
23 /**\file algorithms/smith-form-valence.h
24 * @example algorithms/smith-form-valence.h
25 \brief Valence of sparse matrix over Z or Zp.
26 \ingroup examples
27 */
28
29 #include <string>
30 #include <givaro/modular.h>
31 #include <givaro/givintnumtheo.h>
32
33 #include <givaro/gf2.h>
34 #include <linbox/field/field-traits.h>
35 #include <linbox/blackbox/transpose.h>
36 #include <linbox/blackbox/compose.h>
37 #include <linbox/matrix/sparse-matrix.h>
38 #include <linbox/solutions/rank.h>
39 #include <linbox/solutions/valence.h>
40 #include <linbox/solutions/smith-form.h>
41 #include <linbox/solutions/valence.h>
42 #include <linbox/algorithms/smith-form-sparseelim-local.h>
43 #include <linbox/util/matrix-stream.h>
44 #include <linbox/util/timer.h>
45 #include <linbox/util/error.h>
46
47 #ifndef __VALENCE_FACTOR_LOOPS__
48 #define __VALENCE_FACTOR_LOOPS__ 50000
49 #endif
50
51 #ifndef __VALENCE_REPORTING__
52 # ifdef _LB_DEBUG
53 # define __VALENCE_REPORTING__ 1
54 # else
55 # define __VALENCE_REPORTING__ 0
56 # endif
57 #endif
58
59 #ifdef __LINBOX_USE_OPENMP
60 #include <omp.h>
61 #define THREAD_NUM omp_get_thread_num()
62 #else
63 #define THREAD_NUM 1
64 #endif
65
66 namespace LinBox {
67
68 template<class Field>
TempLRank(size_t & r,const char * filename,const Field & F)69 size_t& TempLRank(size_t& r, const char * filename, const Field& F)
70 {
71 std::ifstream input(filename);
72 MatrixStream< Field > msf( F, input );
73 SparseMatrix<Field,SparseMatrixFormat::SparseSeq> FA(msf);
74 input.close();
75 Timer tim; tim.start();
76 rankInPlace(r, FA);
77 tim.stop();
78 if (__VALENCE_REPORTING__)
79 F.write(std::clog << "Rank over ") << " is " << r << ' ' << tim << " on T" << THREAD_NUM << std::endl;
80 return r;
81 }
82
TempLRank(size_t & r,const char * filename,const GF2 & F2)83 size_t& TempLRank(size_t& r, const char * filename, const GF2& F2)
84 {
85 std::ifstream input(filename);
86 ZeroOne<GF2> A;
87 A.read(input);
88 input.close();
89
90 Timer tim; tim.start();
91 rankInPlace(r, A, Method::SparseElimination() );
92 tim.stop();
93 if (__VALENCE_REPORTING__)
94 F2.write(std::clog << "Rank over ") << " is " << r << ' ' << tim << " on T" << THREAD_NUM << std::endl;
95 return r;
96 }
97
LRank(size_t & r,const char * filename,Givaro::Integer p)98 size_t& LRank(size_t& r, const char * filename,Givaro::Integer p)
99 {
100
101 Givaro::Integer maxmod16; FieldTraits<Givaro::Modular<int16_t> >::maxModulus(maxmod16);
102 Givaro::Integer maxmod32; FieldTraits<Givaro::Modular<int32_t> >::maxModulus(maxmod32);
103 Givaro::Integer maxmod53; FieldTraits<Givaro::Modular<double> >::maxModulus(maxmod53);
104 Givaro::Integer maxmod64; FieldTraits<Givaro::Modular<int64_t> >::maxModulus(maxmod64);
105 if (p == 2) {
106 GF2 F2;
107 return TempLRank(r, filename, F2);
108 }
109 else if (p <= maxmod16) {
110 typedef Givaro::Modular<int16_t> Field;
111 Field F(p);
112 return TempLRank(r, filename, F);
113 }
114 else if (p <= maxmod32) {
115 typedef Givaro::Modular<int32_t> Field;
116 Field F(p);
117 return TempLRank(r, filename, F);
118 }
119 else if (p <= maxmod53) {
120 typedef Givaro::Modular<double> Field;
121 Field F(p);
122 return TempLRank(r, filename, F);
123 }
124 else if (p <= maxmod64) {
125 typedef Givaro::Modular<int64_t> Field;
126 Field F(p);
127 return TempLRank(r, filename, F);
128 }
129 else {
130 typedef Givaro::Modular<Givaro::Integer> Field;
131 Field F(p);
132 return TempLRank(r, filename, F);
133 }
134 return r;
135 }
136
PRank(std::vector<size_t> & ranks,size_t & effective_exponent,const char * filename,Givaro::Integer p,size_t e,size_t intr)137 std::vector<size_t>& PRank(std::vector<size_t>& ranks, size_t& effective_exponent, const char * filename,Givaro::Integer p, size_t e, size_t intr)
138 {
139 effective_exponent = e;
140 Givaro::Integer maxmod;
141 FieldTraits<Givaro::Modular<int64_t> >::maxModulus(maxmod);
142 if (p <= maxmod) {
143 typedef Givaro::Modular<int64_t> Ring;
144 int64_t lp(p);
145 Givaro::Integer q = pow(p,uint64_t(e)); int64_t lq(q);
146 if (q > Ring::maxCardinality()) {
147 if (__VALENCE_REPORTING__)
148 std::clog << "Power rank might need extra large composite (" << p << '^' << e << ")." << std::endl;
149 q = p;
150 for(effective_exponent=1; q <= Ring::maxCardinality(); ++effective_exponent) {
151 q *= p;
152 }
153 q/=p; --effective_exponent;
154 lq = (int64_t)q;
155 if (effective_exponent <= 1) {
156 // Not able to use Modular<int64_t>,
157 // modulus is ok, but is already too large when squared
158 // Return that no prime power was performed
159 if (__VALENCE_REPORTING__)
160 std::clog << "Exceeding int64_t ... power rank useless, nothing done." << std::endl;
161 effective_exponent=1;
162 return ranks;
163 }
164 if (__VALENCE_REPORTING__)
165 std::clog << "First trying: " << lq << " (=" << p << '^' << effective_exponent << ", without further warning this will be sufficient)." << std::endl;
166 }
167 Ring F(lq);
168 std::ifstream input(filename);
169 MatrixStream<Ring> ms( F, input );
170 SparseMatrix<Ring,SparseMatrixFormat::SparseSeq > A (ms);
171 input.close();
172
173 PowerGaussDomain< Ring > PGD( F );
174 Permutation<Ring> Q(F,A.coldim());
175
176 Timer tim; tim.clear(); tim.start();
177 PGD.prime_power_rankin( lq, lp, ranks, A, Q, A.rowdim(), A.coldim(), std::vector<size_t>());
178 tim.stop();
179 if (__VALENCE_REPORTING__) {
180 F.write(std::clog << "Ranks over ") << " are " ;
181 for(auto const & rit:ranks) std::clog << rit << ' ';
182 std::clog << ' ' << tim << " on T" << THREAD_NUM << std::endl;
183 }
184 } else {
185 // Not able to use Modular<int64_t>, modulus too large
186 // Return that no prime power was performed
187 if (__VALENCE_REPORTING__)
188 std::clog << "Exceeding int64_t ... even for the prime, nothing done." << std::endl;
189 effective_exponent=1;
190 }
191 return ranks;
192 }
193 }
194
195 #include <linbox/field/gf2.h>
196 #include <linbox/algorithms/smith-form-sparseelim-poweroftwo.h>
197
198 namespace LinBox {
199
PRankPowerOfTwo(std::vector<size_t> & ranks,size_t & effective_exponent,const char * filename,size_t e,size_t intr)200 std::vector<size_t>& PRankPowerOfTwo(std::vector<size_t>& ranks, size_t& effective_exponent, const char * filename, size_t e, size_t intr)
201 {
202 effective_exponent = e;
203 if (e > 63) {
204 if (__VALENCE_REPORTING__)
205 std::clog << "Power rank power of two might need extra large composite (2^" << e << ")." << std::endl;
206 if (__VALENCE_REPORTING__)
207 std::clog << "First trying: 63, without further warning this will be sufficient)." << std::endl;
208 effective_exponent = 63;
209 }
210
211
212 typedef Givaro::ZRing<int64_t> Ring;
213 Ring F;
214 std::ifstream input(filename);
215 MatrixStream<Ring> ms( F, input );
216 SparseMatrix<Ring,SparseMatrixFormat::SparseSeq > A (ms);
217 input.close();
218 PowerGaussDomainPowerOfTwo< uint64_t > PGD;
219 GF2 F2;
220 Permutation<GF2> Q(F2,A.coldim());
221
222 Timer tim; tim.clear(); tim.start();
223 PGD.prime_power_rankin( effective_exponent, ranks, A, Q, A.rowdim(), A.coldim(), std::vector<size_t>());
224 tim.stop();
225 if (__VALENCE_REPORTING__) {
226 F.write(std::clog << "Ranks over ") << " modulo 2^" << effective_exponent << " are " ;
227 for(auto const& rit: ranks) std::clog << rit << ' ';
228 std::clog << ' ' << tim << " on T" << THREAD_NUM << std::endl;
229 }
230 return ranks;
231 }
232
PRankInteger(std::vector<size_t> & ranks,const char * filename,Givaro::Integer p,size_t e,size_t intr)233 std::vector<size_t>& PRankInteger(std::vector<size_t>& ranks, const char * filename,Givaro::Integer p, size_t e, size_t intr)
234 {
235 typedef Givaro::Modular<Givaro::Integer> Ring;
236 Givaro::Integer q = pow(p,uint64_t(e));
237 Ring F(q);
238 std::ifstream input(filename);
239 MatrixStream<Ring> ms( F, input );
240 SparseMatrix<Ring,SparseMatrixFormat::SparseSeq > A (ms);
241 input.close();
242 PowerGaussDomain< Ring > PGD( F );
243 Permutation<Ring> Q(F,A.coldim());
244
245 Timer tim; tim.clear(); tim.start();
246 PGD.prime_power_rankin( q, p, ranks, A, Q, A.rowdim(), A.coldim(), std::vector<size_t>());
247 tim.stop();
248 if (__VALENCE_REPORTING__) {
249 F.write(std::clog << "Ranks over ") << " are " ;
250 for(auto const & rit: ranks) std::clog << rit << ' ';
251 std::clog << ' ' << tim << std::endl;
252 }
253 return ranks;
254 }
255
PRankIntegerPowerOfTwo(std::vector<size_t> & ranks,const char * filename,size_t e,size_t intr)256 std::vector<size_t>& PRankIntegerPowerOfTwo(std::vector<size_t>& ranks, const char * filename, size_t e, size_t intr)
257 {
258 typedef Givaro::ZRing<Givaro::Integer> Ring;
259 Ring ZZ;
260 std::ifstream input(filename);
261 MatrixStream<Ring> ms( ZZ, input );
262 SparseMatrix<Ring,SparseMatrixFormat::SparseSeq > A (ms);
263 input.close();
264 PowerGaussDomainPowerOfTwo< Givaro::Integer > PGD;
265 Permutation<Ring> Q(ZZ, A.coldim());
266
267 Timer tim; tim.clear(); tim.start();
268 PGD.prime_power_rankin( e, ranks, A, Q, A.rowdim(), A.coldim(), std::vector<size_t>());
269 tim.stop();
270 if (__VALENCE_REPORTING__) {
271 ZZ.write(std::clog << "Ranks over ") << " modulo 2^" << e << " are " ;
272 for(auto const& rit: ranks) std::clog << rit << ' ';
273 std::clog << ' ' << tim << std::endl;
274 }
275 return ranks;
276 }
277
278
279 typedef std::pair<Givaro::Integer,size_t> PairIntRk;
280
281
AllPowersRanks(std::vector<size_t> & ranks,const Givaro::Integer & squarefreePrime,const size_t & squarefreeRank,const size_t & exponentBound,const size_t & coprimeRank,const char * filename)282 std::vector<size_t>& AllPowersRanks(
283 std::vector<size_t>& ranks,
284 const Givaro::Integer& squarefreePrime,// smith[j].first
285 const size_t& squarefreeRank,// smith[j].second
286 const size_t& exponentBound, // exponents[j]
287 const size_t& coprimeRank, // coprimeR
288 const char * filename) { // argv[1]
289
290 if (squarefreeRank != coprimeRank) {
291
292 ranks.push_back(squarefreeRank);
293 size_t effexp;
294 if (exponentBound > 1) {
295 // See if a not too small, not too large exponent would work
296 // Usually, closest to word size
297 if (squarefreePrime == 2)
298 PRankPowerOfTwo(ranks, effexp, filename, exponentBound, coprimeRank);
299 else
300 PRank(ranks, effexp, filename, squarefreePrime, exponentBound, coprimeRank);
301 } else {
302 // Square does not divide valence
303 // Try first with the smallest possible exponent: 2
304 if (squarefreePrime == 2)
305 PRankPowerOfTwo(ranks, effexp, filename, 2, coprimeRank);
306 else
307 PRank(ranks, effexp, filename, squarefreePrime, 2, coprimeRank);
308 }
309
310 if (effexp < exponentBound) {
311 // Above report shows that more powers are needed,
312 // try successive doublings Over abitrary precision
313 for(size_t expo = effexp<<1; ranks.back() < coprimeRank; expo<<=1) {
314 if (squarefreePrime == 2)
315 PRankIntegerPowerOfTwo(ranks, filename, expo, coprimeRank);
316 else
317 PRankInteger(ranks, filename, squarefreePrime, expo, coprimeRank);
318 }
319 } else {
320 // Larger exponents are needed
321 // Try first small precision, then arbitrary
322 for(size_t expo = (exponentBound)<<1; ranks.back() < coprimeRank; expo<<=1) {
323 if (squarefreePrime == 2)
324 PRankPowerOfTwo(ranks, effexp, filename, expo, coprimeRank);
325 else
326 PRank(ranks, effexp, filename, squarefreePrime, expo, coprimeRank);
327 if (ranks.size() < expo) {
328 if (__VALENCE_REPORTING__)
329 std::clog << "It seems we need a larger prime power, it will take longer ..." << std::endl;
330 // break;
331 if (squarefreePrime == 2)
332 PRankIntegerPowerOfTwo(ranks, filename, expo, coprimeRank);
333 else
334 PRankInteger(ranks, filename, squarefreePrime, expo, coprimeRank);
335 }
336 }
337 }
338 }
339
340 return ranks;
341 }
342
populateSmithForm(std::vector<Givaro::Integer> & SmithDiagonal,const std::vector<size_t> & ranks,const Givaro::Integer & squarefreePrime,const size_t & squarefreeRank,const size_t & coprimeRank)343 std::vector<Givaro::Integer>& populateSmithForm(
344 std::vector<Givaro::Integer>& SmithDiagonal,
345 const std::vector<size_t>& ranks,
346 const Givaro::Integer& squarefreePrime,// smith[j].first
347 const size_t& squarefreeRank,// smith[j].second
348 const size_t& coprimeRank) { // coprimeR
349
350
351 SmithDiagonal.resize(coprimeRank, Givaro::Integer(1) );
352
353 for(size_t i=squarefreeRank; i < coprimeRank; ++i) {
354 SmithDiagonal[i] *= squarefreePrime;
355 }
356 auto rit=ranks.begin(); for(++rit; rit!= ranks.end(); ++rit) {
357 if ((*rit)>= coprimeRank) break;
358 for(size_t i=(*rit); i < coprimeRank; ++i)
359 SmithDiagonal[i] *= squarefreePrime;
360 }
361
362 return SmithDiagonal;
363 }
364
365 template<class Blackbox>
366 std::vector<Givaro::Integer>& smithValence(std::vector<Givaro::Integer>& SmithDiagonal,
367 Givaro::Integer& valence,
368 const Blackbox& A,
369 const std::string& filename,
370 Givaro::Integer& coprimeV,
371 size_t method=0) {
372 // method for valence squarization:
373 // 0 for automatic, 1 for aat, 2 for ata
374 // Blackbox provides the Integer matrix rereadable from filename
375 // if valence != 0:
376 // then the valence is not computed and the parameter is used
377 // if coprimeV != 1:
378 // then this value is supposed to be coprime with the valence
379
380 if (__VALENCE_REPORTING__)
381 std::clog << "sV threads: " << NUM_THREADS << std::endl;
382
383 if (valence == 0) {
384 squarizeValence(valence, A, method);
385 }
386
387 if (__VALENCE_REPORTING__)
388 std::clog << "Valence is " << valence << std::endl;
389
390 std::vector<Givaro::Integer> Moduli;
391 std::vector<size_t> exponents;
392
393 if (__VALENCE_REPORTING__)
394 std::clog << "Some factors (" << __VALENCE_FACTOR_LOOPS__ << " factoring loop bound): ";
395
396 Givaro::IntFactorDom<> FTD;
397 FTD.set(Moduli, exponents, valence, __VALENCE_FACTOR_LOOPS__);
398
399 if (__VALENCE_REPORTING__) {
400 auto eit=exponents.begin();
401 for(auto const &mit: Moduli) std::clog << mit << '^' << *eit++ << ' ';
402 std::clog << std::endl;
403 }
404
405 std::vector< size_t > smith(Moduli.size());
406
407 if (coprimeV == 1) {
408 coprimeV=2;
409 while ( gcd(valence,coprimeV) > 1 ) {
410 FTD.nextprimein(coprimeV);
411 }
412 }
413
414 size_t coprimeR;
415 std::vector<std::vector<size_t> > AllRanks(Moduli.size());
416
417 for(size_t j=0; j<Moduli.size(); ++j) {
418 { TASK(MODE(CONSTREFERENCE(Moduli,smith,filename) WRITE(smith[j]) ),
419 {
420 LRank(smith[j], filename.c_str(), Moduli[j]);
421 })}
422 }
423
424 // { TASK(MODE(CONSTREFERENCE(coprimeV,filename) WRITE(coprimeR) ),
425 // {
426 LRank(coprimeR, filename.c_str(), coprimeV);
427 // })}
428
429 WAIT;
430
431 SYNCH_GROUP(
432 for(size_t j=0; j<Moduli.size(); ++j) {
433 { TASK(MODE(CONSTREFERENCE(smith,Moduli,AllRanks,filename,coprimeR,exponents)
434 WRITE(AllRanks[j])),
435 {
436 AllPowersRanks(AllRanks[j], Moduli[j], smith[j], exponents[j],
437 coprimeR, filename.c_str());
438 })}
439 }
440 )
441
442 for(size_t j=0; j<Moduli.size(); ++j) {
443 if (smith[j] != coprimeR) {
444 populateSmithForm(SmithDiagonal, AllRanks[j], Moduli[j], smith[j], coprimeR);
445 }
446 }
447
448 return SmithDiagonal;
449 }
450
451 template<class Blackbox>
452 std::vector<Givaro::Integer>& smithValence(
453 std::vector<Givaro::Integer>& SmithDiagonal,
454 const Blackbox& A,
455 const std::string& filename,
456 size_t method=0)
457 {
458 Givaro::Integer valence(0);
459 Givaro::Integer coprimeV(1);
460 return smithValence(SmithDiagonal, valence, A, filename, coprimeV, method);
461 }
462
463
464 template<class PIR>
writeCompressedSmith(std::ostream & out,const std::vector<typename PIR::Element> & SmithDiagonal,const PIR & ZZ,const size_t m,const size_t n)465 std::ostream& writeCompressedSmith(
466 std::ostream& out,
467 const std::vector<typename PIR::Element>& SmithDiagonal,
468 const PIR& ZZ,
469 const size_t m, const size_t n) {
470 // Output is a list of pairs (integral value, number of repetitions)
471
472 SmithList<PIR> tempSL;
473 compressedSmith(tempSL, SmithDiagonal, ZZ, m, n);
474 out << '(';
475 for( auto sit : tempSL )
476 out << '[' << sit.first << ',' << sit.second << "] ";
477 return out << ')';
478 }
479
480 }
481
482
483 // Local Variables:
484 // mode: C++
485 // tab-width: 4
486 // indent-tabs-mode: nil
487 // c-basic-offset: 4
488 // End:
489 // vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
490