IMPLEMENTING HIGH-PERFORMANCE COMPLEX MATRIX

MULTIPLICATION VIA THE 1M METHOD

FIELD G. VAN ZEE∗

Abstract. Almost all efforts to optimize high-performance matrix-matrix multiplication have

been focused on the case where matrices contain real elements. The community’s collective assump-
tion appears to have been that the techniques and methods developed for the real domain carry
over directly to the complex domain. As a result, implementors have mostly overlooked a class of
methods that compute complex matrix multiplication using only real matrix products. This is the
second in a series of articles that investigate these so-called induced methods. In the previous arti-
cle, we found that algorithms based on the more generally applicable of the two methods—the 4m
method—lead to implementations that, for various reasons, often underperform their real domain
counterparts. To overcome these limitations, we derive a superior 1m method for expressing complex
matrix multiplication, one which addresses virtually all of the shortcomings inherent in 4m. We show
that the method is actually a special case of a larger family of algorithms based on a 2m method,
which is generally well-suited for storage formats that organize real and imaginary parts into separate
matrices. Implementations are developed within the BLIS framework, and testing on a recent Intel
microarchitecture confirms that the 1m method yields performance that is competitive with solutions
based on conventionally implemented complex kernels.

Key words. high-performance, complex, matrix, multiplication, microkernel, kernel, BLAS,

BLIS, 1m, 2m, 4m, induced, linear algebra, DLA

AMS subject classifications. 65Y04

1. Introduction. Over the last several decades, matrix multiplication research

has resulted in methods and implementations that primarily target the real domain.
Recent trends in implementation efforts have condensed virtually all matrix product
computation into relatively small kernels—building blocks of highly optimized code
(typically written in assembly language) upon which more generalized functionality
is constructed via various levels of nested loops [23, 6, 4, 22, 3]. Because most effort
is focused on the real domain, the complex domain is either left as an unimplemented
afterthought—perhaps because the product is merely a proof-of-concept or proto-
type [6], or because the project primarily targets applications and uses cases that
require only real computation [3]—or it is implemented in a manner that mimics the
real domain down to the level of the assembly kernel [23, 4, 5].1 Most modern mi-
croarchitectures lack machine instructions for directly computing complex arithmetic
on complex numbers, and so when the effort to implement these kernels is undertaken,
kernel developers encounter additional programming challenges that do not manifest
in the real domain. Specifically, these kernel developers must explicitly orchestrate
computation on the real and imaginary components in order to implement multi-
plication and addition on complex scalars, and they must do so in terms of vector
instructions to ensure high performance is achievable.
Pushing the nuances and complexities of complex arithmetic down to the level of
the kernel allows the higher-level loop infrastructure within the matrix multiplication
∗ Oden Institute for Computational Engineering & Sciences, The University of Texas at Austin,

Austin, TX (field@[Link])
Funding: This research was partially sponsored by grants from Intel Corporation and the
National Science Foundation (Award ACI-1550493). Any opinions, findings and conclusions or
recommendations expressed in this material are those of the author(s) and do not necessarily reflect
the views of the National Science Foundation (NSF).
1 Because they exhibit slightly less favorable numerical properties, we exclude Strassen-like efforts

from this characterization.

1
2 FIELD G. VAN ZEE

5th loop around micro-kernel

nC nC

Cj += A Bj

4th loop around micro-kernel

kC Bp
Cj Ap
+=
kC ~
Pack Bp → Bp
3rd loop around micro-kernel
~
Ci mC mC Ai Bp
+= ~
Pack Ai → Ai

2nd loop around micro-kernel

nR ~ nR ~
Ci Ai Bp
mR

+=
kC
1st loop around micro-kernel
nR

mR
+= kC

micro-kernel
main memory
1
L3 cache
L2 cache
L1 cache += 1
registers

Fig. 1.1. An illustration of the algorithm for computing high-performance matrix multiplication,
taken from [21], which expresses computation in terms of a so-called “block-panel” subproblem.

to remain largely the same as its real domain counterpart. (See Figure 1.1.) However,
this approach also doubles the number of assembly kernels that must be written in
order to fully support computation in either domain (real or complex) for the desired
floating-point precisions. And while computation in the complex domain may not be
of interest to all developers, it is absolutely essential for many fields and applications
in part because of complex numbers’ unique ability to encode both the phase and
magnitude of a wave. Thus, the maintainers of general-purpose matrix libraries—such
as those that export the Basic Linear Algebra Subprograms (BLAS) [2]—are typically
compelled by their diverse user bases to support complex matrix multiplication despite
the implementation and maintenance costs it may impose.
Life would be simpler for matrix library developers if complex matrix multipli-
cation was not necessary. Of course, complex matrix multiplication will always be
necessary. But what if complex matrix multiplication kernels were found to be un-
necessary? To certain actors, particularly matrix library developers, such a finding
would carry non-trivial consequence.
The predecessor to the current article investigates whether (and to what degree
of effectiveness) real domain matrix multiplication kernels can be repurposed and
leveraged toward the implementation of complex matrix multiplication [21]. In that
article, the authors develop a new class of algorithms that implement these so-called
“induced methods” for matrix multiplication in the complex domain. Instead of re-
lying on an assembly-coded complex kernel, as a conventional implementation would,
these algorithms express complex matrix multiplication only in terms of real domain
 COMPLEX MATRIX MULTIPLICATION VIA THE 1M METHOD 3

primitives.2 The authors first investigate the 4m method, which is merely the matrix
analogue of the classic definition of complex scalar multiplication and addition. They
show how this method’s implementation is facilitated by reordering real and imaginary
elements within the internal storage format used when making temporary contiguous
“packed” copies of the current matrix blocks. The second3 approach, the 3m method,
benefits from requiring only three (instead of four) matrix products at the cost of an
additional three accumulations. This Strassen-like method, which is capable of faster
runtimes but exhibits slightly lower numerical accuracy, is implemented by employing
a special data format similar to that of 4m, except that it also serves as temporary
workspace for storing the results of intermediate computations.4
We consider the current article a companion and follow-up to that previous
work [21]. Here, we will once again consider a new method for emulating complex
matrix multiplication using only real domain building blocks, and we will once again
show that a rearrangement of the real and imaginary elements within the packed ma-
trices is key to facilitating an induced method capable of delivering performance that
rivals solutions based on the conventional approach.
1.1. Audience. It is also worth emphasizing that this article’s primary audience
consists of dense linear algebra (DLA) library implemenators and related experts.
While this audience is rather narrow, the contributions we present have the potential
for wide indirect impact since the community relies heavily on this relatively small
group of implementors to develop and maintain DLA kernels for the benefit of the
entire high-performance computing (HPC) ecosystem.
1.2. Contributions. This article makes the following contributions:
• It introduces a new induced method—the 1m method—that replaces each
complex matrix multiplication with only a single real matrix multiplication.
We introduce two algorithms and, as with the previous article, analyze issues
germane to their high-performance implementations, including workspace,
packing formats, cache behavior, multithreadability, and programming effort.
A detailed review shows how 1m avoids virtually all of the challenges observed
of the 4m method.
• It promotes code reuse and portability by continuing the previous article’s
focus on solutions which may be cast in terms of real matrix multiplication
kernels. Such solutions have clear implications for developer productivity, as
they allow kernel authors to focus their efforts on fewer and simpler kernels.
• It builds on the theme of the BLIS framework as a productivity multiplier [22],
further demonstrating how complex matrix multiplication may be imple-
mented with relatively minor modifications to the source code and in such a
way that results in immediate instantiation of complex implementations for
all level-3 BLAS-like operations.
• It demonstrates performance of 1m implementations that is not only superior
to the previous effort based on the 4m method but also competitive with

2 In [21], the authors use the term “primitive” to refer to a functional abstraction that implements

a single real matrix multiplication. Such primitives are often not general purpose, and may come
with significant prerequisites to facilitate their use.
3 While this was the first method investigated, the authors of [21] actually choose to present

them in reverse order.

4 Others have exploited the careful design of packing and computational primitives in an effort to

improve performance, including in the context of Strassen’s algorithm [8, 10, 11, 12], the computation
of the K-Nearest Neighbors [24], and tensor contraction [9].
4 FIELD G. VAN ZEE

solutions based on complex matrix kernels.

• It serves as a reference guide to the 1m implementations for complex matrix
multiplication found within the BLIS framework, which is available to the
community under an open source software license.5
We believe that these contributions are consequential because the 1m method effec-
tively obviates the previous state-of-the-art established via the 4m method. Further-
more, we believe the thorough treatment of induced methods encompassed by the
present article and its predecessor will have lasting archival as well as pedagogical
value, to say nothing of the potential impact on developer productivity.
1.3. Notation. In this article, we continue the notation established in [21].
Specifically, we use uppercase Roman letters (e.g. A, B, and C) to refer to ma-
trices, lowercase Roman letters (e.g. x, y, and z) to refer to vectors, and lowercase
Greek letters (e.g. χ, ψ, and ζ) to refer to scalars. Subscripts are used typically

to
denote sub-matrices within a larger matrix (e.g. A = A0 A1 · · · An−1 ) or
scalars within a larger matrix or vector.
We make extensive use of superscripts to denote the real and imaginary compo-
nents of a scalar, vector, or (sub-)matrix. For example, αr , αi ∈ R denote the real
and imaginary parts, respectively, of a scalar α ∈ C. Similarly, Ar and Ai refer to the
real and imaginary parts of a complex matrix A, where Ar and Ai are real matrices
with dimensions identical to A. Note that while this notation for real, imaginary, and
complex matrices encodes information about content and origin, it does not encode
how the matrices are actually stored. We will explicitly address storage details as
implementation issues are discussed.
Also, at times we find it useful to refer to the real and imaginary elements of
a complex object indistinguishably as fundamental elements (or f.e.). We also ab-
breviate floating-point operations as “flops” and memory operations as “memops”.
We define the former to be a Multiply or Add (or Subtract) operation whose
operands are fundamental elements and the latter to be a load or store operation on
a single fundamental element. These definitions allow for a consistent accounting of
complex computation relative to the real domain.
We also discuss cache and register blocksizes that are key features of the matrix
multiplication algorithm discussed elsewhere [22, 20, 21]. Unless otherwise noted,
blocksizes nC , mC , kC , mR , and nR refer to those appropriate for computation in the
real domain. Complex domain blocksizes will be denoted with a superscript z .
This article discusses and references several hypothetical algorithms and func-
tions. Unless otherwise noted, a call to function func that implements C := C + AB
appears as [ C ] := func( A, B, C ). We will also reference functions that access
properties of matrices. For example, m(A) and n(A) would return the m and n di-
mensions of a matrix A, while rs(B) and cs(B) would return the row and column
strides of B.
2. Background and review.
2.1. Motivation. In [21], the authors list three primary motivating factors be-
hind their effort to seek out methods for inducing complex matrix multiplication via
real domain kernels:
• Productivity. By inducing complex matrix multiplication from real domain
kernels, the number of kernels that must be supported would be halved.
5 The BLIS framework is available under the so-called “new” or “modified” or “3-clause” BSD
license.
 COMPLEX MATRIX MULTIPLICATION VIA THE 1M METHOD 5

