From 50b37e2f2609887e4d8217490d1ec3904ea852d9 Mon Sep 17 00:00:00 2001 From: Pascal Giorgi Date: Sun, 25 Jun 2023 22:59:08 +0200 Subject: [PATCH 1/8] add arm64 support for simd128 -> TODO must ne changer to use a better coverage of 64-bits arch --- autotune/winograd.C | 2 +- benchmarks/benchmark-dtrsm.C | 2 +- benchmarks/benchmark-fgemm-mp.C | 14 +++---- benchmarks/benchmark-fgemm.C | 16 +++++--- benchmarks/benchmark-ftrsm.C | 10 ++--- benchmarks/benchmark-pluq.C | 2 +- examples/matmul.C | 2 +- fflas-ffpack/config-blas.h | 32 ++++++++++----- .../fflas/fflas_fgemm/fgemm_classical_mp.inl | 7 ++-- fflas-ffpack/fflas/fflas_helpers.inl | 2 +- fflas-ffpack/fflas/fflas_simd/simd128.inl | 2 +- fflas-ffpack/field/rns-integer.h | 41 ++++++++++++------- tests/test-fgemm.C | 2 +- 13 files changed, 82 insertions(+), 52 deletions(-) diff --git a/autotune/winograd.C b/autotune/winograd.C index 303f5d7e5..2ab303ebb 100644 --- a/autotune/winograd.C +++ b/autotune/winograd.C @@ -65,7 +65,7 @@ int main () { typedef FIELD Field; Field F(17); typedef Field::Element Element ; - size_t n=512, nmax=4000, prec=512, nbest=0, count=0; + size_t n=512, nmax=10000, prec=512, nbest=0, count=0; TTimer chrono; bool bound=false; diff --git a/benchmarks/benchmark-dtrsm.C b/benchmarks/benchmark-dtrsm.C index 576835957..5ffd80ce6 100644 --- a/benchmarks/benchmark-dtrsm.C +++ b/benchmarks/benchmark-dtrsm.C @@ -104,7 +104,7 @@ int main(int argc, char** argv) { // ----------- // Standard output for benchmark - Alexis Breust 2014/11/14 std::cout << "Time: " << time / double(iter) - << " Gfops: " << (2.*double(n)/1000.*double(n)/1000.*double(n)/1000.0) / time * double(iter) / 3.; + << " Gfops: " << (1.*double(n)/1000.*double(n)/1000.*double(n)/1000.0) / time * double(iter) ; FFLAS::writeCommandString(std::cout, as) << std::endl; return 0; diff --git a/benchmarks/benchmark-fgemm-mp.C b/benchmarks/benchmark-fgemm-mp.C index 4d349cc06..90058dcab 100644 --- a/benchmarks/benchmark-fgemm-mp.C +++ b/benchmarks/benchmark-fgemm-mp.C @@ -28,7 +28,8 @@ // everywhere in the call stack #define __FFLASFFPACK_OPENBLAS_NT_ALREADY_SET 1 -#define INTEGER_NO_RNS 1 +//#define INTEGER_NO_RNS 1 +#define PROFILE_FGEMM_MP #if not defined(MG_DEFAULT) #define MG_DEFAULT MG_ACTIVE @@ -88,7 +89,7 @@ static Argument as[] = { }; template -int tmain(){ +int tmain(){ srand( (int)seed); srand48(seed); Givaro::Integer::seeding(seed); @@ -183,17 +184,16 @@ int tmain(){ fmpz_mat_clear(BB); #endif //END FLINT CODE // - using FFLAS::CuttingStrategy::Recursive; - using FFLAS::StrategyParameter::TwoDAdaptive; + //using FFLAS::CuttingStrategy::Recursive; + //using FFLAS::StrategyParameter::TwoDAdaptive; // RNS MUL_LA chrono.clear();chrono.start(); // PAR_BLOCK{ // FFLAS::fgemm(F,FFLAS::FflasNoTrans,FFLAS::FflasNoTrans,m,n,k,alpha,A,lda,B,ldb,beta,C,ldc, SPLITTER(NUM_THREADS,Recursive,TwoDAdaptive) ); // } { - FFLAS::fgemm(F,FFLAS::FflasNoTrans,FFLAS::FflasNoTrans,m,n,k,alpha,A,lda,B,ldb,beta,C,ldc,FFLAS::ParSeqHelper::Sequential()); + FFLAS::fgemm(F,FFLAS::FflasNoTrans,FFLAS::FflasNoTrans,m,n,k,alpha,A,lda,B,ldb,beta,C,ldc,FFLAS::ParSeqHelper::Sequential()); } - chrono.stop(); time+=chrono.realtime(); @@ -234,7 +234,7 @@ int tmain(){ int main(int argc, char** argv){ #ifdef __FFLASFFPACK_OPENBLAS_NUM_THREADS - openblas_set_num_threads(__FFLASFFPACK_OPENBLAS_NUM_THREADS); + openblas_set_num_threads(__FFLASFFPACK_OPENBLAS_NUM_THREADS); #endif FFLAS::parseArguments(argc,argv,as); diff --git a/benchmarks/benchmark-fgemm.C b/benchmarks/benchmark-fgemm.C index 91ec24f94..063c831ad 100644 --- a/benchmarks/benchmark-fgemm.C +++ b/benchmarks/benchmark-fgemm.C @@ -28,8 +28,8 @@ //#define WINO_PARALLEL_TMPS //#define __FFLASFFPACK_FORCE_SEQ //#define PFGEMM_WINO_SEQ 32 -//#define CLASSIC_SEQ -#define CLASSIC_HYBRID +#define CLASSIC_SEQ +//#define CLASSIC_HYBRID //#define WINO_SEQ //#define FFT_PROFILER //#define PROFILE_FGEMM_MP @@ -82,15 +82,18 @@ int main(int argc, char** argv) { // typedef Givaro::Modular Field; // typedef Givaro::Modular Field; typedef Givaro::Modular Field; - // typedef Givaro::Modular Field; + //typedef Givaro::Modular Field; // typedef Givaro::ModularBalanced Field; - // typedef Givaro::ModularBalanced Field; + //typedef Givaro::ModularBalanced Field; // typedef Givaro::ModularBalanced Field; // typedef Givaro::Modular Field; + //typedef Givaro::DoubleDomain Field; + + typedef Field::Element Element; Field F(q); - if (q > F.maxCardinality()) return 1; + //if (q > F.maxCardinality()) return 1; Timer chrono, TimFreivalds; double timev=0.0; @@ -99,6 +102,8 @@ int main(int argc, char** argv) { Element * A, * B, * C; Field::RandIter G(F); + //G.setBitsize(20); + A = fflas_new(F,m,k,Alignment::CACHE_PAGESIZE); //#pragma omp parallel for collapse(2) schedule(runtime) PAR_BLOCK { pfrand(F,G, m,k,A,m/size_t(NBK)); } @@ -112,6 +117,7 @@ int main(int argc, char** argv) { PAR_BLOCK { pfzero(F, m,n,C,m/NBK); } + for (size_t i=0;i<=iter;++i){ chrono.clear(); diff --git a/benchmarks/benchmark-ftrsm.C b/benchmarks/benchmark-ftrsm.C index 44d85255e..e239b955e 100644 --- a/benchmarks/benchmark-ftrsm.C +++ b/benchmarks/benchmark-ftrsm.C @@ -81,7 +81,7 @@ int main(int argc, char** argv) { FFLAS::ReadMatrix (file1.c_str(),F,m,m,A); } else{ - A = FFLAS::fflas_new (F,m,m,Alignment::CACHE_PAGESIZE); + A = FFLAS::fflas_new (F,m,m,Alignment::CACHE_PAGESIZE); PAR_BLOCK{ FFLAS::pfrand(F,G,m,m,A,m/NBK); } for (size_t k=0;k<(size_t)m;++k) @@ -102,10 +102,10 @@ for (size_t i=0;i<=iter;++i){ if (i) chrono.start(); if (!p){ - FFLAS::ParSeqHelper::Sequential H; - FFLAS::ftrsm (F, FFLAS::FflasLeft, FFLAS::FflasLower, - FFLAS::FflasNoTrans, FFLAS::FflasNonUnit, - m,n, F.one, A, m, B, n, H); + FFLAS::ParSeqHelper::Sequential H; + FFLAS::ftrsm (F, FFLAS::FflasLeft, FFLAS::FflasLower, + FFLAS::FflasNoTrans, FFLAS::FflasNonUnit, + m,n, F.one, A, m, B, n, H); } else{ FFLAS::ParSeqHelper::Parallel PSH(t); diff --git a/benchmarks/benchmark-pluq.C b/benchmarks/benchmark-pluq.C index 07f19979a..19516a4b1 100644 --- a/benchmarks/benchmark-pluq.C +++ b/benchmarks/benchmark-pluq.C @@ -269,7 +269,7 @@ int main(int argc, char** argv) { if (slab) R = FFPACK::LUdivine (F, diag, FFLAS::FflasNoTrans, m, n, A, n, P, Q); else - R = FFPACK::PLUQ(F, diag, m, n, A, n, P, Q); + R = FFPACK::PLUQ(F, diag, m, n, A, n, P, Q); } if (i) {chrono.stop(); time[i-1]=chrono.realtime();} diff --git a/examples/matmul.C b/examples/matmul.C index 598cb957b..6dc076ef6 100644 --- a/examples/matmul.C +++ b/examples/matmul.C @@ -31,7 +31,7 @@ int main(int argc, char** argv) { if (argc != 4) { - std::cerr << "Usage: matmul

" << std::endl; + std::cerr << "Usage: matmul

" << std::endl; return -1; } diff --git a/fflas-ffpack/config-blas.h b/fflas-ffpack/config-blas.h index e503ab382..ce4932b19 100644 --- a/fflas-ffpack/config-blas.h +++ b/fflas-ffpack/config-blas.h @@ -425,7 +425,8 @@ static const char* EXT_BLAS_DIAG (CBLAS_DIAG t) { if (t == CblasUni // define external link to LAPACK routines extern "C" { //!@bug we should also allow lapacke from MLK - void dgetrf_ (const CBLAS_INT *, const CBLAS_INT *, double *, const CBLAS_INT *, CBLAS_INT *, CBLAS_INT *); + void dgetrf_ (const CBLAS_INT *, const CBLAS_INT *, double *, const CBLAS_INT *, CBLAS_INT *, CBLAS_INT *); + void sgetrf_ (const CBLAS_INT *, const CBLAS_INT *, float *, const CBLAS_INT *, CBLAS_INT *, CBLAS_INT *); void dgetri_ (const CBLAS_INT *, double *, const CBLAS_INT *, const CBLAS_INT *, double *, const CBLAS_INT *, CBLAS_INT *); void dtrtri_ (const char *, const char *, const CBLAS_INT *, double *, const CBLAS_INT *, CBLAS_INT *); void dswap_ (const CBLAS_INT *, double *, const CBLAS_INT *, double *, const CBLAS_INT *); @@ -447,13 +448,23 @@ extern "C" { // return A=P.L.U (L unitary) with ColMajor // return A=L.U.P (U unitary) with RowMajor //! @bug Order is not used. we should use ATLAS/interfaces/lapack/C/src/clapack_dgetrf.c or similar - inline CBLAS_INT clapack_dgetrf(const blas_enum CBLAS_ORDER, const CBLAS_INT M, const CBLAS_INT N, - double *A, const CBLAS_INT lda, CBLAS_INT *ipiv) - { - CBLAS_INT info; - dgetrf_ ( &M, &N, A, &lda, ipiv, &info); - return info; - } + inline CBLAS_INT clapack_dgetrf(const blas_enum CBLAS_ORDER, const CBLAS_INT M, const CBLAS_INT N, + double *A, const CBLAS_INT lda, CBLAS_INT *ipiv) + { + CBLAS_INT info; + dgetrf_ ( &M, &N, A, &lda, ipiv, &info); + return info; + } + + inline CBLAS_INT clapack_sgetrf(const blas_enum CBLAS_ORDER, const CBLAS_INT M, const CBLAS_INT N, + float *A, const CBLAS_INT lda, CBLAS_INT *ipiv) + { + CBLAS_INT info; + sgetrf_ ( &M, &N, A, &lda, ipiv, &info); + return info; + } + + inline CBLAS_INT clapack_dgetri(const blas_enum CBLAS_ORDER, const CBLAS_INT N, double *A, const CBLAS_INT lda, const CBLAS_INT *ipiv) @@ -517,13 +528,12 @@ extern "C" { extern "C" { // LAPACK routines - CBLAS_INT clapack_dgetrf(const blas_enum CBLAS_ORDER Order, const CBLAS_INT M, const CBLAS_INT N, - double *A, const CBLAS_INT lda, CBLAS_INT *ipiv); + CBLAS_INT clapack_dgetrf(const blas_enum CBLAS_ORDER Order, const CBLAS_INT M, const CBLAS_INT N, + double *A, const CBLAS_INT lda, CBLAS_INT *ipiv); CBLAS_INT clapack_dgetri(const blas_enum CBLAS_ORDER Order, const CBLAS_INT N, double *A, const CBLAS_INT lda, const CBLAS_INT *ipiv); CBLAS_INT clapack_dtrtri(const blas_enum CBLAS_ORDER Order,const blas_enum CBLAS_UPLO Uplo, const blas_enum CBLAS_DIAG Diag,const CBLAS_INT N, double *A, const CBLAS_INT lda); - } #endif // CLAPACK ? diff --git a/fflas-ffpack/fflas/fflas_fgemm/fgemm_classical_mp.inl b/fflas-ffpack/fflas/fflas_fgemm/fgemm_classical_mp.inl index 14a2ac3a1..2cd44decb 100644 --- a/fflas-ffpack/fflas/fflas_fgemm/fgemm_classical_mp.inl +++ b/fflas-ffpack/fflas/fflas_fgemm/fgemm_classical_mp.inl @@ -322,7 +322,8 @@ namespace FFLAS { Givaro::Integer* C, const size_t ldc, MMHelper, MMHelperAlgo::Classic, ModeCategories::ConvertTo, ParSeq > & H) { - //std::cerr<<"Entering fgemm> ParSeq"<> ParSeq"<, MMHelperAlgo::Classic, ModeCategories::ConvertTo > & H) { // compute the product over Z - //std::cerr<<"Entering fgemm>"<>"< IntegerDomain; Givaro::Integer p; F.cardinality(p); diff --git a/fflas-ffpack/fflas/fflas_helpers.inl b/fflas-ffpack/fflas/fflas_helpers.inl index 85478d4d2..15c56a31d 100644 --- a/fflas-ffpack/fflas/fflas_helpers.inl +++ b/fflas-ffpack/fflas/fflas_helpers.inl @@ -195,7 +195,7 @@ namespace FFLAS { } bool Aunfit(){ return Protected::unfit(std::max(static_cast(-Amin),Amax));} bool Bunfit(){ return Protected::unfit(std::max(static_cast(-Bmin),Bmax));} - void setOutBounds(const size_t k, const DFElt alpha, const DFElt beta) + void setOutBounds(const uint64_t k, const DFElt alpha, const DFElt beta) { if (beta<0){ Outmin = beta*Cmax; diff --git a/fflas-ffpack/fflas/fflas_simd/simd128.inl b/fflas-ffpack/fflas/fflas_simd/simd128.inl index df87a9ae3..09c14904a 100644 --- a/fflas-ffpack/fflas/fflas_simd/simd128.inl +++ b/fflas-ffpack/fflas/fflas_simd/simd128.inl @@ -105,7 +105,7 @@ Simd128_impl::value, std::is_integral::value, std::is_s #include "simd128_int16.inl" #include "simd128_int32.inl" -#ifdef __x86_64__ +#if defined( __x86_64__) or defined(__aarch64__) #include "simd128_int64.inl" #endif #endif //#ifdef SIMD_INT diff --git a/fflas-ffpack/field/rns-integer.h b/fflas-ffpack/field/rns-integer.h index 27782c83f..95ada64d8 100644 --- a/fflas-ffpack/field/rns-integer.h +++ b/fflas-ffpack/field/rns-integer.h @@ -111,7 +111,7 @@ namespace FFPACK { } Element& init(Element& x, const Givaro::Integer& y) const{ init(x); - size_t k =(y.bitsize())/16+((y.bitsize())%16?1:0); + uint64_t k =(y.bitsize())/16+((y.bitsize())%16?1:0); _rns->init(1,1,x._ptr,x._stride, &y,1,k); return x; } @@ -153,19 +153,32 @@ namespace FFPACK { namespace FFLAS { - // specialization for the fflas alloc function - template<> - inline FFPACK::rns_double_elt_ptr - fflas_new(const FFPACK::RNSInteger &F, const size_t m, const Alignment align){ - double *ptr=FFLAS::fflas_new(m*F.size(), align); - return FFPACK::rns_double_elt_ptr(ptr,m); - } - - template<> - inline FFPACK::rns_double_elt_ptr - fflas_new(const FFPACK::RNSInteger &F, const size_t m, const size_t n, const Alignment align){ - return fflas_new(F, m*n, align); - } + // specialization for the fflas alloc function + template<> + inline FFPACK::rns_double_elt_ptr + fflas_new(const FFPACK::RNSInteger &F, const size_t m, const Alignment align){ + double *ptr=FFLAS::fflas_new(m*F.size(), align); + return FFPACK::rns_double_elt_ptr(ptr,m); + } + + template<> + inline FFPACK::rns_double_elt_ptr + fflas_new(const FFPACK::RNSInteger &F, const size_t m, const size_t n, const Alignment align){ + return fflas_new(F, m*n, align); + } + // specialization for the fflas alloc function + template<> + inline FFPACK::rns_double_elt_ptr + fflas_new(const FFPACK::RNSInteger &F, const size_t m, const Alignment align){ + double *ptr=FFLAS::fflas_new(m*F.size(), align); + return FFPACK::rns_double_elt_ptr(ptr,m); + } + + template<> + inline FFPACK::rns_double_elt_ptr + fflas_new(const FFPACK::RNSInteger &F, const size_t m, const size_t n, const Alignment align){ + return fflas_new(F, m*n, align); + } // function to convert from integer to RNS (note: this is not the finit function from FFLAS, extra k) template diff --git a/tests/test-fgemm.C b/tests/test-fgemm.C index 443e51d5f..9131be3da 100644 --- a/tests/test-fgemm.C +++ b/tests/test-fgemm.C @@ -206,7 +206,7 @@ bool launch_MM(const Field & F, fgemm (F, ta, tb,m,n,k,alpha, A,lda, B,ldb, beta,C,ldc,WH); } ok = ok && check_MM(F, D, ta, tb,m,n,k,alpha, A,lda, B,ldb, beta,C,ldc); - +prln fflas_delete(A); fflas_delete(B); From 5416218560a125fec18606dc300d048883780a06 Mon Sep 17 00:00:00 2001 From: Pascal Giorgi Date: Wed, 13 Sep 2023 17:51:59 +0200 Subject: [PATCH 2/8] update rns-double with some comments and retrieve rns-double-extended code --- benchmarks/benchmark-fgemm.C | 2 +- fflas-ffpack/field/Makefile.am | 2 + fflas-ffpack/field/rns-double-extended.h | 402 ++++++++ fflas-ffpack/field/rns-double-extended.inl | 1015 ++++++++++++++++++++ fflas-ffpack/field/rns-double.h | 228 +---- fflas-ffpack/field/rns-double.inl | 93 +- fflas-ffpack/field/rns-integer.h | 17 +- 7 files changed, 1447 insertions(+), 312 deletions(-) create mode 100644 fflas-ffpack/field/rns-double-extended.h create mode 100644 fflas-ffpack/field/rns-double-extended.inl diff --git a/benchmarks/benchmark-fgemm.C b/benchmarks/benchmark-fgemm.C index 063c831ad..3378b1537 100644 --- a/benchmarks/benchmark-fgemm.C +++ b/benchmarks/benchmark-fgemm.C @@ -211,7 +211,7 @@ int main(int argc, char** argv) { } else{ - MMHelper//, + MMHelper//, //typename FieldTraits::value, //ParSeqHelper::Sequential> WH (F, nbw, ParSeqHelper::Sequential()); diff --git a/fflas-ffpack/field/Makefile.am b/fflas-ffpack/field/Makefile.am index a1ac51f67..8f115d761 100644 --- a/fflas-ffpack/field/Makefile.am +++ b/fflas-ffpack/field/Makefile.am @@ -27,9 +27,11 @@ pkgincludesubdir=$(pkgincludedir)/field RNS=rns.h \ rns.inl \ rns-double.h \ +rns-double-extended.h \ rns-double-elt.h \ rns-double.inl \ rns-double-recint.inl \ + rns-double-extended.inl \ rns-integer.h \ rns-integer-mod.h diff --git a/fflas-ffpack/field/rns-double-extended.h b/fflas-ffpack/field/rns-double-extended.h new file mode 100644 index 000000000..e1b581c25 --- /dev/null +++ b/fflas-ffpack/field/rns-double-extended.h @@ -0,0 +1,402 @@ +/* + * Copyright (C) 2014 the FFLAS-FFPACK group + * + * Written by Pascal Giorgi + * + * + * ========LICENCE======== + * This file is part of the library FFLAS-FFPACK. + * + * FFLAS-FFPACK is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with this library; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + * ========LICENCE======== + *. + */ + +/*! @file field/rns-double-extended.h + * @ingroup field + * @brief rns structure with double support (full precision) + */ + +#ifndef __FFPACK_rns_double_extended_H +#define __FFPACK_rns_double_extended_H + +// Bigger multiple of s lesser or equal than x, s must be a power of two +#ifndef ROUND_DOWN +#define ROUND_DOWN(x, s) ((x) & ~((s)-1)) +#endif + +#include +#include +#include +#include +#include "givaro/modular-extended.h" +#include +#include "fflas-ffpack/config-blas.h" +#include "fflas-ffpack/utils/fflas_memory.h" +#include "fflas-ffpack/utils/align-allocator.h" +#include "fflas-ffpack/field/rns-double-elt.h" + +namespace FFPACK { + + // forward declaration of different variants + struct rns_double_extended_V1 ; + struct rns_double_extended_V2 ; + + // set default value to one variant + using rns_double_extended = rns_double_extended_V2; + + +#define RALIGN setw(40)< ModField; + + std::vector> _basis; // the rns moduli (mi) + std::vector> _basisMax; // (mi-1) + std::vector> _negbasis; // (-mi) + std::vector> _invbasis; // the inverse of rns moduli (1/mi) + std::vector _field_rns; // the associated prime field for each mi + integer _M; // the product of the mi's + std::vector _Mi; // _M/mi + std::vector _MMi; // (_Mi)^(-1) mod mi + std::vector _crt_in[6]; // 2^(16*j) mod mi + std::vector _crt_out[6]; // (_Mi._MMi) written in base 2^16 + size_t _size; // the size of the rns basis (number of mi's) + size_t _pbits; // the size in bit of the mi's + size_t _ldm; // log[2^16](_M) + + typedef double BasisElement; + typedef rns_double_elt Element; + typedef rns_double_elt_ptr Element_ptr; + typedef rns_double_elt_cstptr ConstElement_ptr; + + rns_double_extended_V1(const integer& bound, size_t pbits, bool rnsmod=false, long seed=time(NULL)) + : _M(1), _size(0), _pbits(pbits) + { + integer::seeding(seed); + integer prime; Givaro::IntPrimeDom IPD; + integer sum=1; + while (_M < bound*sum) { + _basis.resize(_size+1); + do { + integer::random_exact_2exp(prime, _pbits-1); + IPD.nextprimein(prime); + } while (_M%prime == 0); + _basis[_size]=prime; + _size++; + _M*=prime; + if (rnsmod) sum+=prime; + } + precompute_cst(); + } + + rns_double_extended_V1(size_t pbits, size_t size, long seed=time(NULL)) + : _M(1), _size(size), _pbits(pbits) + { + integer::seeding(seed); + integer prime; Givaro::IntPrimeDom IPD; + integer sum=1; + _basis.resize(size); + _negbasis.resize(size); + _basisMax.resize(size); + for(size_t i = 0 ; i < _size ; ++i){ + integer::random_exact_2exp(prime, _pbits-1); + IPD.nextprimein(prime); + _basis[i]=prime; + _basisMax[i] = prime-1; + _negbasis[i] = 0-prime; + _M*=prime; + } + precompute_cst(); + } + + template + rns_double_extended_V1(const Vect& basis, bool rnsmod=false, long seed=time(NULL)) + : _basis(basis.begin(),basis.end()), _basisMax(basis.size()), _negbasis(basis.size()), _M(1), _size(basis.size()), _pbits(0) + { + for(size_t i=0;i<_size;i++){ + _M*=_basis[i]; + _pbits=std::max(_pbits, integer(_basis[i]).bitsize()); + } + precompute_cst(); + } + + + void precompute_cst(){ + _ldm = (_M.bitsize()/48) + ((_M.bitsize()%48)?1:0) ; + _invbasis.resize(_size); + _basisMax.resize(_size); + _negbasis.resize(_size); + _field_rns.resize(_size); + _Mi.resize(_size); + _MMi.resize(_size); + for (size_t i=0;i<6;i++){ + _crt_in[i] .resize(_size*_ldm); + _crt_out[i].resize(_size*_ldm); + } + const unsigned int MASK=0xFFFF; + for (size_t i=0;i<_size;i++){ + _invbasis[i] = 1./_basis[i]; + _basisMax[i] = _basis[i]-1; + _negbasis[i] = 0-_basis[i]; + _field_rns[i] = ModField(_basis[i]); + _Mi[i] = _M/(uint64_t)_basis[i]; + _field_rns[i].init(_MMi[i], _Mi[i] % (double)_basis[i]); + _field_rns[i].invin(_MMi[i]); + integer tmp= _Mi[i]*(uint64_t)_MMi[i]; + double a0,a1,a2; + for(size_t j=0;j<_ldm;j++){ + uint64_t limb= tmp[0]; + a0= double(limb&MASK); + a1= double((limb>>16)&MASK); + a2= double((limb>>32)&MASK); + _crt_out[0][j+i*_ldm]=a0; + _crt_out[1][j+i*_ldm]=a1; + _crt_out[2][j+i*_ldm]=a2; + _crt_out[3][j+i*_ldm]=a0+a1; + _crt_out[4][j+i*_ldm]=a1+a2; + _crt_out[5][j+i*_ldm]=a0+a1+a2; + tmp>>=48; + } + double beta=double(1UL<<48); + double acc=1; + for(size_t j=0;j<_ldm;j++){ + uint64_t limb= acc; + a0= double(limb&MASK); + a1= double((limb>>16)&MASK); + a2= double((limb>>32)&MASK); + _crt_in[0][j+i*_ldm]=a0; + _crt_in[1][j+i*_ldm]=a1; + _crt_in[2][j+i*_ldm]=a2; + _crt_in[3][j+i*_ldm]=a0+a1; + _crt_in[4][j+i*_ldm]=a1+a2; + _crt_in[5][j+i*_ldm]=a0+a1+a2; + _field_rns[i].mulin(acc,beta); + } + } +#ifdef RNS_DEBUG + std::cout<<"basis= ["; + for(size_t i=0;i<_size;i++) + std::cout<<(int64_t)_basis[i]<<(i!=_size-1?",":"]\n"); + Givaro::ModularExtended ZZ(2UL<<48);; + for(size_t l=0;l<6;l++){ + std::cout<<"crtin"< ModField; + + std::vector> _basis; // the rns moduli (mi) + std::vector> _basisMax; // (mi-1) + std::vector> _negbasis; // (-mi) + std::vector> _invbasis; // the inverse of rns moduli (1/mi) + std::vector _field_rns; // the associated prime field for each mi + integer _M; // the product of the mi's + std::vector _Mi; // _M/mi + std::vector _MMi; // (_Mi)^(-1) mod mi + std::vector _crt_in; // 2^(16*j) mod mi -> each entry are splitted into two line (a0 + 2^27 a1) + std::vector _crt_out; // (_Mi._MMi) written in base 2^16 + size_t _size; // the size of the rns basis (number of mi's) + size_t _pbits; // the size in bit of the mi's + size_t _ldm; // log[2^16](_M) + uint64_t _shift; + + typedef double BasisElement; + typedef rns_double_elt Element; + typedef rns_double_elt_ptr Element_ptr; + typedef rns_double_elt_cstptr ConstElement_ptr; + + rns_double_extended_V2 (const integer& bound, size_t pbits, bool rnsmod=false, long seed=time(NULL)) + : _M(1), _size(0), _pbits(pbits) + { + if (pbits > 52){ + std::cerr<<"FFLAS- RNS EXTENDED ERROR: prime bitsize above 52 bits ... log(p)="<= (integer(1)<<37) ){ + std::cerr<<"FFLAS- RNS EXTENDED ERROR: prime bitsize ("<(&tmp); + const uint16_t* m0_ptr = reinterpret_cast(m0[0]->_mp_d); + size_t maxs=std::min(_ldm,(tmp.size())*sizeof(mp_limb_t)/2);// to ensure 32 bits portability +#ifdef BENCH_RNS_PRECOMP + chrono.stop(); + t1+=chrono.usertime(); + chrono.start(); +#endif + size_t l=0; + for(;l> _shift); + + _field_rns[i].mulin(acc,beta); + //std::cout<<"RNS precomp ("< "<< (int64_t)_crt_in[j+i*_ldm]<<" "<< (int64_t)_crt_in[j+(i+_size)*_ldm]< ZZ(2UL<<48);; + std::cout<<"CRTmat:="; + //write_field(ZZ,std::cout,_crt_out.data(), _size, _ldm,_ldm,true); + FFLAS::WriteMatrix(std::cout, ZZ, _size, _ldm, _crt_out.data(), _ldm); +#endif + } + + // Arns must be an array of m*n*_size + // abs(||A||) <= maxA + template + void init(size_t m, size_t n, double* Arns, size_t rda, const T* A, size_t lda, const integer& maxA, bool RNS_MAJOR=false) const + { + init(m,n,Arns,rda,A,lda, maxA.bitsize()/16 + (maxA.bitsize()%16?1:0),RNS_MAJOR); + } + + void init(size_t m, size_t n, double* Arns, size_t rda, const integer* A, size_t lda, size_t k, bool RNS_MAJOR=false) const; + // void init_transpose(size_t m, size_t n, double* Arns, size_t rda, const integer* A, size_t lda, size_t k, bool RNS_MAJOR=false) const; + void convert(size_t m, size_t n, integer gamma, integer* A, size_t lda, const double* Arns, size_t rda, bool RNS_MAJOR=false) const; + // void convert_transpose(size_t m, size_t n, integer gamma, integer* A, size_t lda, const double* Arns, size_t rda, bool RNS_MAJOR=false) const; + + // reduce entries of Arns to be less than the rns basis elements + void reduce(size_t n, double* Arns, size_t rda, bool RNS_MAJOR=false) const; + + }; + + + + +} // end of namespace FFPACK + +#include "rns-double-extended.inl" + +#endif //__FFPACK_rns_double_extended_H +/* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ +// vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s diff --git a/fflas-ffpack/field/rns-double-extended.inl b/fflas-ffpack/field/rns-double-extended.inl new file mode 100644 index 000000000..9a63436c4 --- /dev/null +++ b/fflas-ffpack/field/rns-double-extended.inl @@ -0,0 +1,1015 @@ +/* + * Copyright (C) 2016 the FFLAS-FFPACK group + * + * Written by Pascal Giorgi + * + * + * ========LICENCE======== + * This file is part of the library FFLAS-FFPACK. + * + * FFLAS-FFPACK is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with this library; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + * ========LICENCE======== + *. + */ + + +#ifndef __FFLASFFPACK_field_rns_double_extended_INL +#define __FFLASFFPACK_field_rns_double_extended_INL +#include "fflas-ffpack/utils/Matio.h" +#include "fflas-ffpack/fflas/fflas_freduce.h" + +namespace FFLAS { + template<> + void fscalin(const Givaro::ModularExtended & F, const size_t n, double a, + double* X, const size_t incX) + { + for(size_t i=0;i; + using vect_t = typename simd::vect_t; + + if(_size % simd::vect_size == 0){ + //#pragma omp parallel for schedule(static, 256) + for(size_t i = 0 ; i < n ; i++){ + vect_t tmp1, tmp2, v, min, max, basis, inv, neg; + min = simd::set1(0.); + for(size_t j = 0 ; j < _size ; j+=simd::vect_size){ + basis = simd::load(_basis.data()+j); + inv = simd::load(_invbasis.data()+j); + max = simd::load(_basisMax.data()+j); + neg = simd::load(_negbasis.data()+j); + v = simd::load(Arns+i*_size+j); + simd::mod(v, basis, inv, neg, min, max, tmp1,tmp2); + simd::store(Arns+i*_size+j, v); + } + } + } else{ + //#pragma omp parallel for schedule(static, 256) + for(size_t i = 0 ; i < n ; i++){ + vect_t tmp1, tmp2, tmp3, v, min, max, basis, inv, neg; + size_t j = 0; + for( ; j < ROUND_DOWN(_size, simd::vect_size) ; j+=simd::vect_size){ + basis = simd::load(_basis.data()+j); + inv = simd::load(_invbasis.data()+j); + max = simd::load(_basisMax.data()+j); + neg = simd::load(_negbasis.data()+j); + v = simd::load(Arns+i*_size+j); + simd::mod(v, basis, inv, neg, min, max, tmp1,tmp2); + simd::store(Arns+i*_size+j, v); + } + for( ; j < _size ; ++j){ + _field_rns[j].reduce(Arns[i*_size+j]); + } + } + } +#else + for(size_t i=0;i<_size;i++) + FFLAS::freduce (_field_rns[i],n,Arns+i*rda,1); + +#endif + + } + + + + void rns_double_extended_V1::init(size_t m, size_t n, double* Arns, size_t rda, const integer* A, size_t lda, size_t k, bool RNS_MAJOR) const + { + if (_size>= (1<<16)){ + std::cerr<<"RNS EXTENDED DOUBLE: init Error -> the nbr of moduli in RNS basis is > 2^16, not implemented. aborting\n";std::terminate(); + } +#ifdef BENCH_RNS + if (m!=1 && n!=1){ + std::cerr< rns size ("<<_size<<") kronecker size ("< Numbit(M)="<<_M.bitsize()<_ldm){ + FFPACK::failure()(__func__,__FILE__,__LINE__,"rns_struct: init (too large entry)"); + std::cerr<<"k="< "< "<<((tmp[0] & 0xFFFF)) << " != "< "<<((tmp[0] & 0xFFFFFFFF)>>16) << " != "< "<<((tmp[0] & 0xFFFFFFFFFFFF)>>32) << " != "<>=48; + } + } + } +#endif + tkr.stop(); +#ifdef BENCH_RNS + if (m!=1 && n!=1) + std::cerr< ZD; + + // Using Helper for potential parallelism -> need to be activated by hand + FFLAS::MMHelper, FFLAS::MMHelperAlgo::Winograd> MMH (ZD, -1, FFLAS::ParSeqHelper::Sequential()); + + if (RNS_MAJOR==false) { + // A_rns = _crt_in x A_beta^T + Givaro::Timer tfgemm; tfgemm.start(); + FFLAS::fgemm (ZD, FFLAS::FflasNoTrans,FFLAS::FflasTrans,_size,mn,k,1.0,_crt_in[0].data(),_ldm,A_beta0,k,0.,A_rns0, mn, MMH); + FFLAS::fgemm (ZD, FFLAS::FflasNoTrans,FFLAS::FflasTrans,_size,mn,k,1.0,_crt_in[1].data(),_ldm,A_beta1,k,0.,A_rns1, mn, MMH); + FFLAS::fgemm (ZD, FFLAS::FflasNoTrans,FFLAS::FflasTrans,_size,mn,k,1.0,_crt_in[2].data(),_ldm,A_beta2,k,0.,A_rns2, mn, MMH); + FFLAS::fgemm (ZD, FFLAS::FflasNoTrans,FFLAS::FflasTrans,_size,mn,k,1.0,_crt_in[3].data(),_ldm,A_beta3,k,0.,A_rns3, mn, MMH); + FFLAS::fgemm (ZD, FFLAS::FflasNoTrans,FFLAS::FflasTrans,_size,mn,k,1.0,_crt_in[4].data(),_ldm,A_beta4,k,0.,A_rns4, mn, MMH); + FFLAS::fgemm (ZD, FFLAS::FflasNoTrans,FFLAS::FflasTrans,_size,mn,k,1.0,_crt_in[5].data(),_ldm,A_beta5,k,0.,A_rns5, mn, MMH); + + tfgemm.stop(); +#ifdef BENCH_RNS + if (m!=1 && n!=1) + std::cerr<=0?0:_basis[i]); + + q1= std::floor(c1*two_16_invmi);//std::cout<<"Q1:="<<(int64_t)q1<<";\n"; + q2= std::floor(c2*two_32_invmi);//std::cout<<"Q2:="<<(int64_t)q2<<";\n"; + q3= std::floor(c3*two_48_invmi);//std::cout<<"Q3:="<<(int64_t)q3<<";\n"; + q4= std::floor(c4*two_64_invmi);//std::cout<<"Q4:="<<(int64_t)q4<<";\n"; + + c1*=two16; + c2*=two32; + c3*=two48; + c4*=two64; + + c1=fma(q1,_negbasis[i],c1); //std::cout<<"D1:="<<(int64_t)c1<<";\n"; + c2=fma(q2,_negbasis[i],c2); //std::cout<<"D2:="<<(int64_t)c2<<";\n"; + c3=fma(q3,_negbasis[i],c3); //std::cout<<"D3:="<<(int64_t)c3<<";\n"; + c4=fma(q4,_negbasis[i],c4); //std::cout<<"D4:="<<(int64_t)c4<<";\n"; + + c0+=c1; + c2+=c3; + c0+=c4+c2; + while(c0<0.) c0+=_basis[i]; + while(c0>_basis[i]) c0-=_basis[i]; + Arns[j+i*rda]= c0; + + // c1+=c2<<16; + // c3+=c4<<16; + // _field_rns[i].reduce(C1,c1); + // _field_rns[i].reduce(C3,c3); + // _field_rns[i].axpy(C1,C1,two16_mod_mi,c0); + // _field_rns[i].mul(C3,C3,two48_mod_mi); + // _field_rns[i].add(C,C1,C3); + // Arns[j+i*rda]= C+(C>=0?0:_basis[i]); + } + ); + //reduce(mn,Arns,rda,RNS_MAJOR); + + tred.stop(); +#ifdef BENCH_RNS + if (m!=1 && n!=1) + std::cerr< "< "<<((A[i*lda+j] % (int64_t) _basis[k])+(A[i*lda+j]<0?(int64_t)_basis[k]:0)) + // <<" == " + // <<(int64_t) Arns[i*n+j+k*rda] + // <<" --> "<= (1<<16)){ + std::cerr<<"RNS EXTENDED DOUBLE: convert Error -> the nbr of moduli in RNS basis is > 2^16, not implemented. aborting\n";std::terminate(); + } +#ifdef BENCH_RNS + if (m!=1 && n!=1) + std::cerr< rns size ("<<_size<<") kronecker size ("<<_ldm<<") data dim ("<>1; + size_t mn= m*n; + size_t mnldm=mn*_ldm; + double *A_beta= FFLAS::fflas_new(6*mnldm); + double *A_beta0=A_beta; + double *A_beta1=A_beta+1*mnldm; + double *A_beta2=A_beta+2*mnldm; + double *A_beta3=A_beta+3*mnldm; + double *A_beta4=A_beta+4*mnldm; + double *A_beta5=A_beta+5*mnldm; + size_t mnsize=mn*_size; + double *A_rns_tmp = FFLAS::fflas_new(mnsize*6); + double *A_rns0=A_rns_tmp; + double *A_rns1=A_rns_tmp+1*mnsize; + double *A_rns2=A_rns_tmp+2*mnsize; + double *A_rns3=A_rns_tmp+3*mnsize; + double *A_rns4=A_rns_tmp+4*mnsize; + double *A_rns5=A_rns_tmp+5*mnsize; + + Givaro::Timer tsplit; + tsplit.start(); + uint64_t tmp,idx; + double aa0,aa1,aa2; + for (size_t i=0;i<_size;i++) + for(size_t j=0;j>16)&0xFFFF; + aa2=tmp>>32; + A_rns0[idx]=aa0; + A_rns1[idx]=aa1; + A_rns2[idx]=aa2; + A_rns3[idx]=aa0+aa1; + A_rns4[idx]=aa1+aa2; + A_rns5[idx]=aa0+aa1+aa2; + } + tsplit.stop(); +#ifdef BENCH_RNS + if (m!=1 && n!=1) + std::cerr<(),FFLAS::FflasTrans, FFLAS::FflasNoTrans,(int) mn,(int) _ldm,(int) _size, 1.0 , A_rns0,(int) rda, _crt_out[0].data(),(int) _ldm, 0., A_beta0,(int)_ldm, FFLAS::ParSeqHelper::Parallel()); + FFLAS::fgemm(Givaro::ZRing(),FFLAS::FflasTrans, FFLAS::FflasNoTrans,(int) mn,(int) _ldm,(int) _size, 1.0 , A_rns1,(int) rda, _crt_out[1].data(),(int) _ldm, 0., A_beta1,(int)_ldm, FFLAS::ParSeqHelper::Parallel()); + FFLAS::fgemm(Givaro::ZRing(),FFLAS::FflasTrans, FFLAS::FflasNoTrans,(int) mn,(int) _ldm,(int) _size, 1.0 , A_rns2,(int) rda, _crt_out[2].data(),(int) _ldm, 0., A_beta2,(int)_ldm, FFLAS::ParSeqHelper::Parallel()); + FFLAS::fgemm(Givaro::ZRing(),FFLAS::FflasTrans, FFLAS::FflasNoTrans,(int) mn,(int) _ldm,(int) _size, 1.0 , A_rns3,(int) rda, _crt_out[3].data(),(int) _ldm, 0., A_beta3,(int)_ldm, FFLAS::ParSeqHelper::Parallel()); + FFLAS::fgemm(Givaro::ZRing(),FFLAS::FflasTrans, FFLAS::FflasNoTrans,(int) mn,(int) _ldm,(int) _size, 1.0 , A_rns4,(int) rda, _crt_out[4].data(),(int) _ldm, 0., A_beta4,(int)_ldm, FFLAS::ParSeqHelper::Parallel()); + FFLAS::fgemm(Givaro::ZRing(),FFLAS::FflasTrans, FFLAS::FflasNoTrans,(int) mn,(int) _ldm,(int) _size, 1.0 , A_rns5,(int) rda, _crt_out[5].data(),(int) _ldm, 0., A_beta5,(int)_ldm, FFLAS::ParSeqHelper::Parallel()); + } + else { + // compute A_beta = Ap x M_Beta + std::cerr<<"NOT YET IMPLEMENTED .... aborting\n"; std::terminate(); + } + tfgemmc.stop(); +#ifdef BENCH_RNS + if (m!=1 && n!=1) + std::cerr< A0_tmp(k64<<2,0),A1_tmp(k64<<2,0),A2_tmp(k64<<2,0),A3_tmp(k64<<2,0),A4_tmp(k64<<2,0); + uint48_t *A0,*A1,*A2,*A3,*A4; + A0= reinterpret_cast(A0_tmp.data()); + A1= reinterpret_cast(A1_tmp.data()+1); + A2= reinterpret_cast(A2_tmp.data()+2); + A3= reinterpret_cast(A3_tmp.data()+3); + A4= reinterpret_cast(A4_tmp.data()+4); + integer a0,a1,a2,a3,a4,res; + mpz_t *m0,*m1,*m2,*m3,*m4; + m0= reinterpret_cast(&a0); + m1= reinterpret_cast(&a1); + m2= reinterpret_cast(&a2); + m3= reinterpret_cast(&a3); + m4= reinterpret_cast(&a4); + mp_limb_t *m0_d,*m1_d,*m2_d,*m3_d,*m4_d; + m0_d = m0[0]->_mp_d; + m1_d = m1[0]->_mp_d; + m2_d = m2[0]->_mp_d; + m3_d = m3[0]->_mp_d; + m4_d = m4[0]->_mp_d; + m0[0]->_mp_alloc = m1[0]->_mp_alloc = m2[0]->_mp_alloc = m3[0]->_mp_alloc = m4[0]->_mp_alloc = (int) (k64*8/sizeof(mp_limb_t)); // to ensure 32 bits portability + m0[0]->_mp_size = m1[0]->_mp_size = m2[0]->_mp_size = m3[0]->_mp_size = m4[0]->_mp_size = (int) (k64*8/sizeof(mp_limb_t)); // to ensure 32 bits portability + // auto sp=SPLITTER(); + // PARFOR1D(i,m,sp, +#ifdef RNS_DEBUG + std::cout<<"M:="<<_M<_mp_d= reinterpret_cast(A0_tmp.data()); + m1[0]->_mp_d= reinterpret_cast(A1_tmp.data()); + m2[0]->_mp_d= reinterpret_cast(A2_tmp.data()); + m3[0]->_mp_d= reinterpret_cast(A3_tmp.data()); + m4[0]->_mp_d= reinterpret_cast(A4_tmp.data()); + res = a0;res+= a1;res+= a2;res+= a3;res+=a4; + res%=_M; + +#ifdef RNS_DEBUG + std::cout<<"a0:="<hM) + res-=_M; + if (gamma==0) + Aiter[j+i*lda]=res; + else + if (gamma==integer(1)) + Aiter[j+i*lda]+=res; + else + if (gamma==integer(-1)) + Aiter[j+i*lda]=res-Aiter[j+i*lda]; + else{ + Aiter[j+i*lda]*=gamma; + Aiter[j+i*lda]+=res; + } + + } + tkroc.stop(); +#ifdef BENCH_RNS + if (m!=1 && n!=1) + std::cerr<_mp_d = m0_d; + m1[0]->_mp_d = m1_d; + m2[0]->_mp_d = m2_d; + m3[0]->_mp_d = m3_d; + m4[0]->_mp_d = m4_d; + m0[0]->_mp_alloc = m1[0]->_mp_alloc = m2[0]->_mp_alloc= m3[0]->_mp_alloc = m4[0]->_mp_alloc = 1; + m0[0]->_mp_size = m1[0]->_mp_size = m2[0]->_mp_size = m3[0]->_mp_size = m4[0]->_mp_size = 0; + FFLAS::fflas_delete( A_beta); + +#ifdef CHECK_RNS + bool ok=true; + for (size_t i=0;i "<; + using vect_t = typename simd::vect_t; + + if(_size % simd::vect_size == 0){ + //#pragma omp parallel for schedule(static, 256) + for(size_t i = 0 ; i < n ; i++){ + vect_t tmp1, tmp2, v, min, max, basis, inv, neg; + min = simd::set1(0.); + for(size_t j = 0 ; j < _size ; j+=simd::vect_size){ + basis = simd::load(_basis.data()+j); + inv = simd::load(_invbasis.data()+j); + max = simd::load(_basisMax.data()+j); + neg = simd::load(_negbasis.data()+j); + v = simd::load(Arns+i*_size+j); + simd::mod(v, basis, inv, neg, min, max, tmp1,tmp2); + simd::store(Arns+i*_size+j, v); + } + } + } else{ + //#pragma omp parallel for schedule(static, 256) + for(size_t i = 0 ; i < n ; i++){ + vect_t tmp1, tmp2, tmp3, v, min, max, basis, inv, neg; + size_t j = 0; + for( ; j < ROUND_DOWN(_size, simd::vect_size) ; j+=simd::vect_size){ + basis = simd::load(_basis.data()+j); + inv = simd::load(_invbasis.data()+j); + max = simd::load(_basisMax.data()+j); + neg = simd::load(_negbasis.data()+j); + v = simd::load(Arns+i*_size+j); + simd::mod(v, basis, inv, neg, min, max, tmp1,tmp2); + simd::store(Arns+i*_size+j, v); + } + for( ; j < _size ; ++j){ + _field_rns[j].reduce(Arns[i*_size+j]); + } + } + } +#else + for(size_t i=0;i<_size;i++) + FFLAS::freduce (_field_rns[i],n,Arns+i*rda,1); + +#endif + + } + + + inline void rns_double_extended_V2::init(size_t m, size_t n, double* Arns, size_t rda, const integer* A, size_t lda, size_t k, bool RNS_MAJOR) const + { + if (k>_ldm){ + FFPACK::failure()(__func__,__FILE__,__LINE__,"rns_struct: init (too large entry)"); + std::cerr<<"k="< rns size ("<<_size<<") kronecker size ("< Numbit(M)="<<_M.bitsize()<(Aiter+j+i*lda); + const uint16_t* m0_ptr = reinterpret_cast(m0[0]->_mp_d); + size_t l=0; + //size_t maxs=std::min(k,(Aiter[j+i*lda].size())<<2); + size_t maxs=std::min(k,(Aiter[j+i*lda].size())*sizeof(mp_limb_t)/2);// to ensure 32 bits portability + + if (m0[0]->_mp_size >= 0) + for (;l(), FFLAS::FflasNoTrans,FFLAS::FflasTrans,2*_size,mn,k,1.0,_crt_in.data(),_ldm,A_beta,k,0.,Arns_tmp,mn, + FFLAS::ParSeqHelper::Parallel()); + tfgemm.stop(); +#ifdef BENCH_RNS + if(m>1 && n>1) std::cerr<<"RNS EXTENDED double (To) - fgemm : "<1 && n>1) std::cerr<<"RNS EXTENDED double (To) - Reduce : "< "< rns size ("<<_size<<") kronecker size ("<<_ldm<<") data dim ("<(mn*_size*2),accTmp; + double* Arns_tmp2= Arns_tmp+mn; + double mask=double((1<<_shift)+1), scaling=1./double(1<<_shift);; + for (size_t i=0;i<_size;i++) + for (size_t j=0;j> _shift); + } + + +#ifdef RNS_DEBUG + Givaro::ModularExtended ZZ(2UL<<48);; + std::cout<<"Arns:="; + FFLAS::WriteMatrix(std::cout, ZZ, _size, mn, Arns, mn); + //write_field(ZZ,std::cout,Arns, _size, mn,mn,true); + std::cout<<"Arns1:="; + FFLAS::WriteMatrix(std::cout, ZZ, _size, mn, Arns_tmp, 2*mn); + //write_field(ZZ,std::cout,Arns_tmp, _size, mn,2*mn,true); + std::cout<<"Arns2:="; + FFLAS::WriteMatrix(std::cout, ZZ, _size, mn, Arns_tmp2, 2*mn); + //write_field(ZZ,std::cout,Arns_tmp2, _size, mn,2*mn,true); +#endif + + integer hM= (_M-1)>>1; + double *A_beta= FFLAS::fflas_new(2*mn*_ldm); + double *A_beta2 = A_beta+mn*_ldm; + tsplit.stop(); +#ifdef BENCH_RNS + if(m>1 && n>1) std::cerr<<"RNS EXTENDED double (From) - split : "<(),FFLAS::FflasTrans, FFLAS::FflasNoTrans,(int) mn*2,(int) _ldm,(int) _size, 1.0 , Arns_tmp,(int) mn*2, _crt_out.data(),(int) _ldm, 0., A_beta,(int)_ldm, + FFLAS::ParSeqHelper::Parallel()); + // FFLAS::ParSeqHelper::Parallel()); + + else // compute A_beta = Ap x M_Beta + cblas_dgemm(CblasRowMajor,CblasNoTrans, CblasNoTrans, (int)mn*2, (int)_ldm, (int)_size, 1.0 , Arns_tmp, (int)_size, _crt_out.data(), (int)_ldm, 0., A_beta,(int)_ldm); + + FFLAS::fflas_delete( Arns_tmp); + tfgemmc.stop(); +#ifdef RNS_DEBUG + std::cout<<"Abeta1:=";write_field(ZZ,std::cout,A_beta, _ldm, mn,mn,true); + std::cout<<"Abeta2:=";write_field(ZZ,std::cout,A_beta2, _ldm, mn,mn,true); +#endif + + +#ifdef BENCH_RNS + if(m>1 && n>1) std::cerr<<"RNS EXTENDED double (From) - fgemm : "<>2)+ (((k+3)%4==0)?0:1); + std::vector A0(k4<<2,0),A1(k4<<2,0),A2(k4<<2,0),A3(k4<<2,0); + integer a0,a1,a2,a3,res,res2; + mpz_t *m0,*m1,*m2,*m3; + m0= reinterpret_cast(&a0); + m1= reinterpret_cast(&a1); + m2= reinterpret_cast(&a2); + m3= reinterpret_cast(&a3); + mp_limb_t *m0_d,*m1_d,*m2_d,*m3_d; + m0_d = m0[0]->_mp_d; + m1_d = m1[0]->_mp_d; + m2_d = m2[0]->_mp_d; + m3_d = m3[0]->_mp_d; + m0[0]->_mp_alloc = m1[0]->_mp_alloc = m2[0]->_mp_alloc = m3[0]->_mp_alloc = (int) (k4*8/sizeof(mp_limb_t)); // to ensure 32 bits portability + m0[0]->_mp_size = m1[0]->_mp_size = m2[0]->_mp_size = m3[0]->_mp_size = (int) (k4*8/sizeof(mp_limb_t)); // to ensure 32 bits portability + Givaro::Timer tkroc; + tkroc.start(); + // auto sp=SPLITTER(); + // PARFOR1D(i,m,sp, + for(size_t i=0;i(&tmp); + A0[l ]= tptr[0]; + A1[l+1]= tptr[1]; + A2[l+2]= tptr[2]; + A3[l+3]= tptr[3]; + } + // see A0,A1,A2,A3 as a the gmp integers a0,a1,a2,a3 + m0[0]->_mp_d= reinterpret_cast(&A0[0]); + m1[0]->_mp_d= reinterpret_cast(&A1[0]); + m2[0]->_mp_d= reinterpret_cast(&A2[0]); + m3[0]->_mp_d= reinterpret_cast(&A3[0]); + res = a0;res+= a1;res+= a2;res+= a3; + + + for (size_t l=0;l(&tmp); + A0[l ]= tptr[0]; + A1[l+1]= tptr[1]; + A2[l+2]= tptr[2]; + A3[l+3]= tptr[3]; + } + // see A0,A1,A2,A3 as a the gmp integers a0,a1,a2,a3 + m0[0]->_mp_d= reinterpret_cast(&A0[0]); + m1[0]->_mp_d= reinterpret_cast(&A1[0]); + m2[0]->_mp_d= reinterpret_cast(&A2[0]); + m3[0]->_mp_d= reinterpret_cast(&A3[0]); + res2 = a0;res2+= a1;res2+= a2;res2+= a3; +#ifdef RNS_DEBUG + std::cout<<"res1:="<hM) + res-=_M; + if (gamma==0) + Aiter[j+i*lda]=res; + else + if (gamma==integer(1)) + Aiter[j+i*lda]+=res; + else + if (gamma==integer(-1)) + Aiter[j+i*lda]=res-Aiter[j+i*lda]; + else{ + Aiter[j+i*lda]*=gamma; + Aiter[j+i*lda]+=res; + } + + } + tkroc.stop(); +#ifdef BENCH_RNS + if(m>1 && n>1) std::cerr<<"RNS EXTENDED double (From) - Convert : "<_mp_d = m0_d; + m1[0]->_mp_d = m1_d; + m2[0]->_mp_d = m2_d; + m3[0]->_mp_d = m3_d; + m0[0]->_mp_alloc = m1[0]->_mp_alloc = m2[0]->_mp_alloc= m3[0]->_mp_alloc = 1; + m0[0]->_mp_size = m1[0]->_mp_size = m2[0]->_mp_size = m3[0]->_mp_size = 0; + FFLAS::fflas_delete( A_beta); + + +#ifdef CHECK_RNS + bool ok=true; + for (size_t i=0;is,f0,{0,g0,(0,\:0,t0,+0,=s diff --git a/fflas-ffpack/field/rns-double.h b/fflas-ffpack/field/rns-double.h index 8286c8768..19afaf2a7 100644 --- a/fflas-ffpack/field/rns-double.h +++ b/fflas-ffpack/field/rns-double.h @@ -52,7 +52,7 @@ namespace FFPACK { - /* Structure that handles rns representation given a bound and bitsize for prime moduli + /* Structure that handles rns representation with moduli stored in double FP numbers * support sign representation (i.e. the bound must be twice larger then ||A||) */ struct rns_double { @@ -79,6 +79,9 @@ namespace FFPACK { typedef rns_double_elt_ptr Element_ptr; typedef rns_double_elt_cstptr ConstElement_ptr; + + // construct an RNS basis with primes of bit-length (pbits) ensuring to represent integers lying in [0, bound[ + // Rmk: when rnsmod is set to true the RNS basis (m1,m2, ..., mk) satisfies that : m1*m2*...*mk >= bound * (m1+m2+...+mk) rns_double(const integer& bound, size_t pbits, bool rnsmod=false, long seed=time(NULL)) : _M(1), _size(0), _pbits(pbits), _mi_sum(1) { @@ -96,15 +99,11 @@ namespace FFPACK { _M*=prime; if (rnsmod) _mi_sum+=prime; } - // std::ostream_iterator out_it (std::cout,", "); - // std::cout<<"RNS basis ="; - // std::copy ( _basis.begin(), _basis.end(), out_it ); - // std::cout< - rns_double(const Vect& basis, bool rnsmod=false, long seed=time(NULL)) + rns_double(const Vect& basis, long seed=time(NULL)) : _basis(basis.begin(),basis.end()), _basisMax(basis.size()), _negbasis(basis.size()), _M(1), _size(basis.size()), _pbits(0), _mi_sum(1) { for(size_t i=0;i<_size;i++){ - //std::cout<<"basis["< hM) - A[i] -= _M; - } - } - - // reduce entries of Arns to be less than the rns basis elements - inline void rns_double_extended::reduce(size_t n, double* Arns, size_t rda, bool RNS_MAJOR) const{ - -#ifdef __FFLASFFPACK_HAVE_SSE4_1_INSTRUCTIONS - using simd = Simd; - using vect_t = typename simd::vect_t; - - if(_size % simd::vect_size == 0){ - //#pragma omp parallel for schedule(static, 256) - for(size_t i = 0 ; i < n ; i++){ - vect_t tmp1, tmp2, tmp3, v, max, basis, inv, neg; - for(size_t j = 0 ; j < _size ; j+=simd::vect_size){ - basis = simd::load(_basis.data()+j); - inv = simd::load(_invbasis.data()+j); - max = simd::load(_basisMax.data()+j); - neg = simd::load(_negbasis.data()+j); - v = simd::load(Arns+i*_size+j); - tmp2 = modSimd(v, basis, inv, neg); - tmp1 = simd::greater(tmp2, max); - tmp3 = simd::lesser(tmp2, simd::zero()); - tmp1 = simd::vand(tmp1, neg); - tmp3 = simd::vand(tmp3, basis); - tmp1 = simd::vor(tmp1, tmp3); - tmp2 = simd::add(tmp2, tmp1); - simd::store(Arns+i*_size+j, tmp2); - } - } - }else{ - //#pragma omp parallel for schedule(static, 256) - for(size_t i = 0 ; i < n ; i++){ - vect_t tmp1, tmp2, tmp3, v, max, basis, inv, neg; - size_t j = 0; - for( ; j < ROUND_DOWN(_size, simd::vect_size) ; j+=simd::vect_size){ - basis = simd::load(_basis.data()+j); - inv = simd::load(_invbasis.data()+j); - max = simd::load(_basisMax.data()+j); - neg = simd::load(_negbasis.data()+j); - v = simd::loadu(Arns+i*_size+j); - tmp2 = modSimd(v, basis, inv, neg); - tmp1 = simd::greater(tmp2, max); - tmp3 = simd::lesser(tmp2, simd::zero()); - tmp1 = simd::vand(tmp1, neg); - tmp3 = simd::vand(tmp3, basis); - tmp1 = simd::vor(tmp1, tmp3); - tmp2 = simd::add(tmp2, tmp1); - simd::storeu(Arns+i*_size+j, tmp2); - } - for( ; j < _size ; ++j){ - _field_rns[j].reduce(Arns[i*_size+j]); - } - } - } -#else - - // TODO : SIMD version - for(size_t i = 0 ; i < n ; i+= _size){ - for(size_t j = 0 ; j < _size ; ++j){ - //_field_rns.reduce(Arns+i*_size+j); - _field_rns[i].reduce(Arns[i*_size+j]); - } - } - -#endif - - } - } // FFPACK #endif // __FFLASFFPACK_field_rns_double_INL /* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ // vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s + diff --git a/fflas-ffpack/field/rns-integer.h b/fflas-ffpack/field/rns-integer.h index 95ada64d8..2ed6bdbac 100644 --- a/fflas-ffpack/field/rns-integer.h +++ b/fflas-ffpack/field/rns-integer.h @@ -34,6 +34,7 @@ #include #include "fflas-ffpack/field/rns-double.h" +#include "fflas-ffpack/field/rns-double-extended.h" namespace FFPACK { @@ -169,14 +170,26 @@ namespace FFLAS { // specialization for the fflas alloc function template<> inline FFPACK::rns_double_elt_ptr - fflas_new(const FFPACK::RNSInteger &F, const size_t m, const Alignment align){ + fflas_new(const FFPACK::RNSInteger &F, const size_t m, const Alignment align){ double *ptr=FFLAS::fflas_new(m*F.size(), align); return FFPACK::rns_double_elt_ptr(ptr,m); } template<> inline FFPACK::rns_double_elt_ptr - fflas_new(const FFPACK::RNSInteger &F, const size_t m, const size_t n, const Alignment align){ + fflas_new(const FFPACK::RNSInteger &F, const size_t m, const size_t n, const Alignment align){ + return fflas_new(F, m*n, align); + } + template<> + inline FFPACK::rns_double_elt_ptr + fflas_new(const FFPACK::RNSInteger &F, const size_t m, const Alignment align){ + double *ptr=FFLAS::fflas_new(m*F.size(), align); + return FFPACK::rns_double_elt_ptr(ptr,m); + } + + template<> + inline FFPACK::rns_double_elt_ptr + fflas_new(const FFPACK::RNSInteger &F, const size_t m, const size_t n, const Alignment align){ return fflas_new(F, m*n, align); } From a93979226eba80f76ff886a4b8700f55e9014b33 Mon Sep 17 00:00:00 2001 From: Pascal Giorgi Date: Wed, 13 Sep 2023 18:47:21 +0200 Subject: [PATCH 3/8] add benchmark for rns-double and rns-double-extended + possibility to bench against FLINT library --- benchmarks/Makefile.am | 3 +- benchmarks/benchmark-rns-double.C | 337 +++++++++++++++++++++++ fflas-ffpack/field/rns-double-extended.h | 17 ++ fflas-ffpack/field/rns-double.h | 8 + 4 files changed, 364 insertions(+), 1 deletion(-) create mode 100644 benchmarks/benchmark-rns-double.C diff --git a/benchmarks/Makefile.am b/benchmarks/Makefile.am index 43324f644..41c7e734c 100755 --- a/benchmarks/Makefile.am +++ b/benchmarks/Makefile.am @@ -33,7 +33,7 @@ endif PERFPUBLISHERFILE=benchmarks-report.xml -FFLA_BENCH = benchmark-fgemm benchmark-fgemm-rns benchmark-wino benchmark-ftrsm benchmark-fgesv benchmark-ftrsv benchmark-ftrtri benchmark-inverse benchmark-fsytrf benchmark-fsyrk benchmark-lqup benchmark-fsyr2k benchmark-pluq benchmark-charpoly benchmark-charpoly-mp benchmark-fgemm-mp benchmark-fgemv-mp benchmark-ftrsm-mp benchmark-lqup-mp benchmark-checkers benchmark-fadd-lvl2 benchmark-fdot benchmark-fgemv benchmark-quasisep benchmark-sss benchmark-storage-transpose benchmark-qscomp +FFLA_BENCH = benchmark-fgemm benchmark-fgemm-rns benchmark-wino benchmark-ftrsm benchmark-fgesv benchmark-ftrsv benchmark-ftrtri benchmark-inverse benchmark-fsytrf benchmark-fsyrk benchmark-lqup benchmark-fsyr2k benchmark-pluq benchmark-charpoly benchmark-charpoly-mp benchmark-fgemm-mp benchmark-fgemv-mp benchmark-ftrsm-mp benchmark-lqup-mp benchmark-checkers benchmark-fadd-lvl2 benchmark-fdot benchmark-fgemv benchmark-quasisep benchmark-sss benchmark-storage-transpose benchmark-qscomp benchmark-rns-double BLAS_BENCH = benchmark-sgemm$(EXEEXT) benchmark-dgemm benchmark-dtrsm LAPA_BENCH = benchmark-dtrtri benchmark-dgetri benchmark-dgetrf benchmark-dsytrf @@ -87,6 +87,7 @@ benchmark_fadd_lvl2_SOURCES = benchmark-fadd-lvl2.C benchmark_fdot_SOURCES = benchmark-fdot.C benchmark_fgemv_SOURCES = benchmark-fgemv.C benchmark_storage_transpose_SOURCES = benchmark-storage-transpose.C +benchmark_rns_double_SOURCES = benchmark-rns-double.C benchmark_sgemm_CXXFLAGS = $(AM_CXXFLAGS) -D__SGEMM__ diff --git a/benchmarks/benchmark-rns-double.C b/benchmarks/benchmark-rns-double.C new file mode 100644 index 000000000..5bdc31532 --- /dev/null +++ b/benchmarks/benchmark-rns-double.C @@ -0,0 +1,337 @@ +/* + * Copyright (C) 2023 the FFLAS-FFPACK group + * + * Written by Pascal Giorgi + * + * This file is Free Software and part of FFLAS-FFPACK. + * + * ========LICENCE======== + * This file is part of the library FFLAS-FFPACK. + * + * FFLAS-FFPACK is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with this library; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + * ========LICENCE======== + *. + */ + +//#define CHECK_RNS +//#define RNS_DEBUG +//#define BENCH_RNS + +#define __FFLASFFPACK_SEQUENTIAL + +#include "fflas-ffpack/fflas-ffpack-config.h" + + +#include +#include +#include +#include +#include +using namespace std; + +#include "fflas-ffpack/utils/timer.h" +#include "fflas-ffpack/fflas/fflas.h" +#include "fflas-ffpack/field/rns-double.h" +#include "fflas-ffpack/field/rns-double-extended.h" +#include "fflas-ffpack/field/rns-integer.h" +#include "fflas-ffpack/utils/args-parser.h" +#include "givaro/givinteger.h" +#include "givaro/modular-integer.h" +#include "givaro/zring.h" +#include + +//#define BENCH_FLINT +#ifdef BENCH_FLINT +#define __GMP_BITS_PER_MP_LIMB 64 +extern "C" { +#include "flint/longlong.h" +#include "flint/long_extras.h" +#include "flint/fmpz_mat.h" +#include "flint/fmpz.h" +#include "flint/flint.h" +} +#endif + + +#define LOOPS_TIME 1 + +template +void run_bench(size_t n, size_t primes_bits, size_t b, size_t iters, size_t seed, bool matmul=false){ + Givaro::Integer p; + FFLAS::Timer chrono; + double timeFFLASToRNS=0.,timeFFLASFromRNS=0.,timeFFLASPrecomp=0.; + double timeFFLASToRNS_ext=0.,timeFFLASFromRNS_ext=0.,timeFFLASPrecomp_ext=0.; +#ifdef BENCH_FLINT + double timeFlintToRNS=0., timeFlintToRNSnaive=0.,timeFlintFromRNS=0.,timeFlintFromRNSnaive=0.,timeFlintPrecomp=0.; +#endif + size_t loop=0; + size_t bits=b; + int logb= log((double)n)/log(2.0); + + size_t fflas_primes_bits = (primes_bits ? primes_bits : std::min(25,41-logb)); + size_t fflas_primes_bits_ext = 2* fflas_primes_bits; + // std::cout<<"FFLAS primes_bitsize= "<(end-start) < std::chrono::seconds(LOOPS_TIME)) ;loop++){ + + Givaro::Integer::random_exact_2exp(p, bits); + //nextprime(p,p); avoid p to b a prime has it is not needed here + Field F(p); + typename Field::RandIter Rand(F,0,(uint64_t)seed); + typename Field::Element_ptr A; + A= FFLAS::fflas_new(F,n,n); + + FFLAS::frand(F,Rand,n,n,A,n); + + if (matmul){ + rns_bound=p*p*(uint64_t)n; + rns_bitsize=rns_bound.bitsize(); + }else { + rns_bound=p; + rns_bitsize=rns_bound.bitsize(); + } + + int logb= log(double(rns_bitsize))/log(2.); // log of the integer bitsize of entry in A + + fflas_primes_bits = (primes_bits ? primes_bits : std::min(25,41-logb)); + fflas_primes_bits_ext = 2* fflas_primes_bits; + + + +#ifdef BENCH_FLINT + // FLINT RNS CODE // + { + fmpz_mat_t AA; + fmpz_mat_init (AA, n, n); + flint_rand_t randstate; + flint_randinit(randstate); + fmpz_mat_randbits(AA, randstate, bits); + + fmpz_comb_t comb; + fmpz_comb_temp_t comb_temp; + nmod_mat_t * mod_A; + mp_limb_t * primes; + mp_limb_t * residues; + size_t num_primes = (rns_bitsize + flint_primes_bits - 1) / flint_primes_bits; + + /* FLINT RNS initialization */ + residues = (mp_limb_t *)flint_malloc(sizeof(mp_limb_t) * num_primes); + mod_A = (nmod_mat_t *)flint_malloc(sizeof(nmod_mat_t) * num_primes); + primes = (mp_limb_t *)flint_malloc(sizeof(mp_limb_t) * num_primes); + primes[0] = n_nextprime(UWORD(1) << flint_primes_bits, 0); + nmod_mat_init(mod_A[0], AA->r, AA->c, primes[0]); + for (size_t i = 1; i < num_primes; i++) + { + primes[i] = n_nextprime(primes[i-1], 0); + nmod_mat_init(mod_A[i], AA->r, AA->c, primes[i]); + } + + chrono.clear();chrono.start(); + fmpz_comb_init(comb, primes, num_primes); + fmpz_comb_temp_init(comb_temp, comb); + chrono.stop(); + timeFlintPrecomp+=chrono.usertime(); + chrono.clear();chrono.start(); + + /* Calculate residues of AA */ + fmpz_mat_multi_mod_ui_precomp(mod_A, num_primes, AA, comb, comb_temp); + + chrono.stop(); + timeFlintToRNS+=chrono.usertime(); + chrono.clear();chrono.start(); + + + /* Calculate residue without COMB approach (NO FAST RNS) */ + //fmpz_mat_multi_mod_ui(mod_A, num_primes, AA); + // for (long i = 0; i < AA->r; i++) + // { for (long j = 0; j < AA->c; j++) + // for (long l = 0; l < num_primes; l++) + // nmod_mat_entry(mod_A[l],i,j) = fmpz_fdiv_ui(fmpz_mat_entry(AA,i,j), primes[l]); + // } + + chrono.stop(); + timeFlintToRNSnaive+=chrono.usertime(); + chrono.clear();chrono.start(); + /* Chinese remaindering */ + for (long i = 0; i < AA->r; i++) + { for (long j = 0; j < AA->c; j++) + { for (size_t k = 0; k < num_primes; k++) + residues[k] = mod_A[k]->rows[i][j]; + fmpz_multi_CRT_ui(&AA->rows[i][j], residues, comb, comb_temp, 1); + } + } + chrono.stop(); + timeFlintFromRNS+=chrono.usertime(); + + fmpz_mat_clear(AA); + flint_free(mod_A); + fmpz_comb_temp_clear(comb_temp); + fmpz_comb_clear(comb); + flint_free(residues); + flint_free(primes); + + } + //END FLINT CODE // +#endif + + // FFLAS RNS DOUBLE CODE + if (fflas_primes_bits<=41-logb){ + chrono.clear(); chrono.start(); + // construct an RNS structure and its associated Domain + FFPACK::rns_double RNS(rns_bound, fflas_primes_bits,false,seed); + typedef FFPACK::RNSInteger RnsDomain; + RnsDomain Zrns(RNS); + typename RnsDomain::Element_ptr mod_A = FFLAS::fflas_new(Zrns,n,n); + chrono.stop(); + timeFFLASPrecomp+=chrono.usertime(); + chrono.clear(); chrono.start(); + RNS.init(n,n,mod_A._ptr,mod_A._stride,A,n,p); + chrono.stop(); + timeFFLASToRNS+=chrono.usertime(); + chrono.clear(); chrono.start(); + RNS.convert(n,n,Givaro::Integer(0), A,n, mod_A._ptr,mod_A._stride); + //FFLAS::fconvert_rns(Zrns,n,n,Givaro::Integer(0),A,n,mod_A); + chrono.stop(); + timeFFLASFromRNS+=chrono.usertime(); + FFLAS::fflas_delete(mod_A); + } + + // FFLAS RNS DOUBLE EXTENDED CODE + if (fflas_primes_bits_ext<=2*(41-logb)){ + chrono.clear(); chrono.start(); + // construct an RNS structure and its associated Domain + typedef FFPACK::rns_double_extended RNSExt; + RNSExt RNS(rns_bound, fflas_primes_bits_ext,false,seed); + typedef FFPACK::RNSInteger RnsDomain; + RnsDomain Zrns(RNS); + typename RnsDomain::Element_ptr mod_A = FFLAS::fflas_new(Zrns,n,n); + chrono.stop(); + timeFFLASPrecomp_ext+=chrono.usertime(); + chrono.clear(); chrono.start(); + + RNS.init(n,n,mod_A._ptr,mod_A._stride,A,n,p); + chrono.stop(); + timeFFLASToRNS_ext+=chrono.usertime(); + chrono.clear(); chrono.start(); + //FFLAS::fconvert_rns(Zrns,n,n,Givaro::Integer(0),A,n,mod_A); + RNS.convert(n,n,Givaro::Integer(0), A,n, mod_A._ptr,mod_A._stride); + chrono.stop(); + timeFFLASFromRNS_ext+=chrono.usertime(); + FFLAS::fflas_delete(mod_A); + } + FFLAS::fflas_delete(A); + end = std::chrono::system_clock::now(); + } + +#define SPC1 10 +#define SPC2 12 +#define SPC3 33 +#define PREC 2 + + cout< Field; + + cout<<"### running RNS conversions benchmark ###"<(m,p,b,iters,seed,matmul); + else { + size_t matdim=m; + for (size_t bits=128; bits<(1<<20);bits<<=1) + run_bench(matdim, p, bits, iters,seed++,matmul); + + } + return 0; +} + + +/* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ +// vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s + diff --git a/fflas-ffpack/field/rns-double-extended.h b/fflas-ffpack/field/rns-double-extended.h index e1b581c25..0241601ba 100644 --- a/fflas-ffpack/field/rns-double-extended.h +++ b/fflas-ffpack/field/rns-double-extended.h @@ -138,6 +138,15 @@ namespace FFPACK { void precompute_cst(){ + + // Check that _pbits <= 48 + if ( _pbits > 48 ){ + std::cout<<"FFLAS Error in rns_double_extended: primes bitsize "<<_pbits<<" is too large ... aborting\n"; + std::terminate(); + } + + + _ldm = (_M.bitsize()/48) + ((_M.bitsize()%48)?1:0) ; _invbasis.resize(_size); _basisMax.resize(_size); @@ -284,6 +293,14 @@ namespace FFPACK { } // can force to reduce integer entries larger than M void precompute_cst(size_t K=0){ + + // Check that _pbits satisfies log(_M)/16* 2^(_pbits/2) * 2^16 < 2^53 (Required for correctness) + // => pbits <= 2*(41 - loglog(_M)) + if ( _pbits > 2*(41 - log(double(_M.bitsize()))/log(2.))){ + std::cout<<"FFLAS Error in rns_double_extended: primes bitsize "<<_pbits<<" is too large ... aborting\n"; + std::terminate(); + } + if (K!=0) _ldm=K; else diff --git a/fflas-ffpack/field/rns-double.h b/fflas-ffpack/field/rns-double.h index 19afaf2a7..5b2d2218b 100644 --- a/fflas-ffpack/field/rns-double.h +++ b/fflas-ffpack/field/rns-double.h @@ -141,6 +141,14 @@ namespace FFPACK { // can force to reduce integer entries larger than M void precompute_cst(size_t K=0){ + + // Check that _pbits satisfies log(_M)/16* 2^(_pbits) * 2^16 < 2^53 (Required for correctness) + // => pbits <= 41 - loglog(_M) + if ( _pbits > 41 - log(double(_M.bitsize()))/log(2.)){ + std::cout<<"FFLAS Error in rns_double: primes bitsize is too large ... aborting\n"; + std::terminate(); + } + if (K!=0) _ldm=K; else From bb261377a4edc0a9b4b873011698ae361bec69b5 Mon Sep 17 00:00:00 2001 From: Pascal Giorgi Date: Thu, 14 Sep 2023 17:27:24 +0200 Subject: [PATCH 4/8] few updates and cleanups --- benchmarks/benchmark-rns-double.C | 24 ++++++++++------------ fflas-ffpack/field/rns-double-extended.h | 12 ++--------- fflas-ffpack/field/rns-double-extended.inl | 3 +-- 3 files changed, 14 insertions(+), 25 deletions(-) diff --git a/benchmarks/benchmark-rns-double.C b/benchmarks/benchmark-rns-double.C index 5bdc31532..20be8f688 100644 --- a/benchmarks/benchmark-rns-double.C +++ b/benchmarks/benchmark-rns-double.C @@ -80,7 +80,7 @@ void run_bench(size_t n, size_t primes_bits, size_t b, size_t iters, size_t seed size_t bits=b; int logb= log((double)n)/log(2.0); - size_t fflas_primes_bits = (primes_bits ? primes_bits : std::min(25,41-logb)); + size_t fflas_primes_bits = (primes_bits ? primes_bits : 41-logb); size_t fflas_primes_bits_ext = 2* fflas_primes_bits; // std::cout<<"FFLAS primes_bitsize= "< Field; - + cout<<"### running RNS conversions benchmark ###"<> _shift); - _field_rns[i].mulin(acc,beta); - //std::cout<<"RNS precomp ("< "<< (int64_t)_crt_in[j+i*_ldm]<<" "<< (int64_t)_crt_in[j+(i+_size)*_ldm]<(mn*_size*2),accTmp; + double *Arns_tmp = FFLAS::fflas_new(mn*_size*2); double* Arns_tmp2= Arns_tmp+mn; - double mask=double((1<<_shift)+1), scaling=1./double(1<<_shift);; for (size_t i=0;i<_size;i++) for (size_t j=0;j Date: Thu, 14 Sep 2023 18:10:35 +0200 Subject: [PATCH 5/8] add test suite for rns-double --- tests/Makefile.am | 2 + tests/test-rns-double.C | 144 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 146 insertions(+) create mode 100644 tests/test-rns-double.C diff --git a/tests/Makefile.am b/tests/Makefile.am index c3c4338e7..ca18a1937 100755 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -77,6 +77,7 @@ BASIC_TESTS = \ test-fgemv \ test-nullspace \ test-storage-transpose \ + test-rns-double \ regression-check if FFLASFFPACK_PRECOMPILED @@ -162,6 +163,7 @@ test_interfaces_c_SOURCES = test-interfaces-c.c test_maxdelayeddim_SOURCES = test-maxdelayeddim.C # test_fspmv_SOURCES = test-fspmv.C test_storage_transpose_SOURCES = test-storage-transpose.C +test_rns_double_SOURCES = test-rns-double.C regression_check_SOURCES = regression-check.C test_solve_SOURCES = test-solve.C diff --git a/tests/test-rns-double.C b/tests/test-rns-double.C new file mode 100644 index 000000000..fb80ba700 --- /dev/null +++ b/tests/test-rns-double.C @@ -0,0 +1,144 @@ +/* + * Copyright (C) 2023 the FFLAS-FFPACK group + * + * Written by Pascal Giorgi + * + * This file is Free Software and part of FFLAS-FFPACK. + * + * ========LICENCE======== + * This file is part of the library FFLAS-FFPACK. + * + * FFLAS-FFPACK is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with this library; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + * ========LICENCE======== + *. + */ + +//#define CHECK_RNS +//#define RNS_DEBUG +//#define BENCH_RNS + +#define __FFLASFFPACK_SEQUENTIAL + +#include "fflas-ffpack/fflas-ffpack-config.h" + + +#include +#include +#include +#include +#include +using namespace std; + +#include "fflas-ffpack/utils/timer.h" +#include "fflas-ffpack/fflas/fflas.h" +#include "fflas-ffpack/field/rns-double.h" +#include "fflas-ffpack/field/rns-double-extended.h" +#include "fflas-ffpack/field/rns-integer.h" +#include "fflas-ffpack/utils/args-parser.h" +#include "givaro/givinteger.h" +#include "givaro/modular-integer.h" +#include "givaro/zring.h" +#include + +template +void print_rns(); + +template<> +void print_rns(){ std::cout< +void print_rns(){ std::cout< +bool run_check(size_t n, size_t primes_bits, size_t b, size_t iters, size_t seed){ + Givaro::Integer p; + + size_t bits=b; + + std::cout<<"Checking "; print_rns(); + std::cout<<" int bits="<; + Givaro::Integer::random_exact_2exp(p, bits); + Field F(p); + Field::RandIter Rand(F,0,(uint64_t)seed); + Field::Element_ptr A, Acopy; + A = FFLAS::fflas_new(F,n,n); + Acopy = FFLAS::fflas_new(F,n,n); + FFLAS::frand(F,Rand,n,n,A,n); + FFLAS::fassign(F,n,n,Acopy,n,A,n); + + Givaro::Integer rns_bound = 3*p; + + // construct an RNS structure and its associated Domain + RNS_t RNS(rns_bound, primes_bits,false,seed); + typedef FFPACK::RNSInteger RnsDomain; + RnsDomain Zrns(RNS); + typename RnsDomain::Element_ptr mod_A = FFLAS::fflas_new(Zrns,n,n); + RNS.init(n,n,mod_A._ptr,mod_A._stride,A,n,p); + RNS.convert(n,n,Givaro::Integer(0), Acopy, n, mod_A._ptr,mod_A._stride); + + + Givaro::ZRing ZZ; + if (!FFLAS::fequal(ZZ,n,n,A,n,Acopy,n)){ + std::cout<<"FAILED\n"; + return false; + } + + FFLAS::fflas_delete(mod_A); + FFLAS::fflas_delete(A); + FFLAS::fflas_delete(Acopy); + } + + std::cout<<"PASSED\n"; + return true; +} + + + +int main(int argc, char** argv){ + static size_t iters = 10 ; + static unsigned long b = 512 ; + static unsigned long p = 20 ; + static size_t m = 16 ; + static size_t seed= time(NULL); + static Argument as[] = { + { 'b', "-b B", "Set the bitsize of the matrix entries.", TYPE_INT , &b }, + { 'p', "-p B", "Set the bitsize of the RNS Prime.", TYPE_INT , &p }, + { 'm', "-m M", "Set the dimension m of the matrix.", TYPE_INT , &m }, + { 'i', "-i R", "Set minimal number of repetitions.", TYPE_INT , &iters }, + { 's', "-s S", "Set the seeding value", TYPE_INT , &seed }, + END_OF_ARGUMENTS + }; + FFLAS::parseArguments(argc,argv,as); + + /* int logb= log((double)n)/log(2.0); */ + /* size_t fflas_primes_bits = (primes_bits ? primes_bits : 41-logb); */ + /* size_t fflas_primes_bits_ext = 2* fflas_primes_bits; */ + + bool ok=true; + ok = ok and run_check (m, p,b,iters,seed); + ok = ok and run_check(m,2*p,b,iters,seed); + + return ok; +} + + +/* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ +// vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s + From af7eeeadbb8dc39a3d447184678df1ca4866ec54 Mon Sep 17 00:00:00 2001 From: Pascal Giorgi Date: Fri, 15 Sep 2023 13:27:29 +0200 Subject: [PATCH 6/8] cleanup --- fflas-ffpack/field/rns-double-extended.inl | 5 ----- 1 file changed, 5 deletions(-) diff --git a/fflas-ffpack/field/rns-double-extended.inl b/fflas-ffpack/field/rns-double-extended.inl index a4305052a..5481957ae 100644 --- a/fflas-ffpack/field/rns-double-extended.inl +++ b/fflas-ffpack/field/rns-double-extended.inl @@ -858,13 +858,10 @@ namespace FFPACK { Givaro::ModularExtended ZZ(2UL<<48);; std::cout<<"Arns:="; FFLAS::WriteMatrix(std::cout, ZZ, _size, mn, Arns, mn); - //write_field(ZZ,std::cout,Arns, _size, mn,mn,true); std::cout<<"Arns1:="; FFLAS::WriteMatrix(std::cout, ZZ, _size, mn, Arns_tmp, 2*mn); - //write_field(ZZ,std::cout,Arns_tmp, _size, mn,2*mn,true); std::cout<<"Arns2:="; FFLAS::WriteMatrix(std::cout, ZZ, _size, mn, Arns_tmp2, 2*mn); - //write_field(ZZ,std::cout,Arns_tmp2, _size, mn,2*mn,true); #endif integer hM= (_M-1)>>1; @@ -881,8 +878,6 @@ namespace FFPACK { // compute A_beta = Ap^T x M_beta FFLAS::fgemm(Givaro::ZRing(),FFLAS::FflasTrans, FFLAS::FflasNoTrans,(int) mn*2,(int) _ldm,(int) _size, 1.0 , Arns_tmp,(int) mn*2, _crt_out.data(),(int) _ldm, 0., A_beta,(int)_ldm, FFLAS::ParSeqHelper::Parallel()); - // FFLAS::ParSeqHelper::Parallel()); - else // compute A_beta = Ap x M_Beta cblas_dgemm(CblasRowMajor,CblasNoTrans, CblasNoTrans, (int)mn*2, (int)_ldm, (int)_size, 1.0 , Arns_tmp, (int)_size, _crt_out.data(), (int)_ldm, 0., A_beta,(int)_ldm); From e5816840d3b1b34397ff1739c825e7fcc6596ec6 Mon Sep 17 00:00:00 2001 From: Pascal Giorgi Date: Mon, 23 Oct 2023 16:47:44 +0200 Subject: [PATCH 7/8] fix all bugs with RNS double extended --- fflas-ffpack/fflas/fflas.h | 2 + fflas-ffpack/fflas/fflas_fadd.h | 1 + fflas-ffpack/fflas/fflas_fadd.inl | 6 +- fflas-ffpack/fflas/fflas_fscal.inl | 4 +- fflas-ffpack/fflas/fflas_level1.inl | 2 +- fflas-ffpack/fflas/fflas_transpose.h | 2 +- fflas-ffpack/field/Makefile.am | 1 + fflas-ffpack/field/rns-double-extended.h | 52 +- fflas-ffpack/field/rns-double-extended.inl | 266 ++++++--- fflas-ffpack/field/rns-double.h | 396 ++++++------- fflas-ffpack/field/rns-double.inl | 632 +++++---------------- tests/test-rns-double.C | 138 ++++- 12 files changed, 692 insertions(+), 810 deletions(-) diff --git a/fflas-ffpack/fflas/fflas.h b/fflas-ffpack/fflas/fflas.h index efd71b7e2..f545d5a84 100644 --- a/fflas-ffpack/fflas/fflas.h +++ b/fflas-ffpack/fflas/fflas.h @@ -103,6 +103,7 @@ #include "fflas_fadd.h" #include "fflas_fscal.h" #include "fflas_fassign.h" +#include "fflas_transpose.h" #include "fflas_fgemm.inl" #include "fflas_pfgemm.inl" @@ -121,6 +122,7 @@ #include "fflas_faxpy.inl" #include "fflas_fdot.inl" + //--------------------------------------------------------------------- // MultiPrecision routines //--------------------------------------------------------------------- diff --git a/fflas-ffpack/fflas/fflas_fadd.h b/fflas-ffpack/fflas/fflas_fadd.h index b693ff13a..cb0416706 100644 --- a/fflas-ffpack/fflas/fflas_fadd.h +++ b/fflas-ffpack/fflas/fflas_fadd.h @@ -64,6 +64,7 @@ namespace FFLAS { typename Field::ConstElement_ptr B, const size_t incb, typename Field::Element_ptr C, const size_t incc) { + details::fadd(F,N,A,inca,B,incb,C,incc , typename FieldTraits::category() ); } diff --git a/fflas-ffpack/fflas/fflas_fadd.inl b/fflas-ffpack/fflas/fflas_fadd.inl index 917cb0207..b7ebb9e67 100644 --- a/fflas-ffpack/fflas/fflas_fadd.inl +++ b/fflas-ffpack/fflas/fflas_fadd.inl @@ -353,7 +353,7 @@ namespace FFLAS { namespace details { ) { if (inca == 1 && incb == 1 && incc == 1) { - for (size_t i=0; i::category()); + details::fscalin(F,n,a,X,incX,typename FieldTraits::category()); } template diff --git a/fflas-ffpack/fflas/fflas_level1.inl b/fflas-ffpack/fflas/fflas_level1.inl index 0e3c94f4c..f1af3bf5f 100644 --- a/fflas-ffpack/fflas/fflas_level1.inl +++ b/fflas-ffpack/fflas/fflas_level1.inl @@ -117,7 +117,7 @@ namespace FFLAS { typename Field::ConstElement_ptr Yi = Y ; for (; Xi < X+n*incX; Xi+=incX, Yi += incY ) F.convert( *Xi , *Yi); - } + } /** fnegin * \f$x \gets - x\f$. diff --git a/fflas-ffpack/fflas/fflas_transpose.h b/fflas-ffpack/fflas/fflas_transpose.h index 3ff5b264a..5c4891655 100644 --- a/fflas-ffpack/fflas/fflas_transpose.h +++ b/fflas-ffpack/fflas/fflas_transpose.h @@ -291,7 +291,7 @@ namespace FFLAS { } } - /**************************************************************************/ + /**************************************************************************/ /* * Perfom transposition on the matrix A and store the result in matrix B. * B[j,i] = A[i,j] for 0 <= i < m, 0 <= j < n diff --git a/fflas-ffpack/field/Makefile.am b/fflas-ffpack/field/Makefile.am index 8f115d761..798740a41 100644 --- a/fflas-ffpack/field/Makefile.am +++ b/fflas-ffpack/field/Makefile.am @@ -27,6 +27,7 @@ pkgincludesubdir=$(pkgincludedir)/field RNS=rns.h \ rns.inl \ rns-double.h \ + rns-common.h\ rns-double-extended.h \ rns-double-elt.h \ rns-double.inl \ diff --git a/fflas-ffpack/field/rns-double-extended.h b/fflas-ffpack/field/rns-double-extended.h index 6ba2d317d..87e881dca 100644 --- a/fflas-ffpack/field/rns-double-extended.h +++ b/fflas-ffpack/field/rns-double-extended.h @@ -145,8 +145,6 @@ namespace FFPACK { std::terminate(); } - - _ldm = (_M.bitsize()/48) + ((_M.bitsize()%48)?1:0) ; _invbasis.resize(_size); _basisMax.resize(_size); @@ -339,9 +337,14 @@ namespace FFPACK { chrono.start(); #endif size_t l=0; - for(;l - void init(size_t m, size_t n, double* Arns, size_t rda, const T* A, size_t lda, const integer& maxA, bool RNS_MAJOR=false) const - { - init(m,n,Arns,rda,A,lda, maxA.bitsize()/16 + (maxA.bitsize()%16?1:0),RNS_MAJOR); - } + // Arns must be an array of m*n*_size + // abs(||A||) <= maxA + + + + void init(size_t m, size_t n, double* Arns, size_t rda, const integer* A, size_t lda, size_t k, bool RNS_MAJOR=false, const FFLAS::FFLAS_TRANSPOSE ta=FFLAS::FflasNoTrans) const; + void convert(size_t m, size_t n, integer gamma, integer* A, size_t lda, const double* Arns, size_t rda, bool RNS_MAJOR=false, const FFLAS::FFLAS_TRANSPOSE ta=FFLAS::FflasNoTrans) const; + void convert_bis(size_t m, size_t n, integer gamma, integer* A, size_t lda, const double* Arns, size_t rda, bool RNS_MAJOR=false, const FFLAS::FFLAS_TRANSPOSE ta=FFLAS::FflasNoTrans) const; + + // Arns must be an array of m*n*_size with abs(||A||) <= maxA + void init(size_t m, size_t n, double* Arns, size_t rda, const integer* A, size_t lda, const integer& maxA, bool RNS_MAJOR=false) const{ + init(m,n,Arns,rda,A,lda, maxA.bitsize()/16 + (maxA.bitsize()%16?1:0),RNS_MAJOR, FFLAS::FflasNoTrans); + } + void init_transpose(size_t m, size_t n, double* Arns, size_t rda, const integer* A, size_t lda, const integer& maxA, bool RNS_MAJOR=false) const{ + init(m,n,Arns,rda,A,lda, maxA.bitsize()/16 + (maxA.bitsize()%16?1:0),RNS_MAJOR, FFLAS::FflasTrans); + } + void init_transpose(size_t m, size_t n, double* Arns, size_t rda, const integer* A, size_t lda, size_t k, bool RNS_MAJOR=false) const{ + init(m,n,Arns,rda,A,lda,k,RNS_MAJOR, FFLAS::FflasTrans); + } + void convert_transpose(size_t m, size_t n, integer gamma, integer* A, size_t lda, const double* Arns, size_t rda, bool RNS_MAJOR=false) const { + convert(m,n,gamma,A,lda,Arns,rda,RNS_MAJOR,FFLAS::FflasTrans); + } + - void init(size_t m, size_t n, double* Arns, size_t rda, const integer* A, size_t lda, size_t k, bool RNS_MAJOR=false) const; - // void init_transpose(size_t m, size_t n, double* Arns, size_t rda, const integer* A, size_t lda, size_t k, bool RNS_MAJOR=false) const; - void convert(size_t m, size_t n, integer gamma, integer* A, size_t lda, const double* Arns, size_t rda, bool RNS_MAJOR=false) const; - // void convert_transpose(size_t m, size_t n, integer gamma, integer* A, size_t lda, const double* Arns, size_t rda, bool RNS_MAJOR=false) const; // reduce entries of Arns to be less than the rns basis elements void reduce(size_t n, double* Arns, size_t rda, bool RNS_MAJOR=false) const; + + template + void init(size_t m, size_t n, double* Arns, size_t rda, const RecInt::ruint* A, size_t lda, size_t k, bool RNS_MAJOR=false) const; + template + void convert(size_t m, size_t n, integer gamma, RecInt::ruint* A, size_t lda, const double* Arns, size_t rda, integer p=0,bool RNS_MAJOR=false) const; }; diff --git a/fflas-ffpack/field/rns-double-extended.inl b/fflas-ffpack/field/rns-double-extended.inl index 5481957ae..67392ffd1 100644 --- a/fflas-ffpack/field/rns-double-extended.inl +++ b/fflas-ffpack/field/rns-double-extended.inl @@ -730,105 +730,199 @@ namespace FFPACK { } - inline void rns_double_extended_V2::init(size_t m, size_t n, double* Arns, size_t rda, const integer* A, size_t lda, size_t k, bool RNS_MAJOR) const + inline void rns_double_extended_V2::init(size_t m, size_t n, double* Arns, size_t rda, const integer* A, size_t lda, size_t k, bool RNS_MAJOR, const FFLAS::FFLAS_TRANSPOSE trans) const { if (k>_ldm){ - FFPACK::failure()(__func__,__FILE__,__LINE__,"rns_struct: init (too large entry)"); - std::cerr<<"k="< rns size ("<<_size<<") kronecker size ("< Numbit(M)="<<_M.bitsize()< rns basis is too small to handle integers with 2^(16*k) values "); + std::cerr<<"with k="<(), FFLAS::FflasNoTrans,FFLAS::FflasTrans,2*_size,mn,k,1.0,_crt_in.data(),_ldm,A_beta,k,0.,Arns_tmp,mn, - FFLAS::ParSeqHelper::Parallel()); - tfgemm.stop(); -#ifdef BENCH_RNS - if(m>1 && n>1) std::cerr<<"RNS EXTENDED double (To) - fgemm : "<(mn*_size*2); + if (RNS_MAJOR==false) // Arns = _crt_in x A_beta^T + FFLAS::fgemm (ZD, FFLAS::FflasNoTrans,FFLAS::FflasTrans,2*_size,mn,k,1.0,_crt_in.data(),_ldm,A_beta,k,0.,Arns_tmp,mn,MatMulHelper); + else // Arns = A_beta x _crt_in^T + FFLAS::fgemm (ZD, FFLAS::FflasNoTrans,FFLAS::FflasTrans,mn,2*_size,k,1.0,A_beta, k, _crt_in.data(),_ldm,0.,Arns_tmp,2*_size,MatMulHelper); + + + // Reduce the pseudo reduce matrix modulo RNS basis + // reduce(mn,Arns,rda,RNS_MAJOR); + /* + NEED TO MODIFY THE CODE FOR TAKING INTO ACCOUNT RNS MAJOR MODE + */ + size_t inc = (RNS_MAJOR) ? 2*_size : 1; + size_t inc_input = (RNS_MAJOR) ? _size : 1; + size_t real_rda = (RNS_MAJOR) ? 1 : mn; + double *Arns_tmp_high = (RNS_MAJOR) ? Arns_tmp+_size: Arns_tmp+ mn*_size; + + + // Givaro::ZRing ZZ; + // int64_t * ArnsConverted = new int64_t[2*mn*_size]; + + + + // std::cout<<"basis=["; + // for(size_t i=0;i<_size-1;i++) + // std::cout<<(int64_t)_basis[i]<<","; + // std::cout<<(int64_t)_basis[_size-1]<<"]\n"; + + // std::cout<<"shift="<(),2*_size*m*n,ArnsConverted,1,Arns_tmp,1); + // FFLAS::WriteMatrix(std::cout<<"AA=",ZZ,m*n,2*_size, ArnsConverted, 2*_size,FFLAS::FflasSageMath)< scaling by 2^shift allow a proper modular reduction + + + + for(size_t i=0;i<_size;i++){ + // PG: it seems that the fscalin is buggy when called below -> use fscal + faddin instead of fscalin + fadd + //FFLAS::fscalin (_field_rns[i],mn, double(1<<_shift), Arns_tmp_high+i*real_rda, inc); + //FFLAS::fadd(_field_rns[i],mn,Arns_tmp+i*real_rda,inc, Arns_tmp_high+i*real_rda,inc, Arns+i*real_rda,inc_input); + + FFLAS::fscal (_field_rns[i],mn, double(1<<_shift), Arns_tmp_high+i*real_rda, inc, Arns+i*real_rda, inc_input); // This two lines cannot be replaced with an axpy + FFLAS::faddin(_field_rns[i],mn, Arns_tmp+i*real_rda,inc, Arns+i*real_rda,inc_input); + + // FFLAS::fconvert(Givaro::ZRing(),_size*m*n,ArnsConverted,1,Arns,1); + // FFLAS::WriteMatrix(std::cout<<"AA"<(),_size*m*n,ArnsConverted,1,Arns,1); + // FFLAS::WriteMatrix(std::cout<<"A=",ZZ,m*n,_size, ArnsConverted, _size,FFLAS::FflasSageMath)<1 && n>1) std::cerr<<"RNS EXTENDED double (To) - Reduce : "< ZD; + Givaro::ZRing ZZ; + if (mn) { +#ifdef CHECK_RNS + integer* Acopy=new integer[m*n]; + FFLAS::fassign(ZZ,m,n,Acopy,n,A,lda); #endif - FFLAS::fflas_delete( A_beta); - FFLAS::fflas_delete( Arns_tmp); + double *A_beta = FFLAS::fflas_new(mn2*_ldm); + + // Split Arns into two matrix (lower half bit and upper half bit) + double *Arns_low = FFLAS::fflas_new(mn2*_size); + + // Using Helper for potential parallelism -> need to be activated by hand + //using ParallelStrategy= FFLAS::ParSeqHelper::Parallel; + using ParallelStrategy= FFLAS::ParSeqHelper::Sequential; + FFLAS::MMHelper, FFLAS::MMHelperAlgo::Winograd> MatMulHelper (ZD, -1, ParallelStrategy()); + + if (RNS_MAJOR==false) {// compute A_beta = Arns^T x M_beta + double* Arns_high = Arns_low+mn; + for (size_t i=0;i<_size;i++) + for (size_t j=0;j> _shift); + } + FFLAS::fgemm(Givaro::ZRing(),FFLAS::FflasTrans, FFLAS::FflasNoTrans, mn2, _ldm, _size, 1.0 , Arns_low, mn2, _crt_out.data(), _ldm, 0., A_beta,_ldm, MatMulHelper); + } + else { // compute A_beta = Arns x M_Beta + double* Arns_high = Arns_low+mn*_size; + for (size_t i=0;i> _shift); + } + FFLAS::fgemm(Givaro::ZRing(),FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, mn2, _ldm, _size, 1.0 , Arns_low, _size, _crt_out.data(), _ldm, 0., A_beta, _ldm, MatMulHelper); + } + + // std::cout<<"****\n RNS partial reconstruction:\n"; + // std::cout<<"shift="<<_shift<<"\n"; + // Givaro::ZRing ZZ; + // int64_t * ArnsInt = new int64_t[mn2*_size]; + // FFLAS::fconvert(Givaro::ZRing(),_size,mn2,ArnsInt,mn2,Arns_low,mn2 ); + // std::cout<<"Arns="; + // FFLAS::WriteMatrix(std::cout,ZZ,_size, mn2, ArnsInt, mn2,FFLAS::FflasSageMath)<(),m, n, A, n)<(),m, n, A, n)< "<> _shift); } - + #ifdef RNS_DEBUG Givaro::ModularExtended ZZ(2UL<<48);; std::cout<<"Arns:="; @@ -889,7 +983,7 @@ namespace FFPACK { #endif -#ifdef BENCH_RNS +#ifdef BENCH_RNS if(m>1 && n>1) std::cerr<<"RNS EXTENDED double (From) - fgemm : "< ModField; - - std::vector> _basis; // the rns moduli (mi) - std::vector> _basisMax; // (mi-1) - std::vector> _negbasis; // (-mi) - std::vector> _invbasis; // the inverse of rns moduli (1/mi) - std::vector _field_rns; // the associated prime field for each mi - integer _M; // the product of the mi's - std::vector _Mi; // _M/mi - std::vector _MMi; // (_Mi)^(-1) mod mi - std::vector _crt_in; // 2^(16*j) mod mi - std::vector _crt_out; // (_Mi._MMi) written in base 2^16 - size_t _size; // the size of the rns basis (number of mi's) - size_t _pbits; // the size in bit of the mi's - size_t _ldm; // log[2^16](_M) - integer _mi_sum; // the product of the mi's - - typedef double BasisElement; - typedef rns_double_elt Element; - typedef rns_double_elt_ptr Element_ptr; - typedef rns_double_elt_cstptr ConstElement_ptr; - - - // construct an RNS basis with primes of bit-length (pbits) ensuring to represent integers lying in [0, bound[ - // Rmk: when rnsmod is set to true the RNS basis (m1,m2, ..., mk) satisfies that : m1*m2*...*mk >= bound * (m1+m2+...+mk) - rns_double(const integer& bound, size_t pbits, bool rnsmod=false, long seed=time(NULL)) - : _M(1), _size(0), _pbits(pbits), _mi_sum(1) - { - integer::seeding(seed); - Givaro::IntPrimeDom IPD; - integer prime; - while (_M < bound*_mi_sum) { - _basis.resize(_size+1); - do { - integer::random_exact_2exp(prime, _pbits-1); - IPD.nextprimein(prime); - } while (_M%prime == 0); - _basis[_size]=prime; - _size++; - _M*=prime; - if (rnsmod) _mi_sum+=prime; - } - precompute_cst(); - } - - // construct an RNS basis with (size) primes of bit-length (pbits) - rns_double(size_t pbits, size_t size, long seed=time(NULL)) - : _M(1), _size(size), _pbits(pbits), _mi_sum(1) - { - integer::seeding(seed); - Givaro::IntPrimeDom IPD; - integer prime; - _basis.resize(size); - _negbasis.resize(size); - _basisMax.resize(size); - for(size_t i = 0 ; i < _size ; ++i){ - integer::random_exact_2exp(prime, _pbits-1); - IPD.nextprimein(prime); - _basis[i]=prime; - _basisMax[i] = prime-1; - _negbasis[i] = 0-prime; - _M*=prime; - } - precompute_cst(); - } - - // construct an RNS basis from a vector of relatively prime numbers (basis) - template - rns_double(const Vect& basis, long seed=time(NULL)) - : _basis(basis.begin(),basis.end()), _basisMax(basis.size()), _negbasis(basis.size()), _M(1), _size(basis.size()), _pbits(0), _mi_sum(1) - { - for(size_t i=0;i<_size;i++){ - _M*=_basis[i]; - _pbits=std::max(_pbits, integer(_basis[i]).bitsize()); - } - precompute_cst(); - } - - rns_double(const RNSIntegerMod& basis, bool rnsmod=false, long seed=time(NULL)) { - - } - - // can force to reduce integer entries larger than M - void precompute_cst(size_t K=0){ - - // Check that _pbits satisfies log(_M)/16* 2^(_pbits) * 2^16 < 2^53 (Required for correctness) - // => pbits <= 41 - loglog(_M) - if ( _pbits > 41 - log(double(_M.bitsize()))/log(2.)){ - std::cout<<"FFLAS Error in rns_double: primes bitsize is too large ... aborting\n"; - std::terminate(); - } - - if (K!=0) - _ldm=K; - else - _ldm = (_M.bitsize()/16) + ((_M.bitsize()%16)?1:0) ; - _invbasis.resize(_size); - _field_rns.resize(_size); - _Mi.resize(_size); - _MMi.resize(_size); - _basisMax.resize(_size); - _negbasis.resize(_size); - _crt_in.resize(_size*_ldm); - _crt_out.resize(_size*_ldm); - for (size_t i=0;i<_size;i++){ - _invbasis[i] = 1./_basis[i]; - _basisMax[i] = _basis[i]-1; - _negbasis[i] = 0-_basis[i]; - _field_rns[i] = ModField(_basis[i]); - _Mi[i] = _M/(uint64_t)_basis[i]; - _field_rns[i].init(_MMi[i], _Mi[i] % (double)_basis[i]); - _field_rns[i].invin(_MMi[i]); - integer tmp= _Mi[i]*(uint64_t)_MMi[i]; - const mpz_t* m0 = reinterpret_cast(&tmp); - const uint16_t* m0_ptr = reinterpret_cast(m0[0]->_mp_d); - size_t maxs=std::min(_ldm,(tmp.size())*sizeof(mp_limb_t)/2);// to ensure 32 bits portability - - size_t l=0; + /* Structure that handles rns representation with moduli stored in double FP numbers + * support sign representation (i.e. the bound must be twice larger then ||A||) + */ + struct rns_double { + typedef Givaro::Integer integer; + typedef Givaro::Modular ModField; + + std::vector> _basis; // the rns moduli (mi) + std::vector> _basisMax; // (mi-1) + std::vector> _negbasis; // (-mi) + std::vector> _invbasis; // the inverse of rns moduli (1/mi) + std::vector _field_rns; // the associated prime field for each mi + integer _M; // the product of the mi's + std::vector _Mi; // _M/mi + std::vector _MMi; // (_Mi)^(-1) mod mi + std::vector _crt_in; // 2^(16*j) mod mi + std::vector _crt_out; // (_Mi._MMi) written in base 2^16 + size_t _size; // the size of the rns basis (number of mi's) + size_t _pbits; // the size in bit of the mi's + size_t _ldm; // log[2^16](_M) + integer _mi_sum; // the product of the mi's + + typedef double BasisElement; + typedef rns_double_elt Element; + typedef rns_double_elt_ptr Element_ptr; + typedef rns_double_elt_cstptr ConstElement_ptr; + + + // construct an RNS basis with primes of bit-length (pbits) ensuring to represent integers lying in [0, bound[ + // Rmk: when rnsmod is set to true the RNS basis (m1,m2, ..., mk) satisfies that : m1*m2*...*mk >= bound * (m1+m2+...+mk) + rns_double(const integer& bound, size_t pbits, bool rnsmod=false, long seed=time(NULL)) + : _M(1), _size(0), _pbits(pbits), _mi_sum(1) + { + integer::seeding(seed); + Givaro::IntPrimeDom IPD; + integer prime; + while (_M < bound*_mi_sum) { + _basis.resize(_size+1); + do { + integer::random_exact_2exp(prime, _pbits-1); + IPD.nextprimein(prime); + } while (_M%prime == 0); + _basis[_size]=prime; + _size++; + _M*=prime; + if (rnsmod) _mi_sum+=prime; + } + precompute_cst(); + } + + // construct an RNS basis with (size) primes of bit-length (pbits) + rns_double(size_t pbits, size_t size, long seed=time(NULL)) + : _M(1), _size(size), _pbits(pbits), _mi_sum(1) + { + integer::seeding(seed); + Givaro::IntPrimeDom IPD; + integer prime; + _basis.resize(size); + _negbasis.resize(size); + _basisMax.resize(size); + for(size_t i = 0 ; i < _size ; ++i){ + integer::random_exact_2exp(prime, _pbits-1); + IPD.nextprimein(prime); + _basis[i]=prime; + _basisMax[i] = prime-1; + _negbasis[i] = 0-prime; + _M*=prime; + } + precompute_cst(); + } + + // construct an RNS basis from a vector of relatively prime numbers (basis) + template + rns_double(const Vect& basis, long seed=time(NULL)) + : _basis(basis.begin(),basis.end()), _basisMax(basis.size()), _negbasis(basis.size()), _M(1), _size(basis.size()), _pbits(0), _mi_sum(1) + { + for(size_t i=0;i<_size;i++){ + _M*=_basis[i]; + _pbits=std::max(_pbits, integer(_basis[i]).bitsize()); + } + precompute_cst(); + } + + rns_double(const RNSIntegerMod& basis, bool rnsmod=false, long seed=time(NULL)) { + + } + + // can force to reduce integer entries larger than M + void precompute_cst(size_t K=0){ + + // Check that _pbits satisfies log(_M)/16* 2^(_pbits) * 2^16 < 2^53 (Required for correctness) + // => pbits <= 41 - loglog(_M) + if ( _pbits > 41 - log(double(_M.bitsize()))/log(2.)){ + std::cout<<"FFLAS Error in rns_double: primes bitsize is too large ... aborting\n"; + std::terminate(); + } + + if (K!=0) + _ldm=K; + else + _ldm = (_M.bitsize()/16) + ((_M.bitsize()%16)?1:0) ; + _invbasis.resize(_size); + _field_rns.resize(_size); + _Mi.resize(_size); + _MMi.resize(_size); + _basisMax.resize(_size); + _negbasis.resize(_size); + _crt_in.resize(_size*_ldm); + _crt_out.resize(_size*_ldm); + for (size_t i=0;i<_size;i++){ + _invbasis[i] = 1./_basis[i]; + _basisMax[i] = _basis[i]-1; + _negbasis[i] = 0-_basis[i]; + _field_rns[i] = ModField(_basis[i]); + _Mi[i] = _M/(uint64_t)_basis[i]; + _field_rns[i].init(_MMi[i], _Mi[i] % (double)_basis[i]); + _field_rns[i].invin(_MMi[i]); + integer tmp= _Mi[i]*(uint64_t)_MMi[i]; + const mpz_t* m0 = reinterpret_cast(&tmp); + const uint16_t* m0_ptr = reinterpret_cast(m0[0]->_mp_d); + size_t maxs=std::min(_ldm,(tmp.size())*sizeof(mp_limb_t)/2);// to ensure 32 bits portability + + size_t l=0; #ifdef __FFLASFFPACK_HAVE_LITTLE_ENDIAN - for(;l - void init(size_t m, size_t n, double* Arns, size_t rda, const T* A, size_t lda, const integer& maxA, bool RNS_MAJOR=false) const - { - init(m,n,Arns,rda,A,lda, maxA.bitsize()/16 + (maxA.bitsize()%16?1:0),RNS_MAJOR); - } - - void init(size_t m, size_t n, double* Arns, size_t rda, const integer* A, size_t lda, size_t k, bool RNS_MAJOR=false) const; - void init_transpose(size_t m, size_t n, double* Arns, size_t rda, const integer* A, size_t lda, size_t k, bool RNS_MAJOR=false) const; - void convert(size_t m, size_t n, integer gamma, integer* A, size_t lda, const double* Arns, size_t rda, bool RNS_MAJOR=false) const; - void convert_transpose(size_t m, size_t n, integer gamma, integer* A, size_t lda, const double* Arns, size_t rda, bool RNS_MAJOR=false) const; - - // reduce entries of Arns to be less than the rns basis elements - void reduce(size_t n, double* Arns, size_t rda, bool RNS_MAJOR=false) const; + // rda is the distance between two consecutive residues : i.e. A_ij mod mk and A_ij mod m_(k+1) -> it is either 1 (RNS_MAJOR) or mxn (NOT RNS_MAJOR) + void init(size_t m, size_t n, double* Arns, size_t rda, const integer* A, size_t lda, size_t k, bool RNS_MAJOR=false, const FFLAS::FFLAS_TRANSPOSE ta=FFLAS::FflasNoTrans) const; + void convert(size_t m, size_t n, integer gamma, integer* A, size_t lda, const double* Arns, size_t rda, bool RNS_MAJOR=false, const FFLAS::FFLAS_TRANSPOSE ta=FFLAS::FflasNoTrans) const; - template - void init(size_t m, size_t n, double* Arns, size_t rda, const RecInt::ruint* A, size_t lda, size_t k, bool RNS_MAJOR=false) const; - template - void convert(size_t m, size_t n, integer gamma, RecInt::ruint* A, size_t lda, const double* Arns, size_t rda, integer p=0,bool RNS_MAJOR=false) const; + // Arns must be an array of m*n*_size with abs(||A||) <= maxA + void init(size_t m, size_t n, double* Arns, size_t rda, const integer* A, size_t lda, const integer& maxA, bool RNS_MAJOR=false) const{ + init(m,n,Arns,rda,A,lda, maxA.bitsize()/16 + (maxA.bitsize()%16?1:0),RNS_MAJOR, FFLAS::FflasNoTrans); + } + void init_transpose(size_t m, size_t n, double* Arns, size_t rda, const integer* A, size_t lda, const integer& maxA, bool RNS_MAJOR=false) const{ + init(m,n,Arns,rda,A,lda, maxA.bitsize()/16 + (maxA.bitsize()%16?1:0),RNS_MAJOR, FFLAS::FflasTrans); + } + void init_transpose(size_t m, size_t n, double* Arns, size_t rda, const integer* A, size_t lda, size_t k, bool RNS_MAJOR=false) const{ + init(m,n,Arns,rda,A,lda,k,RNS_MAJOR, FFLAS::FflasTrans); + } + void convert_transpose(size_t m, size_t n, integer gamma, integer* A, size_t lda, const double* Arns, size_t rda, bool RNS_MAJOR=false) const { + convert(m,n,gamma,A,lda,Arns,rda,RNS_MAJOR,FFLAS::FflasTrans); + } - }; // end of struct rns_double + // reduce entries of Arns to be less than the rns basis elements + void reduce(size_t n, double* Arns, size_t rda, bool RNS_MAJOR=false) const; + template + void init(size_t m, size_t n, double* Arns, size_t rda, const RecInt::ruint* A, size_t lda, size_t k, bool RNS_MAJOR=false) const; + template + void convert(size_t m, size_t n, integer gamma, RecInt::ruint* A, size_t lda, const double* Arns, size_t rda, integer p=0,bool RNS_MAJOR=false) const; - template - class rnsRandIter { - std::vector _RNS_rand; - const RNS& _domain; - - public: - rnsRandIter(const RNS& R, uint64_t seed=0) - : _domain(R) { - for(const auto& F : R._field_rns) - _RNS_rand.emplace_back(F,seed); - } + }; // end of struct rns_double - /** RNS ring Element random assignement. - * Element is supposed to be initialized - * @return random ring Element - */ - typename RNS::Element& random(typename RNS::Element& elt) const { - auto coefficient(elt._ptr); - for(auto & iter : _RNS_rand) { - iter.random( *coefficient ); - coefficient += elt._stride; - } - return elt; - } - typename RNS::Element& operator()(typename RNS::Element& elt) const { - return this->random(elt); - } - typename RNS::Element operator()() const { - typename RNS::Element tmp; _domain.init(tmp); - return this->operator()(tmp); - } - typename RNS::Element random() const { - return this->operator()(); - } + template + class rnsRandIter { + std::vector _RNS_rand; + const RNS& _domain; - const RNS& ring() const { return _domain; } + public: + rnsRandIter(const RNS& R, uint64_t seed=0) + : _domain(R) { + for(const auto& F : R._field_rns) + _RNS_rand.emplace_back(F,seed); + } - }; + /** RNS ring Element random assignement. + * Element is supposed to be initialized + * @return random ring Element + */ + typename RNS::Element& random(typename RNS::Element& elt) const { + auto coefficient(elt._ptr); + for(auto & iter : _RNS_rand) { + iter.random( *coefficient ); + coefficient += elt._stride; + } + return elt; + } + + typename RNS::Element& operator()(typename RNS::Element& elt) const { + return this->random(elt); + } + + typename RNS::Element operator()() const { + typename RNS::Element tmp; _domain.init(tmp); + return this->operator()(tmp); + } + typename RNS::Element random() const { + return this->operator()(); + } + + const RNS& ring() const { return _domain; } + + }; } // end of namespace FFPACK @@ -266,10 +272,10 @@ namespace FFPACK { #include "rns-double-recint.inl" namespace FFLAS { - template<> - inline void fflas_delete (FFPACK::rns_double_elt_ptr A) {FFLAS::fflas_delete( A._ptr);} - template<> - inline void fflas_delete (FFPACK::rns_double_elt_cstptr A) {delete[] A._ptr;} + template<> + inline void fflas_delete (FFPACK::rns_double_elt_ptr A) {FFLAS::fflas_delete( A._ptr);} + template<> + inline void fflas_delete (FFPACK::rns_double_elt_cstptr A) {delete[] A._ptr;} } diff --git a/fflas-ffpack/field/rns-double.inl b/fflas-ffpack/field/rns-double.inl index 7f8cd841a..deec21b4f 100644 --- a/fflas-ffpack/field/rns-double.inl +++ b/fflas-ffpack/field/rns-double.inl @@ -29,516 +29,182 @@ #define __FFLASFFPACK_field_rns_double_INL #include "fflas-ffpack/fflas/fflas_freduce.h" +#include "fflas-ffpack/field/rns-common.h" -namespace FFPACK { - // Arns must be an array of m*n*_size - // abs(||A||) < 2^(16k) - inline void rns_double::init(size_t m, size_t n, double* Arns, size_t rda, const integer* A, size_t lda, size_t k, bool RNS_MAJOR) const - { - if (k>_ldm){ - FFPACK::failure()(__func__,__FILE__,__LINE__,"rns_double [init] -> rns basis is too small to handle integers with 2^(16*k) values "); - std::cerr<<"with k="<(), FFLAS::FflasNoTrans,FFLAS::FflasTrans,_size,mn,k,1.0,_crt_in.data(),_ldm,A_beta,k,0.,Arns,rda, - FFLAS::ParSeqHelper::Parallel()); -#else - cblas_dgemm(CblasRowMajor,CblasNoTrans,CblasTrans,(int)_size,(int)mn,(int)k,1.0,_crt_in.data(),(int)_ldm,A_beta,(int)k,0.,Arns,(int)rda); -#endif - tfgemm.stop(); - //if(m>1 && n>1) std::cerr<<"fgemm : "<(), FFLAS::FflasNoTrans,FFLAS::FflasTrans,mn,_size,k,1.0,A_beta, k, _crt_in.data(),_ldm,0.,Arns,_size, - FFLAS::ParSeqHelper::Parallel()); -#else - cblas_dgemm(CblasRowMajor,CblasNoTrans,CblasTrans,(int)mn,(int)_size,(int)k,1.0,A_beta,(int)k,_crt_in.data(),(int)_ldm,0.,Arns,(int)_size); -#endif - } - Givaro::Timer tred; tred.start(); - - reduce(mn,Arns,rda,RNS_MAJOR); - tred.stop(); - //if(m>1 && n>1) std::cerr<<"Reduce : "<_ldm){ + FFPACK::failure()(__func__,__FILE__,__LINE__,"rns_double [init] -> rns basis is too small to handle integers with 2^(16*k) values "); + std::cerr<<"with k="<(mn*k); - const integer* Aiter=A; - // split A into A_beta according to a Kronecker transform in base 2^16 - for(size_t j=0;j(Aiter+j+i*lda); - const uint16_t* m0_ptr = reinterpret_cast(m0[0]->_mp_d); - size_t l=0; - //size_t maxs=std::min(k,(Aiter[j+i*lda].size())<<2); - size_t maxs=std::min(k,(Aiter[j+i*lda].size())*sizeof(mp_limb_t)/2); // to ensure 32 bits portability -#ifdef __FFLASFFPACK_HAVE_LITTLE_ENDIAN - if (m0[0]->_mp_size >= 0) - for (;l_mp_size >= 0) - for (;l(), FFLAS::FflasNoTrans,FFLAS::FflasTrans,_size,mn,k,1.0,_crt_in.data(),_ldm, A_beta, k, 0.,Arns,rda, - FFLAS::ParSeqHelper::Parallel()); -#else - cblas_dgemm(CblasRowMajor,CblasNoTrans,CblasTrans,(int)_size,(int)mn,(int)k,1.0,_crt_in.data(),(int)_ldm,A_beta,(int)k,0.,Arns,(int)rda); -#endif - } - else { - // Arns = A_beta x _crt_in^T -#ifndef ENABLE_CHECKER_fgemm - FFLAS::fgemm (Givaro::ZRing(), FFLAS::FflasNoTrans,FFLAS::FflasTrans,mn,_size,k,1.0,A_beta, k, _crt_in.data(),_ldm,0.,Arns,_size, - FFLAS::ParSeqHelper::Parallel()); -#else - cblas_dgemm(CblasRowMajor,CblasNoTrans,CblasTrans,(int)mn,(int)_size,(int)k,1.0,A_beta,(int)k,_crt_in.data(),(int)_ldm,0.,Arns,(int)_size); -#endif - } - reduce(mn,Arns,rda,RNS_MAJOR); - - FFLAS::fflas_delete( A_beta); + inline void rns_double::convert(size_t m, size_t n, integer gamma, integer* A, size_t lda, + const double* Arns, size_t rda, bool RNS_MAJOR, const FFLAS::FFLAS_TRANSPOSE trans) const + { + const size_t mn= m*n; + Givaro::ZRing ZD; + Givaro::ZRing ZZ; + if (mn) { #ifdef CHECK_RNS - bool ok=true; - for (size_t i=0;i(mn*_ldm); -#endif + // Using Helper for potential parallelism -> need to be activated by hand + //using ParallelStrategy= FFLAS::ParSeqHelper::Parallel; + using ParallelStrategy= FFLAS::ParSeqHelper::Sequential; + FFLAS::MMHelper, FFLAS::MMHelperAlgo::Winograd> MatMulHelper (ZD, -1, ParallelStrategy()); - integer hM= (_M-1)>>1; - double *A_beta= FFLAS::fflas_new(mn*_ldm); - Givaro::Timer tfgemmc;tfgemmc.start(); - if (RNS_MAJOR==false) {// compute A_beta = Ap^T x M_beta -#ifndef ENABLE_CHECKER_fgemm - FFLAS::fgemm(Givaro::ZRing(),FFLAS::FflasTrans, FFLAS::FflasNoTrans, mn, _ldm, _size, 1.0 , Arns, rda, _crt_out.data(), _ldm, 0., A_beta,_ldm, FFLAS::ParSeqHelper::Parallel()); -#else - cblas_dgemm(CblasRowMajor,CblasTrans, CblasNoTrans, (int)mn, (int)_ldm, (int)_size, 1.0 , Arns, (int)rda, _crt_out.data(), (int)_ldm, 0., A_beta,(int)_ldm); -#endif - } - else {// compute A_beta = Ap x M_Beta -#ifndef ENABLE_CHECKER_fgemm - FFLAS::fgemm(Givaro::ZRing(),FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, mn, _ldm, _size, 1.0 , Arns, _size, _crt_out.data(), _ldm, 0., A_beta, _ldm, FFLAS::ParSeqHelper::Parallel()); -#else - cblas_dgemm(CblasRowMajor,CblasNoTrans, CblasNoTrans, (int)mn, (int)_ldm, (int)_size, 1.0 , Arns, (int)_size, _crt_out.data(), (int)_ldm, 0., A_beta,(int)_ldm); -#endif - } - tfgemmc.stop(); - //if(m>1 && n>1) std::cerr<<"fgemm Convert : "<>2)+ (((k+3)%4==0)?0:1); - std::vector A0(k4<<2,0),A1(k4<<2,0),A2(k4<<2,0),A3(k4<<2,0); - integer a0,a1,a2,a3,res; - mpz_t *m0,*m1,*m2,*m3; - m0= reinterpret_cast(&a0); - m1= reinterpret_cast(&a1); - m2= reinterpret_cast(&a2); - m3= reinterpret_cast(&a3); - mp_limb_t *m0_d,*m1_d,*m2_d,*m3_d; - m0_d = m0[0]->_mp_d; - m1_d = m1[0]->_mp_d; - m2_d = m2[0]->_mp_d; - m3_d = m3[0]->_mp_d; - m0[0]->_mp_alloc = m1[0]->_mp_alloc = m2[0]->_mp_alloc = m3[0]->_mp_alloc = (int) (k4*8/sizeof(mp_limb_t)); // to ensure 32 bits portability - m0[0]->_mp_size = m1[0]->_mp_size = m2[0]->_mp_size = m3[0]->_mp_size = (int) (k4*8/sizeof(mp_limb_t)); // to ensure 32 bits portability - Givaro::Timer tkroc; - tkroc.start(); - // auto sp=SPLITTER(); - // PARFOR1D(i,m,sp, - for(size_t i=0;i(&tmp); -#ifdef __FFLASFFPACK_HAVE_LITTLE_ENDIAN - A0[l ]= tptr[0]; - A1[l+1]= tptr[1]; - A2[l+2]= tptr[2]; - A3[l+3]= tptr[3]; -#else - A0[l ^ ((sizeof(mp_limb_t)/2U) - 1U)] = tptr[3]; - A1[(l+1) ^ ((sizeof(mp_limb_t)/2U) - 1U)] = tptr[2]; - A2[(l+2) ^ ((sizeof(mp_limb_t)/2U) - 1U)] = tptr[1]; - A3[(l+3) ^ ((sizeof(mp_limb_t)/2U) - 1U)] = tptr[0]; -#endif - } - // see A0,A1,A2,A3 as a the gmp integers a0,a1,a2,a3 - m0[0]->_mp_d= reinterpret_cast(&A0[0]); - m1[0]->_mp_d= reinterpret_cast(&A1[0]); - m2[0]->_mp_d= reinterpret_cast(&A2[0]); - m3[0]->_mp_d= reinterpret_cast(&A3[0]); - res = a0;res+= a1;res+= a2;res+= a3; - res%=_M; + if (RNS_MAJOR==false) // compute A_beta = Ap^T x M_beta + FFLAS::fgemm(Givaro::ZRing(),FFLAS::FflasTrans, FFLAS::FflasNoTrans, mn, _ldm, _size, 1.0 , Arns, rda, _crt_out.data(), _ldm, 0., A_beta,_ldm, MatMulHelper); + else // compute A_beta = Ap x M_Beta + FFLAS::fgemm(Givaro::ZRing(),FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, mn, _ldm, _size, 1.0 , Arns, _size, _crt_out.data(), _ldm, 0., A_beta, _ldm, MatMulHelper); - // get the correct result according to the expected sign of A - if (res>hM) - res-=_M; - if (gamma==0) - Aiter[j+i*lda]=res; - else - if (gamma==integer(1)) - Aiter[j+i*lda]+=res; - else - if (gamma==integer(-1)) - Aiter[j+i*lda]=res-Aiter[j+i*lda]; - else{ - Aiter[j+i*lda]*=gamma; - Aiter[j+i*lda]+=res; - } - } - tkroc.stop(); - //if(m>1 && n>1) std::cerr<<"Kronecker Convert : "<_mp_d = m0_d; - m1[0]->_mp_d = m1_d; - m2[0]->_mp_d = m2_d; - m3[0]->_mp_d = m3_d; - m0[0]->_mp_alloc = m1[0]->_mp_alloc = m2[0]->_mp_alloc= m3[0]->_mp_alloc = 1; - m0[0]->_mp_size = m1[0]->_mp_size = m2[0]->_mp_size = m3[0]->_mp_size = 0; - FFLAS::fflas_delete( A_beta); + //FFLAS::WriteMatrix(std::cout,Givaro::ZRing(),m,n, A, lda, FFLAS::FflasSageMath)<(),1,_size, _basis.data(), _size, FFLAS::FflasSageMath)<(),mn,_size, Arns, rda,FFLAS::FflasSageMath)<(),_size, mn, Arns, rda,FFLAS::FflasSageMath)<>1; - double *A_beta= FFLAS::fflas_new(mn*_ldm); - - if (RNS_MAJOR==false){ - // compute A_beta = Ap^T x M_beta -#ifndef ENABLE_CHECKER_fgemm - FFLAS::fgemm(Givaro::ZRing(),FFLAS::FflasTrans, FFLAS::FflasNoTrans, mn, _ldm, _size, 1.0 , Arns, rda, _crt_out.data(), _ldm, 0., A_beta, _ldm, FFLAS::ParSeqHelper::Parallel()); -#else - cblas_dgemm(CblasRowMajor,CblasTrans, CblasNoTrans,(int) mn,(int) _ldm,(int) _size, 1.0 , Arns,(int) rda, _crt_out.data(),(int) _ldm, 0., A_beta,(int)_ldm); -#endif - } - else { // compute A_beta = Ap x M_Beta -#ifndef ENABLE_CHECKER_fgemm - FFLAS::fgemm(Givaro::ZRing(),FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, mn, _ldm, _size, 1.0 , Arns, rda, _crt_out.data(), _ldm, 0., A_beta, _ldm, FFLAS::ParSeqHelper::Parallel()); -#else - cblas_dgemm(CblasRowMajor,CblasNoTrans, CblasNoTrans, (int)mn, (int)_ldm, (int)_size, 1.0 , Arns, (int)rda, _crt_out.data(), (int)_ldm, 0., A_beta,(int)_ldm); + RNS_COMMON::check_rns_to_integers(trans, m, n, Arns, rda, A, lda, gamma, Acopy, n, _basis, RNS_MAJOR); + FFLAS::fflas_delete(Acopy); #endif - } - // compute A using inverse Kronecker transform of A_beta expressed in base 2^log_beta - integer* Aiter= A; - size_t k=_ldm; - size_t k4=((k+3)>>2)+ (((k+3)%4==0)?0:1); - std::vector A0(k4<<2,0),A1(k4<<2,0),A2(k4<<2,0),A3(k4<<2,0); - integer a0,a1,a2,a3,res; - mpz_t *m0,*m1,*m2,*m3; - m0= reinterpret_cast(&a0); - m1= reinterpret_cast(&a1); - m2= reinterpret_cast(&a2); - m3= reinterpret_cast(&a3); - mp_limb_t *m0_d,*m1_d,*m2_d,*m3_d; - m0_d = m0[0]->_mp_d; - m1_d = m1[0]->_mp_d; - m2_d = m2[0]->_mp_d; - m3_d = m3[0]->_mp_d; - m0[0]->_mp_alloc = m1[0]->_mp_alloc = m2[0]->_mp_alloc = m3[0]->_mp_alloc = (int32_t)(k4*8/sizeof(mp_limb_t)); // to ensure 32 bits portability - m0[0]->_mp_size = m1[0]->_mp_size = m2[0]->_mp_size = m3[0]->_mp_size = (int32_t)(k4*8/sizeof(mp_limb_t)); // to ensure 32 bits portability - for (size_t j=0;j(&tmp); -#ifdef __FFLASFFPACK_HAVE_LITTLE_ENDIAN - A0[l ]= tptr[0]; - A1[l+1]= tptr[1]; - A2[l+2]= tptr[2]; - A3[l+3]= tptr[3]; -#else - A0[l ^ ((sizeof(mp_limb_t)/2U) - 1U)] = tptr[3]; - A1[(l+1) ^ ((sizeof(mp_limb_t)/2U) - 1U)] = tptr[2]; - A2[(l+2) ^ ((sizeof(mp_limb_t)/2U) - 1U)] = tptr[1]; - A3[(l+3) ^ ((sizeof(mp_limb_t)/2U) - 1U)] = tptr[0]; -#endif - } - // see A0,A1,A2,A3 as a the gmp integers a0,a1,a2,a3 - m0[0]->_mp_d= reinterpret_cast(&A0[0]); - m1[0]->_mp_d= reinterpret_cast(&A1[0]); - m2[0]->_mp_d= reinterpret_cast(&A2[0]); - m3[0]->_mp_d= reinterpret_cast(&A3[0]); - res = a0;res+= a1;res+= a2;res+= a3; - res%=_M; +} + } - // get the correct result according to the expected sign of A - if (res>hM) - res-=_M; - if (gamma==0) - Aiter[j+i*lda]=res; - else - if (gamma==integer(1)) - Aiter[j+i*lda]+=res; - else - if (gamma==integer(-1)) - Aiter[j+i*lda]=res-Aiter[j+i*lda]; - else{ - Aiter[j+i*lda]*=gamma; - Aiter[j+i*lda]+=res; - } - } - m0[0]->_mp_d = m0_d; - m1[0]->_mp_d = m1_d; - m2[0]->_mp_d = m2_d; - m3[0]->_mp_d = m3_d; - m0[0]->_mp_alloc = m1[0]->_mp_alloc = m2[0]->_mp_alloc= m3[0]->_mp_alloc = 1; - m0[0]->_mp_size = m1[0]->_mp_size = m2[0]->_mp_size = m3[0]->_mp_size = 0; - FFLAS::fflas_delete( A_beta); -#ifdef CHECK_RNS - bool ok=true; - for (size_t i=0;i; - using vect_t = typename simd::vect_t; - - if(_size % simd::vect_size == 0){ - for(size_t i = 0 ; i < n ; i++){ - vect_t tmp1, tmp2, tmp3, v, max, basis, inv, neg; - for(size_t j = 0 ; j < _size ; j+=simd::vect_size){ - basis = simd::load(_basis.data()+j); - inv = simd::load(_invbasis.data()+j); - max = simd::load(_basisMax.data()+j); - neg = simd::load(_negbasis.data()+j); - v = simd::load(Arns+i*_size+j); - tmp1 = simd::floor(simd::mul(v, inv)); - tmp2 = simd::fnmadd(v, tmp1, basis); - tmp1 = simd::greater(tmp2, max); - tmp3 = simd::lesser(tmp2, simd::zero()); - tmp1 = simd::vand(tmp1, neg); - tmp3 = simd::vand(tmp3, basis); - tmp1 = simd::vor(tmp1, tmp3); - tmp2 = simd::add(tmp2, tmp1); - simd::store(Arns+i*_size+j, tmp2); - } - } - }else{ - for(size_t i = 0 ; i < n ; i++){ - vect_t tmp1, tmp2, tmp3, v, max, basis, inv, neg; - size_t j = 0; - for( ; j < ROUND_DOWN(_size, simd::vect_size) ; j+=simd::vect_size){ - basis = simd::load(_basis.data()+j); - inv = simd::load(_invbasis.data()+j); - max = simd::load(_basisMax.data()+j); - neg = simd::load(_negbasis.data()+j); - v = simd::loadu(Arns+i*_size+j); - tmp1 = simd::floor(simd::mul(v, inv)); - tmp2 = simd::fnmadd(v, tmp1, basis); - tmp1 = simd::greater(tmp2, max); - tmp3 = simd::lesser(tmp2, simd::zero()); - tmp1 = simd::vand(tmp1, neg); - tmp3 = simd::vand(tmp3, basis); - tmp1 = simd::vor(tmp1, tmp3); - tmp2 = simd::add(tmp2, tmp1); - simd::storeu(Arns+i*_size+j, tmp2); - } - for( ; j < _size ; ++j){ - // std::cout << j << std::endl; - // auto x = std::floor(Arns[i*_size+j] * _invbasis[j]); - Arns[i*_size+j] -= std::floor(Arns[i*_size+j]*_invbasis[j])*_basis[j]; - // Arns[i*_size+j] = std::fma(Arns[i*_size+j], -x, _basis[j]); - if(Arns[i*_size+j] >= _basis[j]){ - Arns[i*_size+j] -= _basis[j]; - }else if(Arns[i*_size+j] < 0){ - Arns[i*_size+j] += _basis[j]; - } - } - } - } -#else - for(size_t i = 0 ; i < n ; i+= _size){ - for(size_t j = 0 ; j < _size ; ++j){ - //_field_rns.reduce(Arns+i*_size+j); - _field_rns[i].reduce(Arns[i*_size+j]); - } + using simd = Simd; + using vect_t = typename simd::vect_t; + + if(_size % simd::vect_size == 0){ + for(size_t i = 0 ; i < n ; i++){ + vect_t tmp1, tmp2, tmp3, v, max, basis, inv, neg; + for(size_t j = 0 ; j < _size ; j+=simd::vect_size){ + basis = simd::load(_basis.data()+j); + inv = simd::load(_invbasis.data()+j); + max = simd::load(_basisMax.data()+j); + neg = simd::load(_negbasis.data()+j); + v = simd::load(Arns+i*_size+j); + tmp1 = simd::floor(simd::mul(v, inv)); + tmp2 = simd::fnmadd(v, tmp1, basis); + tmp1 = simd::greater(tmp2, max); + tmp3 = simd::lesser(tmp2, simd::zero()); + tmp1 = simd::vand(tmp1, neg); + tmp3 = simd::vand(tmp3, basis); + tmp1 = simd::vor(tmp1, tmp3); + tmp2 = simd::add(tmp2, tmp1); + simd::store(Arns+i*_size+j, tmp2); + } + } + }else{ + for(size_t i = 0 ; i < n ; i++){ + vect_t tmp1, tmp2, tmp3, v, max, basis, inv, neg; + size_t j = 0; + for( ; j < ROUND_DOWN(_size, simd::vect_size) ; j+=simd::vect_size){ + basis = simd::load(_basis.data()+j); + inv = simd::load(_invbasis.data()+j); + max = simd::load(_basisMax.data()+j); + neg = simd::load(_negbasis.data()+j); + v = simd::loadu(Arns+i*_size+j); + tmp1 = simd::floor(simd::mul(v, inv)); + tmp2 = simd::fnmadd(v, tmp1, basis); + tmp1 = simd::greater(tmp2, max); + tmp3 = simd::lesser(tmp2, simd::zero()); + tmp1 = simd::vand(tmp1, neg); + tmp3 = simd::vand(tmp3, basis); + tmp1 = simd::vor(tmp1, tmp3); + tmp2 = simd::add(tmp2, tmp1); + simd::storeu(Arns+i*_size+j, tmp2); + } + for( ; j < _size ; ++j){ + // std::cout << j << std::endl; + // auto x = std::floor(Arns[i*_size+j] * _invbasis[j]); + Arns[i*_size+j] -= std::floor(Arns[i*_size+j]*_invbasis[j])*_basis[j]; + // Arns[i*_size+j] = std::fma(Arns[i*_size+j], -x, _basis[j]); + if(Arns[i*_size+j] >= _basis[j]){ + Arns[i*_size+j] -= _basis[j]; + }else if(Arns[i*_size+j] < 0){ + Arns[i*_size+j] += _basis[j]; } -#endif + } } - else { // NOT IN RNS MAJOR - // #ifndef __FFLASFFPACK_SEQUENTIAL - // auto sp=SPLITTER(MAX_THREADS); - // #else - // auto sp=SPLITTER(1); - // #endif - PARFOR1D(i,_size,SPLITTER(NUM_THREADS), - //for(size_t i=0;i<_size;i++) - FFLAS::freduce (_field_rns[i],n,Arns+i*rda,1); - ); + } +#else + for(size_t i = 0 ; i < n ; i++){ + for(size_t j = 0 ; j < _size ; ++j){ + //_field_rns.reduce(Arns+i*_size+j); + _field_rns[j].reduce(Arns[i*_size+j]); } - + } +#endif } + else { // NOT IN RNS MAJOR + // #ifndef __FFLASFFPACK_SEQUENTIAL + // auto sp=SPLITTER(MAX_THREADS); + // #else + // auto sp=SPLITTER(1); + // #endif + PARFOR1D(i,_size,SPLITTER(NUM_THREADS), + //for(size_t i=0;i<_size;i++) + FFLAS::freduce (_field_rns[i],n,Arns+i*rda,1); + ); + } + + } - } // FFPACK +} // FFPACK #endif // __FFLASFFPACK_field_rns_double_INL - /* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ - // vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s +/* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ +// vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s diff --git a/tests/test-rns-double.C b/tests/test-rns-double.C index fb80ba700..7879871b0 100644 --- a/tests/test-rns-double.C +++ b/tests/test-rns-double.C @@ -31,6 +31,7 @@ #define __FFLASFFPACK_SEQUENTIAL + #include "fflas-ffpack/fflas-ffpack-config.h" @@ -42,12 +43,13 @@ using namespace std; #include "fflas-ffpack/utils/timer.h" +#include "fflas-ffpack/utils/fflas_io.h" #include "fflas-ffpack/fflas/fflas.h" #include "fflas-ffpack/field/rns-double.h" #include "fflas-ffpack/field/rns-double-extended.h" #include "fflas-ffpack/field/rns-integer.h" #include "fflas-ffpack/utils/args-parser.h" -#include "givaro/givinteger.h" +#include "givaro/givinteger.h" #include "givaro/modular-integer.h" #include "givaro/zring.h" #include @@ -67,52 +69,142 @@ bool run_check(size_t n, size_t primes_bits, size_t b, size_t iters, size_t seed Givaro::Integer p; size_t bits=b; + bool ok=true; std::cout<<"Checking "; print_rns(); - std::cout<<" int bits="< ZZ; - if (!FFLAS::fequal(ZZ,n,n,A,n,Acopy,n)){ - std::cout<<"FAILED\n"; - return false; - } +#endif FFLAS::fflas_delete(mod_A); + FFLAS::fflas_delete(mod_At); FFLAS::fflas_delete(A); - FFLAS::fflas_delete(Acopy); + FFLAS::fflas_delete(At); + FFLAS::fflas_delete(B); } - std::cout<<"PASSED\n"; - return true; + std::cout<<"Test suite for rns-double ...."<<(ok?"PASSED\n":"FAILED\n"); + if (!ok) std::cout<<"seed was "< Date: Mon, 23 Oct 2023 16:54:25 +0200 Subject: [PATCH 8/8] test-rns-double must return false when succeed --- tests/test-fgemm.C | 2 +- tests/test-rns-double.C | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test-fgemm.C b/tests/test-fgemm.C index 9131be3da..443e51d5f 100644 --- a/tests/test-fgemm.C +++ b/tests/test-fgemm.C @@ -206,7 +206,7 @@ bool launch_MM(const Field & F, fgemm (F, ta, tb,m,n,k,alpha, A,lda, B,ldb, beta,C,ldc,WH); } ok = ok && check_MM(F, D, ta, tb,m,n,k,alpha, A,lda, B,ldb, beta,C,ldc); -prln + fflas_delete(A); fflas_delete(B); diff --git a/tests/test-rns-double.C b/tests/test-rns-double.C index 7879871b0..507345439 100644 --- a/tests/test-rns-double.C +++ b/tests/test-rns-double.C @@ -227,7 +227,7 @@ int main(int argc, char** argv){ ok = ok and run_check (m, p,b,iters,seed); ok = ok and run_check(m,2*p,b,iters,seed); - return ok; + return !ok; }