Tí ulo a ículo / Tí ol a icle:
De i ing dense linea algeb a lib a ies
Au o es / Au o s
Bien inesi, Paolo ; Gunnels, John A. ; Mye s,
Ma ga e E. ; Quin ana O í, En ique S. ; Rhodes,
Tyle ; Van de Geijn, Robe A. ; Van Zee, Field G.
Re is a:
Fo mal Aspec s o Compu ing,
Ve sión / Ve sió:
Ve sió pos -p in
Ci a bibliog á ica / Ci a
bibliog à ica (ISO 690):
BIENTINESI, Paolo, e al. De i ing dense linea
algeb a lib a ies. Fo mal Aspec s o Compu ing,
2013, ol. 25, no 6, p. 933-945.
u l Reposi o i UJI:
h p://hdl.handle.ne /10234/113024
Unde conside a ion o publica ion in Fo mal Aspec s o Compu ing
De i ing Dense Linea Algeb a
Lib a ies
Paolo Bien inesi1, John Gunnels2, Ma ga e Mye s3, En ique Quin ana-O ´ı4,
Tyle Rhodes3, Robe an de Geijn3, and Field G. Van Zee3
1RWTH Aachen Uni e si y, Aachen, Ge many,
2IBM T.J. Wa son Resea ch Cen e , Yo k own Heigh s, NY,
3The Uni e si y o Texas a Aus in, Aus in, TX,
4Uni e sidad Jaime I, Cas ell´on, Spain
Abs ac . S a ing in he la e 1960s compu e scien is s including Dijks a and Hoa e ad oca ed goal-
o ien ed p og amming and he o mal de i a ion o algo i hms. The chie impedimen o ealizing his o
loop-based p og ams was ha a p io i de e mina ion o loop-in a ian s, a p e equisi e o de eloping loops,
was a ask oo complex o any bu he simples o ope a ions. A ound 2000, hese echniques we e o he
i s ime success ully applied o he domain o high-pe o mance dense linea algeb a lib a ies. This has led
o a mul i ude o pape s, mos ly published in he ACM T ansac ions o Ma hema ical So wa e, a sys em
o he mechanical de i a ion o algo i hms, and a high-pe o mance linea algeb a lib a y, lib lame, ha
includes mo e han a housand a ian s o algo i hms o mo e han a hund ed linea algeb a ope a ions.
To ou knowledge, his success s o y has un olded wi h limi ed awa eness on he pa he o mal me hods
communi y. This pape epo s on en yea s o expe ience and is mean o aise ha awa eness.
Keywo ds: Fo mal de i a ion; Linea algeb a lib a ies; Scien i ic compu ing
1. In oduc ion
Linea algeb a lib a ies eside a he bo om o he scien i ic compu ing ood chain. While mos p ac ical
applica ions gi e ise o spa se linea algeb a p oblems, a signi ican numbe o hem spends mos compu a-
ional ime sol ing dense ma ix p oblems. E en spa se linea algeb a p oblems o en ha e dense subp oblems
o be sol ed. As a esul , LAPACK [ABB+99], a package o dense ma ix ope a ions, de eloped in he la e
1980s and ea ly 1990s, is undoub edly he mos commonly used lib a y in his ield.
Since 2000, he FLAME p ojec a The Uni e si y o Texas a Aus in, Uni e sidad Jaume I (Spain),
and RWTH Aachen Uni e si y (Ge many) has been pu suing a eplacemen o LAPACK, lib lame [Zee09].
The domain poses a ew in e es ing challenges: scien i ic compu ing ends o exploi he la es a chi ec u es,
Co espondence and o p in eques s o: Robe an de Geijn, Depa men o Compu e Science, The Uni e si y o Texas a
Aus in, 1 Uni e si y S a ion, Aus in, TX 78712, USA. e-mail: [email p o ec ed]
2 P. Bien inesi e al.
demanding he highes possible pe o mance. This equi es he design o loop-based algo i hms ha cas mos
compu a ion in e ms o high-pe o ming ma ix-ma ix ope a ions (like ma ix-ma ix mul iplica ion). The
loop s eps h ough ma ices wi h block sizes chosen so as o op imize he euse o da a in caches. Fo a
speci ic ope a ion, he e a e o en mul iple algo i hmic a ian s, wi h one algo i hmic a ian ma ching a
gi en a chi ec u e be e han he o he s, yielding highe pe o mance. Thus, a lib a y should inco po a e
mul iple loop-based algo i hms o a gi en ope a ion so ha he bes one can be chosen. The LAPACK lib a y
does no include such a wide a ie y o algo i hms. This b ings up he ques ion o how o sys ema ically
ind algo i hmic a ian s o a gi en ope a ion. The o mal de i a ion o loop-based algo i hms u ns ou o
be he answe [Bie06, BGM+05, BGdG, Gun01, GGH dG01, G dG01, QO dG03, dGQO08], as we will
illus a e.
This pape does no p o ide a schola ly ea men o he ield o de i a ion o algo i hms. All we ha e
e e needed o de elop he desc ibed echniques is gi en in he ex by G ies [G i81], which i sel is based
on he wo ks o Dijks a [Dij68, Dij76] and Hoa e [Hoa69]. Wha he pape does p o ide is wha we belie e
o be an excellen p ac ical example o he applica ion o o mal de i a ion o loops o he domain o dense
linea algeb a.
2. De i a ion o Linea Algeb a Algo i hms
In his sec ion, we walk he eade h ough he de i a ion p ocess. The p ocedu e is comple ely ou ine o
us, as we ha e applied i o mo e han a hund ed ope a ions, yielding mo e han a housand ou ines ha
a e pa o he lib lame lib a y. We use he solu ion o he iangula Lyapuno equa ion as ou mo i a ing
example. A eade who is no well- e sed in linea algeb a should no wo y: he me hology is sys ema ic o
he poin whe e one does no need o be an expe in o de o apply i . A mo e basic ea men ha a ge s
no ices and has been used a he unde g adua e le el can be ound in [ dGQO08].
2.1. The FLAME me hodology
A undamen al insigh in ou p ojec was he ealiza ion ha he Fundamen al In a iance Theo em [G i81],
used o p o e he co ec ness o a loop in a p og am, can be o mula ed as a wo kshee ha is sys ema ically
illed ou , i s wi h asse ions (p edica es) and hen wi h commands (impe a i e s a emen s) [BGM+05].
The wo kshee , gi en in Figu e 1, will be illed ou wi h p edica es ha indica e p esc ibed s a es and
commands ha achie e hose s a es. I is illed ou in he o de indica ed in he column ma ked by S ep.
The p edica es Pp e ,Ppos ,Pin , and G ep esen he p econdi ion, pos condi ion, loop-in a an , and loop-
gua d, espec i ely. The loop-in a ian has o be ue in ou di e en places: be o e and a e he loop, and
a he op and bo om o he loop body. The o he pa s o he wo kshee will become ob ious as we ill i
ou o ou example.
2.2. Example: he solu ion o he iangula Lyapuno equa ion
We now show how he me hodology is applied o a p o o ypical example: he solu ion o he iangula
Lyapuno equa ion gi en by UTX+XU +C= 0 (o , al e na i ely, UTX+XU =−C), whe e Uis an
uppe iangula ma ix and Cand Xa e symme ic ma ices. He e he supe sc ip “T” indica es ma ix
ansposi ion. The solu ion Xis o be compu ed. Because o symme y, only he uppe iangula pa o Cis
s o ed and is o e w i en wi h he uppe iangula pa o X. This ope a ion is p eceded by a p e-p ocessing
phase (no discussed in his pape ) ha ans o ms he gene al (non- iangula ) ime-con inuous Lyapuno
equa ion (an ope a ion encoun e ed in con ol heo y [Kha02]) o he gi en iangula Lyapuno o m.
2.3. Filling ou he wo kshee
A his poin , he eade should imagine he wo kshee in Figu e 2 as being emp y and he s eps de ailed
below as illing ou he wo kshee in he indica ed o de .
S ep 1: The p econdi ion and pos condi ion.
De i ing Dense Linea Algeb a Lib a ies 3
S ep Anno a ed Algo i hm: [D,E,F,...] := op(A, B, C, D, . . .)
1a {Pp e }
4Pa i ion
whe e
2{Pin }
3while Gdo
2,3 {(Pin )∧(G)}
5a Repa i ion
whe e
6{Pbe o e }
8SU
5b Con inue wi h
7{Pa e }
2{Pin }
endwhile
2,3 {(Pin )∧ ¬ (G)}
1b {Ppos }
Fig. 1. Blank FLAME wo kshee used o de i e algo i hms.
The p econdi ion is gi en by C=ˆ
Cwhile he pos condi ion is C=X∧UTX+XU =−ˆ
C. He e
heˆ is needed o be able o eason abou he cu en con en s (s a e) o a iable C ela i e o he ini ial
con en s, ˆ
C. Fo b e i y, p ope ies o he ma ices (like iangula s uc u e and sizes) a e no exp essed in
he p econdi ion no in o he p edica es.
S ep 2: De i ing he loop-in a ian s. A undamen al insigh is ha many algo i hms sweep h ough
ma ices (a ays) in a sys ema ic and p edic able ashion. In his example, all ma ices a e ei he iangula
o symme ic and a e pa i ioned in o quad an s since his exposes egions o he ma ices ha a e ei he
(implici ly o explici ly) ze o, o he iangula ma ix, o no s o ed, o he symme ic ma ices:
U= UT L UT R
0UBR !, X = XT L XT R
⋆ XBR !, C = CT L CT R
⋆ CBR !,and ˆ
C= ˆ
CT L ˆ
CT R
⋆ˆ
CBR !,(1)
whe e UT L,XT L,CT L, and ˆ
CT L a e all con o mal (o he same size) and squa e. He e T L,T R,BL, and
BR s and o “ op-le ”, “ op- igh ”, “bo om-le ”, and “bo om- igh ”, espec i ely. The 0 and ⋆indica e
subma ices ha a e en i ely ze o o no s o ed, espec i ely. As he algo i hm p og esses, he op-le (TL)
quad an s will g ow om emp y (0 ×0) o encompassing he en i e ma ix.
Subs i u ing he pa i ioned ma ices in (1) in o he pos condi ion yields he exp ession in Figu e 3,
which we call he Pa i ioned Ma ix Exp ession [Bie06, dGQO08] (PME). I is a ecu si e de ini ion o he
ope a ion in e ms o he exposed subma ices.
The PME exp esses he compu a ion o be pe o med (which mus make he pos condi ion ue) in e ms
o he quad an s. Obse e ha as long as he loop has no inished, only pa o he compu a ion exp essed
by he PME is sa is ied: o come up wi h po en ial loop-in a ian s, one dele es some o he subexp essions
in he PME, as illus a ed in Figu e 4. As long as he exp ession ha is le is a alid exp ession ha has
he co ec size, i is a candida e. Some po en ial loop-in a ian s a e such ha subsequen s eps canno be
pe o med, which means ha hey do no yield an (admissible) algo i hm. Fo example, i only he o iginal
con en s o he ma ix a e le a e dele ing subexp essions om he PME, he loop can clea ly no comple e
4 P. Bien inesi e al.
S ep Anno a ed Algo i hm: C:= lyap unb(U, C)
1a nC=ˆ
Co
4Pa i ion U→ UT L UT R
0UBR !,X→ XT L XT R
⋆ XBR !,C→ CT L CT R
⋆ CBR !
whe e UT L is 0 ×0, XT L is 0 ×0, CT L is 0 ×0, ˆ
CT L is 0 ×0
28
>
<
>
: CT L CT R
⋆ CBR != XT L XT R
⋆ XBR !∧8
>
<
>
:
UT
T LXT L +XT LUT L =−ˆ
CT L
UT
T LXT R +XT R UBR =−ˆ
CT R −XT LUT R
XBR =ˆ
CBR
9
>
=
>
;
3while m(UT L )< m(U)do
2,3 8
>
>
>
<
>
>
>
:
0
B
@ CT L CT R
⋆ CBR != XT L XT R
⋆ XBR !∧8
>
<
>
:
UT
T LXT L +XT LUT L =−ˆ
CT L
UT
T LXT R +XT R UBR =−ˆ
CT R −XT LUT R
XBR =ˆ
CBR
1
C
A
∧(m(UT L )< m(U))
9
>
>
>
=
>
>
>
;
5a Repa i ion
UT L UT R
0UBR !→0
B
@
U00 u01 U02
0υ11 uT
12
0 0 U22
1
C
A, XT L XT R
⋆ XBR !→0
B
@
X00 x01 X02
⋆ χ11 xT
12
⋆ ⋆ X22
1
C
A,· · ·
whe e υ11,χ11,γ11 a e scala s
68
>
>
>
<
>
>
>
:
0
B
@
C00 c01 C02
⋆ γ11 cT
12
⋆ ⋆ C22
1
C
A=0
B
@
X00 x01 X02
⋆ χ11 xT
12
⋆ ⋆ X22
1
C
A
∧8
>
>
>
<
>
>
>
:
UT
00X00 +X00U00 =−ˆ
C00
UT
00x01 +X00u01 +x01υ11 =−ˆc01
UT
00X02 +X00U02 +x01uT
12 +X02U22 =−ˆ
C02
χ11 = ˆγ11 ∧xT
12 = ˆcT
12 ∧X22 =ˆ
C22
9
>
>
>
=
>
>
>
;
8
γ11 := (−γ11 −2uT
01c01)/(2υ11 )
cT
12 := −cT
12 −uT
01C02 −cT
01U02 −γ11uT
12
Sol e υ11xT
12 +xT
12U22 =cT
12 o e w i ing cT
12 wi h xT
12
5b Con inue wi h
UT L UT R
0UBR !←0
B
@
U00 u01 U02
0υ11 uT
12
0 0 U22
1
C
A, XT L XT R
⋆ XBR !←0
B
@
X00 x01 X02
⋆ χ11 xT
12
⋆ ⋆ X22
1
C
A,· · ·
7
8
>
>
>
>
>
>
>
>
>
>
>
>
>
>
<
>
>
>
>
>
>
>
>
>
>
>
>
>
>
:
0
B
@
C00 c01 C02
⋆ γ11 cT
12
⋆ ⋆ C22
1
C
A=0
B
@
X00 x01 X02
⋆ χ11 xT
12
⋆ ⋆ X22
1
C
A
∧
8
>
>
>
>
>
>
>
<
>
>
>
>
>
>
>
:
UT
00X00 +X00U00 =−ˆ
C00
UT
00x01 +X00u01 +x01υ11 =−ˆc01
UT
00X02 +X00U02 +x01uT
12 +X02U22 =−ˆ
C02
2uT
01x01 + 2υ11χ11 =−ˆγ11
uT
01X02 +υ11xT
12 +xT
01U02 +χ11uT
12 +xT
12U22 =−ˆcT
12
X22 =ˆ
C22
9
>
>
>
>
>
>
>
>
>
>
>
>
>
>
=
>
>
>
>
>
>
>
>
>
>
>
>
>
>
;
28
>
<
>
: CT L CT R
⋆ CBR != XT L XT R
⋆ XBR !∧8
>
<
>
:
UT
T LXT L +XT LUT L =−ˆ
CT L
UT
T LXT R +XT R UBR =−ˆ
CT R −XT LUT R
XBR =ˆ
CBR
9
>
=
>
;
endwhile
2,3 8
>
>
>
<
>
>
>
:
0
B
@ CT L CT R
⋆ CBR != XT L XT R
⋆ XBR !∧8
>
<
>
:
UT
T LXT L +XT LUT L =−ˆ
CT L
UT
T LXT R +XT RUBR =−ˆ
CT R −XT LUT R
XBR =ˆ
CBR
1
C
A
∧¬ (m(UT L )< m(U))
9
>
>
>
=
>
>
>
;
1b ˘UTX+XU =−C¯
Fig. 2. Wo kshee o de i ing he unblocked algo i hm o sol ing he iangula Lyapuno equa ion co esponding o Loop-
in a ian 3.

