Program f12jtfe
! F12JTF Example Program Text
! Mark 29.0 Release. NAG Copyright 2023.
! Intrinsic Procedures ..
! .. Use Statements ..
Use, Intrinsic :: iso_c_binding, Only: c_ptr
Use nag_library, Only: f11dnf, f11dqf, f11xnf, f11znf, f12jaf, f12jbf, &
f12jgf, f12jtf, 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) :: ze
Real (Kind=nag_wp) :: eps, rnorm, tol
Integer :: ccn, 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 :: a(:), az(:), azh(:), d(:), w(:), &
work(:), x(:,:), y(:,:), z(:,:), &
zedge(:)
Real (Kind=nag_wp), Allocatable :: resid(:)
Integer, Allocatable :: icol(:), icolz(:), icolzh(:), &
idiag(:), idiagh(:), ipiv(:), &
ipivp(:), ipivph(:), ipivq(:), &
ipivqh(:), irow(:), irowz(:), &
irowzh(:), istr(:), istrh(:), &
iwork(:), nedge(:), tedge(:)
! .. Intrinsic Procedures ..
Intrinsic :: conjg, max, min, sqrt
! .. Executable Statements ..
Write (nout,*) 'F12JTF 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,2*m0),y(ldy,m0),z(ldz,2*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=handle,ifail=ifail)
! Set the integration type
Call f12jbf(handle=handle,str='Integration Type = Trapezoidal', &
ifail=ifail)
! Define arrays for the custom contour
ccn = 4
Allocate (nedge(ccn),tedge(ccn),zedge(ccn))
zedge(1) = (0.0_nag_wp,1.0_nag_wp)
zedge(2) = (0.0_nag_wp,-1.0_nag_wp)
zedge(3) = (-1.0_nag_wp,-1.0_nag_wp)
zedge(4) = (-1.0_nag_wp,1.0_nag_wp)
tedge(1) = 50
tedge(2) = 0
tedge(3) = 0
tedge(4) = 0
nedge(1) = 20
nedge(2) = 1
nedge(3) = 1
nedge(4) = 1
! Generate the contour
Call f12jgf(handle=handle,ccn=ccn,nedge=nedge,tedge=tedge,zedge=zedge, &
ifail=ifail)
irevcm = 0
revcomm: Do
Call f12jtf(handle=handle,irevcm=irevcm,ze=ze,n=n,x=x,ldx=ldx,y=y, &
ldy=ldy,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 sparse matrix ze I - A
nnzz = nnz + n
az(1:nnz) = -a(1:nnz)
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)y = w, 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) = -conjg(a(1:nnz))
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 y = w, 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 f11xnf('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^H z
Do i = 1, m0
Call f11xnf('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 .Or. ifail==14) Then
! Print solution
Write (nout,*) 'Eigenvalues'
Write (nout,99999) d(1:nconv)
Flush (nout)
Call x04daf('General',' ',n,nconv,z,ldz,'Right eigenvectors',ifail)
Call x04daf('General',' ',n,nconv,z(1,m0+1),ldz,'Left eigenvectors', &
ifail)
Else
Write (nout,99998) 'Failure in f12jtf 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 f12jtfe