This allows the DLA library developers to focus on a smaller and simpler
set of real domain kernels. This benefit would manifest most obviously when
instantiating BLAS-like functionality on new hardware [20].
• Portability. Induced methods avoid dependence on complex domain kernels
because they encode the idea of complex matrix product at a higher level.
This would naturally allow us to encode such methods portably within a
framework such as BLIS [22]. Once integrated into the framework, developers
and users would benefit from the immediate availability of complex matrix
multiplication implementations whenever real matrix kernels were present.
• Performance. Implementations of complex matrix multiplication that rely
on real domain kernels would likely inherit the high-performance properties
of those kernels. Any improvement to the real kernels would benefit both real
and complex domains.
Thus, it is clear that finding a suitable induced method would carry significant benefit
to DLA library and kernel developers.
2.2. The 3m and 4m methods. The authors of [21] investigated two general
ways of inducing complex matrix multiplication: the 3m method and the 4m method.
These methods are then contrasted to the conventional approach, whereby a blocked
matrix multiplication algorithm is executed with a complex domain kernel—one that
implements complex arithmetic at the scalar level, in assembly language.
The 4m method begins with the classic definition of complex scalar multiplication
and addition in terms of real and imaginary components of α, β, γ ∈ C:

γ r := γ r + αr β r − αi β i
(2.1) γ i := γ i + αi β r + αr β i

We then observe that we can apply such a definition to complex matrices A ∈ Cm×k ,
B ∈ Ck×n , and C ∈ Cm×n , provided that we can reference the real and imaginary
parts as logically separate submatrices:

C r := C r + Ar B r − Ai B i
(2.2) C i := C i + Ai B r + Ar B i

This definition expresses a complex matrix multiplication in terms of four matrix

products (hence the name 4m) and four matrix accumulations (i.e., additions or sub-
tractions).
The 3m method relies on a Strassen-like algebraic equivalent of Eq. 2.2:

C r := C r + Ar B r − Ai B i
C i := C i + Ar + Ai B r + B i − Ar B r − Ai B i

This re-expression reduces the number of matrix products to three at the expense of
increasing the number of accumulations from four to seven. However, when the cost
of a matrix product greatly exceeds that of an accumulation, this trade-off can result
in a net reduction in computational runtime.
The authors of [21] observe that both methods may be applied to any particular
level of a blocked matrix multiplication algorithm, resulting in several algorithms,
each exhibiting somewhat different properties. The blocked algorithm used in that
article is shown in Figure 1.1 and revisited in Section 2.4 of the present article.
6 FIELD G. VAN ZEE

Algorithms that implement the 3m method were found to yield “effective flops
per second” performance that not only exceeded that of 4m, but also approached or
exceeded the theoretical peak rate of the hardware.6 Unfortunately, these compelling
results come at a cost: the numerical properties of implementations based on 3m are
slightly less robust than that of algorithms based on the conventional approach or 4m.
And although the author of [7] found that 3m was stable enough for most practical
purposes, many applications simply will not be willing to stray from the numerical
expectations implicit in conventional matrix multiplication. Thus, going forward, we
will turn our attention away from 3m and instead focus on the 4m as the standard
reference method against which we will compare.
2.3. Previous findings. For the reader’s convenience, we will now summarize
the key findings, observations, and other highlights from the previous article regarding
algorithms and implementations based on the 4m method [21].
• Since all algorithms in the 4m family execute the same number of flops, the
algorithms’ relative performance depends entirely on (1) the number of mem-
ops executed and (2) the level of cache from which fundamental elements of
matrices A and B (or rather, f.e. of the packed copies of these matrices, Ãi
and B̃p ) are reused7 . The number of memops is affected only by a halving of
certain cache blocksize needed in order to leave cache footprints of Ãi and B̃p
unchanged. The level from which f.e. are reused is determined by the level
of the matrix multiplication algorithm to which the 4m method was originally
applied. The lower the 4m method is applied, the higher the efficiency of data
reuse from and movement through the cache hierarchy.
• The lowest-level application, algorithm 4m 1a, efficiently moves f.e. of A,
B, and C from main memory to the L1 cache only once per rank-kC update,
with virtually no excess movement due to incidental cache line proximity,
and reuses f.e. from the L1 cache. It relies on a relatively simple packing
format in which complex micro-panels are stored with real and imaginary f.e.
separated into two consecutive real micro-panels, each with identical register
blocksize and k dimensions. Algorithm 4m 1a requires negligible workspace
(limited to the storage capacity of the vector register set), is well-suited for
multithreading, and is minimally disruptive to the encoding within the BLIS
framework. And while algorithm 4m 1b—a slightly higher-level application—
very narrowly outperforms 4m 1a by trading away the most optimal cache
reuse behavior for an unreduced kC cache blocksize, 4m 1a proves to be more
versatile and can be extended relatively easily to all other level-3 operations.
• The conventional assembly-based approach to complex matrix multiplication
can be viewed as a special case of 4m in which f.e. are reused from registers
rather than cache. In this way, a conventional implementation embodies the
lowest-level application of 4m possible, in which the method is applied to
individual scalars (and, typically, then optimally encoded via vector instruc-
tions).
• The way complex numbers are stored has a significant effect on performance.
6 Note that 3m and other Strassen-like algorithms are able to exceed the hardware’s theoretical

peak performance when measured in effective flops per second: that is, the 3m implementation’s
wall clock time—now shorter because of avoided matrix products—divided into the flop count of a
conventional algorithm.
7 Here, the term “reuse” refers to the reuse of f.e. that corresponds to the recurrence of Ar ,

Ai , B r , and B i in Eq. 2.2, not the reuse of whole (complex) elements that naturally occurs in the
execution of the matrix multiplication algorithm in Figure 1.1.
 COMPLEX MATRIX MULTIPLICATION VIA THE 1M METHOD 7

Interleaved pair-wise storage of real and imaginary values naturally favors

implementations that reuse f.e. from vector registers, as is common with
conventional implementations.8 However, this storage is awkward for algo-
rithms based on 4m (and 3m), largely because it prevents the use of vector
instructions for loading and storing f.e. of C r and C i910 . The 4m 1a algo-
rithm already suffers from a quadrupling11 of the number of memops on C, in
addition to being forced to access these f.e. in a non-contiguous manner. If,
however, applications stored complex matrices with real and imaginary parts
separated, the penalty paid by 4m (and 3m) would be partially mitigated.
• While observed performance of low-level applications of 4m is decent, far
exceeding an unoptimized reference implementation, it not only falls short of a
comparable conventional solution based on a complex kernel, it appears to fall
short of its real domain “benchmark”—that is, the performance of a similar
problem size in the real domain computed by an optimized algorithm using the
same real domain kernel. This level of performance may be disappointing for
some, even if it achieves 90-95% of what is possible with a complex kernel. The
authors conclude that its somewhat attenuated performance would relegate
4m, in practice, to serving mostly as a placeholder, to be used when complex
kernels have not yet been written, rather than a competitive replacement that
makes complex kernels unnecessary over a longer time horizon.
2.4. Revisiting the matrix multiplication algorithm. In this section, we
review a common algorithm for high-performance matrix multiplication on conven-
tional microprocessor architectures. This algorithm was first reported on in [4] and
further refined in [22]. Figure 1.1 illustrates the key features of this algorithm.
The current state-of-the-art formulation of the matrix multiplication algorithm
consists of six loops, the last of which resides within a micro-kernel that is typically
highly optimized for the target hardware. These loops partition the matrix operands
using carefully chosen cache (nC , kC , and mC ) and register (mR and nR ) blocksizes
that result in submatrices residing favorably at various levels of the cache hierarchy
so as to allow data to be reused many times. In addition, submatrices of A and B
are copied (“packed”) to temporary workspace matrices (Ãi and B̃p , respectively)
in such a way that allows the micro-kernel to subsequently access matrix elements
contiguously in memory, which improves cache and TLB performance. The cost of this
packing is amortized over enough computation that its impact on overall performance
is negligible for all but the smallest problems. At the lowest level, within the micro-
kernel loop, an mR × 1 micro-column and a 1 × nR micro-row are loaded from the
current micro-panels of Ãi and B̃p , respectively, so that the outer product of these
vectors may be computed to update the corresponding mR × nR submatrix, or micro-
tile, of C. The individual floating-point operations that constitute these tiny rank-1
updates are oftentimes executed via vector instructions (if the architecture supports
them) in order to maximize utilization of the floating-point unit(s).

8 The advantages of interleaving data in advance of computation via vector instructions is also

discussed by the authors of [1] in the context of performing “batched” level-3 operations.
9 The traditional pair-wise storage is also awkward for 4m algorithms during the packing of data

from A and B, but this effect is not nearly as dramatic.

10 This awkwardness will persist unless and until hardware architectures begin providing vector

instructions for efficiently loading and storing non-contiguous elements in memory.

11 A factor of two comes from the fact that, as shown in Eq. 2.2, 4m touches C r and C i twice

each, while another factor of two comes from the cache blocksize scaling required on kC in order to
maintain the cache footprints of micro-panels of Ãi and B̃p .
8 FIELD G. VAN ZEE

The algorithm captured by Figure 1.1 forms the basis for all level-3 implementa-
tions found in the BLIS framework (as of this writing). This algorithm is based on a
so-called block-panel matrix multiplication.12 The register (mR , nR ) and cache (mC ,
kC , nC ) blocksizes labeled in the algorithmic diagram are typically chosen by the
kernel developer as a function of hardware characteristics, such as the vector register
set, cache sizes, and cache associativity. The authors of [16] present an analytical
model for identifying suitable (if not optimal) values for these blocksizes.
3. 1m method. The primary motivation for seeking a better induced method
comes from the observation that 4m inherently must update real and imaginary f.e.
of C: (1) in separate steps, and may not use vector instructions to do so (due to the
traditional pair-wise storage format); and (2) twice as frequently, in the case of 4m 1a,
due to the algorithm’s half-of-optimal cache blocksize kC . As reviewed in Section 2.3,
this imposes a significant drag on performance. If there existed an induced method
that could update real and imaginary elements in one step, it may conveniently avoid
both issues.
3.1. Derivation. Consider the classic definition of complex scalar multiplication
and accumulation, shown in Eq. 2.1, refactored and expressed in terms of matrix and
vector notation:
 r  r