De i ing Dense Linea Algeb a Lib a ies 5
Subs i u ing he pa i ioned ope ands in (1) in o UTX+XU =−ˆ
C( he pos condi ion), yields
CT L CT R
⋆ CBR != XT L XT R
⋆ XBR !
∧ UT
T L 0
UT
T R UT
BR ! XT L XT R
⋆ XBR !+ XT L XT R
⋆ XBR ! UT L UT R
0UBR != −ˆ
CT L −ˆ
CT R
⋆−ˆ
CBR !
which (by linea algeb a manipula ion) is equi alen o
CT L CT R
⋆ CBR != XT L XT R
⋆ XBR !∧




UT
T LXT L +XT LUT L =−ˆ
CT L
UT
T LXT R +XT RUBR =−ˆ
CT R −XT LUT R
UT
BRXBR +XBRUBR =−ˆ
CBR −(UT
T RXT R +XT
T RUT R)
Fig. 3. The Pa i ioned Ma ix Exp ession (PME) ( ecu si e de ini ion o he ope a ion) o o e w i ing Cwi h he solu ion
o he iangula Lyapuno equa ion.
Loop-in a ian 1:
CT L CT R
⋆ CBR != XT L XT R
⋆ XBR !∧




UT
T LXT L +XT LUT L =−ˆ
CT L
UT
T LXT R+XT RUBR = + ˆ
CT R−XT LUT R
UT
BRXBR+XBRUBR = + ˆ
CBR−(UT
T RXT R +XT
T RUT R)
Loop-in a ian 2:
CT L CT R
⋆ CBR != XT L XT R
⋆ XBR !∧




