! F02FJF Example Program Text
! Mark 26.1 Release. NAG Copyright 2016.
Module f02fjfe_mod
! F02FJF Example Program Module:
! Parameters and User-defined Routines
! .. Use Statements ..
Use nag_library, Only: nag_wp
! .. Implicit None Statement ..
Implicit None
! .. Accessibility Statements ..
Private
Public :: dot, image, monit
! .. Parameters ..
Real (Kind=nag_wp), Parameter :: half = 0.5_nag_wp
Real (Kind=nag_wp), Parameter, Public :: zero = 0.0_nag_wp
Integer, Parameter, Public :: nin = 5, nout = 6
Contains
Function dot(iflag,n,z,w,ruser,lruser,iuser,liuser)
! This function implements the dot product - transpose(W)*B*Z.
! DOT assumes that N is at least 3.
! .. Function Return Value ..
Real (Kind=nag_wp) :: dot
! .. Scalar Arguments ..
Integer, Intent (Inout) :: iflag
Integer, Intent (In) :: liuser, lruser, n
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Inout) :: ruser(lruser)
Real (Kind=nag_wp), Intent (In) :: w(n), z(n)
Integer, Intent (Inout) :: iuser(liuser)
! .. Local Scalars ..
Real (Kind=nag_wp) :: s
Integer :: i
! .. Executable Statements ..
s = zero
s = s + (z(1)-half*z(2))*w(1)
s = s + (-half*z(n-1)+z(n))*w(n)
Do i = 2, n - 1
s = s + (-half*z(i-1)+z(i)-half*z(i+1))*w(i)
End Do
dot = s
! Set iflag negative to terminate execution for any reason.
iflag = 0
Return
End Function dot
Subroutine image(iflag,n,z,w,ruser,lruser,iuser,liuser)
! This routine solves A*W = B*Z for W.
! The routine assumes that N is at least 3.
! The data A, NNZ, LA, IROW, ICOL, IPIV and ISTR on exit from
! F11JAF have been packed into the xUSER communication arrays in
! the following way:
! IUSER(1:2) = (/NNZ, LA/)
! RUSER(1:LA) = A
! IUSER(3:(2*LA+2*N+3)) = (/IROW, ICOL, IPIV, ISTR/)
! We'll also use RUSER((LA+1):(LA+N)) as space for F11JCF's dummy
! arg. B, and RUSER((LA+N+1):(LA+7*N+120)) as space for F11JCF's
! dummy arg. WORK
! .. Use Statements ..
Use nag_library, Only: f11jcf, x02ajf
! .. Scalar Arguments ..
Integer, Intent (Inout) :: iflag
Integer, Intent (In) :: liuser, lruser, n
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Inout) :: ruser(lruser)
Real (Kind=nag_wp), Intent (Out) :: w(n)
Real (Kind=nag_wp), Intent (In) :: z(n)
Integer, Intent (Inout) :: iuser(liuser)
! .. Local Scalars ..
Real (Kind=nag_wp) :: rnorm, tol
Integer :: ifail, itn, j, la, lwork, maxitn, &
nnz
Character (2) :: method
! .. Executable Statements ..
nnz = iuser(1)
la = iuser(2)
! Form B*Z in RUSER((LA+1):(LA+N)) and initialize W to
! zero.
w(1:n) = zero
ruser(la+1) = z(1) - half*z(2)
Do j = 2, n - 1
ruser(la+j) = -half*z(j-1) + z(j) - half*z(j+1)
End Do
ruser(la+n) = -half*z(n-1) + z(n)
! Call F11JCF to solve the equations A*W = B*Z.
method = 'CG'
tol = x02ajf()
maxitn = 100
lwork = 6*n + 120
! ifail: behaviour on error exit
! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
ifail = 1
Call f11jcf(method,n,nnz,ruser,la,iuser(3),iuser(la+3),iuser(2*la+3), &
iuser(2*la+n+3),ruser(la+1),tol,maxitn,w,rnorm,itn,ruser(la+n+1), &
lwork,ifail)
If (ifail>0) Then
iflag = -ifail
End If
Return
End Subroutine image
Subroutine monit(istate,nextit,nevals,nevecs,k,f,d)
! Monitoring routine for F02FJF.
! .. Parameters ..
Integer, Parameter :: nout = 6
! .. Scalar Arguments ..
Integer, Intent (In) :: istate, k, nevals, nevecs, nextit
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (In) :: d(k), f(k)
! .. Local Scalars ..
Integer :: i
! .. Executable Statements ..
If (istate/=0) Then
Write (nout,*)
Write (nout,99999) ' ISTATE = ', istate, ' NEXTIT = ', nextit
Write (nout,99999) ' NEVALS = ', nevals, ' NEVECS = ', nevecs
Write (nout,*) ' F D'
Write (nout,99998)(f(i),d(i),i=1,k)
End If
Return
99999 Format (1X,A,I4,A,I4)
99998 Format (1X,1P,E11.3,3X,E11.3)
End Subroutine monit
End Module f02fjfe_mod
Program f02fjfe
! F02FJF Example Main Program
! .. Use Statements ..
Use f02fjfe_mod, Only: dot, image, monit, nin, nout, zero
Use nag_library, Only: f02fjf, f06fef, f11jaf, nag_wp, x04cbf
! .. Implicit None Statement ..
Implicit None
! .. Local Scalars ..
Real (Kind=nag_wp) :: dscale, dtol, tol
Integer :: i, ifail, k, l, la, ldx, lfill, &
liuser, lruser, lwork, m, n, nnz, &
nnzc, noits, novecs, npivm
Character (1) :: mic, pstrat
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: a(:), d(:), ruser(:), work(:), &
x(:,:)
Integer, Allocatable :: icol(:), ipiv(:), irow(:), istr(:), &
iuser(:)
Character (1) :: clabs(1), rlabs(1)
! .. Executable Statements ..
Write (nout,*) 'F02FJF Example Program Results'
Write (nout,*)
Flush (nout)
! Skip heading in data file
Read (nin,*)
Read (nin,*) n, m, k, tol
la = 10*n
ldx = n
liuser = 2*la + 2*n + 3
lruser = la + 7*n + 120
lwork = 5*k + 2*n
Allocate (a(la),d(n),ruser(lruser),work(lwork),x(ldx,k),icol(la), &
ipiv(n),irow(la),istr(n+1),iuser(liuser))
! Set up the sparse symmetric coefficient matrix A.
l = 0
Do i = 1, n
If (i>=5) Then
l = l + 1
a(l) = -0.25_nag_wp
irow(l) = i
icol(l) = i - 4
End If
If (i>=2) Then
l = l + 1
a(l) = -0.25_nag_wp
irow(l) = i
icol(l) = i - 1
End If
l = l + 1
a(l) = 1.0_nag_wp
irow(l) = i
icol(l) = i
End Do
nnz = l
! Call F11JAF to find an incomplete Cholesky factorization of A.
lfill = 2
dtol = zero
mic = 'M'
dscale = zero
pstrat = 'M'
! ifail: behaviour on error exit
! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
ifail = 0
Call f11jaf(n,nnz,a,la,irow,icol,lfill,dtol,mic,dscale,pstrat,ipiv,istr, &
nnzc,npivm,iuser,liuser,ifail)
! Call F02FJF to find eigenvalues and eigenvectors.
noits = 1000
novecs = 0
! Communicate A, NNZ, LA, IROW, ICOL, IPIV and ISTR to IMAGE
! thread-safely using RUSER and IUSER.
! In addition to using RUSER for storing A, we'll also use
! 7*N+120 elements of RUSER in place of local arrays in IMAGE.
! Initialized A goes into ruser.
ruser(1:nnz+nnzc) = a(1:nnz+nnzc)
! NNZ, LA, IROW, ICOL, IPIV and ISTR go into IUSER, in that order.
! Only the first NNZ+NNZC elements of IROW and ICOL will have been
! initialized:
iuser(1) = nnz
iuser(2) = la
iuser(3:2+nnz+nnzc) = irow(1:nnz+nnzc)
iuser(la+3:la+2+nnz+nnzc) = icol(1:nnz+nnzc)
iuser(2*la+3:2*la+2+n) = ipiv(1:n)
iuser(liuser-n:liuser) = istr(1:n+1)
ifail = -1
Call f02fjf(n,m,k,noits,tol,dot,image,monit,novecs,x,ldx,d,work,lwork, &
ruser,lruser,iuser,liuser,ifail)
If (ifail>=0) Then
If (ifail/=1 .And. ifail<=4 .And. m>=1) Then
Do i = 1, m
d(i) = 1.0_nag_wp/d(i)
End Do
Write (nout,*) 'Final results'
Write (nout,*)
Write (nout,*) ' Eigenvalues'
Write (nout,99999) d(1:m)
Write (nout,*)
Flush (nout)
! Normalize eigenvectors
Do i = 1, m
Call f06fef(n,x(1,i),x(1,i),1)
End Do
Call x04cbf('General',' ',n,m,x,ldx,'1P,E12.3',' Eigenvectors','N', &
rlabs,'N',clabs,80,0,ifail)
End If
End If
99999 Format (1X,1P,4E12.3)
End Program f02fjfe