α −αi
 r 
γ β
(3.1) +=
γi αi αr βi

Here, we have a singleton complex matrix multiplication problem that can naturally
be expressed as a tiny real matrix multiplication where m = k = 2 and n = 1.
Let us assume we implement this very small matrix multiplication according to the
high-performance algorithm discussed in Section 2.4.
From this, we make the following key observation: If we pack α to Ãi in such a
way that duplicates αr and αi to the second column of the micro-panel (while also
swapping the placement of the duplicates and negating the duplicated αi ), and if we
pack β to B̃p such that β i is stored to the second row of the micro-panel (which,
granted, only has one column), then a real domain gemm micro-kernel executed on
those micro-panels will compute the correct result in the complex domain and do so
with a single invocation of that micro-kernel.
Thus, Eq. 3.1 serves as a packing template that hints at how the data must be
stored. Furthermore, this template can be generalized. We augment α, β, γ with
conventional row and column indices to denote the complex elements of matrices A,
B, and C, respectively. Also, let us apply the Eq. 3.1 to the special case of m = 3,
n = 4, and k = 2 to better observe the general pattern.
 r r r r
  r i r i

γ00 γ01 γ02 γ03 α00 −α00 α01 −α01
i i i i  i r i r 
 r r r r

 γ00
 r γ01 γ02 γ03  α00 α00 α01 α01 β00 β01 β02 β03
r r

r 
 
 γ10 γ11 γ12 γ13  r −α10 i r
−α11i  i i i i 
(3.2)  i i i i
 +=  α10 i r
α11
i r
  β00
r
β01
r
β02
r
β03
r 

 γ10 γ11 γ12 γ13   α10 α10 α11 α11   β10 β11 β12 β13
 r r r

r 
 r i r

i  i i i i
 γ20 γ21 γ22 γ23  α20 −α20 α21 −α21 β10 β11 β12 β13
i i i i i r i r
γ20 γ21 γ22 γ23 α20 α20 α21 α21

From this, we can make the following observations:

12 This terminology describes the shape of the typical problem computed by the macro-kernel, i.e.

the second loop around the micro-kernel. An alternative algorithm that casts its largest cache-bound
subproblem in terms of panel-block matrix multiplication is discussed in [19].
 COMPLEX MATRIX MULTIPLICATION VIA THE 1M METHOD 9

• The complex matrix multiplication C := C + AB with m = 3, n = 4, and

k = 2 becomes a real matrix multiplication with m = 6, n = 4, and k = 4.
In other words, the m and k dimensions are doubled for the purposes of the
real gemm primitive.
• If the primitive is the real gemm micro-kernel, and we assume that matrices
A and B above represent column-stored and row-stored micro-panels from
Ãi and B̃p , respectively, and also that the dimensions are conformal to the
register blocksizes of this micro-kernel (i.e., m = mR and n = nR ) then the
micro-panels of Ãi are packed from a 12 mR × 12 kC submatrix of A, which, when
expanded in the special packing format, appears as the mR × kC micro-panel
that the real gemm micro-kernel expects.
• Similarly, the micro-panels of B̃p are packed from a 12 kC × nR submatrix of
B, which, when reordered into a second special packing format, appears as
the kC × nR micro-panel that the real gemm micro-kernel expects.
It is easy to see by inspection that the real matrix multiplication implied by
Eq. 3.2 induces the desired complex matrix multiplication. We will refer to the packing
format used on matrix A above as the 1e format, since the f.e. are “expanded”
(i.e., duplicated to the next column, with the duplicates swapped and the imaginary
duplicate negated). Similarly, we will refer to the packing format used on matrix B
above as the 1r format, since the f.e. are merely reordered (i.e., imaginary elements
moved to the next row).
3.2. Two variants. Notice that implicit in the 1m method suggested by Eq. 3.2
is the fact that matrix C is stored by columns. This assumption is important; when A
and B are packed according to the 1e and 1r formats, respectively, C must be stored
by columns in order to allow the real domain primitive (or micro-kernel) to correctly
update the individual real and imaginary f.e. of C with the corresponding f.e. from
the matrix product AB.
Suppose that we instead refactored and expressed Eq. 2.1 as follows:
 r i
r i

r i

β β
(3.3) γ γ += α α
−β i β r
This gives us a different template, one that implies different packing formats for
matrices A and B. Applying Eq. 3.3 to the special case of m = 4, n = 3, and k = 2,
yields:
(3.4)
 r i r i r i   r i r i  r i r i r i

γ00 γ00 γ01 γ01 γ02 γ02 α00 α00 α01 α01 β00 β00 β01 β01 β02 β02
 γr γi γr γi γr γi   αr αi αr αi   −β i β r −β i β r −β i β r 
 10 10 11 11 12 12 
 r i r i r i  +=  10 10 11 11  
 00 00 01 01 02 02 

 γ20 γ20 γ21 γ21 γ22 γ22  r
 α20 i r i  r i r i r i 
α20 α21 α21   β10 β10 β11 β11 β12 β12 
r i r i r i r i r i i r i r i r
γ30 γ30 γ31 γ31 γ32 γ32 α30 α30 α31 α31 −β10 β10 −β11 β11 −β12 β12
In this variant, we see that matrix B, not A, is stored according to the 1e format
(where columns become rows), while matrix A is stored according to 1r (where rows
become columns). Also, we can see that matrix C must be stored by rows in order to
allow the real gemm micro-kernel to correctly update its f.e. with the corresponding
values from the matrix product AB.
Henceforth, we will refer to the 1m variant exemplified in Eq. 3.2 as 1m c since
it is predicated on column storage of the output matrix C, and we will refer to the
variant depicted in Eq. 3.4 as 1m r since it assumes the output matrix is stored by
rows.
10 FIELD G. VAN ZEE

Table 3.1
1m complex domain blocksizes as a function of real domain blocksizes

Blocksizes, in terms of real domain

Variant values, required for . . .
z
kC mzC nzC mzR mzP nzR nzP
1 1 1
1m c 2 kC 2 mC nC 2 mR mP nR nP
1 1 1
1m r 2 kC mC 2 nC mR mP 2 nR nP

Note: Blocksizes mP and nP represent the so-called “packing dimensions” for the micro-panels of
Ãi and B̃p , respectively. These values are analogous to the leading dimensions of matrices stored
by columns or rows. In BLIS micro-kernels, typically mR = mP and nR = nP , but sometimes the
kernel author may find it useful for mR < mP or nR < nP .

3.3. Determining complex blocksizes. As we alluded in Section 3.1, the

appropriate blocksizes to use with 1m are a function of the real domain blocksizes.
This makes sense, since the idea is to fool the real gemm micro-kernel, and the
various loops for register and cache blocking around the micro-kernel, into thinking
that it is computing a real domain matrix multiplication. Which blocksizes must be
modified (halved) and which are used unchanged depends on the variant of 1m being
executed—or, more specifically, which matrix is packed according to the 1e format.
Table 3.1 summarizes the complex domain blocksizes prescribed for 1m c and
1m r as a function of the real domain values. This is somewhat analogous to the
blocksize scaling described in Tables II and III in [21]. However, in that article the
scaling was optional in the sense that different scaling factors would still work, albeit
perhaps with a performance penalty. Here, in the case of Table 3.1, some scaling
factors (namely, on mR or nR ) are required in order for the 1m algorithm to function
properly.
Those familiar with the matrix multiplication algorithm implemented by the BLIS
framework, as depicted in Figure 1.1, may be unfamiliar with mP and nP , the so-
called packing dimensions. These values are, effectively, the leading dimensions of the
micro-panels. For most architectures, mP and nP are almost always equal to mR and
nR , respectively. However, in some situations, it may be convenient (or necessary) to
use mR < mP or nR < nP . In any case, these packing dimensions are never scaled,
even when their corresponding register blocksizes are scaled to accommodate the 1e
format, because the halving that would otherwise be called for is cancelled out by the
doubling of f.e. that manifests in the 1e format.
3.4. Algorithms.
3.4.1. General algorithm. Before investigating 1m method algorithms, we will
first provide algorithms for computing real matrix multiplication to serve as a reference
for the reader. Specifically, we provide pseudo-code, targeting the real domain, for
the block-panel algorithm depicted in Figure 1.1. This algorithm is shown as rmmbp
in Figure 3.1.
3.4.2. 1m-specific algorithm. When applied to the block-panel algorithm de-
picted in Figure 1.1, 1m c and 1m r yield nearly identical algorithms whose differences
can be encoded within a few conditional statements within key parts of the high and
low levels of code. We will refer to these specific algorithms as 1m c bp and 1m r bp,
 COMPLEX MATRIX MULTIPLICATION VIA THE 1M METHOD 11

Algorithm: [ C ] := rmmbp( A, B, C )
for ( j = 0 : n − 1 : nC )
Identify Bj , Cj from B, C
for ( p = 0 : k − 1 : kC )
Identify Ap , Bjp from A, Bj
Pack Bjp → B̃p
for ( i = 0 : m − 1 : mC )
Identify Api , Cji from Ap , Cj
Pack Api → Ãi
for ( h = 0 : nC − 1 : nR )
Identify B̃ph , Cjih from B̃p , Cji
for ( l = 0 : mC − 1 : mR )
Identify Ãil , Cjihl from Ãi , Cjih
Cjihl := rkern( Ãil , B̃ph , Cjihl )

Fig. 3.1. Abbreviated pseudo-code for implementing the general matrix multiplication algorithm
depicted in Figure 1.1. Here, rkern calls a real domain gemm micro-kernel.

Algorithm: [ C ] := 1m ? bp( A, B, C ) [ C ] := vk1m( A, B, C )

Set bool colStore if rs( C ) = 1 Acquire workspace W
for ( j = 0 : n − 1 : nC ) Determine if using W ; set usew
Identify Bj , Cj from B, C if ( usew )
for ( p = 0 : k − 1 : kC ) Alias Cuse ← W , Cin ← 0
Identify Ap , Bjp from A, Bj else
if colStore Pack1r Bjp → B̃p Alias Cuse ← C, Cin ← C
else Pack1e Bjp → B̃p Set bool colStore if rs(Cuse ) = 1
for ( i = 0 : m − 1 : mC ) if ( colStore ) cs(Cuse ) ×= 2
Identify Api , Cji from Ap , Cj else rs(Cuse ) ×= 2
if colStore Pack1e Api → Ãi n(A) ×= 2; m(B) ×= 2
else Pack1r Api → Ãi Cuse := rkern( A, B, Cin )
for ( h = 0 : nC − 1 : nR ) if ( usew )
Identify B̃ph , Cjih from B̃p , Cji C := W
for ( l = 0 : mC − 1 : mR )
Identify Ãil , Cjihl from Ãi , Cjih
Cjihl := vk1m( Ãil , B̃ph , Cjihl )