UT
T LXT L +XT LUT L =−ˆ
CT L
UT
T LXT R+XT RUBR =−ˆ
CT R −XT LUT R
UT
BRXBR+XBRUBR = + ˆ
CBR−(UT
T RXT R +XT
T RUT R)
Loop-in a ian 3:
CT L CT R
⋆ CBR != XT L XT R
⋆ XBR !∧




UT
T LXT L +XT LUT L =−ˆ
CT L
UT
T LXT R +XT RUBR =−ˆ
CT R −XT LUT R
UT
BRXBR+XBRUBR = + ˆ
CBR−(UT
T RXT R +XT
T RUT R)
Loop-in a ian 4:
CT L CT R
⋆ CBR != XT L XT R
⋆ XBR !∧




UT
T LXT L +XT LUT L =−ˆ
CT L
UT
T LXT R +XT RUBR =−ˆ
CT R −XT LUT R
UT
BRXBR+XBRUBR =−ˆ
CBR −(UT
T RXT R +XT
T RUT R)
Fig. 4. Loop-in a ian s o he iangula Lyapuno equa ion.
in a s a e whe e he pos condi ion holds; his exhibi s i sel when no loop-gua d can be ound in S ep 3.
In a ian s a e hus sys ema ically de i ed om he PME.
We will now ocus on one loop-in a ian (Loop-in a ian 3) as we ill ou he emainde o he wo kshee
in Figu e 2. The me hodology yields algo i hms co esponding o he o he loop-in a ian s in an analogous
ashion.
S ep 3: Loop gua d G.We know ha a e he loop comple es, {Pin ∧ ¬G}is ue. No commands exis s
be ween his and he pos condi ion {Ppos }. Thus, Gmus be chosen so ha
CT L CT R
⋆ CBR != XT L XT R
⋆ XBR !∧






