NAGLibraryRoutineDocumentF08KBFDGESVDbeforeusingthisroutinepleasere - PDF document

1 NAGLibraryRoutineDocumentF08KBF(DGESVD)b
NAGLibraryRoutineDocumentF08KBF(DGESVD)beforeusingthisroutine,pleasereadtheUsersNoteforyourimplementationtochecktheinterpretationofbolditalicisedtermsandotherimplementation-dependentdetails.1Purpose F08…Least-squaresandEigenvalueProblems(LAPACK)Mark25F08KBF.1 Nocolumnsof(noleftsingularvectors)arecomputed.:JOBUA,S,OorN.JOBVT…CHARACTER(1)Onentry:speci“esoptionsforcomputingallorpartofthematrixrowsofarereturnedinthearrayVT.The“rstminm;nrowsof(therightsingularvectors)arereturnedinthearrayVT.The“rstminm;nrowsof(therightsingularvectors)areoverwrittenonthearrayA.Norowsof(norightsingularvectors)arecomputed.A,S,OorN;JOBVTandJOBUcannotbothbeO.M…INTEGEROnentry,thenumberofrowsofthematrixN…INTEGEROnentry,thenumberofcolumnsofthematrix…REAL(KIND=nag_wp)array:theseconddimensionofthearrayAmustbeatleastmax1Onentry:theOnexit:ifJOBUO,Aisoverwrittenwiththe“rstminm;ncolumnsof(theleftsingularvectors,storedcolumn-wise).IfJOBVTO,Aisoverwrittenwiththe“rstminm;nrowsof(therightsingularvectors,storedrow-wise).IfJOBUOandJOBVTO,thecontentsofAaredestroyed.LDA…INTEGEROnentry:the“rstdimensionofthearrayAasdeclaredinthe(sub)programfromwhichF08KBF(DGESVD)isca

