This wiki describes the computational kernels used by the BLIS framework.
One of the primary features of BLIS is that it provides a large set of dense linear algebra functionality while simultaneously minimizing the amount of kernel code that must be optimized for a given architecture. BLIS does this by isolating a handful of kernels which, when implemented, facilitate functionality and performance of several of the higher-level operations.
Presently, BLIS supports several groups of operations:
- Level-1v: Operations on vectors:
- Level-1d: Element-wise operations on matrix diagonals:
- Level-1m: Element-wise operations on matrices:
- Level-1f: Fused operations on multiple vectors:
- Level-2: Operations with one matrix and (at least) one vector operand:
- Level-3: Operations with matrices that are multiplication-like:
- Utility: Miscellaneous operations on matrices and vectors:
Most of the interest with BLAS libraries centers around level-3 operations because they exhibit favorable ratios of floating-point operations (flops) to memory operations (memops), which allows high performance. Some applications also require level-2 computation; however, these operations are at an inherent disadvantage on modern architectures due to their less favorable flop-to-memop ratio. The BLIS framework allows developers to quickly and easily build high performance level-3 operations, as well as relatively well-performing level-2 operations, simply by optimizing a small set of kernels. These kernels, and their relationship to the other higher-level operations supported by BLIS, are the subject of this wiki.
Some level-1v, level-1m, and level-1d operations may also be accelerated, but since they are memory-bound, optimization typically yields minor performance improvement.
This section lists and briefly describes each of the main computational kernels supported by the BLIS framework. (Other kernels are supported, but they are not of interest to most developers.)
BLIS supports the following three level-3 microkernels. These microkernels are used to implement optimized level-3 operations.
- gemm: The
gemm
microkernel performs a small matrix multiplication and is used by every level-3 operation. - trsm: The
trsm
microkernel performs a small triangular solve with multiple right-hand sides. It is not required for optimal performance and in fact is only needed when the developer opts to not implement the fusedgemmtrsm
kernel. - gemmtrsm: The
gemmtrsm
microkernel implements a fused operation whereby agemm
and atrsm
subproblem are fused together in a single routine. This avoids redundant memory operations that would otherwise be incurred if the operations were executed separately.
The following shows the steps one would take to optimize, to varying degrees, the level-3 operations supported by BLIS:
- By implementing and optimizing the
gemm
microkernel, all level-3 operations excepttrsm
are fully optimized. In this scenario, thetrsm
operation may achieve 60-90% of attainable peak performance, depending on the architecture and problem size. - If one goes further and implements and optimizes the
trsm
microkernel, this kernel, when paired with an optimizedgemm
microkernel, results in atrsm
implementation that is accelerated (but not optimized). - Alternatively, if one implements and optimizes the fused
gemmtrsm
microkernel, this kernel, when paired with an optimizedgemm
microkernel, enables a fully optimizedtrsm
implementation.
BLIS supports the following five level-1f (fused) kernels. These kernels are used to implement optimized level-2 operations (as well as self-similar level-1f operations; that is, the axpyf
kernel can be invoked indirectly via the axpyf
operation).
- axpy2v: Performs and fuses two axpyv operations, accumulating to the same output vector.
- dotaxpyv: Performs and fuses a dotv followed by an axpyv operation with x.
- axpyf: Performs and fuses some implementation-dependent number of axpyv operations, accumulating to the same output vector. Can also be expressed as a gemv operation where matrix A is m x nf, where nf is the number of fused operations (fusing factor).
- dotxf: Performs and fuses some implementation-dependent number of dotxv operations, reusing the
y
vector for each dotxv. - dotxaxpyf: Performs and fuses a dotxf and axpyf in which the matrix operand is reused.
BLIS supports the following 14 level-1v kernels. These kernels are used primarily to implement their self-similar operations. However, they are occasionally used to handle special cases of level-1f kernels or in situations where level-2 operations are partially optimized.
- addv: Performs a vector addition operation.
- amaxv: Performs a search for the index of the element with the largest absolute value (or complex modulus).
- axpyv: Performs a vector scale-and-accumulate operation.
- axpbyv: Performs an extended vector scale-and-accumulate operation similar to axpyv except that the output vector is scaled by a second scalar.
- copyv: Performs a vector copy operation
- dotv: Performs a dot product where the output scalar is overwritten.
- dotxv: Performs an extended dot product operation where the dot product is first scaled and then accumulated into a scaled output scalar.
- invertv: Performs an element-wise vector inversion operation.
- scalv: Performs an in-place (destructive) vector scaling operation.
- scal2v: Performs an out-of-place (non-destructive) vector scaling operation.
- setv: Performs a vector broadcast operation.
- subv: Performs a vector subtraction operation.
- swapv: Performs a vector swap operation.
- xpbyv: Performs a alternate vector scale-and-accumulate operation.
The table below shows dependencies between level-2 operations and each of the level-1v and level-1f kernels.
Kernels marked with a "1" for a given level-2 operation are preferred for optimization because they facilitate an optimal implementation on most architectures. Kernels marked with a "2", "3", or "4" denote those which need to be optimized for alternative implementations that would typically be second, third, or fourth choices, respectively, if the preferred kernels are not optimized.
operation / kernel | effective storage | axpyv |
dotxv |
axpy2v |
dotaxpyv |
axpyf |
dotxf |
dotxaxpyf |
---|---|---|---|---|---|---|---|---|
gemv, trmv, trsv |
row-wise | 2 | 1 | |||||
column-wise | 2 | 1 | ||||||
hemv, symv |
row- or column-wise | 4 | 4 | 3 | 2 | 2 | 1 | |
ger, her, syr |
row- or column-wise | 1 | ||||||
her2, syr2 |
row- or column-wise | 2 | 1 |
Note: The "effective storage" column reflects the orientation of the matrix operand after transposition via the corresponding trans_t
parameter (if applicable). For example, calling gemv
with a column-stored matrix A
and the transa
parameter equal to BLIS_TRANSPOSE
would be effectively equivalent to row-wise storage.
Note that all kernels, whether they be reference implementations or based on fully optimized assembly code, use names that are architecture- and implementation-specific. (This appears as a <suffix>
in the kernel reference below.) Therefore, the easiest way to call the kernel is by querying a pointer from a valid context.
The first step is to obtain a valid context. Contexts store all of the information specific to a particular sub-configuration (usually loosely specific to a microarchitecture or group of closely-related microarchitectuers). If a context is not already available in your current scope, a default context for the hardware for which BLIS was configured (or, in the case of multi-configuration builds, the hardware on which BLIS is currently running) may be queried via:
cntx_t* bli_gks_query_cntx( void );
Once this cntx_t*
pointer is obtained, you may call one of three functions to query any of the computation kernels described in this document:
void* bli_cntx_get_l3_nat_ukr_dt
(
num_t dt,
l3ukr_t ker_id,
cntx_t* cntx
);
void* bli_cntx_get_l1f_ker_dt
(
num_t dt,
l1fkr_t ker_id,
cntx_t* cntx
);
void* bli_cntx_get_l1v_ker_dt
(
num_t dt,
l1vkr_t ker_id,
cntx_t* cntx
);
The dt
and ker_id
parameters specify the floating-point datatype and the
kernel operation you wish to query, respectively.
Valid values for dt
are BLIS_FLOAT
, BLIS_DOUBLE
, BLIS_SCOMPLEX
, and
BLIS_DCOMPLEX
for single- and double-precision real, and single- and
double-precision complex, respectively.
Valid values for ker_id
are given in the tables below.
Also, note that the return values of bli_cntx_get_l1v_ker_dt
bli_cntx_get_l1f_ker_dt()
, and bli_cntx_get_l3_nat_ukr_dt()
,
will be void*
and must be typecast to typed function pointers before being called.
As a convenience, BLIS defines function pointer types appropriate for usage in these
situations. The function pointer type for each operation is given in the third
columns of each table, with the ?
taking the place of one of the supported
datatype characters.
kernel operation | l3ukr_t | function pointer type |
---|---|---|
gemm | BLIS_GEMM |
?gemm_ukr_ft |
trsm_l | BLIS_TRSM_L_UKR |
?trsm_ukr_ft |
trsm_u | BLIS_TRSM_U_UKR |
?trsm_ukr_ft |
gemmtrsm_l | BLIS_GEMMTRSM_L_UKR |
?gemmtrsm_ukr_ft |
gemmtrsm_u | BLIS_GEMMTRSM_U_UKR |
?gemmtrsm_ukr_ft |
kernel operation | l1fkr_t | function pointer type |
---|---|---|
axpy2v | BLIS_AXPY2V_KER |
?axpy2v_ft |
dotaxpyv | BLIS_DOTAXPYV_KER |
?dotaxpyv_ft |
axpyf | BLIS_AXPYF_KER |
?axpyf_ft |
dotxf | BLIS_DOTXF_KER |
?dotxf_ft |
dotxaxpyf | BLIS_DOTXAXPYF_KER |
?dotxaxpyf_ft |
kernel operation | l1vkr_t | function pointer type |
---|---|---|
addv | BLIS_ADDV_KER |
?addv_ft |
amaxv | BLIS_AMAXV_KER |
?amaxv_ft |
axpyv | BLIS_AXPYV_KER |
?axpyv_ft |
axpbyv | BLIS_AXPBYV_KER |
?axpbyv_ft |
dotaxpyv | BLIS_DOTAXPYV_KER |
?dotaxpyv_ft |
copyv | BLIS_COPYV_KER |
?copyv_ft |
dotxv | BLIS_DOTXV_KER |
?dotxv_ft |
invertv | BLIS_INVERTV_KER |
?invertv_ft |
scalv | BLIS_SCALV_KER |
?scalv_ft |
scal2v | BLIS_SCAL2V_KER |
?scal2v_ft |
setv | BLIS_SETV_KER |
?setv_ft |
subv | BLIS_SUBV_KER |
?subv_ft |
swapv | BLIS_SWAPV_KER |
?swapv_ft |
xpybv | BLIS_XPBYV_KER |
?xpbyv_ft |
The specific information behind a queried function pointer is not typically available. However, it is guaranteed that the function pointer will always be valid (usually either an optimized assembly implementation or a reference implementation).
This section seeks to provide developers with a complete reference for each of the following BLIS kernels, including function prototypes, parameter descriptions, implementation notes, and diagrams:
The function prototypes in this section follow the same guidelines as those listed in the BLIS typed API reference. Namely:
- Any occurrence of
?
should be replaced withs
,d
,c
, orz
to form an actual function name. - Any occurrence of
ctype
should be replaced with the actual C type corresponding to the datatype instance in question. - Some matrix arguments have associated row and column strides arguments that proceed them, typically listed as
rsX
andcsX
for a given matrixX
. Row strides are always listed first, and column strides are always listed second. The semantic meaning of a row stride is "the distance, in units of elements, from any given element to the corresponding element (within the same column) of the next row," and the meaning of a column stride is "the distance, in units of elements, from any given element to the corresponding element (within the same row) of the next column." Thus, unit row stride implies column-major storage and unit column stride implies row-major storage. - All occurrences of
alpha
andbeta
parameters are scalars.
This section describes in detail the various level-3 microkernels supported by BLIS:
void bli_?gemm_<suffix>
(
dim_t k,
ctype* restrict alpha,
ctype* restrict a1,
ctype* restrict b1,
ctype* restrict beta,
ctype* restrict c11, inc_t rsc, inc_t csc,
auxinfo_t* restrict data,
cntx_t* restrict cntx
);
where <suffix>
is implementation-dependent. (Recall that the precise <suffix>
associated with the microkernel along with the rest of the function name doesn't matter if you are querying the function address from the context. See section on calling kernels for details.) The following (more portable) wrapper is also defined:
void bli_?gemm_ukernel
(
dim_t k,
ctype* restrict alpha,
ctype* restrict a1,
ctype* restrict b1,
ctype* restrict beta,
ctype* restrict c11, inc_t rsc, inc_t csc,
auxinfo_t* restrict data,
cntx_t* restrict cntx
);
The gemm
microkernel, sometimes simply referred to as "the BLIS microkernel" or "the microkernel", performs the following operation:
C11 := beta * C11 + A1 * B1
where A1
is an MR x k "micropanel" matrix stored in packed (column-wise) format, B1
is a k x NR "micropanel" matrix stored in packed (row-wise) format, C11
is an MR x NR general matrix stored according to its row and column strides rsc
and csc
, and alpha
and beta are scalars.
MR and NR are the register blocksizes associated with the microkernel. They are chosen by the developer when the microkernel is written and then encoded into a BLIS configuration, which will reference the microkernel when the BLIS framework is instantiated into a library. For more information on setting register blocksizes and related constants, please see the BLIS developer configuration guide.
Parameters:
k
: The number of columns ofA1
and rows ofB1
.alpha
: The address of a scalar to theA1 * B1
product.a1
: The address of a micropanel of matrixA
of dimension MR x k, stored by columns with leading dimension PACKMR, where typically PACKMR = MR. (See Implementation Notes for gemm for a discussion of PACKMR.)b1
: The address of a micropanel of matrixB
of dimension k x NR, stored by rows with leading dimension PACKNR, where typically PACKNR = NR. (See Implementation Notes for gemm for a discussion of PACKNR.)beta
: The address of a scalar to the input value of matrixC11
.c11
: The address of a matrixC11
of dimension MR x NR, stored according torsc
andcsc
.rsc
: The row stride of matrixC11
(ie: the distance to the next row, in units of matrix elements).csc
: The column stride of matrixC11
(ie: the distance to the next column, in units of matrix elements).data
: The address of anauxinfo_t
object that contains auxiliary information that may be useful when optimizing thegemm
microkernel implementation. (See Using the auxinfo_t object for a discussion of the kinds of values available viaauxinfo_t
.)cntx
: The address of the runtime context. The context can be queried for implementation-specific values such as cache and register blocksizes. However, most microkernels intrinsically "know" these values already, and thus thecntx
argument usually can be safely ignored.
The diagram below shows the packed micropanel operands and how elements of each would be stored when MR = NR = 4. The hex digits indicate the layout and order (but NOT the numeric contents) of the elements in memory. Note that the storage of C11
is not shown since it is determined by the row and column strides of C11
.
c11: a1: b1:
_______ ______________________ _______
| | |0 4 8 C | |0 1 2 3|
MR | | |1 5 9 D . . . | |4 5 6 7|
| | += |2 6 A E | |8 9 A B|
|_______| |3_7_B_F_______________| |C D E F|
| . |
NR k | . | k
| . |
| |
| |
|_______|
NR
- Register blocksizes. The register blocksizes
MR
andNR
, corresponding to the number of logical rows ina1
and columns inb1
, respectively, are defined in the context and may be queried viabli_cntx_get_blksz_def_dt()
. However, you shouldn't need to query these values since the implementation inherently "knows" them already. - Leading dimensions of
a1
andb1
: PACKMR and PACKNR. The packed micropanelsa1
andb1
are simply stored in column-major and row-major order, respectively. Usually, the width of either micropanel (ie: the number of logical rows ofa1
, or MR, and the number of columns ofb1
, or NR) is equal to that micropanel's so-called "leading dimension", or number of physical rows. Sometimes, it may be beneficial to specify a leading dimension that is larger than the panel width. This may be desirable because it allows each column ofa1
or row ofb1
to maintain a certain alignment in memory that would not otherwise be maintained by MR and/or NR. In this case, you should index througha1
andb1
using the values PACKMR and PACKNR, respectively (which are stored in the context as the blocksize "maximums" associated with thebszid_t
valuesBLIS_MR
andBLIS_NR
). These values are defined in the context and may be queried viabli_cntx_get_blksz_max_dt()
. However, you shouldn't need to query these values since the implementation inherently "knows" them already. - Storage preference of
c11
. Usually, an optimizedgemm
microkernel will have a "preferred" storage format forC11
--typically either contiguous row-storage (i.e.cs_c
= 1) or contiguous column-storage (i.e.rs_c
= 1). This preference comes from how the microkernel is most efficiently able to load/store elements ofC11
from/to memory. Most microkernels use vector instructions to access contiguous columns (or column segments) ofC11
. However, the developer may decide that accessing contiguous rows (or row segments) is more desirable. If this is the case, this preference should be indicated via thebool_t
argument when registering microkernels viabli_cntx_set_l3_nat_ukrs()
--TRUE
indicating a row preference andFALSE
indicating a column preference. Properly setting this property allows the framework to perform a runtime optimization that will ensure the microkernel preference is honored, if at all possible. - Edge cases in MR, NR dimensions. Sometimes the microkernel will be called with micropanels
a1
andb1
that correspond to edge cases, where only partial results are needed. Zero-padding is handled automatically by the packing function to facilitate reuse of the same microkernel. Similarly, the logic for computing to temporary storage and then saving only the elements that correspond to elements ofC11
that exist (at the edges) is handled automatically within the macrokernel. - Alignment of
a1
andb1
. By default, the alignment of addressesa1
andb1
are aligned only tosizeof(type)
. IfBLIS_POOL_ADDR_ALIGN_SIZE
is set to some larger multiple ofsizeof(type)
, such as the page size, then the firsta1
andb1
micropanels will be aligned to that value, but subsequent micropanels will only be aligned tosizeof(type)
, or, ifBLIS_POOL_ADDR_ALIGN_SIZE
is a multiple ofPACKMR
andPACKNR
, then subsequent micropanelsa1
andb1
will be aligned toPACKMR * sizeof(type)
andPACKNR * sizeof(type)
, respectively. - Unrolling loops. As a general rule of thumb, the loop over k is sometimes moderately unrolled; for example, in our experience, an unrolling factor of u = 4 is fairly common. If unrolling is applied in the k dimension, edge cases must be handled to support values of k that are not multiples of u. It is nearly universally true that there should be no loops in the MR or NR directions; in other words, iteration over these dimensions should always be fully unrolled (within the loop over k).
- Zero
beta
. Ifbeta
= 0.0 (or 0.0 + 0.0i for complex datatypes), then the microkernel should NOT use it explicitly, asC11
may contain uninitialized memory (including elements containingNaN
orInf
). This case should be detected and handled separately by overwritingC11
with thealpha * A1 * B1
product.
Each microkernel (gemm, trsm, and gemmtrsm) takes as its last argument a pointer of type auxinfo_t
. This BLIS-defined type is defined as a struct
whose fields contain auxiliary values that may be useful to some microkernel authors, particularly when implementing certain optimization techniques. BLIS provides kernel authors access to the fields of the auxinfo_t
object via the following function-like preprocessor macros. Each macro takes a single argument, the auxinfo_t
pointer, and returns one of the values stored within the object.
bli_auxinfo_next_a()
. Returns the address (void*
) of the micropanel ofA
that will be used the next time the microkernel will be called.bli_auxinfo_next_b()
. Returns the address (void*
) of the micropanel ofB
that will be used the next time the microkernel will be called.bli_auxinfo_ps_a()
. Returns the panel stride (inc_t
) of the current micropanel ofA
.bli_auxinfo_ps_b()
. Returns the panel stride (inc_t
) of the current micropanel ofB
.
The addresses of the next micropanels of A
and B
may be used by the microkernel to perform prefetching, if prefetching is supported by the architecture. Similarly, it may be useful to know the precise distance in memory to the next micropanel. (Note that sometimes the next micropanel to be used is not the same as the next micropanel in memory.)
Any and all of these values may be safely ignored; they are completely optional. However, BLIS guarantees that all values accessed via the macros listed above will always be initialized and meaningful, for every invocation of each microkernel (gemm
, trsm
, and gemmtrsm
).
An example implementation of the gemm
microkernel may be found in the template
configuration directory in:
Note that this implementation is coded in C99 and lacks several kinds of optimization that are typical of real-world optimized microkernels, such as vector instructions (or intrinsics) and loop unrolling in MR or NR. It is meant to serve only as a starting point for a microkernel developer.
void bli_?trsm_l_<suffix>
(
ctype* restrict a11,
ctype* restrict b11,
ctype* restrict c11, inc_t rsc, inc_t csc,
auxinfo_t* restrict data,
cntx_t* restrict cntx
);
void bli_?trsm_u_<suffix>
(
ctype* restrict a11,
ctype* restrict b11,
ctype* restrict c11, inc_t rsc, inc_t csc,
auxinfo_t* restrict data,
cntx_t* restrict cntx
);
where <suffix>
is implementation-dependent. (Recall that the precise <suffix>
associated with the microkernel along with the rest of the function name doesn't matter if you are querying the function address from the context. See section on calling kernels for details.) The following (more portable) wrappers are also defined:
void bli_?trsm_l_ukernel
(
ctype* restrict a11,
ctype* restrict b11,
ctype* restrict c11, inc_t rsc, inc_t csc,
auxinfo_t* restrict data,
cntx_t* restrict cntx
);
void bli_?trsm_u_ukernel
(
ctype* restrict a11,
ctype* restrict b11,
ctype* restrict c11, inc_t rsc, inc_t csc,
auxinfo_t* restrict data,
cntx_t* restrict cntx
);
The trsm_l
and trsm_u
microkernels perform the following operation:
C11 := inv(A11) * B11
where A11
is MR x MR and lower (trsm_l
) or upper (trsm_u
) triangular, B11
is MR x NR, and C11
is MR x NR.
MR and NR are the register blocksizes associated with the microkernel. They are chosen by the developer when the microkernel is written and then encoded into a BLIS configuration, which will reference the microkernel when the BLIS framework is instantiated into a library. For more information on setting register blocksizes and related constants, please see the BLIS developer configuration guide.
Parameters:
a11
: The address ofA11
, which is the MR x MR lower (trsm_l
) or upper (trsm_u
) triangular submatrix within the packed micropanel of matrixA
.A11
is stored by columns with leading dimension PACKMR, where typically PACKMR = MR. (See Implementation Notes for gemm for a discussion of PACKMR.) Note thatA11
contains elements in both triangles, though elements in the unstored triangle are not guaranteed to be zero and thus should not be referenced.b11
: The address ofB11
, which is an MR x NR submatrix of the packed micropanel ofB
.B11
is stored by rows with leading dimension PACKNR, where typically PACKNR = NR. (See Implementation Notes for gemm for a discussion of PACKNR.)c11
: The address ofC11
, which is an MR x NR submatrix of matrixC
, stored according torsc
andcsc
.C11
is the submatrix withinC
that corresponds to the elements which were packed intoB11
. Thus,C
is the original input matrixB
to the overalltrsm
operation.rsc
: The row stride of matrixC11
(ie: the distance to the next row, in units of matrix elements).csc
: The column stride of matrixC11
(ie: the distance to the next column, in units of matrix elements).data
: The address of anauxinfo_t
object that contains auxiliary information that may be useful when optimizing thetrsm
microkernel implementation. (See Using the auxinfo_t object for a discussion of the kinds of values available viaauxinfo_t
, and also Implementation Notes for trsm for caveats.)cntx
: The address of the runtime context. The context can be queried for implementation-specific values such as cache and register blocksizes. However, most microkernels intrinsically "know" these values already, and thus thecntx
argument usually can be safely ignored.
Please see the diagram for gemmtrsm_l and gemmtrsm_u to see depictions of the trsm_l
and trsm_u
microkernel operations and where they fit in with their preceding gemm
subproblems.
- Register blocksizes. See Implementation Notes for gemm.
- Leading dimensions of
a11
andb11
: PACKMR and PACKNR. See Implementation Notes for gemm. - Edge cases in MR, NR dimensions. See Implementation Notes for gemm.
- Alignment of
a11
andb11
. The addressesa11
andb11
are aligned according toPACKMR * sizeof(type)
andPACKNR * sizeof(type)
, respectively. - Unrolling loops. Most optimized implementations should unroll all three loops within the
trsm
microkernel. - Prefetching next micropanels of
A
andB
. We advise against using thebli_auxinfo_next_a()
andbli_auxinfo_next_b()
macros from within thetrsm_l
andtrsm_u
microkernels, since the values returned usually only make sense in the context of the overallgemmtrsm
subproblem. - Diagonal elements of
A11
. At the time this microkernel is called, the diagonal entries of triangular matrixA11
contain the inverse of the original elements. This inversion is done during packing so that we can avoid expensive division instructions within the microkernel itself. If thediag
parameter to the higher leveltrsm
operation was equal toBLIS_UNIT_DIAG
, the diagonal elements will be explicitly unit. - Zero elements of
A11
. SinceA11
is lower triangular (fortrsm_l
), the strictly upper triangle implicitly contains zeros. Similarly, the strictly lower triangle ofA11
implicitly contains zeros whenA11
is upper triangular (fortrsm_u
). However, the packing function may or may not actually write zeros to this region. Thus, the implementation should not reference these elements. - Output. This microkernel must write its result to two places: the submatrix
B11
of the current packed micropanel ofB
and the submatrixC11
of the output matrixC
.
Example implementations of the trsm
microkernels may be found in the template
configuration directory in:
Note that these implementations are coded in C99 and lack several kinds of optimization that are typical of real-world optimized microkernels, such as vector instructions (or intrinsics) and loop unrolling in MR or NR. They are meant to serve only as a starting point for a microkernel developer.
void bli_?gemmtrsm_l_<suffix>
(
dim_t k,
ctype* restrict alpha,
ctype* restrict a10,
ctype* restrict a11,
ctype* restrict b01,
ctype* restrict b11,
ctype* restrict c11, inc_t rsc, inc_t csc,
auxinfo_t* restrict data,
cntx_t* restrict cntx
);
void bli_?gemmtrsm_u_<suffix>
(
dim_t k,
ctype* restrict alpha,
ctype* restrict a12,
ctype* restrict a11,
ctype* restrict b21,
ctype* restrict b11,
ctype* restrict c11, inc_t rsc, inc_t csc,
auxinfo_t* restrict data,
cntx_t* restrict cntx
);
where <suffix>
is implementation-dependent. (Recall that the precise <suffix>
associated with the microkernel along with the rest of the function name doesn't matter if you are querying the function address from the context. See section on calling kernels for details.) The following (more portable) wrappers are also defined:
void bli_?gemmtrsm_l_ukernel
(
dim_t k,
ctype* restrict alpha,
ctype* restrict a10,
ctype* restrict a11,
ctype* restrict b01,
ctype* restrict b11,
ctype* restrict c11, inc_t rsc, inc_t csc,
auxinfo_t* restrict data,
cntx_t* restrict cntx
);
void bli_?gemmtrsm_u_ukernel
(
dim_t k,
ctype* restrict alpha,
ctype* restrict a12,
ctype* restrict a11,
ctype* restrict b21,
ctype* restrict b11,
ctype* restrict c11, inc_t rsc, inc_t csc,
auxinfo_t* restrict data,
cntx_t* restrict cntx
);
The gemmtrsm_l
microkernel performs the following compound operation:
B11 := alpha * B11 - A10 * B01
B11 := inv(A11) * B11
C11 := B11
where A11
is MR x MR and lower triangular, A10
is MR x k, and B01
is k x NR.
The gemmtrsm_u
microkernel performs:
B11 := alpha * B11 - A12 * B21
B11 := inv(A11) * B11
C11 := B11
where A11
is MR x MR and upper triangular, A12
is MR x k, and B21
is k x NR.
In both cases, B11
is MR x NR and alpha
is a scalar. Here, inv()
denotes matrix inverse.
MR and NR are the register blocksizes associated with the microkernel. They are chosen by the developer when the microkernel is written and then encoded into a BLIS configuration, which will reference the microkernel when the BLIS framework is instantiated into a library. For more information on setting register blocksizes and related constants, please see the BLIS developer configuration guide.
Parameters:
k
: The number of columns ofA10
and rows ofB01
(trsm_l
); the number of columns ofA12
and rows ofB21
(trsm_u
).alpha
: The address of a scalar to be applied toB11
.a10
,a12
: The address ofA10
orA12
, which is the MR x k submatrix of the packed micropanel ofA
that is situated to the left (trsm_l
) or right (trsm_u
) of the MR x MR triangular submatrixA11
.A10
andA12
are stored by columns with leading dimension PACKMR, where typically PACKMR = MR. (See Implementation Notes for gemm for a discussion of PACKMR.)a11
: The address ofA11
, which is the MR x MR lower (trsm_l
) or upper (trsm_u
) triangular submatrix within the packed micropanel of matrixA
that is situated to the right ofA10
(trsm_l
) or the left ofA12
(trsm_u
).A11
is stored by columns with leading dimension PACKMR, where typically PACKMR = MR. (See Implementation Notes for gemm for a discussion of PACKMR.) Note thatA11
contains elements in both triangles, though elements in the unstored triangle are not guaranteed to be zero and thus should not be referenced.b01
,b21
: The address ofB01
andB21
, which is the k x NR submatrix of the packed micropanel ofB
that is situated above (trsm_l
) or below (trsm_u
) the MR x NR blockB11
.B01
andB21
are stored by rows with leading dimension PACKNR, where typically PACKNR = NR. (See Implementation Notes for gemm for a discussion of PACKNR.)b11
: The address ofB11
, which is the MR x NR submatrix of the packed micropanel ofB
, situated belowB01
(trsm_l
) or aboveB21
(trsm_u
).B11
is stored by rows with leading dimension PACKNR, where typically PACKNR = NR. (See Implementation Notes for gemm for a discussion of PACKNR.)c11
: The address ofC11
, which is an MR x NR submatrix of matrixC
, stored according torsc
andcsc
.C11
is the submatrix withinC
that corresponds to the elements which were packed intoB11
. Thus,C
is the original input matrixB
to the overalltrsm
operation.rsc
: The row stride of matrixC11
(ie: the distance to the next row, in units of matrix elements).csc
: The column stride of matrixC11
(ie: the distance to the next column, in units of matrix elements).data
: The address of anauxinfo_t
object that contains auxiliary information that may be useful when optimizing thegemmtrsm
microkernel implementation. (See Using the auxinfo_t object for a discussion of the kinds of values available viaauxinfo_t
, and also Implementation Notes for gemmtrsm for caveats.)cntx
: The address of the runtime context. The context can be queried for implementation-specific values such as cache and register blocksizes. However, most microkernels intrinsically "know" these values already, and thus thecntx
argument usually can be safely ignored.
The diagram below shows the packed micropanel operands for trsm_l
and how elements of each would be stored when MR = NR = 4. (The hex digits indicate the layout and order (but NOT the numeric contents) in memory. Here, matrix A11
(referenced by a11
) is lower triangular. Matrix A11
does contain elements corresponding to the strictly upper triangle, however, they are not guaranteed to contain zeros and thus these elements should not be referenced.
NR
_______
b01:|0 1 2 3|
|4 5 6 7|
|8 9 A B|
|C D E F|
k | . |
| . |
a10: a11: | . |
___________________ _______ |_______|
|0 4 8 C |`. | b11:| |
MR |1 5 9 D . . . | `. | | |
|2 6 A E | `. | MR | |
|3_7_B_F____________|______`.| |_______|
k MR
The diagram below shows the packed micropanel operands for trsm_u
and how elements of each would be stored when MR = NR = 4. (The hex digits indicate the layout and order (but NOT the numeric contents) in memory. Here, matrix A11
(referenced by a11
) is upper triangular. Matrix A11
does contain elements corresponding to the strictly lower triangle, however, they are not guaranteed to contain zeros and thus these elements should not be referenced.
a11: a12: NR
________ ___________________ _______
|`. |0 4 8 | b11:|0 1 2 3|
MR | `. |1 5 9 . . . | |4 5 6 7|
| `. |2 6 A | MR |8 9 A B|
|______`.|3_7_B______________| |___.___|
b21:| . |
MR k | . |
| |
| |
NOTE: Storage digits are shown k | |
starting with a12 to avoid | |
obscuring triangular structure | |
of a11. |_______|
- Register blocksizes. See Implementation Notes for gemm.
- Leading dimensions of
a1
andb1
: PACKMR and PACKNR. See Implementation Notes for gemm. - Edge cases in MR, NR dimensions. See Implementation Notes for gemm.
- Alignment of
a1
andb1
. See Implementation Notes for gemm. - Unrolling loops. Most optimized implementations should unroll all three loops within the
trsm
subproblem ofgemmtrsm
. See Implementation Notes for gemm for remarks on unrolling thegemm
subproblem. - Prefetching next micropanels of
A
andB
. When invoked from within agemmtrsm_l
microkernel, the addresses accessible viabli_auxinfo_next_a()
andbli_auxinfo_next_b()
refer to the next invocation'sa10
andb01
, respectively, while ingemmtrsm_u
, the_next_a()
and_next_b()
macros return the addresses of the next invocation'sa11
andb11
(since those submatrices precedea12
andb21
). - Zero
alpha
. The microkernel can safely assume thatalpha
is non-zero; "alpha equals zero" handling is performed at a much higher level, which means that, in such a scenario, the microkernel will never get called. - Diagonal elements of
A11
. See Implementation Notes for trsm. - Zero elements of
A11
. See Implementation Notes for trsm. - Output. See Implementation Notes for trsm.
- Optimization. Let's assume that the gemm microkernel has already been optimized. You have two options with regard to optimizing the fused
gemmtrsm
microkernels:- Optimize only the trsm microkernels. This will result in the
gemm
andtrsm_l
microkernels being called in sequence. (Likewise forgemm
andtrsm_u
.) - Fuse the implementation of the
gemm
microkernel with that of thetrsm
microkernels by inlining both into thegemmtrsm_l
andgemmtrsm_u
microkernel definitions. This option is more labor-intensive, but also more likely to yield higher performance because it avoids redundant memory operations on the packed MR x NR submatrixB11
.
- Optimize only the trsm microkernels. This will result in the
Example implementations of the gemmtrsm
microkernels may be found in the template
configuration directory in:
- config/template/kernels/3/bli_gemmtrsm_l_opt_mxn.c
- config/template/kernels/3/bli_gemmtrsm_u_opt_mxn.c
Note that these implementations are coded in C99 and lack several kinds of optimization that are typical of real-world optimized microkernels, such as vector instructions (or intrinsics) and loop unrolling in MR or NR. They are meant to serve only as a starting point for a microkernel developer.
void bli_?axpy2v_<suffix>
(
conj_t conjx,
conj_t conjy,
dim_t n,
ctype* restrict alphax,
ctype* restrict alphay,
ctype* restrict x, inc_t incx,
ctype* restrict y, inc_t incy,
ctype* restrict z, inc_t incz,
cntx_t* restrict cntx
)
This kernel performs the following operation:
z := z + alphax * conjx(x) + alphay * conjy(y)
where x
, y
, and z
are vectors of length n stored with strides incx
, incy
, and incz
, respectively. This kernel is typically implemented as the fusion of two axpyv
operations on different input vectors x
and y
and with different scalars alphax
and alpay
to update the same output vector z
.
void bli_?dotaxpyv_<suffix>
(
conj_t conjxt,
conj_t conjx,
conj_t conjy,
dim_t n,
ctype* restrict alpha,
ctype* restrict x, inc_t incx,
ctype* restrict y, inc_t incy,
ctype* restrict rho,
ctype* restrict z, inc_t incz,
cntx_t* restrict cntx
)
This kernel performs the following operation:
rho := conjxt(x)^T * conjy(y)
z := z + alpha * conjx(x)
where x
, y
, and z
are vectors of length n stored with strides incx
, incy
, and incz
, respectively, and rho
is a scalar. This kernel is typically implemented as a dotv
operation fused with an axpyv
operation.
void bli_?axpyf_<suffix>
(
conj_t conja,
conj_t conjx,
dim_t m,
dim_t b,
ctype* restrict alpha,
ctype* restrict a, inc_t inca, inc_t lda,
ctype* restrict x, inc_t incx,
ctype* restrict y, inc_t incy,
cntx_t* restrict cntx
)
This kernel performs the following operation:
y := y + alpha * conja(a) * conjy(x)
where a
is an m x b matrix, x
is a vector of length b, and y
is a vector of length m. Vectors x
and y
are stored with strides incx
and incy
, respectively. Matrix a
is stored with row stride inca
and column stride lda
, though inca
is most often (in practice) unit. This kernel is typically implemented as a fused series of b axpyv
operations updating the same vector y
(with the elements of x
serving as the scalars and the columns of a
serving as the vectors to be scaled).
void bli_?dotxf_<suffix>
(
conj_t conjat,
conj_t conjx,
dim_t m,
dim_t b,
ctype* restrict alpha,
ctype* restrict a, inc_t inca, inc_t lda,
ctype* restrict x, inc_t incx,
ctype* restrict beta,
ctype* restrict y, inc_t incy,
cntx_t* restrict cntx
)
This kernel performs the following operation:
y := beta * y + alpha * conjat(a)^T conjx(x)
where a
is an m x b matrix, where w
is a vector of length m, y
is a vector of length b, and alpha
is a scalar.
Vectors x
and y
are stored with strides incx
and incy
, respectively. Matrix a
is stored with row stride inca
and column stride lda
, though inca
is most often (in practice) unit.
This kernel is typically implemented as a series of b dotxv
operations with the same right-hand operand vector x
(contracted with the rows of a^T
and accumulating to the corresponding elements of vector y
).
void bli_?dotxaxpyf_<suffix>
(
conj_t conjat,
conj_t conja,
conj_t conjw,
conj_t conjx,
dim_t m,
dim_t b,
ctype* restrict alpha,
ctype* restrict a, inc_t inca, inc_t lda,
ctype* restrict w, inc_t incw,
ctype* restrict x, inc_t incx,
ctype* restrict beta,
ctype* restrict y, inc_t incy,
ctype* restrict z, inc_t incz,
cntx_t* restrict cntx
)
This kernel performs the following operation:
y := beta * y + alpha * conjat(a)^T conjw(w)
z := z + alpha * conja(a) conjx(x)
where a
is an m x b matrix, w
and z
are vectors of length m, x
and y
are vectors of length b, and alpha
and beta
are scalars.
Vectors w
, z
, x
and y
are stored with strides incw
, incz
, incx
, and incy
, respectively. Matrix a
is stored with row stride inca
and column stride lda
, though inca
is most often (in practice) unit.
This kernel is typically implemented as a series of b dotxv
operations with the same right-hand operand vector w
fused with a series of b axpyv
operations updating the same vector z
.
void bli_?addv_<suffix>
(
conj_t conjx,
dim_t n,
ctype* restrict x, inc_t incx,
ctype* restrict y, inc_t incy,
cntx_t* restrict cntx
)
This kernel performs the following operation:
y := y + conjx(x)
where x
and y
are vectors of length n stored with strides incx
and incy
, respectively.
void bli_?amaxv_<suffix>
(
dim_t n,
ctype* restrict x, inc_t incx,
dim_t* restrict index,
cntx_t* restrict cntx
)
Given a vector of length n, this kernel returns the zero-based index index
of the element of vector x
that contains the largest absolute value (or, in the complex domain, the largest complex modulus).
If NaN
is encountered, it is treated as if it were a valid value that was smaller than any other value in the vector. If more than one element contains the same maximum value, the index of the latter element is returned via index
.
void bli_?axpyv_<suffix>
(
conj_t conjx,
dim_t n,
ctype* restrict alpha,
ctype* restrict x, inc_t incx,
ctype* restrict y, inc_t incy,
cntx_t* restrict cntx
)
This kernel performs the following operation:
y := y + alpha * conjx(x)
where x
and y
are vectors of length n stored with strides incx
and incy
, respectively, and alpha
is a scalar.
void bli_?axpbyv_<suffix>
(
conj_t conjx,
dim_t n,
ctype* restrict alpha,
ctype* restrict x, inc_t incx,
ctype* restrict beta,
ctype* restrict y, inc_t incy,
cntx_t* restrict cntx
)
This kernel performs the following operation:
y := beta * y + alpha * conjx(x)
where x
and y
are vectors of length n stored with strides incx
and incy
, respectively, and alpha
and beta
are scalars.
void bli_?copyv_<suffix>
(
conj_t conjx,
dim_t n,
ctype* restrict x, inc_t incx,
ctype* restrict y, inc_t incy,
cntx_t* restrict cntx
)
This kernel performs the following operation:
y := conjx(x)
where x
and y
are vectors of length n stored with strides incx
and incy
, respectively.
void bli_?dotv_<suffix>
(
conj_t conjx,
conj_t conjy,
dim_t n,
ctype* restrict x, inc_t incx,
ctype* restrict y, inc_t incy,
ctype* restrict rho,
cntx_t* restrict cntx
)
This kernel performs the following operation:
rho := conjxt(x)^T * conjy(y)
where x
and y
are vectors of length n stored with strides incx
and incy
, respectively, and rho
is a scalar.
void bli_?dotxv_<suffix>
(
conj_t conjx,
conj_t conjy,
dim_t n,
ctype* restrict alpha,
ctype* restrict x, inc_t incx,
ctype* restrict y, inc_t incy,
ctype* restrict beta,
ctype* restrict rho,
cntx_t* restrict cntx
)
This kernel performs the following operation:
rho := beta * rho + alpha * conjxt(x)^T * conjy(y)
where x
and y
are vectors of length n stored with strides incx
and incy
, respectively, and alpha
, beta
, and rho
are scalars.
void bli_?invertv_<suffix>
(
dim_t n,
ctype* restrict x, inc_t incx,
cntx_t* restrict cntx
)
This kernel inverts all elements of an n-length vector x
.
void bli_?scalv_<suffix>
(
conj_t conjalpha,
dim_t n,
ctype* restrict alpha,
ctype* restrict x, inc_t incx,
cntx_t* restrict cntx
)
This kernel performs the following operation:
x := conjalpha(alpha) * x
where x
is a vector of length n stored with stride incx
and alpha
is a scalar.
void bli_?scal2v_<suffix>
(
conj_t conjx,
dim_t n,
ctype* restrict alpha,
ctype* restrict x, inc_t incx,
ctype* restrict y, inc_t incy,
cntx_t* restrict cntx
)
This kernel performs the following operation:
y := alpha * conjx(x)
where x
and y
are vectors of length n stored with strides incx
and incy
, respectively, and alpha
is a scalar.
void bli_?setv_<suffix>
(
conj_t conjalpha,
dim_t n,
ctype* restrict alpha,
ctype* restrict x, inc_t incx,
cntx_t* restrict cntx
)
This kernel performs the following operation:
x := conjalpha(alpha)
where x
is a vector of length n stored with stride incx
and alpha
is a scalar. Note that here, the :=
operator represents a broadcast of conjalpha(alpha)
to every element in x
.
void bli_?subv_<suffix>
(
conj_t conjx,
dim_t n,
ctype* restrict x, inc_t incx,
ctype* restrict y, inc_t incy,
cntx_t* restrict cntx
)
This kernel performs the following operation:
y := y - conjx(x)
where x
and y
are vectors of length n.
void bli_?swapv_<suffix>
(
dim_t n,
ctype* restrict x, inc_t incx,
ctype* restrict y, inc_t incy,
cntx_t* restrict cntx
)
This kernel swaps corresponding elements of two n-length vectors x
and y
stored with strides incx
and incy
, respectively.
void bli_?xpbyv_<suffix>
(
conj_t conjx,
dim_t n,
ctype* restrict x, inc_t incx,
ctype* restrict beta,
ctype* restrict y, inc_t incy,
cntx_t* restrict cntx
)
This kernel performs the following operation:
y := beta * y + conjx(x)
where x
and y
are vectors of length n stored with strides incx
and incy
, respectively, and beta
is a scalar.