大型稀疏矩阵的LU分解及特征值求解

2013.7.20大型稀疏矩阵的LU分解及特征值求解陈英时2016.1.9稀疏矩阵求解的广泛应用•矩阵求解是数值计算的核心[1]•稀疏矩阵求解是数值计算的关键之一偏微分方程，积分方程，特征值，优化…万阶以上densematrix不可行•稀疏矩阵求解往往是资源瓶颈时间瓶颈，内存，外存等瓶颈•直接法基于高斯消元法，即计算A的LU分解。A通常要经过一系列置换排序，可归并为左置换矩阵P，右置换矩阵Q。基本步骤如下：1）符号分析：得到置换排序矩阵PQ2）数值分解：3）前代回代：I.S.Duff,A.M.Erisman,andJ.K.Reid.DirectMethodsforSparseMatrices.London:OxfordUniv.Press,1986.J.J.Dongarra,I.S.Duff,...NumericalLinearAlgebraforHigh-PerformanceComputers.G.H.Golob,C.F.Vanloan.MatrixComputations.TheJohnsHopkinsUniversityPress.1996稀疏矩阵复杂、多变•基本参数对称性，稀疏性，非零元分布•敏感性，病态矩阵条件数•格式多变Harwell-BoeingExchangeFormat。。。•测试集Harwell-BoeingSparseMatrixCollectionUFsparsematrixcollection求解器的飞速发展•BBMAT阶，分解后元素超过四千万.1988巨型机cray-2上1000秒20034Gumfpack432.6秒[4]20063.0GGSS1.215秒20123.0G4核GSS2.34秒2014i78核GSS2.41秒•硬件的发展CPU，内存等•稀疏算法逐渐成熟稀疏排序，密集子阵multifrontal,supernodal…•数学库BLAS，LAPACK分解算法的关键根据符号分析，数值分解算法的不同，直接法有以下几类[15]：1）generaltechnique(通用方法)：主要采用消去树等结构进行LU分解。适用于非常稀疏的非结构化矩阵。2）frontalscheme(波前法)LU分解过程中，将连续多行合并为一个密集子块(波前)，对这个子块采用BLAS等高效数学库进行分解。3）multifrontalmethod(多波前法)将符号结构相同的多行(列)合并为多个密集子块，以这些密集子块为单位进行LU分解。这些子块的生成，消去，装配，释放等都需要特定的数据结构及算法。1零元是动态的概念，需要稀疏排序,减少注入元(fill-in)2密集子阵稀疏LU分解（不考虑零元的）LU分解稀疏排序排序算法的作用是减少矩阵LU分解过程中产生的注入元，求解矩阵的最优排序方法是个NP完全问题（不能够在合理的时间内进行求解），对具体矩阵而言，目前也没有方法或指标来判定哪种算法好。因此实测不同的算法，对比产生的注入元，是确定排序算法的可靠依据。主要的排序算法有最小度排序（MMD，minimumdegreeorderingalgorithm）和嵌套排序（nesteddissection）两种。矩阵排序方面相应的成熟软件库有AMD[12]、COLAMD、METIS[13]等。YannokakisM.Computingtheminimumfill-ininNP-Complete.SIAMJ.AlgebraicDiscreteMethods,1981,2:77~79近似最小度排序算法–商图近似最小度排序（AMD，ApproximateMinimumDegreeOrderingAlgorithm）算法于1996年左右由PatrickR.Amestoy等学者提出AMESTOY,P.R.,DAVIS,...1996a.Anapproximateminimumdegreeorderingalgorithm.SIAMJ.MatrixAnal.Applic.17,4,886–905.为何需要密集子块(DenseMatrix)=Benchmark多波前法(multifrontal)简介•发展DuffandRaid[2]J.W.H.Liu等分析，改进[3]T.A.Davis开发UMFPACK[4]•基本算法利用稀疏矩阵的特性，得到一系列密集子阵（波前）。将LU分解转化为对这些波前的装配，消去，更新等操作。•多波前法的优点波前是densematrix,可直接调用高性能库（BLAS等）密集子阵可以节省下标存储提高并行性•目前主要的求解器UMFPACK,,GSS,MUMPS,PARDISO,WSMP,HSLMA41等LU分解形成frontal10阶矩阵。蓝点代表非零元。红点表示分解产生的注入元(fill-in)Frontal划分{a},{b}{c}{d}{e}{f,g}{h,i,j}Frontal的装配，消去，更新过程achc··h··{c}{f,g}{b}{e}{h,i,j}{a}c,g,hg··h··beje··j··f,g,hgg·h··e,i,ji··j··h,i,jii·j·j{d}d,i,ji··j··消去树消去树消去树是符号分析的关键结构，其每个节点对应于矩阵的一列（行），该节点只与其父节点相连，父节点定义如下：J.W.H.Liu.TheRoleofEliminationTreesinSparseFactorization.SIAMJ.MatrixAnal.Applic.,11(1):134-172,1990.J.W.H.Liu.EliminationstructuresforunsymmetricsparseLUfactors.SIAMJ.MatrixAnal.Appl.14,no.2,334--352,1993.GSS简介•标准C开发，适用于各种平台•比INTELPARDISO更快，更稳定•CPU/GPU混合计算•突破32位Windows内存限制•32个用户参数•支持用户定制模块高校，研究所中国电力科学研究院香港大学南京大学河海大学中国石油大学电子科技大学三峡大学山东大学userdetailWhytheychooseGSScrosslightIndustryleaderinTCADsimulationHybridGPU/CPUversion,morethan2timesfasterthanPARDISO,MUMPSandothersparsesolvers.soilvisionThemosttechnicallyadvancedsuiteof1D/2D/3DgeotechnicalsoftwareMuchfasterthantheirownsparsesolver.FEMconsultingTheleadingresearchteamsintheareaoftheFiniteElementMethodsince1967GSSisfasterthanPARDISOandprovidemanycustommodule.GSCADLeadingbuildingsoftwareinChinaGSSprovideauser-specificmoduletodealwithill-conditionedmatrix.ICAROSAglobalturnkeygeospatialmappingserviceproviderandstateoftheartphotogrammetrictechnologiesdeveloper.GSSisfasterthanPARDISO.Alsoprovidesometechnicalhelp.EPRIChinaElectricPowerResearchInstitute3-4timesfasterthanKLU他们为什么选择GSS?GSS–加权消去树工作量消去树基于消去树结构来计算数值分解的工作量，将LU分解划分为多个独立的任务，为高效并行计算奠定基础。GSS-双阈值列选主元算法GSS-CPU/GPU混合计算1)AfterdividesLUfactorizationintomanyparalleltasks,GSSwilluseadaptivestrategytorunthesetasksindifferenthardware(CPU,GPU…).Thatis,ifGPUhavehighcomputingpower,thenitwillrunmoretasksautomatically.IfCPUismorepowerful,thenGSSwillgiveitmoretasks.2)Andfurthermore,ifCPUisfree(havefinishedalltasks)andGPUstillrunatask,thenGSSwilldividethistasktosomesmalltasksthenassignsomechild-tasktoCPU,thenCPUdocomputingagain.Sogethigherparallelefficiency.3)GSSwillalsodosometestingtogetbestparametersfordifferenthardware.GSS–求解频域谱元方法生成的矩阵•矩阵较稠密约40万阶，15亿个非零元GSS约需15G内存•需要求解多个右端项32个右端项需80秒？•进一步优化CPU/GPU混合计算数值分解约35秒重复利用符号分析结果根据矩阵的特殊结构来进一步减少非零元估计80/4=20秒•一次LU分解符号分解时间50秒数值分解时间46秒回代2.5秒对比测试ThetestmatricesareallfromtheUFsparsematrixcollectionPARDISOisfromIntelComposerXE2013SP1.GSS2.4useCPU-GPUhybridcomputing.ThetestingCPUisINTELCorei7-4770(3.4GHz)with24Gmemory.ThegraphicscardisASUSGTX780(withcomputecapability3.5).NVIDIACUDAToolkitis5.5.TheoperatingsystemisWindows764.Bothsolversusedefaultparameters.Forlargematricesneedlongtimecomputing,GSS2.4isNearly3timesfasterthanPARDISO.Formatricesneedshorttimecomputing,PARDISOisfasterthanGSS.OnereasonisthatcomplexsynchronizationbetweenCPU/GPUdoneedsomeextratime.大型稀疏矩阵的特征值求解重视xAxA为全过程动态仿真程序中大规模稀疏矩阵为特征值，x为对应的特征向量150Yearsoldandstillalive:eigenproblems1997---byHenkA.vanderVorst,GeneH.Golub稀疏LU分解-理论上即高斯消元法稀疏特征值–趋近于纯数学代数特征值问题J.H.WilkinsonMatrixAlgorithms||EigenSystems.G.W.StewartG.H.Golob,C.F.Vanloan.MatrixComputations.TheJohnsHopkinsUniversityPress.1996TemplatesfortheSolutionofAlgebraicEigenvalueProblems:APracticalGui