From 6d171624ca6125677f7df2cfa4bf94454e9b71d1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Cl=C3=A9ment=20Pernet?= Date: Tue, 27 Jun 2023 15:34:03 +0200 Subject: [PATCH 1/9] enable Strassen-Winograd for non delayed fields WIP --- .../fflas/fflas_fgemm/fgemm_winograd.inl | 58 ++++--- .../fflas/fflas_fgemm/schedule_winograd.inl | 139 ++++++++++++++- .../fflas_fgemm/schedule_winograd_acc.inl | 158 ++++++++++++++++-- fflas-ffpack/fflas/fflas_helpers.inl | 67 ++++++-- fflas-ffpack/field/field-traits.h | 12 +- tests/test-fgemm.C | 44 ++--- 6 files changed, 405 insertions(+), 73 deletions(-) diff --git a/fflas-ffpack/fflas/fflas_fgemm/fgemm_winograd.inl b/fflas-ffpack/fflas/fflas_fgemm/fgemm_winograd.inl index 37a794afd..0573756f5 100644 --- a/fflas-ffpack/fflas/fflas_fgemm/fgemm_winograd.inl +++ b/fflas-ffpack/fflas/fflas_fgemm/fgemm_winograd.inl @@ -197,9 +197,9 @@ namespace FFLAS { namespace Protected { typename Field::ConstElement_ptr B, const size_t ldb, 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 & H ) + // const typename MMHelper::DelayedField::Element Cmin, + // const typename MMHelper::DelayedField::Element Cmax) { size_t mkn =(size_t)( (bool)(nr > 0)+ ((bool)(kr > 0) << 1)+ ((bool)(mr > 0) << 2)); if (mkn == 0) return; @@ -226,12 +226,15 @@ namespace FFLAS { namespace Protected { MMHelper HModd(H); MMHelper HNodd(H); - Hacc.Cmin = H.Outmin; Hacc.Cmax = H.Outmax; + //Hacc.Cmin = H.Outmin; Hacc.Cmax = H.Outmax; + copyAccumulator(H, Hacc); + 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; + + //HModd.Cmin = Cmin; HModd.Cmax = Cmax; + //HNodd.Cmin = Cmin; HNodd.Cmax = Cmax; + // HModd.Amax = H.Bmax; HModd.Amin = H.Bmin; + // HModd.Bmax = H.Amax; HModd.Bmin = H.Amin; switch (mkn) { case 1: // n oddsized @@ -270,9 +273,10 @@ 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); + // 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); } // #define NEWIP @@ -398,17 +402,18 @@ 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); @ todo should replace the following 2 lines. + // // Not ready yet, because of API of DynamicPeeling + copyOutBounds(HC,H); + // H.Outmax = HC.Outmax; + // H.Outmin = HC.Outmin; return C; } // Then w >0 - typedef typename MMHelper::DelayedField::Element DFElt; - DFElt Cmin = H.Cmin; - DFElt Cmax = H.Cmax; + // typedef typename MMHelper::DelayedField::Element DFElt; + // DFElt Cmin = H.Cmin; + // DFElt Cmax = H.Cmax; #ifdef OLD_DYNAMIC_PEELING @@ -435,7 +440,8 @@ 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, Cmin, Cmax); + Protected::DynamicPeeling2 (F, ta, tb, m, n, k, mr, nr, kr, alpha, A, lda, B, ldb, beta, C, ldc, H); // Let's see if Cmin, Cmax in H are still valid #endif return C; } // fgemm @@ -497,17 +503,18 @@ namespace FFLAS{ // 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); + // H.Outmax = HC.Outmax; + // H.Outmin = HC.Outmin; return C; } /// // Then w >0 - typedef typename MMHelper::DelayedField::Element DFElt; - DFElt Cmin = H.Cmin; - DFElt Cmax = H.Cmax; + // 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); @@ -533,7 +540,8 @@ 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, Cmin, Cmax); + Protected::DynamicPeeling2 (F, ta, tb, m, n, k, mr, nr, kr, alpha, A, lda, B, ldb, beta, C, ldc, HC); #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..99a980928 100644 --- a/fflas-ffpack/fflas/fflas_fgemm/schedule_winograd.inl +++ b/fflas-ffpack/fflas/fflas_fgemm/schedule_winograd.inl @@ -341,7 +341,8 @@ namespace FFLAS { namespace BLAS3 { template < class Field, class FieldTrait > - inline void Winograd (const Field& F, + 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, @@ -535,6 +536,142 @@ namespace FFLAS { namespace BLAS3 { } // Winograd + + 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..a9e8db7bb 100644 --- a/fflas-ffpack/fflas/fflas_fgemm/schedule_winograd_acc.inl +++ b/fflas-ffpack/fflas/fflas_fgemm/schedule_winograd_acc.inl @@ -185,17 +185,17 @@ namespace FFLAS { namespace BLAS3 { // 3 temps and 21 ops 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; @@ -388,7 +388,143 @@ 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 + 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..9394e88cc 100644 --- a/fflas-ffpack/fflas/fflas_helpers.inl +++ b/fflas-ffpack/fflas/fflas_helpers.inl @@ -360,21 +360,58 @@ 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 + 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){} + + template + void mergeOutBounds (const MMHelper &H1, + const MMHelper &H2, + const MMHelper &H3, + MMHelper & Dest){ + Dest.Outmax = max4 (H1.Outmax, H2.Outmax, H3.Outmax, Dest.Outmax); + Dest.Outmin = max4 (H1.Outmin, H2.Outmin, H3.Outmin, Dest.Outmin); + } + template + void mergeOutBounds (const MMHelper &H1, + const MMHelper &H2, + const MMHelper &H3, + MMHelper & Dest){ + Dest.Outmax = max4 (H1.Outmax, H2.Outmax, H3.Outmax, Dest.Outmax); + Dest.Outmin = max4 (H1.Outmin, H2.Outmin, H3.Outmin, Dest.Outmin); + } + template + void mergeOutBounds (const MMH1& H1, const MMH2& H2, const MMH3& H3, MMH4& Dest){} + + template + void copyAccumulator (const MMHelper& Source, + MMHelper& Dest){ + Dest.Cmin = Source.Outmin; + Dest.Cmax = Source.Outmax; + } + template + void copyAccumulator( const MMHelper& Source, + MMHelper& Dest){ + Dest.Cmin = Source.Outmin; + Dest.Cmax = Source.Outmax; + } + + template + void copyAccumulator(const MMH1& Source, MMH2& Dest){} + + //Hacc.Cmin = H.Outmin; Hacc.Cmax = H.Outmax; + /*! StructureHelper for ftrsm */ namespace StructureHelper { struct Recursive{}; diff --git a/fflas-ffpack/field/field-traits.h b/fflas-ffpack/field/field-traits.h index a713978d8..76abcf950 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,22 @@ namespace FFLAS { /* Traits */ /*! ModeTraits - */ + */ + template + class isDelayed : public std::false_type{}; + + template <> + class isDelayed : public std::true_type{}; + template <> + class isDelayed : 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::DefaultTag value;}; template struct ModeTraits > {typedef typename ModeCategories::ConvertTo value;}; template struct ModeTraits > {typedef typename ModeCategories::ConvertTo value;}; diff --git a/tests/test-fgemm.C b/tests/test-fgemm.C index 443e51d5f..409e3a19a 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 @@ -395,25 +396,28 @@ int main(int argc, char** argv) bool ok = true; srand(seed); do{ - 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); - 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); -#ifndef __FFLASFFPACK_HAVE_AVX512F_INSTRUCTIONS - // int32_t simd not yet fully implemented over AVX512 - 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: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); - ok = ok && run_with_field,RecInt::ruint<8> > >(q,b?b:127_ui64,m,n,k,nbw,iters, p, seed); - ok = ok && run_with_field >(q,(b?b:512_ui64),m,n,k,nbw,iters,p, seed); - ok = ok && run_with_field >(0,(b?b:512_ui64),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,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,m,n,k,nbw,iters,p, seed); +// #ifndef __FFLASFFPACK_HAVE_AVX512F_INSTRUCTIONS +// // int32_t simd not yet fully implemented over AVX512 +// 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); + // ok = ok && run_with_field,RecInt::ruint<8> > >(q,b?b:127_ui64,m,n,k,nbw,iters, p, seed); + // ok = ok && run_with_field >(q,(b?b:512_ui64),m,n,k,nbw,iters,p, seed); + // ok = ok && run_with_field >(0,(b?b:512_ui64),m,n,k,nbw,iters,p, seed); seed++; } while (loop && ok); From 2f79a03497e4695af2743c52359ff4fee29a57fb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Cl=C3=A9ment=20Pernet?= Date: Wed, 28 Jun 2023 09:46:49 +0200 Subject: [PATCH 2/9] done --- .../fflas/fflas_fgemm/fgemm_classical.inl | 1 + .../fflas/fflas_fgemm/fgemm_winograd.inl | 17 ++++--- fflas-ffpack/fflas/fflas_helpers.inl | 23 ++++------ tests/test-fgemm.C | 44 +++++++++---------- 4 files changed, 43 insertions(+), 42 deletions(-) diff --git a/fflas-ffpack/fflas/fflas_fgemm/fgemm_classical.inl b/fflas-ffpack/fflas/fflas_fgemm/fgemm_classical.inl index db1d79282..b50fb777d 100644 --- a/fflas-ffpack/fflas/fflas_fgemm/fgemm_classical.inl +++ b/fflas-ffpack/fflas/fflas_fgemm/fgemm_classical.inl @@ -226,6 +226,7 @@ namespace FFLAS { for (size_t l = 0; l < k; ++l) F.axpyin (*(C+i*ldc+j), *(A+l*lda+i), *(B+j*ldb+l)); fscalin(F,m,n,alpha,C,ldc); + } template < class Field> inline void fgemm (const Field& F, diff --git a/fflas-ffpack/fflas/fflas_fgemm/fgemm_winograd.inl b/fflas-ffpack/fflas/fflas_fgemm/fgemm_winograd.inl index 0573756f5..8941e36dc 100644 --- a/fflas-ffpack/fflas/fflas_fgemm/fgemm_winograd.inl +++ b/fflas-ffpack/fflas/fflas_fgemm/fgemm_winograd.inl @@ -197,7 +197,8 @@ namespace FFLAS { namespace Protected { typename Field::ConstElement_ptr B, const size_t ldb, const typename Field::Element beta, typename Field::Element_ptr C, const size_t ldc, - MMHelper & H ) + MMHelper & H, + MMHelper & Hsave) // const typename MMHelper::DelayedField::Element Cmin, // const typename MMHelper::DelayedField::Element Cmax) { @@ -228,11 +229,11 @@ namespace FFLAS { namespace Protected { //Hacc.Cmin = H.Outmin; Hacc.Cmax = H.Outmax; copyAccumulator(H, Hacc); - + Hacc.recLevel=-1;HModd.recLevel=-1;HNodd.recLevel=-1; - //HModd.Cmin = Cmin; HModd.Cmax = Cmax; - //HNodd.Cmin = Cmin; HNodd.Cmax = Cmax; + HModd.Cmin = Hsave.Cmin; HModd.Cmax = Hsave.Cmax; + HNodd.Cmin = Hsave.Cmin; HNodd.Cmax = Hsave.Cmax; // HModd.Amax = H.Bmax; HModd.Amin = H.Bmin; // HModd.Bmax = H.Amax; HModd.Bmin = H.Amin; @@ -430,6 +431,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; @@ -441,7 +444,7 @@ namespace FFLAS{ 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); // Let's see if Cmin, Cmax in H are still valid + Protected::DynamicPeeling2 (F, ta, tb, m, n, k, mr, nr, kr, alpha, A, lda, B, ldb, beta, C, ldc, H, Hsave); // Let's see if Cmin, Cmax in Hsave are still valid #endif return C; } // fgemm @@ -530,6 +533,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; @@ -541,7 +546,7 @@ namespace FFLAS{ 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); + 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_helpers.inl b/fflas-ffpack/fflas/fflas_helpers.inl index 9394e88cc..b177d0527 100644 --- a/fflas-ffpack/fflas/fflas_helpers.inl +++ b/fflas-ffpack/fflas/fflas_helpers.inl @@ -360,19 +360,14 @@ 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 + typename std::enable_if::value && FFLAS::isDelayed::value, void>::type + copyOutBounds(const MMHelper &Source, + MMHelper & Dest){ + Dest.Outmax = Source.Outmax; + Dest.Outmin = Source.Outmin; } - template + template void copyOutBounds(const MMH1 &Source, MMH2 & Dest){} template @@ -381,7 +376,7 @@ namespace FFLAS { const MMHelper &H3, MMHelper & Dest){ Dest.Outmax = max4 (H1.Outmax, H2.Outmax, H3.Outmax, Dest.Outmax); - Dest.Outmin = max4 (H1.Outmin, H2.Outmin, H3.Outmin, Dest.Outmin); + Dest.Outmin = min4 (H1.Outmin, H2.Outmin, H3.Outmin, Dest.Outmin); } template void mergeOutBounds (const MMHelper &H1, @@ -389,7 +384,7 @@ namespace FFLAS { const MMHelper &H3, MMHelper & Dest){ Dest.Outmax = max4 (H1.Outmax, H2.Outmax, H3.Outmax, Dest.Outmax); - Dest.Outmin = max4 (H1.Outmin, H2.Outmin, H3.Outmin, Dest.Outmin); + Dest.Outmin = min4 (H1.Outmin, H2.Outmin, H3.Outmin, Dest.Outmin); } template void mergeOutBounds (const MMH1& H1, const MMH2& H2, const MMH3& H3, MMH4& Dest){} diff --git a/tests/test-fgemm.C b/tests/test-fgemm.C index 409e3a19a..60e80b1b7 100644 --- a/tests/test-fgemm.C +++ b/tests/test-fgemm.C @@ -396,28 +396,28 @@ int main(int argc, char** argv) bool ok = true; srand(seed); do{ -// 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); -// 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); -// #ifndef __FFLASFFPACK_HAVE_AVX512F_INSTRUCTIONS -// // int32_t simd not yet fully implemented over AVX512 -// 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); - // ok = ok && run_with_field,RecInt::ruint<8> > >(q,b?b:127_ui64,m,n,k,nbw,iters, p, seed); - // ok = ok && run_with_field >(q,(b?b:512_ui64),m,n,k,nbw,iters,p, seed); - // ok = ok && run_with_field >(0,(b?b:512_ui64),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,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,m,n,k,nbw,iters,p, seed); +#ifndef __FFLASFFPACK_HAVE_AVX512F_INSTRUCTIONS + // int32_t simd not yet fully implemented over AVX512 + 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); + ok = ok && run_with_field,RecInt::ruint<8> > >(q,b?b:127_ui64,m,n,k,nbw,iters, p, seed); + ok = ok && run_with_field >(q,(b?b:512_ui64),m,n,k,nbw,iters,p, seed); + ok = ok && run_with_field >(0,(b?b:512_ui64),m,n,k,nbw,iters,p, seed); seed++; } while (loop && ok); From b4929357aae12f2dc5acc12f32976806ce602bed Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Cl=C3=A9ment=20Pernet?= Date: Wed, 28 Jun 2023 10:40:30 +0200 Subject: [PATCH 3/9] cleaning up comments --- .../fflas/fflas_fgemm/fgemm_winograd.inl | 22 +--------------- .../fflas/fflas_fgemm/schedule_winograd.inl | 26 +++++++++---------- .../fflas_fgemm/schedule_winograd_acc.inl | 2 ++ 3 files changed, 16 insertions(+), 34 deletions(-) diff --git a/fflas-ffpack/fflas/fflas_fgemm/fgemm_winograd.inl b/fflas-ffpack/fflas/fflas_fgemm/fgemm_winograd.inl index 8941e36dc..67725cb1c 100644 --- a/fflas-ffpack/fflas/fflas_fgemm/fgemm_winograd.inl +++ b/fflas-ffpack/fflas/fflas_fgemm/fgemm_winograd.inl @@ -199,8 +199,6 @@ namespace FFLAS { namespace Protected { typename Field::Element_ptr C, const size_t ldc, MMHelper & H, MMHelper & Hsave) - // const typename MMHelper::DelayedField::Element Cmin, - // const typename MMHelper::DelayedField::Element Cmax) { size_t mkn =(size_t)( (bool)(nr > 0)+ ((bool)(kr > 0) << 1)+ ((bool)(mr > 0) << 2)); if (mkn == 0) return; @@ -227,15 +225,12 @@ namespace FFLAS { namespace Protected { MMHelper HModd(H); MMHelper HNodd(H); - //Hacc.Cmin = H.Outmin; Hacc.Cmax = H.Outmax; copyAccumulator(H, Hacc); - Hacc.recLevel=-1;HModd.recLevel=-1;HNodd.recLevel=-1; + Hacc.recLevel=-1; HModd.recLevel=-1; HNodd.recLevel=-1; HModd.Cmin = Hsave.Cmin; HModd.Cmax = Hsave.Cmax; HNodd.Cmin = Hsave.Cmin; HNodd.Cmax = Hsave.Cmax; - // HModd.Amax = H.Bmax; HModd.Amin = H.Bmin; - // HModd.Bmax = H.Amax; HModd.Bmin = H.Amin; switch (mkn) { case 1: // n oddsized @@ -275,9 +270,6 @@ namespace FFLAS { namespace Protected { break; } mergeOutBounds (HModd, HNodd, Hacc, H); - // 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); } // #define NEWIP @@ -403,18 +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 copyOutBounds(HC,H); - // H.Outmax = HC.Outmax; - // H.Outmin = HC.Outmin; return C; } // Then w >0 - // typedef typename MMHelper::DelayedField::Element DFElt; - // DFElt Cmin = H.Cmin; - // DFElt Cmax = H.Cmax; #ifdef OLD_DYNAMIC_PEELING @@ -443,7 +428,6 @@ 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); // Let's see if Cmin, Cmax in Hsave are still valid #endif return C; @@ -503,12 +487,9 @@ 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); copyOutBounds(H,HC); - // H.Outmax = HC.Outmax; - // H.Outmin = HC.Outmin; return C; } @@ -545,7 +526,6 @@ 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; diff --git a/fflas-ffpack/fflas/fflas_fgemm/schedule_winograd.inl b/fflas-ffpack/fflas/fflas_fgemm/schedule_winograd.inl index 99a980928..9c6b7815d 100644 --- a/fflas-ffpack/fflas/fflas_fgemm/schedule_winograd.inl +++ b/fflas-ffpack/fflas/fflas_fgemm/schedule_winograd.inl @@ -339,21 +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 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 - ) + 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)); @@ -537,6 +536,7 @@ 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, diff --git a/fflas-ffpack/fflas/fflas_fgemm/schedule_winograd_acc.inl b/fflas-ffpack/fflas/fflas_fgemm/schedule_winograd_acc.inl index a9e8db7bb..f30614352 100644 --- a/fflas-ffpack/fflas/fflas_fgemm/schedule_winograd_acc.inl +++ b/fflas-ffpack/fflas/fflas_fgemm/schedule_winograd_acc.inl @@ -184,6 +184,7 @@ 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 typename std::enable_if::value, void>::type WinogradAcc_3_21 (const Field& F, @@ -390,6 +391,7 @@ namespace FFLAS { namespace BLAS3 { } // 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, From 4eea199576d373a4da3ab821c785798d398107d3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Cl=C3=A9ment=20Pernet?= Date: Wed, 28 Jun 2023 10:41:05 +0200 Subject: [PATCH 4/9] cleaning up comments --- fflas-ffpack/fflas/fflas_fgemm/fgemm_classical.inl | 1 - 1 file changed, 1 deletion(-) diff --git a/fflas-ffpack/fflas/fflas_fgemm/fgemm_classical.inl b/fflas-ffpack/fflas/fflas_fgemm/fgemm_classical.inl index b50fb777d..db1d79282 100644 --- a/fflas-ffpack/fflas/fflas_fgemm/fgemm_classical.inl +++ b/fflas-ffpack/fflas/fflas_fgemm/fgemm_classical.inl @@ -226,7 +226,6 @@ namespace FFLAS { for (size_t l = 0; l < k; ++l) F.axpyin (*(C+i*ldc+j), *(A+l*lda+i), *(B+j*ldb+l)); fscalin(F,m,n,alpha,C,ldc); - } template < class Field> inline void fgemm (const Field& F, From 9b498defa0b026e58bdf8429a3bf173c36db8814 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Cl=C3=A9ment=20Pernet?= Date: Wed, 28 Jun 2023 10:43:41 +0200 Subject: [PATCH 5/9] cleaning up comments --- fflas-ffpack/fflas/fflas_fgemm/fgemm_winograd.inl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/fflas-ffpack/fflas/fflas_fgemm/fgemm_winograd.inl b/fflas-ffpack/fflas/fflas_fgemm/fgemm_winograd.inl index 67725cb1c..64bd8ecba 100644 --- a/fflas-ffpack/fflas/fflas_fgemm/fgemm_winograd.inl +++ b/fflas-ffpack/fflas/fflas_fgemm/fgemm_winograd.inl @@ -428,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, Hsave); // Let's see if Cmin, Cmax in Hsave are still valid + 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 From f1188c86bebc51ab65018493ac973ec1969ccb80 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Cl=C3=A9ment=20Pernet?= Date: Wed, 28 Jun 2023 11:49:51 +0200 Subject: [PATCH 6/9] isDelayed->hasBounds is a better name --- fflas-ffpack/fflas/fflas_fgemm/fgemm_winograd.inl | 3 --- .../fflas/fflas_fgemm/schedule_winograd.inl | 4 ++-- .../fflas/fflas_fgemm/schedule_winograd_acc.inl | 4 ++-- fflas-ffpack/fflas/fflas_helpers.inl | 13 ++++++++----- fflas-ffpack/field/field-traits.h | 8 +++++--- 5 files changed, 17 insertions(+), 15 deletions(-) diff --git a/fflas-ffpack/fflas/fflas_fgemm/fgemm_winograd.inl b/fflas-ffpack/fflas/fflas_fgemm/fgemm_winograd.inl index 64bd8ecba..8829d0514 100644 --- a/fflas-ffpack/fflas/fflas_fgemm/fgemm_winograd.inl +++ b/fflas-ffpack/fflas/fflas_fgemm/fgemm_winograd.inl @@ -496,9 +496,6 @@ namespace FFLAS{ /// // 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); diff --git a/fflas-ffpack/fflas/fflas_fgemm/schedule_winograd.inl b/fflas-ffpack/fflas/fflas_fgemm/schedule_winograd.inl index 9c6b7815d..24909a0da 100644 --- a/fflas-ffpack/fflas/fflas_fgemm/schedule_winograd.inl +++ b/fflas-ffpack/fflas/fflas_fgemm/schedule_winograd.inl @@ -341,7 +341,7 @@ namespace FFLAS { namespace BLAS3 { // Specialization for Delayed or Lazy fields, requiring management of bounds template < class Field, class FieldTrait > - inline typename std::enable_if::value, void>::type + inline typename std::enable_if::value, void>::type Winograd (const Field& F, const FFLAS_TRANSPOSE ta, const FFLAS_TRANSPOSE tb, @@ -538,7 +538,7 @@ namespace FFLAS { namespace BLAS3 { // Case for fields not requiring management of bounds template < class Field, class FieldTrait > - inline typename std::enable_if::value, void>::type + inline typename std::enable_if::value, void>::type Winograd (const Field& F, const FFLAS_TRANSPOSE ta, const FFLAS_TRANSPOSE tb, diff --git a/fflas-ffpack/fflas/fflas_fgemm/schedule_winograd_acc.inl b/fflas-ffpack/fflas/fflas_fgemm/schedule_winograd_acc.inl index f30614352..2aada96b4 100644 --- a/fflas-ffpack/fflas/fflas_fgemm/schedule_winograd_acc.inl +++ b/fflas-ffpack/fflas/fflas_fgemm/schedule_winograd_acc.inl @@ -186,7 +186,7 @@ namespace FFLAS { namespace BLAS3 { // 3 temps and 21 ops // Specialization for Delayed or Lazy fields, requiring management of bounds template < class Field, class FieldTrait> - inline typename std::enable_if::value, void>::type + inline typename std::enable_if::value, void>::type WinogradAcc_3_21 (const Field& F, const FFLAS_TRANSPOSE ta, const FFLAS_TRANSPOSE tb, @@ -393,7 +393,7 @@ namespace FFLAS { namespace BLAS3 { // 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 + inline typename std::enable_if::value, void>::type WinogradAcc_3_21 (const Field& F, const FFLAS_TRANSPOSE ta, const FFLAS_TRANSPOSE tb, diff --git a/fflas-ffpack/fflas/fflas_helpers.inl b/fflas-ffpack/fflas/fflas_helpers.inl index b177d0527..8c9392040 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,15 +363,15 @@ namespace FFLAS { // to be used in the future, when Winograd's algorithm will be made generic wrt the ModeTrait - template - typename std::enable_if::value && FFLAS::isDelayed::value, void>::type - copyOutBounds(const MMHelper &Source, - MMHelper & Dest){ + template + 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 - void copyOutBounds(const MMH1 &Source, MMH2 & Dest){} + typename std::enable_if::value && FFLAS::hasBounds::value, void>::type + copyOutBounds(const MMH1& Source, MMH2& Dest){} template void mergeOutBounds (const MMHelper &H1, diff --git a/fflas-ffpack/field/field-traits.h b/fflas-ffpack/field/field-traits.h index 76abcf950..611428b4a 100644 --- a/fflas-ffpack/field/field-traits.h +++ b/fflas-ffpack/field/field-traits.h @@ -166,12 +166,14 @@ namespace FFLAS { /* Traits */ /*! ModeTraits */ template - class isDelayed : public std::false_type{}; + class hasBounds : public std::false_type{}; template <> - class isDelayed : public std::true_type{}; + class hasBounds : public std::true_type{}; template <> - class isDelayed : public std::true_type{}; + class hasBounds : public std::true_type{}; + template <> + class hasBounds : public std::true_type{}; template struct ModeTraits {typedef typename ModeCategories::DefaultTag value;}; From 2707919b556a12459ca950cec6a146baa4dfa25a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Cl=C3=A9ment=20Pernet?= Date: Wed, 28 Jun 2023 14:11:33 +0200 Subject: [PATCH 7/9] fixes and factorisation of copyAcc --- .../fflas/fflas_fgemm/fgemm_winograd.inl | 6 ++-- fflas-ffpack/fflas/fflas_helpers.inl | 36 ++++++++++--------- tests/test-fgemm.C | 2 +- 3 files changed, 23 insertions(+), 21 deletions(-) diff --git a/fflas-ffpack/fflas/fflas_fgemm/fgemm_winograd.inl b/fflas-ffpack/fflas/fflas_fgemm/fgemm_winograd.inl index 8829d0514..b69a42dd9 100644 --- a/fflas-ffpack/fflas/fflas_fgemm/fgemm_winograd.inl +++ b/fflas-ffpack/fflas/fflas_fgemm/fgemm_winograd.inl @@ -225,12 +225,12 @@ namespace FFLAS { namespace Protected { MMHelper HModd(H); MMHelper HNodd(H); - copyAccumulator(H, Hacc); + copyAccumulator (true, H, Hacc); // H.Out -> Hacc.C Hacc.recLevel=-1; HModd.recLevel=-1; HNodd.recLevel=-1; - HModd.Cmin = Hsave.Cmin; HModd.Cmax = Hsave.Cmax; - HNodd.Cmin = Hsave.Cmin; HNodd.Cmax = Hsave.Cmax; + copyAccumulator (false, Hsave, HModd); // Hsave.C -> HModd.C + copyAccumulator (false, Hsave, HNodd); // Hsave.C -> HNodd.C switch (mkn) { case 1: // n oddsized diff --git a/fflas-ffpack/fflas/fflas_helpers.inl b/fflas-ffpack/fflas/fflas_helpers.inl index 8c9392040..d9b8a7b05 100644 --- a/fflas-ffpack/fflas/fflas_helpers.inl +++ b/fflas-ffpack/fflas/fflas_helpers.inl @@ -364,13 +364,15 @@ namespace FFLAS { // to be used in the future, when Winograd's algorithm will be made generic wrt the ModeTrait template - typename std::enable_if::value && FFLAS::hasBounds::value, void>::type + 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 - typename std::enable_if::value && FFLAS::hasBounds::value, void>::type + typename std::enable_if::value || + !FFLAS::hasBounds::value, void>::type copyOutBounds(const MMH1& Source, MMH2& Dest){} template @@ -392,23 +394,23 @@ namespace FFLAS { template void mergeOutBounds (const MMH1& H1, const MMH2& H2, const MMH3& H3, MMH4& Dest){} - template - void copyAccumulator (const MMHelper& Source, - MMHelper& Dest){ - Dest.Cmin = Source.Outmin; - Dest.Cmax = Source.Outmax; - } - template - void copyAccumulator( const MMHelper& Source, - MMHelper& Dest){ - Dest.Cmin = Source.Outmin; - Dest.Cmax = Source.Outmax; + 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) {} - template - void copyAccumulator(const MMH1& Source, MMH2& Dest){} - - //Hacc.Cmin = H.Outmin; Hacc.Cmax = H.Outmax; /*! StructureHelper for ftrsm */ namespace StructureHelper { diff --git a/tests/test-fgemm.C b/tests/test-fgemm.C index 60e80b1b7..6c5377a93 100644 --- a/tests/test-fgemm.C +++ b/tests/test-fgemm.C @@ -410,7 +410,7 @@ 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?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: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); From 8313f132d4a7bd676048e31fcb00cd0831df773e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Cl=C3=A9ment=20Pernet?= Date: Wed, 28 Jun 2023 14:33:57 +0200 Subject: [PATCH 8/9] factorise code --- fflas-ffpack/fflas/fflas_helpers.inl | 33 +++++++++++++--------------- 1 file changed, 15 insertions(+), 18 deletions(-) diff --git a/fflas-ffpack/fflas/fflas_helpers.inl b/fflas-ffpack/fflas/fflas_helpers.inl index d9b8a7b05..f48aa7f1d 100644 --- a/fflas-ffpack/fflas/fflas_helpers.inl +++ b/fflas-ffpack/fflas/fflas_helpers.inl @@ -364,35 +364,32 @@ namespace FFLAS { // to be used in the future, when Winograd's algorithm will be made generic wrt the ModeTrait template - typename std::enable_if::value && - FFLAS::hasBounds::value, void>::type + 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 - typename std::enable_if::value || - !FFLAS::hasBounds::value, void>::type + inline typename std::enable_if::value || + !FFLAS::hasBounds::value, void>::type copyOutBounds(const MMH1& Source, MMH2& Dest){} - template - void mergeOutBounds (const MMHelper &H1, - const MMHelper &H2, - const MMHelper &H3, - MMHelper & Dest){ - Dest.Outmax = max4 (H1.Outmax, H2.Outmax, H3.Outmax, Dest.Outmax); - Dest.Outmin = min4 (H1.Outmin, H2.Outmin, H3.Outmin, Dest.Outmin); - } - template - void mergeOutBounds (const MMHelper &H1, - const MMHelper &H2, - const MMHelper &H3, - MMHelper & 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 - void mergeOutBounds (const MMH1& H1, const MMH2& H2, const MMH3& H3, MMH4& Dest){} + 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 && From aa894698592bb54742c3e05bebfb0e28e8e79ee1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Cl=C3=A9ment=20Pernet?= Date: Mon, 3 Jul 2023 23:32:18 +0200 Subject: [PATCH 9/9] fixed --- benchmarks/benchmark-fgemm.C | 185 +++++++++--------- fflas-ffpack/fflas/fflas_fadd.h | 37 ++-- fflas-ffpack/fflas/fflas_fgemm.inl | 32 ++- .../fflas/fflas_fgemm/fgemm_classical.inl | 8 +- .../fflas_fgemm/schedule_winograd_acc.inl | 40 ++-- fflas-ffpack/fflas/fflas_level2.inl | 6 +- fflas-ffpack/field/field-traits.h | 7 +- tests/test-fgemm.C | 28 +-- 8 files changed, 191 insertions(+), 152 deletions(-) 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/schedule_winograd_acc.inl b/fflas-ffpack/fflas/fflas_fgemm/schedule_winograd_acc.inl index 2aada96b4..253f1e151 100644 --- a/fflas-ffpack/fflas/fflas_fgemm/schedule_winograd_acc.inl +++ b/fflas-ffpack/fflas/fflas_fgemm/schedule_winograd_acc.inl @@ -206,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; @@ -222,7 +221,6 @@ namespace FFLAS { namespace BLAS3 { } else { DF.assign(betadf, beta); } - size_t ldX3; if (ta == FflasTrans) { @@ -272,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, 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 611428b4a..7977ce4ef 100644 --- a/fflas-ffpack/field/field-traits.h +++ b/fflas-ffpack/field/field-traits.h @@ -181,7 +181,7 @@ namespace FFLAS { /* Traits */ template struct ModeTraits >{typedef typename ModeCategories::DelayedTag value;}; template<> struct ModeTraits > {typedef typename ModeCategories::DefaultTag 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;}; @@ -339,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 6c5377a93..67a9a5d7f 100644 --- a/tests/test-fgemm.C +++ b/tests/test-fgemm.C @@ -86,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); @@ -244,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) @@ -314,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); @@ -338,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; @@ -405,11 +404,12 @@ 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);