/*
 * Copyright (c) IBM Corporation 2020.
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are
 * met:
 *
 *    1. Redistributions of source code must retain the above copyright
 *       notice, this list of conditions and the following disclaimer.
 *
 *    2. Redistributions in binary form must reproduce the above copyright
 *       notice, this list of conditions and the following disclaimer in
 *       the documentation and/or other materials provided with the
 *       distribution.
 *    3. Neither the name of the OpenBLAS project nor the names of
 *       its contributors may be used to endorse or promote products
 *       derived from this software without specific prior written
 *       permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
 * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
 * USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */
#include "common.h"
#include "vector-common.h"

#include <stdbool.h>
#include <stdio.h>
#include <stdlib.h>


#ifdef COMPLEX
#error "Handling for complex numbers is not supported in this kernel"
#endif

#ifdef DOUBLE
#define UNROLL_M DGEMM_DEFAULT_UNROLL_M
#define UNROLL_N DGEMM_DEFAULT_UNROLL_N
#else
#define UNROLL_M SGEMM_DEFAULT_UNROLL_M
#define UNROLL_N SGEMM_DEFAULT_UNROLL_N
#endif

static const size_t unroll_m = UNROLL_M;
static const size_t unroll_n = UNROLL_N;

/* Handling of triangular matrices */
#ifdef TRMMKERNEL
static const bool trmm = true;
static const bool left =
#ifdef LEFT
	true;
#else
	false;
#endif

static const bool backwards =
#if defined(LEFT) != defined(TRANSA)
	true;
#else
	false;
#endif

#else
static const bool trmm = false;
static const bool left = false;
static const bool backwards = false;
#endif /* TRMMKERNEL */

/*
 * Background:
 *
 * The algorithm of GotoBLAS / OpenBLAS breaks down the matrix multiplication
 * problem by splitting all matrices into partitions multiple times, so that the
 * submatrices fit into the L1 or L2 caches. As a result, each multiplication of
 * submatrices can stream data fast from L1 and L2 caches. Inbetween, it copies
 * and rearranges the submatrices to enable contiguous memory accesses to
 * improve locality in both caches and TLBs.
 *
 * At the heart of the algorithm is this kernel, which multiplies, a "Block
 * matrix" A (small dimensions) with a "Panel matrix" B (number of rows is
 * small) and adds the result into a "Panel matrix" C; GotoBLAS calls this
 * operation GEBP. This kernel further partitions GEBP twice, such that (1)
 * submatrices of C and B fit into the L1 caches (GEBP_column_block) and (2) a
 * block of C fits into the registers, while multiplying panels from A and B
 * streamed from the L2 and L1 cache, respectively (GEBP_block).
 *
 *
 * Algorithm GEBP(A, B, C, m, n, k, alpha):
 *
 * The problem is calculating C += alpha * (A * B)
 * C is an m x n matrix, A is an m x k matrix, B is an k x n matrix.
 *
 * - C is in column-major-order, with an offset of ldc to the element in the
 *   next column (same row).
 * - A is in row-major-order yet stores SGEMM_UNROLL_M elements of each column
 *   contiguously while walking along rows.
 * - B is in column-major-order but packs SGEMM_UNROLL_N elements of a row
 *   contiguously.
 * If the numbers of rows and columns are not multiples of SGEMM_UNROLL_M or
 * SGEMM_UNROLL_N, the remaining elements are arranged in blocks with power-of-2
 * dimensions (e.g., 5 remaining columns would be in a block-of-4 and a
 * block-of-1).
 *
 * Note that packing A and B into that form is taken care of by the caller in
 * driver/level3/level3.c (actually done by "copy kernels").
 *
 * Steps:
 * - Partition C and B into blocks of n_r (SGEMM_UNROLL_N) columns, C_j and B_j.
 *   Now, B_j should fit into the L1 cache.
 * - For each partition, calculate C_j += alpha * (A * B_j) by
 *     (1) Calculate C_aux := A * B_j (see below)
 *     (2) unpack C_j = C_j + alpha * C_aux
 *
 *
 * Algorithm for Calculating C_aux:
 *
 * - Further partition C_aux and A into groups of m_r (SGEMM_UNROLL_M) rows,
 *   such that the m_r x n_r-submatrix of C_aux can be held in registers. Each
 *   submatrix of C_aux can be calculated independently, and the registers are
 *   added back into C_j.
 *
 * - For each row-block of C_aux:
 *   (uses a row block of A and full B_j)
 *    - stream over all columns of A, multiply with elements from B and
 *      accumulate in registers. (use different inner-kernels to exploit
 *      vectorization for varying block sizes)
 *    - add alpha * row block of C_aux back into C_j.
 *
 * Note that there are additional mechanics for handling triangular matrices,
 * calculating B := alpha (A * B) where either of the matrices A or B can be
 * triangular. In case of A, the macro "LEFT" is defined. In addition, A can
 * optionally be transposed.
 * The code effectively skips an "offset" number of columns in A and rows of B
 * in each block, to save unnecessary work by exploiting the triangular nature.
 * To handle all cases, the code discerns (1) a "left" mode when A is triangular
 * and (2) "forward" / "backwards" modes where only the first "offset"
 * columns/rows of A/B are used or where the first "offset" columns/rows are
 * skipped, respectively.
 *
 * Reference:
 *
 * The summary above is based on staring at various kernel implementations and:
 * K. Goto and R. A. Van de Geijn, Anatomy of High-Performance Matrix
 * Multiplication, in ACM Transactions of Mathematical Software, Vol.  34, No.
 * 3, May 2008.
 */