UT
T LXT L +XT LUT L =−ˆ
CT L
UT
T LXT R +XT RUBR =−ˆ
CT R −XT LUT R
XBR =ˆ
CBR


∧ ¬G
implies UTX+XU =−C. This yields he (nonunique) choice o he loop-gua d: G= (m(UT L)< m(U)),
6 P. Bien inesi e al.
whe e m(·) e u ns he ow dimension o he a gumen . (An implici assump ion he e is ha ma ices UT L,
XT L,CT L, and ˆ
CT L a e always kep con o mal, i.e., o he same size and squa e.)
S ep 4: Ini ializa ion. The ini ializa ion is an indexing s ep: he ma ices a e pa i ioned as in (1). The
ac ha his mus place he a iables in a s a e whe e Pin holds dic a es he choice whe ein he op-le
quad an s a e 0 ×0 (emp y). The e could be o he choices o he ini ializa ion, bu hose would in a iably
in ol e pe o ming compu a ion wi h he ma ices, al e ing hei con en s.
S ep 5: Mo ing h ough he ma ices. In S eps 5a and 5b, subma ices a e exposed so ha we make
p og ess h ough he ma ices. He e, hick lines ha e seman ic meaning: a new ow and column a e exposed.
Upda es will happen in he loop body, and hen ha ow and column a e mo ed ac oss he hick line o
cap u e he a e sal h ough he ma ices.
In exposing subma ices, we use no a ional con en ions ha allow special p ope ies o he subma ices
o be easily ecognized: lowe case G eek le e s deno e scala s, lowe case Roman le e s deno e column
ec o s, and uppe case Roman le e deno e ma ices. Subma ices like uT
12 can be easily ecognized as
being pa o a ow and hence a ow ec o ( ansposed column ec o ).
The ac ha he op-le quad an is ini ially emp y (0×0) and mus e en ually en elop he en i e ma ix
dic a es how he algo i hm a e ses h ough he ma ices. The a e sal h ough he ma ix, oge he wi h
he ini e size o he ope ands, means ha he e is a na u al, mono onically dec easing loop-bound unc ion,
=n−m(UT L), ha is bounded below. This ensu es ha he loop e mina es.
S ep 6: S a e be o e he upda e. The commands in S ep 5a a e me ely indexing ope a ions. Since no
compu a ion occu s be ween he op o he loop and S ep 6 in he wo kshee , he s a e o he subma ices ha
a e exposed by S ep 5a can be de e mined by ex ual subs i u ion and linea algeb a manipula ion, as illus-
a ed in Figu e 5. This yields he s a e desc ibed by S ep 6. The in a ian oge he wi h he epa i ioning
in S ep 5a dic a es he p edica e in S ep 6.
S ep 7: S a e a e he upda e. Simila ly, he commands in S ep 5b a e me ely indexing ope a ions. Since
he in a ian mus again hold, he s a e in S ep 7 can be sys ema ically de i ed by ex ual subs i u ion o
he subma ices in S ep 5b in o he loop-in a ian and linea algeb a manipula ion, as illus a ed in Figu e 6.
The loop-in a ian oge he wi h he ede ini ion o he quad an s in S ep 5b dic a e he p edica e in S ep 7.
S ep 8: Upda e. The upda e in S ep 8 is now dic a ed by he s a e ha he a iables a e in a S ep 6 and
he s a e ha hey mus be in a S ep 7:
•C00 al eady con ains X00, he solu ion o UT
00X00 +X00U00 =−ˆ
C00, and hence is no upda ed.
•c01 al eady con ains x01, he solu ion o UT
00x01 +X00u01 +x01υ11 =−ˆc01, and hence is no upda ed.
•C02 al eady con ains X02, he solu ion o UT
00X02 +X00U02 +x01uT
12 +X02U22 =−ˆ
C02, and hence is no
upda ed.
•γ11 con ains χ11 = ˆγ11 and needs o be o e w i en by he solu ion, χ11, o 2uT
01x01 + 2υ11χ11 =−ˆγ11.
Recognizing ha a his poin c01 con ains x01 and γ11 con ains ˆγ11, his can be accomplished by upda ing
γ11 wi h γ11 := (−γ11 −2uT
01c01)/(2υ11), whe e := is used o assignmen .
•cT
12 holds xT
12 = ˆcT
12 and needs o be o e w i en wi h he solu ion, xT
12, o uT
01X02 +υ11xT
12 +xT
01U02 +
χ11uT
12 +xT
12U22 =−ˆcT
12. Recognizing ha X02 has o e w i en C02, e c., his can be accomplished by
upda ing i wi h he solu ion, xT
12, o υ11xT
12 +xT
12U22 =−cT
12 −uT
01C02 −cT
01U02 −γ11uT
12:
–Fi s , we upda e cT
12 := −cT
12 −uT
01C02 −cT
01U02 −γ11uT
12.
–Nex , we compu e cT
12 := cT
12 (υ11I+U22)−1. This equi es he solu ion o a iangula sys em o
equa ions, since i is equi alen o compu ing he solu ion o (υ11I+U22)Tx12 =c12 and o e w i ing
cT
12 wi h xT
12.
We ei e a e ha he s a e in S ep 6 and he desi ed s a e in S ep 7 dic a e how he a iables (subma ices
o C) mus be upda ed.
Resul ing algo i hm. The esul ing algo i hm, s ipped o he anno a ions ha we e used o de i e i , is
gi en in Figu e 7 (le ), execu ing only he commands indica ed unde Va ian 3.
De i ing Dense Linea Algeb a Lib a ies 7
S ep 5a, gi en by
UT L UT R
0UBR !→


