1 #include "linbox/linbox-config.h"
2
3 #include <iostream>
4 #include <fstream>
5
6 #include <omp.h>
7 #include <time.h>
8 #include <set>
9 #include <utility>
10 #include <sstream>
11
12 #include "linbox/ring/modular.h"
13 #include "linbox/matrix/sparse-matrix.h"
14 #include "linbox/blackbox/transpose.h"
15 #include "linbox/vector/vector-domain.h"
16 #include "linbox/matrix/dense-matrix.h"
17 #include "linbox/matrix/matrix-domain.h"
18
19 #include "linbox/algorithms/blackbox-block-container.h"
20 #include "linbox/algorithms/block-coppersmith-domain.h"
21
22 #include "linbox/solutions/det.h"
23 #include "linbox/solutions/rank.h"
24 #include "linbox/solutions/methods.h"
25
26 #include "linbox/algorithms/wiedemann.h"
27
28 #include "examples/map-sparse.h"
29
30 // Generates random dense matrices U and V and sparse matrix M for use
31 // with block-coppersmith-benchmark and invariant-factors-benchmark.
32 // Also computes the min-poly and saves that as well
33
34 using namespace std;
35 using namespace LinBox;
36
37 typedef Givaro::Modular<double> Field;
38 typedef typename Field::Element Element;
39 typedef SparseMatrix<Field> SparseMat;
40 typedef MatrixDomain<Field> Domain;
41 typedef typename Domain::OwnMatrix Block;
42
randRange(int start,int end)43 int randRange(int start, int end)
44 {
45 double rval = rand();
46 static const double NORMALIZING_CONSTANT = 1.0/(1.0+RAND_MAX);
47 double normedRVal = rval*NORMALIZING_CONSTANT;
48 double rangeSize = end-start;
49 int offset = rangeSize*normedRVal;
50 return start+offset;
51 }
52
randomBlock(Block & block,Field & field,int q,int m,int n)53 void randomBlock(Block& block, Field& field, int q, int m, int n)
54 {
55 Element d;
56 for (int i=0;i<m;++i) {
57 for (int j=0;j<n;++j) {
58 field.init(d,randRange(0,q));
59 block.setEntry(i,j,d);
60 }
61 }
62 }
63
randomNonSingular(Block & block,Field & field,int q,int m,int n)64 void randomNonSingular(Block& block, Field& field, int q, int m, int n)
65 {
66 randomBlock(block,field,q,m,n);
67 long unsigned r;
68 long unsigned b=(long unsigned)(m<n?m:n);
69 std::cerr << "starting" << std::endl;
70 while (LinBox::rank(r,block,Method::DenseElimination()) < b) {
71 randomBlock(block,field,q,m,n);
72 std::cerr << "loop" << std::endl;
73 }
74 }
75
fileDesc(int n,int b,double s,int q)76 std::string fileDesc(int n, int b, double s, int q)
77 {
78 std::stringstream ss;
79 ss << "n" << n;
80 ss << "b" << b;
81 int sInt=10000*(s-floor(s));
82 ss << "s" << setfill('0') << setw(4) << sInt << setw(0);
83 ss << "q" << q;
84 ss << ".txt";
85
86 return ss.str();
87 }
88
computeMinPoly(Domain & MD,SparseMat & M,Field & F,Block & U,Block & V,std::vector<Block> & gen,std::vector<size_t> & deg,int t)89 void computeMinPoly(Domain& MD,
90 SparseMat& M,
91 Field& F,
92 Block& U,
93 Block& V,
94 std::vector<Block>& gen,
95 std::vector<size_t>& deg,
96 int t)
97 {
98 BlackboxBlockContainer<Field,SparseMat> blockseq(&M,F,U,V);
99 BlockCoppersmithDomain<Domain,BlackboxBlockContainer<Field,SparseMat> > BCD(MD,&blockseq,t);
100
101 deg=BCD.right_minpoly(gen);
102 }
103
main(int argc,char ** argv)104 int main(int argc, char** argv) {
105 int seed=501;
106 int n;
107 int b;
108 int nnz;
109 double sparsity;
110 int t;
111 int q=3;
112
113 srand(seed);
114
115 static Argument args[] = {
116 { 'n', "-n N", "Set row/col dimension of test matrix n.", TYPE_INT, &n },
117 { 'b', "-b B", "Set block size b.", TYPE_INT, &b },
118 { 's', "-s S", "Set the sparsity [0-1).", TYPE_DOUBLE, &sparsity },
119 { 'q', "-q Q", "Set the field GF(p)", TYPE_INT, &q},
120 { 't', "-t T", "Early term threshold", TYPE_INT, &t},
121 END_OF_ARGUMENTS
122 };
123
124 parseArguments (argc, argv, args);
125
126 nnz=(int)((((double)(n))*n)*sparsity);
127 std::cout << nnz << std::endl;
128
129 Field F(q);
130 Domain MD(F);
131 SparseMat M(F,n,n);
132 Block U(F,b,n),V(F,n,b);
133
134 std::vector<Block> gen;
135 std::vector<size_t> deg;
136
137 randomNonSingular(U,F,q,b,n);
138 randomNonSingular(V,F,q,n,b);
139
140 std::cerr << "done1" << std::endl;
141 MapSparse<Field> sparse(F,n,n);
142 MapSparse<Field>::generateSparseNonSingular(sparse,nnz,seed);
143 sparse.copy(M);
144 std::cerr << "done2" << std::endl;
145
146 {
147 ofstream oF("U"+fileDesc(n,b,sparsity,q));
148 U.write(oF);
149 oF.close();
150 }
151
152 {
153 ofstream oF("V"+fileDesc(n,b,sparsity,q));
154 V.write(oF);
155 oF.close();
156 }
157
158 {
159 ofstream oF("M"+fileDesc(n,b,sparsity,q));
160 M.write(oF);
161 oF.close();
162 }
163
164 std::cerr << "done3" << std::endl;
165 computeMinPoly(MD,M,F,U,V,gen,deg,t);
166
167 {
168 ofstream oF("MP"+fileDesc(n,b,sparsity,q));
169 for (int i=0;i<deg.size()-1;++i) {
170 oF << deg[i] << " ";
171 }
172 if (deg.size() > 0) {
173 oF << deg[deg.size()-1];
174 }
175 oF << std::endl;
176
177 for (int i=0;i<gen.size();++i) {
178 oF << gen[i] << std::endl;
179 }
180 }
181
182 std::cout << fileDesc(n,b,sparsity,q) << std::endl;
183 return 0;
184 }
185
186 // Local Variables:
187 // mode: C++
188 // tab-width: 4
189 // indent-tabs-mode: nil
190 // c-basic-offset: 4
191 // End:
192 // vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
193