Program f12jvfe
! F12JVF Example Program Text
! Mark 29.2 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, f12jgf, &
f12jvf, 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, deg, i, ifail, ind, irevcm, &
iter, itn, j, k, la, ldx, ldy, ldz, &
liwork, lwork, m, m0, n, nconv, &
nnzc, nnzch, nnzmax, nnzsum, 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(:), nnza(:), &
tedge(:)
! .. Intrinsic Procedures ..
Intrinsic :: conjg, max, maxval, min, sqrt, sum
! .. Executable Statements ..
Write (nout,*) 'F12JVF Example Program Results'
Write (nout,*)
! Skip heading in data file
Read (nin,*)
Read (nin,*) n, deg
Allocate (nnza(deg+1))
Read (nin,*) nnza(1:deg+1)
nnzmax = maxval(nnza)
nnzsum = sum(nnza)
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*(nnzsum)
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(nnzmax,deg+1),w(n),icol(nnzmax,deg+1),irow(nnzmax,deg+1),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 matrices from data file
Do j = 1, deg + 1
Do i = 1, nnza(j)
Read (nin,*) a(i,j), irow(i,j), icol(i,j)
End Do
End Do
! Initialize the data handle
ifail = 0
Call f12jaf(handle,ifail)
! Define arrays for the custom contour
ccn = 2
Allocate (nedge(ccn),tedge(ccn),zedge(ccn))
zedge(1) = (-1.5_nag_wp,-0.5_nag_wp)
zedge(2) = (0.5_nag_wp,1.5_nag_wp)
tedge(1) = 100
tedge(2) = 0
nedge(1) = 20
nedge(2) = 10
! Generate the contour
Call f12jgf(handle,ccn,nedge,tedge,zedge,ifail)
irevcm = 0
revcomm: Do
Call f12jvf(handle,irevcm,deg,ze,n,x,ldx,y,ldy,k,m0,nconv,d,z,ldz,eps, &
iter,resid,ifail)
Select Case (irevcm)
Case (1)
! Form the sparse matrix \sum ze^i A_i
nnzz = nnzsum
ind = 1
Do i = 0, deg
az(ind:ind+nnza(i+1)-1) = (ze**i)*a(1:nnza(i+1),i+1)
irowz(ind:ind+nnza(i+1)-1) = irow(1:nnza(i+1),i+1)
icolz(ind:ind+nnza(i+1)-1) = icol(1:nnza(i+1),i+1)
ind = ind + nnza(i+1)
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 \sum ze^i A_i
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 (\sum ze^i A_i) 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 (\sum ze^i A_i)^H
nnzzh = nnzsum
ind = 1
Do i = 0, deg
azh(ind:ind+nnza(i+1)-1) = conjg((ze**i)*a(1:nnza(i+1),i+1))
irowzh(ind:ind+nnza(i+1)-1) = icol(1:nnza(i+1),i+1)
icolzh(ind:ind+nnza(i+1)-1) = irow(1:nnza(i+1),i+1)
ind = ind + nnza(i+1)
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 (\sum ze^i A_i)^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 (\sum ze^i A_i)^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 <- A_k z
Do i = 1, m0
Call f11xnf('N',n,nnza(k+1),a(1,k+1),irow(1,k+1),icol(1,k+1),'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_k^H z
Do i = 1, m0
Call f11xnf('T',n,nnza(k+1),a(1,k+1),irow(1,k+1),icol(1,k+1),'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)
Call x04daf('General',' ',n,nconv,z(1,m0+1),ldz,'Left eigenvectors', &
ifail)
Else
Write (nout,99998) 'Failure in f12jvf IFAIL =', ifail
End If
! Destroy data handle
Call f12jzf(handle,ifail)
99999 Format ((3X,4(' (',F7.4,',',F7.4,')',:)))
99998 Format (1X,A,I4)
End Program f12jvfe