root / tests / test-fgesv.C

Revision 51, 3.8 kB (checked in by pernet, 1 year ago)

Make test-fgesv work. No bugs known.

Line 
1/* -*- mode: C++; tab-width: 8; indent-tabs-mode: t; c-basic-offset: 8 -*- */
2//--------------------------------------------------------------------------
3//                        Test for fgesv : 1 computation
4//                 
5//--------------------------------------------------------------------------
6// Clement Pernet
7//-------------------------------------------------------------------------
8
9#define DEBUG 1
10#define TIME 1
11
12#include <iomanip>
13#include <iostream>
14using namespace std;
15
16#include "fflas-ffpack/modular-balanced.h"
17#include "timer.h"
18#include "Matio.h"
19#include "fflas-ffpack/ffpack.h"
20
21
22
23typedef Modular<double> Field;
24
25int main(int argc, char** argv){
26
27        int n,m,mb,nb;
28        cerr<<setprecision(10);
29        Field::Element zero, one;
30
31        if (argc != 6)  {
32                cerr<<"Usage : test-fgesv <p> <A> <b> <iter> <left/right>"
33                    <<endl;
34                exit(-1);
35        }
36        int nbit=atoi(argv[4]); // number of times the product is performed
37        Field F(atoi(argv[1]));
38        F.init(zero,0.0);
39        F.init(one,1.0);
40        Field::Element * A, *B, *B2, *X=NULL;
41        A = read_field(F,argv[2],&m,&n);
42        B = read_field(F,argv[3],&mb,&nb);
43
44        FFLAS::FFLAS_SIDE side = (atoi(argv[5])) ? FFLAS::FflasRight :  FFLAS::FflasLeft;
45               
46        size_t ldx=0;
47        size_t rhs = (side == FFLAS::FflasLeft) ? nb : mb;
48        if (m != n)
49                if (side == FFLAS::FflasLeft){
50                        X = new Field::Element[n*nb];
51                        ldx = nb;
52                } else {
53                        X = new Field::Element[mb*m];
54                        ldx = m;
55                }
56       
57        if ( ((side == FFLAS::FflasRight) && (n != nb))
58             || ((side == FFLAS::FflasLeft)&&(m != mb)) ) {
59                cerr<<"Error in the dimensions of the input matrices"<<endl;
60                exit(-1);
61        }
62        int info=0;
63        Timer t; t.clear();
64        double time=0.0;
65        //write_field(F, cerr<<"A="<<endl, A, k,k,k);
66        size_t R=0;
67        for (int i = 0;i<nbit;++i){
68                t.clear();
69                t.start();
70                if (m == n)
71                        R = FFPACK::fgesv (F, side, mb, nb, A, n, B, nb, &info);
72                else
73                        R = FFPACK::fgesv (F, side, m, n, rhs, A, n, X, ldx, B, nb, &info);
74                if (info > 0){
75                        std::cerr<<"System is inconsistent"<<std::endl;
76                        exit(-1);
77                }
78                       
79                t.stop();
80                time+=t.usertime();
81                if (i+1<nbit){
82                        delete[]A;
83                        A = read_field(F,argv[2],&m,&n);
84                        delete[] B;
85                        B = read_field(F,argv[3],&mb,&nb);
86                }
87        }
88
89#if DEBUG
90        delete[] A;
91
92        if (info > 0){
93                std::cerr<<"System inconsistent"<<std::endl;
94                exit (-1);
95        }
96               
97        A = read_field(F,argv[2],&m,&n);
98                       
99        B2 = new Field::Element[mb*nb];
100
101       
102        if (m==n)
103                if (side == FFLAS::FflasLeft)
104                        FFLAS::fgemm (F, FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, m, nb, n,
105                                      one, A, n, B, nb, zero, B2, nb);
106                else
107                        FFLAS::fgemm (F, FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, mb, n, m,
108                                      one, B, nb, A, n, zero, B2, nb);
109        else
110                if (side == FFLAS::FflasLeft)
111                        FFLAS::fgemm (F, FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, m, nb, n,
112                                      one, A, n, X, ldx, zero, B2, nb);
113                else
114                        FFLAS::fgemm (F, FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, mb, n, m,
115                                      one, X, ldx, A, n, zero, B2, nb);
116        delete[] B;
117        delete[] X;
118
119        B = read_field(F,argv[3],&mb,&nb);
120
121        bool wrong = false;
122        for (int i=0;i<mb;++i)
123                for (int j=0;j<nb;++j)
124                        if ( !F.areEqual(*(B2+i*nb+j), *(B+i*nb+j))){
125                                cerr<<"B2 ["<<i<<", "<<j<<"] = "<<(*(B2+i*nb+j))
126                                    <<" ; B ["<<i<<", "<<j<<"] = "<<(*(B+i*nb+j))
127                                    <<endl;
128                                wrong = true;
129                        }
130       
131        if (wrong) {
132                cerr<<"FAIL"<<endl;
133                //write_field (F,cerr<<"B2="<<endl,B2,m,n,n);
134                //write_field (F,cerr<<"B="<<endl,B,m,n,n);
135        }else{
136       
137                cerr<<"PASS"<<endl;
138        }
139
140
141#endif
142
143        delete[] A;
144        delete[] B;
145        delete[] B2;
146#if TIME
147        double mflops;
148        if (side == FFLAS::FflasLeft)
149                mflops = (n*m*m-m*m*m/3+2*R*R*n)/1000000.0*nbit/time;
150        else
151                mflops = (n*m*m-m*m*m/3+2*R*R*m)/1000000.0*nbit/time;
152        cerr<<"m,n,mb,nb = "<<m<<" "<<n<<" "<<mb<<" "<<nb<<". fgesv "
153            <<((side == FFLAS::FflasLeft)?" Left ":" Right ")
154            <<"over Z/"<<atoi(argv[1])<<"Z :"
155            <<endl
156            <<"t= "
157            << time/nbit 
158            << " s, Mffops = "<<mflops
159            << endl;
160       
161        cout<<m<<" "<<n<<" "<<mflops<<" "<<time/nbit<<endl;
162#endif
163}
Note: See TracBrowser for help on using the browser.