Fig. 3.2. Left: Pseudo-code for Algorithms 1m c bp and 1m r bp, which result from applying
1m c and 1m r algorithmic variants to the block-panel algorithm depicted in Figure 1.1. Here,
Pack1e and Pack1r pack matrices into the 1e and 1r formats, respectively. Right: Pseudo-code
for a virtual micro-kernel used by all 1m algorithms.

respectively. Figure 3.2 shows a hybrid algorithm that encompasses both, supporting
row- and column-stored matrices C.
As with the 3m and 4m algorithms in [21], we have separated the so-called virtual
micro-kernel into a separate function, shown in Figure 3.2 (right). This 1m-specific
virtual micro-kernel, vk1m, largely consists of a call to the real domain micro-kernel
12 FIELD G. VAN ZEE

rkern, with some added special case handling and book-keeping needed to properly
induce complex matrix multiplication. Some of the details of the virtual micro-kernel
will be addressed later.
3.5. Performance properties. Table 3.2 tallies the total number of f.e. mem-
ops required by the block-panel algorithm for both variants of 1m (1m c and 1m r).
For comparison, we also include the corresponding memop counts for a selection of
4m algorithms as well as a conventional assembly-based solution, as first published in
Table III in [21].
Notice that 1m c bp (and 1m r bp) incur additional memops relative to a con-
ventional assembly-based solution. This stems from the fact that, unlike the assembly-
based solution, 1m implementations cannot reuse13 all real and imaginary f.e. from
vector registers.
We can hypothesize that the observed performance signatures of 1m c bp and
1m r bp may be slightly different, because each places the additional memop overhead
that is unique to 1m on different parts of the computation. This stems from the fact
that there exists an asymmetry in the assignment of packing formats to matrices in
each 1m variant. Specifically, 50% more memops—relative to a conventional assembly
solution—are required during the initial packing and the movement between caches
for the matrix packed according to 1e, since that format writes four f.e. for every
two that it reads from the source operand. (Packing to 1r incurs the same number
of memops as an assembly-based solution.) Also, if 1m c bp and 1m r bp use real
micro-kernels with different micro-tile shapes (i.e., different values of mR and nR ),
those micro-kernels’ differing performance properties will likely cause the performance
signatures of 1m c bp and 1m r bp to deviate further.
Table 3.3 summarizes Table 3.2 and also lists the level of the memory hierarchy
from which each matrix operand is reused as well as a measure of memory movement
efficiency. The information listed for 4m and assembly algorithms is reproduced from
Table IV of the previous article.
3.6. Algorithm details. This section lays out important details that must be
handled when implementing the 1m method.
3.6.1. Micro-kernel I/O preference. Within the BLIS framework, micro-
kernels are registered with a property that describes their input/output preference.
The I/O preference describes whether the micro-kernel is set up to ideally use vector
instructions to load and store elements of the micro-tile by rows or by columns. This
property typically originates from the semantic orientation of vector registers used to
accumulate the mR × nR micro-panel product. Whenever possible, the BLIS frame-
work will perform logical transpositions14 so that the apparent storage of C matches
the preference property of the micro-kernel being used. This guarantees that the
micro-kernel will be able to load and store f.e. of C using vector instructions.
This preference property is merely an interesting performance detail for conven-
tional implementations (real and complex). However, in the case of 1m, it becomes
important for constructing a correctly-functioning implementation. Specifically, the
micro-kernel’s I/O preference determines whether the 1m c or 1m r algorithm is pre-
scribed. Generally speaking, a 1m c algorithmic variant must employ a micro-kernel
that prefers to access C by columns, while a 1m r algorithmic variant must use a
13
Here, the term “reuse” refers to the same reuse described in Footnote 7.
14
This amounts to a swapping of the row and column strides and a swapping of the m and n
dimensions.
 COMPLEX MATRIX MULTIPLICATION VIA THE 1M METHOD 13
Table 3.2
f.e. memops incurred by various algorithms, broken down by stage of computation

f.e. memops required to . . . b

a
Algorithm update micro- pack move Ãi from pack move B̃p from
tilesc C r , C i Ãi L2 to L1 cache B̃p L3 to L1 cache
4m h 8mn kkC 8mk nnC 4mk nnR 8kn 4kn mmC
4m 1b 8mn kkC 8mk n2nC 4mk nnR 8kn 2m
4kn mC

4m 1a 8mn k2kC 8mk nnC 4mk nnR 8kn 4kn mmC

assembly 4mn kkC 4mk nnC 2mk nnR 4kn 2kn mmC
1m c bp 6mk nnC 4mk nnR 4kn 2m
2kn m
4mn k2kC C

1m r bp 4mk nnC 2mk n2nR 6kn 4kn mmC

a Algorithms 4m h, 4m 1b, 4m 1a, and the assembly implementation employ a block-panel algo-
rithm, and therefore their memop counts more closely resemble those of 1m c bp and 1m r bp.
b We express the number of iterations executed in the 5th, 4th, 3rd, and 2nd loops as nn , kk , mm ,
C C C
and nn . The precise number of iterations along a dimension x using a cache blocksize xC would
R
actually be d xx e. Similarly, when blocksize scaling of 21 is required, the precise value x x/2
 
C C
is expressed as x2x . These simplifications allow easier comparison between algorithms while still
C
providing meaningful approximations.
c As described in Section 3.6.2, mR × nR workspace sometimes becomes mandatory, such as when
β i 6= 0. When workspace is employed in a 4m-based algorithm, the number of f.e. memops
incurred updating the micro-tile typically doubles from the values shown here.

micro-kernel that prefers to access C by rows.

3.6.2. Workspace. In some cases, a small amount of mR × nR workspace is
needed. These cases fall into one of four scenarios: (1) C is row-stored and the real
micro-kernel rkern has a column preference; (2) C is column-stored and rkern has
a row preference; (3) C is general-stored (i.e., neither rs(C) nor cs(C) is unit); and
(4) β has a non-zero imaginary component. If any of these situations apply, then the
1m virtual micro-kernel will need to use workspace. This corresponds to the setting
of usew in vk1m (in Figure 3.2). The idea is simply that rkern will be called to
compute the micro-panel product and store it to the workspace W . Subsequently, the
result in W can be accumulated back to C.
Cases (1) and (2), while supported, actually never occur in practice because BLIS
will perform (at a high level within the framework) a logical transposition whenever
necessary. The net effect is that the storage of C will always appear to match the I/O
preference of the micro-kernel.
Case (3) is needed because the real micro-kernel is programmed to support the
updating of real matrices stored with general stride, which cannot be spoofed to
match the updating of complex matrices stored with general stride. The reason is
even when stored with general stride, complex matrices store real and imaginary f.e.
in contiguous pairs. There is no way to coax this pattern of data access from a real
domain micro-kernel, given its existing API. Thus, general stride support must be
implemented outside rkern, within vk1m.
Case (4) is needed because real domain micro-kernels are not capable of scaling
14 FIELD G. VAN ZEE

Table 3.3
Performance properties of various algorithms

Level from which f.e. of ma-

trix X are reused, and lL1 :
Algorithm

Total f.e. memops required # of times each cache line is

(Sum of columns of Table 3.2) moved into the L1 cache (per
rank-kC update).
C A B
C lL1 A lL1 B lL1
     
k 2n n 2m
4m h 8mn kC
+ 4mk nC
+ nR
+ 2kn 4 + mC
Mem 4 Mem 4 Mem 4
     
4m 1b 8mn k
kC
+ 4mk 4n
nC
+ n
nR
+ 2kn 4 + 4m
mC
L2 2a L2 1 L1 1
     
2k 2n n 2m a
4m 1a 8mn kC
+ 4mk nC
+ nR
+ 2kn 4 + mC
L1 1 L1 1 L1 1
     
k 2n n m
assembly 4mn kC
+ 2mk nC
+ nR
+ 2kn 2+ mC
Reg 1 Reg 1 Reg 1
     
1m c bp 4mn 2k
kC
+ 2mk 3n
nC
+ 2n
nR
+ 2kn 2 + 2m
mC
Reg 1 L2b 1 Reg 1
     
2k 2n 2n 2m b
1m r bp 4mn kC
+ 2mk nC
+ nR
+ 2kn 3 + mC
Reg 1 Reg 1 L1 1

a This assumes that the micro-tile is not evicted from the L1 cache during the next call to rkern.
b In the case of 1m algorithms, we consider f.e. of A and B to be “reused” from the level of cache
in which the 1e-formatted matrix resides.

C by complex scalars β (that is, β such that β i 6= 0).

3.6.3. Handling alpha and beta scalars. As in the previous article, we have
simplified the general matrix multiplication to C := C+AB. In practice, the operation
is implemented as C := βC + αAB, where α, β ∈ C. Let us use Algorithms 1m c bp
and 1m r bp in Figure 3.2 as an example of how to support arbitrary values of α and
β.
If no workspace is needed (because none of the four situations described in Sec-
tion 3.6.2 apply), we can simply pass β r into the rkern call. However, if β is complex,
or (regardless of whether β is real or complex) if any of the other three workspace cases
apply, then we must pass in a local βuse = 0 to rkern, compute to local workspace
W , and then apply β at the end of vk1m, when W is accumulated to C.
When α is real, the scaling may be performed directly by rkern. This situation
is ideal since it almost always incurs no additional costs (since many micro-kernels
multiply their intermediate AB product by α unconditionally). Scaling by α with non-
zero imaginary components can be performed by the packing function when either Ãi
or B̃p are packed. Though somewhat less than ideal, the overhead incurred by this
treatment of α is minimal since packing is a memory-bound operation.
3.6.4. Multithreading. As with Algorithm 4m 1a in the previous article, Al-
gorithms 1m c bp and 1m r bp parallelize in a straightforward manner for multicore
and many-core environments. Because those algorithms encode the 1m method en-
tirely within the packing functions and the virtual micro-kernel, all other levels of
code are completely oblivious to, and therefore unaffected by, the specifics of the new
algorithms. Therefore, we expect that 1m c bp and 1m r bp will yield multithreaded
performance that is on-par with that of the corresponding real domain matrix multi-
plication function, rmmbp.
 COMPLEX MATRIX MULTIPLICATION VIA THE 1M METHOD 15

3.6.5. Bypassing the virtual micro-kernel. Because the 1m virtual micro-