U00 u01 U02
0υ11 uT
12
0 0 U22


, XT L XT R
⋆ XBR !→


X00 x01 X02
⋆ χ11 xT
12
⋆ ⋆ X22


,
CT L CT R
⋆ CBR !→


C00 c01 C02
⋆ γ11 cT
12
⋆ ⋆ C22


, ˆ
CT L ˆ
CT R
⋆ˆ
CBR !→


ˆ
C00 ˆc01 ˆ
C02
⋆ˆγ11 ˆcT
12
⋆ ⋆ ˆ
C22


,
exp esses



UT L →U00 UT R →u01 U02 
0UBR → υ11 uT
12
0U22 !

,


XT L →X00 XT R →x01 X02 
⋆ XBR → χ11 xT
12
⋆ X22 !

,



CT L →C00 CT R →c01 C02 
⋆ CBR → γ11 cT
12
⋆ C22 !

,



ˆ
CT L →ˆ
C00 ˆ
CT R →ˆc01 ˆ
C02 
⋆ˆ
CBR → ˆγ11 ˆcT
12
⋆ˆ
C22 !



.
Subs i u ing hese in o he s a e o he a iables a he op o he loop ( he in a ian ) gi en by
CT L CT R
⋆ CBR != XT L XT R
⋆ XBR !∧




