Program f12jsfe
! F12JSF Example Program Text
! Mark 29.2 Release. NAG Copyright 2023.
! .. Use Statements ..
Use, Intrinsic :: iso_c_binding, Only: c_ptr
Use nag_library, Only: f12jaf, f12jbf, f12jff, f12jsf, 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, ze
Real (Kind=nag_wp) :: eps, r
Integer :: i, ifail, irevcm, iter, lda, ldb, &
ldx, ldy, ldz, lwork, m0, n, nconv
! .. Local Arrays ..
Complex (Kind=nag_wp), Allocatable :: a(:,:), az(:,:), b(:,:), d(:), &
work(:), x(:,:), y(:,:), z(:,:)
Real (Kind=nag_wp), Allocatable :: resid(:)
Integer, Allocatable :: ipiv(:)
! .. Executable Statements ..
Write (nout,*) 'F12JSF Example Program Results'
Write (nout,*)
! Skip heading in data file
Read (nin,*)
Read (nin,*) n
lda = n
ldb = 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), &
b(ldb,n),az(n,n),ipiv(n),work(lwork))
! Read the upper triangular part of the matrix A from data file
Read (nin,*)(a(i,i:n),i=1,n)
! Read the upper triangular part of the matrix B from data file
Read (nin,*)(b(i,i:n),i=1,n)
! Initialize the data handle
ifail = 0
Call f12jaf(handle,ifail)
! Set the ellipse rotation angle and minor/major axis ratio
Call f12jbf(handle,'Ellipse Contour Ratio = 10',ifail)
Call f12jbf(handle,'Ellipse Rotation Angle = 45',ifail)
emid = (0.0_nag_wp,0.0_nag_wp)
r = 1.0_nag_wp
! Generate the contour
Call f12jff(handle,emid,r,ifail)
irevcm = 0
revcomm: Do
Call f12jsf(handle,irevcm,ze,n,x,ldx,y,ldy,m0,nconv,d,z,ldz,eps,iter, &
resid,ifail)
Select Case (irevcm)
Case (1)
! Form the matrix ze B - A
Do i = 1, n
az(1:i,i) = ze*b(1:i,i) - cone*a(1:i,i)
End Do
! Compute a Bunch-Kaufman factorization of ze B -A
! 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 (ze I -A)w = y
! 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 <- Az
Call zsymm('Left','Upper',n,m0,cone,a,lda,z,ldz,czero,x,ldx)
Cycle revcomm
Case (4)
! Compute x <- Bz
Call zsymm('Left','Upper',n,m0,cone,b,ldb,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)
Call x04daf('General',' ',n,nconv,z,ldz,'Eigenvectors',ifail)
Else
Write (nout,99998) 'Failure in f12jsf IFAIL =', ifail
End If
! Destroy the data handle
Call f12jzf(handle,ifail)
99999 Format ((3X,4(' (',F7.4,',',F7.4,')',:)))
99998 Format (1X,A,I4)
End Program f12jsfe