kernel serves as a wrapper to the real domain micro-kernel, it would seem at first glance
that there exists the potential for additional overhead in 1m algorithms, particularly
from the extra function calls. However, there are a few things to consider.
First, we should consider that a conventional solution would implement matrix
multiplication using a complex micro-kernel, which usually has a smaller micro-tile
footprint (i.e., fewer f.e.). But, a complex micro-kernel that updates fewer f.e. would
need to be called more times per rank-kC update in order to fully update the output
matrix C. Thus, the function call overhead incurred by 1m algorithms may already
be at or near parity with that of a conventional implementation.
Secondly, even if the complex micro-kernel updates the same number of f.e. as
its real domain counterpart, there exists a simple optimization that can be applied as
long as β i = 0 and C is either row- or column-stored (but not general-stored). If these
conditions are met and detected by the implementation, the real domain macro-kernel
can be called with modified parameters to induce the equivalent complex domain
subproblem. This simple optimization avoids all overhead introduced by the virtual
micro-kernel, including (but not limited to) the cost incurred by additional function
calls.
Finally, we suspect that, even if this optimization cannot be applied, the slowdown
that results from the additional overhead may be tolerable to many applications.

3.7. 2m. The previous article noted that because of its expression in terms of
real and imaginary matrices, 4m would perform more favorably on complex matrices
that stored real and imaginary values separately, as two real matrices. This would not
benefit the accesses on A and B since those matrices must almost always be packed to
achieve high performance and can easily be separated during that process. However,
it would benefit the accesses on C. Referring back to Eq. 3.2, we can see that storing
the real and imaginary f.e. of C separately is equivalent to a permutation of the
rows of C. In order to keep the computation expressed unchanged, this permutation
would need to be applied to matrix A as well since they both share an m dimension.
Applying such a permutation to Eq. 3.2 yields:

r r r r r i r i
   
γ00 γ01 γ02 γ03 α00 −α00 α01 −α01
r r r r r i r i
 r r r r

 γ10 γ11 γ12 γ13   α10 −α10 α11 −α11  β00 β01 β02 β03
 r r r r
  r i r i
 i i i i 
 γ20 γ21 γ22 γ23  +=  α20
  −α20 α21 −α21   β00 β01 β02 β03
(3.5) 
 γi i i i i r i r
 r
  β10 r r r 

 00 γ01 γ02 γ03 

α
 00 α00 α01 α01  β11 β12 β13
i i i i i r i r  β10i i i i
 γ10 γ11 γ12 γ13   α10 α10 α11 α11 β11 β12 β13
i i i i i r i r
γ20 γ21 γ22 γ23 α20 α20 α21 α21

If we also permute even-indexed columns of A to be consecutive with one another

(and grouped together) while odd-numbered columns are permuted to immediately
follow them, and also permute rows of B accordingly, we have:

r r r r r r i i
   
γ00 γ01 γ02 γ03 α00 α01 −α00 −α01
r r r r r r i i
 r r r r

 γ10
 r γ11 γ12 γ13 

 α10
 r α11 −α10 −α11  β00
 β01 β02 β03
r r r r i i r r r r 
 γ20
 i γ21 γ22 γ23   α20
 +=  i α21 −α20 −α21   β10
 i β11 β12 β13 
γ i i i i r r i i i 
 00 γ01 γ02 γ03 

α
 00 α01 α00 α01  β
 00 β01 β02 β03
i
 γ10 i i i i i r r i i i i
γ11 γ12 γ13   α10 α11 α10 α11  β10 β11 β12 β13
i i i i i i r r
γ20 γ21 γ22 γ23 α20 α21 α20 α21
16 FIELD G. VAN ZEE

Or, more generally:

Cr Ar −Ai Br
    
(3.6) +=
Ci Ai Ar Bi

This matrix multiplication can be computed via just two calls to a real domain matrix
multiplication primitive:

Br Br
   
C r += Ar −Ai , C i += Ai Ar



Bi Bi

Br
 
Ar −Ai , Ai Ar , and



provided that can each be stored so that they
Bi
can be referenced as single matrices. We call this the 2m method and refer to the
specific instance derived from Eq. 3.6 as 2m c.
Applying similar permutations to the n and k dimensions to Eq. 3.4 yields:

Br Bi
 
r i r i



(3.7) C C += A A
−B i B r

which can also be broken down in two instances of real matrix multiplication:

Br Bi
   
r r i i r i



C += A A , C += A A
−B i Br

which corresponds to 2m r.
We can now make a few observations about 2m:
• Eqs. 3.6 and 3.7 are identical to Eqs. 3.1 and 3.3, respectively, except that
scalars are replaced with matrices.
• The permutation along the k dimension is actually unnecessary and so the
formatting captured by Eq. 3.5 also falls under 2m. This permutation (or lack
thereof) only changes the order in which intermediate terms are accumulated.
• The storage of C r and C i in both Eqs. 3.6 and 3.7 is unspecified and does not
depend on which matrix—A or B—was originally formatted with 1e (prior
to permutation). Either C r or C i may be stored by rows, columns, or a more
exotic storage scheme, and their storage formats need not even be identical.
Thus, neither 2m c nor 2m r implies the storage of C; rather, they only
imply how matrices A and B are formatted and stored—that is, which one
contains duplicated (and negated) f.e..
• The 2m method can be applied at an arbitrary level of matrix multiplication.
For example, if we assume from Eq. 3.6 that input matrices Ar −Ai ,


 r
B
Ai Ar , and


are each stored as micro-panels (stored by columns,
Bi
columns, and rows, respectively), then the implied primitive is the real gemm
micro-kernel, which will update mR × nR micro-tiles of C.
4. Performance. In this section we present performance results for implemen-
tations of 1m algorithms on a recent Intel architecture. For comparison, we include
results for a 4m algorithm as well as conventional assembly-based approaches in the
real and complex domains.
 COMPLEX MATRIX MULTIPLICATION VIA THE 1M METHOD 17
Table 4.1
Register and cache blocksizes used by various BLIS implementations of matrix multiplication,
as configured for an Intel Xeon E5-2690 v3 “Haswell” processor

Precision/Domain Implementation mzR nzR mzC z

kC nzC
BLIS 1m c 16/2 6 144/2 256/2 4080
BLIS 1m r 6 16/2 144 256/2 4080/2
single complex
BLIS assembly (c) 8 3 56 256 4080
BLIS assembly (r) 3 8 75 256 4080
BLIS 1m c 8/2 6 72/2 256/2 4080
BLIS 1m r 6 8/2 72 256/2 4080/2
double complex
BLIS assembly (c) 4 3 44 256 4080
BLIS assembly (r) 3 4 192 256 4080

Note: For 1m implementations, division by 2 is made explicit to allow the reader to quickly see both
the complex blocksize values as well as the values that would be used by the underlying real domain
micro-kernels when performing real matrix multiplication. The I/O preference of the assembly-based
implementations is indicated by a “(c)” or “(r)” (for column- or row-preferring).

4.1. Platform and implementation details. Results presented in this section

were gathered on a single Cray XC40 compute node consisting of two 12-core Intel
Xeon E5-2690 v3 processors featuring the “Haswell” microarchitecture. Each core,
running at a clock rate of 3.2 GHz15 , provides a single-core peak performance of 51.2
gigaflops (GFLOPS) in double precision and 102.4 GFLOPS in single precision.16
Each socket has a 30MB L3 cache that is shared among cores, and each core has a
private 256KB L2 cache and 32KB L1 (data) cache. Performance experiments were
gathered under the Cray Linux Environment 6 operating system running the Linux
4.4.103 (x86 64) kernel. Source code was compiled by the GNU C compiler (gcc)
version [Link] The version of BLIS used in these tests was not officially released at
the time of this writing, and was adapted from version 0.6.0-11.18
Algorithms 1m c bp and 1m r bp were implemented in the BLIS framework,
as described in Section 3.4. We also refer to results based on existing conventional
assembly-based micro-kernels written by hand for the Haswell microarchitecture via
GNU extended inline assembly syntax.
All experiments were performed on randomized, column-stored matrices with
gemm scalars held constant: α = β = 1. In all performance graphs, each data
point represents the best of three trials.
Blocksizes for each of the BLIS implementations tested are provided in Table 4.1.
In all graphs presented in this section the x-axes denote the problem size, the
y-axes show observed floating-point performance in units of GFLOPS per core, and
15 This system uses Intel’s Turbo Boost 2.0 dynamic frequency throttling technology. According

to [15], the maximum the clock frequency when executing AVX instructions is 3.2 GHz when utilizing
one or two cores, and 3.0 GHz when utilizing three or more cores.
16 Accounting for the reduced AVX clock frequency, the peak performance when utilizing 24 cores

is 48 GFLOPS/core in double precision and 96 GFLOPS/core in single precision.

17 The following optimization flags were used during compilation of BLIS and its test drivers: -O3

-mavx2 -mfma -mfpmath=sse -march=haswell.

18 Despite not yet having an official version number, this version of BLIS may be uniquely

identified, with high probability, by the first 10 digits of its git “commit” (SHA1 hash) number:
ceee2f973e.
18 FIELD G. VAN ZEE

gemm, single-precision (1 thread) gemm, single-precision (1 thread)

100 100

80 80
GFLOPS

GFLOPS
60 60

40 40
cgemm assembly (c) cgemm assembly (r)
sgemm assembly (c) sgemm assembly (r)
20 20
cgemm 1m_c_bp cgemm 1m_r_bp
cgemm 4m_1a (c) cgemm 4m_1a (r)
0 0
0 500 1000 1500 2000 0 500 1000 1500 2000
problem size (m = n = k) problem size (m = n = k)
gemm, double-precision (1 thread) gemm, double-precision (1 thread)
50 50

40 40
GFLOPS

GFLOPS

30 30

20 20
zgemm assembly (c) zgemm assembly (r)
dgemm assembly (c) dgemm assembly (r)
10 10
zgemm 1m_c_bp zgemm 1m_r_bp
zgemm 4m_1a (c) zgemm 4m_1a (r)
0 0
0 500 1000 1500 2000 0 500 1000 1500 2000
problem size (m = n = k) problem size (m = n = k)

Fig. 4.1. Single-threaded performance of various implementations of single-precision (top) and

double-precision (bottom) complex gemm on a single core of an Intel Xeon E5-2690 v3 “Haswell”
processor. The left and right graphs differ in which 1m implementation they report, with the left
graphs reporting 1m c bp (which employs a column-preferring micro-kernel) and the right graphs
reporting 1m r bp (which employs a row-preferring micro-kernel). The graphs also contain three
reference curves for comparison: an assembly-coded complex gemm, an assembly-coded real gemm,
and the 4m 1a implementation found in BLIS (with the latter two using the same micro-kernel
as the 1m implementation shown in the same graph). For consistency with the 1m curves, these
reference implementations differ from left to right graphs in the I/O preference of their underlying
micro-kernel, indicated by a “(c)” or “(r)” (for column- or row-preferring) in the legends. The
theoretical peak performance coincides with the top of each graph.

