1 /*
2 * Copyright (C) FFLAS-FFPACK
3 * Written by Clément Pernet
4 * This file is Free Software and part of FFLAS-FFPACK.
5 *
6 * ========LICENCE========
7 * This file is part of the library FFLAS-FFPACK.
8 *
9 * FFLAS-FFPACK 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 //--------------------------------------------------------------------------
27 // Test for the computations of rank profiles
28 //--------------------------------------------------------------------------
29 #define __FFLASFFPACK_SEQUENTIAL
30 #include "fflas-ffpack/fflas-ffpack-config.h"
31 #include "fflas-ffpack/ffpack/ffpack.h"
32 #include "fflas-ffpack/utils/args-parser.h"
33 #include "fflas-ffpack/utils/test-utils.h"
34 #include <givaro/modular.h>
35
36 #include <iostream>
37 #include <iomanip>
38 #include <random>
39 #include <chrono>
40
41 using namespace FFPACK;
42 using namespace FFLAS;
43
44
45 template<class Field>
run_with_field(Givaro::Integer q,uint64_t b,size_t m,size_t n,size_t r,size_t iters,uint64_t seed,bool par)46 bool run_with_field(Givaro::Integer q, uint64_t b, size_t m, size_t n, size_t r, size_t iters, uint64_t seed, bool par){
47 bool ok = true ;
48 int nbit=(int)iters;
49
50 while (ok && nbit){
51 // choose Field
52 Field* F= chooseField<Field>(q,b,seed);
53 if (F==nullptr)
54 return true;
55
56 typename Field::RandIter G (*F,b,seed++);
57 std::ostringstream oss;
58 F->write(oss);
59
60 std::cout.fill('.');
61 std::cout<<"Checking ";
62 std::cout.width(40);
63 std::cout<<oss.str();
64 std::cout<<" ... ";
65
66
67 size_t lda = n;
68 typename Field::Element_ptr A=fflas_new (*F, m,lda);
69 typename Field::Element_ptr B=fflas_new (*F, m,lda);
70 RandomMatrixWithRankandRandomRPM(*F,m,n,r,A,lda,G);
71 fassign (*F, m, n, A, lda, B, lda);
72
73 {
74 // Testing if LUdivine and PLUQ return the same result
75 size_t* RP1, * RP2;
76 if(!par)
77 FFPACK::RowRankProfile (*F, m, n, A, lda, RP1, FFPACK::FfpackSlabRecursive);
78 else
79 FFPACK::pRowRankProfile (*F, m, n, A, lda, RP1, 0, FFPACK::FfpackSlabRecursive);
80 fassign (*F, m, n, B, lda, A, lda);
81 if(!par)
82 FFPACK::RowRankProfile (*F, m, n, A, lda, RP2, FFPACK::FfpackTileRecursive);
83 else
84 FFPACK::pRowRankProfile (*F, m, n, A, lda, RP2, 0, FFPACK::FfpackTileRecursive);
85 for (size_t i=0; i<r; i++)
86 ok = ok && (RP1[i] == RP2[i]);
87 fflas_delete(RP1);
88 fflas_delete(RP2);
89
90 fassign (*F, m, n, B, lda, A, lda);
91 if(!par)
92 FFPACK::ColumnRankProfile (*F, m, n, A, lda, RP1, FFPACK::FfpackSlabRecursive);
93 else
94 FFPACK::pColumnRankProfile (*F, m, n, A, lda, RP1, 0, FFPACK::FfpackSlabRecursive);
95 fassign (*F, m, n, B, lda, A, lda);
96 if(!par)
97 FFPACK::ColumnRankProfile (*F, m, n, A, lda, RP2, FFPACK::FfpackTileRecursive);
98 else
99 FFPACK::pColumnRankProfile (*F, m, n, A, lda, RP2, 0, FFPACK::FfpackTileRecursive);
100
101 for (size_t i=0; i<r; i++)
102 ok = ok && (RP1[i] == RP2[i]);
103 fflas_delete(RP1);
104 fflas_delete(RP2);
105 }
106 {
107 // Testing if 1 PLUQ computes the rank profiles of all leading submatrices
108 size_t* RP1, * RP2;
109 size_t * P = fflas_new<size_t>(m);
110 size_t * Q = fflas_new<size_t>(n);
111 fassign (*F, m, n, B, lda, A, lda);
112 PLUQ(*F, FflasNonUnit, m, n, A, lda, P, Q);
113
114 for (size_t i=0; i<3;i++){
115 size_t mm = 1 + (rand() % m);
116 size_t nn = 1 + (rand() % n);
117 fassign (*F, m, n, B, lda, A, lda);
118 size_t rr;
119 if(!par)
120 rr = FFPACK::ColumnRankProfile (*F, mm, nn, A, lda, RP1, FFPACK::FfpackSlabRecursive);
121 else
122 rr = FFPACK::pColumnRankProfile (*F, mm, nn, A, lda, RP1, 0, FFPACK::FfpackSlabRecursive);
123 fassign (*F, m, n, B, lda, A, lda);
124 if(!par)
125 FFPACK::RowRankProfile (*F, mm, nn, A, lda, RP2, FFPACK::FfpackSlabRecursive);
126 else
127 FFPACK::pRowRankProfile (*F, mm, nn, A, lda, RP2, 0, FFPACK::FfpackSlabRecursive);
128 size_t* RRP = fflas_new<size_t>(r);
129 size_t* CRP = fflas_new<size_t>(r);
130
131 LeadingSubmatrixRankProfiles (m,n,r,mm,nn,P,Q,RRP,CRP);
132 for (size_t ii=0; ii<rr; ii++)
133 ok = ok && (RP1[ii] == CRP[ii]) && (RP2[ii] == RRP[ii]);
134
135 fflas_delete(CRP);
136 fflas_delete(RRP);
137 fflas_delete(RP1);
138 fflas_delete(RP2);
139
140 }
141 fflas_delete(P);
142 fflas_delete(Q);
143 }
144 {
145 // Testing PLUQ and LUDivine return a specified rank profile
146 size_t* RRP = fflas_new<size_t>(r);
147 size_t* CRP = fflas_new<size_t>(r);
148 size_t* RRPLUD, * RRPPLUQ, *CRPLUD, *CRPPLUQ;
149
150 RandomRankProfileMatrix (m, n, r, RRP, CRP);
151 RandomMatrixWithRankandRPM(*F,m,n,r,A,lda, RRP, CRP, G);
152 size_t cs, ct;
153 fassign (*F, m, n, A, lda, B, lda);
154 if(!par)
155 cs = FFPACK::ColumnRankProfile (*F, m, n, A, lda, CRPLUD, FFPACK::FfpackSlabRecursive);
156 else
157 cs = FFPACK::pColumnRankProfile (*F, m, n, A, lda, CRPLUD, 0, FFPACK::FfpackSlabRecursive);
158 fassign (*F, m, n, B, lda, A, lda);
159 if(!par)
160 ct = FFPACK::ColumnRankProfile (*F, m, n, A, lda, CRPPLUQ, FFPACK::FfpackTileRecursive);
161 else
162 ct = FFPACK::pColumnRankProfile (*F, m, n, A, lda, CRPPLUQ, 0, FFPACK::FfpackTileRecursive);
163 fassign (*F, m, n, B, lda, A, lda);
164 size_t rs, rt;
165 if(!par)
166 rs = FFPACK::RowRankProfile (*F, m, n, A, lda, RRPLUD, FFPACK::FfpackSlabRecursive);
167 else
168 rs = FFPACK::pRowRankProfile (*F, m, n, A, lda, RRPLUD, 0, FFPACK::FfpackSlabRecursive);
169 fassign (*F, m, n, B, lda, A, lda);
170 if(!par)
171 rt = FFPACK::RowRankProfile (*F, m, n, A, lda, RRPPLUQ, FFPACK::FfpackTileRecursive);
172 else
173 rt = FFPACK::pRowRankProfile (*F, m, n, A, lda, RRPPLUQ, 0, FFPACK::FfpackTileRecursive);
174 std::sort(CRP,CRP+r);
175 std::sort(RRP,RRP+r);
176 ok = ok && (cs==ct)&&(cs==rs)&&(cs==rt)&&(cs==r);
177 for (size_t i=0; i<r; i++)
178 ok = ok && (CRPLUD[i] == CRP[i]) && (CRPPLUQ[i] == CRP[i]) &&
179 (RRPLUD[i] == RRP[i]) && (RRPPLUQ[i] == RRP[i]);
180 fflas_delete(CRP);
181 fflas_delete(RRP);
182 fflas_delete(CRPLUD);
183 fflas_delete(RRPLUD);
184 fflas_delete(CRPPLUQ);
185 fflas_delete(RRPPLUQ);
186 }
187 fflas_delete(A);
188 fflas_delete(B);
189 {
190 // Testing PLUQ and LUDivine return a specified rank profile (symmetric case)
191 size_t* RRP = fflas_new<size_t>(r);
192 size_t* CRP = fflas_new<size_t>(r);
193 size_t* RRPLUD, * RRPPLUQ, *CRPLUD, *CRPPLUQ;
194
195 typename Field::Element_ptr A=fflas_new (*F, n,lda);
196 typename Field::Element_ptr B=fflas_new (*F, n,lda);
197 RandomSymmetricRankProfileMatrix (n, r, RRP, CRP);
198 RandomSymmetricMatrixWithRankandRPM(*F,n,r,A,lda, RRP, CRP, G);
199 // typename Field::Element_ptr RPM=FFLAS::fflas_new(*F,n*n);
200 // fzero(*F,n,n,RPM,n);
201 // for (size_t i=0; i<r; i++)
202 // F->assign(RPM[RRP[i]*n+CRP[i]],F->one);
203 size_t cs, ct;
204 fassign (*F, n, n, A, lda, B, lda);
205 if(!par)
206 cs= FFPACK::ColumnRankProfile (*F, n, n, A, lda, CRPLUD, FFPACK::FfpackSlabRecursive);
207 else
208 cs = FFPACK::pColumnRankProfile (*F, n, n, A, lda, CRPLUD, 0, FFPACK::FfpackSlabRecursive);
209 fassign (*F, n, n, B, lda, A, lda);
210 if(!par)
211 ct = FFPACK::ColumnRankProfile (*F, n, n, A, lda, CRPPLUQ, FFPACK::FfpackTileRecursive);
212 else
213 ct = FFPACK::pColumnRankProfile (*F, n, n, A, lda, CRPPLUQ, 0, FFPACK::FfpackTileRecursive);
214 size_t rs, rt;
215 fassign (*F, n, n, B, lda, A, lda);
216 if(!par)
217 rs = FFPACK::RowRankProfile (*F, n, n, A, lda, RRPLUD, FFPACK::FfpackSlabRecursive);
218 else
219 rs = FFPACK::pRowRankProfile (*F, n, n, A, lda, RRPLUD, 0, FFPACK::FfpackSlabRecursive);
220 fassign (*F, n, n, B, lda, A, lda);
221 if(!par)
222 rt = FFPACK::RowRankProfile (*F, n, n, A, lda, RRPPLUQ, FFPACK::FfpackTileRecursive);
223 else
224 rt = FFPACK::pRowRankProfile (*F, n, n, A, lda, RRPPLUQ, 0, FFPACK::FfpackTileRecursive);
225 std::sort(CRP,CRP+r);
226 std::sort(RRP,RRP+r);
227 ok = ok && (cs==ct) && (cs==rs) && (cs==rt) && (cs==r);
228 for (size_t i=0; i<r; i++)
229 ok = ok && (CRPLUD[i] == CRP[i]) && (CRPPLUQ[i] == CRP[i]) &&
230 (RRPLUD[i] == RRP[i]) && (RRPPLUQ[i] == RRP[i]);
231 fflas_delete(CRP);
232 fflas_delete(RRP);
233 fflas_delete(CRPLUD);
234 fflas_delete(RRPLUD);
235 fflas_delete(CRPPLUQ);
236 fflas_delete(RRPPLUQ);
237 fflas_delete(A);
238 fflas_delete(B);
239 }
240 delete F;
241
242 nbit--;
243 if (!ok)
244 //std::cout << "\033[1;31mFAILED\033[0m "<<std::endl;
245 std::cout << "FAILED "<<std::endl;
246 else
247 //std::cout << "\033[1;32mPASSED\033[0m "<<std::endl;
248 std::cout << "PASSED "<<std::endl;
249 }
250 return ok;
251 }
252
main(int argc,char ** argv)253 int main(int argc, char** argv){
254 std::cerr<<std::setprecision(20);
255
256 Givaro::Integer q = -1;
257 size_t b = 0;
258 size_t m = 150;
259 size_t n = 280;
260 size_t r = 85;
261 size_t iters = 6 ;
262 bool loop=false;
263 uint64_t seed = getSeed();
264 Argument as[] = {
265 { 'q', "-q Q", "Set the field cardinality.", TYPE_INTEGER , &q },
266 { 'b', "-b B", "Set the bitsize of the field characteristic.", TYPE_INT , &b },
267 { 'n', "-n N", "Set the number of cols in the matrix.", TYPE_INT , &n },
268 { 'm', "-m N", "Set the number of rows in the matrix.", TYPE_INT , &m },
269 { 'r', "-r r", "Set the rank of the matrix." , TYPE_INT , &r },
270 { 'i', "-i R", "Set number of repetitions.", TYPE_INT , &iters },
271 { 'l', "-loop Y/N", "run the test in an infinite loop.", TYPE_BOOL , &loop },
272 { 's', "-s seed", "Set seed for the random generator", TYPE_UINT64, &seed },
273 // { 'f', "-f file", "Set input file", TYPE_STR, &file },
274 END_OF_ARGUMENTS
275 };
276
277 parseArguments(argc,argv,as);
278
279 srand(seed);
280
281 if (r > std::min (m,n))
282 r = std::min (m, n);
283
284 bool ok=true;
285 do{
286 ok = ok &&run_with_field<Givaro::Modular<float> > (q,b,m,n,r,iters,seed,false);
287 ok = ok &&run_with_field<Givaro::Modular<double> > (q,b,m,n,r,iters,seed,false);
288 ok = ok &&run_with_field<Givaro::ModularBalanced<float> > (q,b,m,n,r,iters,seed,false);
289 ok = ok &&run_with_field<Givaro::ModularBalanced<double> > (q,b,m,n,r,iters,seed,false);
290 ok = ok &&run_with_field<Givaro::Modular<int32_t> > (q,b,m,n,r,iters,seed,false);
291 ok = ok &&run_with_field<Givaro::ModularBalanced<int32_t> > (q,b,m,n,r,iters,seed,false);
292 ok = ok &&run_with_field<Givaro::Modular<int64_t> > (q,b,m,n,r,iters,seed,false);
293 ok = ok &&run_with_field<Givaro::ModularBalanced<int64_t> > (q,b,m,n,r,iters,seed,false);
294 ok = ok &&run_with_field<Givaro::Modular<Givaro::Integer> >(q,(b?b:128),m/4+1,n/4+1,r/4+1,iters,seed,false);
295
296 ok = ok &&run_with_field<Givaro::Modular<float> > (q,b,m,n,r,iters,seed,true);
297 ok = ok &&run_with_field<Givaro::Modular<double> > (q,b,m,n,r,iters,seed,true);
298 ok = ok &&run_with_field<Givaro::ModularBalanced<float> > (q,b,m,n,r,iters,seed,true);
299 ok = ok &&run_with_field<Givaro::ModularBalanced<double> > (q,b,m,n,r,iters,seed,true);
300 ok = ok &&run_with_field<Givaro::Modular<int32_t> > (q,b,m,n,r,iters,seed,true);
301 ok = ok &&run_with_field<Givaro::ModularBalanced<int32_t> > (q,b,m,n,r,iters,seed,true);
302 ok = ok &&run_with_field<Givaro::Modular<int64_t> > (q,b,m,n,r,iters,seed,true);
303 ok = ok &&run_with_field<Givaro::ModularBalanced<int64_t> > (q,b,m,n,r,iters,seed,true);
304 ok = ok &&run_with_field<Givaro::Modular<Givaro::Integer> >(q,(b?b:128),m/4+1,n/4+1,r/4+1,iters,seed,true);
305 } while (loop && ok);
306
307 return !ok;
308 }
309 /* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
310 // vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
311