批量厄米矩阵特征值分解的GPU算法_黄荣锋.pdf

h t t p:/ww w.j s j k x.c o mD O I:1 0.1 1 8 9 6/j s j k x.2 2 0 1 0 0 2 3 21)h t t p s:/d o c s.n v i d i a.c o m/c u d a/p d f/c u b l a s _L i b r a r y.p d f2)h t t p s:/g i t h u b.c o m/R O C m S o f t w a r e P l a t f o r m/r o c B L A S3)h t t p:/i c l.c s.u t k.e d u/m a g m a/4)h t t p s:/g i t h u b.c o m/e c r c/k b l a s-g p u到稿日期:2 0 2 2-0 1-2 5 返修日期:2 0 2 2-0 8-2 4基金项目:国家重点研发计划(2 0 1 7 Y F B 0 2 0 2 2 0 2);中国科学院战略性先导科技专项(X D C 0 5 0 0 0 0 0 0)T h i sw o r kw a ss u p p o r t e db yt h eN a t i o n a lK e yR e s e a r c ha n dD e v e l o p m e n tP r o g r a mo fC h i n a(2 0 1 7 Y F B 0 2 0 2 2 0 2)a n dS t r a t e g i cP r i o r i t yR e s e a r c hP r o g r a mo fC h i n e s eA c a d e m yo fS c i e n c e s(X D C 0 5 0 0 0 0 0 0).通信作者:赵永华(y h z h a o s c c a s.c n)批量厄米矩阵特征值分解的G P U算法黄荣锋1,2刘世芳1,2赵永华11中国科学院计算机网络信息中心 北京1 0 0 0 8 02中国科学院大学 北京1 0 0 0 8 0(h r f s c c a s.c n)摘 要 批量矩阵计算问题广泛存在于科学计算与工程应用领域。随着性能的快速提升,G P U已成为解决这类问题的重要工具之一。矩阵特征值分解属于双边分解,需要使用迭代算法进行求解,不同矩阵的迭代次数可能不同,因此,在G P U上设计批量矩阵的特征值分解算法比设计L U分解等单边分解算法更具挑战性。文中针对不同规模的矩阵,基于J a c o b i算法设计了相应的批量厄米矩阵特征值分解G P U算法。对于共享内存无法存储的矩阵,采用矩阵“块”操作技术提升计算强度,从而提高G P U的资源利用率。所提算法完全在G P U上运行,避免了C P U与G P U之间的通信。在算法实现上,通过k e r n e l融合减少了k e r n e l启动负载和全局内存访问。在V 1 0 0G P U上的实验结果表明,所提算法优于已有工作。R o o f l i n e性能分析模型表明,文中给出的实现已接近理论上限,达到了4.1 1 T F L O P S。关键词:厄米矩阵;特征值分解;批量计算;R o o f l i n e模型;性能分析中图法分类号 T P 3 0 1 B a t c h e dE i g e n v a l u eD e c o m p o s i t i o nA l g o r i t h m s f o rH e r m i t i a nM a t r i c e so nG P UHUAN GR o n g f e n g1,2,L I US h i f a n g1,2a n dZ HA OY o n g h u a11C o m p u t e rN e t w o r kI n f o r m a t i o nC e n t e r,C h i n e s eA c a d e m yo fS c i e n c e s,B e i j i n g1 0 0 0 8 0,C h i n a2U n i v e r s i t yo fC h i n e s eA c a d e m yo fS c i e n c e s,B e i j i n g1 0 0 0 8 0,C h i n a A b s t r a c t B a t c h e dm a t r i xc o m p u t i n gp r o b l e m sa r ew i d e l ye x i s t e d i ns c i e n t i f i cc o m p u t i n ga n de n g i n e e r i n ga p p l i c a t i o n s.W i t hr a p i dp e r f o r m a n c e i m p r o v e m e n t s,G P Uh a sb e c o m ea n i m p o r t a n t t o o l t os o l v es u c hp r o b l e m s.T h ee i g e n v a l u ed e c o m p o s i t i o nb e l o n g s t ot h et w o-s i d e dd e c o m p o s i t i o na n dm u s t b e s o l v e db y t h e i t e r a t i v e a l g o r i t h m.I t e r a t i v en u m b e r s f o r d i f f e r e n tm a t r i c e s c a nb ev a r i e d.T h e r e f o r e,d e s i g n i n ge i g e n v a l u ed e c o m p o s i t i o na l g o r i t h m sf o rb a t c h e dm a t r i c e so nt h eG P Ui sm o r ec h a l l e n g i n gt h a nd e s i g n i n gb a t c h e da l g o r i t h m s f o r t h eo n e-s i d e dd e c o m p o s i t i o n,s u c ha sL Ud e c o m p o s i t i o n.T h i sp a p e rp r o p o s e sb a t c h e da l g o r i t h m sb a s e do nt h eJ a c o b ia l g o r i t h m sf o re i g e n v a l u ed e c o m p o s i t i o no fH e r m i t i a n m a t r i c e s.F o rm a t r i c e st h a tc a n n o tr e s i d ei ns h a r e d m e m o r yw h o l l y,t h eb l o c kt e c h n i q u ei su s e dt oi m p r o v et h ea r i t h m e t i ci n t e n s i t y,t h u si m p r o v i n gt h eu s eo fG P Ur e s o u r c e s.A l g o r i t h m sp r e s e n t e d i nt h i sp a p e rr u nc o m p l e t e l yo nt h eG P U,a v o i d i n gt h ec o mm u n i c a t i o nb e t w e e nt h eC P Ua n dG P U.K e r n e l f u s i o ni sa d o p t e dt od e c r e a s e t h eo v e r h e a do f l a u n c h i n gk e r n e l a n dg l o b a lm e m o r ya c c e s s.E x p e r i m e n t a l r e s u l t so nV 1 0 0G P Us h o wt h a to u ra l g o r i t h m sa r eb e t t e r t h a ne x i s t i n gw o r k s.P e r f o r m a n c e e v a l u a t i o nr e s u l t so f t h eR o o f l i n em o d e l i n d i c a t e t h a t o u r i m p l e m e n t a-t i o n sa r ec l o s e t ot h eu p p e rb o u n d,a p p r o a c h i n g4.1 1 T F L O P S.K e y w o r d s H e r m i t i a nm a t r i x,E i g e n v a l u ed e c o m p o s i t i o n,B a t c hc o m p u t i n g,R o o f l i n em o d e l,P e r f o r m a n c ee v a l u a t i o n 1 引言批量矩阵计算问题广泛存在于机器学习、计算机视觉、天体物理学等领域1-4。对于多核C P U,这类问题的求解相对容易。例如结合O p e n MP和高度优化的L A P A C K/B L A S库(如MK L,o p e n B L A S)即可获得较好的性能,因为许多计算可以在C P U的高速缓存中进行。然 而,由 于 缓 存 较 小,在G P U上求解批量矩阵计算问题十分困难。B L A S库是科学计算领域的基础函数库之一,库中的矩阵乘(G EMM)函数占据着核心地位。为此,许多G P U供应商提供批量G EMM的高效实现1),2),开源软件包MAGMA3)和K B L A S4)也提供批量G EMM在G P U上的实现。L A P A-C K库中函数的批量算法获得了广泛关注。例如,D o n g等5研究了批量C h o l e s k y分解;A b d e l f a t t a h等6研究了部分选主元的批量L U分解;文献7 对批量矩阵计算问题的研究进行了较为全面的总结。与C h o l e s k y分解等单边分解不同,矩阵特征值分解属于双边分解,必须使用迭代算法进行求解。不同矩阵的迭代次数可能不同,因此进行批量计算的矩阵几乎不能同时收敛。当一些矩阵收敛后,计算量变小,剩余的矩阵难以充分利用G P U上的众多计算核心。因此,在G P U上设计批量矩阵的特征值分解算法更具挑战性。1)h t t p s:/d o c s.n v i d i a.c o m/c u d a/c u s o l v e r/i n d e x.h t m l特征值分解算法主要分为两类:三对角化方法和J a c o b i方法。三对角化方法主要包括三步:第一步使用H o u s e h o l d e r正交变换将原始矩阵转化为三对角矩阵;第二步使用Q R迭代算法8-9、D&C算法1 0-1 1或MR R R算法1 2-1 3计算三对角矩阵的特征值分解;最后使用H o u s e h o l d e r正交变换将三对角矩阵的特征向量恢复成原始矩阵的特征向量。与之不同,J a-c o b i方法1 4-1 5无须将原始矩阵转化为三对角