diff --git a/benchmarks/benchmark-fgemm.C b/benchmarks/benchmark-fgemm.C index 91ec24f94..21a336b33 100644 --- a/benchmarks/benchmark-fgemm.C +++ b/benchmarks/benchmark-fgemm.C @@ -82,6 +82,7 @@ int main(int argc, char** argv) { // typedef Givaro::Modular Field; // typedef Givaro::Modular Field; typedef Givaro::Modular Field; +// typedef Givaro::ModularExtended Field; // typedef Givaro::Modular Field; // typedef Givaro::ModularBalanced Field; // typedef Givaro::ModularBalanced Field; @@ -115,95 +116,95 @@ int main(int argc, char** argv) { for (size_t i=0;i<=iter;++i){ chrono.clear(); - if (p && p!=7){ - // CuttingStrategy meth = RECURSIVE; - // StrategyParameter strat = THREADS; - - typedef CuttingStrategy::Block block; - typedef CuttingStrategy::Recursive rec; - typedef StrategyParameter::Threads threads; - typedef StrategyParameter::TwoD twod; - typedef StrategyParameter::TwoDAdaptive twoda; - typedef StrategyParameter::ThreeD threed; - typedef StrategyParameter::ThreeDAdaptive threeda; - typedef StrategyParameter::ThreeDInPlace threedip; - PAR_BLOCK{ - if (i) { chrono.start(); } - - switch (p){ - case 1:{ - MMHelper::value, ParSeqHelper::Parallel > WH(F,nbw, SPLITTER(t,block,threads)); - fgemm (F, FflasNoTrans, FflasNoTrans, m,n,k, F.one, A, k, B, n, F.zero, C,n, WH); - break;} - case 2:{ - MMHelper::value, ParSeqHelper::Parallel > WH(F,nbw, SPLITTER(t,rec,twod)); - fgemm (F, FflasNoTrans, FflasNoTrans, m,n,k, F.one, A, k, B, n, F.zero, C,n, WH); - break; - } - case 3:{ - MMHelper::value, ParSeqHelper::Parallel > WH(F,nbw, SPLITTER(t,rec,twoda)); - fgemm (F, FflasNoTrans, FflasNoTrans, m,n,k, F.one, A, k, B, n, F.zero, C,n, WH); - break; - } - case 4:{ - MMHelper::value, ParSeqHelper::Parallel > WH(F,nbw, SPLITTER(t,rec,threedip)); - fgemm (F, FflasNoTrans, FflasNoTrans, m,n,k, F.one, A, k, B, n, F.zero, C,n, WH); - break; - } - case 5:{ - MMHelper::value, ParSeqHelper::Parallel > WH(F,nbw, SPLITTER(t,rec,threed)); - fgemm (F, FflasNoTrans, FflasNoTrans, m,n,k, F.one, A, k, B, n, F.zero, C,n, WH); - break; - } - case 6:{ - MMHelper::value, ParSeqHelper::Parallel > WH(F,nbw, SPLITTER(t,rec,threeda)); - fgemm (F, FflasNoTrans, FflasNoTrans, m,n,k, F.one, A, k, B, n, F.zero, C,n, WH); - break; - } - default:{ - MMHelper::value, ParSeqHelper::Parallel > WH(F,nbw, SPLITTER(t,block,threads)); - fgemm (F, FflasNoTrans, FflasNoTrans, m,n,k, F.one, A, k, B, n, F.zero, C,n, WH); - break; - } - } - } - if (i) {chrono.stop(); time[i-1]=chrono.realtime();} - }else{ - if(p==7){ - - int nrec = 0; - int dim = m; - // if(dim < 19000) - nrec--; - while(dim >= __FFLASFFPACK_WINOTHRESHOLD*2){ - dim=dim/2; - nrec++; - } - nrec=std::max(1,nrec); - // std::cout<<" WINO_THREShold"<<__FFLASFFPACK_WINOTHRESHOLD<<" nrec = "<::value, ParSeqHelper::Parallel > WH(F,nbw, SPLITTER(t,block,threads)); + // fgemm (F, FflasNoTrans, FflasNoTrans, m,n,k, F.one, A, k, B, n, F.zero, C,n, WH); + // break;} + // case 2:{ + // MMHelper::value, ParSeqHelper::Parallel > WH(F,nbw, SPLITTER(t,rec,twod)); + // fgemm (F, FflasNoTrans, FflasNoTrans, m,n,k, F.one, A, k, B, n, F.zero, C,n, WH); + // break; + // } + // case 3:{ + // MMHelper::value, ParSeqHelper::Parallel > WH(F,nbw, SPLITTER(t,rec,twoda)); + // fgemm (F, FflasNoTrans, FflasNoTrans, m,n,k, F.one, A, k, B, n, F.zero, C,n, WH); + // break; + // } + // case 4:{ + // MMHelper::value, ParSeqHelper::Parallel > WH(F,nbw, SPLITTER(t,rec,threedip)); + // fgemm (F, FflasNoTrans, FflasNoTrans, m,n,k, F.one, A, k, B, n, F.zero, C,n, WH); + // break; + // } + // case 5:{ + // MMHelper::value, ParSeqHelper::Parallel > WH(F,nbw, SPLITTER(t,rec,threed)); + // fgemm (F, FflasNoTrans, FflasNoTrans, m,n,k, F.one, A, k, B, n, F.zero, C,n, WH); + // break; + // } + // case 6:{ + // MMHelper::value, ParSeqHelper::Parallel > WH(F,nbw, SPLITTER(t,rec,threeda)); + // fgemm (F, FflasNoTrans, FflasNoTrans, m,n,k, F.one, A, k, B, n, F.zero, C,n, WH); + // break; + // } + // default:{ + // MMHelper::value, ParSeqHelper::Parallel > WH(F,nbw, SPLITTER(t,block,threads)); + // fgemm (F, FflasNoTrans, FflasNoTrans, m,n,k, F.one, A, k, B, n, F.zero, C,n, WH); + // break; + // } + // } + // } + // if (i) {chrono.stop(); time[i-1]=chrono.realtime();} + // }else{ + // if(p==7){ + + // int nrec = 0; + // int dim = m; + // // if(dim < 19000) + // nrec--; + // while(dim >= __FFLASFFPACK_WINOTHRESHOLD*2){ + // dim=dim/2; + // nrec++; + // } + // nrec=std::max(1,nrec); + // // std::cout<<" WINO_THREShold"<<__FFLASFFPACK_WINOTHRESHOLD<<" nrec = "<//, //typename FieldTraits::value, @@ -212,10 +213,10 @@ int main(int argc, char** argv) { if (i) chrono.start(); fgemm (F, FflasNoTrans, FflasNoTrans, m,n,k, F.one, A, k, B, n, F.zero, C,n,WH); if (i) {chrono.stop(); time[i-1]=chrono.realtime();} - } - } + //} + - TimFreivalds.clear(); + TimFreivalds.clear(); TimFreivalds.start(); bool pass = freivalds(F, FflasNoTrans, FflasNoTrans, m,n,k, F.one, A, k, B, n, C,n); diff --git a/fflas-ffpack/fflas/fflas_fadd.h b/fflas-ffpack/fflas/fflas_fadd.h index b693ff13a..59d083a3b 100644 --- a/fflas-ffpack/fflas/fflas_fadd.h +++ b/fflas-ffpack/fflas/fflas_fadd.h @@ -106,12 +106,11 @@ namespace FFLAS { // C = A + a B template - void - fadd (const Field& F, const size_t N, - typename Field::ConstElement_ptr A, const size_t inca, - const typename Field::Element alpha, - typename Field::ConstElement_ptr B, const size_t incb, - typename Field::Element_ptr C, const size_t incc) + void fadd (const Field& F, const size_t N, + typename Field::ConstElement_ptr A, const size_t inca, + const typename Field::Element alpha, + typename Field::ConstElement_ptr B, const size_t incb, + typename Field::Element_ptr C, const size_t incc) { if (C == A && inca == incc) return faxpy(F,N,alpha,B,incb,C,incc); @@ -125,9 +124,10 @@ namespace FFLAS { if (inca == 1 && incb == 1 && incc == 1) { for (size_t i = 0 ; i < N ; ++i) { - //!@todo optimise here - F.mul(C[i],alpha,B[i]); - F.addin(C[i],A[i]); + //!@todo optimise here + F.axpy (C[i], alpha, B[i], A[i]); + // F.mul(C[i],alpha,B[i]); + // F.addin(C[i],A[i]); } return; } @@ -302,21 +302,26 @@ namespace FFLAS { return faxpy(F,M,N,alpha,B,ldb,C,ldc); if (F.isOne(alpha)) return fadd(F,M,N,A,lda,B,ldb,C,ldc); - if (F.isMOne(alpha)) + if (F.isMOne(alpha)){ return fsub(F,M,N,A,lda,B,ldb,C,ldc); + } if (F.isZero(alpha)) return fassign(F,M,N,A,lda,C,ldc); if (N == lda && N == ldb && N == ldc) return fadd(F,M*N,A,1,alpha,B,1,C,1); - +//TODO: use level 1 row by row when lda!=n typename Field::ConstElement_ptr Ai = A, Bi = B; typename Field::Element_ptr Ci = C; - for (; Ai < A+M*lda; Ai+=lda, Bi+=ldb, Ci+=ldc) - for (size_t i=0; i - inline bool NeedDoublePreAddReduction (Element& Outmin, Element& Outmax, + inline int NeedDoublePreAddReduction (Element& Outmin, Element& Outmax, Element& Op1min, Element& Op1max, Element& Op2min, Element& Op2max, Element beta, MMHelper& WH) @@ -168,27 +171,40 @@ namespace FFLAS{ namespace Protected{ // Testing if P5 need to be reduced Outmin = std::min(beta*Op2min,beta*Op2max); Outmax = std::max(beta*Op2min,beta*Op2max); - if (Op1max > WH.MaxStorableValue-Outmax || - -Op1min > WH.MaxStorableValue+Outmin){ + // std::cerr<? "<< WH.MaxStorableValue <<" - "<? "<< WH.MaxStorableValue <<" + "< WH.MaxStorableValue-Outmax || -Op1min > WH.MaxStorableValue+Outmin){ Outmin += WH.FieldMin; Outmax += WH.FieldMax; - return true; + return 1; } else{ Outmin += Op1min; Outmax += Op1max; - return false; + return 0; } } template - inline bool NeedDoublePreAddReduction (Element& Outmin, Element& Outmax, + inline int NeedDoublePreAddReduction (Element& Outmin, Element& Outmax, Element& Op1min, Element& Op1max, Element& Op2min, Element& Op2max, Element beta, MMHelper& WH) { Outmin = WH.FieldMin; Outmax = WH.FieldMax; - return false; + return 0; + } + template + inline int NeedDoublePreAddReduction (Element& Outmin, Element& Outmax, + Element& Op1min, Element& Op1max, + Element& Op2min, Element& Op2max, Element beta, + MMHelper, AlgoT, ModeCategories::LazyTag, ParSeqTrait >& WH){ + Outmin = WH.FieldMin; + Outmax = WH.FieldMax; + return 2; } template diff --git a/fflas-ffpack/fflas/fflas_fgemm/fgemm_classical.inl b/fflas-ffpack/fflas/fflas_fgemm/fgemm_classical.inl index db1d79282..88ff68c8a 100644 --- a/fflas-ffpack/fflas/fflas_fgemm/fgemm_classical.inl +++ b/fflas-ffpack/fflas/fflas_fgemm/fgemm_classical.inl @@ -80,9 +80,10 @@ namespace FFLAS { if (F.isMOne(betadf)) betadf = -F.one; size_t kmax = H.MaxDelayedDim (betadf); + H.checkA(F,ta, m,k,A,lda); H.checkB(F,tb, k,n,B,ldb); - if (kmax <= k/2 || H.Aunfit() || H.Bunfit() ){ + if (kmax <= k/2 || H.Aunfit() || H.Bunfit() ){ // TODO: K/2 may be too big // Might as well reduce inputs if (H.Amin < H.FieldMin || H.Amax>H.FieldMax){ H.initA(); @@ -207,8 +208,9 @@ namespace FFLAS { if (tb == FflasNoTrans) for (size_t i = 0; i < m; ++i) for (size_t l = 0; l < k; ++l) - for (size_t j = 0; j < n; ++j) - F.axpyin (*(C+i*ldc+j), *(A+i*lda+l), *(B+l*ldb+j)); + for (size_t j = 0; j < n; ++j) { + F.axpyin (*(C+i*ldc+j), *(A+i*lda+l), *(B+l*ldb+j)); + } else for (size_t i = 0; i < m; ++i) for (size_t j = 0; j < n; ++j) diff --git a/fflas-ffpack/fflas/fflas_fgemm/fgemm_winograd.inl b/fflas-ffpack/fflas/fflas_fgemm/fgemm_winograd.inl index 37a794afd..b69a42dd9 100644 --- a/fflas-ffpack/fflas/fflas_fgemm/fgemm_winograd.inl +++ b/fflas-ffpack/fflas/fflas_fgemm/fgemm_winograd.inl @@ -198,8 +198,7 @@ namespace FFLAS { namespace Protected { const typename Field::Element beta, typename Field::Element_ptr C, const size_t ldc, MMHelper & H, - const typename MMHelper::DelayedField::Element Cmin, - const typename MMHelper::DelayedField::Element Cmax) + MMHelper & Hsave) { size_t mkn =(size_t)( (bool)(nr > 0)+ ((bool)(kr > 0) << 1)+ ((bool)(mr > 0) << 2)); if (mkn == 0) return; @@ -226,12 +225,12 @@ namespace FFLAS { namespace Protected { MMHelper HModd(H); MMHelper HNodd(H); - Hacc.Cmin = H.Outmin; Hacc.Cmax = H.Outmax; - Hacc.recLevel=-1;HModd.recLevel=-1;HNodd.recLevel=-1; - HModd.Cmin = Cmin; HModd.Cmax = Cmax; - HModd.Amax = H.Bmax; HModd.Amin = H.Bmin; - HModd.Bmax = H.Amax; HModd.Bmin = H.Amin; - HNodd.Cmin = Cmin; HNodd.Cmax = Cmax; + copyAccumulator (true, H, Hacc); // H.Out -> Hacc.C + + Hacc.recLevel=-1; HModd.recLevel=-1; HNodd.recLevel=-1; + + copyAccumulator (false, Hsave, HModd); // Hsave.C -> HModd.C + copyAccumulator (false, Hsave, HNodd); // Hsave.C -> HNodd.C switch (mkn) { case 1: // n oddsized @@ -270,9 +269,7 @@ namespace FFLAS { namespace Protected { fgemm (F, ta, tb, m, nr, k, alpha, A, lda, b12, ldb, beta, C+(n-nr), ldc, HNodd); break; } - H.Outmin = min4(HModd.Outmin,HNodd.Outmin, Hacc.Outmin, H.Outmin); - H.Outmax = max4(HModd.Outmax,HNodd.Outmax, Hacc.Outmax, H.Outmax); - H.checkOut(F, m,n, C, ldc); + mergeOutBounds (HModd, HNodd, Hacc, H); } // #define NEWIP @@ -398,17 +395,11 @@ namespace FFLAS{ if (H.recLevel == 0){ MMHelper HC(H); fgemm (F, ta, tb, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc, HC); - // copyOutBounds (HC, H); @ todo should replace the following 2 lines. - // Not ready yet, because of API of DynamicPeeling - H.Outmax = HC.Outmax; - H.Outmin = HC.Outmin; + copyOutBounds(HC,H); return C; } // Then w >0 - typedef typename MMHelper::DelayedField::Element DFElt; - DFElt Cmin = H.Cmin; - DFElt Cmax = H.Cmax; #ifdef OLD_DYNAMIC_PEELING @@ -425,6 +416,8 @@ namespace FFLAS{ size_t n2 = (n >> ww) << (ww-1) ; size_t k2 = (k >> ww) << (ww-1) ; + MMHelper Hsave (H); + Protected::WinogradCalc (F, ta, tb, m2, n2, k2, alpha, A, lda, B, ldb, beta, C, ldc, H); size_t mr = m -2*m2; @@ -435,7 +428,7 @@ namespace FFLAS{ FFLASFFPACK_check(n == n2*2+nr); FFLASFFPACK_check(k == k2*2+kr); - Protected::DynamicPeeling2 (F, ta, tb, m, n, k, mr, nr, kr, alpha, A, lda, B, ldb, beta, C, ldc, H, Cmin, Cmax); + Protected::DynamicPeeling2 (F, ta, tb, m, n, k, mr, nr, kr, alpha, A, lda, B, ldb, beta, C, ldc, H, Hsave); #endif return C; } // fgemm @@ -494,20 +487,15 @@ namespace FFLAS{ FFLAS::ParSeqHelper::Parallel > HC (F, 0, ParSeqHelper::Parallel(NUM_THREADS)); #endif - // MMHelper HC(H); fgemm (F, ta, tb, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc, HC); - H.Outmax = HC.Outmax; - H.Outmin = HC.Outmin; + copyOutBounds(H,HC); return C; } /// // Then w >0 - typedef typename MMHelper::DelayedField::Element DFElt; - DFElt Cmin = H.Cmin; - DFElt Cmax = H.Cmax; #ifdef OLD_DYNAMIC_PEELING BLAS3::WinoPar (F, ta, tb, m/2, n/2, k/2, alpha, A, lda, B, ldb, beta, C, ldc, H); @@ -523,6 +511,8 @@ namespace FFLAS{ size_t n2 = (n >> ww) << (ww-1) ; size_t k2 = (k >> ww) << (ww-1) ; + MMHelper Hsave(H); + BLAS3::WinoPar (F, ta, tb, m2, n2, k2, alpha, A, lda, B, ldb, beta, C, ldc, H); size_t mr = m -2*m2; @@ -533,7 +523,7 @@ namespace FFLAS{ FFLASFFPACK_check(n == n2*2+nr); FFLASFFPACK_check(k == k2*2+kr); MMHelper HC(H); - Protected::DynamicPeeling2 (F, ta, tb, m, n, k, mr, nr, kr, alpha, A, lda, B, ldb, beta, C, ldc, HC, Cmin, Cmax); + Protected::DynamicPeeling2 (F, ta, tb, m, n, k, mr, nr, kr, alpha, A, lda, B, ldb, beta, C, ldc, HC, Hsave); #endif return C; } // fgemm diff --git a/fflas-ffpack/fflas/fflas_fgemm/schedule_winograd.inl b/fflas-ffpack/fflas/fflas_fgemm/schedule_winograd.inl index fbaffa805..24909a0da 100644 --- a/fflas-ffpack/fflas/fflas_fgemm/schedule_winograd.inl +++ b/fflas-ffpack/fflas/fflas_fgemm/schedule_winograd.inl @@ -339,20 +339,20 @@ namespace FFLAS { namespace BLAS3 { return C; } //wino parallel - + // Specialization for Delayed or Lazy fields, requiring management of bounds template < class Field, class FieldTrait > - inline void Winograd (const Field& F, - const FFLAS_TRANSPOSE ta, - const FFLAS_TRANSPOSE tb, - const size_t mr, const size_t nr, const size_t kr, - const typename Field::Element alpha, - typename Field::ConstElement_ptr A,const size_t lda, - typename Field::ConstElement_ptr B,const size_t ldb, - const typename Field::Element beta, - typename Field::Element_ptr C, const size_t ldc, - // const size_t kmax, const size_t w, const FFLAS_BASE base - MMHelper & WH - ) + inline typename std::enable_if::value, void>::type + Winograd (const Field& F, + const FFLAS_TRANSPOSE ta, + const FFLAS_TRANSPOSE tb, + const size_t mr, const size_t nr, const size_t kr, + const typename Field::Element alpha, + typename Field::ConstElement_ptr A,const size_t lda, + typename Field::ConstElement_ptr B,const size_t ldb, + const typename Field::Element beta, + typename Field::Element_ptr C, const size_t ldc, + // const size_t kmax, const size_t w, const FFLAS_BASE base + MMHelper & WH) { FFLASFFPACK_check(F.isZero(beta)); @@ -535,6 +535,143 @@ namespace FFLAS { namespace BLAS3 { } // Winograd + + // Case for fields not requiring management of bounds + template < class Field, class FieldTrait > + inline typename std::enable_if::value, void>::type + Winograd (const Field& F, + const FFLAS_TRANSPOSE ta, + const FFLAS_TRANSPOSE tb, + const size_t mr, const size_t nr, const size_t kr, + const typename Field::Element alpha, + typename Field::ConstElement_ptr A,const size_t lda, + typename Field::ConstElement_ptr B,const size_t ldb, + const typename Field::Element beta, + typename Field::Element_ptr C, const size_t ldc, + // const size_t kmax, const size_t w, const FFLAS_BASE base + MMHelper & WH + ) + { + FFLASFFPACK_check(F.isZero(beta)); + + typedef MMHelper MMH_t; + MMH_t H = WH ; + H.recLevel--; + + size_t lb, cb, la, ca, ldX2; + // size_t x3rd = std::max(mr,kr); + typename Field::ConstElement_ptr A11=A, A12, A21, A22; + typename Field::ConstElement_ptr B11=B, B12, B21, B22; + + typename Field::Element_ptr C11=C, C12=C+nr, C21=C+mr*ldc, C22=C21+nr; + + size_t x1rd = std::max(nr,kr); + size_t ldX1; + if (ta == FflasTrans) { + A21 = A + mr; + A12 = A + kr*lda; + A22 = A12 + mr; + la = kr; + ca = mr; + ldX1 = mr; + } else { + A12 = A + kr; + A21 = A + mr*lda; + A22 = A21 + kr; + la = mr; + ca = kr; + ldX1 = x1rd; + } + if (tb == FflasTrans) { + B21 = B + kr; + B12 = B + nr*ldb; + B22 = B12 + kr; + lb = nr; + cb = kr; + ldX2 = kr; + } else { + B12 = B + nr; + B21 = B + kr*ldb; + B22 = B21 + nr; + lb = kr; + ldX2 = cb = nr; + } + // Two temporary submatrices are required + typename Field::Element_ptr X2 = fflas_new (F, kr, nr); + + // T3 = B22 - B12 in X2 + fsub(F,lb,cb, B22,ldb, B12,ldb, X2,ldX2); + + // S3 = A11 - A21 in X1 + typename Field::Element_ptr X1 = fflas_new (F,mr,x1rd); + fsub(F,la,ca,A11,lda,A21,lda,X1,ldX1); + + // P7 = alpha . S3 * T3 in C21 + fgemm (F, ta, tb, mr, nr, kr, alpha, X1, ldX1, X2, ldX2, F.zero, C21, ldc, H); + + // T1 = B12 - B11 in X2 + fsub(F,lb,cb,B12,ldb,B11,ldb,X2,ldX2); + + // S1 = A21 + A22 in X1 + fadd(F,la,ca,A21,lda,A22,lda,X1,ldX1); + + // P5 = alpha . S1*T1 in C22 + fgemm (F, ta, tb, mr, nr, kr, alpha, X1, ldX1, X2, ldX2, F.zero, C22, ldc, H); + + // T2 = B22 - T1 in X2 + fsub(F,lb,cb,B22,ldb,X2,ldX2,X2,ldX2); + + // S2 = S1 - A11 in X1 + fsubin(F,la,ca,A11,lda,X1,ldX1); + + // P6 = alpha . S2 * T2 in C12 + fgemm (F, ta, tb, mr, nr, kr, alpha, X1, ldX1, X2, ldX2, F.zero, C12, ldc, H); + + // S4 = A12 -S2 in X1 + fsub(F,la,ca,A12,lda,X1,ldX1,X1,ldX1); + + // P3 = alpha . S4*B22 in C11 + fgemm (F, ta, tb, mr, nr, kr, alpha, X1, ldX1, B22, ldb, F.zero, C11, ldc, H); + + // P1 = alpha . A11 * B11 in X1 + fgemm (F, ta, tb, mr, nr, kr, alpha, A11, lda, B11, ldb, F.zero, X1, nr, H); + + // U2 = P1 + P6 in C12 and + faddin(F,mr,nr,X1,nr,C12,ldc); + + // U3 = P7 + U2 in C21 and + faddin(F,mr,nr,C12,ldc,C21,ldc); + + + // U4 = P5 + U2 in C12 and + faddin(F,mr,nr,C22,ldc,C12,ldc); + + // U7 = P5 + U3 in C22 and + faddin(F,mr,nr,C21,ldc,C22,ldc); + + // U5 = P3 + U4 in C12 + faddin(F,mr,nr,C11,ldc,C12,ldc); + + // T4 = T2 - B21 in X2 + fsubin(F,lb,cb,B21,ldb,X2,ldX2); + + // P4 = alpha . A22 * T4 in C11 + fgemm (F, ta, tb, mr, nr, kr, alpha, A22, lda, X2, ldX2, F.zero, C11, ldc, H); + + fflas_delete (X2); + + // U6 = U3 - P4 in C21 + fsubin(F,mr,nr,C11,ldc,C21,ldc); + + // P2 = alpha . A12 * B21 in C11 + fgemm (F, ta, tb, mr, nr, kr, alpha, A12, lda, B21, ldb, F.zero, C11, ldc, H); + + // U1 = P2 + P1 in C11 + faddin(F,mr,nr,X1,nr,C11,ldc); + + fflas_delete (X1); + } // Winograd + } // BLAS3 diff --git a/fflas-ffpack/fflas/fflas_fgemm/schedule_winograd_acc.inl b/fflas-ffpack/fflas/fflas_fgemm/schedule_winograd_acc.inl index 4d0f7c263..253f1e151 100644 --- a/fflas-ffpack/fflas/fflas_fgemm/schedule_winograd_acc.inl +++ b/fflas-ffpack/fflas/fflas_fgemm/schedule_winograd_acc.inl @@ -184,18 +184,19 @@ namespace FFLAS { namespace BLAS3 { } // WinogradAccOld // 3 temps and 21 ops + // Specialization for Delayed or Lazy fields, requiring management of bounds template < class Field, class FieldTrait> - inline void WinogradAcc_3_21 (const Field& F, - const FFLAS_TRANSPOSE ta, - const FFLAS_TRANSPOSE tb, - const size_t mr, const size_t nr, const size_t kr, - const typename Field::Element alpha, - typename Field::ConstElement_ptr A,const size_t lda, - typename Field::ConstElement_ptr B,const size_t ldb, - const typename Field::Element beta, - typename Field::Element_ptr C, const size_t ldc, - MMHelper & WH - ) + inline typename std::enable_if::value, void>::type + WinogradAcc_3_21 (const Field& F, + const FFLAS_TRANSPOSE ta, + const FFLAS_TRANSPOSE tb, + const size_t mr, const size_t nr, const size_t kr, + const typename Field::Element alpha, + typename Field::ConstElement_ptr A,const size_t lda, + typename Field::ConstElement_ptr B,const size_t ldb, + const typename Field::Element beta, + typename Field::Element_ptr C, const size_t ldc, + MMHelper & WH) { typedef MMHelper MMH_t; typedef typename MMH_t::DelayedField::Element_ptr DFEptr; @@ -205,7 +206,6 @@ namespace FFLAS { namespace BLAS3 { const typename MMH_t::DelayedField & DF = WH.delayedField; FFLASFFPACK_check(!DF.isZero(beta)); - size_t lb, cb, la, ca; size_t x3rd = std::max(mr,kr); typename Field::ConstElement_ptr A11=A, A12, A21, A22; @@ -221,7 +221,6 @@ namespace FFLAS { namespace BLAS3 { } else { DF.assign(betadf, beta); } - size_t ldX3; if (ta == FflasTrans) { @@ -271,21 +270,33 @@ namespace FFLAS { namespace BLAS3 { WH.Bmax-WH.Bmin, 0, 0); fgemm (F, ta, tb, mr, nr, kr, alpha, X2, ca, X3, ldX3, F.zero, X1, nr, H5); - - DFElt C22Min, C22Max; - DFElt C12Min, C12Max; + + DFElt C22Min=0, C22Max=0; + DFElt C12Min=0, C12Max=0; // This test will be optimized out - if (Protected::NeedDoublePreAddReduction (C12Min, C12Max, H5.Outmin, H5.Outmax, WH.Cmin, WH.Cmax, betadf, WH)){ - freduce(F,mr,nr,X1,nr); - H5.initOut(); + switch (Protected::NeedDoublePreAddReduction (C12Min, C12Max, H5.Outmin, H5.Outmax, WH.Cmin, WH.Cmax, betadf, WH)){ + case 2: + freduce(F,mr,nr,X1,nr); + H5.initOut(); + C22Min = C12Min; C22Max = C12Max; + // C22 = P5 + beta C22 in C22 + fadd(F,mr,nr,X1,nr,beta,C22,ldc,C22,ldc); + + // C12 = P5 + beta C12 in C12 + fadd(F,mr,nr,X1,nr,beta,C12,ldc,C12,ldc); + break; + case 1: + freduce(F,mr,nr,X1,nr); + H5.initOut(); + case 0: + C22Min = C12Min; C22Max = C12Max; + // C22 = P5 + beta C22 in C22 + fadd(DF,mr,nr,(DFCEptr)X1,nr,betadf,(DFCEptr)C22,ldc,(DFEptr)C22,ldc); + + // C12 = P5 + beta C12 in C12 + fadd(DF,mr,nr,(DFCEptr)X1,nr,betadf,(DFCEptr)C12,ldc,(DFEptr)C12,ldc); + break; } - C22Min = C12Min; C22Max = C12Max; - - // C22 = P5 + beta C22 in C22 - fadd(DF,mr,nr,(DFCEptr)X1,nr,betadf,(DFCEptr)C22,ldc,(DFEptr)C22,ldc); - - // C12 = P5 + beta C12 in C12 - fadd(DF,mr,nr,(DFCEptr)X1,nr,betadf,(DFCEptr)C12,ldc,(DFEptr)C12,ldc); // P1 = alpha . A11 * B11 in X1 MMH_t H1(F, WH.recLevel-1, @@ -388,7 +399,144 @@ namespace FFLAS { namespace BLAS3 { WH.Outmin = min4 (U1Min, H3.Outmin, U6Min, U7Min); WH.Outmax = max4 (U1Max, H3.Outmax, U6Max, U7Max); } // WinogradAcc + + // 3 temps and 21 ops + // Specialization other fields not requiring management of bounds + template < class Field, class FieldTrait> + inline typename std::enable_if::value, void>::type + WinogradAcc_3_21 (const Field& F, + const FFLAS_TRANSPOSE ta, + const FFLAS_TRANSPOSE tb, + const size_t mr, const size_t nr, const size_t kr, + const typename Field::Element alpha, + typename Field::ConstElement_ptr A,const size_t lda, + typename Field::ConstElement_ptr B,const size_t ldb, + const typename Field::Element beta, + typename Field::Element_ptr C, const size_t ldc, + MMHelper & WH) + { + typedef MMHelper MMH_t; + + FFLASFFPACK_check(!F.isZero(beta)); + MMH_t H = WH ; + H.recLevel--; + + typename Field::Element mbeta; + F.neg(mbeta,beta); + + size_t lb, cb, la, ca; + size_t x3rd = std::max(mr,kr); + typename Field::ConstElement_ptr A11=A, A12, A21, A22; + typename Field::ConstElement_ptr B11=B, B12, B21, B22; + typename Field::Element_ptr C11=C, C12=C+nr, C21=C+mr*ldc, C22=C21+nr; + + size_t ldX3; + + if (ta == FflasTrans) { + A21 = A + mr; + A12 = A + kr*lda; + A22 = A12 + mr; + la = kr; + ca = mr; + } else { // ta == FflasNoTrans + A12 = A + kr; + A21 = A + mr*lda; + A22 = A21 + kr; + la = mr; + ca = kr; + } + if (tb == FflasTrans) { + B21 = B + kr; + B12 = B + nr*ldb; + B22 = B12 + kr; + lb = nr; + cb = kr; + ldX3 = x3rd; + } else { // ta == FflasNoTrans + B12 = B + nr; + B21 = B + kr*ldb; + B22 = B21 + nr; + lb = kr; + ldX3 = cb = nr; + } + + // Three temporary submatrices are required + typename Field::Element_ptr X3 = fflas_new (F, x3rd, nr); + + // T1 = B12 - B11 in X3 + fsub(F,lb,cb,B12,ldb,B11,ldb,X3,ldX3); + + typename Field::Element_ptr X2 = fflas_new(F,mr,kr); + + // S1 = A21 + A22 in X2 + fadd(F,la,ca,A21,lda,A22,lda,X2,ca); + + typename Field::Element_ptr X1 = fflas_new(F,mr,nr); + // P5 = alpha . S1*T1 in X1 + fgemm (F, ta, tb, mr, nr, kr, alpha, X2, ca, X3, ldX3, F.zero, X1, nr, H); + + // C22 = P5 + beta C22 in C22 + fadd(F,mr,nr,X1,nr,beta,C22,ldc,C22,ldc); + + // C12 = P5 + beta C12 in C12 + fadd(F,mr,nr,X1,nr,beta,C12,ldc,C12,ldc); + + // P1 = alpha . A11 * B11 in X1 + fgemm (F, ta, tb, mr, nr, kr, alpha, A11, lda, B11, ldb, F.zero, X1, nr, H); + + // P2 = alpha . A12 * B21 + beta . C11 in C11 + fgemm (F, ta, tb, mr, nr, kr, alpha, A12, lda, B21, ldb, beta, C11, ldc, H); + + // U1 = P2 + P1 in C11 + faddin(F,mr,nr,X1,nr,C11,ldc); + + // T2 = B22 - T1 in X3 + fsub(F,lb,cb,B22,ldb,X3,ldX3,X3,ldX3); + + // S2 = S1 - A11 in X2 + fsubin(F,la,ca,A11,lda,X2,ca); + + // U2 = P6 + P1 = alpha . S2 * T2 + P1 in X1 + fgemm (F, ta, tb, mr, nr, kr, alpha, X2, ca, X3, ldX3, F.one, X1, nr, H); + + // U4 = U2 + C12 in C12 + faddin(F,mr,nr,X1,nr,C12,ldc); + + // T4 = T2 - B21 in X3 + fsubin(F,lb,cb,B21,ldb,X3,ldX3); + + // S4 = A12 -S2 in X2 + fsub(F,la,ca,A12,lda,X2,ca,X2,ca); + + // P4 = alpha . A22 * T4 - beta . C21 in C21 + fgemm (F, ta, tb, mr, nr, kr, alpha, A22, lda, X3, ldX3, mbeta, C21, ldc, H); + + // U5 = P3 + U4 = alpha . S4*B22 + U4 in C12 + fgemm (F, ta, tb, mr, nr, kr, alpha, X2, ca, B22, ldb, F.one, C12, ldc, H); + + // T3 = B22 - B12 in X3 + fsub(F,lb,cb,B22,ldb,B12,ldb,X3,ldX3); + + // S3 = A11 - A21 in X2 + fsub(F,la,ca,A11,lda,A21,lda,X2,ca); + + // U3 = P7 + U2 = alpha . S3 * T3 + U2 in X1 + fgemm (F, ta, tb, mr, nr, kr, alpha, X2, ca, X3, ldX3, F.one, X1, nr, H); + + fflas_delete (X2); + fflas_delete (X3); + + // U7 = U3 + C22 in C22 + faddin(F,mr,nr,X1,nr,C22,ldc); + + // U6 = U3 - P4 in C21 + fsub(F,mr,nr,X1,nr,C21,ldc,C21,ldc); + + fflas_delete (X1); + + } // WinogradAcc + // 2 temps and 24 ops // TODO: Add check for modular reductions before final additions diff --git a/fflas-ffpack/fflas/fflas_helpers.inl b/fflas-ffpack/fflas/fflas_helpers.inl index 85478d4d2..f48aa7f1d 100644 --- a/fflas-ffpack/fflas/fflas_helpers.inl +++ b/fflas-ffpack/fflas/fflas_helpers.inl @@ -100,6 +100,7 @@ namespace FFLAS { struct MMHelper { typedef MMHelper Self_t; + typedef ModeCategories::DefaultTag Mode_t; int recLevel ; ParSeqTrait parseq; @@ -127,6 +128,7 @@ namespace FFLAS { struct MMHelper, ParSeqTrait> { typedef MMHelper,ParSeqTrait> Self_t; + typedef ModeCategories::ConvertTo Mode_t; int recLevel ; ParSeqTrait parseq; @@ -154,6 +156,7 @@ namespace FFLAS { typename ParSeqTrait> struct MMHelper { typedef MMHelper Self_t; + typedef ModeTrait Mode_t; typedef typename associatedDelayedField::type DelayedField_t; typedef typename associatedDelayedField::field DelayedField; typedef typename DelayedField::Element DFElt; @@ -360,21 +363,52 @@ namespace FFLAS { // to be used in the future, when Winograd's algorithm will be made generic wrt the ModeTrait - // template - // void copyOutBounds(const MMHelper &Source, - // MMHelper & Dest){ - // Dest.Outmax = Source.Outmax; - // Dest.Outmin = Source.Outmin; - // } - // template - // void copyOutBounds(const MMHelper &Source, - // MMHelper & Dest){ - // Dest.Outmax = Source.Outmax; - // Dest.Outmin = Source.Outmin; - // } - // template - // void copyOutBounds(const MMH1 &Source, MMH2 & Dest){} - /*! StructureHelper for ftrsm + template + inline typename std::enable_if::value && + FFLAS::hasBounds::value, void>::type + copyOutBounds(const MMH1& Source, MMH2& Dest){ + Dest.Outmax = Source.Outmax; + Dest.Outmin = Source.Outmin; + } + template + inline typename std::enable_if::value || + !FFLAS::hasBounds::value, void>::type + copyOutBounds(const MMH1& Source, MMH2& Dest){} + + template + inline typename std::enable_if::value && + FFLAS::hasBounds::value && + FFLAS::hasBounds::value && + FFLAS::hasBounds::value, void>::type + mergeOutBounds (const MMH1& H1, const MMH2& H2, const MMH3& H3, MMH4& Dest){ + Dest.Outmax = max4 (H1.Outmax, H2.Outmax, H3.Outmax, Dest.Outmax); + Dest.Outmin = min4 (H1.Outmin, H2.Outmin, H3.Outmin, Dest.Outmin); + } + template + inline typename std::enable_if::value || + !FFLAS::hasBounds::value || + !FFLAS::hasBounds::value || + !FFLAS::hasBounds::value, void>::type + mergeOutBounds (const MMH1& H1, const MMH2& H2, const MMH3& H3, MMH4& Dest){} + + template + inline typename std::enable_if::value && + FFLAS::hasBounds::value, void>::type + copyAccumulator (const bool fromOut, const MMH1& Source, MMH2& Dest){ + if (fromOut){ + Dest.Cmin = Source.Outmin; + Dest.Cmax = Source.Outmax; + } else { + Dest.Cmin = Source.Cmin; + Dest.Cmax = Source.Cmax; + } + } + template + inline typename std::enable_if::value || + !FFLAS::hasBounds::value, void>::type + copyAccumulator (const bool fromOut, const MMH1& Source, MMH2& Dest) {} + + /*! StructureHelper for ftrsm */ namespace StructureHelper { struct Recursive{}; diff --git a/fflas-ffpack/fflas/fflas_level2.inl b/fflas-ffpack/fflas/fflas_level2.inl index 791005012..8d753fe56 100644 --- a/fflas-ffpack/fflas/fflas_level2.inl +++ b/fflas-ffpack/fflas/fflas_level2.inl @@ -362,14 +362,14 @@ namespace FFLAS { typename Field::ConstElement_ptr A, const size_t lda, typename Field::Element_ptr B, const size_t ldb); - /** \brief faxpy : \f$y \gets \alpha \cdot x + y\f$. + /** \brief faxpy : \f$Y \gets \alpha \cdot X + Y\f$. * @param F field * @param m row dimension * @param n column dimension * @param alpha scalar - * \param[in] X vector in \p F + * \param[in] X matrix in \p F * \param ldx leading dimension of \p X - * \param[in,out] Y vector in \p F + * \param[in,out] Y matrix in \p F * \param ldy leading dimension of \p Y */ template diff --git a/fflas-ffpack/field/field-traits.h b/fflas-ffpack/field/field-traits.h index a713978d8..7977ce4ef 100644 --- a/fflas-ffpack/field/field-traits.h +++ b/fflas-ffpack/field/field-traits.h @@ -39,6 +39,7 @@ #include "recint/rmint.h" #include "givaro/modular-general.h" +#include "givaro/modular-extended.h" #include "givaro/zring.h" namespace RecInt { @@ -163,13 +164,24 @@ namespace FFLAS { /* Traits */ /*! ModeTraits - */ + */ + template + class hasBounds : public std::false_type{}; + + template <> + class hasBounds : public std::true_type{}; + template <> + class hasBounds : public std::true_type{}; + template <> + class hasBounds : public std::true_type{}; + template struct ModeTraits {typedef typename ModeCategories::DefaultTag value;}; template struct ModeTraits >{typedef typename ModeCategories::DelayedTag value;}; template<> struct ModeTraits > {typedef typename ModeCategories::DefaultTag value;}; + template<> struct ModeTraits > {typedef typename ModeCategories::DelayedTag value;}; template struct ModeTraits > {typedef typename ModeCategories::ConvertTo value;}; template struct ModeTraits > {typedef typename ModeCategories::ConvertTo value;}; @@ -327,6 +339,11 @@ namespace FFLAS { /* associatedDelayedField */ typedef Givaro::ZRing type; }; template + struct associatedDelayedField> { + typedef Givaro::ZRing field; + typedef Givaro::ZRing type; + }; + template struct associatedDelayedField> { typedef Givaro::ZRing field; typedef Givaro::ZRing type; diff --git a/tests/test-fgemm.C b/tests/test-fgemm.C index 443e51d5f..67a9a5d7f 100644 --- a/tests/test-fgemm.C +++ b/tests/test-fgemm.C @@ -41,7 +41,7 @@ #include #include -#include +#include #include @@ -59,6 +59,7 @@ using namespace FFPACK; using Givaro::Modular; using Givaro::ModularBalanced; +using Givaro::ModularExtended; // checks that D = alpha . C + beta . A ^ta * B ^tb @@ -85,7 +86,6 @@ bool check_MM(const Field & F, ConstElement_ptr ail,blj; Element_ptr D = fflas_new (F,m,n); fassign(F,m,n,Cd,n,D,n); - for (size_t i = 0; i < m; ++i) for (size_t j = 0; j < n; ++j){ F.mulin(*(D+i*n+j),beta); @@ -243,8 +243,8 @@ bool launch_MM_dispatch(const Field &F, FFLAS_TRANSPOSE ta = FflasNoTrans ; FFLAS_TRANSPOSE tb = FflasNoTrans ; if (! par) { - if (random()%2) ta = FflasTrans ; - if (random()%2) tb = FflasTrans ; + // if (random()%2) ta = FflasTrans ; + // if (random()%2) tb = FflasTrans ; } if (mm<0) @@ -313,21 +313,21 @@ bool run_with_field (Givaro::Integer q, uint64_t b, int m, int n, int k, int nbw ok = ok && launch_MM_dispatch(*F,m,n,k,F->one,F->zero,iters,nbw, par, R); //std::cout << k << "/24" << std::endl; ++k; ok = ok && launch_MM_dispatch(*F,m,n,k,F->zero,F->zero,iters,nbw, par, R); - //std::cout << k << "/24" << std::endl; ++k; + // //std::cout << k << "/24" << std::endl; ++k; ok = ok && launch_MM_dispatch(*F,m,n,k,F->mOne,F->zero,iters,nbw, par, R); - //std::cout << k << "/24" << std::endl; ++k; + // //std::cout << k << "/24" << std::endl; ++k; ok = ok && launch_MM_dispatch(*F,m,n,k,F->one ,F->one,iters,nbw, par, R); - //std::cout << k << "/24" << std::endl; ++k; + // //std::cout << k << "/24" << std::endl; ++k; ok = ok && launch_MM_dispatch(*F,m,n,k,F->zero,F->one,iters,nbw, par, R); - //std::cout << k << "/24" << std::endl; ++k; + // //std::cout << k << "/24" << std::endl; ++k; ok = ok && launch_MM_dispatch(*F,m,n,k,F->mOne,F->one,iters,nbw, par, R); - //std::cout << k << "/24" << std::endl; ++k; + // //std::cout << k << "/24" << std::endl; ++k; ok = ok && launch_MM_dispatch(*F,m,n,k,F->one ,F->mOne,iters,nbw, par, R); - //std::cout << k << "/24" << std::endl; ++k; + // //std::cout << k << "/24" << std::endl; ++k; ok = ok && launch_MM_dispatch(*F,m,n,k,F->zero,F->mOne,iters,nbw, par, R); - //std::cout << k << "/24" << std::endl; ++k; + // //std::cout << k << "/24" << std::endl; ++k; ok = ok && launch_MM_dispatch(*F,m,n,k,F->mOne,F->mOne,iters,nbw, par, R); - //std::cout << k << "/24" << std::endl; ++k; + // //std::cout << k << "/24" << std::endl; ++k; Element alpha,beta ; NZR.random(alpha); @@ -337,10 +337,10 @@ bool run_with_field (Givaro::Integer q, uint64_t b, int m, int n, int k, int nbw //std::cout << k << "/24" << std::endl; ++k; ok = ok && launch_MM_dispatch(*F,m,n,k,F->mOne,alpha,iters,nbw, par, R); //std::cout << k << "/24" << std::endl; ++k; - ok = ok && launch_MM_dispatch(*F,m,n,k,alpha,F->one ,iters,nbw, par, R); + ok = ok && launch_MM_dispatch(*F,m,n,k,alpha,F->one ,iters,nbw, par, R); //std::cout << k << "/24" << std::endl; ++k; ok = ok && launch_MM_dispatch(*F,m,n,k,alpha,F->zero,iters,nbw, par, R); - //std::cout << k << "/24" << std::endl; ++k; + //std::cout << k << "/24" << std::endl; ++k; ok = ok && launch_MM_dispatch(*F,m,n,k,alpha,F->mOne,iters,nbw, par, R); //std::cout << k << "/24" << std::endl; ++k; @@ -404,10 +404,14 @@ int main(int argc, char** argv) ok = ok && run_with_field >(q,b,m,n,k,nbw,iters,p, seed); ok = ok && run_with_field >(q,b,m,n,k,nbw,iters,p, seed); #endif + ok = ok && run_with_field >(q,b,m,n,k,nbw,iters, p, seed); ok = ok && run_with_field >(q,b?b:25,m,n,k,nbw,iters, p, seed); ok = ok && run_with_field >(q,b,m,n,k,nbw,iters, p, seed); ok = ok && run_with_field >(q,b?b:25,m,n,k,nbw,iters, p, seed); + + ok = ok && run_with_field >(q,b?b:45,m,n,k,nbw,iters, p, seed); + ok = ok && run_with_field > >(q,b?b:63_ui64,m,n,k,nbw,iters, p, seed); ok = ok && run_with_field > >(q,b?b:63_ui64,m,n,k,nbw,iters, p, seed); ok = ok && run_with_field > >(q,b?b:127_ui64,m,n,k,nbw,iters, p, seed);