1 /*
2 * Copyright (C) the FFLAS-FFPACK group
3 * Written by Clément Pernet
4 * Brice Boyer (briceboyer) <boyer.brice@gmail.com>
5 * This file is Free Software and part of FFLAS-FFPACK.
6 *
7 * ========LICENCE========
8 * This file is part of the library FFLAS-FFPACK.
9 *
10 * FFLAS-FFPACK is free software: you can redistribute it and/or modify
11 * it under the terms of the GNU Lesser General Public
12 * License as published by the Free Software Foundation; either
13 * version 2.1 of the License, or (at your option) any later version.
14 *
15 * This library is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 * Lesser General Public License for more details.
19 *
20 * You should have received a copy of the GNU Lesser General Public
21 * License along with this library; if not, write to the Free Software
22 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23 * ========LICENCE========
24 *.
25 */
26
27 // #ifndef NEWINO
28 // #define NEWWINO
29 // #endif
30
31 // #define WINOTHRESHOLD 100
32 // #define OLD_DYNAMIC_PEELING
33
34
35
36 #include "fflas-ffpack/fflas-ffpack-config.h"
37
38 #include <iomanip>
39 #include <iostream>
40
41 #include <givaro/modular.h>
42
43 #include <recint/rint.h>
44
45 #include "fflas-ffpack/utils/timer.h"
46 #include "fflas-ffpack/fflas/fflas.h"
47
48 #include "fflas-ffpack/utils/args-parser.h"
49 #include "fflas-ffpack/utils/test-utils.h"
50
51 using namespace std;
52 using namespace FFPACK;
53 using namespace FFLAS;
54
55 using Givaro::Modular;
56 using Givaro::ModularBalanced;
57
58
59 // checks that D = beta . Y + alpha . A ^ta * X
60 template<class Field>
check_MV(const Field & F,const typename Field::Element_ptr Cd,enum FFLAS_TRANSPOSE & ta,const size_t m,const size_t k,const typename Field::Element & alpha,const typename Field::Element_ptr A,size_t lda,const typename Field::Element_ptr X,size_t incX,const typename Field::Element & beta,const typename Field::Element_ptr Y,size_t incY)61 bool check_MV(const Field & F,
62 const typename Field::Element_ptr Cd, // c0
63 enum FFLAS_TRANSPOSE & ta,
64 const size_t m,
65 const size_t k,
66 const typename Field::Element & alpha,
67 const typename Field::Element_ptr A, size_t lda,
68 const typename Field::Element_ptr X, size_t incX,
69 const typename Field::Element & beta,
70 const typename Field::Element_ptr Y, size_t incY)
71 {
72 bool wrong = false;
73 typename Field::Element_ptr D;
74 if (ta == FflasNoTrans){
75 D = fflas_new(F,m);
76 fassign (F, m, Cd, 1, D, 1);
77 for (size_t i=0; i<m; i++){
78 F.mulin (D[i], beta);
79 typename Field::Element tmp;
80 F.init(tmp);
81 for (size_t j=0; j<k; j++){
82 F.axpyin (tmp, A[i*lda+j], X[j*incX]);
83 }
84 F.axpyin(D[i],alpha,tmp);
85 }
86 wrong = !fequal(F, m, D, 1, Y, incY);
87 } else {
88 D = fflas_new(F,k);
89 fassign (F, k, Cd, 1, D, 1);
90 for (size_t i=0; i<k; i++){
91 F.mulin (D[i], beta);
92 typename Field::Element tmp;
93 F.init(tmp);
94 for (size_t j=0; j<m; j++){
95 F.axpyin (tmp, A[i+j*lda], X[j*incX]);
96 }
97 F.axpyin(D[i],alpha,tmp);
98 }
99 wrong = !fequal(F, k, D, 1, Y, incY);
100 }
101 size_t Ydim = (ta==FflasNoTrans)? m : k;
102 if ( wrong ){
103 size_t canprint = 20 ;
104 std::cerr<<"FAIL"<<std::endl;
105 std::cerr << "alpha :" << alpha<<", beta : " << beta << std::endl;
106 std::cerr << "m :" << m << ", k : " << k << std::endl;
107 std::cerr << "ldA :" << lda << ", incX : " << incX << ", incY : " << incY << std::endl;
108 for (size_t i=0; i<Ydim && canprint; ++i){
109 if (!F.areEqual( Y[i*incY], D[i] ) ) {
110 std::cerr<<"Error Y["<<i<<"]="<<Y[i*incY]<<" D["<<i<<"]="<<D[i]<<std::endl;
111 canprint--;
112 }
113 }
114 if (Ydim<80) {
115 for (size_t i=0; i<Ydim ; ++i){
116 if (!F.areEqual( Y[i*incY], D[i] ) )
117 std::cout << 'X' ;
118 else
119 std::cout << '.' ;
120 }
121 std::cout << std::endl;
122 }
123 }
124
125 fflas_delete (D);
126
127 return !wrong ;
128
129 }
130
131
132 template<class Field, class RandIter>
launch_MV(const Field & F,const size_t m,const size_t k,const typename Field::Element alpha,const typename Field::Element beta,const size_t lda,enum FFLAS_TRANSPOSE ta,const size_t incX,const size_t incY,size_t iters,bool par,RandIter & G)133 bool launch_MV(const Field & F,
134 const size_t m,
135 const size_t k,
136 const typename Field::Element alpha,
137 const typename Field::Element beta,
138 const size_t lda,
139 enum FFLAS_TRANSPOSE ta,
140 const size_t incX,
141 const size_t incY,
142 size_t iters,
143 bool par,
144 RandIter& G)
145 {
146
147 bool ok = true;
148
149 typedef typename Field::Element_ptr Element_ptr;
150 Element_ptr A ;
151 for(size_t i = 0;i<iters;++i){
152 FFLASFFPACK_check(lda >= k);
153 A = fflas_new (F, m, lda);
154 fzero(F,m,lda,A,lda);
155 RandomMatrix(F, m, k, A, lda, G);
156 size_t Xdim = (ta == FflasNoTrans)? k : m;
157 size_t Ydim = (ta == FflasNoTrans)? m : k;
158 Element_ptr X = fflas_new (F, Xdim, incX);
159 Element_ptr Y = fflas_new (F, Ydim, incY);
160 fzero (F, Xdim, incX, X, incX);
161 fzero (F, Ydim, incY, Y, incY);
162 Element_ptr D = fflas_new (F, Ydim);
163
164 RandomMatrix (F, Xdim, 1, X, incX, G);
165 RandomMatrix (F, Ydim, 1, Y, incY, G);
166 fassign (F, Ydim, Y, incY, D, 1);
167
168 //Y will be modified so keep a copy of Y as Y2
169 Element_ptr Y2 = fflas_new (F, Ydim, incY);
170 fassign (F, Ydim, Y, incY, Y2, incY);
171
172
173 if (par){
174 {
175 MMHelper<Field, MMHelperAlgo::Classic, ModeTraits<Field>, ParSeqHelper::Parallel<CuttingStrategy::Recursive, StrategyParameter::Threads> > WH;
176
177 PAR_BLOCK{
178 fgemv(F, ta, m,k,alpha, A,lda, X, incX, beta, Y, incY, WH);
179 }
180 }
181 }else{
182 //MMHelper<Field,MMHelperAlgo::Auto,typename ModeTraits<Field>::value> WH(F,nbw,ParSeqHelper::Sequential());
183 fgemv(F, ta, m, k,alpha, A,lda, X, incX, beta, Y, incY);
184 }
185
186 ok = ok && check_MV(F, D, ta, m, k,alpha, A, lda, X, incX, beta, Y, incY);
187
188
189
190 if (!ok){
191 fflas_delete (A, X, Y, Y2, D);
192 break;
193 }
194
195
196
197 fassign (F, Ydim, Y2, incY, D, 1);
198 fassign (F, Ydim, Y2, incY, Y, incY);
199 if (par){
200 {
201 MMHelper<Field, MMHelperAlgo::Classic, ModeTraits<Field>, ParSeqHelper::Parallel<CuttingStrategy::Row, StrategyParameter::Threads> > WH;
202
203 PAR_BLOCK{
204 fgemv(F, ta, m,k,alpha, A,lda, X, incX, beta, Y, incY, WH);
205 }
206 }
207 }else{
208 //MMHelper<Field,MMHelperAlgo::Auto,typename ModeTraits<Field>::value> WH(F,nbw,ParSeqHelper::Sequential());
209 fgemv(F, ta, m, k,alpha, A,lda, X, incX, beta, Y, incY);
210 }
211
212 ok = ok && check_MV(F, D, ta, m, k,alpha, A, lda, X, incX, beta, Y, incY);
213
214
215 if (!ok){
216 fflas_delete (A, X, Y, Y2, D);
217 break;
218 }
219
220
221
222
223 fassign (F, Ydim, Y2, incY, D, 1);
224 fassign (F, Ydim, Y2, incY, Y, incY);
225 if (par){
226 {
227 ParSeqHelper::Parallel<CuttingStrategy::Row,StrategyParameter::Grain> WH(4);
228
229 PAR_BLOCK{
230 fgemv(F, ta, m,k,alpha, A,lda, X, incX, beta, Y, incY, WH);
231 }
232 }
233 }else{
234 //MMHelper<Field,MMHelperAlgo::Auto,typename ModeTraits<Field>::value> WH(F,nbw,ParSeqHelper::Sequential());
235 fgemv(F, ta, m, k,alpha, A,lda, X, incX, beta, Y, incY);
236 }
237
238 ok = ok && check_MV(F, D, ta, m, k,alpha, A, lda, X, incX, beta, Y, incY);
239
240 fflas_delete (A, X, Y, Y2, D);
241 if (!ok){
242
243 break;
244 }
245
246
247
248 }
249 return ok ;
250 }
251
252
253 template<class Field, class RandIter>
launch_MV_dispatch(const Field & F,const int mm,const int kk,const typename Field::Element alpha,const typename Field::Element beta,const size_t iters,const bool par,RandIter & G)254 bool launch_MV_dispatch(const Field &F,
255 const int mm,
256 const int kk,
257 const typename Field::Element alpha,
258 const typename Field::Element beta,
259 const size_t iters,
260 const bool par,
261 RandIter& G)
262 {
263 bool ok = true;
264 size_t m,k;
265 size_t lda,incX, incY;
266 size_t ld = 13 ;
267 {
268 //FFLAS_TRANSPOSE ta = FflasNoTrans ;
269 //if (! par) {
270 //if (random()%2) ta = FflasTrans ;
271 //}
272
273 if (mm<0)
274 m = 1+(size_t)random() % -mm;
275 else m = mm;
276 if (kk<0)
277 k = 1+(size_t)random() % -kk;
278 else k = kk;
279
280 lda = k+(size_t)random()%ld;
281 incX = 1+(size_t)random()%ld;
282 incY = 1+(size_t)random()%ld;
283
284 ok = ok && launch_MV (F, m, k, alpha,beta, lda, FflasNoTrans, incX, incY, iters, par, G);
285 ok = ok && launch_MV (F, m, k, alpha,beta, lda, FflasTrans, incX, incY, iters, par, G);
286 }
287 return ok ;
288 }
289
290 template <class Field>
run_with_field(Givaro::Integer q,uint64_t b,int m,int k,size_t iters,bool par,uint64_t seed)291 bool run_with_field (Givaro::Integer q, uint64_t b, int m, int k, size_t iters, bool par, uint64_t seed){
292 bool ok = true ;
293
294 int nbit=(int)iters;
295
296 while (ok && nbit){
297 typedef typename Field::Element Element ;
298 // choose Field
299 Field* F= chooseField<Field>(q,b,seed);
300 if (F==nullptr)
301 return true;
302
303 std::ostringstream oss;
304 F->write(oss);
305 std::cout.fill('.');
306 std::cout<<"Checking ";
307 std::cout.width(50);
308 std::cout<<oss.str();
309 std::cout<<" ... ";
310
311 #ifdef __FFLASFFPACK_DEBUG
312 F->write(std::cerr) << std::endl;
313 #endif
314 typedef typename Field::Element Element ;
315 typename Field::RandIter R(*F,b,seed++);
316 typename Field::NonZeroRandIter NZR(R);
317
318 //size_t k = 0 ;
319 //std::cout << k << "/24" << std::endl; ++k;
320 ok = ok && launch_MV_dispatch<Field>(*F,m,k,F->one,F->zero,iters, par, R);
321 //std::cout << k << "/24" << std::endl; ++k;
322 ok = ok && launch_MV_dispatch<Field>(*F,m,k,F->zero,F->zero,iters, par, R);
323 //std::cout << k << "/24" << std::endl; ++k;
324 ok = ok && launch_MV_dispatch<Field>(*F,m,k,F->mOne,F->zero,iters, par, R);
325 //std::cout << k << "/24" << std::endl; ++k;
326 ok = ok && launch_MV_dispatch<Field>(*F,m,k,F->one ,F->one,iters, par, R);
327 //std::cout << k << "/24" << std::endl; ++k;
328 ok = ok && launch_MV_dispatch<Field>(*F,m,k,F->zero,F->one,iters, par, R);
329 //std::cout << k << "/24" << std::endl; ++k;
330 ok = ok && launch_MV_dispatch<Field>(*F,m,k,F->mOne,F->one,iters, par, R);
331 //std::cout << k << "/24" << std::endl; ++k;
332 ok = ok && launch_MV_dispatch<Field>(*F,m,k,F->one ,F->mOne,iters, par, R);
333 //std::cout << k << "/24" << std::endl; ++k;
334 ok = ok && launch_MV_dispatch<Field>(*F,m,k,F->zero,F->mOne,iters, par, R);
335 //std::cout << k << "/24" << std::endl; ++k;
336 ok = ok && launch_MV_dispatch<Field>(*F,m,k,F->mOne,F->mOne,iters, par, R);
337 //std::cout << k << "/24" << std::endl; ++k;
338
339 Element alpha,beta ;
340 NZR.random(alpha);
341 ok = ok && launch_MV_dispatch<Field>(*F,m,k,F->one ,alpha,iters, par, R);
342 //std::cout << k << "/24" << std::endl; ++k;
343 ok = ok && launch_MV_dispatch<Field>(*F,m,k,F->zero,alpha,iters, par, R);
344 //std::cout << k << "/24" << std::endl; ++k;
345 ok = ok && launch_MV_dispatch<Field>(*F,m,k,F->mOne,alpha,iters, par, R);
346 //std::cout << k << "/24" << std::endl; ++k;
347 ok = ok && launch_MV_dispatch<Field>(*F,m,k,alpha,F->one ,iters, par, R);
348 //std::cout << k << "/24" << std::endl; ++k;
349 ok = ok && launch_MV_dispatch<Field>(*F,m,k,alpha,F->zero,iters, par, R);
350 //std::cout << k << "/24" << std::endl; ++k;
351 ok = ok && launch_MV_dispatch<Field>(*F,m,k,alpha,F->mOne,iters, par, R);
352 //std::cout << k << "/24" << std::endl; ++k;
353
354 for (size_t j = 0 ; j < 3 ; ++j) {
355 R.random(alpha);
356 R.random(beta);
357 ok = ok && launch_MV_dispatch<Field>(*F,m,k,alpha,beta,iters, par, R);
358 //std::cout << k << "/24" << std::endl; ++k;
359 }
360 //std::cout<<std::endl;
361 nbit--;
362 if ( !ok )
363 //std::cout << "\033[1;31mFAILED\033[0m "<<std::endl;
364 std::cout << "FAILED "<<std::endl;
365 else
366 //std::cout << "\033[1;32mPASSED\033[0m "<<std::endl;
367 std::cout << "PASSED "<<std::endl;
368 delete F;
369 }
370 return ok;
371 }
main(int argc,char ** argv)372 int main(int argc, char** argv)
373 {
374 std::cout<<setprecision(17);
375 std::cerr<<setprecision(17);
376
377 uint64_t seed = getSeed();
378 size_t iters = 3 ;
379 Givaro::Integer q = -1 ;
380 uint64_t b = 0 ;
381 int m = -50 ;
382 int k = -50 ;
383 int nbw = -1 ;
384 bool loop = false;
385 bool p = false;
386 Argument as[] = {
387 { 'q', "-q Q", "Set the field characteristic (-1 for random).", TYPE_INTEGER , &q },
388 { 'b', "-b B", "Set the bitsize of the random characteristic.", TYPE_INT , &b },
389 { 'm', "-m M", "Set the dimension of the matrix (negative values, mean, any random value between 0 and |n|).", TYPE_INT , &m },
390 { 'k', "-k K", "Set the dimension of the matrix (negative values, mean, any random value between 0 and |k|).", TYPE_INT , &k },
391 { 'w', "-w N", "Set the number of winograd levels (-1 for random).", TYPE_INT , &nbw },
392 { 'i', "-i R", "Set number of repetitions.", TYPE_INT , &iters },
393 { 'l', "-l Y/N", "run the test in an infinte loop.", TYPE_BOOL , &loop },
394 { 'p', "-p Y/N", "run the parallel fgemv.", TYPE_BOOL , &p },
395 { 's', "-s seed", "Set seed for the random generator", TYPE_UINT64, &seed },
396 END_OF_ARGUMENTS
397 };
398
399 parseArguments(argc,argv,as);
400
401 bool ok = true;
402 srand(seed);
403 do{
404 ok = ok && run_with_field<Modular<double> >(q,b,m,k,iters,p, seed);
405 ok = ok && run_with_field<ModularBalanced<double> >(q,b,m,k,iters,p, seed);
406 ok = ok && run_with_field<Modular<float> >(q,b,m,k,iters,p, seed);
407 ok = ok && run_with_field<ModularBalanced<float> >(q,b,m,k,iters,p, seed);
408 ok = ok && run_with_field<Modular<int32_t> >(q,b,m,k,iters,p, seed);
409 ok = ok && run_with_field<ModularBalanced<int32_t> >(q,b,m,k,iters,p, seed);
410 ok = ok && run_with_field<Modular<int64_t> >(q,b,m,k,iters, p, seed);
411 ok = ok && run_with_field<ModularBalanced<int64_t> >(q,b,m,k,iters, p, seed);
412
413 ok = ok && run_with_field<Modular<RecInt::rint<8> > >(q,b?b:127_ui64,m,k,iters, p, seed);
414 ok = ok && run_with_field<Modular<RecInt::ruint<7>,RecInt::ruint<8> > >(q,b?b:127_ui64,m,k,iters, p, seed);
415 ok = ok && run_with_field<Modular<Givaro::Integer> >(q,(b?b:512_ui64),m,k,iters,p, seed);
416 ok = ok && run_with_field<Givaro::ZRing<Givaro::Integer> >(0,(b?b:512_ui64),m,k,iters,p, seed);
417 } while (loop && ok);
418
419
420
421
422 return !ok ;
423 }
424 /* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
425 // vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
426