1 /* tests/test-blas-domain.C
2 * Copyright (C) 2004 Pascal Giorgi
3 *
4 * Written by Pascal Giorgi <pascal.giorgi@ens-lyon.fr>
5 *
6 * ---------------------------------------------------------
7 *
8 *
9 * ========LICENCE========
10 * This file is part of the library LinBox.
11 *
12 * LinBox is free software: you can redistribute it and/or modify
13 * it under the terms of the GNU Lesser General Public
14 * License as published by the Free Software Foundation; either
15 * version 2.1 of the License, or (at your option) any later version.
16 *
17 * This library is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
21 *
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with this library; if not, write to the Free Software
24 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
25 * ========LICENCE========
26 *.
27 *
28 *
29 */
30
31
32 /*! @file tests/test-blas-domain.C
33 * @ingroup tests
34 * @brief no doc
35 * @test NO DOC
36 */
37
38
39 #include <linbox/linbox-config.h>
40 #include <iostream>
41 #include <string>
42
43 #include "linbox/integer.h"
44 //#include "linbox/field/gf2.h"
45 //#include "linbox/domain/blas_matrix_domain_gf2.h"
46 // defines template<> class BlasMatrixDomain<GF2> { ... };
47 // ... and template<> BlasMatrix<GF2> {}
48 #if 0
49 #include "linbox/field/gf3.h"
50 //#include "linbox/domain/blas_matrix_domain_gf3.h"
51 // defines template<> class BlasMatrixDomain<GF3> { ... };
52 // ... and template<> BlasMatrix<GF3> {}
53 #endif
54 #include "linbox/ring/modular.h"
55 #include <givaro/modular-balanced.h>
56
57 #ifdef __LINBOX_HAVE_NTL
58 #include "linbox/ring/ntl.h"
59 #endif
60 #include "linbox/matrix/dense-matrix.h"
61 #include "linbox/matrix/matrix-domain.h"
62 #include <givaro/givranditer.h>
63 #include "linbox/util/commentator.h"
64 #include "givaro/zring.h"
65 // #include "linbox/algorithms/matrix-hom.h"
66
67 #include "linbox/matrix/random-matrix.h"
68 #include "linbox/blackbox/scalar-matrix.h"
69
70
71
72 #include "test-common.h"
73
74 using namespace LinBox;
75
76 const int maxpretty = 35;
77
78 string blank;
79
80 const char* pretty(string a)
81 ;
82
83 template <class Vector>
84 bool localAreEqual(
85 const Vector& a,
86 const Vector& b)
87 ;
88
89 template <class Field>
90 static bool testMulAdd (const Field& F, size_t n, int iterations)
91 ;
92 template <class Field>
93 static bool testMulAddAgain (const Field& Zp, size_t n, int iterations)
94 ;
95 template <class Field>
96 static bool testMulAddShapeTrans (const Field &F, size_t m, size_t n, size_t k, int iterations)
97 ;
98 template<class Field, bool LeftSide, bool UnitDiag>
99 static bool testTriangMulShapeTrans (const Field &F, size_t m, size_t n, int iterations)
100 ;
101 template <class Field>
102 static bool testRank (const Field& F,size_t n, int iterations)
103 ;
104 template <class Field>
105 static bool testDet (const Field& F,size_t n, int iterations)
106 ;
107 template <class Field>
108 static bool testInv (const Field& F,size_t n, int iterations)
109 ;
110 template <class Field>
111 static bool testTriangularSolve (const Field& F, size_t m, size_t n, int iterations)
112 ;
113 template <class Field>
114 static bool testSolve (const Field& F, size_t m, size_t n, int iterations)
115 ;
116 template <class Field>
117 static bool testPermutation (const Field& F, size_t m, int iterations)
118 ;
119 template <class Field>
120 static bool testPLUQ (const Field& F, size_t m, size_t n, int iterations)
121 ;
122 template <class Field>
123 static bool testMinPoly (const Field& F, size_t n, int iterations)
124 ;
125 template <class Field>
126 static bool testCharPoly (const Field& F, size_t n, int iterations)
127 ;
128 template<class Field>
129 static bool testBlasMatrixConstructors(const Field& Fld, size_t m, size_t n)
130 ;
131 template<class Field>
132 int launch_tests(Field & F, size_t n, int iterations)
133 ;
134 bool launch_gf2_tests(GF2 & F, size_t n)
135 ;
136 #if 0
137 bool launch_gf3_tests(GF3 & F, size_t n)
138 ;
139 #endif
140
141 // test BlasMatrixDomain<Field> for various fields
main(int argc,char ** argv)142 int main(int argc, char **argv)
143 {
144
145 static size_t n = 3;
146 //static integer q = 1000003U;
147 static integer q = 65521;
148 static int iterations = 1;
149
150 static Argument args[] = {
151 { 'n', "-n N", "Set dimension of test matrices to NxN", TYPE_INT, &n },
152 { 'q', "-q Q", "Operate over the \"field\" GF(Q) [1]", TYPE_INTEGER, &q },
153 { 'i', "-i I", "Perform each test for I iterations", TYPE_INT, &iterations },
154 END_OF_ARGUMENTS
155 };
156
157 parseArguments (argc, argv, args);
158
159 // BB. (Blas)MatrixDomain are not very generic...
160
161 bool pass = true;
162 srand ((unsigned)time (NULL));
163
164 commentator().getMessageClass (INTERNAL_DESCRIPTION).setMaxDepth (3);
165 commentator().getMessageClass (INTERNAL_DESCRIPTION).setMaxDetailLevel (Commentator::LEVEL_NORMAL);
166 commentator().start("BlasMatrixDomain test suite", "BlasMatrixDomain");
167
168 Givaro::Modular<double> F1 (q);
169 GF2 F2 ;
170 //GF3 F3 ;
171 Givaro::ModularBalanced<double> F4(q);
172 Givaro::Modular<float> F5(2011);
173 Givaro::Modular<uint32_t> F6(1009); // (2011);
174
175 pass &= launch_tests(F1,n,iterations);
176
177 pass &= launch_gf2_tests(F2,n);
178
179 #pragma message "#warning GF3 -> working on sliced wrapper"
180 //pass &= launch_gf3_tests(F3,n);
181
182 pass &= launch_tests(F4,n,iterations);
183
184 pass &= launch_tests(F5,n,iterations);
185
186 pass &= launch_tests(F6,n,iterations);
187
188 //#pragma message "#warning Givaro::Modular<bool> is not working"
189 //pass &= launch_tests(F7,n,iterations);
190
191 #ifdef __LINBOX_HAVE_NTL
192 #pragma message "#warning NTL_ZZp is not working at all"
193 NTL::ZZ_p::init(NTL::to_ZZ((size_t)q));
194 NTL_ZZ_p F8;
195 // pass &= launch_tests(F8,n,iterations);
196 #endif
197
198 commentator().stop(MSG_STATUS (pass), (const char *) 0,"BlasMatrixDomain test suite");
199 return pass ? 0 : -1;
200 }
201
pretty(string a)202 const char* pretty(string a)
203 {
204
205 blank = " " + a;
206 int msgsize= maxpretty - (int)blank.size();
207 string dot(".");
208 for (int i=0;i<msgsize ;++i)
209 blank+=dot;
210 return blank.c_str();
211 }
212 #define mycommentator commentator
213 ostream & report = mycommentator().report();
214
215 template<class Vector>
localAreEqual(const Vector & a,const Vector & b)216 bool localAreEqual(
217 const Vector& a,
218 const Vector& b)
219 {
220 if ( a.size() != b.size() )
221 return false;
222
223 const typename Vector::Field & F = a.field();
224
225 typename Vector::const_iterator it = a.begin();
226 typename Vector::const_iterator jt = b.begin();
227 for ( ; it != a.end(); ++it, ++jt)
228 if (! F.areEqual(*it,*jt))
229 return false;
230 return true;
231 }
232
233 template <class Field>
testMulAdd(const Field & F,size_t n,int iterations)234 static bool testMulAdd (const Field& F, size_t n, int iterations)
235 {
236
237 typedef typename Field::Element Element;
238 typedef typename Field::RandIter RandIter;
239 typedef BlasMatrix<Field> Matrix;
240
241 //Commentator mycommentator;
242 //mycommentator().getMessageClass (INTERNAL_DESCRIPTION).setMaxDepth (3);
243 //mycommentator().getMessageClass (INTERNAL_DESCRIPTION).setMaxDetailLevel (Commentator::LEVEL_NORMAL);
244 mycommentator().start (pretty("Testing muladd"),"testMulAdd",(unsigned int)iterations);
245
246 RandIter G(F);
247 bool ret = true;
248 BlasMatrixDomain<Field> BMD(F);
249
250 Matrix A(F,n,n);
251 for (int k=0;k<iterations; ++k) {
252
253 mycommentator().progress(k);
254 // A.init(F, n, n);
255 Matrix /*A(F, n,n),*/B(F, n,n),C(F, n,n),D(F, n,n),T(F, n,n),R(F, n,n);
256 BlasVector<Field> x(F,n),y(F,n),z(F,n),t(F,n);
257
258 Element alpha, beta,malpha;
259
260
261 // Create 3 random n*n matrices
262 A.random();
263 B.random();
264 C.random();
265
266 // Create 2 random vectors
267 for (size_t i=0;i<n;++i) {
268 G.random(x[i]);
269 G.random(y[i]);
270 }
271
272 // create 2 random element
273 G.random(alpha);
274 G.random(beta);
275
276 F.neg(malpha,alpha);
277
278 // compute D = -alpha.(A*C+B*C) + alpha.(A+B)*C
279
280 BMD.mul(D,A,C);
281 BMD.mul(T,B,C);
282 BMD.addin(D,T);
283
284 BMD.add(T,A,B);
285 BMD.muladd(R,malpha,D,alpha,T,C);
286
287 if (!BMD.isZero(R)) {
288 ret=false;
289 }
290
291 // compute z = beta.y + alpha.A*x
292 // report << "beta := " << beta << ";" << std::endl;
293 // report << "y := " <<y << ";" << std::endl;
294 // report << "alpha := " << alpha << ";" << std::endl;
295 // A.write(report << "A :=",Tag::FileFormat::Maple) << ";" << std::endl;
296 // report << "x := " << x << std::endl;
297
298 BMD.muladd(z,beta,y,alpha,A,x);
299 // report << "z := "<< z << std::endl;
300
301 // report << "#apply" <<std::endl;
302 A.apply(t, x);
303 //MD.vectorMul(t,A,x);
304 for (size_t i=0;i<n;++i){
305 F.mulin(t[i],alpha);
306 F.axpyin(t[i],beta,y[i]);
307 }
308 // report << "t := "<< t << std::endl;
309
310 if (!localAreEqual(t,z)){
311 exit(-1);
312 report << "2 alpha = " << alpha << "mod " << F.characteristic() << std::endl;
313 //report << "2 alpha = " << alpha << "mod " << F.characteristic() << std::endl;
314 ret=false;
315 }
316 }
317
318 mycommentator().stop(MSG_STATUS (ret), (const char *) 0, "testMulAdd");
319
320 return ret;
321 }
322
323 // computes D = alpha A B + beta C on integers and check the result is ok mod p.
324 // actually we check the mod p muladd here...
325 template <class Field>
CheckMulAdd(const Field & Zp,const Integer & alpha,const BlasMatrix<Givaro::ZRing<Integer>> & A,const BlasMatrix<Givaro::ZRing<Integer>> & B,const Integer & beta,const BlasMatrix<Givaro::ZRing<Integer>> & C)326 bool CheckMulAdd(const Field& Zp, const Integer & alpha ,
327 const BlasMatrix<Givaro::ZRing<Integer> > & A ,
328 const BlasMatrix<Givaro::ZRing<Integer> > & B ,
329 const Integer & beta ,
330 const BlasMatrix<Givaro::ZRing<Integer> > & C)
331 {
332
333 size_t M = C.rowdim();
334 size_t N = C.coldim();
335
336 //typedef Givaro::Modular<double> Field ;
337 typedef typename Field::Element Element ;
338
339 Givaro::ZRing<Integer> ZZ ;
340 // compiles, but wrong if BlasMatrixDomain is used
341 MatrixDomain<Givaro::ZRing<Integer> > ZMD(ZZ);
342
343 BlasMatrix<Givaro::ZRing<Integer> > D(ZZ,M,N);
344
345 /*
346 Integer p = Integer::random_between(10,12) ;
347 nextprime(p,p); //!@bug si p n'est pas premier, fgemm fait n'importe quoi (division par alpha)
348 Field Zp (p);
349 */
350
351 BlasMatrixDomain<Field> BMD (Zp);
352
353 // Ep = b C + a A B
354 ZMD.muladd(D,beta,C,alpha,A,B);
355 BlasMatrix<Field> Dp(D,Zp); // D mod p
356
357 BlasMatrix<Field> Ap(A,Zp);
358 BlasMatrix<Field> Bp(B,Zp);
359 BlasMatrix<Field> Cp(C,Zp);
360 // BlasMatrix<Field> Ap(A.rowdim(),A.coldim());
361 // BlasMatrix<Field> Bp(B.rowdim(),B.coldim());
362 // BlasMatrix<Field> Cp(C.rowdim(),C.coldim());
363 // MatrixHom::map(Ap,A,Zp);
364 // MatrixHom::map(Bp,B,Zp);
365 // MatrixHom::map(Cp,C,Zp);
366 BlasMatrix<Field> Ep(Zp,M,N); // D mod p
367
368 Element ap, bp ;
369 Zp.init(ap,alpha);
370 Zp.init(bp,beta);
371
372 // Ep = bp Cp + ap Ap Bp mod p
373 BMD.muladd(Ep,bp,Cp,ap,Ap,Bp);
374
375 bool pass = BMD.areEqual(Ep,Dp);
376 if (!pass) {
377 #if 0 /* maple check on stdout */
378 report << "#########################################" << std::endl;
379 report << "p := " << p << ';' << std::endl;
380 report << "ap,bp := " << ap << ',' << bp << ';' << std::endl;
381 Ap.write(report << "Ap :=", Tag::FileFormat::Maple) << ";" << std::endl;
382 Bp.write(report << "Bp :=", Tag::FileFormat::Maple) << ";" << std::endl;
383 Cp.write(report << "Cp :=", Tag::FileFormat::Maple) << ";" << std::endl;
384 Dp.write(report << "Dp :=", Tag::FileFormat::Maple) << ";" << std::endl;
385 Ep.write(report << "Ep :=", Tag::FileFormat::Maple) << ";" << std::endl;
386 report << "alpha,beta := " << alpha << ',' << beta << ';' << std::endl;
387 A.write(report << "A :=",Tag::FileFormat::Maple) << ';' << std::endl;
388 B.write(report << "B :=",Tag::FileFormat::Maple) << ';' << std::endl;
389 C.write(report << "C :=",Tag::FileFormat::Maple) << ';' << std::endl;
390 D.write(report << "E :=",Tag::FileFormat::Maple) << ';' << std::endl;
391 report << "evalm(E-alpha*A.B-beta*C);" << std::endl;
392 report << "#########################################" << std::endl;
393 #endif
394 mycommentator().report() << " *** BMD ERROR (" << alpha << ',' << beta << ") *** " << std::endl;
395 }
396
397 // Ep = bp Cp + ap Ap Bp mod p
398 BMD.muladd(Ep,bp,Cp,ap,Ap,Bp);
399 bool all = BMD.areEqual(Ep,Dp);
400 if (!all) {
401 mycommentator().report() << " *** MD ERROR *** " << std::endl;
402 }
403
404 return pass&all ;
405
406 }
407
408 // tests MulAdd for various parameters alpha and beta.
409 // This test ignores Field F. It tests only integer matrices.
410 // It should not really be included in launch_tests<Field>.
411 // For now we'll call it when characteristic is zero.
412 template <class Field>
testMulAddAgain(const Field & F,size_t n,int iterations)413 static bool testMulAddAgain (const Field& F, size_t n, int iterations)
414 {
415
416 Givaro::ZRing<Integer> ZZ ;
417 typedef BlasMatrix<Givaro::ZRing<Integer> > IMatrix;
418
419 //Commentator mycommentator;
420 mycommentator().getMessageClass (INTERNAL_DESCRIPTION).setMaxDepth (3);
421 mycommentator().getMessageClass (INTERNAL_DESCRIPTION).setMaxDetailLevel (Commentator::LEVEL_NORMAL);
422 mycommentator().start (pretty("Testing muladd again"),"testMulAddAgain",(unsigned int)iterations);
423
424 bool ret = true;
425
426 size_t ll = 17 ;
427 Integer::seeding();
428 #if 0
429 size_t lA = 15 ;
430 size_t lB = 18 ;
431 size_t lC = 19 ;
432
433 Givaro::ZRing<Integer> ZZ ;
434
435 typedef RandomIntegerIterator<> IntRandIter ;
436 typedef RandomDenseMatrix<IntRandIter, Givaro::ZRing<Integer> > IntRand_t;
437
438 IntRandIter RA(ZZ,lA);
439 IntRandIter RB(ZZ,lB);
440 IntRandIter RC(ZZ,lC);
441
442 IntRand_t Arand (ZZ,RA);
443 IntRand_t Brand (ZZ,RB);
444 IntRand_t Crand (ZZ,RC);
445 #endif
446
447 for (int k=0;k<iterations; ++k) {
448
449 mycommentator().progress(k);
450 IMatrix A(ZZ, n,n),B(ZZ, n,n),C(ZZ, n,n),D(ZZ, n,n),E(ZZ, n,n);
451
452 Integer a , b;
453 #if 0
454 Arand.random<IMatrix>(A);
455 Brand.random<IMatrix>(B);
456 Crand.random<IMatrix>(C);
457 #endif
458 for (size_t i=0;i<n;++i)
459 for (size_t j=0;j<n;++j){
460 A.setEntry(i,j,Integer::random());
461 B.setEntry(i,j,Integer::random());
462 C.setEntry(i,j,Integer::random());
463 }
464
465 a = 1 ; b = 1 ;
466 if (!CheckMulAdd(F,a,A,B,b,C)) ret = false ;
467 a = 1 ; b = -1 ;
468 if (!CheckMulAdd(F,a,A,B,b,C)) ret = false ;
469 a = -1 ; b = 1 ;
470 if (!CheckMulAdd(F,a,A,B,b,C)) ret = false ;
471 a = -1 ; b = -1 ;
472 if (!CheckMulAdd(F,a,A,B,b,C)) ret = false ;
473 a = 0 ; b = 1 ;
474 if (!CheckMulAdd(F,a,A,B,b,C)) ret = false ;
475 a = 1 ; b = 0 ;
476 if (!CheckMulAdd(F,a,A,B,b,C)) ret = false ;
477 a = 0 ; b = -1 ;
478 if (!CheckMulAdd(F,a,A,B,b,C)) ret = false ;
479 a = -1 ; b = 0 ;
480 if (!CheckMulAdd(F,a,A,B,b,C)) ret = false ;
481 a = Integer::random<false>(ll) ; b = 1 ;
482 if (!CheckMulAdd(F,a,A,B,b,C)) ret = false ;
483 a =Integer::random<false>(ll) ; b = -1 ;
484 if (!CheckMulAdd(F,a,A,B,b,C)) ret = false ;
485 a = 1 ; b = Integer::random<false>(ll) ;
486 if (!CheckMulAdd(F,a,A,B,b,C)) ret = false ;
487 a = -1 ; b = Integer::random<false>(ll) ;
488 if (!CheckMulAdd(F,a,A,B,b,C)) ret = false ;
489 a = 0 ; b = Integer::random<false>(ll) ;
490 if (!CheckMulAdd(F,a,A,B,b,C)) ret = false ;
491 a = Integer::random<false>(ll) ; b = 0 ;
492 if (!CheckMulAdd(F,a,A,B,b,C)) ret = false ;
493 a = Integer::random<false>(ll) ; b = Integer::random<false>(ll) ;
494 if (!CheckMulAdd(F,a,A,B,b,C)) ret = false ;
495
496 }
497
498
499 mycommentator().stop(MSG_STATUS (ret), (const char *) 0, "testMulAddAgain");
500
501 return ret;
502 }
503
504 // tests MulAdd for various shapes and values of transposition.
505 template <class Field>
testMulAddShapeTrans(const Field & F,size_t m,size_t n,size_t k,int iterations)506 static bool testMulAddShapeTrans (const Field &F, size_t m, size_t n, size_t k, int iterations)
507 {
508 bool ret = true ;
509 mycommentator().getMessageClass (INTERNAL_DESCRIPTION).setMaxDepth (3);
510 mycommentator().getMessageClass (INTERNAL_DESCRIPTION).setMaxDetailLevel (Commentator::LEVEL_NORMAL);
511 mycommentator().start (pretty("Testing muladd for shapes and transposition"),"testMulAddShapeTrans",(unsigned long)iterations);
512
513
514 typedef typename Field::Element Element;
515 typedef BlasMatrix<Field> Matrix ;
516 typedef TransposedBlasMatrix<Matrix> TransposedMatrix ;
517 typedef typename Field::RandIter Randiter ;
518 Randiter R(F) ;
519 RandomDenseMatrix<Randiter,Field> RandMat(F,R);
520
521 BlasMatrixDomain<Field> BMD (F);
522
523 // input matrix
524 Matrix A(F, m,k);
525 Matrix B(F, k,n);
526 Matrix C(F, m,n);
527 // result matrix
528 Matrix D(F, m,n);
529 Matrix E(F, m,n);
530
531 // random A,B
532 RandMat.random(A);
533 RandMat.random(B);
534 RandMat.random(C);
535
536 // hard tranpose A,B
537 Matrix A1 (F, k,m) ;
538 A.transpose(A1) ;
539 Matrix B1 (F, n,k) ;
540 B.transpose(B1) ;
541 TransposedMatrix tA(A1); // t(tA)=A
542 TransposedMatrix tB(B1); // t(tB)=B
543
544 // random alpha, beta
545 Element alpha ;
546 do {R.random(alpha);} while (F.isZero(alpha)); // nonzerorandom
547
548 Element beta ;
549 R.random(beta);
550
551 // témoin.
552 BMD.muladd(D,beta,C,alpha,A,B);
553
554 // A,B
555 BMD.muladd(E,beta,C,alpha,A,B);
556 if (!BMD.areEqual(E,D)) {
557 ret = false ;
558 mycommentator().report() << " *** BMD ERROR (" << alpha << ',' << beta << ") (noTrans, noTrans) *** " << std::endl;
559 }
560
561 BMD.muladd(E,beta,C,alpha,A,tB);
562 if (!BMD.areEqual(E,D)) {
563 ret = false ;
564 mycommentator().report() << " *** BMD ERROR (" << alpha << ',' << beta << ") (noTrans, Trans) *** " << std::endl;
565 }
566
567 BMD.muladd(E,beta,C,alpha,tA,B);
568 if (!BMD.areEqual(E,D)) {
569 ret = false ;
570 mycommentator().report() << " *** BMD ERROR (" << alpha << ',' << beta << ") (Trans, noTrans) *** " << std::endl;
571 }
572
573 BMD.muladd(E,beta,C,alpha,tA,tB);
574 if (!BMD.areEqual(E,D)) {
575 ret = false ;
576 mycommentator().report() << " *** BMD ERROR (" << alpha << ',' << beta << ") (Trans, Trans) *** " << std::endl;
577 }
578
579
580 mycommentator().stop(MSG_STATUS (ret), (const char *) 0, "testMulAddShapeTrans");
581 return ret ;
582 }
583
584 // tests MulAdd for various shapes and values of transposition.
585 template<class Field, bool LeftSide, bool UnitDiag>
testTriangMulShapeTrans(const Field & F,size_t m,size_t n,int iterations)586 static bool testTriangMulShapeTrans (const Field &F, size_t m, size_t n, int iterations)
587 {
588 bool ret = true ;
589 mycommentator().getMessageClass (INTERNAL_DESCRIPTION).setMaxDepth (3);
590 mycommentator().getMessageClass (INTERNAL_DESCRIPTION).setMaxDetailLevel (Commentator::LEVEL_NORMAL);
591 mycommentator().start (pretty("Testing triangular matmul for shapes and transposition"),"testTriangMulShapeTrans",(unsigned int)iterations);
592
593
594 // typedef typename Field::Element Element;
595 typedef BlasMatrix<Field> Matrix ;
596 typedef TriangularBlasMatrix<Field> TriangularMatrix ;
597 // typedef TransposedBlasMatrix<Matrix> TransposedMatrix ;
598 typedef TransposedBlasMatrix<TriangularMatrix > TransposedTriangular ;
599 typedef typename Field::RandIter Randiter ;
600 Randiter R(F) ;
601 RandomDenseMatrix<Randiter,Field> RandMat(F,R);
602
603
604 BlasMatrixDomain<Field> BMD (F);
605
606 int k =(int) (LeftSide?m:n) ;
607 // input matrix
608 Matrix A(F, k,k); // A = L+U-I. Either L or U is unit.
609 Matrix B(F, m,n);
610 // result matrix
611 Matrix D(F, m,n);
612 Matrix E(F, m,n);
613
614 // random A,B
615 RandMat.random(A);
616 RandMat.random(B);
617
618 // hard tranpose A,B
619 Matrix A1 (F, k,k) ;
620 A.transpose(A1) ;
621
622
623
624 /* test (L+U-I) B+B = LB+UB */
625 if (LeftSide)
626 BMD.muladd(D,F.one,B,F.one,A,B);
627 else
628 BMD.muladd(D,F.one,B,F.one,B,A);
629
630 /**** DIRECT ****/
631 {
632 /* L */
633 TriangularMatrix L (A, Tag::Shape::Lower,
634 (UnitDiag?Tag::Diag::Unit:Tag::Diag::NonUnit));
635
636 /* U */
637 TriangularMatrix U (A, Tag::Shape::Upper,
638 (!UnitDiag?Tag::Diag::Unit:Tag::Diag::NonUnit));
639
640 /* make product */
641 E = B ;
642 // Matrix G(m,n);
643 // G = E ;
644 Matrix G((const Matrix&)E); //!@warning on n'oublie pas l'esperluette !!!
645 if(LeftSide) {
646 BMD.mulin_right(L,E) ; // B <- AB
647 BMD.mulin_right(U,G) ;
648 }
649 else {
650 BMD.mulin_left(G,L) ; // B <- BA
651 BMD.mulin_left(E,U) ;
652 }
653 BMD.addin(E,G);
654
655 /* check equality */
656 if (!BMD.areEqual(E,D)) {
657 ret = false ;
658 mycommentator().report() << " *** BMD ERROR (" << (LeftSide?"left":"right") << ',' << (UnitDiag?" L":" U") << " is unit) *** " << std::endl;
659 }
660 else {
661 mycommentator().report() << " direct triangular multiplication ok." << std::endl;
662 }
663 }
664 /**** Transpose ****/
665 {
666 /* L */
667 TriangularMatrix L1 (A1, Tag::Shape::Lower,
668 (UnitDiag?Tag::Diag::Unit:Tag::Diag::NonUnit));
669
670 /* U */
671 TriangularMatrix U1 (A1, Tag::Shape::Upper,
672 (!UnitDiag?Tag::Diag::Unit:Tag::Diag::NonUnit));
673
674 TransposedTriangular L(L1);
675 TransposedTriangular U(U1);
676 /* make product */
677 E = B ;
678 // Matrix G(m,n);
679 // G = E ;
680 Matrix G((const Matrix&)E); //!@warning on n'oublie pas l'esperluette !!!
681 if(LeftSide) {
682 BMD.mulin_right(L,E) ; // B <- AB
683 BMD.mulin_right(U,G) ;
684 }
685 else {
686 BMD.mulin_left(G,L) ; // B <- BA
687 BMD.mulin_left(E,U) ;
688 }
689 BMD.addin(E,G);
690
691 /* check equality */
692 if (!BMD.areEqual(E,D)) {
693 ret = false ;
694 mycommentator().report() << " *** BMD ERROR Transpose (" << (LeftSide?"left":"right") << ',' << (UnitDiag?" L":" U") << " is unit) *** " << std::endl;
695 }
696 else {
697 mycommentator().report() << " transposed triangular multiplication ok." << std::endl;
698 }
699 }
700
701 mycommentator().stop(MSG_STATUS (ret), (const char *) 0, "testMulAddShapeTrans");
702 return ret ;
703 }
704
705
706 /*
707 * Testing the rank of dense matrices using BlasDomain
708 * construct a n*n matrices of rank r and compute the rank
709 */
710 template <class Field>
testRank(const Field & F,size_t n,int iterations)711 static bool testRank (const Field& F,size_t n, int iterations)
712 {
713
714 typedef typename Field::Element Element;
715 typedef typename Field::RandIter RandIter;
716
717 //Commentator mycommentator;
718 mycommentator().getMessageClass (INTERNAL_DESCRIPTION).setMaxDepth (3);
719 mycommentator().getMessageClass (INTERNAL_DESCRIPTION).setMaxDetailLevel (Commentator::LEVEL_NORMAL);
720 mycommentator().start (pretty("Testing rank"),"testRank",(unsigned int)iterations);
721
722 RandIter G(F);
723 Givaro::GeneralRingNonZeroRandIter<Field> Gn(G);
724 Element tmp;
725 bool ret = true;
726 BlasMatrixDomain<Field> BMD(F);
727
728 for (int k=0;k<iterations; ++k) {
729 unsigned int r;
730
731 mycommentator().progress(k);
732 BlasMatrix<Field> A(F,n,n),S(F,n,n), L(F,n,n);
733
734 r = (unsigned int)((size_t)random() % n);
735 // create S as an upper triangular matrix with r nonzero rows
736 for (size_t i=0;i<r;++i){
737 S.setEntry(i,i,Gn.random(tmp));
738 for (size_t j=i+1;j<n;++j)
739 S.setEntry(i,j,G.random(tmp));
740 }
741 BMD.write(commentator().report(), S) << std::endl;
742
743
744 // create L as a lower triangular matrix with nonzero elements on the diagonal
745 for (size_t i=0;i<n;++i){
746 for (size_t j=0;j<i;++j)
747 L.setEntry(i,j,G.random(tmp));
748 L.setEntry(i,i,Gn.random(tmp));
749 }
750 BMD.write(commentator().report(), L) << std::endl;
751
752 // compute A=LS
753 BMD.mul(A,L,S);
754 BMD.write(commentator().report(), A) << std::endl;
755
756 // compute the rank of A
757 unsigned int rank= BMD.rankInPlace(A);
758 commentator().report() << "Rank " << rank << " should be " << r << std::endl;
759
760 if (rank!=r)
761 ret=false;
762 }
763
764 mycommentator().stop(MSG_STATUS (ret), (const char *) 0, "testRank");
765
766 return ret;
767 }
768
769
770 /*
771 * Testing the determinant of dense matrices using BlasDomain
772 * construct a n*n matrices of determinant d and compute the determinant
773 */
774 template <class Field>
testDet(const Field & F,size_t n,int iterations)775 static bool testDet (const Field& F,size_t n, int iterations)
776 {
777
778 typedef typename Field::Element Element;
779 typedef typename Field::RandIter RandIter;
780
781 //Commentator mycommentator;
782 mycommentator().getMessageClass (INTERNAL_DESCRIPTION).setMaxDepth (3);
783 mycommentator().getMessageClass (INTERNAL_DESCRIPTION).setMaxDetailLevel (Commentator::LEVEL_NORMAL);
784 mycommentator().start (pretty("Testing determinant"),"testDet",(unsigned int)iterations);
785
786 RandIter G(F);
787 Givaro::GeneralRingNonZeroRandIter<Field> Gn(G);
788 Element tmp,d;
789
790 bool ret = true;
791 BlasMatrixDomain<Field> BMD(F);
792
793 for (int k=0;k<iterations;++k) {
794
795 mycommentator().progress(k);
796
797 G.random(d);
798
799 BlasMatrix<Field> A(F,n,n),S(F,n,n), L(F,n,n);
800
801 // create S as an upper triangular matrix of full rank
802 // with diagonal's element equal to 1 except the first entry wich equals to d
803 for (size_t i=0;i<n;++i){
804 S.setEntry(i,i,F.one);
805 for (size_t j=i+1;j<n;++j)
806 S.setEntry(i,j,G.random(tmp));
807 }
808 S.setEntry(0,0,d);
809
810 // create L as a lower triangular matrix with only 1's on diagonal
811 for (size_t i=0;i<n;++i){
812 for (size_t j=0;j<i;++j)
813 L.setEntry(i,j,G.random(tmp));
814 L.setEntry(i,i,F.one);
815 }
816
817
818 // compute A=LS
819 BMD.mul(A,L,S);
820
821 // compute the determinant of A
822 Element det= BMD.detInPlace(A);
823
824 if (!F.areEqual(det,d))
825 ret=false;
826 }
827
828 mycommentator().stop(MSG_STATUS (ret), (const char *) 0, "testDet");
829
830 return ret;
831 }
832
833 /*
834 * Testing the inverse of dense matrices using BlasDomain
835 * construct a non-singular n*n matrices
836 */
837 template <class Field>
testInv(const Field & F,size_t n,int iterations)838 static bool testInv (const Field& F,size_t n, int iterations)
839 {
840
841 typedef typename Field::Element Element;
842 typedef typename Field::RandIter RandIter;
843 typedef BlasMatrix<Field> Matrix;
844
845 //Commentator mycommentator;
846 mycommentator().getMessageClass (INTERNAL_DESCRIPTION).setMaxDepth (3);
847 mycommentator().getMessageClass (INTERNAL_DESCRIPTION).setMaxDetailLevel (Commentator::LEVEL_NORMAL);
848 mycommentator().start (pretty("Testing inverse"),"testInv",(unsigned int)iterations);
849
850 RandIter G(F);
851 Givaro::GeneralRingNonZeroRandIter<Field> Gn(G);
852 Element tmp;
853
854 bool ret = true;
855 BlasMatrixDomain<Field> BMD(F);
856
857 Matrix Id(F, n,n);
858 for (size_t i=0;i<n;++i)
859 Id.setEntry(i,i,F.one);
860
861 if (n < 10)
862 Id.write(mycommentator().report() << "Id" << std::endl) << std::endl;
863 for (int k=0;k<iterations;++k) {
864
865 mycommentator().progress(k);
866
867
868 Matrix A(F, n,n),S(F, n,n), L(F, n,n); // , invA(F, n,n);
869
870 // create S as an upper triangular matrix of full rank
871 // with nonzero random diagonal's element
872 for (size_t i=0;i<n;++i){
873 S.setEntry(i,i,Gn.random(tmp));
874 for (size_t j=i+1;j<n;++j)
875 S.setEntry(i,j,G.random(tmp));
876 }
877
878 // create L as a lower triangular matrix
879 // with only 1's on diagonal
880 for (size_t i=0;i<n;++i){
881 for (size_t j=0;j<i;++j)
882 L.setEntry(i,j,G.random(tmp));
883 L.setEntry(i,i,F.one);
884 }
885
886 // compute A=LS
887 BMD.mul(A,L,S);
888 //for (size_t i=0;i<n;++i){A.setEntry(i,i,F.one); }
889 //A.setEntry(0, n-1, F.mOne);
890
891 if (n < 10)
892 A.write(mycommentator().report() << "A") << std::endl;
893
894 // compute the inverse of A
895 Matrix invA(A);
896 if (n < 10)
897 invA.write(mycommentator().report() << "before inversion, invA") << std::endl;
898 //int nullity;
899 //FFPACK::Invert2 (F, invA.rowdim(), A.getPointer(), A.getStride(), invA.getPointer(), invA.getStride(), nullity);
900 BMD.invin(invA);
901 if (n < 10)
902 invA.write(mycommentator().report() << "invA") << std::endl;
903
904 // compute Ainv*A and A*Ainv
905 BMD.mul(L,invA,A);
906 if (n < 10)
907 L.write(mycommentator().report() << "invA*A") << std::endl;
908 BMD.mul(S,A,invA);
909 if (n < 10)
910 S.write(mycommentator().report() << "A*invA") << std::endl;
911
912 if (!BMD.areEqual(L,Id) || !BMD.areEqual(S,Id))
913 ret=false;
914 }
915
916 mycommentator().stop(MSG_STATUS (ret), (const char *) 0, "testInv");
917
918 return ret;
919 }
920
921
922 /*
923 * Test resolution of linear system with a triangular matrix
924 */
925 template <class Field>
testTriangularSolve(const Field & F,size_t m,size_t n,int iterations)926 static bool testTriangularSolve (const Field& F, size_t m, size_t n, int iterations)
927 {
928
929 typedef typename Field::Element Element;
930 typedef BlasMatrix<Field> Matrix;
931 typedef TriangularBlasMatrix<Field> TriangularMatrix;
932 typedef typename Field::RandIter RandIter;
933
934 //Commentator mycommentator;
935 mycommentator().getMessageClass (INTERNAL_DESCRIPTION).setMaxDepth (3);
936 mycommentator().getMessageClass (INTERNAL_DESCRIPTION).setMaxDetailLevel (Commentator::LEVEL_NORMAL);
937 mycommentator().start (pretty("Testing triangular solver"),"testTriangularSolve",(unsigned int)iterations);
938
939 RandIter G(F);
940 Givaro::GeneralRingNonZeroRandIter<Field> Gn(G);
941 Element tmp;
942
943 bool ret = true;
944 BlasMatrixDomain<Field> BMD(F);
945
946 for (int k=0;k<iterations;++k) {
947
948 mycommentator().progress(k);
949
950 Matrix Al(F, m,m),Au(F, m,m);
951 Matrix X(F, m,n), B(F, m,n), C(F, m,n);
952
953 BlasVector<Field> b(F,m),x(F,m),c(F,m);
954
955 // Create B a random matrix
956 for (size_t i=0;i<m;++i)
957 for (size_t j=0;j<m;++j)
958 B.setEntry(i,j,G.random(tmp));
959
960 //create random vector b
961 for( size_t i=0;i<m;++i)
962 F.init(b[i],G.random(tmp));
963
964 // Create Au a random full rank upper triangular matrix
965 for (size_t i=0;i<m;++i){
966 Au.setEntry(i,i,Gn.random(tmp));
967 for (size_t j=i+1;j<m;++j)
968 Au.setEntry(i,j,G.random(tmp));
969 }
970
971 // Create Al a random full rank lower triangular matrix
972 for (size_t i=0;i<m;++i){
973 for (size_t j=0;j<i;++j)
974 Al.setEntry(i,j,G.random(tmp));
975 Al.setEntry(i,i,Gn.random(tmp));
976 }
977
978 // Create 2 trinagular matrix as view of matrix
979 TriangularMatrix TAl(Al,Tag::Shape::Lower,Tag::Diag::NonUnit), TAu(Au,Tag::Shape::Upper,Tag::Diag::NonUnit);
980
981 // testing solver with matrix right hand side
982 BMD.left_solve(X,TAl,B);
983 BMD.mul(C,Al,X);
984 if (!BMD.areEqual(C,B))
985 ret=false;
986
987 BMD.left_solve(X,TAu,B);
988 BMD.mul(C,Au,X);
989 if (!BMD.areEqual(C,B))
990 ret=false;
991
992 // testing solver with matrix left hand side
993 BMD.right_solve(X,TAl,B);
994 BMD.mul(C,X,Al);
995 if (!BMD.areEqual(C,B))
996 ret=false;
997
998 BMD.right_solve(X,TAu,B);
999 BMD.mul(C,X,Au);
1000 if (!BMD.areEqual(C,B))
1001 ret=false;
1002
1003
1004 // testing solver with vector right hand side
1005 BMD.left_solve(x,TAl,b);
1006 BMD.mul(c,Al,x);
1007 if (!localAreEqual(c,b))
1008 ret=false;
1009
1010 BMD.left_solve(x,TAu,b);
1011 BMD.mul(c,Au,x);
1012 if (!localAreEqual(c,b))
1013 ret=false;
1014
1015 // testing solver with vector left hand side
1016 BMD.right_solve(x,TAl,b);
1017 BMD.mul(c,x,Al);
1018 if (!localAreEqual(c,b))
1019 ret=false;
1020
1021 BMD.right_solve(x,TAu,b);
1022 BMD.mul(c,x,Au);
1023 if (!localAreEqual(c,b))
1024 ret=false;
1025 }
1026
1027 mycommentator().stop(MSG_STATUS (ret), (const char *) 0, "testTriangularSolve");
1028
1029 return ret;
1030 }
1031
1032 /*
1033 * Test resolution of linear system with a matrix
1034 */
1035 template <class Field>
testSolve(const Field & F,size_t m,size_t n,int iterations)1036 static bool testSolve (const Field& F, size_t m, size_t n, int iterations)
1037 {
1038
1039 typedef typename Field::Element Element;
1040 typedef BlasMatrix<Field> Matrix;
1041 typedef typename Field::RandIter RandIter;
1042
1043 //Commentator mycommentator;
1044 mycommentator().getMessageClass (INTERNAL_DESCRIPTION).setMaxDepth (3);
1045 mycommentator().getMessageClass (INTERNAL_DESCRIPTION).setMaxDetailLevel (Commentator::LEVEL_NORMAL);
1046 mycommentator().start (pretty("Testing solver"),"testTriangularSolve",(unsigned int)iterations);
1047
1048 RandIter G(F);
1049 Givaro::GeneralRingNonZeroRandIter<Field> Gn(G);
1050 Element tmp;
1051
1052 bool ret = true;
1053 BlasMatrixDomain<Field> BMD(F);
1054
1055 for (int k=0;k<iterations;++k) {
1056
1057 mycommentator().progress(k);
1058
1059 Matrix A(F, m,m),L(F, m,m),S(F, m,m);
1060 Matrix X(F, m,n), B(F, m,n), C(F, m,n);
1061
1062 BlasVector<Field> b(F,m),x(F,m),c(F,m);
1063
1064 // Create B a random matrix
1065 for (size_t i=0;i<m;++i)
1066 for (size_t j=0;j<n;++j)
1067 B.setEntry(i,j,G.random(tmp));
1068
1069 //create a random vector b
1070 for( size_t i=0;i<m;++i)
1071 F.init(b[i],G.random(tmp));
1072
1073
1074 // create S as an upper triangular matrix of full rank
1075 // with nonzero random diagonal's element
1076 for (size_t i=0;i<m;++i){
1077 S.setEntry(i,i,Gn.random(tmp));
1078 for (size_t j=i+1;j<n;++j)
1079 S.setEntry(i,j,G.random(tmp));
1080 }
1081
1082 // create L as a lower triangular matrix
1083 // with only 1's on diagonal
1084 for (size_t i=0;i<m;++i){
1085 for (size_t j=0;j<i;++j)
1086 L.setEntry(i,j,G.random(tmp));
1087 L.setEntry(i,i,F.one);
1088 }
1089
1090 // compute A=LS
1091 BMD.mul(A,L,S);
1092
1093
1094 // testing solver with matrix right hand side
1095 BMD.left_solve(X,A,B);
1096 BMD.mul(C,A,X);
1097 if (!BMD.areEqual(C,B))
1098 ret=false;
1099
1100 BMD.left_solve(X,A,B);
1101 BMD.mul(C,A,X);
1102 if (!BMD.areEqual(C,B))
1103 ret=false;
1104
1105 // testing solver with matrix left hand side
1106 BMD.right_solve(X,A,B);
1107 BMD.mul(C,X,A);
1108 if (!BMD.areEqual(C,B))
1109 ret=false;
1110
1111 BMD.right_solve(X,A,B);
1112 BMD.mul(C,X,A);
1113 if (!BMD.areEqual(C,B))
1114 ret=false;
1115
1116
1117 // testing solver with vector right hand side
1118 BMD.left_solve(x,A,b);
1119 BMD.mul(c,A,x);
1120 if (!localAreEqual(c,b))
1121 ret=false;
1122
1123 BMD.left_solve(x,A,b);
1124 BMD.mul(c,A,x);
1125 if (!localAreEqual(c,b))
1126 ret=false;
1127
1128 // testing solver with vector left hand side
1129 BMD.right_solve(x,A,b);
1130 BMD.mul(c,x,A);
1131 if (!localAreEqual(c,b))
1132 ret=false;
1133
1134 BMD.right_solve(x,A,b);
1135 BMD.mul(c,x,A);
1136 if (!localAreEqual(c,b))
1137 ret=false;
1138
1139 }
1140
1141 mycommentator().stop(MSG_STATUS (ret), (const char *) 0, "testTriangularSolve");
1142
1143 return ret;
1144 }
1145
1146 /*
1147 * Test of the BlasPermutations
1148 */
1149 template <class Field>
testPermutation(const Field & F,size_t m,int iterations)1150 static bool testPermutation (const Field& F, size_t m, int iterations)
1151 {
1152
1153 typedef typename Field::Element Element;
1154 typedef BlasMatrix<Field> Matrix;
1155 typedef typename Field::RandIter RandIter;
1156
1157 //Commentator mycommentator;
1158 mycommentator().getMessageClass (INTERNAL_DESCRIPTION).setMaxDepth (3);
1159 mycommentator().getMessageClass (INTERNAL_DESCRIPTION).setMaxDetailLevel (Commentator::LEVEL_NORMAL);
1160 mycommentator().start (pretty("Testing permutations"),"testPermutation",(unsigned int)iterations);
1161
1162 RandIter G(F);
1163 Givaro::GeneralRingNonZeroRandIter<Field> Gn(G);
1164 Element tmp;
1165
1166 bool ret = true;
1167 BlasMatrixDomain<Field> BMD(F);
1168
1169 for (int k=0;k<iterations;++k) {
1170
1171 mycommentator().progress(k);
1172
1173 // std::vector<size_t> P(m);
1174
1175 // Field Z2(2);
1176 // RandIter G2(Z2);
1177
1178 // for (size_t i=0; i<m; ++i){
1179 // G.random(tmp);
1180 // if ( Z2.isZero(G2.random(tmp2) ) )
1181 // P[i] = i + ( (size_t) random() % (m-i) );
1182 // else
1183 // P[i] = i;
1184 // }
1185
1186 //report<<P<<std::endl;
1187 Matrix A(F, m,m), Abis(F, m,m), B(F, m,m), C(F, m,m), D(F, m,m);
1188 BlasVector<Field> a(F,m),abis(F,m),b(F,m),c(F,m), d(F,m);
1189 BlasPermutation<size_t> Perm(m);
1190 RandomBlasPermutation(Perm);
1191
1192 // Create A a random matrix
1193 for (size_t i=0;i<m;++i)
1194 for (size_t j=0;j<m;++j)
1195 A.setEntry(i,j,Gn.random(tmp));
1196 // Create a a random vector
1197 for (size_t i=0;i<m;++i)
1198 F.assign(a[i],Gn.random(tmp));
1199
1200 /*
1201 * Test A.P.P^t == A
1202 */
1203
1204 // B = A.P
1205 BMD.mul( B, A, Perm);
1206 // C = B.P^t
1207 BMD.mul( C, B, TransposedBlasMatrix<BlasPermutation<size_t> >(Perm) );
1208 // Test C==A
1209 if (!BMD.areEqual(A,C))
1210 ret=false;
1211 /*
1212 * Test A.P^t.P == A
1213 */
1214
1215 // B = A.P^t
1216 BMD.mul( B, A, TransposedBlasMatrix<BlasPermutation<size_t> >(Perm));
1217 // C = B.P
1218 BMD.mul( C, B, Perm );
1219 // Test C==A
1220 if (!BMD.areEqual(A,C))
1221 ret=false;
1222 /*
1223 * Test P.P^t.A == A
1224 */
1225
1226 // B = P.A
1227 BMD.mul( B, Perm, A);
1228 // C = P^t.B
1229 BMD.mul( C, TransposedBlasMatrix<BlasPermutation<size_t> >(Perm) , B);
1230 // Test C==A
1231 if (!BMD.areEqual(A,C))
1232 ret=false;
1233 /*
1234 * Test P^t.P.A == A
1235 */
1236
1237 // B = P^t.A
1238 BMD.mul( B, TransposedBlasMatrix<BlasPermutation<size_t> >(Perm), A);
1239 // C = P.B
1240 BMD.mul( C, Perm, B);
1241 // Test C==A
1242 if (!BMD.areEqual(A,C))
1243 ret=false;
1244
1245 /*
1246 * Test a.P.P^t == a
1247 */
1248
1249 // b = a.P
1250 BMD.mul( b, a, Perm);
1251 // c = b.P^t
1252 BMD.mul( c, b, TransposedBlasMatrix<BlasPermutation<size_t> >(Perm) );
1253 // Test c==a
1254 if (!localAreEqual(a,c))
1255 ret=false;
1256
1257 /*
1258 * Test a.P^t.P == a
1259 */
1260
1261 // b = a.P^t
1262 BMD.mul( b, a, TransposedBlasMatrix<BlasPermutation<size_t> >(Perm));
1263 // c = B.P
1264 BMD.mul( c, b, Perm );
1265 // Test c==a
1266 if (!localAreEqual(a,c))
1267 ret=false;
1268 /*
1269 * Test P.P^t.a == a
1270 */
1271
1272 // b = P.a
1273 BMD.mul( b, Perm, a);
1274 // c = P^t.b
1275 BMD.mul( c, TransposedBlasMatrix<BlasPermutation<size_t> >(Perm) , b);
1276 // Test c==a
1277 if (!localAreEqual(a,c))
1278 ret=false;
1279
1280 /*
1281 * Test P^t.P.a == a
1282 */
1283
1284 // b = P^t.a
1285 BMD.mul( b, TransposedBlasMatrix<BlasPermutation<size_t> >(Perm), a);
1286 // c = P.b
1287 BMD.mul( c, Perm, b);
1288 // Test c==a
1289 if (!localAreEqual(a,c))
1290 ret=false;
1291
1292
1293 /*
1294 * Test P^t.A.(P.A)^-1.B == B
1295 */
1296 // Create B a random matrix
1297 for (size_t i=0;i<m;++i)
1298 for (size_t j=0;j<m;++j)
1299 B.setEntry(i,j,G.random(tmp));
1300
1301 // Abis = P.A
1302 BMD.mul( Abis, Perm, A);
1303 // C = (P.A)^-1.B
1304 BMD.left_solve( C, Abis, B);
1305 // D = A.C (= P^-1.B)
1306 BMD.mul(D, A, C);
1307 // D = P.D
1308 BMD.mulin_right( Perm,D);
1309 if (!BMD.areEqual(D,B))
1310 ret=false;
1311 /*
1312 * Test A.P^t.(A.P)^-1.B == B
1313 */
1314 // Create B a random matrix
1315 for (size_t i=0;i<m;++i)
1316 for (size_t j=0;j<m;++j)
1317 B.setEntry(i,j,G.random(tmp));
1318
1319 // Abis = A.P
1320 BMD.mul( Abis, A, Perm);
1321 // C = (A.P)^-1.B
1322 BMD.left_solve( C, Abis, B);
1323 // C = P.C
1324 BMD.mulin_right( Perm,C);
1325 // D = A.C (= P^-1.B)
1326 BMD.mul(D, A, C);
1327
1328 if (!BMD.areEqual(D,B))
1329 ret=false;
1330 /*
1331 * Test B.P^t.A.(P.A)^-1 == B
1332 */
1333 // Create B a random matrix
1334 for (size_t i=0;i<m;++i)
1335 for (size_t j=0;j<m;++j)
1336 B.setEntry(i,j,G.random(tmp));
1337
1338 // Abis = P.A
1339 BMD.mul( Abis, Perm, A);
1340 // C = B.(P.A)^-1
1341 BMD.right_solve( C, Abis, B);
1342 // C = C.P
1343 BMD.mulin_left( C,Perm);
1344 // D = C.A (=B)
1345 BMD.mul(D, C, A);
1346 if (!BMD.areEqual(D,B))
1347 ret=false;
1348
1349 /*
1350 * Test B.A.P^t.(A.P)^-1 == B
1351 */
1352 // Create B a random matrix
1353 for (size_t i=0;i<m;++i)
1354 for (size_t j=0;j<m;++j)
1355 B.setEntry(i,j,G.random(tmp));
1356
1357 // Abis = A.P
1358 BMD.mul( Abis, A, Perm);
1359 // C = B.(A.P)^-1
1360 BMD.right_solve( C, Abis, B);
1361 // D = C.A (= B.P^t)
1362 BMD.mul(D, C, A);
1363 // C = C.P
1364 BMD.mulin_left( D, Perm);
1365
1366 if (!BMD.areEqual(D,B))
1367 ret=false;
1368 /*
1369 * Test P.A.(P^t.A)^-1.B == B
1370 */
1371 // Create B a random matrix
1372 for (size_t i=0;i<m;++i)
1373 for (size_t j=0;j<m;++j)
1374 B.setEntry(i,j,G.random(tmp));
1375
1376 // Abis = P^t.A
1377 BMD.mul( Abis, TransposedBlasMatrix<BlasPermutation<size_t> >(Perm), A);
1378 // C = (P^t.A)^-1.B
1379 BMD.left_solve( C, Abis, B);
1380 // D = A.C (= P.B)
1381 BMD.mul(D, A, C);
1382 // D = P^t.D
1383 BMD.mulin_right( TransposedBlasMatrix<BlasPermutation<size_t> >(Perm),D);
1384 if (!BMD.areEqual(D,B))
1385 ret=false;
1386 /*
1387 * Test A.P.(A.P^t)^-1.B == B
1388 */
1389 // Create B a random matrix
1390 for (size_t i=0;i<m;++i)
1391 for (size_t j=0;j<m;++j)
1392 B.setEntry(i,j,G.random(tmp));
1393
1394 // Abis = A.P^t
1395 BMD.mul( Abis, A, TransposedBlasMatrix<BlasPermutation<size_t> >(Perm));
1396 // C = (A.P^t)^-1.B
1397 BMD.left_solve( C, Abis, B);
1398 // C = P^t.C
1399 BMD.mulin_right( TransposedBlasMatrix<BlasPermutation<size_t> >(Perm),C);
1400 // D = A.C (= P.B)
1401 BMD.mul(D, A, C);
1402
1403 if (!BMD.areEqual(D,B))
1404 ret=false;
1405 /*
1406 * Test B.P.A.(P^t.A)^-1 == B
1407 */
1408 // Create B a random matrix
1409 for (size_t i=0;i<m;++i)
1410 for (size_t j=0;j<m;++j)
1411 B.setEntry(i,j,G.random(tmp));
1412
1413 // Abis = P^t.A
1414 BMD.mul( Abis, TransposedBlasMatrix<BlasPermutation<size_t> >(Perm), A);
1415 // C = B.(P^t.A)^-1
1416 BMD.right_solve( C, Abis, B);
1417 // C = C.P^t
1418 BMD.mulin_left( C,TransposedBlasMatrix<BlasPermutation<size_t> >(Perm));
1419 // D = C.A (=B)
1420 BMD.mul(D, C, A);
1421 if (!BMD.areEqual(D,B))
1422 ret=false;
1423
1424 /*
1425 * Test B.A.P.(A.P^t)^-1 == B
1426 */
1427 // Create B a random matrix
1428 for (size_t i=0;i<m;++i)
1429 for (size_t j=0;j<m;++j)
1430 B.setEntry(i,j,G.random(tmp));
1431
1432 // Abis = A.P^t
1433 BMD.mul( Abis, A, TransposedBlasMatrix<BlasPermutation<size_t> >(Perm));
1434 // C = B.(A.P^t)^-1
1435 BMD.right_solve( C, Abis, B);
1436 // D = C.A (= B.P)
1437 BMD.mul(D, C, A);
1438 // C = C.P^t
1439 BMD.mulin_left( D, TransposedBlasMatrix<BlasPermutation<size_t> >(Perm));
1440
1441 if (!BMD.areEqual(D,B))
1442 ret=false;
1443 }
1444 mycommentator().stop(MSG_STATUS (ret), (const char *) 0, "testPLUQ");
1445
1446 return ret;
1447 }
1448
1449 /*
1450 * Test of the PLUQMatrix class
1451 */
1452 template <class Field>
testPLUQ(const Field & F,size_t m,size_t n,int iterations)1453 static bool testPLUQ (const Field& F, size_t m, size_t n, int iterations)
1454 {
1455
1456 typedef typename Field::Element Element;
1457 typedef BlasMatrix<Field> Matrix;
1458 typedef typename Field::RandIter RandIter;
1459
1460 //Commentator mycommentator;
1461 mycommentator().getMessageClass (INTERNAL_DESCRIPTION).setMaxDepth (3);
1462 mycommentator().getMessageClass (INTERNAL_DESCRIPTION).setMaxDetailLevel (Commentator::LEVEL_NORMAL);
1463 mycommentator().start (pretty("Testing PLUQ factorization"),"testPLUQ",(unsigned int)iterations);
1464
1465 RandIter G(F);
1466 Givaro::GeneralRingNonZeroRandIter<Field> Gn(G);
1467 Element tmp;
1468
1469 bool ret = true;
1470 BlasMatrixDomain<Field> BMD(F);
1471
1472 for (int k=0;k<iterations;++k) {
1473
1474 mycommentator().progress(k);
1475
1476 Matrix A(F, m,n), Abis(F, m,n), B(F, m,m), C(F, m,n);
1477
1478
1479 // Create B a random matrix of rank n/2
1480 for (size_t j=0;j<m;++j)
1481 if ( j % 2 )
1482 for (size_t i=0;i<m;++i)
1483 B.setEntry(i,j,G.random(tmp));
1484 else
1485 for (size_t i=0;i<m;++i)
1486 B.setEntry(i,j,F.zero);
1487 // Create C a random matrix of rank n/2
1488 for (size_t i=0;i<m;++i)
1489 if ( i % 2 )
1490 for (size_t j=0;j<n;++j)
1491 C.setEntry(i,j,G.random(tmp));
1492 else
1493 for (size_t j=0;j<n;++j)
1494 C.setEntry(i,j,F.zero);
1495
1496 // A = B*C
1497 BMD.mul(A, B, C);
1498
1499 Abis = A;
1500
1501 BlasPermutation<size_t> P(A.rowdim()),Q(A.coldim());
1502 PLUQMatrix<Field> X(A,P,Q);
1503
1504 TriangularBlasMatrix<Field> L(F,m,m,Tag::Shape::Lower,Tag::Diag::Unit);
1505 TriangularBlasMatrix<Field> U(F,m,n,Tag::Shape::Upper,Tag::Diag::NonUnit);
1506 X.getL(L);
1507 X.getU(U);
1508 P=X.getP();
1509
1510 Q=X.getQ();
1511
1512 // C = U*Q
1513 BMD.mul( C, U, Q);
1514 // A = L*C
1515 BMD.mul( A, L, C);
1516 // C = P*C
1517 BMD.mulin_right( P, C);
1518
1519 if (!BMD.areEqual(A,Abis))
1520 ret=false;
1521
1522 // Second pass
1523 // A = B*C
1524 BMD.mul(A, B, C);
1525
1526 Abis = A;
1527
1528 PLUQMatrix<Field> Y(A,P,Q);
1529
1530 TriangularBlasMatrix<Field> L2(F,m,m,Tag::Shape::Lower,Tag::Diag::Unit);
1531 TriangularBlasMatrix<Field> U2(F,m,n,Tag::Shape::Upper,Tag::Diag::NonUnit);
1532 Y.getL(L2);
1533 Y.getU(U2);
1534 P=Y.getP();
1535
1536 Q=Y.getQ();
1537
1538 // C = L2*U2
1539 BMD.mul( C,L2,U2);
1540 // C = C*Q
1541 BMD.mulin_left( C,Q);
1542 // A = P*C
1543 BMD.mul( A, P, C);
1544
1545 if (!BMD.areEqual(A,Abis))
1546 ret=false;
1547 }
1548
1549 mycommentator().stop(MSG_STATUS (ret), (const char *) 0, "testPLUQ");
1550
1551 return ret;
1552 }
1553
1554 template <class Field>
testMinPoly(const Field & F,size_t n,int iterations)1555 static bool testMinPoly (const Field& F, size_t n, int iterations)
1556 {
1557 typedef typename Field::Element Element;
1558 typedef BlasMatrix<Field> Matrix;
1559 typedef typename Field::RandIter RandIter;
1560 typedef BlasVector<Field> Polynomial;
1561 //Commentator mycommentator;
1562 mycommentator().getMessageClass (INTERNAL_DESCRIPTION).setMaxDepth (3);
1563 mycommentator().getMessageClass (INTERNAL_DESCRIPTION).setMaxDetailLevel (Commentator::LEVEL_NORMAL);
1564 mycommentator().start (pretty("Testing minpoly"),"testMinPoly",(unsigned int)iterations);
1565 Element tmp ;
1566 RandIter G(F);
1567 Givaro::GeneralRingNonZeroRandIter<Field> Gn(G);
1568 bool ret = true;
1569 BlasMatrixDomain<Field> BMD(F);
1570
1571 for (int k=0;k<iterations;++k) {
1572
1573 mycommentator().progress(k);
1574
1575 Matrix A(F,n,n);
1576 Polynomial P(F);
1577 // Test MinPoly(In) = X-1
1578 for (size_t i=0;i<n;++i)
1579 A.setEntry(i,i,F.one);
1580
1581 BMD.minpoly( P, A );
1582
1583 if ( P.size() !=2 )
1584 ret = false;
1585 if ( !F.areEqual(P[0], F.mOne) )
1586 ret = false;
1587 if ( !F.areEqual(P[1], F.one) )
1588 ret = false;
1589
1590 // Test MinPoly(a*In) = X-a
1591 G.random(tmp);
1592
1593 for (size_t i=0;i<n;++i)
1594 A.setEntry(i,i,tmp);
1595 F.negin(tmp);
1596 BMD.minpoly( P, A );
1597 if ( P.size() !=2 )
1598 ret = false;
1599 if ( !F.areEqual(P[0], tmp) )
1600 ret = false;
1601 if ( !F.areEqual(P[1], F.one) )
1602 ret = false;
1603
1604 for (size_t i=0;i<n-1;++i){
1605 for (size_t j=0; j<i+1; ++j)
1606 A.setEntry(i,j,F.zero);
1607 A.setEntry(i,i+1,F.one);
1608 if (i<n-2)
1609 for (size_t j=i+2; j<n; ++j)
1610 A.setEntry(i,j,F.zero);
1611 }
1612 for (size_t j=0;j<n;++j)
1613 A.setEntry(n-1,j,F.zero);
1614
1615 BMD.minpoly( P, A );
1616 if ( P.size() !=n+1 )
1617 ret = false;
1618 for (size_t i=0; i<n;++i)
1619 if ( !F.areEqual(P[i], F.zero) )
1620 ret = false;
1621 if ( !F.areEqual(P[n], F.one) )
1622 ret = false;
1623
1624 }
1625
1626 mycommentator().stop(MSG_STATUS (ret), (const char *) 0, "testMinPoly");
1627
1628 return ret;
1629 }
1630
1631 template <class Field>
testCharPoly(const Field & F,size_t n,int iterations)1632 static bool testCharPoly (const Field& F, size_t n, int iterations)
1633 {
1634 typedef typename Field::Element Element;
1635 typedef BlasMatrix<Field> Matrix;
1636 typedef typename Field::RandIter RandIter;
1637 // typedef BlasVector<Field> Polynomial;
1638 typedef typename Givaro::Poly1Dom<Field> PolRing;
1639 typedef typename PolRing::Element Polynomial;
1640 //Commentator mycommentator;
1641 mycommentator().getMessageClass (INTERNAL_DESCRIPTION).setMaxDepth (3);
1642 mycommentator().getMessageClass (INTERNAL_DESCRIPTION).setMaxDetailLevel (Commentator::LEVEL_NORMAL);
1643 mycommentator().start (pretty("Testing charpoly"),"testCharPoly",(unsigned int)iterations);
1644 Element tmp ;
1645 PolRing R(F);
1646 RandIter G(F);
1647 Givaro::GeneralRingNonZeroRandIter<Field> Gn(G);
1648 bool ret = true;
1649 BlasMatrixDomain<Field> BMD(F);
1650
1651 for (int k=0;k<iterations;++k) {
1652
1653 mycommentator().progress(k);
1654
1655 Matrix A(F,n,n);
1656 list<Polynomial> P;
1657 // Test CharPoly(In) = (X-1)^n
1658 for (size_t i=0;i<n;++i){
1659 for (size_t j=0;j<i;++j)
1660 A.setEntry(i,j,F.zero);
1661 A.setEntry(i,i,F.one);
1662 for (size_t j=i+1;j<n;++j)
1663 A.setEntry(i,j,F.zero);
1664 }
1665 P.clear();
1666
1667 BMD.charpoly( P, A );
1668
1669 auto P_it = P.begin();
1670 while (P_it != P.end()){
1671 if ( P_it->size() !=2 )
1672 ret = false;
1673 if ( !F.areEqual(P_it->operator[](0), F.mOne) )
1674 ret = false;
1675 if ( !F.areEqual(P_it->operator[](1), F.one) )
1676 ret = false;
1677
1678 ++P_it;
1679 }
1680
1681 // Test CharPoly(a*In) = X-a
1682 G.random(tmp);
1683
1684 for (size_t i=0;i<n;++i)
1685 A.setEntry(i,i,tmp);
1686 F.negin(tmp);
1687 P.clear();
1688 BMD.charpoly( P, A );
1689 P_it = P.begin();
1690 while (P_it != P.end()){
1691 if ( P_it->size() !=2 )
1692 ret = false;
1693 if ( !F.areEqual(P_it->operator[](0), tmp) )
1694 ret = false;
1695 if ( !F.areEqual(P_it->operator[](1), F.one) )
1696 ret = false;
1697 ++P_it;
1698 }
1699 }
1700
1701 mycommentator().stop(MSG_STATUS (ret), (const char *) 0, "testCharPoly");
1702
1703 return ret;
1704 }
1705
1706 template<class Field>
testBlasMatrixConstructors(const Field & Fld,size_t m,size_t n)1707 static bool testBlasMatrixConstructors(const Field& Fld, size_t m, size_t n)
1708 {
1709 bool pass = true;
1710 // typedef typename Field::Element Element;
1711 BlasMatrixDomain<Field> BMD(Fld);
1712 //BlasMatrix<Field> A; // nowhere to go
1713
1714 BlasMatrix<Field> B(Fld);
1715 B.resize(m, n, Fld.zero);
1716 // don't understand the variations on integer types for the indices...
1717
1718 BlasMatrix<Field> C(Fld,m,n);
1719 pass = pass and BMD.areEqual(B, C);
1720
1721 // MatrixStream<Field> ms; ...
1722 // BlasMatrix<Field> D(ms);
1723 // pass = pass and BMD.areEqual(B, D);
1724
1725 ScalarMatrix<Field> Eo(Fld, n, n, Fld.zero);
1726 BlasMatrix<Field> E(Eo); // copy a bb
1727 pass = pass and BMD.areEqual(B, E);
1728
1729 #if 0
1730 ScalarMatrix<Field> Fo(Fld, 2*m, 2*m, Fld.zero);
1731 BlasMatrix F(Fo, m, 0, m, n) ; // copy subm of a bb
1732 pass = pass and BMD.areEqual(B, F);
1733
1734 BlasMatrix<Field> G(Eo, Fld); // other field?
1735 pass = pass and BMD.areEqual(B, G);
1736
1737 BlasMatrix<Field> H(B); // copy cstor
1738 pass = pass and BMD.areEqual(B, H);
1739
1740 BlasVector<Field> v(Fld,m*n, Fld.zero);
1741 BlasMatrix<Field> I(Fld,v,m,n);
1742 pass = pass and BMD.areEqual(B, I);
1743
1744 Element *p = &v[0];
1745 BlasMatrix<Field> J(Fld,p,m,n);
1746 pass = pass and BMD.areEqual(B, J);
1747 #endif
1748
1749 return pass;
1750 }
1751
1752 // returns true if ok, false if not.
1753 template<class Field>
launch_tests(Field & F,size_t n,int iterations)1754 int launch_tests(Field & F, size_t n, int iterations)
1755 {
1756 bool pass = true ;
1757 //report << "no blas tests for now" << std::endl;
1758 // no slow test while I work on io
1759 if (!testBlasMatrixConstructors(F, n, n)) pass=false;
1760 if (!testMulAdd (F,n,iterations)) pass=false;
1761 if (0==F.characteristic()) if
1762 (!testMulAddAgain (F,n,iterations)) pass=false;
1763 size_t m = n+n/2 ; size_t k = 2*n+1 ;
1764 if (!testMulAddShapeTrans (F,n,m,k,iterations)) pass=false;
1765 if (!testMulAddShapeTrans (F,n,k,m,iterations)) pass=false;
1766 if (!testMulAddShapeTrans (F,m,n,k,iterations)) pass=false;
1767 if (!testMulAddShapeTrans (F,m,k,n,iterations)) pass=false;
1768 if (!testMulAddShapeTrans (F,k,n,m,iterations)) pass=false;
1769 if (!testMulAddShapeTrans (F,k,m,n,iterations)) pass=false;
1770 if (!testTriangMulShapeTrans<Field,true,true> (F,m,n,iterations)) pass=false;
1771 if (!testTriangMulShapeTrans<Field,true,true> (F,n,m,iterations)) pass=false;
1772 if (!testTriangMulShapeTrans<Field,false,true> (F,m,n,iterations)) pass=false;
1773 if (!testTriangMulShapeTrans<Field,false,true> (F,n,m,iterations)) pass=false;
1774 if (!testTriangMulShapeTrans<Field,true,false> (F,m,n,iterations)) pass=false;
1775 if (!testTriangMulShapeTrans<Field,true,false> (F,n,m,iterations)) pass=false;
1776 if (!testTriangMulShapeTrans<Field,false,false> (F,m,n,iterations)) pass=false;
1777 if (!testTriangMulShapeTrans<Field,false,false> (F,n,m,iterations)) pass=false;
1778 if (!testRank (F, n, iterations)) pass=false;
1779 if (!testDet (F, n, iterations)) pass=false;
1780 if (!testInv (F, n, iterations)) pass=false;
1781 if (!testTriangularSolve (F,n,n,iterations)) pass=false;
1782 if (!testSolve (F,n,n,iterations)) pass=false;
1783 if (!testPermutation (F,n,iterations)) pass=false;
1784 if (!testPLUQ (F,n,n,iterations)) pass=false;
1785 if (!testMinPoly (F,n,iterations)) pass=false;
1786 if (!testCharPoly (F,n,iterations)) pass=false;
1787 //
1788 //
1789 if (not pass) F.write(report) << endl;
1790 return pass ;
1791
1792 }
1793
launch_gf2_tests(GF2 & F,size_t n)1794 bool launch_gf2_tests(GF2 & F, size_t n)
1795 {
1796 bool pass = true;
1797 //pass = pass and testRank(F, n, 1);
1798 //pass = pass and testDet(F, n, 1);
1799 //pass = pass and testMulAdd(F, n, 1);
1800 return pass ;
1801 }
1802
1803 #if 0
1804 bool launch_gf3_tests(GF3 & F, size_t n)
1805 {
1806 bool pass = true;
1807 //pass = pass and testRank(F, n, 1);
1808 //pass = pass and launch_tests(F, n, 1); // fails with std BlasMatrixDomain
1809 return pass ;
1810 }
1811 #endif
1812
1813 // Local Variables:
1814 // mode: C++
1815 // tab-width: 4
1816 // indent-tabs-mode: nil
1817 // c-basic-offset: 4
1818 // End:
1819 // vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
1820