/**
 * Calculate for a row-block in C_i of size ROWSxCOLS using vector intrinsics.
 *
 * @param[in] 	A	Pointer current block of input matrix A.
 * @param[in]	k	Number of columns in A.
 * @param[in]	B	Pointer current block of input matrix B.
 * @param[inout] C	Pointer current block of output matrix C.
 * @param[in]	ldc	Offset between elements in adjacent columns in C.
 * @param[in]	alpha	Scalar factor.
 */
#define VECTOR_BLOCK(ROWS, COLS)                                              \
	static inline void GEBP_block_##ROWS##_##COLS(                        \
	    FLOAT const *restrict A, BLASLONG bk, FLOAT const *restrict B,    \
	    FLOAT *restrict C, BLASLONG ldc, FLOAT alpha) {                   \
		_Static_assert(                                               \
		    ROWS % VLEN_FLOATS == 0,                                  \
		    "rows in block must be multiples of vector length");      \
		vector_float Caux[ROWS / VLEN_FLOATS][COLS];                  \
                                                                              \
		for (BLASLONG i = 0; i < ROWS / VLEN_FLOATS; i++) {           \
			vector_float A0 =                                     \
			    vec_load_hinted(A + i * VLEN_FLOATS);             \
			for (BLASLONG j = 0; j < COLS; j++)                   \
				Caux[i][j] = A0 * B[j];                       \
		}                                                             \
                                                                              \
		/*                                                            \
		 * Stream over the row-block of A, which is packed            \
		 * column-by-column, multiply by coefficients in B and add up \
		 * into temporaries Caux (which the compiler will hold in     \
		 * registers). Vectorization: Multiply column vectors from A  \
		 * with scalars from B and add up in column vectors of Caux.  \
		 * That equates to unrolling the loop over rows (in i) and    \
		 * executing each unrolled iteration as a vector element.     \
		 */                                                           \
		for (BLASLONG k = 1; k < bk; k++) {                           \
			for (BLASLONG i = 0; i < ROWS / VLEN_FLOATS; i++) {   \
				vector_float Ak = vec_load_hinted(            \
				    A + i * VLEN_FLOATS + k * ROWS);          \
                                                                              \
				for (BLASLONG j = 0; j < COLS; j++)           \
					Caux[i][j] += Ak * B[j + k * COLS];   \
			}                                                     \
		}                                                             \
                                                                              \
		/*                                                            \
		 * Unpack row-block of C_aux into outer C_i, multiply by      \
		 * alpha and add up.                                          \
		 */                                                           \
		for (BLASLONG j = 0; j < COLS; j++) {                         \
			for (BLASLONG i = 0; i < ROWS / VLEN_FLOATS; i++) {   \
				vector_float *C_ij =                          \
				    (vector_float *)(C + i * VLEN_FLOATS +    \
						     j * ldc);                \
				if (trmm) {                                   \
					*C_ij = alpha * Caux[i][j];           \
				} else {                                      \
					*C_ij += alpha * Caux[i][j];          \
				}                                             \
			}                                                     \
		}                                                             \
	}


