Program f12jkfe
! F12JKF Example Program Text
! Mark 30.1 Release. NAG Copyright 2024.
! Intrinsic Procedures ..
! .. Use Statements ..
Use, Intrinsic :: iso_c_binding, Only: c_ptr
Use nag_library, Only: f11dnf, f11dqf, f11xaf, f11znf, f12jaf, f12jbf, &
f12jff, f12jkf, f12jzf, nag_wp, x02ajf, x04daf
! .. Implicit None Statement ..
Implicit None
! .. Parameters ..
Complex (Kind=nag_wp), Parameter :: cone = (1.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, rnorm, tol
Integer :: i, ifail, irevcm, iter, itn, la, &
ldx, ldy, ldz, liwork, lwork, m, m0, &
n, nconv, nnz, nnzc, nnzch, nnzz, &
nnzzh, npivm
! .. Local Arrays ..
Complex (Kind=nag_wp), Allocatable :: az(:), azh(:), d(:), w(:), &
work(:), y(:,:)
Real (Kind=nag_wp), Allocatable :: a(:), resid(:), x(:,:), z(:,:)
Integer, Allocatable :: icol(:), icolz(:), icolzh(:), &
idiag(:), idiagh(:), ipiv(:), &
ipivp(:), ipivph(:), ipivq(:), &
ipivqh(:), irow(:), irowz(:), &
irowzh(:), istr(:), istrh(:), &
iwork(:)
! .. Intrinsic Procedures ..
Intrinsic :: cmplx, conjg, max, min, sqrt
! .. Executable Statements ..
Write (nout,*) 'F12JKF Example Program Results'
Write (nout,*)
! Skip heading in data file
Read (nin,*)
Read (nin,*) n
Read (nin,*) nnz
m0 = n
ldx = n
ldy = n
ldz = n
m = min(n,50)
lwork = max(64*n,4*n+m*(m+n+5)+121)
la = 2*(nnz+n)
liwork = 7*n + 2
tol = sqrt(x02ajf())
Allocate (x(ldx,m0),y(ldy,2*m0),z(ldz,m0),resid(2*m0),d(m0),a(nnz),w(n), &
icol(nnz),irow(nnz),az(la),azh(la),icolz(la),irowz(la),icolzh(la), &
irowzh(la),ipiv(n),iwork(liwork),idiag(n),ipivp(n),ipivq(n),ipivph(n), &
ipivqh(n),istrh(n+1),idiagh(n),istr(n+1),work(lwork))
! Read the matrix A from data file
Do i = 1, nnz
Read (nin,*) a(i), irow(i), icol(i)
End Do
! Initialize the data handle
ifail = 0
Call f12jaf(handle,ifail)
! Set the ellipse rotation angle, and the minor/major axis ratio
Call f12jbf(handle,'Ellipse Rotation Angle = 42',ifail)
Call f12jbf(handle,'Ellipse Contour Ratio = 20',ifail)
emid = (2.5_nag_wp,1.2_nag_wp)
r = 1.4_nag_wp
! Generate the contour
Call f12jff(handle,emid,r,ifail)
irevcm = 0
revcomm: Do
Call f12jkf(handle,irevcm,ze,n,x,ldx,y,ldy,m0,nconv,d,z,ldz,eps,iter, &
resid,ifail)
Select Case (irevcm)
Case (1)
! Form the sparse matrix ze I - A
nnzz = nnz + n
az(1:nnz) = cmplx(-a(1:nnz),kind=nag_wp)
irowz(1:nnz) = irow(1:nnz)
icolz(1:nnz) = icol(1:nnz)
Do i = 1, n
irowz(nnz+i) = i
icolz(nnz+i) = i
az(nnz+i) = ze
End Do
! Sort the elements into correct coordinate storage format
Call f11znf(n,nnzz,az,irowz,icolz,'S','R',istr,iwork,ifail)
! Form incomplete LU factorization of ze I - A
Call f11dnf(n,nnzz,az,la,irowz,icolz,0,0.0_nag_wp,'P','N',ipivp, &
ipivq,istr,idiag,nnzc,npivm,iwork,liwork,ifail)
If (ifail/=0) Then
Exit revcomm
End If
Cycle revcomm
Case (2)
! Solve the linear system (ze I - A)w = y, with m0 righthand sides
Do i = 1, m0
! We will need to store the final result in y
w(1:n) = y(1:n,i)
! Initial guess
y(1:n,i) = cone
Call f11dqf('RGMRES',n,nnzz,az,la,irowz,icolz,ipivp,ipivq,istr, &
idiag,w,m,tol,500,y(1,i),rnorm,itn,work,lwork,ifail)
End Do
If (ifail/=0) Then
Exit revcomm
End If
Cycle revcomm
Case (3)
! Form the sparse matrix (ze I - A)^H
nnzzh = nnz + n
azh(1:nnz) = cmplx(-a(1:nnz),kind=nag_wp)
irowzh(1:nnz) = icol(1:nnz)
icolzh(1:nnz) = irow(1:nnz)
Do i = 1, n
irowzh(nnz+i) = i
icolzh(nnz+i) = i
azh(nnz+i) = conjg(ze)
End Do
! Sort the elements into correct coordinate storage format
Call f11znf(n,nnzzh,azh,irowzh,icolzh,'S','R',istrh,iwork,ifail)
! Form incomplete LU factorization of (ze I - A)^H
Call f11dnf(n,nnzzh,azh,la,irowzh,icolzh,0,0.0_nag_wp,'P','N', &
ipivph,ipivqh,istrh,idiagh,nnzch,npivm,iwork,liwork,ifail)
If (ifail/=0) Then
Exit revcomm
End If
Cycle revcomm
Case (4)
! Solve the linear system (ze I - A)^H w = y, with m0 righthand sides
Do i = 1, m0
! We will need to store the final result in y
w(1:n) = y(1:n,i)
! Initial guess
y(1:n,i) = cone
Call f11dqf('RGMRES',n,nnzzh,azh,la,irowzh,icolzh,ipivph,ipivqh, &
istrh,idiagh,w,m,tol,500,y(1,i),rnorm,itn,work,lwork,ifail)
End Do
If (ifail/=0) Then
Exit revcomm
End If
Cycle revcomm
Case (5)
! Compute x <- Az
Do i = 1, m0
Call f11xaf('N',n,nnz,a,irow,icol,'N',z(1,i),x(1,i),ifail)
End Do
If (ifail/=0) Then
Exit revcomm
End If
Cycle revcomm
Case (6)
! Compute x <- A^T z
Do i = 1, m0
Call f11xaf('T',n,nnz,a,irow,icol,'N',z(1,i),x(1,i),ifail)
End Do
If (ifail/=0) Then
Exit revcomm
End If
Cycle revcomm
Case (7)
! Since we are not solving a generalized eigenvalue problem set x = z
x(1:n,1:m0) = z(1:n,1:m0)
Cycle revcomm
Case (8)
! Since we are not solving a generalized eigenvalue problem set x = z
x(1:n,1:m0) = z(1:n,1:m0)
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,y,ldz,'Right eigenvectors',ifail)
Call x04daf('General',' ',n,nconv,y(1,m0+1),ldz,'Left eigenvectors', &
ifail)
Else
Write (nout,99998) 'Failure in f12jkf 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 f12jkfe