the theoretical peak performance coincides with the top of each graph.
4.2. Sequential results. Figure 4.1 reports performance results for various im-
plementations of double- and single-precision complex matrix multiplication on a
single core of the Haswell processor. For these results, all matrix dimensions were
equal (e.g. m = n = k). Results for 1m c bp (which uses a column-preferring
micro-kernel) appears on the left of Figure 4.1 while those of 1m r bp (which uses a
row-preferring micro-kernel) appears on the right.
Each graph in Figure 4.1 also contains three reference implementations: BLIS’s
complex gemm based on conventional assembly-coded kernels (e.g. “cgemm assem-
bly”); BLIS’s real gemm (e.g. “sgemm assembly”); and the 4m 1a implementation
found in BLIS.19 We configured all three of these reference codes to use column-

19 Within any given graph of Figures 4.1 and 4.2, the 1m and 4m 1a implementations use the same
 COMPLEX MATRIX MULTIPLICATION VIA THE 1M METHOD 19

preferential micro-kernels on the left and row-preferential micro-kernels on the right,

as indicated by a “(c)” or “(r)” in the legends, in order to provide consistency with
the 1m results.
As predicted in Section 3.5, we find that the performance signatures of the
1m c bp and 1m r bp algorithms differ slightly. This was expected given that the
1e and 1r packing formats place different memory access burdens on different packed
matrices, Ãi and B̃p , which reside in different levels of cache. It was not previously
clear, however, which would be superior over the other. It seems that, at least in
the sequential case, the difference is somewhat more noticeable in double-precision,
though even there it is quite subtle. This difference is almost certainly due to the
individual performance characteristics of the underlying row- and column-preferential
micro-kernels. We find evidence of this in the 4m 1a results, which was also affected
by the change in micro-kernel I/O preference.
In all cases, the 1m implementations outperform 4m 1a, with the margin some-
what larger in single-precision.
The 1m implementations match or exceed the performance of their real domain
gemm benchmarks (the dotted lines in each graph), and are quite competitive with
assembly-coded complex gemm (the solid lines) regardless of the algorithm employed.

4.3. Multithreaded results. Figure 4.2 shows single- and double-precision per-
formance using 24 threads, with one thread bound to each physical core of the proces-
sor. Performance is presented in units of gigaflops per core to facilitate visual assess-
ment of scalability. For all BLIS implementations, we employed 4-way parallelism
within the 5th loop, 3-way parallelism within the 3rd loop, and 2-way parallelism in
the 2nd loop for a total of 24 threads. This parallelization scheme was chosen in a
manner consistent with that of the previous article using a strategy set forth in [18].
Compared to the single-threaded case, we find a more noticable difference in
multithreaded performance between the 1m algorithms. Specifically, the 1m r bp
implementation (based on a row-preferring micro-kernel) outperforms that of 1m c bp
(based on a column-preferring micro-kernel), with the difference more pronounced in
single-precision. We suspect this is rooted not in the algorithms, per se, but in the
differing micro-kernel implementations used by each 1m algorithm. The 1m r bp
algorithm is implemented with a real micro-kernel that is 6 × 16 and 6 × 8 in the
single- and double-precision cases, respectively, while 1m c bp uses 16 × 6 and 8 ×
6 micro-kernels for single- and double-precision implementations, respectively. The
observed difference in performance between the 1m algorithms is likely attributable
to the fact that the micro-kernels’ different values for mR and nR place different
bandwidth requirements when reading f.e. from the caches (primarily L1 and L2).
More specifically, larger values of mR place a heavier burden on loading elements from
the L2 cache, which is usually disadvantageous since that cache can typically sustain
lower bandwidth. By contrast, a micro-kernel with larger nR loads more elements
(per mR × nR rank-1 update) from the L1 cache, which resides closer to the processor
and can typically sustain higher bandwidth than the L2 cache.
The multithreaded 1m implementation approximately matches or exceeds its real
domain counterpart in all cases.
The 1m algorithm based on a row-preferential micro-kernel, 1m r bp, outper-
forms 4m 1a, especially in single-precision where the margin is quite wide. The
1m algorithm based on column-preferential micro-kernels, 1m c bp, performs more

real-domain micro-kernel that of the real gemm (e.g. “sgemm assembly” or “dgemm assembly”).
20 FIELD G. VAN ZEE

gemm, single-precision (24 threads) gemm, single-precision (24 threads)

80 80
GFLOPS/core

GFLOPS/core
60 60

40 40
cgemm assembly (c) cgemm assembly (r)
20 sgemm assembly (c) 20 sgemm assembly (r)
cgemm 1m_c_bp cgemm 1m_r_bp
cgemm 4m_1a (c) cgemm 4m_1a (r)
0 0
0 2000 4000 6000 0 2000 4000 6000
problem size (m = n = k) problem size (m = n = k)
gemm, double-precision (24 threads) gemm, double-precision (24 threads)

40 40
GFLOPS/core

GFLOPS/core

30 30

20 20
zgemm assembly (c) zgemm assembly (r)
10 dgemm assembly (c) 10 dgemm assembly (r)
zgemm 1m_c_bp zgemm 1m_r_bp
zgemm 4m_1a (c) zgemm 4m_1a (r)
0 0
0 2000 4000 6000 0 2000 4000 6000
problem size (m = n = k) problem size (m = n = k)

Fig. 4.2. Multithreaded performance of various implementations of single-precision (top) and

double-precision (bottom) complex gemm on two Intel Xeon E5-2690 v3 “Haswell” processors, each
with 12 cores. All data points reflect the use of 24 threads. The left and right graphs differ in which
1m implementation they report, with the left graphs reporting 1m c bp (which employs a column-
preferring micro-kernel) and the right graphs reporting 1m r bp (which employs a row-preferring
micro-kernel). The graphs also contain three reference curves for comparison: an assembly-coded
complex gemm, an assembly-coded real gemm, and the 4m 1a implementation found in BLIS (with
the latter two using the same micro-kernel as the 1m implementation shown in the same graph).
For consistency with the 1m curves, these reference implementations differ from left to right graphs
in the I/O preference of their underlying micro-kernel, indicated by a “(c)” or “(r)” (for column-
or row-preferring) in the legends. The theoretical peak performance coincides with the top of each
graph.

poorly, barely ediging out 4m 1a in single precision and tracking closely with 4m 1a
in double precision. We suspect that 4m 1a is more resilient to the lower-performing
column-preferential micro-kernel by virtue of the fact that the algorithm’s virtual
micro-kernel leans heavily on the L1 cache, which on this architecture is capable of
being read from and written to at relatively high bandwidth (64 bytes/cycle and 32
bytes/cycle, respectively) [14].
4.4. Comparing to other implementations. While our primary goal is not
to compare the performance of the newly developed 1m implementations with that
of other established BLAS solutions, some basic comparison is merited and thus we
have included Figure 4.3 (left). These graphs are similar to those in Figure 4.1,
except that: we show only implementations based on row-preferential micro-kernels;
we omit 4m 1a; and we include results for complex gemm implementations provided
 COMPLEX MATRIX MULTIPLICATION VIA THE 1M METHOD 21

gemm, single-precision (1 thread) gemm, single-precision (24 threads)

100

80
80

GFLOPS/core
60
GFLOPS

60

40 cgemm MKL
40 cgemm MKL
cgemm assembly (r) cgemm assembly (r)
sgemm assembly (r) 20 sgemm assembly (r)
20
cgemm 1m_r_bp cgemm 1m_r_bp
cgemm OpenBLAS cgemm OpenBLAS
0 0
0 500 1000 1500 2000 0 2000 4000 6000
problem size (m = n = k) problem size (m = n = k)
gemm, double-precision (1 thread) gemm, double-precision (24 threads)
50

40
40
GFLOPS/core

30
GFLOPS

30

20 zgemm MKL 20 zgemm MKL

zgemm assembly (r) zgemm assembly (r)
dgemm assembly (r) 10 dgemm assembly (r)
10
zgemm 1m_r_bp zgemm 1m_r_bp
zgemm OpenBLAS zgemm OpenBLAS
0 0
0 500 1000 1500 2000 0 2000 4000 6000
problem size (m = n = k) problem size (m = n = k)

Fig. 4.3. Single-threaded (left) and multithreaded (right) performance of various implementa-
tions of single-precision (top) and double-precision (bottom) complex gemm on a single core (left)
or 12 cores (right) of an Intel Xeon E5-2690 v3 “Haswell” processor. All multithreaded data points
reflect the use of 24 threads. The 1m curves are identical from those shown in Figures 4.1 and 4.2.
The theoretical peak performance coincides with the top of each graph.

by OpenBLAS 0.3.6 [17] and Intel MKL 2019 Update 4 [13].

Figure 4.3 (right) shows multithreaded performance of the same implementations
running with 24 threads.
These graphs show that BLIS’s complex assembly-based and 1m implementations
typically outperform OpenBLAS while falling short in most (but not all) cases when
compared to Intel’s MKL library.
4.5. Additional results. Additional performance results were gathered on a 56-
core Marvell ThunderX2 compute server. For brevity, we present and discuss that data
in Appendix A. Those results reinforce the narrative provided here, lending even more
evidence that the 1m method is capable of yielding high-performance implementations
of complex matrix multiplication that are competitive with (and often outperform)
other leading library solutions.
5. Observations.
5.1. 4m limitations circumvented. The previous article concluded by iden-
tifying a number of limitations inherent in the 4m method that collectively prevent
the approach from becoming both a feasible and competitive alternative to matrix
22 FIELD G. VAN ZEE

