! F02FJF Example Program Text ! Mark 23 Release. NAG Copyright 2011. 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 ! .. Parameters .. REAL (KIND=nag_wp), PARAMETER :: half = 0.5_nag_wp REAL (KIND=nag_wp), PARAMETER :: zero = 0.0_nag_wp INTEGER, PARAMETER :: 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. ! .. Implicit None Statement .. IMPLICIT NONE ! .. 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 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 ! .. Implicit None Statement .. IMPLICIT NONE ! .. 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) iflag = -ifail RETURN END SUBROUTINE image SUBROUTINE monit(istate,nextit,nevals,nevecs,k,f,d) ! Monitoring routine for F02FJF. ! .. Implicit None Statement .. IMPLICIT NONE ! .. 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 nag_library, ONLY : f02fjf, f02fjz, f11jaf, x04cbf USE f02fjfe_mod, ONLY : dot, image, nag_wp, nin, nout, zero ! .. 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 factorisation 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) ! * To obtain monitoring information from the supplied subroutine MONIT, ! replace the name F02FJZ by MONIT in the call to F02FJF, and USE MONIT ! from f02fjfe_mod * ifail = -1 CALL f02fjf(n,m,k,noits,tol,dot,image,f02fjz,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 x(1:n,i) = x(1:n,i)/x(1,i) 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