Polynomial filtered Lanczos iterations with applic

Polynomial¯lteredLanczositerationswithapplicationsinDensityFunctionalTheory¤C.Bekas,E.Kokiopoulou,YousefSaadComputerScience&EngineeringDept.UniversityofMinnesota,TwinCities.fbekas,kokiopou,[e-mail]15,2005AbstractThemostexpensivepartofallElectronicStructureCalculationsbasedonDensityFunctionalThe-ory,liesinthecomputationofaninvariantsubspaceassociatedwithsomeofthesmallesteigenvaluesofadiscretizedHamiltonianoperator.Thedimensionofthissubspacetypicallydependsonthetotalnumberofvalenceelectronsinthesystem,andcaneasilyreachhundredsoreventhousandswhenlargesystemswithmanyatomsareconsidered.Atthesametime,thediscretizationofHamiltoniansassociatedwithlargesystemsyieldsverylargematrices,whetherwithplane-waveorrealspacedis-cretizations.Thecombinationofthesetwofactorsresultsinoneofthemostsigni¯cantbottlenecksinComputationalMaterialsScience.Inthispaperweshowhowtoe±cientlycomputealargeinvariantsubspaceassociatedwiththesmallesteigenvaluesofaHermitianmatrixusingpolynomially¯lteredLanczositerations.Theproposedmethoddoesnottrytoextractindividualeigenvaluesandeigen-vectors.Instead,itconstructsanorthogonalbasisoftheinvariantsubspacebycombiningtwomainingredients.The¯rstisa¯lteringtechniquetodampentheundesirablecontributionofthelargesteigenvaluesateachmatrix-vectorproductintheLanczosalgorithm.Thistechniquesemploysawell-selectedlowpass¯lterpolynomial,obtainedviaaconjugateresidual-typealgorithminpolynomialspace.ThesecondingredientistheLanczosalgorithmwithpartialreorthogonalization.Experimentsarereportedtoillustratethee±ciencyoftheproposedschemecomparedtostate-of-the-artimplicitlyrestartedtechniques.Keywords:PolynomialFiltering,ConjugateResidual,LanczosAlgorithm,DensityFunctionalTheory1IntroductionAnumberofscienti¯candengineeringapplicationsrequirethecomputationofalargenumberofeigen-vectorsassociatedwiththesmallesteigenvaluesofalargesymmetric/hermitianmatrix.Animportantexampleisinelectronicstructurecalculationswherethechargedensity%(r)atapointrinspaceiscalculatedfromtheeigenfunctionsªioftheHamiltonianAviatheformula%(r)=noXi=1jªi(r)j2;(1)wherethesummationistakenoveralloccupiedstates(no)ofthesystemunderstudy.ThisisacrucialcalculationinDensityFunctionalTheorysincethepotentialVoftheHamiltonianA=r2+V,depends¤WorksupportedbyNSFgrantsNSF/ITR-0325218andNSF/ITR-0428774,byDOEunderGrantsDE-FG02-03ER25585,DE-FG02-03ER15491,andbytheMinnesotaSupercomputingInstitute1onthechargedensity%,whichinturndependsontheeigenvectorsÃiofA(see(1)),andasaresult,aniterativeloopisrequiredtoachieveself-consistence.Computingthechargedensity%(r)via(1),requireseigenvectors,thoughitismoreaccuratetosaythatwhatisneededisanorthogonalbasisoftheinvariantsubspaceassociatedwiththenoalgebraicallysmallesteigenvalues.Thisisbecause%(r)isinvariantunderanorthogonaltransformationsofthebasisofeigenfunctionsfªig.IfthesymmetricmatrixAisthediscretizationoftheHamiltonianAandthevectorsÃiarethecorrespondingdiscretizationsoftheeigenfunctionsªi(r)withrespecttor,then,thechargedensitiesarethediagonalentriesofthe\functionaldensitymatrixP=QnoQnowithQno=[Ã1;:::;Ãno]:(2)Speci¯cally,thechargedensityatthej-thpointrjisthej-thdiagonalentryofP.Infact,anyorthogonalbasisQwhichspansthesamesubspaceastheeigenvectorsÃi;i=1;:::;nocanbeused.Thisobservationhasleadtoimprovedschemeswhichdonotfocusonextractingindividualeigenvectors.Forexample,[4]showedthatthesemi-orthogonalbasiscomputedbytheLanczosalgorithmwithpartialreorthogonalizationcanbeusedinordertoextractaccurateapproximationstothechargedensity.ThisschemeresultsinsubstantialsavingsrelativetoschemeswhichrelyonthefullreorthogonalizationoftheLanczosvectors.Insimpleterms,theproblemconsideredinthispapercanbestatedasfollows.Givenarealsymmetric(orcomplexHermitian)matrixA2Rn£nwitheigenvalues¸1·¸2·:::·¸n,computetheinvari-antsubspaceSnoassociatedwiththeeigenvalueswhichdonotexceedacertainlimit°.Inelectronicstructures,°istheFermienergylevelandtheinterval[a;°]containsthe(algebraically)smallestoccu-piedeigenstateseigenvalues¸1;:::;¸no.Weassumethatwearegivenaninterval[®;¯]which(tightly)containsthespectrumofA.Thenatureofthealgorithmsusedalsorequiresthat®¸0.Ifthisisnotsatis¯ed,weshiftmatrixAbyascalar¾sothatA+¾Idoesnothaveanynegativeeigenvalues.Methodsforcomputinganinterval[®;¯]whenthisisnotreadilyavailablearediscussedinsection3.1.Thiswellknownproblemreceivedmuchattentioninthepastfewdecades,ascanbeseenfromthesurvey[2]1whichprovidesagoodillustrationofthewealthofalgorithmicresearchinthis¯eld.ForacomprehensivetheoreticaldiscussionoftheproblemonecantorefertoParlett'sclassicbook[24]andthereferencestherein.Whenthenumberofdesiredeigenvaluesisrathersmall,intheorderofafewdozens,avarietyofalgorithmssuccessfullyaddresstheproblem.ThemostgeneralpurposeandextensivelyusedmethodisbasedonimplicitlyrestartedLanczositerations[36]andisimplementedinthesoftwarepackageARPACK[16].However,theproblembecomesparticularlydemandinginthecaseswhenweseektocomputealargenumberofeigenvaluesthatreachdeepintotheinteriorofthespectrumofthematrixathand.Indeed,inelectronicstructurecalculations,thedimension