multiplication via conventional assembly-based kernels. We now revisit this list and
briefly discuss whether, to what degree, and how those limitations are overcome by
algorithms based on the 1m method.
Number of calls to primitive. The most versatile 4m algorithm, 4m 1a, incurs
up to a four-fold increase in function call overhead over a comparable assembly-based
implementation. By comparison, 1m algorithms require at most a doubling of micro-
kernel function call overhead, and in certain common cases (e.g., when β ∈ R and C
is row- or column-stored), this overhead can be avoided completely. The 1m method
is clearly an improvement over 4m due to its reliance on a single invocation of the
matrix multiplication primitive.
Inefficient reuse of input data from A, B, and C. The most cache-efficient
application of 4m is the lowest level algorithm, 4m 1a, which reuses f.e. of A, B, and
C from the L1 cache. But, as shown in Table 3.3, both 1m r and 1m c variants reuse
f.e. of two of the three matrices from registers, with 1m r bp reusing f.e. of the third
matrix from the L1 cache. This would seem to be a significant improvement, though
observed performance improvement over 4m 1a will depend on properties specific to
the hardware.
Non-contiguous output to C. Algorithms based on the 4m method must up-
date only the real and then only the imaginary parts of the output matrix, twice each.
Since C is typically stored (by rows or columns) with real and imaginary f.e. inter-
leaved, this piecemeal approach prevents the real micro-kernel from using vector load
and store instructions on C during those four updates. The 1m method avoids this
issue altogether by packing A and B to formats that allow the real micro-kernel to
update contiguous real and imaginary f.e. of C simultaneously within a single invo-
cation. We suspect that this is, perhaps, the largest contributor to 1m’s performance
superiority over 4m.
Reduction of kC . Algorithm 4m 1a requires that the real micro-kernel’s pre-
ferred kC blocksize be halved in the complex algorithm in order to maintain proper
cache footprints of Ãi and B̃p as well the footprints of their constituent micro-panels.20
Using such sub-optimally sized micro-panels can noticeably hobble the performance
of 4m 1a. Looking back at Table 3.1, it may seem like 1m suffers a similar handicap;
however, the reason for halving kC and its effect are both completely different. In the
case of 1m, the use of kCz
= 12 kC is simply a conversion of units (complex elements to
real f.e.) for the purposes of identifying the size of the complex submatrices to be
packed that will induce the optimal kC value from the perspective of the real micro-
kernel, not a reduction in the f.e. footprint of the micro-panels operated upon by
that real micro-kernel. Indeed, the ability of 1m to achieve high performance when
k = 21 kC is a strength in the context of certain higher-level applications, such as
Cholesky, LU, and QR factorizations based on rank-k update. Those operations tend
to perform better when the algorithmic blocksize (corresponding to kC ) is as nar-
row as possible in order to limit the amount of computation in the lower-performing
unblocked subproblem.
Framework accommodation. The 1m algorithms are no more disruptive to the
BLIS framework than the most accommodating of 4m algorithms, 4m 1a, and much
20 Recall that the halving of k
C for 4m 1a was motivated by the desire to keep not just two, but
four real micro-panels in the L1 cache simultaneously. These correspond to the real and imaginary
parts of the current micro-panels of Ãi and B̃p .
 COMPLEX MATRIX MULTIPLICATION VIA THE 1M METHOD 23

less disruptive than the remaining algorithms. This is because, like with 4m 1a, al-
most all of the 1m implementation details are sequestered within the packing routines
and the virtual micro-kernel.
Interference with multithreading. Because the 1m algorithms are imple-
mented entirely within the packing facility and virtual micro-kernel, they parallelize
just as easily as the most thread-friendly of the 4m algorithms, 4m 1a, and entirely
avoid the threading difficulties of higher-level 4m algorithms.21
Non-applicability to two-operand operations. Certain higher-level appli-
cations of 4m are inherently incompatible with two-operand operations because they
would overwrite the original contents of the input/output operand even though subse-
quent stages of computation depend on that original input. 1m avoids this limitation
entirely. Like 4m 1a, 1m can easily be applied to two-operand level-3 operations such
as trmm and trsm.22
Summary. The analysis above suggests that the 1m method solves or avoids
most of the performance-degrading weaknesses of 4m and in the remaining cases is
no worse off than the best 4m algorithm. Thus, its observed performance superiority
was predictable.
5.2. Further discussion. Before concluding, here we offer some final thoughts
on the 1m method and its place in the larger spectrum of approaches to implementing
complex matrix multiplication.
5.2.1. Geometric interpretation. Matrix multiplication is sometimes thought
of as a three-dimensional operation with a contraction (accumulation) over the k di-
mension. This interpretation carries into the complex domain as well. However, when
each complex element is viewed in terms of its real and imaginary components, we
find that a fourth pseudo-dimension of computation (of fixed size 2) emerges, one
which also involves a contraction. The 1m method reorders and duplicates elements
of A and B in such a way that exposes and “flattens” this extra dimension of com-
putation. This, combined with the exposed treatment of real and imaginary f.e.,
causes the resulting floating-point operations to appear indistinguishable from a real
domain matrix multiplication with m and k dimensions (for column-stored C) or k
and n dimensions (for row-stored C) that are twice as long.
5.2.2. Data reuse: efficiency vs. programmability. Both the conventional
approach and 1m move data efficiently through the memory hierarchy.23 However,
once in registers, a conventional complex micro-kernel reuses those loaded values
to perform twice as many flops as 1m. The previous article observes that all 4m
algorithms make different variations of the same tradeoff: by forgoing the reuse of
f.e. from registers and instead reusing those data from some level of cache, the
algorithms avoid the need to explicitly encode complex arithmetic at the assembly
level. As it turns out, 1m makes a similar tradeoff, but gives up less while gaining
more: it is able to effectively reuse f.e. from two of the three matrix operands from
registers while still avoiding the need for a complex micro-kernel, and it manages

21 This thread-friendly property holds even when the virtual micro-kernel is bypassed altogether

as discussed in Section 3.6.5

22 As with 4m 1a, 1m support for trsm requires a separate pair of virtual micro-kernels that fuse

a matrix multiplication with a triangular solve with nR right-hand sides.

23 This is in contrast to, for example, Algorithm 4m hw, which the previous article showed makes

rather inefficient use of cache lines as they travel through the L3, L2, and L1 caches.
24 FIELD G. VAN ZEE

to replace that kernel operation with a single real matrix multiplication. And we
would argue that increasing programmability and productivity by forfeiting a modest
performance advantage is a good trade to make under almost any circumstance.
5.2.3. Micro-kernel bandwidth. Because 1m algorithms do not explicitly re-
use f.e. of Ãi and B̃p from registers, they require higher memory bandwidth during
the micro-kernel computation than a conventional assembly-based solution.24 Specif-
ically, increased bandwidth is utilized when reading from the copy of the matrix that
is packed into the 1e format, which resides in either the L2 or L1 cache, depending
on which algorithm is being employed.25 In practice, this potential weakness of 1m
is not a concern. Yes, in principle, one could design an architecture with sufficiently
low memory bandwidth from L2 or L1 cache that a conventional complex matrix
multiplication achieves high performance while a 1m-based implementation struggles.
However, this would imply a corresponding performance shortfall in the underlying
real domain matrix kernel. Given the ubiquity and importance of real matrix multi-
plication in the scientific community, hardware vendors have great incentive to design
architectures that allow sgemm and dgemm to achieve high performance. Thus, we
would expect that 1m will remain a viable alternative for the foreseeable future.
5.2.4. Storage. The supremacy of the 1m method is closely tied to the inter-
leaved, pairwise storage of real and imaginary values—specifically, of the output ma-
trix C. If users and applications decide to begin storing complex matrices as two real
matrices (traditionally-stored, by rows or columns), one each for real and imaginary
components, the 2m approach (for numerically sensitive settings) as well as low-level
applications of 3m (for numerically insensitive settings) become more appropriate.
6. Conclusions. We began the article by reviewing the general motivations for
induced methods for complex matrix multiplication as well as the specific methods, 3m
and 4m, studied in the previous article. Then, we recast complex scalar multiplication
(and accumulation) in such a way that revealed a template that could be used to
fashion a new induced method, one that casts complex matrix multiplication in terms
of a single real matrix product. The key is the application of two new packing formats
on the left- and right-hand matrix product operands that allows us to disguise the
complex matrix multiplication as a real matrix multiplication with slightly modified
input parameters. This 1m method is shown to have two variants, depending on
whether the output matrix is stored by rows or columns. We also briefly contemplated
a related 2m method that would be applicable to more exotic storage arrangements
that separate the real and imaginary components of the matrix operands. When
implemented in the BLIS framework, competitive performance was observed for 1m
algorithms on a recent Intel microarchitecture. Finally, we reviewed the limitations of
the 4m method that are overcome by 1m and concluded by discussing a few high-level
observations.

24 While this bandwidth distinction actually holds for all induced methods, different algorithms

will place differing degrees of bandwidth pressure on the memory hierarchy, depending on the level(s)
of cache from which they reuse f.e. of A and B.
25 The bandwidth from the 1r-formatted matrix is unaffected since it only reorders (rather than

duplicates) its f.e.. And since 1m uses half the kC that is optimal in the real domain (and therefore
would perform roughly twice as many rank-kC updates), the bandwidth required for accessing f.e.
of the output matrix may also be higher, though the precise level of increase will depend on the value
of kCz that would be used by a comparable assembly-based implementation. However, bandwidth

requirements on C are already quite low, so we would not expect this difference to measurably impact
performance.
 COMPLEX MATRIX MULTIPLICATION VIA THE 1M METHOD 25

The key takeaway from our study of induced methods is that the real and imag-
inary elements of complex matrices can always be reordered to accommodate the
desired fundamental primitives, whether those primitives are defined to be various
forms of real matrix multiplication (as is the case for the 4m, 3m, 2m, and 1m meth-
ods), or vector instructions (as is the case for micro-kernels that implement complex
arithmetic in assembly code). Indeed, even in the real domain, the classic matrix
multiplication algorithm’s packing format is simply a reordering of data that targets
the fundamental primitive implicit in the micro-kernel—namely, an mR × nR rank-1
update. The family of induced methods presented here and in the previous article ex-
pand upon this basic reordering so that the mathematics of complex arithmetic can be
expressed at different levels of the algorithm and of its corresponding implementation,
each yielding different benefits, costs, and performance.
Appendix A. Additional Performance Results.
In this section we present performance results for implementations of 1m method
on a recent ARM-based architecture. The primary purpose of gathering these results
was to confirm 1m performance on an non-x86 64 architecture, and also to assess the
relative performance of 1m on a platform for which Intel’s MKL was not available.
A.1. Platform and implementation details. Results presented in this sec-
tion were gathered on a single compute node consisting of two 28-core Marvell Thun-
derX2 CN9975 processors.26 Each core, running at a clock rate of 2.2 GHz, provides
a single-core peak performance of 17.6 gigaflops (GFLOPS) in double precision and
35.2 GFLOPS in single precision. Each socket has a 32MB L3 cache that is shared
among cores, and each core has a private 256KB L2 cache and 32KB L1 (data) cache.
Performance experiments were gathered under the Ubuntu 16.04 operating system
running the Linux 4.15.0 kernel. Source code was compiled by the GNU C compiler
(gcc) version [Link] The version of BLIS used in these tests was version 0.5.0-1.28
In this section, we show results for only 1m Algorithms 1m c bp, omitting re-
sults for the other three algorithms. We chose to omit 1m r bp because we did not
develop a row-preferential micro-kernel for this architecture. (Recall that 1m r vari-
ant algorithms require a micro-kernel that reads and writes elements of C in a row
orientation.) And unlike the results shown in Section 4, we did not develop conven-
tional assembly-based micro-kernels with which to compare. For reference, we also
measured performance for the complex gemm implementations found in OpenBLAS29
and ARMPL 18.4.0.
All other parameters, such as values of α and β, and the number of trials per-
formed for each problem size, as well as graphing conventions, such as scaling of the
y-axis, remain identical to those of Section 4.
A.2. Analysis. Figure A.1 contains single-threaded (left) and multithreaded
(right) performance of single-precision (top) and double-precision (bottom) complex

