root / tests / testeur_lqup.C

Revision 62, 7.1 kB (checked in by pernet, 6 months ago)

* Change the design of fflas-bounds
* Modular<double> -> ModularBalanced?<double>
* Fix the winograd recursion issue (when some steps are done over the field)
* Fix a bunch of bugs
* Remove the template specialization by the Element (incompatible with the soon coming compressed representations over small fields)
* Create a randiter file, generic wrt the modular field

Line 
1/* -*- mode: C++; tab-width: 8; indent-tabs-mode: t; c-basic-offset: 8 -*- */
2//--------------------------------------------------------------------------
3//                        Test for the  lqup decomposition
4//                 
5//--------------------------------------------------------------------------
6// Clement Pernet
7//-------------------------------------------------------------------------
8
9
10#include <iostream>
11#include <iomanip>
12using namespace std;
13//#include "fflas-ffpack/modular-int.h"
14//#include "fflas-ffpack/modular-positive.h"
15#include "fflas-ffpack/modular-balanced.h"
16#include "timer.h"
17#include "Matio.h"
18#include "fflas-ffpack/ffpack.h"
19#include "givaro/givintprime.h"
20
21
22
23//typedef Modular<double> Field;
24typedef ModularBalanced<double> Field;
25//typedef Modular<float> Field;
26//typedef ModularBalanced<float> Field;
27//typedef Modular<int> Field;
28//typedef GivaroZpz<Std32> Field;
29//typedef GivaroGfq Field;
30
31int main(int argc, char** argv){
32        Timer tim;
33        IntPrimeDom IPD;
34        unsigned long p;
35        size_t M, N ;
36        bool keepon = true;
37        Integer _p,tmp;
38        Field::Element zero,one;
39        cerr<<setprecision(10);
40        size_t TMAX = 100;
41        size_t PRIMESIZE = 23;
42       
43        if (argc > 1 )
44                TMAX = atoi(argv[1]);
45        if (argc > 2 )
46                PRIMESIZE = atoi(argv[2]);
47
48        FFLAS::FFLAS_TRANSPOSE ta;
49        FFLAS::FFLAS_DIAG diag;
50        size_t lda;
51
52        Field::Element * A, *Abis, *X,* U, *L;
53        size_t *P, *Q;
54        while (keepon){
55                srandom(_p);
56                do{
57                        //              max = Integer::random(2);
58                        _p = random();//max % (2<<30);
59                        IPD.prevprime( tmp, (_p% (1<<PRIMESIZE)) );
60                        p =  tmp;
61                       
62                }while( (p <= 2) );
63               
64                Field F( p); 
65                F.init(zero,0.0);
66                F.init(one,1.0);
67                Field::RandIter RValue( F );
68               
69                do{
70                        M = (size_t)  random() % TMAX;
71                        N = (size_t)  random() % TMAX;
72                } while ((M == 0) || (N == 0));
73                lda = N;
74                if (random()%2)
75                        diag = FFLAS::FflasUnit;
76                else
77                        diag = FFLAS::FflasNonUnit;
78               
79
80                if (random()%2){
81                        ta = FFLAS::FflasTrans;
82                        L = new Field::Element[M*N];
83                        U = new Field::Element[N*N];
84                        P = new size_t[M];
85                        Q = new size_t[N];
86                        for (size_t i=0; i<M; ++i) P[i] = 0;
87                        for (size_t i=0; i<N; ++i) Q[i] = 0;
88                }
89                else{
90                        ta = FFLAS::FflasNoTrans;
91                        L = new Field::Element[M*M];
92                        U = new Field::Element[M*N];
93                        P = new size_t[N];
94                        Q = new size_t[M];
95                        for (size_t i=0; i<N; ++i) P[i] = 0;
96                        for (size_t i=0; i<M; ++i) Q[i] = 0;
97                }
98               
99                size_t R=0;
100                Field::Element * G = new Field::Element[M*M];
101                Field::Element * H = new Field::Element[M*N];
102                size_t t;
103                do{
104                        t = (size_t) random() % 10;
105                } while ((!t)||(t==1));
106                for (size_t i=0; i<M; ++i)
107                        if (!(random() % t))
108                                for (size_t j=0; j < M; ++j)
109                                        RValue.random (*(G+i*M+j));
110                        else
111                                for (size_t j=0; j < M; ++j)
112                                        F.assign(*(G+i*M+j), zero);
113                       
114               
115
116                for (size_t j=0; j < N; ++j)
117                        if (!(random() % t))
118                                for (size_t i=0; i<M; ++i)
119                                        RValue.random (*(H+i*N+j));
120                        else
121                                for (size_t i=0; i<M; ++i)
122                                        F.assign(*(H+i*N+j), zero);
123
124//              write_field(F,cerr<<"G = "<<endl,G,M,M,M);
125//              write_field(F,cerr<<"H = "<<endl,H,M,N,N);
126                A = new Field::Element[M*N];
127                FFLAS::fgemm (F, FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, M, N, M, one, G, M, H, N, zero, A, N);
128                delete[] G;
129                delete[] H;
130               
131                Abis = new Field::Element[M*N];
132                for (size_t i=0; i<M*N; ++i)
133                        *(Abis+i) = *(A+i);
134               
135                X = new Field::Element[M*N];
136               
137               
138                cout <<"p = "<<(size_t)p<<" M = "<<M
139                     <<" N = "<<N
140                     <<((diag==FFLAS::FflasUnit)?" Unit ":" Non Unit ")
141                     <<((ta==FFLAS::FflasNoTrans)?"LQUP ( A ) ":"LQUP ( A^T ) ")
142                     <<"...."; 
143
144                       
145                tim.clear();
146                tim.start();
147                R = FFPACK::LUdivine (F, diag, ta, M, N, A, lda, P, Q);
148                tim.stop();
149
150
151                //write_field(F,cerr<<"Result = "<<endl,Abis,M,N,lda);
152               
153                if (ta == FFLAS::FflasNoTrans){
154                       
155                        for (size_t i=0; i<R; ++i){
156                                for (size_t j=0; j<i; ++j)
157                                        F.assign ( *(U + i*N + j), zero);
158                                for (size_t j=i+1; j<N; ++j)
159                                        F.assign (*(U + i*N + j), *(A+ i*N+j));
160                        }
161                        for (size_t i=R;i<M; ++i)
162                                for (size_t j=0; j<N; ++j)
163                                        F.assign(*(U+i*N+j), zero);
164                        for ( size_t i=0; i<M; ++i ){
165                                size_t j=0;
166                                for (; j< ((i<R)?i:R) ; ++j )
167                                        F.assign( *(L + i*M+j), *(A+i*N+j));
168                                for (; j<M; ++j )
169                                        F.assign( *(L+i*M+j), zero);
170                        }
171                       
172                        //write_field(F,cerr<<"L = "<<endl,L,M,M,M);
173                        //write_field(F,cerr<<"U = "<<endl,U,M,N,N);
174                        FFPACK::applyP( F, FFLAS::FflasRight, FFLAS::FflasNoTrans, M,0,R, L, M, Q);
175                        for ( size_t  i=0; i<M; ++i )
176                                F.assign(*(L+i*(M+1)), one);
177                       
178                        if (diag == FFLAS::FflasNonUnit)
179                                for ( size_t  i=0; i<R; ++i )
180                                        F.assign (*(U+i*(N+1)), *(A+i*(lda+1)));
181                       
182                        else{
183                                for (size_t i=0; i<R; ++i ){
184                                        *(L+Q[i]*(M+1)) = *(A+Q[i]*lda+i);
185                                        F.assign (*(U+i*(N+1)),one);
186                                }
187                        }
188                       
189                        FFPACK::applyP (F, FFLAS::FflasRight, FFLAS::FflasNoTrans, M,0,R, U, N, P);
190                        FFPACK::applyP (F, FFLAS::FflasLeft, FFLAS::FflasTrans, N,0,R, U, N, Q);
191                        FFLAS::fgemm (F, FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, M,N,M, 1.0, L,M, U,N, 0.0, X,N);
192                        //delete[] A;
193                } else {
194                       
195                        for (size_t i=0; i<R; ++i){
196                                for (size_t j=0; j<i; ++j)
197                                        F.assign ( *(L + i + j*N), zero);
198                                for (size_t j=i+1; j<M; ++j)
199                                        F.assign (*(L + i + j*N), *(A+ i+j*N));
200                        }
201                       
202                        for (size_t i=R;i<N; ++i)
203                                for (size_t j=0; j<M; ++j)
204                                        F.assign(*(L+i+j*N), zero);
205                        for ( size_t i=0; i<N; ++i ){
206                                size_t j=0;
207                                for (;  j< ((i<R)?i:R) ; ++j )
208                                        F.assign( *(U + i+j*N), *(A+i+j*N));
209                                for (; j<N; ++j )
210                                        F.assign( *(U+i+j*N), zero);
211                        }
212                       
213                        FFPACK::applyP( F, FFLAS::FflasLeft, FFLAS::FflasTrans, N,0,R, U, N, Q);
214                        for (size_t i=0; i<N; ++i)
215                                F.assign (*(U+i*(N+1)),one);
216                        if (diag == FFLAS::FflasNonUnit)
217                                for ( size_t i=0; i<R; ++i )
218                                        F.assign (*(L+i*(N+1)), *(A+i*(lda+1)));
219                        else{
220                                for ( size_t i=0; i<R; ++i ){
221                                        *(U+Q[i]*(N+1)) = *(A+Q[i]+i*N);
222                                        F.assign (*(L+i*(N+1)),one);
223                                }
224                        }
225                        // write_field(F,cerr<<"L = "<<endl,L,M,N,N);
226//                      write_field(F,cerr<<"U = "<<endl,U,N,N,N);
227                       
228                        FFPACK::applyP (F, FFLAS::FflasLeft, FFLAS::FflasTrans, N,0,R, L, N, P);
229                        FFPACK::applyP (F, FFLAS::FflasRight, FFLAS::FflasNoTrans, M,0,R, L, N, Q);
230                        FFLAS::fgemm (F, FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, M,N,N, 1.0, L,N, U,N, 0.0, X,N);
231                }
232                for (size_t i=0; i<M; ++i)
233                        for (size_t j=0; j<N; ++j)
234                                if (!F.areEqual (*(Abis+i*N+j), *(X+i*N+j))){
235                                        cerr<<"error for i,j="<<i<<" "<<j<<" "<<*(Abis+i*N+j)<<" "<<*(X+i*N+j)<<endl;
236                                        keepon = false;
237                                }
238               
239                //write_field(F,cerr<<"X = "<<endl,X,m,n,n);
240                //write_field(F,cerr<<"B = "<<endl,B,m,n,n);
241               
242                if (keepon){
243                        cout<<"R = "<<R
244                            <<" Passed "
245                            <<(M*M/1000.0*(N-M/3.0)/tim.usertime()/1000.0)<<"Mfops"<<endl; 
246                        delete[] A;
247                        delete[] L;
248                        delete[] U;
249                        delete[] Abis;
250                        delete[] X;
251                        delete[] P;
252                        delete[] Q;
253                }
254                else{
255                        cerr<<"Abis = "<<endl;
256                        write_field( F, cerr, Abis, M, N, N );
257                        cerr<<"X = "<<endl;
258                        write_field( F, cerr, X, M, N, N );
259                }
260        }
261        cout<<endl;
262        cerr<<"FAILED with p = "<<(size_t)p<<" M = "<<M<<" N = "<<N
263            <<" trans = "<<ta<<" diag = "<<diag<<endl;
264
265        cerr<<"A:"<<endl;
266        cerr<<M<<" "<<N<<" M"<<endl;
267        for (size_t i=0; i<M; ++i)
268                for (size_t j=0; j<N; ++j)
269                        if (*(Abis+i*lda+j))
270                                cerr<<i+1<<" "<<j+1<<" "<<((int) *(Abis+i*lda+j) )<<endl;
271        cerr<<"0 0 0"<<endl<<endl;
272       
273        delete[] A;
274        delete[] Abis;
275        delete[] L;
276        delete[] U;
277        delete[] X;
278        delete[] P;
279        delete[] Q;
280}
281
282
283
284
285
286
287
288
289
290
291
292
293
Note: See TracBrowser for help on using the browser.