Program f12apfe ! F12APF Example Program Text ! Mark 24 Release. NAG Copyright 2012. ! .. Use Statements .. Use nag_library, Only: dznrm2, f12anf, f12apf, f12aqf, f12arf, f12asf, & nag_wp, zgttrf, zgttrs ! .. Implicit None Statement .. Implicit None ! .. Parameters .. Complex (Kind=nag_wp), Parameter :: one = (1.0_nag_wp,0.0_nag_wp) Complex (Kind=nag_wp), Parameter :: two = (2.0_nag_wp,0.0_nag_wp) Integer, Parameter :: imon = 0, nerr = 6, nin = 5, nout = 6 ! .. Local Scalars .. Complex (Kind=nag_wp) :: h, h2, rho, s, s1, s2, s3, sigma Integer :: ifail, info, irevcm, j, lcomm, ldv, & licomm, n, nconv, ncv, nev, niter, & nshift, nx ! .. Local Arrays .. Complex (Kind=nag_wp), Allocatable :: comm(:), d(:,:), dd(:), dl(:), & du(:), du2(:), mx(:), resid(:), & v(:,:), x(:) Integer, Allocatable :: icomm(:), ipiv(:) ! .. Intrinsic Procedures .. Intrinsic :: cmplx, int, max ! .. Executable Statements .. Write (nout,*) 'F12APF Example Program Results' Write (nout,*) ! Skip heading in data file Read (nin,*) Read (nin,*) nx, nev, ncv n = nx*nx ! Initialize communication arrays. ! Query the required sizes of the communication arrays. licomm = -1 lcomm = -1 Allocate (icomm(max(1,licomm)),comm(max(1,lcomm))) ifail = 0 Call f12anf(n,nev,ncv,icomm,licomm,comm,lcomm,ifail) licomm = icomm(1) lcomm = int(comm(1)) Deallocate (icomm,comm) Allocate (icomm(max(1,licomm)),comm(max(1,lcomm))) ifail = 0 Call f12anf(n,nev,ncv,icomm,licomm,comm,lcomm,ifail) ! Set the mode. ifail = 0 Call f12arf('SHIFTED INVERSE',icomm,comm,ifail) sigma = cmplx(0,kind=nag_wp) rho = (10.0_nag_wp,0.0_nag_wp) h = one/cmplx(n+1,kind=nag_wp) h2 = h*h s = rho/two s1 = -one/h2 - s/h s2 = two/h2 - sigma s3 = -one/h2 + s/h Allocate (dl(n-1),dd(n),du(n-1),du2(n-2),ipiv(n)) dl(1:n-1) = s1 dd(1:n-1) = s2 du(1:n-1) = s3 dd(n) = s2 ! The NAG name equivalent of zgttrf is f07crf Call zgttrf(n,dl,dd,du,du2,ipiv,info) If (info/=0) Then Write (nerr,99999) info Go To 110 End If ldv = n Allocate (resid(n),v(ldv,ncv),x(n),mx(n),d(ncv,2)) irevcm = 0 ifail = -1 100 Continue Call f12apf(irevcm,resid,v,ldv,x,mx,nshift,comm,icomm,ifail) If (irevcm/=5) Then If (irevcm==-1 .Or. irevcm==1) Then ! Perform x <--- OP*x = inv[A-SIGMA*I]*x ! The NAG name equivalent of zgttrs is f07csf Call zgttrs('N',n,1,dl,dd,du,du2,ipiv,x,n,info) If (info/=0) Then Write (nerr,99998) info Go To 110 End If Else If (irevcm==4 .And. imon/=0) Then ! Output monitoring information Call f12asf(niter,nconv,d,d(1,2),icomm,comm) ! The NAG name equivalent of dznrm2 is f06jjf Write (6,99997) niter, nconv, dznrm2(nev,d(1,2),1) End If Go To 100 End If If (ifail==0) Then ! Post-Process using F12AQF to compute eigenvalues/vectors. ifail = 0 Call f12aqf(nconv,d,v,ldv,sigma,resid,v,ldv,comm,icomm,ifail) Write (nout,99996) nconv Do j = 1, nconv Write (nout,99995) j, d(j,1) End Do End If 110 Continue 99999 Format (1X,'** Error status returned by ZGTTRF, INFO =',I12) 99998 Format (1X,'** Error status returned by ZGTTRS, INFO =',I12) 99997 Format (1X,'Iteration',1X,I3,', No. converged =',1X,I3,', norm o', & 'f estimates =',E16.8) 99996 Format (1X/' The ',I4,' Ritz values of smallest magnitude are:'/) 99995 Format (1X,I8,5X,'( ',F12.4,' , ',F12.4,' )') End Program f12apfe