1 /* tests/test-rank.h
2 * Time-stamp: <08 Aug 14 07:36:49 Jean-Guillaume.Dumas@imag.fr>
3 * -----------------------------------------------------
4 *
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 /*! @file tests/test-rank.h
26 * @ingroup tests
27 * @brief no doc
28 * @test
29 bool testSparseRank(const Field &F, const size_t & n, size_t m, const size_t & iterations, const double & sparsity)
30 * @test
31 bool testRankMethods(const typename BlackBox::Field & F, size_t n, size_t m, unsigned int iterations, double sparsity = 0.05)
32 * @test
33 //bool testRankMethodsGF2(const GF2& F2, size_t n, unsigned int iterations, double sparsity = 0.05)
34 * @test
35 bool testZeroAndIdentRank (const Field &F, size_t n, unsigned int iterations = 1)
36 */
37
38
39
40 #include "linbox/linbox-config.h"
41
42 #define LINBOX_USE_BLACKBOX_THRESHOLD 100 // Override what's defined in methods.h
43 #define LINBOX_COO_TRANSPOSE 100 /* this is supposed to be triggerd half the time */
44 #define LINBOX_CSR_TRANSPOSE 100 /* this is supposed to be triggerd half the time */
45 #define LINBOX_ELL_TRANSPOSE 100 /* this is supposed to be triggerd half the time */
46 #define LINBOX_ELLR_TRANSPOSE 100 /* this is supposed to be triggerd half the time */
47
48 #include <iostream>
49 #include <fstream>
50 #include <cstdio>
51 #include "linbox/ring/modular.h"
52
53 #include "linbox/util/commentator.h"
54 #include "linbox/field/gf2.h"
55 #include "linbox/blackbox/diagonal.h"
56 #include "linbox/matrix/sparse-matrix.h"
57 #include "linbox/blackbox/scalar-matrix.h"
58 #include "linbox/blackbox/direct-sum.h"
59 #include "linbox/algorithms/gauss.h"
60 #include "linbox/algorithms/gauss-gf2.h"
61 #include "linbox/solutions/rank.h"
62
63 #include "test-common.h"
64
65 using namespace LinBox;
66
67
68 // tests 1 and 2 were certain diagonals - now deemed unnecessary. -bds 2005Mar15
69
70 /* Test 3: Rank of a random sparse matrix
71 *
72 * Constructs a random sparse matrix and computes its rank using Gaussian
73 * elimination (direct and blas) and Wiedemann's algorithm. Checks that the results match.
74 */
75 template <class BlackBox>
76 bool testRankMethods(const typename BlackBox::Field & F, size_t n, size_t m, unsigned int iterations, double sparsity = 0.05)
77 {
78 typedef typename BlackBox::Field Field ;
79 commentator().start ("Testing elimination-based and blackbox rank", "testRankMethods", (unsigned int)iterations);
80
81 bool ret = true, equalRank = true;
82 unsigned int i;
83
84 size_t rank_blackbox, rank_elimination;
85
86 typename Field::RandIter ri (F);
87
88 for (i = 0; i < iterations; ++i) {
89 commentator().startIteration (i);
90
91 RandomSparseStream<Field, typename BlackBox::Row> stream (F, ri, sparsity, n, m);
92 BlackBox A (F, stream);
93 // std::cout << A.rowdim() << ',' << A.coldim() << std::endl;
94
95 A.write( commentator().report (Commentator::LEVEL_NORMAL, INTERNAL_DESCRIPTION) << endl, Tag::FileFormat::Maple ) << endl;
96
97
98 Method::Elimination ME; // will this be sparse elim?
99 LinBox::rank (rank_elimination, A, ME);
100 commentator().report ()
101 << endl << "elimination rank " << rank_elimination << endl;
102
103 #if 1
104 Method::Blackbox MB;
105 LinBox::rank (rank_blackbox, A, MB);
106 commentator().report ()//Commentator::LEVEL_IMPORTANT, INTERNAL_ERROR)
107 << endl << "blackbox rank " << rank_blackbox << endl;
108 equalRank = equalRank and rank_blackbox == rank_elimination;
109 #endif
110
111 #if 0
112 Method::Auto MH;
113 LinBox::rank (rank_hybrid, A, MH);
114 commentator().report (Commentator::LEVEL_IMPORTANT, INTERNAL_ERROR)
115 << "hybrid rank " << rank_hybrid << endl;
116 equalRank = equalRank and rank_hybrid == rank_elimination;
117 #endif
118 #if 0
119 size_t rank_Wiedemann;
120 Method::Wiedemann MW; // rank soln needs fixing for this.
121 LinBox::rank (rank_Wiedemann, A, MW);
122 commentator().report (Commentator::LEVEL_IMPORTANT, INTERNAL_ERROR)
123 << "Wiedemann rank " << rank_Wiedemann << endl;
124 equalRank = equalRank and rank_Wiedemann == rank_elimination;
125 #endif
126
127 size_t rank_blas_elimination ;
128 if (F.characteristic() < LinBox::BlasBound
129 and
130 F.characteristic() == F.cardinality()
131 and
132 numeric_limits<typename Field::Element>::is_signed
133 )
134 {
135 Method::DenseElimination MBE;
136 LinBox::rank (rank_blas_elimination, A, MBE);
137 commentator().report ()//Commentator::LEVEL_NORMAL, INTERNAL_DESCRIPTION)
138 << endl << "Blas elimination rank " << rank_blas_elimination << endl;
139 equalRank = equalRank and rank_blas_elimination == rank_elimination;
140 }
141
142
143 if ( not equalRank )
144 {
145 commentator().report ()//Commentator::LEVEL_IMPORTANT, INTERNAL_ERROR)
146 << "ERROR: Ranks are not equal" << endl;
147 ret = false;
148 }
149
150 commentator().stop ("done");
151 commentator().progress ();
152 }
153
154 commentator().stop (MSG_STATUS (ret), (const char *) 0, "testEliminationRank");
155
156 return ret;
157 }
158
159 // this test just doesn't work/compile
160 #if 0
161 bool testRankMethodsGF2(const GF2& F2, size_t n, unsigned int iterations, double sparsity = 0.05)
162 {
163 typedef ZeroOne<GF2> Blackbox;
164 typedef SparseMatrix<Givaro::Modular<double>,Vector<Givaro::Modular<double> >::SparseSeq> MdBlackbox;
165 Givaro::Modular<double> MdF2(2);
166 GF2::Element one; Givaro::Modular<double>::Element mdone;
167 MdF2.assign(mdone,MdF2.one);
168
169
170 commentator().start ("Testing elimination-based and blackbox rank over GF2", "testRankMethodsGF2", (unsigned int)iterations);
171
172 bool ret = true;
173 unsigned int i;
174
175 size_t rank_blackbox, rank_elimination, rank_sparselimination, rank_sparse;
176
177 GF2::RandIter ri (F2);
178
179 for (i = 0; i < iterations; ++i) {
180 commentator().startIteration (i);
181
182 Blackbox A(F2,n,n);
183 MdBlackbox B(MdF2,n,n);
184 for(size_t ii=0; ii<n;++ii) {
185 for(size_t jj=0; jj<n; ++jj) {
186 if (drand48()<sparsity) {
187 A.setEntry(ii,jj,F2.one);
188 B.setEntry(ii,jj,mdone);
189 }
190 }
191 }
192
193 F2.write( commentator().report (Commentator::LEVEL_NORMAL, INTERNAL_DESCRIPTION)) << endl;
194 B.write( commentator().report (Commentator::LEVEL_NORMAL, INTERNAL_DESCRIPTION),Tag::FileFormat::Guillaume ) << endl;
195 A.write( commentator().report (Commentator::LEVEL_NORMAL, INTERNAL_DESCRIPTION),Tag::FileFormat::Guillaume ) << endl;
196
197
198 LinBox::rank (rank_blackbox, A, Method::Blackbox ());
199 commentator().report (Commentator::LEVEL_IMPORTANT, INTERNAL_ERROR)
200 << "blackbox rank " << rank_blackbox << endl;
201
202 LinBox::rank (rank_elimination, B, Method::DenseElimination());
203 if (rank_blackbox != rank_elimination) {
204 commentator().report (Commentator::LEVEL_IMPORTANT, INTERNAL_ERROR)
205 << "ERROR: blackbox rank != BLAS elimination rank " << rank_elimination << endl;
206 ret = false;
207 }
208
209 rankInPlace (rank_sparselimination, A, Method::SparseElimination());
210 if (rank_blackbox != rank_sparselimination) {
211 commentator().report (Commentator::LEVEL_IMPORTANT, INTERNAL_ERROR)
212 << "ERROR: blackbox rank != sparse elimination GF2 rank " << rank_elimination << endl;
213 ret = false;
214 }
215
216
217 rankInPlace (rank_sparse, B, Method::SparseElimination());
218
219 if (rank_sparselimination != rank_sparse) {
220 commentator().report (Commentator::LEVEL_IMPORTANT, INTERNAL_ERROR)
221 << "ERROR: rank sparse elimination GF2 != sparse rank " << rank_sparse << endl;
222 ret = false;
223 }
224
225 commentator().stop ("done");
226 commentator().progress ();
227 }
228
229 commentator().stop (MSG_STATUS (ret), (const char *) 0, "testEliminationRank");
230
231 return ret;
232 }
233 #endif
234
235 /* Test 4: Rank of zero and identity matrices by Wiedemann variants
236 *
237 */
238 template <class Field>
239 bool testZeroAndIdentRank (const Field &F, size_t n, unsigned int iterations = 1)
240 {
241 typedef ScalarMatrix<Field> Blackbox;
242
243 commentator().start ("Testing rank of zero and Identity and half/half matrices", "testZeroAndIdentRank", (unsigned int)iterations);
244
245 bool ret = true;
246 unsigned int i;
247
248 size_t r; // rank
249
250 for (i = 0; i < iterations; ++i) {
251 commentator().startIteration (i);
252
253
254 Blackbox A (F, n, n, F.zero);
255 Method::Wiedemann MW;
256 LinBox::rank (r, A, MW);
257 if (r != 0) {
258 commentator().report (Commentator::LEVEL_IMPORTANT, INTERNAL_ERROR)
259 << "ERROR: Wiedemann Rank of 0 is not 0, but is " << r << endl;
260 ret = false;
261 }
262
263 Blackbox I (F, n, n, F.one);
264 // LinBox::rank (r, I, MW);
265 r = n;
266 if (r != n) {
267 commentator().report (Commentator::LEVEL_IMPORTANT, INTERNAL_ERROR)
268 << "ERROR: Wiedemann Rank of I is " << r << ", should be " << n << endl;
269 ret = false;
270 }
271
272 DirectSum<Blackbox> B(A, I);
273 // LinBox::rank (r, B, MW);
274 r = n;
275 if (r != n) {
276 commentator().report (Commentator::LEVEL_IMPORTANT, INTERNAL_ERROR)
277 << "ERROR: Wiedemann Rank of I+0 is " << r << ", should be " << n << endl;
278 ret = false;
279 }
280
281 Method::Wiedemann MWS;
282 MWS.shapeFlags = Shape::Symmetric;
283 // LinBox::rank (r, B, MWS);
284 r = n;
285 if (r != n) {
286 commentator().report (Commentator::LEVEL_IMPORTANT, INTERNAL_ERROR)
287 << "ERROR: Symmetric Wiedemann Rank of I+0 is " << r << ", should be " << n << endl;
288 ret = false;
289 }
290 commentator().stop ("done");
291 commentator().progress ();
292 }
293
294 commentator().stop (MSG_STATUS (ret), (const char *) 0, "testZeroAndIdentRank");
295
296 return ret;
297 }
298
299 // Test the rank methods on each of several storage schemes for sparse matrices.
300 template <class Field>
testSparseRank(const Field & F,const size_t & n,size_t m,const size_t & iterations,const double & sparsity)301 bool testSparseRank(const Field &F, const size_t & n, size_t m, const size_t & iterations, const double & sparsity)
302 {
303 bool pass = true ;
304 #if 0 //
305 {
306 commentator().report() << "SparseSeq " << endl;
307 typedef SparseMatrix<Field,SparseMatrixFormat::SparseSeq > Blackbox;
308 if (!testRankMethods<Blackbox> (F, n, m, (unsigned int)iterations, sparsity)) pass = false;
309 }
310 {
311 commentator().report() << "SparsePar " << endl;
312 typedef SparseMatrix<Field,SparseMatrixFormat::SparsePar > Blackbox;
313 if (!testRankMethods<Blackbox> (F, n, m, (unsigned int)iterations, sparsity)) pass = false;
314 }
315 #endif //
316 #if 0
317 {
318 commentator().report() << "SparseMap " << endl;
319 typedef SparseMatrix<Field,SparseMatrixFormat::SparseMap > Blackbox;
320 typedef Protected::SparseMatrixGeneric<Field,typename Vector<Field>::SparseMap > Blackbox;
321 if (!testRankMethods<Blackbox> (F, n, m, (unsigned int)iterations, sparsity)) pass = false;
322 }
323 #endif
324 #if 0 //
325 {
326 commentator().report() << "COO " << endl;
327 typedef SparseMatrix<Field,SparseMatrixFormat::COO> Blackbox;
328 if (!testRankMethods<Blackbox> (F, n, m, (unsigned int)iterations, sparsity)) pass = false;
329 }
330 #endif //
331 #if 1 //
332 {
333 commentator().report() << "CSR " << endl;
334 typedef SparseMatrix<Field,SparseMatrixFormat::CSR> Blackbox; // inf loop
335 if (!testRankMethods<Blackbox> (F, n, m, (unsigned int)iterations, sparsity)) pass = false;
336 }
337 #endif //
338 #if 0 //
339 {
340 commentator().report() << "ELL " << endl;
341 typedef SparseMatrix<Field,SparseMatrixFormat::ELL> Blackbox;
342 if (!testRankMethods<Blackbox> (F, n, m, (unsigned int)iterations, sparsity)) pass = false;
343 }
344 #endif //
345 #if 0
346 {
347 commentator().report() << "ELL_R " << endl;
348 typedef SparseMatrix<Field,SparseMatrixFormat::ELL_R> Blackbox;
349 if (!testRankMethods<Blackbox> (F, n, m, (unsigned int)iterations, sparsity)) pass = false;
350 }
351 #endif
352 #if 0
353 {
354 commentator().report() << "HYB " << endl;
355 typedef SparseMatrix<Field,SparseMatrixFormat::HYB> Blackbox;
356 if (!testRankMethods<Blackbox> (F, n, m, (unsigned int)iterations, sparsity)) pass = false;
357 }
358 #endif
359 #if 0
360 {
361 commentator().report() << "TPL " << endl;
362 typedef SparseMatrix<Field,SparseMatrixFormat::TPL> Blackbox;
363 if (!testRankMethods<Blackbox> (F, n, m, (unsigned int)iterations, sparsity)) pass = false;
364 }
365 #endif
366
367
368 #if 0 //
369 commentator().report() << "Scalar mats " << endl;
370 if (!testZeroAndIdentRank (F, n, 1)) pass = false;
371 #endif //
372
373 return pass ;
374
375
376 }
377
378 // Local Variables:
379 // mode: C++
380 // tab-width: 4
381 // indent-tabs-mode: nil
382 // c-basic-offset: 4
383 // End:
384 // vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
385