UT
T LXT L +XT LUT L =−ˆ
CT L
UT
T LXT R +XT RUBR =−ˆ
CT R −XT LUT R
XBR =ˆ
CBR
yields he s a e o he exposed subma ices:



C00 c01 C02 
⋆
⋆! γ11 cT
12
⋆ C22 !

=


X00 x01 X02 
⋆
⋆! χ11 xT
12
⋆ X22 !


∧















UT
00X00 +X00U00 =−ˆ
C00
UT
00 x01 X02 +x01 X02  υ11 uT
12
0U22 !=−ˆc01 ˆ
C02 −X00 u01 U02 
χ11 xT
12
⋆ X22 != ˆγ11 ˆcT
12
⋆ˆ
C22 !.
Algeb aic manipula ion o he abo e exp ession yields S ep 6 in he wo kshee .
Fig. 5. Sys ema ic de i a ion o he s a e o he a iables a S ep 6 in Figu e 2.
2.4. O he algo i hms
O he algo i hms a e de i ed om he o he loop-in a ian s. In addi ion, blocked algo i hms, which cas
mos compu a ion in e ms o ma ix-ma ix ope a ions and hence can a ain highe pe o mance, can be
de i ed by mo ing h ough he ma ix se e al ows and columns a a ime. Resul ing algo i hms a e gi en
in Figu e 7. In he blocked algo i hms, he ope a ion
Sol e UT
II XIJ +XIJ UJJ =CIJ
equi es he solu ion o he Syl es e equa ion. Algo i hms o ha ope a ion we e de i ed in [QO dG03].
8 P. Bien inesi e al.
Subs i u ing he ede ini ion o quad an s in S ep 5b
UT L UT R
0UBR !←


