Hacia una biblio eca BLAS ealmen e po able
en e di e en es ipos de acele ado es
Edua do Rod iguez-Gu iez,1A u o Gonzalez-Esc ibano,2Diego R. Llanos 3
Resumen— Las u inas de ´algeb a lineal BLAS son
ampliamen e u ilizadas en aplicaciones cien ´ı icas de
odo ipo. Exis en implemen aciones espec´ı icamen e
op imizadas pa a di e en es ipos de pla a o mas de
c´ompu o incluyendo acele ado es. Po ejemplo, la im-
plemen aci´on con enida en la biblio eca In el MKL,
apa e de ejecu a se en CPUs, incluye e siones pa-
a Xeon Phi, mien as que la biblio eca cuBLAS es ´a
especialmen e dise˜nada pa a GPUs de NVIDIA.
Sin emba go, los mecanismos pa a ges iona la me-
mo ia u ilizada po las es uc u as de da os sob e las
que se ealiza el compu o son di e en es en cada imple-
men aci´on, as´ı como algunos mecanismos elacionados
con las llamadas y el paso de pa ´ame os.
En es e a ´ıculo p esen amos una in e az ´unica
pa a BLAS, in eg ada en un modelo de p og ama-
ci´on he e og´enea (Con olle s) que sopo a g upos de
n´ucleos de CPU, acele ado es Xeon Phi o GPUs de
NVIDIA de o ma anspa en e pa a el p og amado .
Con es a p opues a es posible cons ui p og amas
po ables basados en u inas BLAS, que se ejecu an
en di e en es ipos de acele ado es cambiando simple-
men e un pa ´ame o de inicializaci´on. Nues a p o-
pues a explo a in e namen e la biblio eca espec´ı ica
pa a cada ipo de disposi i o. Las di e encias en sus
in e aces y en los mecanismos ex e nos pa a ges iona
la memo ia de los disposi i os, minimizando ans e-
encias, son anspa en es pa a el p og amado . Los
esul ados expe imen ales mues an que nues a abs-
acci´on no in oduce p´e didas de endimien o signi-
ica i as.
Palab as cla e— BLAS, P og amaci´on he e og´enea,
Con olle s, Acele ado es, GPU, Xeon Phi, MIC, CU-
DA.
I. In oducci´
on
La especi icaci´on BLAS, muy conocida en el ´ambi-
o de la compu aci´on cien ´ı ica, con iene la de ini-
ci´on de es conjun os (denominados ni eles) de u i-
nas ma em´a icas de bajo ni el que pe mi en ealiza
ope aciones ec o - ec o , ec o -ma iz y ma iz-
ma iz, espec i amen e. Exis en cie as implemen-
aciones de BLAS que pe mi en ob ene un mejo
endimien o cuando se ejecu an sob e un ha dwa e
espec´ı ico. Algunos ejemplos de ´es as son cuBLAS[1],
cuyas unciones es ´an op imizadas pa a GPUs que
sopo en CUDA, y la implemen aci´on que se encuen-
a como pa e de In el MKL[2], que pe mi e acele a
su ejecuci´on bajo pla a o mas In el, en e las que se
incluyen los cop ocesado es Xeon Phi.
Uno de los p oblemas asociados al uso de cop oce-
sado es es ´a en el ele ado cos e de mo e los da os
necesa ios desde el hos al disposi i o. Pa a una se-
cuencia conc e a de llamadas a u inas BLAS, esco-
1Dp o. de In o m´a ica, Uni . Valladolid, e-mail:
[email p o ec ed]
2Dp o. de In o m´a ica, Uni . Valladolid, e-mail:
[email p o ec ed]
3Dp o. de In o m´a ica, Uni . Valladolid, e-mail:
[email p o ec ed]
ge los momen os en los que se sinc onizan y mue en
da os en e el hos y el disposi i o es impo an e. Los
mecanismos pa a ealiza es as ans e encias a ´ıan
en e las di e sas implemen aciones de BLAS, e in-
cluso den o de cada implemen aci´on exis en di e en-
es opciones con di e en es g ados de anspa encia y
esponsabilidad de ca a al p og amado . Las biblio-
ecas MKL de In el y cuBLAS/n BLAS de NVIDIA
p o een dis in os mecanismos de o loading, ya sea
o almen e au om´a ico, asis ido po el compilado
(Compile -Assis ed O loading, CAO), o con ans e-
encias p og amadas manualmen e. A mayo es, aun-
que la in e az de las unciones BLAS pueda pa ece
uni o me en e las di e en es implemen aciones, en
la p ´ac ica exis en peque˜nas di e encias que con iene
ene en cuen a. Po ejemplo, los mismos pa ´ame os
pueden apa ece po alo o po e e encia en di e-
en es biblio ecas, di e en es nomb es o cons an es se
usan pa a indica el mismo ipo ope aciones o da os,
e incluso algunas biblio ecas pueden ene pa ´ame-
os ex a espec´ı icos de la p opia implemen aci´on,
ales como una es uc u a de da os que ep esen a el
con ex o de la biblio eca en cuBLAS.
En abajos p e ios p esen amos el modelo deno-
minado Con olle [3], que pe mi e la implemen aci´on
de algo i mos conc e os en unidades de c´ompu o o
ke nels. Cada ke nel se iden i ica po un nomb e e in-
e az ´unicos, pe o se pueden decla a di e en es im-
plemen aciones del mismo ke nel. Desde una gen´e i-
ca que u iliza una abs acci´on que pe mi e ejecu a lo
en cualquie disposi i o, has a una implemen aci´on
espec´ı icamen e op imizada pa a cada ipo de pla a-
o ma, o incluso pa a cada amilia de a qui ec u a.
Un con olado se asocia a un disposi i o en su c ea-
ci´on. Pe mi e asocia es uc u as de da os a su con-
ex o e in oca la ejecuci´on de ke nels po nomb e.
El con olado escoge la implemen aci´on del ke nel
in ocado m´as adecuada pa a la pla a o ma o dispo-
si i o al que es ´a asociado. La p ime a p opues a so-
po aba como disposi i o an o un g upo de n´ucleos
de CPU, como GPUs de NVIDIA[3]. Pos e io men-
e se a˜nadi´o el sopo e pa a acele ado es In el MIC
(Xeon Phi)[4]. In e namen e los Con olle s u ilizan
OpenMP y CUDA pa a la coo dinaci´on de a eas.
En es e abajo se p esen a una capa so wa e
de abs acci´on en e di e sas implemen aciones de
BLAS, u ilizando como sopo e la biblio eca de Con-
olle s. Es a p opues a pe mi e al usua io p og a-
ma c´odigos con secuencias de llamadas a las u inas
de BLAS con una in e az uni o me y selecciona con
un ´unico pa ´ame o el ha dwa e des ino en el que se
ejecu a ´a el c´odigo. In e namen e, la biblio eca es-
coge ´a la mejo implemen aci´on disponible pa a el
ha dwa e seleccionado, u ilizando llamadas a biblio-
eca BLAS espec´ı icamen e op imizadas pa a dicho
ha dwa e, p o is as po el ab ican e o po e ce as
pa es. Nues a p opues a combina po abilidad de
c´odigo y de endimien o.
El p esen e a ´ıculo es ´a es uc u ado en a ias
pa es. La secci´on 2 p esen a algunos abajos e-
lacionados. En la secci´on 3 se epasan algunas ca ac-
e ´ıs icas de las he amien as u ilizadas. La secci´on 4
de alla las soluciones u ilizadas en nues a p opues a.
En la secci´on 5 se desc ibe el abajo expe imen al
pa a alida la p opues a. Finalmen e, en la secci´on
6 se discu en las conclusiones y el abajo u u o.
II. T abajo elacionado
La idea de ejecu a las unciones de BLAS en un
en o no he e og´eneo no es nue a. Po ejemplo, exis-
e una biblio eca BLAS pensada pa a se ejecu ada
de o ma pa alela denominada PBLAS, p esen ada
po Choi e al.[5] y que puede se desca gada des-
de[6]. Una a ian e de ´es a, llamada He e oPBLAS[7]
y c eada po Manumachu e al., es capaz de ejecu-
a los algo i mos de BLAS explo ando el pa alelis-
mo median e memo ia dis ibuida con He e oMPI en
cl´us e es he e og´eneos, lo que a˜nade la capacidad de
balanceo de ca ga. Sin emba go, He e oPBLAS no es
capaz de explo a las capacidades de c´ompu o de los
cop ocesado es conec ados a los nodos y ´unicamen e
u iliza los co es de CPU.
Una soluci´on al e na i a a la p esen ada aqu´ı es
el conjun o de biblio ecas MAGMA[8], que con iene
e siones op imizadas de las unciones de BLAS y
LAPACK, y es ´a disponible an o pa a CPUs, CU-
DA, OpenCL[9] y Xeon Phi[10]. La p incipal di e-
encia que p esen a la soluci´on desc i a en es e a-
bajo espec o a MAGMA es que no p e ende gene a
una implemen aci´on e icien e po cada cop ocesado
exis en e, sino que u iliza en su luga las biblio e-
cas p opo cionadas po las p opias compa˜n´ıas que
desa ollan el ha dwa e. De es a o ma se e i an los
p oblemas de i ados de ac ualiza pe manen emen e
los algo i mos de BLAS pa a cada nue a a qui ec u-
a[11], o de ene que inclui nue os algo i mos que
puedan apa ece pa a a qui ec u as an iguas.
Finalmen e, dado que es posible p og ama an o
pa a CPUs, GPUs y Xeon Phi u ilizando OpenCL,
ambi´en cabe ci a la biblio eca clBLAS[12], que o -
ma pa e de clMa hLib a ies, que a su ez es ´a siendo
desa ollada po AMD.
III. BLAS y Con olle s
A. Implemen aciones especializadas de BLAS
La especi icaci´on BLAS, muy conocida en el ´ambi-
o de la compu aci´on cien ´ı ica, con iene la de inici´on
de es conjun os (denominados ni eles) de u inas
ma em´a icas b´asicas que pe mi en ealiza ope acio-
nes ec o - ec o , ec o -ma iz y ma iz-ma iz, es-
pec i amen e. Exis en cie as implemen aciones de
BLAS que pe mi en ob ene un mejo endimien o
cuando se ejecu an sob e un ha dwa e espec´ı ico. Al-
gunos ejemplos de ´es as son cuBLAS[1], cuyas un-
ciones es ´an op imizadas pa a GPUs que sopo en
CUDA, y la implemen aci´on que se encuen a como
pa e de In el MKL[2], que pe mi e acele a su ejecu-
ci´on bajo pla a o mas In el, en e las que se incluyen
los cop ocesado es In el MIC Xeon Phi.
Cada u ina iene no malmen e cua o e siones
pa a cua o ipos de da os di e en es. El ipo de da-
os se e leja en el nomb e de la unci´on con una spa-
a eales en pun o lo an e de p ecisi´on simple (p.ej.
SGEMM); d, pa a eales en pun o lo an e de p eci-
si´on doble (DGEMM); c, pa a n´ume os complejos en
pun o lo an e de p ecisi´on simple (CGEMM) y inal-
men e z, pa a complejos en pun o lo an e de p eci-
si´on doble (ZGEMM). No malmen e se suele emplea
un as e isco *o alg´un o o s´ımbolo en luga de la le-
a co espondien e al ipo de da os cuando se ci a
el nomb e de una u ina de o ma gen´e ica.
En cada ni el de BLAS se de inen di e en es u i-
nas. En ni el 1 encon amos u inas que hacen ope-
aciones escala - ec o o ec o - ec o . Po ejemplo,
la u ina *AXPY, que mul iplica los elemen os de un
ec o po un escala , o I*AMAX, que pa a un ec o
de uel e el ´ındice de la celda con el alo m´as al o
del mismo. En ni el 2 encon amos u inas con ope-
aciones ec o -ma iz, como po ejemplo *GEMV,
que de uel e el esul ado de mul iplica una ma iz
po un ec o . Finalmen e, en ni el 3 hay ope aciones
ma iz-ma iz, como *GEMM que pe mi e ob ene el
esul ado de la mul iplicaci´on de dos ma ices.
En cada implemen aci´on especializada pa a cop o-
cesado es encon amos di e sos mecanismos de o -
loading, que son las ope aciones de mo imien o de
da os desde la je a qu´ıa de memo ia del hos al dis-
posi i o, o ice e sa. La biblio eca MKL, po ejem-
plo, pe mi e hace o loading au om´a ico[13] de cie -
as unciones a Xeon Phi. En es e caso, una llamada
a unci´on de BLAS implica de o ma anspa en e el
mo imien o de da os de en ada al cop ocesado y la
ecupe aci´on de los esul ados mo i´endolos al hos .
Es e mecanismo, si bien debe se habili ado po el
usua io median e la llamada a una unci´on o la de ini-
ci´on de una a iable de en o no (MKL MIC ENABLE=1),
su ac i aci´on pa a una llamada a BLAS es decisi´on
de la implemen aci´on. S´olo se ealiza cuando ´es a
juzga que el ama˜no de da os a p ocesa es lo su i-
cien emen e g ande como pa a que el mo imien o al
cop ocesado compense[14]. En caso con a io la u-
ina se ejecu a en el hos . Aunque es una uncionali-
dad au om´a ica y anspa en e, s´olo es ´a disponible
pa a pa a un conjun o muy educido de unciones de
BLAS (*GEMM, *TRSM, y *STRMM), las uncio-
nes de ac o izaci´on LU (*GETRF, *GETRFNPI),
QR (*POTRF) y Cholesky (*GEQRF) y la educ-
ci´on a o ma idiagonal (SYRDB) de LAPACK[15],
[16], [17]. El encadenamien o de secuencias de lla-
madas a es as unciones pueden p oduci m´ul iples
mo imien os de da os innecesa ios, ya que aunque
los esul ados de una u ina se ayan a u iliza como
en adas en o a pos e io , el o loading au om´a ico
siemp e ans ie e las en adas y salidas de cada u-
ina de o ma independien e.
O a opci´on consis e en el uso de In el Langua-
ge Ex ensions o O load (LEO). Se a a de unas
ano aciones del compilado de In el (p agmas) que
pe mi en ealiza comunicaciones expl´ıci as de da os
en e hos y disposi i o en los momen os deseados.
Es o pe mi e el mayo g ado de con ol al p og ama-
do , a cos a de un mayo es ue zo po su pa e pa a
de e mina los mo imien os ´op imos y codi ica los,
apa e de gene a un c´odigo menos po able.
Como pa e de las posibilidades que o ece LEO
es ´a el modo de o loading asis ido po el compilado
(CAO), en el que el p opio usua io es esponsable de
de ini qu´e unciones o g upos de ´es as se ejecu a ´an
en el cop ocesado y cu´ales en el hos , quedando los
de alles del mo imien o pa cialmen e ocul os. Es, po
an o, una soluci´on m´as gene al que la an e io . To-
das las unciones BLAS de MKL pueden se ejecu a-
das median e es a ´ecnica y pe mi e eu iliza da os
en e llamadas a u inas sin mo imien os innecesa-
ios Sin emba go, las decisiones son esponsabilidad
del p og amado , y los mo imien os se sinc onizan
con los comienzos o inales de u inas.
NVIDIA iene una soluci´on simila al o loading
au om´a ico de MKL, median e una biblio eca llama-
da n BLAS, que se ejecu a sob e su implemen aci´on
especializada de BLAS pa a sus GPUs, denomina-
da cuBLAS. Al igual que en el caso an e io , e i ica
si la ans e encia median e el bus PCI compensa el
iempo es imado de ejecuci´on. De nue o al igual que
en MKL, s´olo es ´a disponible pa a un conjun o edu-
cido de unciones (*GEMM, *SYRK, (C,Z)HERK,
*SYR2K, (C,Z)HER2K, *TRSM, *TRMM, *SYMM
y (C,Z)HEMM)[18]. Si no se u iliza es a opci´on, las
ans e encias de memo ia en e hos y GPU deben
se ges ionadas con las llamadas habi uales a la bi-
blio eca de CUDA. Es o es equi alen e a la u iliza-
ci´on de In el LEO en cuan o a mayo con ol, mayo
es ue zo de desa ollo y meno po abilidad.
En cuan o a di e encias de in e az podemos des-
aca que en la segunda e si´on del API de cuBLAS
se u iliza un nue o pa ´ame o handle que es espec´ı i-
co de es a implemen aci´on. Se a a de una es uc-
u a que con iene in o maci´on sob e el con ex o de
la biblio eca, y que debe se c eada, pasada como
p ime pa ´ame o en cada llamada a es a biblio e-
ca, y pos e io men e des uida. O as di e encias me-
no es apa ecen po ejemplo en la o ma de pasa
cie os pa ´ame os. Po ejemplo, los escala es al a,
be a, e c., se pasan como pun e o en odas las lla-
madas a cuBLAS, mien as que en In el MKL se
pasan po alo . O o p oblema eside en los di-
e en es alo es que cada implemen aci´on asigna a
las cons an es o ipos de da o que de inen ope acio-
nes, ales como po ejemplo cublasOpe a ion y
CBLAS TRANSPOSE, que especi ican si se debe ealiza
o no la anspues a de una de las ma ices de en ada
an es de ealiza la ope aci´on. Relacionado con ´es a,
algunas biblio ecas como CBLAS pe mi en de ini si
la ma iz de en ada es ´a en ow-majo o de (com´un
en lenguajes como C) o en column-majo o de (em-
pleado po p og amas en lenguaje Fo an).
B. Con olle s
La biblio eca de Con olle s[3], [4] in oduce una
en idad abs ac a, el con olado , que se asocia en su
c eaci´on a un disposi i o. Es e puede se una GPU,
una acele ado In el MIC Xeon Phi, o un g upo de
n´ucleos de CPU que se manejan de o ma anspa-
en e como un acele ado .
El con olado sopo a uncionalidades pa a aso-
cia /desasocia una es uc u a de da os a su con ex-
o, c ea /des ui es uc u as de da os den o de su
con ex o que no ienen e lejo en el hos , y lanza
ke nels ( u inas de abajo) en el disposi i o.
Las llamadas a ke nels eciben como pa ´ame os
alo es, o e e encias a es uc u as de da os de su
con ex o (asociadas o c eadas). La asociaci´on de una
es uc u a de da os al con ex o de un con olado
implica que an es de se u ilizada po un ke nel los
da os del hos se copia an a la imagen que se c ea en
el con ex o del disposi i o. En el momen o de desaso-
cia la, si la es uc u a ha sido modi icada en el dis-
posi i o, los da os ac uales se copian en la imagen
o iginal del hos . Las copias hacia el disposi i o pue-
den se eage olazy, ealiz´andose en el momen o de
la asociaci´on o e as´andose has a que son necesa ias.
El con olado es esponsable de ges iona la cola
de pe iciones de ejecuci´on de ke nels, y maneja las
ans e encias de memo ia de o ma anspa en e. In-
e namen e, los con olado es u ilizan las biblio ecas
o uncionalidades na i as que p o een los endedo es
de las pla a o mas, o las biblio ecas de e ce os m´as
e icien es pa a cada pla a o ma. Po an o p o een
de una capa de abs acci´on pa a desa olla p og a-
mas po ables, consiguiendo la m´axima e iciencia en
la ges i´on del disposi i o g acias a explo a los me-
canismos na i os de cada uno.
Pa a de ini ke nels, la biblio eca de con olado es
u iliza un sis ema de decla aci´on gen´e ica de in e az
donde se de alla el nomb e, n´ume o de pa ´ame os,
ipos, nomb es y ol de en ada/salida. Un ke nel con
un nomb e ´unico puede se decla ado a ias eces con
un iden i icado de pla a o ma. Exis e un iden i ica-
do gen´e ico que si e pa a decla a ke nels que se
pueden ejecu a en cualquie pla a o ma. El con o-
lado p o ee de un espacio de iden i icaci´on de´ındices
de h ead simila al de CUDA u OpenCL, pe o que
asocia los h eads en o den ow-majo , lo que pe -
mi e po abilidad en e disposi i os. Se p o ee de un
in e az de acceso a los elemen os de es uc u as de
da os ´unico. Po an o el p og amado puede de ini
ke nels sencillos y po ables. La especi icaci´on se ha-
ce en g ano ino. El con olado es el esponsable de
ag upa los h eads en bloques (pa a GPUs) o de ge-
ne a a eas de g ano g ueso que eco en el espacio
de´ındices con bucles (en OpenMP) pa a su ejecuci´on
e icien e.
Sin emba go, ambi´en se pe mi e decla a ke nels
espec´ı icos pa a pla a o mas y a qui ec u as, donde
se pueden inclui es uc u as y p imi i as p opias del
modelo de p og amaci´on in e no. Po ejemplo, los
ke nels decla ados pa a GPUs pueden u iliza sin-
axis y p imi i as de CUDA pa a decla a y usa
memo ia compa ida, sinc oniza los h eads del blo-
que, e c. Es o pe mi e i cons uyendo una biblio eca
de ke nels especializados que el con olado u iliza ´a
p io iz´andolos sob e los ke nels gen´e icos cuando la
a qui ec u a asociada sea la adecuada. Pa ´ame os
de lanzamien o ales como los ama˜nos de bloque y
g id pa a GPUs, se seleccionan au om´a icamen e en
unci´on de una decla aci´on de p opiedades del c´odi-
go suminis ada po el p og amado y un modelo de
p edicci´on in e no. En caso de no ene una ca ac e-
izaci´on espec´ı ica se u ilizan alo es po de ec o que
asegu an m´axima ocupaci´on de los co es de la GPU.
Pa a consegui una uni o midad en la ges i´on de
las es uc u as de da os y en el paso de pa ´ame os,
el modelo de con olado es se apoya en la biblio eca
Hi map. Es a biblio eca pe mi e de ini es uc u as
de da os con dominios de´ındices, ales como a ays o
g a os. Cada es uc u a se ep esen a con una es uc-
u a denominada Hi Tile que con iene me a-da os y
pun e os a la zona de memo ia con los da os. La bi-
blio eca incluye un in e az uni icado pa a la decla a-
ci´on de espacios de ´ındices, c eaci´on, des ucci´on, y
acceso a los da os de un Hi Tile. Los pa ´ame os de
los ke nels en el modelo de con olado es son siem-
p e ipos b´asicos pasados po alo , o Hi Tiles de
cualquie ipo base.
IV. P opues a de soluci´
on
En es a secci´on se de alla la soluci´on p opues a
pa a consegui una biblio eca po able y e icien e de
BLAS a a ´es del modelo de con olado es. La bi-
blio eca de con olado es nos acili a una o ma de
cons ui un in e az ´unico pa a ke nels con di e en-
es implemen aciones especializadas pa a di e en es
pla a o mas. El con olado es el esponsable de es-
coge la implemen aci´on adecuada pa a la a qui ec-
u a a la que se asocia en su cons ucci´on. En caso
de no exis i una ap opiada, el con olado u iliza ´a
la implemen aci´on gen´e ica no especializada.
Nues a p opues a se basa en ex ende el meca-
nismo de los con olado es pa a implemen a ke nels
con c´odigo que se ejecu a en el hos , pa a a˜nadi es
nue os ipos de implemen aciones de ke nels, uno pa-
a cada ipo de disposi i o sopo ado. Los ipos se
denominan libCPU,libXPhi ylibGPU.
Es os nue os ipos de ke nels se decla an como
cualquie o o ke nel, pe o usando es os nomb es co-
mo iden i icado es de a qui ec u a. En luga de im-
plemen a un c´odigo especializado que se lanza y eje-
cu a en el disposi i o, implemen an llamadas a un-
ciones de biblio ecas ex e nas. Po an o su c´odigo
se ejecu a en el hos , pe o la llamada a la biblio eca
na i a dispa a la ejecuci´on de una u ina en el dis-
posi i o. En el caso de g upos de CPUs (libCPU) o
In el MIC XeonPhi (libXPhi), las llamadas se ´an a
la biblio eca MKL. En el caso de GPUs (libGPU) las
llamadas se ´an a cuBlas.
Se ha desa ollado una colecci´on de ke nels. Uno
po cada unci´on de BLAS y ipo de da os, que
ac ´uan c´omo w appe s sob e las biblio ecas o igina-
les. Hemos denominado a es a biblio eca de ke nels
cuBLAS In el MKL O os
Biblio eca Hi C lBLAS
Hi Map
Con olle s
P og ama del Usua io
Fig. 1
Capas de abs acci´
on y he amien as u ilizadas en la
soluci´
on p opues a.
Hi C lBLAS ( e Fig. 1). En es a biblio eca, cada
unci´on de BLAS es un ke nel del modelo de con o-
lado es, con es implemen aciones, una pa a cada
ipo de disposi i o (CPUs, XeonPhi, GPU). El in e -
az p opues o pa a cada unci´on de BLAS es com´un
y uni o me. Den o de cada implemen aci´on se ges-
ionan los de alles necesa ios pa a hace la llamada
adecuada a la biblio eca co espondien e. Po ejem-
plo, las implemen aciones libGPU con ienen in e na-
men e el handle necesa io pa a las llamadas a cuBlas,
que se c ea y des uye con el con olado cuando es e
es ´a asociado a una GPU. O as di e encias espec o
de las API de las biblio ecas u ilizadas como base se
esuel en en la implemen aci´on de cada ke nel.
Toda la ges i´on de ans e encias de da os ya es ´a
esuel a en la biblio eca de con olado es, que u iliza
CUDA, o CAO seg´un el ipo de disposi i o asocia-
do al con olado . Las es uc u as de ipo Hi Tile de
la biblio eca Hi map, u ilizadas en los con olado es,
p o een ambi´en de un in e az uni icado pa a la ges-
i´on de las es uc u as de da os.
V. Es udio expe imen al
A con inuaci´on se p esen a un abajo expe imen-
al pa a comp oba que la implemen aci´on de la p o-
pues a de la biblio eca de ke nels no impone una pe-
nalizaci´on demasiado ele ada con espec o al iempo
que a da ´ıa una aplicaci´on implemen ada di ec a-
men e usando una de las biblio ecas o iginales.
Como caso de p ueba se ha seleccionado un c´odigo
que pe mi e calcula el esul ado de mul iplica una
ma iz Apo o a B(ambas cuad adas) epe idas
eces po la de echa. Dicho de o o modo, pe mi e
ob ene el esul ado de C=A×Bn. En es e caso se
u iliza una implemen aci´on ecu si a: Ck=Ck−1∗B,
pa a k∈ {1,2, ..., n}, donde C0=A. Algunas imple-
men aciones como In el MKL no pe mi en especi i-
ca el mismo a gumen o como pa ´ame o de en ada
y salida en la misma llamada, po lo que el esul ado
de mul iplica Ck×Bno puede almacena se de uel a
en Cpa a mul iplica de nue o po Bpo la de echa
en la siguien e i e aci´on. Es o obliga al uso de una
e ce a ma iz ese ada en la unidad de c´ompu o
donde almacena el esul ado pa cial. U ilizamos un
segundo ke nel pa a hace la copia de la ma iz in-
e media a la que se usa como p ime ope ando en
cada i e aci´on. Es e segundo ke nel no iene ca ga
compu acional apa e de la simple ans e encia de
memo ia, lo que pe mi e e mejo si las p opias lla-
/∗CAL BLAS XPHI. cu ∗/
// C a a c e i z a c i ´on d e l k e ne l : Pa ´a me os po d e e c o
CAL KERNEL GPU CHAR STATIC( scopy , 2 , de , de , de ) ;
// D ec la ac i ´on de l a implemen aci ´on pa a GPU
CAL KERNEL LIBGPU( scopy , 5 ,
IVAL, in , n ,
IN, Hi Tile loa ∗, x ,
IVAL, in , incx ,
IO, Hi Tile loa ∗, y ,
IVAL, in , incy )
{
cublasScopy(comm−>handleCUBLAS , n , x . kda a , incx , y . kda a , incy ) ;
}
/∗CAL BLAS XPHI. cpp ∗/
// D ec la ac i ´on de l a implemen aci ´on pa a XeonPhi
CAL KERNEL XPHI( scopy , 5 ,
IVAL, in , n ,
IN, Hi Tile loa ∗, x ,
IVAL, in , incx ,
IO, Hi Tile loa ∗, y ,
IVAL, in , incy )
{
cb l as sco p y (n , x . da a , incx , y . da a , incy ) ;
}
Fig. 2
Ejemplo de decla aci´
on de dos implemen aciones (pa a GPU-CUDA y pa a In el MIC XeonPhi), del ke nel scopy
de BLAS, en la biblio eca Hi C lBLAS.
madas a los ke nels suponen un o e head ap eciable.
La implemen aci´on pe mi e juga con dos pa ´ame-
os, el n´ume o de eces que se mul iplica la ma iz
B y el ama˜no de la ma iz. An es de comenza la
i e aci´on p incipal, las ma ices iniciales AyBse
copian en la memo ia de la unidad de c´ompu o.
La implemen aci´on del c´odigo pa a las p uebas se
ealiz´o en cua o e apas y dio como esul ado es
c´odigos di e en es adem´as de la p opia biblio eca de
ke nels. En una p ime a ase, la especi icaci´on de la
aplicaci´on se adujo a dos a ian es, una u ilizando
In el MKL y a˜nadiendo los p agmas necesa ios pa-
a ejecu a las unciones de BLAS en una XeonPhi
median e O loading asis ido po compilado y op i-
miz´andola pa a educi las ans e encias de da os,
y una segunda a ian e en la que se u iliza cuBLAS
y el compilado de NVIDIA pa a ealiza el mismo
p ocesamien o en una GPU habili ada pa a CUDA.
En ambos casos se ha e i icado que las unciones de
BLAS se ejecu an ealmen e en sus espec i os co-
p ocesado es y que en ning´un caso el compilado o
el sis ema de ejecuci´on oman la decisi´on de ejecu-
a la unci´on en el hos . Po ejemplo en el caso de la
Xeon Phi u ilizamos a iables de en o no, ales como
OFFLOAD REPORT=3, an es de la ejecuci´on del p og a-
ma[19]. Es as dos a ian es (ma ixABk XeonPhi,
ma ixABk CUDA, no u ilizan Con olle s y se em-
plean como e e encias.
Pos e io men e, se ealiz´o la e si´on que u iliza
con olado es, empezando po la c eaci´on de la capa
in e media de ke nels. Como ya se ha explicado, po
cada unci´on de BLAS se implemen a on dos ke nels
con el mismo nomb e, uno pa a libXPhi que u ili-
za su equi alen e en In el MKL y una segunda que
u iliza la unci´on co espondien e en cuBLAS. Todo
es o se ealiza cua o eces po cada g upo de uncio-
nes (eg. *GEMM →SGEMM, DGEMM, CGEMM,
ZGEMM), y los ke nels de CUDA y los de In el MKL
se c ean en dos iche os sepa ados pa a pode compi-
la cada uno de ellos con su compilado espec´ı ico. A
con inuaci´on se c ea on sendos iche os de cabece a,
y inalmen e el p og ama p incipal.
El aspec o del p og ama p incipal se mues a en
la Fig. 3. En ´el se pueden e las di e en es e apas
de las que suele cons a no malmen e el desa ollo
de p og amas usando Con olle s. En p ime luga
se inicializa el sis ema de con olado es y se decla an
las a iables, ales como ma ices y ec o es, que se
ayan a emplea , u ilizando Hi Tiles. Pa a las ma i-
ces AyBse decla a y ese a espacio de memo ia en
el hos en la misma llamada. Pa a la ma iz Cs´olo
se decla an sus ca ac e ´ıs icas pe o no se ese a me-
mo ia pa a los da os, ya que se u iliza ´a s´olo como
a iable empo al en el disposi i o.
El siguien e paso es decla a los con olado es que
deseemos con el ipo CALCn l. Median e la unci´on
de c eaci´on el con olado queda ´a asociado con una
unidad de c´ompu o (como po ejemplo una GPU,
una Xeon Phi espec´ı ica o un conjun o de uno o m´as
co es CPU). Dado que odas las unciones de con-
olado es eciben como pa ´ame o un CALCn l, es
posible selecciona din´amicamen e (en iempo de eje-
cuci´on) qu´e unciones se ´an ejecu adas po cada uni-
dad de c´ompu o (con olado ). Es o ambi´en pe mi e
que las mismas unciones puedan se ejecu adas po
di e en es con olado es, o bien que cada con olado
ejecu e un conjun o conc e o de unciones.
A con inuaci´on, en el c´odigo, se de ine una con igu-
aci´on que incluye el n´ume o de hilos (l´ogicos o i -
uales) que se u iliza ´an en la unidad de c´ompu o, y
se copian las ma ices desde la CPU a la misma usan-
do CAL Cn lA ach. Como se expuso an e io men e,

in main ( in a gc , cha ∗a g [ ] )
{
CAL Cn lIni ( 1 ) ;
// D ecl a ac i ´on de dominio de i n d i c e s de cada ec o / m a i z
H i T i l e l o a ma ixA , ma ixB , ma ixC ;
h i il e Do m ai n (&ma ixC , loa , 2 , ows , columns ) ;
// Rese a de esp acio ( s ol o pa a e c o e s / ma ices en CPU)
hi ileDomainAlloc(&ma ixA , loa , 2 , ows , columns ) ;
hi ileDomainAlloc(&ma ixB , loa , 2 , ows , columns ) ;
// I n i c i a l i z a da os en hos
...
// C ea con o l a d o
CALCn l commGPU;
CAL Cn lC ea e(&commGPU, CAL CNTRL LIBGPU, 0 ) ;
// N´ume o de h eads ( ea le s en GPU, i u a l e s en CPU)
CALTh ead h eads ;
CALTh eadIni ( h eads , 2 , ows , columns ) ;
// Asocia ma ices que e s i d i ´an an o en hos como en a c e l e a d o
CAL Cn lA ach(&commGPU, &ma ixA ) ;
CAL Cn lA ach(&commGPU, &ma ixB ) ;
// C ea ma ices que e s i d i ´an s ´o l o en e l a cele ado
CAL Cn lIn e nal(&commGPU, &ma ixC ) ;
...
// Lanzamien o de k e n e l s
in imes ;
o ( im es = 0 ; ime s < i n a l p o w e ; imes++)
{
CAL Cn lLaunch (commGPU, sgemm , h eads , 13 , & i l e o p n o a ,
& i l e o p n o a , ows , columns , columns ,
1 . 0 , &ma ixA , columns , &ma ixB , columns ,
0 . 0 , &ma ixC , ows ) ;
CAL Cn lLaunch (commGPU, scopy , h eads , 5 , ows∗columns ,
&ma ixC , 1 , &ma ixA , 1 ) ;
}
...
// L ibe a e s u c u a s
CAL Cn lDe ach(&commGPU, &ma ixA ) ;
CAL Cn lDe ach(&commGPU, &ma ixB ) ;
CAL Cn lDes oyIn e nal(&commGPU, &ma ixC ) ;
CAL Cn lDes oy(&commGPU) ;
h i i l e F e e ( ma ixA ) ;
h i i l e F e e ( ma ixB ) ;
}
Fig. 3
Ex ac o del p og ama p incipal del ejemplo C=A×Bku ilizando la biblio eca de con olado es y
Hi C lBLAS. En es e caso el con olado u ilizado es ´
a asociado a la p ime a GPU del sis ema.
puede se necesa io ese a espacio den o del dis-
posi i o pa a una es uc u a empo al, sin memo ia
ese ada en en la CPU. En el ejemplo se u iliza la
unci´on CAL Cn lIn e nal pa a ealiza es a ese -
a pa a la ma iz C. Seguidamen e, el c´odigo con-
iene un bucle que lanza los ke nels co espondien es
a SGEMM, pasando los pa ´ame os necesa ios. En
unci´on del ipo de disposi i o al que es ´e asociado
el CALCn l que se pasa como pa ´ame o, el ame-
wo k de con olado es elegi ´a en e la implemen a-
ci´on de cuBLAS o la de In el MKL pa a Xeon Phi.
Du an e las i e aciones no se p oduce mo imien o de
da os en e el cop ocesado y el hos , lo que aumen a
la e iciencia. Finalmen e, al ealiza la ope aci´on de
de ach de las ma ices AyB, los da os de la ma iz
esul ado almacenada en Ase copian de uel a a la
CPU. En el caso de B, el con olado de ec a que
no ha sido modi icada y no ealiza ans e encia de
memo ia a la CPU. La ma iz C ese ada s´olo en el
disposi i o se des uye. Se des uye el con olado en
uso y se libe a la memo ia usada en el hos .
Las p uebas con es e c´odigo ue on ejecu adas en
una m´aquina denominada ‘Chime a’, que con iene
las ca ac e ´ıs icas especi icadas en la Tabla I. El mo-
i o pa a su selecci´on, como se puede e en dicha
abla, es que la misma dispone an o de un cop o-
cesado Xeon Phi como de una GPU Ti an Black,
lo que pe mi e hace p uebas u ilizando cualquie a
de los dos disposi i os, asegu ando que el es o del
ha dwa e de la pla a o ma es exac amen e el mismo.
Cada expe imen o se epi e al menos cinco eces, ob-
se ´andose una al a es abilidad en las medidas y e-
p oducibilidad de los esul ados.
Las Fig. 4 y 5 mues an la elaci´on en e los iem-
pos de ejecuci´on del p og ama que ob iene la mul-
iplicaci´on de una ma iz po o a epe idas eces
(C=A×Bk), pa a un k= 3, en unci´on del ama˜no
de la ma iz que se le pasa como pa ´ame o.
La p ime a igu a, Fig. 4, mues a la di e encia en-
e la e si´on que hace llamadas a cuBLAS di ec a-
men e ( ep esen ada median e pun os cuad ados) y
la e si´on que u iliza llamadas a las unciones del
con olado que ac ´uan como w appe s de cada un-
ci´on de cuBLAS. En es a g ´a ica se puede e que el
uso de Con olle s (y po an o de la biblio eca que se
desc ibe en es e abajo). La di e encia de endimien-
o pa a la e si´on con Con olle s comienza con una
pe dida de has a un 79 % pa a ma ices peque˜nas de
1000 ×1000, que se educe ´apidamen e al aumen a
el ama˜no de la ma iz, has a algo menos de un 4 %
de penalizaci´on pa a ma ices de 20000 ×20000.
La segunda g ´a ica, Fig. 5 con iene la misma com-
pa aci´on pa a Xeon Phi. Una e si´on u iliza llamadas
di ec as a la implemen aci´on CBLAS de In el MKL,
y la o a u iliza el mismo c´odigo de con olado u ili-
zado en el caso de la GPU, pe o es a ez asociando el
con olado a la unidad Xeon Phi en el momen o de
su c eaci´on. Se puede e que en es e caso los iempos
de ejecuci´on de la a ian e con con olado , con ma-
ices g andes a pa i de 6500 ×6500, la des en aja
onda siemp e en e el 10 % y el 19 %. Sin emba go,
con ama˜nos peque˜nos de has a 6000×6000, el com-
po amien o depende mucho del ama˜no en conc e o.
En cie os casos se ob ienen penalizaciones de has a
el 19 %. En o os, sin emba go, la e si´on con con-
olado ob iene en ajas de endimien o de has a un
14 %, o incluso en un caso conc e o con una en aja
de un 46 %. Pues o que es os esul ados son consis-
en es a lo la go de a ias ejecuciones, ac ualmen e
se es ´a abajando en ealiza un p o iling de allado
pa a a e igua las azones de es as a iaciones seg´un
el ama˜no de las ma ices,
O a en aja que se ha podido demos a median e
la implemen aci´on de es e so wa e es que el uso de
las nue as biblio ecas de ke nels pe mi en que la apli-
caci´on pueda cambia en e las unciones de CBLAS
MKL y las de cuBLAS an s´olo cambiando un ´uni-
co pa ´ame o al c ea el Con olle , consiguiendo la
po abilidad deseada.
VI. Conclusiones y abajo u u o
En es e abajo se p esen a una implemen aci´on de
una biblio eca de unciones BLAS he e og´enea, de al
o ma que es posible de ini algo i mos gen´e icos in-
dependien es de la pla a o ma, u ilizando de o ma
anspa en e an o implemen aciones p opias como
las unciones BLAS que p o een las biblio ecas de
los p opios ab ican es de los di e en es acele ado es
y cop ocesado es. Es a soluci´on pe mi e po a ´acil-
men e c´odigos que p e iamen e ealiza an llamadas a
biblio ecas espec´ı icas BLAS de o ma muy sencilla,
a la ez que habili a el uso de m´ul iples cop ocesa-
do es de ipos di e en es sin obliga al p og amado
Sis ema ope a i o Cen OS 7.2.1511
P ocesado In el Xeon E5-2620 3
Velocidad de eloj 2.40GHz
Memo ia p incipal 32 GB DDR3 1333 MHz
Memo ia cach´e 15 MB
Cop ocesado In el Xeon Phi (KNC) A3120
Cop ocesado NVIDIA GTX Ti an Black GK110B
En o no In el Pa allel S . XE 2017.0.035
En o no NVIDIA CUDA (n cc) 8.0.44
TABLA I
Ca ac e ´
ıs icas ha dwa e y so wa e del en o no de
p uebas.
a esol e las di e encias en e di e en es in e aces o
duplica el c´odigo.
A mayo es, g acias al sis ema de con olado es, es-
a soluci´on e i a el mo imien o innecesa io de los da-
os que esiden en los di e en es ipos de unidades de
c´ompu o has a que los da os esul ado son expl´ıci a-
men e sinc onizados con el hos po el p og amado .
Al u iliza es e sis ema, el o e head in oducido con
espec o a una implemen aci´on di ec a del p og a-
ma p incipal u ilizando las biblio ecas p opo ciona-
das po los ab ican es es peque˜no.
Uno de los obje i os a co o plazo es la educci´on
del o e head que exis e en las ejecuciones pa a la e -
si´on libXPhi con ama˜nos de ma ices g andes, como
se ha podido ap ecia median e el caso de es udio
u ilizado. Adem´as, pa a pode ene una base expe-
imen al m´as amplia, se es ´an implemen ando y p o-
bando nue os ejemplos, incluyendo casos que pe mi-
an p oba di e sas op imizaciones o solapa c´alculo
y comunicaci´on. Se plan ea ambi´en como abajo
u u o compa a con o as lib e ´ıas de BLAS o ien-
adas a la po abilidad en sis emas he e og´eneos.
Ag adecimien os
Es e abajo ha sido pa cialmen e inanciado po el
Minis e io de Ciencia e Inno aci´on (MICINN) y po
el p og ama ERDF de la Uni´on Eu opea: P oyec o
HomP og-He Sys (TIN2014-58876-P), CAPAP-H6
(TIN2016-81840-REDT), y COST P og am Ac ion
IC1305: Ne wo k o Sus ainable Ul ascale Compu-
ing (NESUS).
Re e encias
[1] NVIDIA Co po a ion, “cuBLAS Lib a y: Use
Guide,” h p://docs.n idia.com/pd /CUBLAS_Lib a y.
pd , Jan. 2017.
[2] In el Co po a ion, “In el R
Ma h Ke nel Lib a y (In el R

MKL)”, h p://so wa e.in el.com/en-us/in el-mkl.
[3] Ana Mo e on-Fe nandez, Hec o O ega-A anz, and A -
u o Gonzalez-Esc ibano, “Con olle s: An abs ac ion o
ease he use o ha dwa e accele a o s,” The In e na ional
Jou nal o High Pe o mance Compu ing Applica ions,
p. 109434201770296, May 2017.
[4] Ana Mo e on-Fe nandez, Edua do Rod iguez-Gu iez, A -
u o Gonzalez-Esc ibano, and Diego R. Llanos, “Suppo -
ing he Xeon Phi cop ocesso in a He e ogeneous P o-
g amming Model,” in Eu o-Pa 2017: Pa allel P oces-
sing, San iago de Compos ela, Galicia, Spain., Aug. 2017,
Sp inge , Cham, In p ess.
[5] Jaeyoung Choi, Jack Donga a, Susan Os oucho , An-
oine Pe i e , Da id Walke , and R. Clin on Whaley, “A
0
2
4
6
8
10
12
14
16
0 5 10 15 20
Tiempo (segundos)
Ca dinalidad de cada dimensión (en miles de elemen os)
C=A×Bn usando cuBLAS en GTX Ti an Black (Chime a)
cuBLAS
cuBLAS + Con olle s
Fig. 4
Resul ados expe imen ales pa a el p og ama de mul iplicaci´
on C=A×Bnejecu ado en una GPU NVIDIA,
an o la e si´
on que ealiza llamadas di ec amen e a cuBLAS, como la e si´
on con Con olle s.
0
1
2
3
4
5
0 2 4 6 8 10
Tiempo (segundos)
Ca dinalidad de cada dimensión (en miles de elemen os)
C=A×Bn usando MKL CBLAS en Xeon Phi 3120 (Chime a)
MKL CBLAS (Xeon Phi)
MKL CBLAS (Xeon Phi) + Con olle s
Fig. 5
Resul ados expe imen ales pa a el p og ama de mul iplicaci´
on C=A×Bn, ejecu ado en un cop ocesado Xeon
Phi, an o la e si´
on que ealiza llamadas di ec amen e a CBLAS de MKL. como la e si´
on con Con olle s.
p oposal o a se o pa allel basic linea algeb a subp o-
g ams,” in Applied Pa allel Compu ing Compu a ions in
Physics, Chemis y and Enginee ing Science. Aug. 1995,
pp. 107–114, Sp inge , Be lin, Heidelbe g.
[6] “PBLAS Home Page,” h p://www.ne lib.o g/
scalapack/pblas_q e .h ml.
[7] Ra i Reddy Manumachu, Alexey Las o e sky, and Pe-
d o Alonso, “He e ogeneous PBLAS: Op imiza ion o
PBLAS o He e ogeneous Compu a ional Clus e s,” in
2008 In e na ional Symposium on Pa allel and Dis ibu-
ed Compu ing, July 2008, pp. 73–80.
[8] S animi e Tomo , Jack Donga a, and Ma c Baboulin,
“Towa ds dense linea algeb a o hyb id GPU accele a-
ed manyco e sys ems,” Pa allel Compu ing, ol. 36, no.
5, pp. 232–240, June 2010.
[9] Peng Du, Rick Webe , Pio Luszczek, S animi e Tomo ,
G ego y Pe e son, and Jack Donga a, “F om CUDA o
OpenCL: Towa ds a pe o mance-po able solu ion o
mul i-pla o m GPU p og amming,” Pa allel Compu-
ing, ol. 38, no. 8, pp. 391–407, Aug. 2012.
[10] Jack Donga a, Ma k Ga es, Azzam Haida , Yulu Jia,
Khai ul Kabi , Pio Luszczek, and S animi e To-
mo , “HPC P og amming on In el Many-In eg a ed-Co e
Ha dwa e wi h MAGMA Po o Xeon Phi,” Scien i ic
P og amming, ol. 2015, pp. e502593, Ap . 2015.
[11] Ma k Ga es, “MAGMA Fo um: Pe o mance is-
sue,” h p://icl.cs.u k.edu/magma/ o um/ iew opic.
php? =2& =1475, Dec. 2016.
[12] Ken Knox, Jian Liu, Da id Tanne , Pa an Yalaman-
chili, Ch is ian Kellne , Hugh Pe kins, Tingxing Tim
Dong, Ga¨e an Lehmann, Ced ic Nug e en, and Benja-
min Coquelle, “clBLAS: a so wa e lib a y con aining
BLAS unc ions w i en in OpenCL,” h ps://gi hub.
com/clMa hLib a ies/clBLAS, May 2017, o iginal-da e:
2013-08-13T15:05:53Z.
[13] In el Co po a ion, “Using In el R
MKL Au oma ic O -
load on In el R
Xeon Phi Cop ocesso s,” h ps://goo.
gl/1kq8GB, 2011.
[14] In el Co po a ion, “In el R
MKL Au oma ic O load
enabled unc ions o In el Xeon Phi cop ocesso s,”
h ps://goo.gl/9jV7PY, Feb. 2013.
[15] Ken Mil eld, “In el Xeon Phi MIC O load P og amming
Models,” h ps://goo.gl/ KxFgx, No . 2014.
[16] In el Co po a ion, “In el R
Ma h Ke nel Lib a y (In el R

MKL) 11.3 Release No es,” h ps://goo.gl/oPU5ay,
Feb. 2016.
[17] In el Co po a ion, “In el R
Ma h Ke nel Lib a y (In el R

MKL) 2017 Release No es,” h ps://goo.gl/82 Mkh,
Sep . 2016.
[18] NVIDIA Co po a ion, “NVBLAS,” h p://docs.
n idia.com/cuda/n blas/index.h ml.
[19] James Je e s, James Reinde s, and A inash Sodani, In-
el Xeon Phi P ocesso High Pe o mance P og amming,
Mo gan Kau mann, Camb idge, MA, edici´on: 2 edi ion,
June 2016.
[20] L. Susan Black o d, An oine Pe i e , Roldan Pozo, Ka-
in Reming on, R. Clin Whaley, James Demmel, Jack
Donga a, Iain Du , S en Hamma ling, G eg Hen y, Mi-
chael He oux, Linda Kau man, and And ew Lumsdai-
ne, “An upda ed se o basic linea algeb a subp og ams
(BLAS),” ACM T ansac ions on Ma hema ical So wa e,
ol. 28, no. 2, pp. 135–151, June 2002.
[21] Ga y W. Howell, James W. Demmel, Cha les T. Ful-
on, S en Hamma ling, and Ka en Ma mol, “Cache e i-
cien bidiagonaliza ion using BLAS 2.5 ope a o s,” ACM
T ansac ions on Ma hema ical So wa e, ol. 34, no. 3,
pp. 1–33, May 2008.
[22] Je emy G. Siek, Ian Ka lin, and E. R. Jessup, “Build o
o de linea algeb a ke nels,” Ap . 2008, pp. 1–8, IEEE.
[23] Ga y W. Howell and Cha les T. Ful on, “Cache E icien
Householde Bidiagonaliza ion,” The College o William
and Ma y, Williamsbu g, Vi ginia, June 2003.