#if UNROLL_M == 16
VECTOR_BLOCK(16, 2)
VECTOR_BLOCK(16, 1)
#endif
#if UNROLL_N == 8
VECTOR_BLOCK(8, 8)
VECTOR_BLOCK(4, 8)
#endif
#ifndef DOUBLE
VECTOR_BLOCK(8, 4)
#endif
VECTOR_BLOCK(8, 2)
VECTOR_BLOCK(8, 1)
VECTOR_BLOCK(4, 4)
VECTOR_BLOCK(4, 2)
VECTOR_BLOCK(4, 1)

/**
 * Calculate for a row-block in C_i of size ROWSxCOLS using scalar operations.
 * Simple implementation for smaller block sizes
 *
 * @param[in] 	A	Pointer current block of input matrix A.
 * @param[in]	k	Number of columns in A.
 * @param[in]	B	Pointer current block of input matrix B.
 * @param[inout] C	Pointer current block of output matrix C.
 * @param[in]	ldc	Offset between elements in adjacent columns in C.
 * @param[in]	alpha	Scalar factor.
 */
#define SCALAR_BLOCK(ROWS, COLS)                                          \
    static inline void GEBP_block_##ROWS##_##COLS(                        \
	FLOAT const *restrict A, BLASLONG k, FLOAT const *restrict B,     \
	FLOAT *restrict C, BLASLONG ldc, FLOAT alpha) {                   \
	FLOAT Caux[ROWS][COLS] __attribute__((aligned(16)));              \
                                                                          \
	/*                                                                \
	 * Peel off first iteration (i.e., column of A) for               \
	 * initializing Caux                                              \
	 */                                                               \
	for (BLASLONG i = 0; i < ROWS; i++)                               \
	    for (BLASLONG j = 0; j < COLS; j++) Caux[i][j] = A[i] * B[j]; \
                                                                          \
	for (BLASLONG kk = 1; kk < k; kk++)                               \
	    for (BLASLONG i = 0; i < ROWS; i++)                           \
		for (BLASLONG j = 0; j < COLS; j++)                       \
		    Caux[i][j] += A[i + kk * ROWS] * B[j + kk * COLS];    \
                                                                          \
	for (BLASLONG i = 0; i < ROWS; i++)                               \
	    for (BLASLONG j = 0; j < COLS; j++)                           \
		if (trmm) {                                               \
		    C[i + j * ldc] = alpha * Caux[i][j];                  \
		} else {                                                  \
		    C[i + j * ldc] += alpha * Caux[i][j];                 \
		}                                                         \
    }

#ifdef DOUBLE
VECTOR_BLOCK(2, 4)
VECTOR_BLOCK(2, 2)
VECTOR_BLOCK(2, 1)
#else
SCALAR_BLOCK(2, 4)
SCALAR_BLOCK(2, 2)
SCALAR_BLOCK(2, 1)
#endif

SCALAR_BLOCK(1, 4)
SCALAR_BLOCK(1, 2)
SCALAR_BLOCK(1, 1)


/**
 * Calculate a row-block that fits 4x4 vector registers using a loop
 * unrolled-by-2 with explicit interleaving to better overlap loads and
 * computation.
 * This function fits 16x4 blocks for SGEMM and 8x4 blocks for DGEMM.
 */
