Program f12jrfe
! F12JRF Example Program Text
! Mark 28.5 Release. NAG Copyright 2022.
! Intrinsic Procedures ..
! .. Use Statements ..
Use, Intrinsic :: iso_c_binding, Only: c_ptr
Use nag_library, Only: f11dnf, f11dqf, f11xsf, f11znf, f12jaf, f12jbf, &
f12jef, f12jrf, 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) :: emax, emin, eps, rnorm, tol
Integer :: i, ifail, irevcm, iter, itn, la, &
ldx, ldy, ldz, liwork, lwork, m, m0, &
n, nconv, nnza, nnzb, nnzc, nnzch, &
nnzz, nnzzh, npivm
! .. Local Arrays ..
Complex (Kind=nag_wp), Allocatable :: a(:), az(:), azh(:), b(:), w(:), &
work(:), x(:,:), y(:,:), z(:,:)
Real (Kind=nag_wp), Allocatable :: d(:), resid(:)
Integer, Allocatable :: icola(:), icolb(:), icolz(:), &
icolzh(:), idiag(:), idiagh(:), &
ipiv(:), ipivp(:), ipivph(:), &
ipivq(:), ipivqh(:), irowa(:), &
irowb(:), irowz(:), irowzh(:), &
istr(:), istrh(:), iwork(:)
! .. Intrinsic Procedures ..
Intrinsic :: conjg, max, min, sqrt
! .. Executable Statements ..
Write (nout,*) 'F12JRF Example Program Results'
Write (nout,*)
! Skip heading in data file
Read (nin,*)
Read (nin,*) n
Read (nin,*) nnza, nnzb
m0 = n
ldx = n
ldy = n
ldz = n
m = min(n,50)
lwork = max(64*n,4*n+m*(m+n+5)+121)
la = 4*(nnza+nnzb)
liwork = 7*n + 2
tol = sqrt(x02ajf())
Allocate (x(ldx,m0),y(ldy,m0),z(ldz,m0),resid(m0),d(m0),a(nnza),b(nnzb), &
w(n),icola(nnza),irowa(nnza),icolb(nnzb),irowb(nnzb),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, nnza
Read (nin,*) a(i), irowa(i), icola(i)
End Do
! Read the matrix B from data file
Do i = 1, nnzb
Read (nin,*) b(i), irowb(i), icolb(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 = Zolotarev', &
ifail=ifail)
emin = -1.0_nag_wp
emax = 0.0_nag_wp
! Generate the contour
Call f12jef(handle=handle,emin=emin,emax=emax,ifail=ifail)
irevcm = 0
revcomm: Do
Call f12jrf(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 B - A
nnzz = 2*(nnza+nnzb)
! We store it in the arrays az, irowz and icolz
az(1:nnza) = -a(1:nnza)
irowz(1:nnza) = irowa(1:nnza)
icolz(1:nnza) = icola(1:nnza)
az(nnza+1:2*nnza) = -a(1:nnza)
irowz(nnza+1:2*nnza) = icola(1:nnza)
icolz(nnza+1:2*nnza) = irowa(1:nnza)
! Add the elements of ze B
irowz(2*nnza+1:2*nnza+nnzb) = irowb(1:nnzb)
icolz(2*nnza+1:2*nnza+nnzb) = icolb(1:nnzb)
az(2*nnza+1:2*nnza+nnzb) = ze*b(1:nnzb)
irowz(2*nnza+nnzb+1:2*nnza+2*nnzb) = icolb(1:nnzb)
icolz(2*nnza+nnzb+1:2*nnza+2*nnzb) = irowb(1:nnzb)
az(2*nnza+nnzb+1:2*nnza+2*nnzb) = conjg(ze)*b(1:nnzb)
! Sort the elements into correct coordinate storage format
Call f11znf(n,nnzz,az,irowz,icolz,'S','R',istr,iwork,ifail)
! We have double counted the diagonal elements, so halve them
Do i = 1, nnzz
If (irowz(i)==icolz(i)) Then
az(i) = 0.5_nag_wp*az(i)
End If
End Do
! Form incomplete LU factorization of ze B - 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 B - A)^H
nnzzh = 2*(nnza+nnzb)
! We store it in the arrays azh, irowzh and icolzh
azh(1:nnza) = -a(1:nnza)
irowzh(1:nnza) = irowa(1:nnza)
icolzh(1:nnza) = icola(1:nnza)
azh(nnza+1:2*nnza) = -a(1:nnza)
irowzh(nnza+1:2*nnza) = icola(1:nnza)
icolzh(nnza+1:2*nnza) = irowa(1:nnza)
! Add the elements of (ze B)^H
irowzh(2*nnza+1:2*nnza+nnzb) = irowb(1:nnzb)
icolzh(2*nnza+1:2*nnza+nnzb) = icolb(1:nnzb)
azh(2*nnza+1:2*nnza+nnzb) = conjg(ze)*b(1:nnzb)
irowzh(2*nnza+nnzb+1:2*nnza+2*nnzb) = icolb(1:nnzb)
icolzh(2*nnza+nnzb+1:2*nnza+2*nnzb) = irowb(1:nnzb)
azh(2*nnza+nnzb+1:2*nnza+2*nnzb) = ze*b(1:nnzb)
! Sort the elements into correct coordinate storage format
Call f11znf(n,nnzzh,azh,irowzh,icolzh,'S','R',istr,iwork,ifail)
! We have double counted the diagonal elements, so halve them
Do i = 1, nnzzh
If (irowzh(i)==icolzh(i)) Then
azh(i) = 0.5_nag_wp*azh(i)
End If
End Do
! Form incomplete LU factorization of (ze B - 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 f11xsf(n,nnza,a,irowa,icola,'N',z(1,i),x(1,i),ifail)
End Do
If (ifail/=0) Then
Exit revcomm
End If
Cycle revcomm
Case (6)
! Compute x <- Bz
Do i = 1, m0
Call f11xsf(n,nnzb,b,irowb,icolb,'N',z(1,i),x(1,i),ifail)
End Do
If (ifail/=0) Then
Exit revcomm
End If
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)
Else
Write (nout,99998) 'Failure in f12jrf IFAIL =', ifail
End If
! Destroy the data handle
Call f12jzf(handle=handle,ifail=ifail)
99999 Format (3X,(8F8.4))
99998 Format (1X,A,I4)
End Program f12jrfe