Program f12jufe
! F12JUF Example Program Text
! Mark 28.4 Release. NAG Copyright 2022.
! .. Use Statements ..
Use, Intrinsic :: iso_c_binding, Only: c_ptr
Use nag_library, Only: f12jaf, f12jbf, f12jff, f12juf, f12jzf, nag_wp, &
x04daf, zsymm, zsytrf, zsytrs
! .. Implicit None Statement ..
Implicit None
! .. Parameters ..
Complex (Kind=nag_wp), Parameter :: cone = (1.0_nag_wp,0.0_nag_wp)
Complex (Kind=nag_wp), Parameter :: czero = (0.0_nag_wp,0.0_nag_wp)
Integer, Parameter :: nin = 5, nout = 6
! .. Local Scalars ..
Type (c_ptr) :: handle
Complex (Kind=nag_wp) :: emid, tmp, ze
Real (Kind=nag_wp) :: eps, r
Integer :: deg, i, ifail, irevcm, iter, j, k, &
lda, ldx, ldy, ldz, lwork, m0, n, &
nconv
! .. Local Arrays ..
Complex (Kind=nag_wp), Allocatable :: a(:,:,:), az(:,:), d(:), work(:), &
x(:,:), y(:,:), z(:,:)
Real (Kind=nag_wp), Allocatable :: resid(:)
Integer, Allocatable :: ipiv(:)
! .. Executable Statements ..
Write (nout,*) 'F12JUF Example Program Results'
Write (nout,*)
! Skip heading in data file
Read (nin,*)
Read (nin,*) n, deg
lda = n
m0 = n
ldx = n
ldy = n
ldz = n
lwork = 64*n
Allocate (x(ldx,m0),y(ldy,m0),z(ldz,m0),resid(m0),d(m0),a(lda,n,deg+1), &
az(n,n),ipiv(n),work(lwork))
! Read the upper triangular part of the matrices from data file
Do j = 1, deg + 1
Read (nin,*)(a(i,i:n,j),i=1,n)
End Do
! Initialize the data handle
ifail = 0
Call f12jaf(handle=handle,ifail=ifail)
! Set the ellipse minor/major axis ratio
Call f12jbf(handle=handle,str='Ellipse Contour Ratio = 80',ifail=ifail)
emid = (-2.0_nag_wp,-2.0_nag_wp)
r = 1.5_nag_wp
! Generate the contour
Call f12jff(handle=handle,emid=emid,r=r,ifail=ifail)
irevcm = 0
revcomm: Do
Call f12juf(handle=handle,irevcm=irevcm,deg=deg,ze=ze,n=n,x=x,ldx=ldx, &
y=y,ldy=ldy,k=k,m0=m0,nconv=nconv,d=d,z=z,ldz=ldz,eps=eps,iter=iter, &
resid=resid,ifail=ifail)
Select Case (irevcm)
Case (1)
! Form the matrix \sum ze^i A_i
Do j = 1, n
az(1:j,j) = a(1:j,j,1)
End Do
Do i = 1, deg
Do j = 1, n
az(1:j,j) = az(1:j,j) + (ze**i)*a(1:j,j,i+1)
End Do
End Do
! Compute a Bunch-Kaufman factorization of \sum ze^i A_i
! The NAG name equivalent of zsytrf is f07nrf
Call zsytrf('Upper',n,az,n,ipiv,work,lwork,ifail)
If (ifail/=0) Then
Exit revcomm
End If
Cycle revcomm
Case (2)
! Solve the linear system (\sum ze^i A_i)w = y, overwriting y with w
! The NAG name equivalent of zsytrs is f07nsf
Call zsytrs('Upper',n,m0,az,n,ipiv,y,ldy,ifail)
If (ifail/=0) Then
Exit revcomm
End If
Cycle revcomm
Case (3)
! Compute x <- A_k z
Call zsymm('Left','Upper',n,m0,cone,a(1,1,k+1),lda,z,ldz,czero,x, &
ldx)
Cycle revcomm
Case Default
Exit revcomm
End Select
End Do revcomm
If (ifail==0) Then
! Print solution
Write (nout,*) 'Eigenvalues'
Write (nout,99999) d(1:nconv)
Flush (nout)
! Normalize the eigenvectors: unit first element
Do i = 1, nconv
tmp = z(1,i)
z(1:n,i) = z(1:n,i)/tmp
End Do
Call x04daf('General',' ',n,nconv,z,ldz,'Eigenvectors',ifail)
Else
Write (nout,99998) 'Failure in f12juf IFAIL =', ifail
End If
! Destroy the data handle
Call f12jzf(handle=handle,ifail=ifail)
99999 Format ((3X,4(' (',F7.4,',',F7.4,')',:)))
99998 Format (1X,A,I4)
End Program f12jufe