U00 u01 U02
0υ11 uT
12
0 0 U22


, XT L XT R
⋆ XBR !←


X00 x01 X02
⋆ χ11 xT
12
⋆ ⋆ X22


,
CT L CT R
⋆ CBR !←


C00 c01 C02
⋆ γ11 cT
12
⋆ ⋆ C22


, ˆ
CT L ˆ
CT R
⋆ˆ
CBR !←


ˆ
C00 ˆc01 ˆ
C02
⋆ˆγ11 ˆcT
12
⋆ ⋆ ˆ
C22



in o he desi ed s a e o he a iables a he bo om o he loop ( he in a ian )
CT L CT R
⋆ CBR != XT L XT R
⋆ XBR !∧




UT
T LXT L +XT LUT L =−ˆ
CT L
UT
T LXT R +XT RUBR =−ˆ
CT R −XT LUT R
XBR =ˆ
CBR
esul s in



C00 c01 C02
⋆ γ11 cT
12
⋆ ⋆ C22


=


X00 x01 X02
⋆ χ11 xT
12
⋆ ⋆ X22



∧

















U00 u01
0υ11 !T X00 x01
⋆ χ11 !+ X00 x01
⋆ χ11 ! U00 u01
0υ11 !=− ˆ
C00 ˆc01
⋆ˆγ11 !
U00 u01
0υ11 !T X02
xT
12 !+ X02
xT
12 !U22 =− ˆ
C02
ˆcT
12 !− X00 x01
⋆ χ11 ! U02
uT
12 !
X22 =ˆ
C22.
Algeb aic manipula ion yields he exp ession in S ep 7.
Fig. 6. Sys ema ic de i a ion o he s a e o he a iables a S ep 7 in Figu e 2.
2.5. Discussion
I is he no a ion we use ha acili a es he de i a ion p ocess: by p esen ing subma ices a he han
index anges, he de i ed algo i hm a oids much o he indexing clu e (bo h in he p esen a ion o he
algo i hm and he de ails o he p oo o co ec ness) ha is ypically ound in con en ional loop-based
algo i hms. Indeed, he algo i hm exposes only one loop e en hough i equi es app oxima ely n3 loa ing
poin ope a ions. The o he loops a e hidden inside o he linea algeb a ope a ions ha o m he body o
he loop. Algo i hms o hese ope a ions hemsel es can be, and ha e been, o mally de i ed, using he same
echniques as desc ibed in his pape .
2.6. F om algo i hm o code
Using a co ec algo i hm as he basis o you implemen a ion does no gua an ee ha he esul ing code
will be co ec . To p ese e he co ec ness o he algo i hm as we ansla e i o code, we de ined APIs o
di e en languages such ha he code closely esembles he algo i hm [BQO dG05, dGQO08, VZC dG+09].
An example o his is gi en in Figu e 8. In ha igu e, unblocked Va ian 3 is coded in M-sc ip , he sc ip ing
language o Ma lab [MLB87].
We c ea ed a hand ul o ou ines ha pa i ion and epa i ion he ma ices. Whi e-space is used o make
he code esemble he algo i hm as closely as possible. The code in Figu e 8 ga e he co ec answe he
i s ime i was un. We ha e de eloped APIs o he C p og amming language [BQO dG05, dGQO08,