root / tests / test-lqup.C

Revision 74, 6.1 kB (checked in by pernet, 3 months ago)

fix minor bug

Line 
1/* -*- mode: C++; tab-width: 8; indent-tabs-mode: t; c-basic-offset: 8 -*- */
2//--------------------------------------------------------------------------
3//          Test for the lsp factorisation
4//--------------------------------------------------------------------------
5// usage: test-lsp p A n, for n lsp factorization 
6// of A over Z/pZ
7//-------------------------------------------------------------------------
8
9//-------------------------------------------------------------------------
10#define DEBUG 0
11// Debug option  0: no debug
12//               1: check A = LQUP
13//-------------------------------------------------------------------------
14using namespace std;
15
16
17//#define __LUDIVINE_CUTOFF 1
18#include <iostream>
19#include <iomanip>
20#include "Matio.h"
21#include "timer.h"
22#include "fflas-ffpack/modular-balanced.h"
23#include "fflas-ffpack/ffpack.h"
24
25//typedef ModularBalanced<double> Field;
26typedef ModularBalanced<float> Field;
27
28int main(int argc, char** argv){
29        cerr<<setprecision(20);
30        int m,n;
31        int R=0;
32
33        if (argc!=4){
34                cerr<<"usage : test-lqup <p> <A>  <i>"<<endl
35                    <<"        to do i LQUP factorisation of A"
36                    <<endl;
37                exit(-1);
38        }
39        Field F(atof(argv[1]));
40        Field::Element * A;
41       
42        A = read_field(F,argv[2],&m,&n);
43       
44        size_t maxP, maxQ;
45                       
46        //      size_t cutoff = atoi(argv[3]);
47        size_t nbf = atoi(argv[3]);
48       
49        Timer tim,timc;
50        timc.clear();
51
52        enum FFLAS::FFLAS_DIAG diag = FFLAS::FflasUnit;
53        enum FFLAS::FFLAS_TRANSPOSE trans = FFLAS::FflasNoTrans;
54        if (trans == FFLAS::FflasTrans){
55                maxP = m;
56                maxQ = n;
57        } else{
58                maxP = n;
59                maxQ = m;
60        }
61        size_t *P = new size_t[maxP];
62        size_t *Q = new size_t[maxQ];
63               
64        //write_field (F,cerr<<"A = "<<endl, A, m,n,n);
65        for ( size_t i=0;i<nbf;i++){
66                if (i) {               
67                        delete[] A;
68                        A = read_field(F,argv[2],&m,&n);
69                }
70                for (size_t j=0;j<maxP;j++)
71                        P[j]=0;
72                for (size_t j=0;j<maxQ;j++)
73                        Q[j]=0;
74                tim.clear();     
75                tim.start();   
76                R = FFPACK::LUdivine (F, diag, trans, m, n, A, n, P, Q,
77                                      FFPACK::FfpackLQUP);
78                tim.stop();
79                timc+=tim;
80        }
81        //write_field (F,cerr<<"Result = "<<endl, A, m,n,n);
82
83//      cerr<<"P = [";
84//      for (size_t i=0; i<maxP; ++i)
85//              cerr<<P[i]<<" ";
86//      cerr<<"]"<<endl;
87//      cerr<<"Q = [";
88//      for (size_t i=0; i<maxQ; ++i)
89//              cerr<<Q[i]<<" ";
90//      cerr<<"]"<<endl;
91#if DEBUG
92        Field::Element * X = new Field::Element[m*n];
93        Field::Element * L, *U;
94        if (trans == FFLAS::FflasNoTrans){
95                L = new Field::Element[m*m];
96                U = new Field::Element[m*n];
97                               
98                Field::Element zero,one;
99                F.init(zero,0.0);
100                F.init(one,1.0);
101                for (int i=0; i<R; ++i){
102                        for (int j=0; j<i; ++j)
103                                F.assign ( *(U + i*n + j), zero);
104                        for (int j=i+1; j<n; ++j)
105                                F.assign (*(U + i*n + j), *(A+ i*n+j));
106                }
107                for (int i=R;i<m; ++i)
108                        for (int j=0; j<n; ++j)
109                                F.assign(*(U+i*n+j), zero);
110                for ( int i=0; i<m; ++i ){
111                        int j=0;
112                        for (; j< ((i<R)?i:R) ; ++j )
113                                F.assign( *(L + i*m+j), *(A+i*n+j));
114                        for (; j<m; ++j )
115                                F.assign( *(L+i*m+j), zero);
116                }
117               
118                // write_field(F,cerr<<"L = "<<endl,L,m,m,m);
119//              write_field(F,cerr<<"U = "<<endl,U,m,n,n);
120                FFPACK::applyP( F, FFLAS::FflasRight, FFLAS::FflasNoTrans, m,0,R, L, m, Q);
121                for ( int i=0; i<m; ++i )
122                        F.assign(*(L+i*(m+1)), one);
123
124                //write_field(F,cerr<<"L = "<<endl,L,m,m,m);
125                //write_field(F,cerr<<"U = "<<endl,U,m,n,n);
126                if (diag == FFLAS::FflasNonUnit)
127                        for ( int i=0; i<R; ++i )
128                                F.assign (*(U+i*(n+1)), *(A+i*(n+1)));
129                       
130                else{
131                        for ( int i=0; i<R; ++i ){
132                                *(L+Q[i]*(m+1)) = *(A+Q[i]*n+i);
133                                F.assign (*(U+i*(n+1)),one);
134                        }
135                }
136//              write_field(F,cerr<<"L = "<<endl,L,m,m,m);
137//              write_field(F,cerr<<"U = "<<endl,U,m,n,n);
138
139                FFPACK::applyP (F, FFLAS::FflasRight, FFLAS::FflasNoTrans, m,0,R, U, n, P);
140                FFPACK::applyP (F, FFLAS::FflasLeft, FFLAS::FflasTrans, n,0,R, U, n, Q);
141                FFLAS::fgemm (F, FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, m,n,m, 1.0, L,m, U,n, 0.0, X,n);
142                //delete[] A;
143        } else {
144
145                L = new Field::Element[m*n];
146                U = new Field::Element[n*n];
147
148               
149                Field::Element zero,one;
150                F.init(zero,0.0);
151                F.init(one,1.0);
152                for (int i=0; i<R; ++i){
153                        for (int j=0; j<i; ++j)
154                                F.assign ( *(L + i + j*n), zero);
155                        for (int j=i+1; j<m; ++j)
156                                F.assign (*(L + i + j*n), *(A+ i+j*n));
157                }
158               
159                for (int i=R;i<n; ++i)
160                        for (int j=0; j<m; ++j)
161                                F.assign(*(L+i+j*n), zero);
162                for ( int i=0; i<n; ++i ){
163                        int j=0;
164                        for (;  j< ((i<R)?i:R) ; ++j )
165                                F.assign( *(U + i+j*n), *(A+i+j*n));
166                        for (; j<n; ++j )
167                                F.assign( *(U+i+j*n), zero);
168                }
169//              write_field(F,cerr<<"L = "<<endl,L,m,n,n);
170//              write_field(F,cerr<<"U = "<<endl,U,n,n,n);
171
172                FFPACK::applyP( F, FFLAS::FflasLeft, FFLAS::FflasTrans, n,0,R, U, n, Q);
173
174
175                for (int i=0; i<n; ++i)
176                        F.assign (*(U+i*(n+1)),one);
177
178                if (diag == FFLAS::FflasNonUnit)
179                        for ( int i=0; i<R; ++i )
180                                F.assign (*(L+i*(n+1)), *(A+i*(n+1)));
181                else{
182                        for ( int i=0; i<R; ++i ){
183                                *(U+Q[i]*(n+1)) = *(A+Q[i]+i*n);
184                                F.assign (*(L+i*(n+1)),one);
185                        }
186                }
187//              write_field(F,cerr<<"L = "<<endl,L,m,n,n);
188//              write_field(F,cerr<<"U = "<<endl,U,n,n,n);
189
190                FFPACK::applyP (F, FFLAS::FflasLeft, FFLAS::FflasTrans, n,0,R, L, n, P);
191                FFPACK::applyP (F, FFLAS::FflasRight, FFLAS::FflasNoTrans, m,0,R, L, n, Q);
192                FFLAS::fgemm (F, FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, m,n,n, 1.0, L,n, U,n, 0.0, X,n);
193        }
194        Field::Element * B =  read_field(F,argv[2],&m,&n);
195        bool fail = false;
196        for (int i=0; i<m; ++i)
197                for (int j=0; j<n; ++j)
198                        if (!F.areEqual (*(B+i*n+j), *(X+i*n+j))){
199                                std::cerr << " B["<<i<<","<<j<<"] = " << (*(B+i*n+j))
200                                          << " X["<<i<<","<<j<<"] = " << (*(X+i*n+j))
201                                          << endl;
202                                fail=true;
203                        }
204//      write_field(F,cerr<<"X = "<<endl,X,m,n,n);
205//      write_field(F,cerr<<"B = "<<endl,B,m,n,n);
206        delete[] B;
207        if (fail)
208                cerr<<"FAIL"<<endl;
209
210
211        else
212                cerr<<"PASS"<<endl;
213        delete[] U;
214        delete[] L;
215        delete[] X;
216#endif
217        delete[] A;
218        delete[] P;
219        delete[] Q;
220       
221        double t = timc.realtime();
222        double numops = m*m/1000.0*(n-m/3.0);
223       
224        cerr<<m<<"x"<< n
225            << " Trans = "<<trans
226            << " Diag = "<<diag
227            << " : rank = " << R << "  ["
228            << ((double)nbf/1000.0*(double)numops / t) 
229            << " MFops "
230            << " in "
231            << t/nbf<<"s"
232            <<"]"<< endl;
233//      cout<<m
234//          <<" "<<((double)nbf/1000.0*(double)numops / t)
235//          <<" "<<t/nbf
236//          <<endl;
237
238        return 0;
239}
Note: See TracBrowser for help on using the browser.