26 While four-way symmetric multithreading is available on this hardware, the feature was disabled

at boot-time so that the operating system detects only one logical core per physical core and schedules
threads accordingly.
27 The following optimization flags were used during compilation of BLIS and its test drivers: -O3

-ftree-vectorize -mtune=cortex-a57. In addition to those flags, the following flags were also used
when compiling assembly kernels: -march=armv8-a+fp+simd -mcpu=cortex-a57.
28 This version of BLIS may be uniquely identified, with high probability, by the first 10 digits of

its git “commit” (SHA1 hash) number: e90e7f309b.

29 This version of OpenBLAS may be uniquely identified, with high probability, by the first 10

digits of its git commit number: 52d3f7af50.

26 FIELD G. VAN ZEE

gemm, single-precision (1 thread) gemm, single-precision (56 threads)

35 35

30 30

25 25

GFLOPS/core
GFLOPS

20 20

15 15

10 cgemm 1m_c_bp 10 cgemm 1m_c_bp

sgemm assembly (c) sgemm assembly (c)
5 cgemm ARMPL 5 cgemm ARMPL
cgemm OpenBLAS cgemm OpenBLAS
0 0
0 500 1000 1500 2000 0 2000 4000 6000 8000 10000
problem size (m = n = k) problem size (m = n = k)
gemm, double-precision (1 thread) gemm, double-precision (56 threads)

15 15
GFLOPS/core
GFLOPS

10 10

5 zgemm 1m_c_bp 5 zgemm 1m_c_bp

dgemm assembly (c) dgemm assembly (c)
zgemm ARMPL zgemm ARMPL
zgemm OpenBLAS zgemm OpenBLAS
0 0
0 500 1000 1500 2000 0 2000 4000 6000 8000 10000
problem size (m = n = k) problem size (m = n = k)

Fig. A.1. Single-threaded (left) and multithreaded (right) performance of various implementa-
tions of single-precision (top) and double-precision (bottom) complex gemm on a single core (left)
or 56 cores (right) of a Marvell ThunderX2 CN9975 processor. All multithreaded data points reflect
the use of 56 threads. The real domain gemm implementation from BLIS uses a column-preferential
microkernel, as indicated the a “(c)” in the legends. (The 1m c bp implementation uses the same
column-preferential microkernel as the real domain gemm implementation.) The theoretical peak
performance coincides with the top of each graph.

gemm implementations. In addition to the 1m c bp implementation within BLIS,

we also show the corresponding real domain gemm implementation and the cgemm or
zgemm found in OpenBLAS and ARMPL.
In Figure A.1 (top-left), single-precision 1m and its corresponding real domain
benchmark track each other closely in the multithreaded configurations tested, as we
would have expected. Somewhat surprisingly, the vendor library, ARMPL, does not
appear to scale well at 56 threads (Figure A.1 (top-right)). Also somewhat surpris-
ingly, OpenBLAS performance is consistently low, even for sequential execution. This
suggests that while parallelism may be well-configured, their kernel is likely under-
performing.
Figure A.1 (bottom) tells a similar story of performance among double-precision
implementations, except that all BLIS implementations are, for reasons not immedi-
ately obvious, somewhat less efficient relative to peak performance than their single-
precision counterparts. ARMPL performance is more competitive for both one and 56
threads, though the single-core graph exposes evidence of a “crossover point” strat-
egy gone awry. ARMPL also seems to exhibit large swings in performance for certain
large, multithreaded problem sizes. Once again, OpenBLAS performance is much
 COMPLEX MATRIX MULTIPLICATION VIA THE 1M METHOD 27

lower, but consistently so.

In summary, BLIS’s 1m implementation performs extremely well on the Marvell
CN9975 when computing in single precision. Performance and scalability in double
precision, while not quite as impressive, is still highly competitive, especially when
compared to OpenBLAS and the ARM Performance Library.
Acknowledgements. We kindly thank Devangi Parikh for gathering the Thun-
derX2 performance results presented in Appendix A. We also thank the Texas Ad-
vanced Computing Center for providing access to the the Intel Xeon “Lonestar5”
compute cluster on which the performance data presented in Section 4 was gathered.

REFERENCES

[1] J. Dongarra, S. Hammarling, N. J. Higham, S. D. Relton, P. Valero-Lara, and

M. Zounon, The design and performance of batched BLAS on modern high-performance
computing systems, Procedia Computer Science, 108 (2017), pp. 495–504, [Link]
org/[Link] [Link]
article/pii/S1877050917307056.
[2] J. J. Dongarra, J. Du Croz, S. Hammarling, and I. Duff, A set of level 3 basic linear
algebra subprograms, ACM Trans. Math. Soft., 16 (1990), pp. 1–17.
[3] G. Frison, D. Kouzoupis, T. Sartor, A. Zanelli, and M. Diehl, BLASFEO: Basic lin-
ear algebra subroutines for embedded optimization, ACM Trans. Math. Soft., 44 (2018),
pp. 42:1–42:30, [Link] [Link]
[4] K. Goto and R. A. van de Geijn, Anatomy of high-performance matrix multiplication, ACM
Trans. Math. Soft., 34 (2008), pp. 12:1–12:25, [Link]
[Link]
[5] K. Goto and R. A. van de Geijn, High-performance implementation of the level-3 BLAS,
ACM Trans. Math. Soft., 35 (2008), pp. 4:1–4:14, [Link]
1377607, [Link]
[6] J. A. Gunnels, G. M. Henry, and R. A. van de Geijn, A family of high-performance matrix
multiplication algorithms, in Proceedings of the International Conference on Computational
Sciences-Part I, ICCS ’01, Berlin, Heidelberg, 2001, Springer-Verlag, pp. 51–60, [Link]
[Link]/[Link]?id=645455.653765.
[7] N. J. Higham, Stability of a method for multiplying complex matrices with three real matrix
multiplications, SIAM J. Matrix Anal. App., 13 (1992), pp. 681–687, [Link]
1137/0613043, [Link]
[8] J. Huang, Practical fast matrix multiplication algorithms, (2018). PhD thesis, The University
of Texas at Austin.
[9] J. Huang, D. A. Matthews, and R. A. van de Geijn, Strassen’s algorithm for tensor contrac-
tion, SIAM Journal on Scientific Computing, 40 (2018), pp. C305–C326, [Link]
10.1137/17M1135578, [Link] [Link]
//[Link]/10.1137/17M1135578.
[10] J. Huang, L. Rice, D. A. Matthews, and R. A. van de Geijn, Generating families of
practical fast matrix multiplication algorithms, in 31th IEEE International Parallel and
Distributed Processing Symposium (IPDPS 2017), May 2017, pp. 656–667, [Link]
org/10.1109/IPDPS.2017.56.
[11] J. Huang, T. M. Smith, G. M. Henry, and R. A. van de Geijn, Strassen’s algorithm reloaded,
in Proceedings of the International Conference for High Performance Computing, Network-
ing, Storage and Analysis, SC ’16, Piscataway, NJ, USA, 2016, IEEE Press, pp. 59:1–59:12,
[Link]
[12] J. Huang, C. D. Yu, and R. A. van de Geijn, Implementing Strassen’s algorithm with CUT-
LASS on NVIDIA Volta GPUs, FLAME Working Note #88, TR-18-08, The University
of Texas at Austin, Department of Computer Science, 2018, [Link]
apps/sites/default/files/tech reports/[Link].
[13] Intel, Math Kernel Library. [Link] 2019.
[14] Intel Corporation, Intel R 64 and IA-32 Architectures Optimization Reference Manual,
no. 248966-033, June 2016.
[15] Intel Corporation, Intel R Xeon R Processor E5 v3 Product Family: Processor Specification
Update, no. 330785-010US, September 2016.
[16] T. M. Low, F. D. Igual, T. M. Smith, and E. S. Quintana-Ortı́, Analytical modeling is
28 FIELD G. VAN ZEE

enough for high-performance BLIS, ACM Trans. Math. Soft., 43 (2016), pp. 12:1–12:18,
[Link] [Link]
[17] OpenBLAS. [Link] 2019.
[18] T. M. Smith, R. A. van de Geijn, M. Smelyanskiy, J. R. Hammond, and F. G. Van
Zee, Anatomy of high-performance many-threaded matrix multiplication, in Proceedings
of the 28th IEEE International Parallel & Distributed Processing Symposium (IPDPS),
IPDPS ’14, Washington, DC, USA, 2014, IEEE Computer Society, pp. 1049–1059, https:
//[Link]/10.1109/IPDPS.2014.110, [Link]
[19] F. G. Van Zee, Inducing complex matrix multiplication via the 1m method, FLAME Working
Note #85 TR-17-03, The University of Texas at Austin, Department of Computer Sciences,
2017.
[20] F. G. Van Zee, T. Smith, F. D. Igual, M. Smelyanskiy, X. Zhang, M. Kistler, V. Austel,
J. Gunnels, T. M. Low, B. Marker, L. Killough, and R. A. van de Geijn, The BLIS
framework: Experiments in portability, ACM Trans. Math. Soft., 42 (2016), pp. 12:1–12:19,
[Link]
[21] F. G. Van Zee and T. M. Smith, Implementing high-performance complex matrix multiplica-
tion via the 3m and 4m methods, ACM Trans. Math. Soft., 44 (2017), pp. 7:1–7:36.
[22] F. G. Van Zee and R. A. van de Geijn, BLIS: A framework for rapidly instantiating BLAS
functionality, ACM Trans. Math. Soft., 41 (2015), pp. 14:1–14:33, [Link]
1145/2764454.
[23] R. C. Whaley, A. Petitet, and J. J. Dongarra, Automated empirical optimization of soft-
ware and the ATLAS project, Parallel Computing, 27 (2001), pp. 3–35, [Link]
[Link] [Link]
article/pii/S0167819100000879. New Trends in High Performance Computing.
[24] C. D. Yu, J. Huang, W. Austin, B. Xiao, and G. Biros, Performance optimization for
the k-nearest neighbors kernel on x86 architectures, in Proceedings of the International
Conference for High Performance Computing, Networking, Storage and Analysis, SC ’15,
New York, NY, USA, 2015, ACM, pp. 7:1–7:12, [Link]
[Link]