#ifdef DOUBLE
static inline void GEBP_block_8_4(
#else // float
static inline void GEBP_block_16_4(
#endif
    FLOAT const *restrict A, BLASLONG bk, FLOAT const *restrict B,
    FLOAT *restrict C, BLASLONG ldc, FLOAT alpha) {
#define VEC_ROWS 4
#define VEC_COLS 4
#define ROWS VEC_ROWS * VLEN_FLOATS
#define COLS (VEC_COLS)

	/*
	 * Hold intermediate results in vector registers.
	 * Since we need to force the compiler's hand in places, we need to use
	 * individual variables in contrast to the generic implementation's
	 * arrays.
	 */
#define INIT_ROW_OF_C(ROW)                                        \
    vector_float A##ROW = vec_load_hinted(A + ROW * VLEN_FLOATS); \
    vector_float C_##ROW##_0 = A##ROW * B[0];                     \
    vector_float C_##ROW##_1 = A##ROW * B[1];                     \
    vector_float C_##ROW##_2 = A##ROW * B[2];                     \
    vector_float C_##ROW##_3 = A##ROW * B[3];

	INIT_ROW_OF_C(0)
	INIT_ROW_OF_C(1)
	INIT_ROW_OF_C(2)
	INIT_ROW_OF_C(3)
#undef INIT_ROW_OF_C

	if (bk > 1) {
		BLASLONG k = 1;
		vector_float Ak[VEC_ROWS], Aknext[VEC_ROWS];
		vector_float Bk[VEC_COLS], Bknext[VEC_COLS];

		/*
		 * Note that in several places, we enforce an instruction
		 * sequence that we identified empirically by utilizing dummy
		 * asm statements.
		 */

		for (BLASLONG j = 0; j < VEC_COLS; j++)
			Bk[j] = vec_splats(B[j + k * COLS]);
		asm("");

		for (BLASLONG i = 0; i < VEC_ROWS; i++)
			Ak[i] = vec_load_hinted(A + i * VLEN_FLOATS + k * ROWS);

		for (; k < (bk - 2); k += 2) {
			/*
			 * Load inputs for (k+1) into registers.
			 * Loading from B first is advantageous.
			 */
			for (BLASLONG j = 0; j < VEC_COLS; j++)
				Bknext[j] = vec_splats(B[j + (k + 1) * COLS]);
			asm("");
			for (BLASLONG i = 0; i < VEC_ROWS; i++)
				Aknext[i] = vec_load_hinted(A + i * VLEN_FLOATS +
						(k + 1) * ROWS);

			/*
			 * To achieve better instruction-level parallelism,
			 * make sure to first load input data for (k+1) before
			 * initiating compute for k. We enforce that ordering
			 * with a pseudo asm statement.
			 * Note that we need to massage this particular "barrier"
			 * depending on the gcc version.
			 */
#if __GNUC__ > 7 || defined(__clang__)
#define BARRIER_READ_BEFORE_COMPUTE(SUFFIX)                                    \
    do {                                                                       \
	asm(""                                                                 \
	    : "+v"(C_0_0), "+v"(C_0_1), "+v"(C_0_2), "+v"(C_0_3), "+v"(C_1_0), \
	      "+v"(C_1_1), "+v"(C_1_2), "+v"(C_1_3)                            \
	    : "v"(B##SUFFIX[0]), "v"(B##SUFFIX[1]), "v"(B##SUFFIX[2]),         \
	      "v"(B##SUFFIX[3]), "v"(A##SUFFIX[0]), "v"(A##SUFFIX[1]),         \
	      "v"(A##SUFFIX[2]), "v"(A##SUFFIX[3]));                           \
	asm(""                                                                 \
	    : "+v"(C_2_0), "+v"(C_2_1), "+v"(C_2_2), "+v"(C_2_3), "+v"(C_3_0), \
	      "+v"(C_3_1), "+v"(C_3_2), "+v"(C_3_3)                            \
	    : "v"(B##SUFFIX[0]), "v"(B##SUFFIX[1]), "v"(B##SUFFIX[2]),         \
	      "v"(B##SUFFIX[3]), "v"(A##SUFFIX[0]), "v"(A##SUFFIX[1]),         \
	      "v"(A##SUFFIX[2]), "v"(A##SUFFIX[3]));                           \
    } while (0)
#else // __GNUC__ <= 7
#define BARRIER_READ_BEFORE_COMPUTE(SUFFIX) \
    do {                                    \
	asm("");                            \
    } while (0)
#endif

			BARRIER_READ_BEFORE_COMPUTE(knext);

			/* Compute for (k) */
			C_0_0 += Ak[0] * Bk[0];
			C_1_0 += Ak[1] * Bk[0];
			C_2_0 += Ak[2] * Bk[0];
			C_3_0 += Ak[3] * Bk[0];

			C_0_1 += Ak[0] * Bk[1];
			C_1_1 += Ak[1] * Bk[1];
			C_2_1 += Ak[2] * Bk[1];
			C_3_1 += Ak[3] * Bk[1];

			C_0_2 += Ak[0] * Bk[2];
			C_1_2 += Ak[1] * Bk[2];
			C_2_2 += Ak[2] * Bk[2];
			C_3_2 += Ak[3] * Bk[2];

			C_0_3 += Ak[0] * Bk[3];
			C_1_3 += Ak[1] * Bk[3];
			C_2_3 += Ak[2] * Bk[3];
			C_3_3 += Ak[3] * Bk[3];

			asm("");

			/*
			 * Load inputs for (k+2) into registers.
			 * First load from B.
			 */
			for (BLASLONG j = 0; j < VEC_COLS; j++)
				Bk[j] = vec_splats(B[j + (k + 2) * COLS]);
			asm("");
			for (BLASLONG i = 0; i < VEC_ROWS; i++)
				Ak[i] = vec_load_hinted(A + i * VLEN_FLOATS + (k + 2) * ROWS);

			/*
			 * As above, make sure to first schedule the loads for (k+2)
			 * before compute for (k+1).
			 */
			BARRIER_READ_BEFORE_COMPUTE(k);

			/* Compute on (k+1) */
			C_0_0 += Aknext[0] * Bknext[0];
			C_1_0 += Aknext[1] * Bknext[0];
			C_2_0 += Aknext[2] * Bknext[0];
			C_3_0 += Aknext[3] * Bknext[0];

			C_0_1 += Aknext[0] * Bknext[1];
			C_1_1 += Aknext[1] * Bknext[1];
			C_2_1 += Aknext[2] * Bknext[1];
			C_3_1 += Aknext[3] * Bknext[1];

			C_0_2 += Aknext[0] * Bknext[2];
			C_1_2 += Aknext[1] * Bknext[2];
			C_2_2 += Aknext[2] * Bknext[2];
			C_3_2 += Aknext[3] * Bknext[2];

			C_0_3 += Aknext[0] * Bknext[3];
			C_1_3 += Aknext[1] * Bknext[3];
			C_2_3 += Aknext[2] * Bknext[3];
			C_3_3 += Aknext[3] * Bknext[3];
		}

		/* Wrapup remaining k's */
		for (; k < bk; k++) {
			vector_float Ak;

#define COMPUTE_WRAPUP_ROW(ROW)                             \
    Ak = vec_load_hinted(A + ROW * VLEN_FLOATS + k * ROWS); \
    C_##ROW##_0 += Ak * B[0 + k * COLS];                    \
    C_##ROW##_1 += Ak * B[1 + k * COLS];                    \
    C_##ROW##_2 += Ak * B[2 + k * COLS];                    \
    C_##ROW##_3 += Ak * B[3 + k * COLS];

			COMPUTE_WRAPUP_ROW(0)
			COMPUTE_WRAPUP_ROW(1)
			COMPUTE_WRAPUP_ROW(2)
			COMPUTE_WRAPUP_ROW(3)
#undef COMPUTE_WRAPUP_ROW
		}
	}

	/*
	 * Unpack row-block of C_aux into outer C_i, multiply by
	 * alpha and add up (or assign for TRMM).
	 */
#define WRITE_BACK_C(ROW, COL)                                   \
    do {                                                         \
	vector_float *Cij =                                      \
	    (vector_float *)(C + ROW * VLEN_FLOATS + COL * ldc); \
	if (trmm) {                                              \
	    *Cij = alpha * C_##ROW##_##COL;                   \
	} else {                                                 \
	    *Cij += alpha * C_##ROW##_##COL;                  \
	}                                                        \
    } while (0)

	WRITE_BACK_C(0, 0); WRITE_BACK_C(0, 1); WRITE_BACK_C(0, 2); WRITE_BACK_C(0, 3);
	WRITE_BACK_C(1, 0); WRITE_BACK_C(1, 1); WRITE_BACK_C(1, 2); WRITE_BACK_C(1, 3);
	WRITE_BACK_C(2, 0); WRITE_BACK_C(2, 1); WRITE_BACK_C(2, 2); WRITE_BACK_C(2, 3);
	WRITE_BACK_C(3, 0); WRITE_BACK_C(3, 1); WRITE_BACK_C(3, 2); WRITE_BACK_C(3, 3);
#undef WRITE_BACK_C

#undef ROWS
#undef VEC_ROWS
#undef COLS
#undef VEC_COLS
#undef BARRIER_READ_BEFORE_COMPUTE
}

/**
 * Handle calculation for row blocks in C_i of any size by dispatching into
 * macro-defined (inline) functions or by deferring to a simple generic
 * implementation. Note that the compiler can remove this awkward-looking
 * dispatching code while inlineing.
 *
 * @param[in]	m	Number of rows in block C_i.
 * @param[in]	n	Number of columns in block C_i.
 * @param[in]	first_row Index of first row of the block C_i (relative to C).
 * @param[in]	A	Pointer to input matrix A (note: all of it).
 * @param[in]	k	Number of columns in A and rows in B.
 * @param[in]	B	Pointer to current column block (panel) of input matrix B.
 * @param[inout] C	Pointer to current column block (panel) of output matrix C.
 * @param[in]	ldc	Offset between elements in adjacent columns in C.
 * @param[in]	alpha	Scalar factor.
 * @param[in]	offset  Number of columns of A and rows of B to skip (for triangular matrices).
 * @param[in]	off	Running offset for handling triangular matrices.
 */
static inline void GEBP_block(BLASLONG m, BLASLONG n,
		       BLASLONG first_row,
		       const FLOAT * restrict A, BLASLONG k,
		       const FLOAT * restrict B,
		       FLOAT *restrict C, BLASLONG ldc,
		       FLOAT alpha,
		       BLASLONG offset, BLASLONG off)
{
	if (trmm && left)
		off = offset + first_row;

	A += first_row * k;
	C += first_row;

	if (trmm) {
		if (backwards) {
			A += off * m;
			B += off * n;
			k -= off;
		} else {
			if (left) {
				k = off + m;
			} else {
				k = off + n;
			}
		}
	}

	/* Dispatch into the implementation for each block size: */

#define BLOCK(bm, bn)                                           \
	if (m == bm && n == bn) {                               \
		GEBP_block_##bm##_##bn(A, k, B, C, ldc, alpha); \
		return;                                         \
	}

#if UNROLL_M == 16
	BLOCK(16, 4); BLOCK(16, 2); BLOCK(16, 1);
#endif
#if UNROLL_N == 8
	BLOCK(8, 8); BLOCK(4, 8);
#endif
	BLOCK(8, 4); BLOCK(8, 2); BLOCK(8, 1);
	BLOCK(4, 4); BLOCK(4, 2); BLOCK(4, 1);

	BLOCK(2, 4); BLOCK(2, 2); BLOCK(2, 1);

	BLOCK(1, 4); BLOCK(1, 2); BLOCK(1, 1);

#undef BLOCK
}

/**
 * Handle a column block (panel) of C and B while calculating C += alpha(A * B).
 *
 * @param[in]	num_cols	Number of columns in the block (in C and B).
 * @param[in]	first_col	First column of the current block (in C and B).
 * @param[in]	A	Pointer to input matrix A.
 * @param[in]	bk	Number of columns in A and rows in B.
 * @param[in]	B	Pointer to input matrix B (note: all of it).
 * @param[in]	bm	Number of rows in C and A.
 * @param[inout] C	Pointer to output matrix C (note: all of it).
 * @param[in]	ldc	Offset between elements in adjacent columns in C.
 * @param[in]	alpha	Scalar factor.
 * @param[in]	offset	Number of columns of A and rows of B to skip (for triangular matrices).
 */
static inline void GEBP_column_block(BLASLONG num_cols, BLASLONG first_col,
			const FLOAT *restrict A, BLASLONG bk,
			const FLOAT *restrict B, BLASLONG bm,
			FLOAT *restrict C, BLASLONG ldc,
			FLOAT alpha,
			BLASLONG const offset) {

	FLOAT *restrict C_i = C + first_col * ldc;
	/*
	 * B is in column-order with n_r packed row elements, which does
	 * not matter -- we always move in full such blocks of
	 * column*pack
	 */
	const FLOAT *restrict B_i = B + first_col * bk;

	BLASLONG off = 0;
	if (trmm) {
		if (left) {
			off = offset;
		} else {
			off = -offset + first_col;
		}
	}

	/*
	 * Calculate C_aux := A * B_j
	 * then unpack C_i += alpha * C_aux.
	 *
	 * For that purpose, further partition C_aux and A into blocks
	 * of m_r (unroll_m) rows, or powers-of-2 if smaller.
	 */
	BLASLONG row = 0;
	for (BLASLONG block_size = unroll_m; block_size > 0; block_size /= 2)
		for (; bm - row >= block_size; row += block_size)
			GEBP_block(block_size, num_cols, row, A, bk, B_i, C_i,
				   ldc, alpha, offset, off);
}

/**
 * Inner kernel for matrix-matrix multiplication. C += alpha (A * B)
 * where C is an m-by-n matrix, A is m-by-k and B is k-by-n. Note that A, B, and
 * C are pointers to submatrices of the actual matrices.
 *
 * For triangular matrix multiplication, calculate B := alpha (A * B) where A
 * or B can be triangular (in case of A, the macro LEFT will be defined).
 *
 * @param[in]	bm	Number of rows in C and A.
 * @param[in]	bn	Number of columns in C and B.
 * @param[in]	bk	Number of columns in A and rows in B.
 * @param[in]	alpha	Scalar factor.
 * @param[in]	ba	Pointer to input matrix A.
 * @param[in]	bb	Pointer to input matrix B.
 * @param[inout] C	Pointer to output matrix C.
 * @param[in]	ldc	Offset between elements in adjacent columns in C.
 * @param[in]	offset	Number of columns of A and rows of B to skip (for triangular matrices).
 * @returns 0 on success.
 */
int CNAME(BLASLONG bm, BLASLONG bn, BLASLONG bk, FLOAT alpha,
	  FLOAT *restrict ba, FLOAT *restrict bb,
	  FLOAT *restrict C, BLASLONG ldc
#ifdef TRMMKERNEL
	  , BLASLONG offset
#endif
	  )
{
	if ( (bm == 0) || (bn == 0) || (bk == 0) || (alpha == ZERO))
		return 0;

	/*
	 * interface code allocates buffers for ba and bb at page
	 * granularity (i.e., using mmap(MAP_ANONYMOUS), so enable the compiler
	 * to make use of the fact in vector load operations.
	 */
	ba = __builtin_assume_aligned(ba, 16);
	bb = __builtin_assume_aligned(bb, 16);

	/*
	 * Use offset and off even when compiled as SGEMMKERNEL to simplify
	 * function signatures and function calls.
	 */
#ifndef TRMMKERNEL
	BLASLONG const offset = 0;
#endif

	/*
	 * Partition B and C into blocks of n_r (unroll_n) columns, called B_i
	 * and C_i. For each partition, calculate C_i += alpha * (A * B_j).
	 *
	 * For remaining columns that do not fill up a block of n_r, iteratively
	 * use smaller block sizes of powers of 2.
	 */
	BLASLONG col = 0;
	for (BLASLONG block_size = unroll_n; block_size > 0; block_size /= 2)
		for (; bn - col >= block_size; col += block_size)
			GEBP_column_block(block_size, col, ba, bk, bb, bm, C, ldc, alpha, offset);

   return 0;
}