Show
Ignore:
Timestamp:
05/01/07 19:23:34 (2 years ago)
Author:
pernet
Message:

Introduction of the new automatic choice of underlying BLAS, for any finite finite field implementation:
float/double is chosen depending on the prime, the dimension, and efficiency considerations.

Files:
1 modified

Legend:

Unmodified
Added
Removed

  • include/fflas-ffpack/fflas_bounds.inl

    r18 r21  
    1919// Computes the maximal dimension k so that the product A*B+beta.C over Z, 
    2020// where A is m*k and B is k*n can be performed correctly with w Winograd 
    21 // recursion levels on the 53 bits of double mantissa 
     21// recursion levels on the number of bits of the floating point mantissa 
    2222//--------------------------------------------------------------------- 
    2323template  <class Field>  
    2424inline size_t FFLAS::DotProdBoundCompute (const Field& F, const size_t w,  
    25                                           const typename Field::Element& beta){ 
    26         return callDotProdBoundCompute<typename Field::Element>() (F, w, beta); 
    27 } 
    28  
    29 template<class Element> 
    30 class FFLAS::callDotProdBoundCompute { 
    31 public: 
    32         template  <class Field>  
    33         size_t operator () (const Field& F, const size_t w,  
    34                             const typename Field::Element& beta) 
    35         { 
    36                 typename Field::Element mone; 
    37                 static FFLAS_INT_TYPE p; 
    38                 F.characteristic(p); 
    39                 F.init (mone, -1.0); 
    40                 size_t kmax; 
    41                 if (p == 0) 
    42                         kmax = 2; 
    43                 else 
    44                         if (w > 0) { 
    45                                 size_t ex=1; 
    46                                 for (size_t i=0; i < w; ++i)    ex *= 3; 
    47                                 //FFLAS_INT_TYPE c = (p-1)*(ex)/2; //bound for a centered representation 
    48                                 long long c = (p-1)*(1+ex)/2; 
    49                                 kmax =  lround(( double(1ULL << DOUBLE_MANTISSA) /double(c*c) + 1)*(1 << w)); 
    50                                 if (kmax ==  ( 1ULL << w)) 
    51                                         kmax = 2; 
    52                         } 
    53                         else{ 
    54                                 long long c = p-1; 
    55                                 long long cplt=0; 
    56                                 if (!F.isZero (beta)) 
    57                                         if (F.isOne (beta) || F.areEqual (beta, mone)) 
    58                                                 cplt = c; 
    59                                         else cplt = c*c; 
    60                                 kmax =  lround(( double((1ULL << DOUBLE_MANTISSA) - cplt)) /double(c*c)); 
    61                                 if (kmax  < 2) 
    62                                         kmax = 2; 
    63                         } 
    64                 return  MIN(kmax,1ULL<<31); 
    65         } 
    66 }; 
    67  
    68 template  <>  
    69 class FFLAS::callDotProdBoundCompute<double> { 
    70 public: 
    71         template <class Field> 
    72         size_t operator() (const Field& F, const size_t w,  
    73                            const double& beta) 
    74         { 
    75                 double mone; 
    76                 static FFLAS_INT_TYPE p; 
    77                 F.characteristic(p); 
    78                 F.init (mone, -1.0); 
    79                 size_t  kmax; 
    80                 if (p == 0) 
    81                         kmax = 2; 
    82                 else 
    83                         if (w > 0) { 
    84                                 size_t ex=1; 
    85                                 for (size_t i=0; i < w; ++i)    ex *= 3; 
    86                                 long long c; 
     25                                          const typename Field::Element& beta, 
     26                                          const FFLAS_BASE base){ 
     27         
     28        typename Field::Element mone; 
     29        static FFLAS_INT_TYPE p; 
     30        F.characteristic(p); 
     31        F.init (mone, -1.0); 
     32        size_t kmax; 
     33 
     34        unsigned long mantissa = (base == FflasDouble)? DOUBLE_MANTISSA : FLOAT_MANTISSA; 
     35         
     36        if (p == 0) 
     37                kmax = 2; 
     38        else 
     39                if (w > 0) { 
     40                        size_t ex=1; 
     41                        for (size_t i=0; i < w; ++i)    ex *= 3; 
     42                        //FFLAS_INT_TYPE c = (p-1)*(ex)/2; //bound for a centered representation 
     43                        long long c; 
    8744#ifndef _LINBOX_CONFIG_H 
    88                                 if (F.balanced) 
    89                                         c = (p-1)*(ex)/2; // balanced representation 
    90                                 else 
     45                        if (F.balanced) 
     46                                c = (p-1)*(ex)/2; // balanced representation 
     47                        else 
    9148#endif 
    92                                         c = (p-1)*(1+ex)/2; // positive representation 
    93                                 kmax =  lround(double(1ULL << DOUBLE_MANTISSA) /(double(c*c) + 1)*(1ULL << w)); 
    94                                 if (kmax ==  ( 1ULL << w)) 
    95                                         kmax = 2; 
    96                         } 
    97                         else{ 
    98                                 long long c = p-1; 
    99                                 long long cplt=0; 
    100                                 if (!F.isZero (beta)) 
    101                                         if (F.isOne (beta) || F.areEqual (beta, mone)) 
    102                                                 cplt = c; 
    103                                         else cplt = c*c; 
    104                                 kmax =  lround( double((1ULL << DOUBLE_MANTISSA) - cplt) /(double(c*c))); 
    105                                 if (kmax  < 2) 
    106                                         kmax = 2; 
    107                         } 
    108                 return (size_t) MIN(kmax,1ULL<<31); 
    109         } 
    110 }; 
    111  
    112 template  <>  
    113 class FFLAS::callDotProdBoundCompute<float> { 
    114 public: 
    115         template <class Field> 
    116         size_t operator () (const Field& F, const size_t w,  
    117                             const float& beta) 
    118         { 
    119                 float mone,one; 
    120                 static FFLAS_INT_TYPE p; 
    121                 F.characteristic(p); 
    122                 F.init (one, 1.0F); 
    123                 F.neg(mone,one); 
    124                 size_t  kmax; 
    125                 if (p == 0) 
    126                         kmax = 2; 
    127                 else 
    128                         if (w > 0) { 
    129                                 size_t ex=1; 
    130                                 for (size_t i=0; i < w; ++i)    ex *= 3; 
    131                                 long long c; 
    132 #ifndef _LINBOX_CONFIG_H 
    133                                 if (F.balanced) 
    134                                         c = (p-1)*(ex)/2; // balanced representation 
    135                                 else 
    136 #endif 
    137                                         c = (p-1)*(1+ex)/2; // positive representation 
    138                                 kmax =  lround(float(1ULL << FLOAT_MANTISSA) /(float(c*c) + 1)*(1ULL << w)); 
    139                                 if (kmax ==  ( 1ULL << w)) 
    140                                         kmax = 2; 
    141                         } 
    142                         else{ 
    143                                 long long c = p-1; 
    144                                 long long cplt=0; 
    145                                 if (!F.isZero (beta)) 
    146                                         if (F.isOne (beta) || F.areEqual (beta, mone)) 
    147                                                 cplt = c; 
    148                                         else cplt = c*c; 
    149                                 kmax =  lround( float((1ULL << FLOAT_MANTISSA) - cplt) /(float(c*c))); 
    150                                 if (kmax  < 2) 
    151                                         kmax = 2; 
    152                         } 
    153                 return (size_t) MIN(kmax,1ULL<<31); 
    154         } 
    155 }; 
     49                                c = (p-1)*(1+ex)/2; 
     50                        kmax =  lround(( double(1ULL << mantissa) /double(c*c) + 1)*(1ULL << w)); 
     51                        if (kmax ==  ( 1ULL << w)) 
     52                                kmax = 2; 
     53                } 
     54                else{ 
     55                        long long c = p-1; 
     56                        long long cplt=0; 
     57                        if (!F.isZero (beta)) 
     58                                if (F.isOne (beta) || F.areEqual (beta, mone)) 
     59                                        cplt = c; 
     60                                else cplt = c*c; 
     61                        kmax =  lround(( double((1ULL << mantissa) - cplt)) /double(c*c)); 
     62                        if (kmax  < 2) 
     63                                kmax = 2; 
     64                } 
     65        //cerr<<"kmax = "<<kmax<<endl; 
     66        return  (size_t) MIN(kmax,1ULL<<31); 
     67} 
     68 
     69 
    15670 
    15771template  < class Field >  
    15872inline size_t FFLAS::DotProdBound (const Field& F, const size_t w,  
    159                                    const typename Field::Element& beta) 
     73                                   const typename Field::Element& beta, 
     74                                   const FFLAS_BASE base) 
    16075{ 
    16176        static Field G = F; 
    … …  
    16681        static typename Field::Element b = beta; 
    16782        static size_t w2 = w; 
    168         static size_t kmax = DotProdBoundCompute (F, w, beta); 
    169         if ((b != beta) || (pif != pig) ||  (w2 != w)) { 
     83        static FFLAS_BASE base2 = base; 
     84        static size_t kmax = DotProdBoundCompute (F, w, beta, base); 
     85        if ((b != beta) || (pif != pig) ||  (w2 != w) || (base2!=base)) { 
    17086                G = F; 
    17187                w2 = w; 
    17288                b = beta; 
    173                 kmax =  DotProdBoundCompute (F, w, beta); 
     89                base2 = base; 
     90                kmax =  DotProdBoundCompute (F, w, beta, base); 
    17491        }        
    17592        return kmax; 
    … …  
    268185                else 
    269186#ifdef _LINBOX_CONFIG_H 
    270                         return nmax= bound_compute (p=pi); //(F.balanced) ? bound_compute_balanced(p=pi) : bound_compute(p=pi); 
     187                        return nmax= bound_compute_double (p=pi); //(F.balanced) ? bound_compute_balanced(p=pi) : bound_compute(p=pi); 
    271188#else 
    272189                return (F.balanced) ? bound_compute_double_balanced(p=pi) : bound_compute_double(p=pi); 
    … …  
    292209                else 
    293210#ifdef _LINBOX_CONFIG_H 
    294                         return nmax= bound_compute (p=pi); //(F.balanced) ? bound_compute_balanced(p=pi) : bound_compute(p=pi); 
     211                        return nmax= bound_compute_float (p=pi); //(F.balanced) ? bound_compute_balanced(p=pi) : bound_compute(p=pi); 
    295212#else 
    296213                return (F.balanced) ? bound_compute_float_balanced(p=pi) : bound_compute_float(p=pi); 
    … …  
    299216}; 
    300217 
     218 
     219 
     220/** Automatic computation of the parameters for Matrix Multiplication 
     221 * 
     222 */ 
     223 
     224// To be tuned up: use statics to perform the computations only once for a given set of parameters 
     225template <class Element> 
     226class FFLAS::setMatMulParam { 
     227 public: 
     228        template <class Field> 
     229                void operator() (const Field& F, const size_t m, 
     230                                 const typename Field::Element& beta, 
     231                                 size_t& w, FFLAS_BASE& base, size_t& kmax){ 
     232 
     233                // Strategy : determine Winograd's recursion first, then choose appropriate 
     234                // floating point representation, and finally the blocking dimension. 
     235                // Can be improved for some cases. 
     236                static size_t ms = m; 
     237                static size_t ks = WinoSteps (m); 
     238 
     239                if (m != ms) 
     240                        ks = WinoSteps (ms = m); 
     241 
     242                w = ks; 
     243                 
     244                base = BaseCompute (F, w); 
     245                kmax = DotProdBound (F, w, beta, base); 
     246        } 
     247}; 
     248 
     249template<> 
     250class FFLAS::setMatMulParam<double> { 
     251public: 
     252        template <class Field> 
     253        void operator () (const Field& F, const size_t m, 
     254                          const typename Field::Element& beta, 
     255                          size_t& w, FFLAS_BASE& base, size_t& kmax){ 
     256 
     257                static size_t ms = m; 
     258                static size_t ks = WinoSteps (m); 
     259 
     260                if (m != ms) 
     261                        ks = WinoSteps (ms = m); 
     262 
     263                w = ks; 
     264                base = FflasDouble; 
     265                kmax = DotProdBound (F, w, beta, base); 
     266        } 
     267}; 
     268 
     269template<> 
     270class FFLAS::setMatMulParam<float> { 
     271public: 
     272        template <class Field> 
     273        void operator () (const Field& F, const size_t m, 
     274                          const typename Field::Element& beta, 
     275                          size_t& w, FFLAS_BASE& base, size_t& kmax){ 
     276 
     277                static size_t ms = m; 
     278                static size_t ks = WinoSteps (m); 
     279 
     280                if (m != ms) 
     281                        ks = WinoSteps (ms = m); 
     282 
     283                w = ks; 
     284                base = FflasFloat; 
     285                kmax = DotProdBound (F, w, beta, base); 
     286        } 
     287}; 
     288 
     289inline size_t FFLAS::WinoSteps (const size_t m) { 
     290 
     291        size_t w=0; 
     292        size_t mt = m; 
     293        while (mt >= WINOTHRESHOLD){ 
     294                w++; 
     295                mt/=2; 
     296        } 
     297        //cerr<<"Winostep = "<<w<<endl; 
     298        return w; 
     299} 
     300 
     301 
     302template <class Field> 
     303inline FFLAS::FFLAS_BASE FFLAS::BaseCompute (const Field& F, const size_t w){ 
     304 
     305        FFLAS_INT_TYPE pi; 
     306        F.characteristic(pi); 
     307        FFLAS_BASE base; 
     308        switch (w) { 
     309        case 0: base = (pi < FLOAT_DOUBLE_THRESHOLD_0)? FflasFloat : FflasDouble; 
     310                break; 
     311        case 1:  base = (pi < FLOAT_DOUBLE_THRESHOLD_1)? FflasFloat : FflasDouble; 
     312                break; 
     313        case 2:  base = (pi < FLOAT_DOUBLE_THRESHOLD_2)? FflasFloat : FflasDouble; 
     314                break; 
     315        default: base = FflasDouble; 
     316                break; 
     317        } 
     318        //cerr<<"base = "<<base<<endl; 
     319        return base; 
     320//      case 0: return (pi < FLOAT_DOUBLE_THRESHOLD_0)? FflasFloat : FflasDouble; 
     321//      case 1: return (pi < FLOAT_DOUBLE_THRESHOLD_1)? FflasFloat : FflasDouble; 
     322//      case 2: return (pi < FLOAT_DOUBLE_THRESHOLD_2)? FflasFloat : FflasDouble; 
     323//      default: return FflasDouble; 
     324         
     325} 
     326