diff --git a/frame/1/bli_l1v_tapi.c b/frame/1/bli_l1v_tapi.c index b49cdc8238..1689bb10da 100644 --- a/frame/1/bli_l1v_tapi.c +++ b/frame/1/bli_l1v_tapi.c @@ -36,6 +36,133 @@ // #included from files that define the typed API macros. #ifdef BLIS_ENABLE_TAPI +// -- Optional OpenMP threading for memory-bound level-1v operations ---------- +// BLIS runs level-1/2 single-threaded, which leaves most of the memory +// bandwidth on the table for large vectors (a single GB10 X925 core sees +// ~45 GB/s vs ~118 GB/s aggregate). When BLIS_ENABLE_L1_OPENMP is defined, +// large contiguous-enough level-1v calls are split across OpenMP threads. +#ifdef BLIS_ENABLE_L1_OPENMP +#include +#ifndef BLIS_L1_MT_THRESHOLD +#define BLIS_L1_MT_THRESHOLD 200000 // min elements to bother threading +#endif +#ifndef BLIS_L1_MT_MAX +#define BLIS_L1_MT_MAX 256 // cap on partial-sum buffer / threads +#endif +#ifndef BLI_L1V_MT_HELPERS +#define BLI_L1V_MT_HELPERS +BLIS_INLINE bool bli_l1v_mt_ok( dim_t n ) +{ + // Large enough, not already nested in a parallel region, and within the + // fixed partial-sum buffer size used by the reduction path. + return ( n >= ( dim_t )BLIS_L1_MT_THRESHOLD ) && + ( omp_get_active_level() == 0 ) && + ( omp_get_max_threads() <= BLIS_L1_MT_MAX ); +} +// Contiguous [start,len) sub-range for the calling thread -- a plain EQUAL +// split. Even on heterogeneous machines (mixed fast/slow cores) equal work per +// thread is the right split for these ops: they are pure streaming with no data +// reuse, so they are memory-bandwidth bound and every core makes memory +// progress at nearly the same rate regardless of its compute width (e.g. on +// GB10 the fast-core:slow-core per-core throughput is ~1.0-1.1 for axpy). This +// is unlike compute-bound level-3 GEMM, whose data reuse would justify giving +// the faster cores proportionally more work. +BLIS_INLINE void bli_l1v_range( dim_t n, dim_t* start, dim_t* len ) +{ + const dim_t nt = omp_get_num_threads(); + const dim_t tid = omp_get_thread_num(); + const dim_t base = n / nt; + const dim_t rem = n % nt; + *start = tid * base + ( tid < rem ? tid : rem ); + *len = base + ( tid < rem ? 1 : 0 ); +} +#endif + +// axpyv/scal2v: disjoint output; split the range and call the kernel per chunk. +#define bli_l1v_axpyv_kercall( ch, ctype, f, conjx, n, alpha, x, incx, y, incy, cntx ) \ +{ \ + if ( bli_l1v_mt_ok( n ) ) \ + { \ + _Pragma( "omp parallel" ) \ + { \ + dim_t s_, l_; bli_l1v_range( (n), &s_, &l_ ); \ + if ( l_ > 0 ) \ + f( conjx, l_, ( ctype* )(alpha), ( ctype* )(x) + s_*(incx), (incx), \ + (y) + s_*(incy), (incy), ( cntx_t* )(cntx) ); \ + } \ + } \ + else \ + f( conjx, (n), ( ctype* )(alpha), ( ctype* )(x), (incx), (y), (incy), ( cntx_t* )(cntx) ); \ +} + +// dotv: reduction; each thread computes a partial dot, combined afterward. +#define bli_l1v_dotv_kercall( ch, ctype, f, conjx, conjy, n, x, incx, y, incy, rho, cntx ) \ +{ \ + if ( bli_l1v_mt_ok( n ) ) \ + { \ + ctype parts_[ BLIS_L1_MT_MAX ]; \ + dim_t nt_ = 1; \ + _Pragma( "omp parallel" ) \ + { \ + const dim_t tid_ = omp_get_thread_num(); \ + dim_t s_, l_; bli_l1v_range( (n), &s_, &l_ ); \ + ctype pr_; bli_tset0s( ch, pr_ ); \ + if ( tid_ == 0 ) nt_ = omp_get_num_threads(); \ + if ( l_ > 0 ) \ + f( conjx, conjy, l_, ( ctype* )(x) + s_*(incx), (incx), \ + ( ctype* )(y) + s_*(incy), (incy), &pr_, ( cntx_t* )(cntx) ); \ + parts_[ tid_ ] = pr_; \ + } \ + ctype acc_; bli_tset0s( ch, acc_ ); \ + for ( dim_t t_ = 0; t_ < nt_; ++t_ ) bli_tadds( ch, ch, ch, parts_[ t_ ], acc_ ); \ + *(rho) = acc_; \ + } \ + else \ + f( conjx, conjy, (n), ( ctype* )(x), (incx), ( ctype* )(y), (incy), (rho), ( cntx_t* )(cntx) ); \ +} + +// copyv/addv/subv: two-vector, disjoint output; split the range. +#define bli_l1v_copyv_kercall( ch, ctype, f, conjx, n, x, incx, y, incy, cntx ) \ +{ \ + if ( bli_l1v_mt_ok( n ) ) \ + { \ + _Pragma( "omp parallel" ) \ + { \ + dim_t s_, l_; bli_l1v_range( (n), &s_, &l_ ); \ + if ( l_ > 0 ) \ + f( conjx, l_, ( ctype* )(x) + s_*(incx), (incx), (y) + s_*(incy), (incy), ( cntx_t* )(cntx) ); \ + } \ + } \ + else \ + f( conjx, (n), ( ctype* )(x), (incx), (y), (incy), ( cntx_t* )(cntx) ); \ +} + +// scalv/invscalv/setv: single in-place vector, disjoint; split the range. +#define bli_l1v_scalv_kercall( ch, ctype, f, conjalpha, n, alpha, x, incx, cntx ) \ +{ \ + if ( bli_l1v_mt_ok( n ) ) \ + { \ + _Pragma( "omp parallel" ) \ + { \ + dim_t s_, l_; bli_l1v_range( (n), &s_, &l_ ); \ + if ( l_ > 0 ) \ + f( conjalpha, l_, ( ctype* )(alpha), (x) + s_*(incx), (incx), ( cntx_t* )(cntx) ); \ + } \ + } \ + else \ + f( conjalpha, (n), ( ctype* )(alpha), (x), (incx), ( cntx_t* )(cntx) ); \ +} +#else +#define bli_l1v_axpyv_kercall( ch, ctype, f, conjx, n, alpha, x, incx, y, incy, cntx ) \ + f( conjx, (n), ( ctype* )(alpha), ( ctype* )(x), (incx), (y), (incy), ( cntx_t* )(cntx) ) +#define bli_l1v_dotv_kercall( ch, ctype, f, conjx, conjy, n, x, incx, y, incy, rho, cntx ) \ + f( conjx, conjy, (n), ( ctype* )(x), (incx), ( ctype* )(y), (incy), (rho), ( cntx_t* )(cntx) ) +#define bli_l1v_copyv_kercall( ch, ctype, f, conjx, n, x, incx, y, incy, cntx ) \ + f( conjx, (n), ( ctype* )(x), (incx), (y), (incy), ( cntx_t* )(cntx) ) +#define bli_l1v_scalv_kercall( ch, ctype, f, conjalpha, n, alpha, x, incx, cntx ) \ + f( conjalpha, (n), ( ctype* )(alpha), (x), (incx), ( cntx_t* )(cntx) ) +#endif + // // Define BLAS-like interfaces with typed operands. // @@ -63,14 +190,7 @@ void PASTEMAC(ch,opname,EX_SUF) \ \ PASTECH(opname,_ker_ft) f = bli_cntx_get_ukr_dt( dt, kerid, cntx ); \ \ - f \ - ( \ - conjx, \ - n, \ - ( ctype* )x, incx, \ - y, incy, \ - ( cntx_t* )cntx \ - ); \ + bli_l1v_copyv_kercall( ch, ctype, f, conjx, n, x, incx, y, incy, cntx ); \ } INSERT_GENTFUNC_BASIC( addv, BLIS_ADDV_KER ) @@ -177,15 +297,7 @@ void PASTEMAC(ch,opname,EX_SUF) \ \ PASTECH(opname,_ker_ft) f = bli_cntx_get_ukr_dt( dt, kerid, cntx ); \ \ - f \ - ( \ - conjx, \ - n, \ - ( ctype* )alpha, \ - ( ctype* )x, incx, \ - y, incy, \ - ( cntx_t* )cntx \ - ); \ + bli_l1v_axpyv_kercall( ch, ctype, f, conjx, n, alpha, x, incx, y, incy, cntx ); \ } INSERT_GENTFUNC_BASIC( axpyv, BLIS_AXPYV_KER ) @@ -217,16 +329,7 @@ void PASTEMAC(ch,opname,EX_SUF) \ \ PASTECH(opname,_ker_ft) f = bli_cntx_get_ukr_dt( dt, kerid, cntx ); \ \ - f \ - ( \ - conjx, \ - conjy, \ - n, \ - ( ctype* )x, incx, \ - ( ctype* )y, incy, \ - rho, \ - ( cntx_t* )cntx \ - ); \ + bli_l1v_dotv_kercall( ch, ctype, f, conjx, conjy, n, x, incx, y, incy, rho, cntx ); \ } INSERT_GENTFUNC_BASIC( dotv, BLIS_DOTV_KER ) @@ -331,14 +434,7 @@ void PASTEMAC(ch,opname,EX_SUF) \ \ PASTECH(opname,_ker_ft) f = bli_cntx_get_ukr_dt( dt, kerid, cntx ); \ \ - f \ - ( \ - conjalpha, \ - n, \ - ( ctype* )alpha, \ - x, incx, \ - ( cntx_t* )cntx \ - ); \ + bli_l1v_scalv_kercall( ch, ctype, f, conjalpha, n, alpha, x, incx, cntx ); \ } INSERT_GENTFUNC_BASIC( invscalv, BLIS_INVSCALV_KER ) diff --git a/frame/2/gemv/bli_gemv_unf_var1.c b/frame/2/gemv/bli_gemv_unf_var1.c index 86cd41c467..416037bc2a 100644 --- a/frame/2/gemv/bli_gemv_unf_var1.c +++ b/frame/2/gemv/bli_gemv_unf_var1.c @@ -34,6 +34,62 @@ #include "blis.h" +// -- Optional OpenMP output-parallel dotxf loop for gemv (transpose) --------- +// var1 computes each output element independently (dotxf: y1 = beta*y1 + +// alpha*A1*x over f output entries), so the output loop is disjoint and can be +// split directly across threads -- no reduction. Each thread owns a contiguous +// band [i0,i1) of output indices; f is clamped to the band. Beta is applied by +// dotxf per output element. Enabled via BLIS_ENABLE_L1_OPENMP. +#ifdef BLIS_ENABLE_L1_OPENMP +#include +#ifndef BLIS_L2_MT_THRESHOLD +#define BLIS_L2_MT_THRESHOLD 262144 +#endif +#define BLI_GEMV_V1_DOTXF_LOOP( kfp_df, conja, conjx, n_elem, n_iter, b_fuse, \ + alpha, a, rs_at, cs_at, x, incx, beta, y, incy, cntx, mn ) \ +{ \ + if ( ( uint64_t )(mn) >= ( uint64_t )BLIS_L2_MT_THRESHOLD && \ + omp_get_active_level() == 0 && omp_get_max_threads() > 1 ) \ + { \ + _Pragma( "omp parallel" ) \ + { \ + const dim_t nt_ = omp_get_num_threads(), tid_ = omp_get_thread_num(); \ + const dim_t bs_ = (n_iter) / nt_, rm_ = (n_iter) % nt_; \ + const dim_t i0_ = tid_ * bs_ + ( tid_ < rm_ ? tid_ : rm_ ); \ + const dim_t i1_ = i0_ + bs_ + ( tid_ < rm_ ? 1 : 0 ); \ + dim_t i_, f_; \ + for ( i_ = i0_; i_ < i1_; i_ += f_ ) { \ + f_ = bli_determine_blocksize_dim_f( i_, i1_, (b_fuse) ); \ + kfp_df( (conja), (conjx), (n_elem), f_, (alpha), \ + (a) + i_*(rs_at), (cs_at), (rs_at), (x), (incx), \ + (beta), (y) + i_*(incy), (incy), (cntx) ); \ + } \ + } \ + } \ + else { \ + dim_t i_, f_; \ + for ( i_ = 0; i_ < (n_iter); i_ += f_ ) { \ + f_ = bli_determine_blocksize_dim_f( i_, (n_iter), (b_fuse) ); \ + kfp_df( (conja), (conjx), (n_elem), f_, (alpha), \ + (a) + i_*(rs_at), (cs_at), (rs_at), (x), (incx), \ + (beta), (y) + i_*(incy), (incy), (cntx) ); \ + } \ + } \ +} +#else +#define BLI_GEMV_V1_DOTXF_LOOP( kfp_df, conja, conjx, n_elem, n_iter, b_fuse, \ + alpha, a, rs_at, cs_at, x, incx, beta, y, incy, cntx, mn ) \ +{ \ + dim_t i_, f_; \ + for ( i_ = 0; i_ < (n_iter); i_ += f_ ) { \ + f_ = bli_determine_blocksize_dim_f( i_, (n_iter), (b_fuse) ); \ + kfp_df( (conja), (conjx), (n_elem), f_, (alpha), \ + (a) + i_*(rs_at), (cs_at), (rs_at), (x), (incx), \ + (beta), (y) + i_*(incy), (incy), (cntx) ); \ + } \ +} +#endif + #undef GENTFUNC #define GENTFUNC( ctype, ch, varname ) \ \ @@ -72,30 +128,11 @@ void PASTEMAC(ch,varname) \ dotxf_ker_ft kfp_df = bli_cntx_get_ukr_dt( dt, BLIS_DOTXF_KER, cntx ); \ b_fuse = bli_cntx_get_blksz_def_dt( dt, BLIS_DF, cntx ); \ \ - for ( i = 0; i < n_iter; i += f ) \ - { \ - f = bli_determine_blocksize_dim_f( i, n_iter, b_fuse ); \ -\ - A1 = a + (i )*rs_at + (0 )*cs_at; \ - x1 = x + (0 )*incy; \ - y1 = y + (i )*incy; \ -\ - /* y1 = beta * y1 + alpha * A1 * x; */ \ - kfp_df \ - ( \ - conja, \ - conjx, \ - n_elem, \ - f, \ - alpha, \ - A1, cs_at, rs_at, \ - x1, incx, \ - beta, \ - y1, incy, \ - cntx \ - ); \ -\ - } \ + /* y = beta*y + alpha * op(A) * x, output-parallel when enabled. */ \ + ( void )A1; ( void )x1; ( void )y1; ( void )i; ( void )f; \ + BLI_GEMV_V1_DOTXF_LOOP( kfp_df, conja, conjx, n_elem, n_iter, b_fuse, \ + alpha, a, rs_at, cs_at, x, incx, beta, y, incy, cntx, \ + ( uint64_t )m * ( uint64_t )n ); \ } INSERT_GENTFUNC_BASIC( gemv_unf_var1 ) diff --git a/frame/2/gemv/bli_gemv_unf_var2.c b/frame/2/gemv/bli_gemv_unf_var2.c index f1be39ca57..95a3439b98 100644 --- a/frame/2/gemv/bli_gemv_unf_var2.c +++ b/frame/2/gemv/bli_gemv_unf_var2.c @@ -34,6 +34,64 @@ #include "blis.h" +// -- Optional OpenMP row-parallel axpyf loop for gemv (no-transpose) --------- +// gemv accumulates A's columns into y, so we cannot split the column loop +// without a reduction. Instead we split y's rows (n_elem): each thread runs +// the full column loop over its own contiguous band of rows, with offset A/y +// pointers -- disjoint outputs, no reduction. Beta-scaling of y already +// happened (serially) before this loop. Enabled via BLIS_ENABLE_L1_OPENMP. +#ifdef BLIS_ENABLE_L1_OPENMP +#include +#ifndef BLIS_L2_MT_THRESHOLD +#define BLIS_L2_MT_THRESHOLD 262144 // min m*n to thread +#endif +#define BLI_GEMV_V2_AXPYF_LOOP( kfp_af, conja, conjx, n_elem, n_iter, b_fuse, \ + alpha, a, rs_at, cs_at, x, incx, y, incy, cntx, mn ) \ +{ \ + if ( ( uint64_t )(mn) >= ( uint64_t )BLIS_L2_MT_THRESHOLD && \ + omp_get_active_level() == 0 && omp_get_max_threads() > 1 ) \ + { \ + _Pragma( "omp parallel" ) \ + { \ + const dim_t nt_ = omp_get_num_threads(), tid_ = omp_get_thread_num(); \ + const dim_t bs_ = (n_elem) / nt_, rm_ = (n_elem) % nt_; \ + const dim_t r0_ = tid_ * bs_ + ( tid_ < rm_ ? tid_ : rm_ ); \ + const dim_t rl_ = bs_ + ( tid_ < rm_ ? 1 : 0 ); \ + if ( rl_ > 0 ) { \ + dim_t i_, f_; \ + for ( i_ = 0; i_ < (n_iter); i_ += f_ ) { \ + f_ = bli_determine_blocksize_dim_f( i_, (n_iter), (b_fuse) ); \ + kfp_af( (conja), (conjx), rl_, f_, (alpha), \ + (a) + r0_*(rs_at) + i_*(cs_at), (rs_at), (cs_at), \ + (x) + i_*(incx), (incx), (y) + r0_*(incy), (incy), (cntx) ); \ + } \ + } \ + } \ + } \ + else { \ + dim_t i_, f_; \ + for ( i_ = 0; i_ < (n_iter); i_ += f_ ) { \ + f_ = bli_determine_blocksize_dim_f( i_, (n_iter), (b_fuse) ); \ + kfp_af( (conja), (conjx), (n_elem), f_, (alpha), \ + (a) + i_*(cs_at), (rs_at), (cs_at), \ + (x) + i_*(incx), (incx), (y), (incy), (cntx) ); \ + } \ + } \ +} +#else +#define BLI_GEMV_V2_AXPYF_LOOP( kfp_af, conja, conjx, n_elem, n_iter, b_fuse, \ + alpha, a, rs_at, cs_at, x, incx, y, incy, cntx, mn ) \ +{ \ + dim_t i_, f_; \ + for ( i_ = 0; i_ < (n_iter); i_ += f_ ) { \ + f_ = bli_determine_blocksize_dim_f( i_, (n_iter), (b_fuse) ); \ + kfp_af( (conja), (conjx), (n_elem), f_, (alpha), \ + (a) + i_*(cs_at), (rs_at), (cs_at), \ + (x) + i_*(incx), (incx), (y), (incy), (cntx) ); \ + } \ +} +#endif + #undef GENTFUNC #define GENTFUNC( ctype, ch, varname ) \ \ @@ -101,28 +159,11 @@ void PASTEMAC(ch,varname) \ axpyf_ker_ft kfp_af = bli_cntx_get_ukr_dt( dt, BLIS_AXPYF_KER, cntx ); \ b_fuse = bli_cntx_get_blksz_def_dt( dt, BLIS_AF, cntx ); \ \ - for ( i = 0; i < n_iter; i += f ) \ - { \ - f = bli_determine_blocksize_dim_f( i, n_iter, b_fuse ); \ -\ - A1 = a + (0 )*rs_at + (i )*cs_at; \ - x1 = x + (i )*incx; \ - y1 = y + (0 )*incy; \ -\ - /* y = y + alpha * A1 * x1; */ \ - kfp_af \ - ( \ - conja, \ - conjx, \ - n_elem, \ - f, \ - alpha, \ - A1, rs_at, cs_at, \ - x1, incx, \ - y1, incy, \ - cntx \ - ); \ - } \ + /* y = y + alpha * A * x, row-parallel when enabled (see macro above). */ \ + ( void )A1; ( void )x1; ( void )y1; ( void )i; ( void )f; \ + BLI_GEMV_V2_AXPYF_LOOP( kfp_af, conja, conjx, n_elem, n_iter, b_fuse, \ + alpha, a, rs_at, cs_at, x, incx, y, incy, cntx, \ + ( uint64_t )m * ( uint64_t )n ); \ } INSERT_GENTFUNC_BASIC( gemv_unf_var2 ) diff --git a/frame/2/ger/bli_ger_unb_var1.c b/frame/2/ger/bli_ger_unb_var1.c index 0b07611302..8edacf02d0 100644 --- a/frame/2/ger/bli_ger_unb_var1.c +++ b/frame/2/ger/bli_ger_unb_var1.c @@ -34,6 +34,55 @@ #include "blis.h" +// Optional OpenMP row-parallel rank-1 update (ger, var1). Each row of A is an +// independent axpyv (A[i,:] += alpha*conj(x[i]) * y), so split the rows across +// threads -- disjoint output, no reduction, raw axpyv kernel (no nesting). +// Enabled via BLIS_ENABLE_L1_OPENMP. +#ifdef BLIS_ENABLE_L1_OPENMP +#include +#ifndef BLIS_L2_MT_THRESHOLD +#define BLIS_L2_MT_THRESHOLD 262144 +#endif +#define BLI_GER_V1_ROWS( ch, ctype, kfp_av, conjx, conjy, m, n, alpha, x, incx, y, incy, a, rs_a, cs_a, cntx ) \ +{ \ + if ( ( uint64_t )(m)*( uint64_t )(n) >= ( uint64_t )BLIS_L2_MT_THRESHOLD && \ + omp_get_active_level() == 0 && omp_get_max_threads() > 1 ) \ + { \ + _Pragma( "omp parallel" ) \ + { \ + const dim_t nt_ = omp_get_num_threads(), tid_ = omp_get_thread_num(); \ + const dim_t bs_ = (m) / nt_, rm_ = (m) % nt_; \ + const dim_t i0_ = tid_*bs_ + ( tid_ < rm_ ? tid_ : rm_ ); \ + const dim_t i1_ = i0_ + bs_ + ( tid_ < rm_ ? 1 : 0 ); \ + for ( dim_t i_ = i0_; i_ < i1_; ++i_ ) { \ + ctype ac_; \ + bli_tcopycjs( ch,ch, (conjx), *((x) + i_*(incx)), ac_ ); \ + bli_tscals( ch,ch,ch, *(alpha), ac_ ); \ + kfp_av( (conjy), (n), &ac_, (y), (incy), (a) + i_*(rs_a), (cs_a), (cntx) ); \ + } \ + } \ + } \ + else { \ + for ( dim_t i_ = 0; i_ < (m); ++i_ ) { \ + ctype ac_; \ + bli_tcopycjs( ch,ch, (conjx), *((x) + i_*(incx)), ac_ ); \ + bli_tscals( ch,ch,ch, *(alpha), ac_ ); \ + kfp_av( (conjy), (n), &ac_, (y), (incy), (a) + i_*(rs_a), (cs_a), (cntx) ); \ + } \ + } \ +} +#else +#define BLI_GER_V1_ROWS( ch, ctype, kfp_av, conjx, conjy, m, n, alpha, x, incx, y, incy, a, rs_a, cs_a, cntx ) \ +{ \ + for ( dim_t i_ = 0; i_ < (m); ++i_ ) { \ + ctype ac_; \ + bli_tcopycjs( ch,ch, (conjx), *((x) + i_*(incx)), ac_ ); \ + bli_tscals( ch,ch,ch, *(alpha), ac_ ); \ + kfp_av( (conjy), (n), &ac_, (y), (incy), (a) + i_*(rs_a), (cs_a), (cntx) ); \ + } \ +} +#endif + #undef GENTFUNC #define GENTFUNC( ctype, ch, varname ) \ \ @@ -61,26 +110,8 @@ void PASTEMAC(ch,varname) \ /* Query the context for the kernel function pointer. */ \ axpyv_ker_ft kfp_av = bli_cntx_get_ukr_dt( dt, BLIS_AXPYV_KER, cntx ); \ \ - for ( i = 0; i < m; ++i ) \ - { \ - a1t = a + (i )*rs_a + (0 )*cs_a; \ - chi1 = x + (i )*incx; \ - y1 = y + (0 )*incy; \ -\ - /* a1t = a1t + alpha * chi1 * y; */ \ - bli_tcopycjs( ch,ch, conjx, *chi1, alpha_chi1 ); \ - bli_tscals( ch,ch,ch, *alpha, alpha_chi1 ); \ -\ - kfp_av \ - ( \ - conjy, \ - n, \ - &alpha_chi1, \ - y1, incy, \ - a1t, cs_a, \ - cntx \ - ); \ - } \ + ( void )a1t; ( void )chi1; ( void )y1; ( void )alpha_chi1; ( void )i; \ + BLI_GER_V1_ROWS( ch, ctype, kfp_av, conjx, conjy, m, n, alpha, x, incx, y, incy, a, rs_a, cs_a, cntx ); \ } INSERT_GENTFUNC_BASIC( ger_unb_var1 ) diff --git a/frame/2/ger/bli_ger_unb_var2.c b/frame/2/ger/bli_ger_unb_var2.c index ee40223ca4..25897087bd 100644 --- a/frame/2/ger/bli_ger_unb_var2.c +++ b/frame/2/ger/bli_ger_unb_var2.c @@ -34,6 +34,56 @@ #include "blis.h" +// Optional OpenMP column-parallel rank-1 update (ger). Each column of A is an +// independent axpyv (A[:,j] += alpha*conj(y[j]) * x), so we split the columns +// across threads -- disjoint output, no reduction. The per-column kernel is the +// raw axpyv micro-kernel (not the threaded tapi wrapper), so there is no +// nesting. Enabled via BLIS_ENABLE_L1_OPENMP. +#ifdef BLIS_ENABLE_L1_OPENMP +#include +#ifndef BLIS_L2_MT_THRESHOLD +#define BLIS_L2_MT_THRESHOLD 262144 +#endif +#define BLI_GER_V2_COLS( ch, ctype, kfp_av, conjx, conjy, m, n, alpha, x, incx, y, incy, a, rs_a, cs_a, cntx ) \ +{ \ + if ( ( uint64_t )(m)*( uint64_t )(n) >= ( uint64_t )BLIS_L2_MT_THRESHOLD && \ + omp_get_active_level() == 0 && omp_get_max_threads() > 1 ) \ + { \ + _Pragma( "omp parallel" ) \ + { \ + const dim_t nt_ = omp_get_num_threads(), tid_ = omp_get_thread_num(); \ + const dim_t bs_ = (n) / nt_, rm_ = (n) % nt_; \ + const dim_t j0_ = tid_*bs_ + ( tid_ < rm_ ? tid_ : rm_ ); \ + const dim_t j1_ = j0_ + bs_ + ( tid_ < rm_ ? 1 : 0 ); \ + for ( dim_t j_ = j0_; j_ < j1_; ++j_ ) { \ + ctype ap_; \ + bli_tcopycjs( ch,ch, (conjy), *((y) + j_*(incy)), ap_ ); \ + bli_tscals( ch,ch,ch, *(alpha), ap_ ); \ + kfp_av( (conjx), (m), &ap_, (x), (incx), (a) + j_*(cs_a), (rs_a), (cntx) ); \ + } \ + } \ + } \ + else { \ + for ( dim_t j_ = 0; j_ < (n); ++j_ ) { \ + ctype ap_; \ + bli_tcopycjs( ch,ch, (conjy), *((y) + j_*(incy)), ap_ ); \ + bli_tscals( ch,ch,ch, *(alpha), ap_ ); \ + kfp_av( (conjx), (m), &ap_, (x), (incx), (a) + j_*(cs_a), (rs_a), (cntx) ); \ + } \ + } \ +} +#else +#define BLI_GER_V2_COLS( ch, ctype, kfp_av, conjx, conjy, m, n, alpha, x, incx, y, incy, a, rs_a, cs_a, cntx ) \ +{ \ + for ( dim_t j_ = 0; j_ < (n); ++j_ ) { \ + ctype ap_; \ + bli_tcopycjs( ch,ch, (conjy), *((y) + j_*(incy)), ap_ ); \ + bli_tscals( ch,ch,ch, *(alpha), ap_ ); \ + kfp_av( (conjx), (m), &ap_, (x), (incx), (a) + j_*(cs_a), (rs_a), (cntx) ); \ + } \ +} +#endif + #undef GENTFUNC #define GENTFUNC( ctype, ch, varname ) \ \ @@ -61,26 +111,8 @@ void PASTEMAC(ch,varname) \ /* Query the context for the kernel function pointer. */ \ axpyv_ker_ft kfp_av = bli_cntx_get_ukr_dt( dt, BLIS_AXPYV_KER, cntx ); \ \ - for ( j = 0; j < n; ++j ) \ - { \ - a1 = a + (0 )*rs_a + (j )*cs_a; \ - x1 = x + (0 )*incx; \ - psi1 = y + (j )*incy; \ -\ - /* a1 = a1 + alpha * psi1 * x; */ \ - bli_tcopycjs( ch,ch, conjy, *psi1, alpha_psi1 ); \ - bli_tscals( ch,ch,ch, *alpha, alpha_psi1 ); \ -\ - kfp_av \ - ( \ - conjx, \ - m, \ - &alpha_psi1, \ - x1, incx, \ - a1, rs_a, \ - cntx \ - ); \ - } \ + ( void )a1; ( void )x1; ( void )psi1; ( void )alpha_psi1; ( void )j; \ + BLI_GER_V2_COLS( ch, ctype, kfp_av, conjx, conjy, m, n, alpha, x, incx, y, incy, a, rs_a, cs_a, cntx ); \ } INSERT_GENTFUNC_BASIC( ger_unb_var2 ) diff --git a/frame/util/bli_util_unb_var1.c b/frame/util/bli_util_unb_var1.c index b0d33e67fc..1476ef5e14 100644 --- a/frame/util/bli_util_unb_var1.c +++ b/frame/util/bli_util_unb_var1.c @@ -40,6 +40,111 @@ // Define BLAS-like interfaces with typed operands. // +// Optional OpenMP threading for the asumv reduction (memory-bound, and BLIS +// otherwise runs it single-threaded). Enabled via BLIS_ENABLE_L1_OPENMP. +#ifdef BLIS_ENABLE_L1_OPENMP +#include +#ifndef BLIS_L1_MT_THRESHOLD +#define BLIS_L1_MT_THRESHOLD 200000 +#endif +#ifndef BLIS_L1_MT_MAX +#define BLIS_L1_MT_MAX 256 +#endif +#ifndef BLI_UTIL_MT_HELPERS +#define BLI_UTIL_MT_HELPERS +BLIS_INLINE bool bli_util_mt_ok( dim_t n ) +{ + return ( n >= ( dim_t )BLIS_L1_MT_THRESHOLD ) && + ( omp_get_active_level() == 0 ) && + ( omp_get_max_threads() <= BLIS_L1_MT_MAX ); +} +BLIS_INLINE void bli_util_range( dim_t n, dim_t* s, dim_t* l ) +{ + const dim_t nt = omp_get_num_threads(), tid = omp_get_thread_num(); + const dim_t b = n / nt, r = n % nt; + *s = tid * b + ( tid < r ? tid : r ); + *l = b + ( tid < r ? 1 : 0 ); +} +#endif +#endif + +// Serial accumulation of sum(|re|+|im|) over [s,e) into acc. +#define BLIS_ASUMV_RANGE( ch, chr, ctype, ctype_r, x, incx, s, e, acc ) \ +{ \ + for ( dim_t i_ = (s); i_ < (e); ++i_ ) { \ + ctype* c_ = (x) + i_*(incx); \ + ctype_r r_, im_; \ + bli_tgets( ch,chr, *c_, r_, im_ ); \ + r_ = bli_fabs( r_ ); im_ = bli_fabs( im_ ); \ + bli_tadds( chr,chr,chr, r_, acc ); \ + bli_tadds( chr,chr,chr, im_, acc ); \ + } \ +} + +#ifdef BLIS_ENABLE_L1_OPENMP +#define BLIS_ASUMV_ACCUM( ch, chr, ctype, ctype_r, n, x, incx, absum ) \ +{ \ + if ( bli_util_mt_ok( n ) ) \ + { \ + ctype_r parts_[ BLIS_L1_MT_MAX ]; dim_t nt_ = 1; \ + _Pragma( "omp parallel" ) \ + { \ + const dim_t tid_ = omp_get_thread_num(); \ + dim_t s_, l_; bli_util_range( (n), &s_, &l_ ); \ + if ( tid_ == 0 ) nt_ = omp_get_num_threads(); \ + ctype_r pa_; bli_tset0s( chr, pa_ ); \ + BLIS_ASUMV_RANGE( ch, chr, ctype, ctype_r, x, incx, s_, s_+l_, pa_ ); \ + parts_[ tid_ ] = pa_; \ + } \ + for ( dim_t t_ = 0; t_ < nt_; ++t_ ) bli_tadds( chr,chr,chr, parts_[ t_ ], absum ); \ + } \ + else \ + BLIS_ASUMV_RANGE( ch, chr, ctype, ctype_r, x, incx, 0, (n), absum ); \ +} +#else +#define BLIS_ASUMV_ACCUM( ch, chr, ctype, ctype_r, n, x, incx, absum ) \ + BLIS_ASUMV_RANGE( ch, chr, ctype, ctype_r, x, incx, 0, (n), absum ) +#endif + +// Threaded sum-of-squares for normfv: run the existing overflow-safe sumsqv +// kernel on each thread's contiguous chunk (each yields a (scale,sumsq) pair +// with chunk-sum-of-squares = scale^2 * sumsq), then merge the pairs with the +// standard LAPACK dlassq combine (rescale to the larger scale -- overflow-safe). +// (scale,sumsq) must be pre-initialized to (0,1) by the caller, matching the +// serial kernel's contract. Reuses bli_util_mt_ok/bli_util_range (above). +#ifdef BLIS_ENABLE_L1_OPENMP +#define BLIS_SUMSQV_MT( ch, chr, ctype_r, kername, n, x, incx, scale, sumsq, cntx, rntm ) \ +{ \ + if ( bli_util_mt_ok( n ) ) \ + { \ + ctype_r sc_[ BLIS_L1_MT_MAX ], sq_[ BLIS_L1_MT_MAX ]; dim_t nt_ = 1; \ + _Pragma( "omp parallel" ) \ + { \ + const dim_t tid_ = omp_get_thread_num(); \ + dim_t s_, l_; bli_util_range( (n), &s_, &l_ ); \ + if ( tid_ == 0 ) nt_ = omp_get_num_threads(); \ + ctype_r lsc_ = ( ctype_r )0, lsq_ = ( ctype_r )1; \ + PASTEMAC(ch,kername)( l_, (x) + s_*(incx), (incx), &lsc_, &lsq_, (cntx), (rntm) ); \ + sc_[ tid_ ] = lsc_; sq_[ tid_ ] = lsq_; \ + } \ + ctype_r cs_ = ( ctype_r )0, cq_ = ( ctype_r )1; \ + for ( dim_t t_ = 0; t_ < nt_; ++t_ ) { \ + const ctype_r as_ = sc_[ t_ ], aq_ = sq_[ t_ ]; \ + if ( as_ != ( ctype_r )0 ) { \ + if ( cs_ >= as_ ) { const ctype_r r_ = as_ / cs_; cq_ = cq_ + aq_ * r_ * r_; } \ + else { const ctype_r r_ = cs_ / as_; cq_ = aq_ + cq_ * r_ * r_; cs_ = as_; } \ + } \ + } \ + (scale) = cs_; (sumsq) = cq_; \ + } \ + else \ + PASTEMAC(ch,kername)( (n), (x), (incx), &(scale), &(sumsq), (cntx), (rntm) ); \ +} +#else +#define BLIS_SUMSQV_MT( ch, chr, ctype_r, kername, n, x, incx, scale, sumsq, cntx, rntm ) \ + PASTEMAC(ch,kername)( (n), (x), (incx), &(scale), &(sumsq), (cntx), (rntm) ) +#endif + #undef GENTFUNCR #define GENTFUNCR( ctype, ctype_r, ch, chr, varname ) \ \ @@ -52,30 +157,10 @@ void PASTEMAC(ch,varname) \ rntm_t* rntm \ ) \ { \ - ctype* chi1; \ - ctype_r chi1_r; \ - ctype_r chi1_i; \ ctype_r absum; \ - dim_t i; \ -\ - /* Initialize the absolute sum accumulator to zero. */ \ bli_tset0s( chr, absum ); \ \ - for ( i = 0; i < n; ++i ) \ - { \ - chi1 = x + (i )*incx; \ -\ - /* Get the real and imaginary components of chi1. */ \ - bli_tgets( ch,chr, *chi1, chi1_r, chi1_i ); \ -\ - /* Replace chi1_r and chi1_i with their absolute values. */ \ - chi1_r = bli_fabs( chi1_r ); \ - chi1_i = bli_fabs( chi1_i ); \ -\ - /* Accumulate the real and imaginary components into absum. */ \ - bli_tadds( chr,chr,chr, chi1_r, absum ); \ - bli_tadds( chr,chr,chr, chi1_i, absum ); \ - } \ + BLIS_ASUMV_ACCUM( ch, chr, ctype, ctype_r, n, x, incx, absum ); \ \ /* Store the final value of absum to the output variable. */ \ bli_tcopys( chr,chr, absum, *asum ); \ @@ -325,15 +410,7 @@ void PASTEMAC(ch,varname) \ bli_tcopys( chr,chr, *one, sumsq ); \ \ /* Compute the sum of the squares of the vector. */ \ - PASTEMAC(ch,kername) \ - ( \ - n, \ - x, incx, \ - &scale, \ - &sumsq, \ - cntx, \ - rntm \ - ); \ + BLIS_SUMSQV_MT( ch, chr, ctype_r, kername, n, x, incx, scale, sumsq, cntx, rntm ); \ \ /* Compute: norm = scale * sqrt( sumsq ) */ \ bli_tsqrt2s( chr,chr,chr, sumsq, sqrt_sumsq ); \ @@ -417,15 +494,7 @@ void PASTEMAC(ch,varname) \ } \ \ /* Compute the sum of the squares of the vector. */ \ - PASTEMAC(ch,kername) \ - ( \ - n, \ - x, incx, \ - &scale, \ - &sumsq, \ - cntx, \ - rntm \ - ); \ + BLIS_SUMSQV_MT( ch, chr, ctype_r, kername, n, x, incx, scale, sumsq, cntx, rntm ); \ \ /* Compute: norm = scale * sqrt( sumsq ) */ \ tsqrt2s( chr, sumsq, sqrt_sumsq ); \ @@ -461,15 +530,7 @@ void PASTEMAC(ch,varname) \ \ /* Compute the sum of the squares of the vector. */ \ \ - PASTEMAC(ch,kername) \ - ( \ - n, \ - x, incx, \ - &scale, \ - &sumsq, \ - cntx, \ - rntm \ - ); \ + BLIS_SUMSQV_MT( ch, chr, ctype_r, kername, n, x, incx, scale, sumsq, cntx, rntm ); \ \ /* Compute: norm = scale * sqrt( sumsq ) */ \ bli_tsqrt2s( chr,chr,chr, sumsq, sqrt_sumsq ); \