2 lled.:LDAmax1…REAL(KIND=nag_wp)array:the
lled.:LDAmax1…REAL(KIND=nag_wp)array:thedimensionofthearraySmustbeatleastmax1minMOnexit:thesingularvaluesof,sortedsothatS NAGLibraryManualF08KBF.2Mark25 …REAL(KIND=nag_wp)array:theseconddimensionofthearrayUmustbeatleastmax1ifJOBUmax1minMifJOBUS,andatleast1otherwise.Onexit:ifJOBUA,UcontainstheorthogonalmatrixIfJOBUS,Ucontainsthe“rstminm;ncolumnsof(theleftsingularvectors,storedIfJOBUNorO,Uisnotreferenced.LDU…INTEGEROnentry:the“rstdimensionofthearrayUasdeclaredinthe(sub)programfromwhichF08KBF(DGESVD)iscalled.ifJOBUAorS,LDUmax1otherwiseLDU10:LDVT…REAL(KIND=nag_wp)array:theseconddimensionofthearrayVTmustbeatleastmax1ifJOBVTAorS,andatleast1otherwise.Onexit:ifJOBVTA,VTcontainstheorthogonalmatrixIfJOBVTS,VTcontainsthe“rstminm;nrowsof(therightsingularvectors,storedIfJOBVTNorO,VTisnotreferenced.11:LDVT…INTEGEROnentry:the“rstdimensionofthearrayVTasdeclaredinthe(sub)programfromwhichF08KBF(DGESVD)iscalled.ifJOBVTA,LDVTmax1ifJOBVTS,LDVTmax1minMotherwiseLDVT12:max1LWORKðÞޅREAL(KIND=nag_wp)arrayWorkspaceOnexit:ifINFO,WORKreturnstheoptimalLWORK.IfINFO,WORKminMðÞÞcontainstheunconvergedsuperdiagonalelementsofanupperbidiagonalmatrixwhosediag

3 onalisinS(notnecessarilysorted).satis“es
onalisinS(notnecessarilysorted).satis“esUBVsoithasthesamesingularvaluesas,andsingularvectorsrelatedby13:LWORK…INTEGEROnentry:thedimensionofthearrayWORKasdeclaredinthe(sub)programfromwhichF08KBF(DGESVD)iscalled.IfLWORK1,aworkspacequeryisassumed;theroutineonlycalculatestheoptimalsizeoftheWORKarray,returnsthisvalueasthe“rstentryoftheWORKarray,andnoerrormessagerelatedtoLWORKisissued.Suggestedvalue:foroptimalperformance,LWORKshouldgenerallybelarger.ConsiderincreasingLWORKbyatleastminM,whereistheoptimalblocksize:LWORKmax1minMðÞþmaxMminM F08…Least-squaresandEigenvalueProblems(LAPACK)Mark25F08KBF.3 14:INFO…INTEGEROnexit:INFO0unlesstheroutinedetectsanerror(seeSection6).6ErrorIndicatorsandWarningsIfINFO,argumenthadanillegalvalue.Anexplanatorymessageisoutput,andexecutionoftheprogramisterminated.IfF08KBF(DGESVD)didnotconverge,INFOspeci“eshowmanysuperdiagonalsofanintermediatebidiagonalformdidnotconvergetozero.7AccuracyThecomputedsingularvaluedecompositionisnearlytheexactsingularvaluedecompositionforanearbymatrix,whereisthemachineprecision.Inaddition,thecomputedsingularvectorsarenearlyorthogonaltoworkingprecision.SeeSection4.9ofAnderson

4 etal.(1999)forfurtherdetails.8Parallelis
etal.(1999)forfurtherdetails.8ParallelismandPerformanceF08KBF(DGESVD)isthreadedbyNAGforparallelexecutioninmultithreadedimplementationsoftheNAGLibrary.F08KBF(DGESVD)makescallstoBLASand/orLAPACKroutines,whichmaybethreadedwithinthevendorlibraryusedbythisimplementation.Consultthedocumentationforthevendorlibraryforfurtherinformation.PleaseconsulttheX06ChapterIntroductionforinformationonhowtocontrolandinterrogatetheOpenMPenvironmentusedwithinthisroutine.PleasealsoconsulttheUsersNoteforyourimplementationforanyadditionalimplementation-speci“cinformation.9FurtherCommentsThetotalnumberof”oating-pointoperationsisapproximatelyproportionaltootherwise.Thesingularvaluesarereturnedindescendingorder.ThecomplexanalogueofthisroutineisF08KPF(ZGESVD).10ExampleThisexample“ndsthesingularvaluesandleftandrightsingularvectorsofthe6by4matrix541670090071220790352452391togetherwithapproximateerrorboundsforthecomputedsingularvaluesandvectors.NAGLibraryManualF08KBF.4Mark25 TheexampleprogramforF08KDF(DGESDD)illustrates“ndingasingularvaluedecompositionforthecase10.1ProgramTextProgramf08kbfe!F08KBFExampleProgramText!Mark25Release.NAGCopyright2014.!..UseSt

5 atements..Usenag_library,Only:dgesvd,nag
atements..Usenag_library,Only:dgesvd,nag_wp!..ImplicitNoneStatement..ImplicitNone!..Parameters..Integer,Parameter::nb=64,nin=5,nout=6,&prerr=0!..LocalScalars..Integer::i,info,lda,ldu,ldvt,lwork,&m,n!..LocalArrays..Real(Kind=nag_wp),Allocatable::a(:,:),a_copy(:,:),b(:),s(:),&u(:,:),vt(:,:),work(:)Real(Kind=nag_wp)::dummy(1,1)!..IntrinsicProcedures..Intrinsic::max,min,nint!..ExecutableStatements..Write(nout,*)F08KBFExampleProgramResultsWrite(nout,*)!SkipheadingindatafileRead(nin,*)Read(nin,*)m,nlda=mldu=mldvt=nAllocate(a(lda,n),a_copy(m,n),s(n),vt(ldvt,n),u(ldu,m),b(m))!ReadthembynmatrixAfromdatafileRead(nin,*)(a(i,1:n),i=1,m)!ReadtherighthandsideofthelinearsystemRead(nin,*)b(1:m)a_copy(1:m,1:n)=a(1:m,1:n)!Useroutineworkspacequerytogetoptimalworkspace.lwork=-1!TheNAGnameequivalentofdgesvdisf08kbfCalldgesvd(A,S,m,n,a,lda,s,u,ldu,vt,ldvt,dummy,lwork,info)!Makesurethatthereisenoughworkspaceforblocksizenb.lwork=max(m+4*n+nb*(m+n),nint(dummy(1,1)))Allocate(work(lwork))!Computethesingularvaluesandleftandrightsingularvectors!ofA.!TheNAGnameequivalentofdgesvdisf08kbfCalldgesvd(A,S,m,n,a,lda,s,u,ldu,vt,ldvt,work,lwork,info)

6 If(info/=0)ThenWrite(nout,99999)Failure
If(info/=0)ThenWrite(nout,99999)FailureinDGESVD.INFO=,info99999Format(1X,A,I4)GoTo100EndIf!PrintthesignificantsingularvaluesofAWrite(nout,*)SingularvaluesofA:F08…Least-squaresandEigenvalueProblems(LAPACK)Mark25F08KBF.5 Write(nout,99998)s(1:min(m,n))99998Format(1X,4(3X,F11.4))If�(prerr0)ThenCallcompute_error_bounds(m,n,s)EndIfIf�(mn)Then!ComputeV*Inv(S)*U^T*btogetleast-squaressolution.Callcompute_least_squares(m,n,a_copy,m,u,ldu,vt,ldvt,s,b)EndIf100ContinueContainsSubroutinecompute_least_squares(m,n,a,lda,u,ldu,vt,ldvt,s,b)!..UseStatements..Usenag_library,Only:dgemv,dnrm2!..ImplicitNoneStatement..ImplicitNone!..ScalarArguments..Integer,Intent(In)::lda,ldu,ldvt,m,n!..ArrayArguments..Real(Kind=nag_wp),Intent(In)::a(lda,n),s(n),u(ldu,m),&Real(Kind=nag_wp),Intent(Inout)::b(m)!..LocalScalars..Real(Kind=nag_wp)::alpha,beta,norm!..LocalArrays..Real(Kind=nag_wp),Allocatable::x(:),y(:)!..IntrinsicProcedures..Intrinsic::allocated!..ExecutableStatements..Allocate(x(n),y(n))!ComputeV*Inv(S)*U^T*btogetleast-squaressolution.!y=U^Tb!TheNAGnameequivelentofdgemvisf06pafalpha=1._nag_wpbeta=0._nag_wpCalldgemv(T,m,n,alpha,u

7 ,ldu,b,1,beta,y,1)y(1:n)=y(1:n)/s(1:n)!x
,ldu,b,1,beta,y,1)y(1:n)=y(1:n)/s(1:n)!x=VyCalldgemv(T,n,n,alpha,vt,ldvt,y,1,beta,x,1)Write(nout,*)Write(nout,*)Leastsquaressolution:Write(nout,99999)x(1:n)!Findnormofresidual||b-Ax||.alpha=-1._nag_wpbeta=1._nag_wpCalldgemv(N,m,n,alpha,a,lda,x,1,beta,b,1)norm=dnrm2(m,b,1)Write(nout,*)Write(nout,*)NormofResidual:Write(nout,99999)normIf(allocated(x))Deallocate(x)If(allocated(y))Deallocate(y)99999Format(1X,4(3X,F11.4))EndSubroutinecompute_least_squaresNAGLibraryManualF08KBF.6Mark25 Subroutinecompute_error_bounds(m,n,s)!Errorestimatesforsingularvaluesandvectorsiscomputed!andprintedhere.!..UseStatements..Usenag_library,Only:ddisna,nag_wp,x02ajf!..ImplicitNoneStatement..ImplicitNone!..ScalarArguments..Integer,Intent(In)::m,n!..ArrayArguments..Real(Kind=nag_wp),Intent(In)::s(n)!..LocalScalars..Real(Kind=nag_wp)::eps,serrbdInteger::i,info!..LocalArrays..Real(Kind=nag_wp),Allocatable::rcondu(:),rcondv(:),&uerrbd(:),verrbd(:)!..ExecutableStatements..Allocate(rcondu(n),rcondv(n),uerrbd(n),verrbd(n))!Getthemachineprecision,EPSandcomputetheapproximate!errorboundforthecomputedsingularvalues.Notethatfor!the2-norm,S(1)=norm(A)eps=

8 x02ajf()serrbd=eps*s(1)!CallDDISNA(F08FL
x02ajf()serrbd=eps*s(1)!CallDDISNA(F08FLF)toestimatereciprocalcondition!numbersforthesingularvectorsCallddisna(Left,m,n,s,rcondu,info)Callddisna(Right,m,n,s,rcondv,info)!ComputetheerrorestimatesforthesingularvectorsDoi=1,nuerrbd(i)=serrbd/rcondu(i)verrbd(i)=serrbd/rcondv(i)EndDo!Printtheapproximateerrorboundsforthesingularvalues!andvectorsWrite(nout,*)Write(nout,*)ErrorestimateforthesingularvaluesWrite(nout,99999)serrbdWrite(nout,*)Write(nout,*)ErrorestimatesfortheleftsingularvectorsWrite(nout,99999)uerrbd(1:n)Write(nout,*)Write(nout,*)ErrorestimatesfortherightsingularvectorsWrite(nout,99999)verrbd(1:n)99999Format(4X,1P,6E11.1)EndSubroutinecompute_error_boundsEndProgramf08kbfe10.2ProgramDataF08KBFExampleProgramData64:ValuesofMandN2.27-1.541.15-1.940.28-1.670.94-0.78-0.48-3.090.99-0.21F08…Least-squaresandEigenvalueProblems(LAPACK)Mark25F08KBF.7 1.071.220.790.63-2.352.93-1.452.300.62-7.391.03-2.57:EndofmatrixA1.01.01.01.01.01.0:RHSb(1:m)10.3ProgramResultsF08KBFExampleProgramResultsSingularvaluesofA:9.99663.68311.35690.5000Leastsquaressolution:-0.0563-0.17000.82020.5545NormofResidual:NAGLibraryManualF08KBF.8(last)Ma