1 /*
2 * Copyright (C) 2016 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 #define ENABLE_ALL_CHECKINGS 1
27
28 #include "fflas-ffpack/fflas-ffpack-config.h"
29
30 #include <iomanip>
31 #include <iostream>
32 #include <random>
33
34 #include "fflas-ffpack/utils/fflas_io.h"
35 #include "fflas-ffpack/utils/timer.h"
36 #include "fflas-ffpack/fflas/fflas.h"
37 #include "fflas-ffpack/utils/args-parser.h"
38 #include "fflas-ffpack/utils/test-utils.h"
39 #include <givaro/modular.h>
40
41
42 using namespace std;
43 using namespace FFLAS;
44 using Givaro::Modular;
45 using Givaro::ModularBalanced;
46
47 template<typename Field, class RandIter>
check_fsyrk(const Field & F,size_t n,size_t k,const typename Field::Element & alpha,const typename Field::Element & beta,FFLAS::FFLAS_UPLO uplo,FFLAS::FFLAS_TRANSPOSE trans,RandIter & Rand)48 bool check_fsyrk (const Field &F, size_t n, size_t k,
49 const typename Field::Element &alpha, const typename Field::Element &beta,
50 FFLAS::FFLAS_UPLO uplo, FFLAS::FFLAS_TRANSPOSE trans, RandIter& Rand){
51
52 typedef typename Field::Element Element;
53 Element * A, *C, *C2;
54 size_t ldc = n+15;
55 size_t Arows = (trans==FFLAS::FflasNoTrans)?n:k;
56 size_t Acols = (trans==FFLAS::FflasNoTrans)?k:n;
57 size_t lda = Acols+13;
58
59 A = FFLAS::fflas_new(F,Arows,lda);
60 C = FFLAS::fflas_new(F,n,ldc);
61 C2 = FFLAS::fflas_new(F,n,ldc);
62
63 FFPACK::RandomTriangularMatrix (F, n, n, uplo, FflasNonUnit, true, C, ldc, Rand);
64 FFPACK::RandomMatrix (F, Arows, Acols, A, lda, Rand);
65 FFLAS::fassign (F, n, n, C, ldc, C2, ldc);
66
67 string ss=string((uplo == FFLAS::FflasLower)?"Lower_":"Upper_")+string((trans == FFLAS::FflasTrans)?"Trans":"NoTrans");
68
69 cout<<std::left<<"Checking FSYRK_";
70 cout.fill('.');
71 cout.width(35);
72 cout<<ss;
73
74 FFLAS::Timer t; t.clear();
75 double time=0.0;
76 t.clear(); t.start();
77
78 fsyrk (F, uplo, trans, n, k, alpha, A, lda, beta, C, ldc);
79
80 t.stop();
81 time+=t.usertime();
82
83 fgemm (F, trans, (trans==FflasNoTrans)?FflasTrans:FflasNoTrans, n, n, k, alpha, A, lda, A, lda, beta, C2, ldc);
84
85 bool ok = true;
86 if (uplo == FflasUpper){
87 for (size_t i=0; i<n; i++)
88 for (size_t j=i; j<n; j++)
89 ok = ok && F.areEqual(C2[i*ldc+j], C[i*ldc+j]);
90 } else {
91 for (size_t i=0; i<n; i++)
92 for (size_t j=0; j<=i; j++)
93 ok = ok && F.areEqual(C2[i*ldc+j], C[i*ldc+j]);
94 }
95 if (ok)
96 //cout << "\033[1;32mPASSED\033[0m ("<<time<<")"<<endl;
97 cout << "PASSED ("<<time<<")"<<endl;
98 //cerr<<"PASSED ("<<time<<")"<<endl;
99 else
100 //cout << "\033[1;31mFAILED\033[0m ("<<time<<")"<<endl;
101 cout << "FAILED ("<<time<<")"<<endl;
102 //cerr<<"FAILED ("<<time<<")"<<endl;
103
104 FFLAS::fflas_delete(A);
105 FFLAS::fflas_delete(C2);
106 FFLAS::fflas_delete(C);
107 return ok;
108 }
109 template<typename Field, class RandIter>
check_fsyrk_diag(const Field & F,size_t n,size_t k,const typename Field::Element & alpha,const typename Field::Element & beta,FFLAS::FFLAS_UPLO uplo,FFLAS::FFLAS_TRANSPOSE trans,RandIter & Rand)110 bool check_fsyrk_diag (const Field &F, size_t n, size_t k,
111 const typename Field::Element &alpha, const typename Field::Element &beta,
112 FFLAS::FFLAS_UPLO uplo, FFLAS::FFLAS_TRANSPOSE trans, RandIter& Rand){
113
114 typedef typename Field::Element Element;
115 Element * A, *B, *C, *C2, *D;
116 size_t ldc = n+(rand()%50);
117 size_t Arows = (trans==FFLAS::FflasNoTrans)?n:k;
118 size_t Acols = (trans==FFLAS::FflasNoTrans)?k:n;
119 size_t lda = Acols+(rand()%50);
120 size_t incD = 1+(rand()%100);
121 A = FFLAS::fflas_new(F,Arows,lda);
122 B = FFLAS::fflas_new(F,Arows,lda);
123 C = FFLAS::fflas_new(F,n,ldc);
124 C2 = FFLAS::fflas_new(F,n,ldc);
125 D = FFLAS::fflas_new(F,k,incD);
126 Givaro::GeneralRingNonZeroRandIter<Field,RandIter> nzRand (Rand);
127 for (size_t i=0; i<k; i++)
128 nzRand.random(D[i*incD]);
129 FFPACK::RandomTriangularMatrix (F, n, n, uplo, FflasNonUnit, true, C, ldc, Rand);
130 FFPACK::RandomMatrix (F, Arows, Acols, A, lda, Rand);
131 FFLAS::fassign (F, n, n, C, ldc, C2, ldc);
132 FFLAS::fassign (F, Arows, Acols, A, lda, B, lda);
133
134 string ss=string((uplo == FFLAS::FflasLower)?"Lower_":"Upper_")+string((trans == FFLAS::FflasTrans)?"Trans":"NoTrans");
135
136 cout<<std::left<<"Checking FSYRK_DIAG_";
137 cout.fill('.');
138 cout.width(30);
139 cout<<ss;
140
141
142 // FFLAS::WriteMatrix ( std::cerr<<"A = "<<std::endl,F,Arows, Acols, A, lda);
143 // FFLAS::WriteMatrix ( std::cerr<<"C = "<<std::endl,F,n,n, C, ldc);
144 // FFLAS::WriteMatrix ( std::cerr<<"D = "<<std::endl,F,k,1,D, incD);
145 FFLAS::Timer t; t.clear();
146 double time=0.0;
147 t.clear(); t.start();
148
149 fsyrk (F, uplo, trans, n, k, alpha, A, lda, D, incD, beta, C, ldc, 13);
150
151 t.stop();
152 time+=t.usertime();
153
154 // std::cerr<<"After fsyrk_diag"<<std::endl;
155 // FFLAS::WriteMatrix (std::cerr<<"A = "<<std::endl,F,Arows, Acols, A, lda);
156 // FFLAS::WriteMatrix (std::cerr<<"C = "<<std::endl,F,n,n,C,ldc);
157
158 bool ok = true;
159
160 typename Field::Element tmp;
161 F.init(tmp);
162 if (trans==FflasNoTrans){
163 // Checking whether A = B x D
164 for (size_t i=0; i < Arows; i++)
165 for (size_t j=0; j < Acols; j++)
166 ok = ok && F.areEqual(A[i*lda+j], F.mul (tmp, B[i*lda+j], D[j*incD]));
167 } else {
168 // Checking whether A = D x B
169 for (size_t i=0; i < Arows; i++)
170 for (size_t j=0; j < Acols; j++){
171 if(!F.areEqual(A[i*lda+j], F.mul (tmp, B[i*lda+j], D[i*incD]))){
172 std::cerr<<"B["<<i<<","<<j<<"] = "<<B[i*lda+j]<<" != "<<D[i*incD]<<" * "<<A[i*lda+j]<<std::endl;
173 ok=false;
174 }
175 }
176 }
177 if (!ok){
178 std::cerr<<"Scaling failed"<<std::endl;
179 return ok;
180 }
181
182 fgemm (F, trans, (trans==FflasNoTrans)?FflasTrans:FflasNoTrans, n, n, k, alpha, A, lda, B, lda, beta, C2, ldc);
183
184 if (uplo == FflasUpper){
185 for (size_t i=0; i<n; i++)
186 for (size_t j=i; j<n; j++){
187 if(!F.areEqual(C2[i*ldc+j], C[i*ldc+j])){
188 std::cerr<<"C2["<<i<<","<<j<<"] = "<<C2[i*ldc+j]<<" != C["<<i<<","<<j<<"] = "<< C[i*ldc+j]<<std::endl;
189 ok=false;
190 }
191 }
192 } else {
193 for (size_t i=0; i<n; i++)
194 for (size_t j=0; j<=i; j++)
195 ok = ok && F.areEqual(C2[i*ldc+j], C[i*ldc+j]);
196 }
197 if (ok)
198 cout << "PASSED ("<<time<<")"<<endl;
199 else
200 cout << "FAILED ("<<time<<")"<<endl;
201
202 FFLAS::fflas_delete(A);
203 FFLAS::fflas_delete(C2);
204 FFLAS::fflas_delete(C);
205 FFLAS::fflas_delete(D);
206 return ok;
207 }
208 template<typename Field, class RandIter>
check_fsyrk_bkdiag(const Field & F,size_t n,size_t k,const typename Field::Element & alpha,const typename Field::Element & beta,FFLAS::FFLAS_UPLO uplo,FFLAS::FFLAS_TRANSPOSE trans,RandIter & Rand)209 bool check_fsyrk_bkdiag (const Field &F, size_t n, size_t k,
210 const typename Field::Element &alpha, const typename Field::Element &beta,
211 FFLAS::FFLAS_UPLO uplo, FFLAS::FFLAS_TRANSPOSE trans, RandIter& Rand){
212
213 typedef typename Field::Element Element;
214 Element * A, *B, *C, *C2, *D;
215 size_t ldc = n+(rand()%50);
216 size_t Arows = (trans==FFLAS::FflasNoTrans)?n:k;
217 size_t Acols = (trans==FFLAS::FflasNoTrans)?k:n;
218 size_t lda = Acols+(rand()%50);
219 size_t incD = 1+(rand()%100);
220 A = FFLAS::fflas_new(F,Arows,lda);
221 B = FFLAS::fflas_new(F,Arows,lda);
222 C = FFLAS::fflas_new(F,n,ldc);
223 C2 = FFLAS::fflas_new(F,n,ldc);
224 D = FFLAS::fflas_new(F,k,incD);
225 Givaro::GeneralRingNonZeroRandIter<Field,RandIter> nzRand (Rand);
226 std::vector<bool> tb(k,false);
227 for (size_t i=0; i<k; i++){
228 if (rand()% 10 != 0 || i==k-1)
229 nzRand.random(D[i*incD]);
230 else{
231 nzRand.random(D[i*incD]);
232 F.assign(D[(i+1)*incD],D[i*incD]);
233 tb[i]=true;
234 i++;
235 }
236 }
237
238 FFPACK::RandomTriangularMatrix (F, n, n, uplo, FflasNonUnit, true, C, ldc, Rand);
239 FFPACK::RandomMatrix (F, Arows, Acols, A, lda, Rand);
240 FFLAS::fassign (F, n, n, C, ldc, C2, ldc);
241 FFLAS::fassign (F, Arows, Acols, A, lda, B, lda);
242
243 string ss=string((uplo == FFLAS::FflasLower)?"Lower_":"Upper_")+string((trans == FFLAS::FflasTrans)?"Trans":"NoTrans");
244
245 cout<<std::left<<"Checking FSYRK_BK_DIAG_";
246 cout.fill('.');
247 cout.width(30);
248 cout<<ss;
249
250
251 // FFLAS::WriteMatrix ( std::cerr<<"A = "<<std::endl,F,Arows, Acols, A, lda);
252 // FFLAS::WriteMatrix ( std::cerr<<"C = "<<std::endl,F,n,n, C, ldc);
253 // FFLAS::WriteMatrix ( std::cerr<<"D = "<<std::endl,F,k,1,D, incD);
254 FFLAS::Timer t; t.clear();
255 double time=0.0;
256 t.clear(); t.start();
257
258 fsyrk (F, uplo, trans, n, k, alpha, A, lda, D, incD, tb, beta, C, ldc, 13);
259
260 t.stop();
261 time+=t.usertime();
262
263 // std::cerr<<"After fsyrk_bk_diag"<<std::endl;
264 // FFLAS::WriteMatrix (std::cerr<<"A = "<<std::endl,F,Arows, Acols, A, lda);
265 // FFLAS::WriteMatrix (std::cerr<<"C = "<<std::endl,F,n,n,C,ldc);
266 bool ok = true;
267
268 typename Field::Element tmp;
269 F.init(tmp);
270 if (trans==FflasNoTrans){
271 // Checking whether A = B x D
272 for (size_t i=0; i < Arows; i++)
273 for (size_t j=0; j < Acols; j++)
274 if (!tb[j])
275 ok = ok && F.areEqual(A[i*lda+j], F.mul (tmp, B[i*lda+j], D[j*incD]));
276 else{
277 ok = ok && F.areEqual(A[i*lda+j], F.mul (tmp, B[i*lda+j+1], D[j*incD]));
278 ok = ok && F.areEqual(A[i*lda+j+1], F.mul (tmp, B[i*lda+j], D[j*incD]));
279 j++;
280 }
281
282 } else {
283 // Checking whether A = D x B
284 for (size_t j=0; j < Acols; j++){
285 for (size_t i=0; i < Arows; i++)
286 if (!tb[i])
287 ok = ok && F.areEqual(A[i*lda+j], F.mul (tmp, B[i*lda+j], D[i*incD]));
288 else{
289 ok = ok && F.areEqual(A[i*lda+j], F.mul (tmp, B[(i+1)*lda+j], D[i*incD]));
290 ok = ok && F.areEqual(A[(i+1)*lda+j], F.mul (tmp, B[i*lda+j], D[i*incD]));
291 i++;
292 }
293 // {
294 // std::cerr<<"B["<<i<<","<<j<<"] = "<<B[i*lda+j]<<" != "<<D[i*incD]<<" * "<<A[i*lda+j]<<std::endl;
295 // ok=false;
296 // }
297 }
298 }
299 if (!ok){
300 std::cerr<<"Scaling failed"<<std::endl;
301 std::cerr<<"alpha = "<<alpha<<" beta="<<beta<<std::endl;
302 std::cerr<<"tb = "<<tb<<std::endl;
303
304 return ok;
305 }
306
307 fgemm (F, trans, (trans==FflasNoTrans)?FflasTrans:FflasNoTrans, n, n, k, alpha, A, lda, B, lda, beta, C2, ldc);
308
309 if (uplo == FflasUpper){
310 for (size_t i=0; i<n; i++)
311 for (size_t j=i; j<n; j++){
312 if(!F.areEqual(C2[i*ldc+j], C[i*ldc+j])){
313 std::cerr<<"C2["<<i<<","<<j<<"] = "<<C2[i*ldc+j]<<" != C["<<i<<","<<j<<"] = "<< C[i*ldc+j]<<std::endl;
314 ok=false;
315 }
316 }
317 } else {
318 for (size_t i=0; i<n; i++)
319 for (size_t j=0; j<=i; j++)
320 ok = ok && F.areEqual(C2[i*ldc+j], C[i*ldc+j]);
321 }
322 if (ok)
323 //cout << "\033[1;32mPASSED\033[0m ("<<time<<")"<<endl;
324 cout << "PASSED ("<<time<<")"<<endl;
325 //cerr<<"PASSED ("<<time<<")"<<endl;
326 else
327 //cout << "\033[1;31mFAILED\033[0m ("<<time<<")"<<endl;
328 cout << "FAILED ("<<time<<")"<<endl;
329 //cerr<<"FAILED ("<<time<<")"<<endl;
330
331 FFLAS::fflas_delete(A);
332 FFLAS::fflas_delete(C2);
333 FFLAS::fflas_delete(C);
334 FFLAS::fflas_delete(D);
335 return ok;
336 }
337
338 template <class Field>
run_with_field(Givaro::Integer q,size_t b,size_t n,size_t k,int a,int c,size_t iters,uint64_t seed)339 bool run_with_field (Givaro::Integer q, size_t b, size_t n, size_t k, int a, int c, size_t iters, uint64_t seed){
340 bool ok = true ;
341 int nbit=(int)iters;
342
343 while (ok && nbit){
344 //typedef typename Field::Element Element ;
345 // choose Field
346 Field* F= FFPACK::chooseField<Field>(q,b,seed);
347 typename Field::RandIter G(*F,0,seed++);
348 if (F==nullptr)
349 return true;
350
351 typename Field::Element alpha, beta;
352 F->init (alpha, (typename Field::Element)a);
353 F->init (beta, (typename Field::Element)c);
354 cout<<"Checking with ";F->write(cout)<<endl;
355
356 ok = ok && check_fsyrk(*F,n,k,alpha,beta,FflasUpper,FflasNoTrans,G);
357 ok = ok && check_fsyrk(*F,n,k,alpha,beta,FflasUpper,FflasTrans,G);
358 ok = ok && check_fsyrk(*F,n,k,alpha,beta,FflasLower,FflasNoTrans,G);
359 ok = ok && check_fsyrk(*F,n,k,alpha,beta,FflasLower,FflasTrans,G);
360 ok = ok && check_fsyrk_diag(*F,n,k,alpha,beta,FflasUpper,FflasNoTrans,G);
361 ok = ok && check_fsyrk_diag(*F,n,k,alpha,beta,FflasUpper,FflasTrans,G);
362 ok = ok && check_fsyrk_diag(*F,n,k,alpha,beta,FflasLower,FflasNoTrans,G);
363 ok = ok && check_fsyrk_diag(*F,n,k,alpha,beta,FflasLower,FflasTrans,G);
364 ok = ok && check_fsyrk_bkdiag(*F,n,k,alpha,beta,FflasUpper,FflasNoTrans,G);
365 ok = ok && check_fsyrk_bkdiag(*F,n,k,alpha,beta,FflasUpper,FflasTrans,G);
366 ok = ok && check_fsyrk_bkdiag(*F,n,k,alpha,beta,FflasLower,FflasNoTrans,G);
367 ok = ok && check_fsyrk_bkdiag(*F,n,k,alpha,beta,FflasLower,FflasTrans,G);
368
369 // checking with k > n (=k+n)
370 ok = ok && check_fsyrk(*F,n,k+n,alpha,beta,FflasUpper,FflasNoTrans,G);
371 ok = ok && check_fsyrk(*F,n,k+n,alpha,beta,FflasUpper,FflasTrans,G);
372 ok = ok && check_fsyrk(*F,n,k+n,alpha,beta,FflasLower,FflasNoTrans,G);
373 ok = ok && check_fsyrk(*F,n,k+n,alpha,beta,FflasLower,FflasTrans,G);
374 ok = ok && check_fsyrk_diag(*F,n,k+n,alpha,beta,FflasUpper,FflasNoTrans,G);
375 ok = ok && check_fsyrk_diag(*F,n,k+n,alpha,beta,FflasUpper,FflasTrans,G);
376 ok = ok && check_fsyrk_diag(*F,n,k+n,alpha,beta,FflasLower,FflasNoTrans,G);
377 ok = ok && check_fsyrk_diag(*F,n,k+n,alpha,beta,FflasLower,FflasTrans,G);
378 ok = ok && check_fsyrk_bkdiag(*F,n,k+n,alpha,beta,FflasUpper,FflasNoTrans,G);
379 ok = ok && check_fsyrk_bkdiag(*F,n,k+n,alpha,beta,FflasUpper,FflasTrans,G);
380 ok = ok && check_fsyrk_bkdiag(*F,n,k+n,alpha,beta,FflasLower,FflasNoTrans,G);
381 ok = ok && check_fsyrk_bkdiag(*F,n,k+n,alpha,beta,FflasLower,FflasTrans,G);
382 nbit--;
383 delete F;
384 }
385 return ok;
386 }
387
main(int argc,char ** argv)388 int main(int argc, char** argv)
389 {
390 cerr<<setprecision(10);
391 Givaro::Integer q=-1;
392 size_t b=0;
393 int k=35;
394 int n=109;
395 int a=-1;
396 int c=1;
397 size_t iters=3;
398 bool loop=false;
399 uint64_t seed = getSeed();
400
401 Argument as[] = {
402 { 'q', "-q Q", "Set the field characteristic (-1 for random).", TYPE_INTEGER , &q },
403 { 'b', "-b B", "Set the bitsize of the field characteristic.", TYPE_INT , &b },
404 { 'k', "-k K", "Set the dimension", TYPE_INT , &k },
405 { 'n', "-n N", "Set the column dimension.", TYPE_INT , &n },
406 { 'a', "-a A", "Set the scaling alpha", TYPE_INT , &a },
407 { 'c', "-c C", "Set the scaling beta", TYPE_INT , &c },
408 { 'i', "-i R", "Set number of repetitions.", TYPE_INT , &iters },
409 { 'l', "-loop Y/N", "run the test in an infinite loop.", TYPE_BOOL , &loop },
410 { 's', "-s seed", "Set seed for the random generator", TYPE_UINT64, &seed },
411 END_OF_ARGUMENTS
412 };
413
414 parseArguments(argc,argv,as);
415 srand(seed);
416 bool ok = true;
417 do{
418 ok = ok && run_with_field<Modular<double> >(q,b,n,k,a,c,iters,seed);
419 ok = ok && run_with_field<ModularBalanced<double> >(q,b,n,k,a,c,iters,seed);
420 ok = ok && run_with_field<Modular<float> >(q,b,n,k,a,c,iters,seed);
421 ok = ok && run_with_field<ModularBalanced<float> >(q,b,n,k,a,c,iters,seed);
422 ok = ok && run_with_field<Modular<int32_t> >(q,b,n,k,a,c,iters,seed);
423 ok = ok && run_with_field<ModularBalanced<int32_t> >(q,b,n,k,a,c,iters,seed);
424 ok = ok && run_with_field<Modular<int64_t> >(q,b,n,k,a,c,iters,seed);
425 ok = ok && run_with_field<ModularBalanced<int64_t> >(q,b,n,k,a,c,iters,seed);
426 ok = ok && run_with_field<Modular<Givaro::Integer> >(q,5,n/4+1,k/4+1,a,c,iters,seed);
427 ok = ok && run_with_field<Modular<Givaro::Integer> >(q,(b?b:512),n/4+1,k/4+1,a,c,iters,seed);
428 } while (loop && ok);
429
430 if (!ok) std::cerr<<"with seed = "<<seed<<std::endl;
431
432 return !ok ;
433 }
434 /* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
435 // vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
436