PROGRAM f04ycfe ! F04YCF Example Program Text ! Mark 23 Release. NAG Copyright 2011. ! .. Implicit None Statement .. IMPLICIT NONE ! .. Parameters .. INTEGER, PARAMETER :: nout = 6 ! .. Executable Statements .. WRITE (nout,*) 'F04YCF Example Program Results' CALL ex1 CALL ex2 CONTAINS SUBROUTINE ex1 ! .. Use Statements .. USE nag_library, ONLY : dasum, dtrsv, f03aff, f04ycf, nag_wp ! .. Implicit None Statement .. IMPLICIT NONE ! .. Parameters .. REAL (KIND=nag_wp), PARAMETER :: zero = 0.0E+0_nag_wp INTEGER, PARAMETER :: nin = 5 ! .. Local Scalars .. REAL (KIND=nag_wp) :: anorm, cond, d1, eps, estnrm INTEGER :: i, icase, id, ifail, j, lda, n ! .. Local Arrays .. REAL (KIND=nag_wp), ALLOCATABLE :: a(:,:), p(:), work(:), x(:) INTEGER, ALLOCATABLE :: iwork(:) ! .. Intrinsic Functions .. INTRINSIC max ! .. Executable Statements .. WRITE (nout,*) WRITE (nout,*) WRITE (nout,*) 'Example 1' WRITE (nout,*) ! Skip heading in data file READ (nin,*) READ (nin,*) READ (nin,*) READ (nin,*) n lda = n ALLOCATE (a(lda,n),p(n),work(n),x(n),iwork(n)) READ (nin,*) (a(i,1:n),i=1,n) ! First compute the norm of A. ! DASUM (f06ekf) returns the sum of the absolute values of a column of A. anorm = zero DO j = 1, n anorm = max(anorm,dasum(n,a(1,j),1)) END DO WRITE (nout,99999) 'Computed norm of A =', anorm ! Next estimate the norm of inverse(A). We do not form the ! inverse explicitly. ! ifail: behaviour on error exit ! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft ifail = 0 ! Factorise A as P*A = L*U using F03AFF. CALL f03aff(n,eps,a,lda,d1,id,p,ifail) icase = 0 LOOP: DO ifail = 0 CALL f04ycf(icase,n,x,estnrm,work,iwork,ifail) IF (icase/=0) THEN IF (icase==1) THEN ! Return the vector inv(P*A)*X by solving the equations ! L*U*Y = X, overwriting Y on X. First solve L*Z = X for Z. ! The NAG name equivalent of dtrsv is f06pjf CALL dtrsv('Lower','No Transpose','Non-Unit',n,a,lda,x,1) ! Then solve U*Y = Z for Y. CALL dtrsv('Upper','No Transpose','Unit',n,a,lda,x,1) ELSE IF (icase==2) THEN ! Return the vector inv(P*A)'*X by solving U'*L'*Y = X, ! overwriting Y on X. First solve U'*Z = X for Z. CALL dtrsv('Upper','Transpose','Unit',n,a,lda,x,1) ! Then solve L'*Y = Z for Y. CALL dtrsv('Lower','Transpose','Non-Unit',n,a,lda,x,1) END IF ! Continue until icase is returned as 0. ELSE WRITE (nout,99999) 'Estimated norm of inverse(A) =', estnrm cond = anorm*estnrm WRITE (nout,99998) 'Estimated condition number of A =', cond WRITE (nout,*) EXIT LOOP END IF END DO LOOP 99999 FORMAT (1X,A,F8.4) 99998 FORMAT (1X,A,F5.1) END SUBROUTINE ex1 SUBROUTINE ex2 ! .. Use Statements .. USE nag_library, ONLY : f01brf, f04axf, f04ycf, nag_wp ! .. Implicit None Statement .. IMPLICIT NONE ! .. Parameters .. REAL (KIND=nag_wp), PARAMETER :: tenth = 0.1E+0_nag_wp REAL (KIND=nag_wp), PARAMETER :: zero = 0.0E+0_nag_wp INTEGER, PARAMETER :: nin = 5 ! .. Local Scalars .. REAL (KIND=nag_wp) :: anorm, cond, estnrm, resid, sum, u INTEGER :: i, icase, ifail, j, licn, lirn, & n, nz LOGICAL :: grow, lblock ! .. Local Arrays .. REAL (KIND=nag_wp), ALLOCATABLE :: a(:), w(:), work1(:), x(:) INTEGER, ALLOCATABLE :: icn(:), ikeep(:), irn(:), iw(:), & iwork(:) INTEGER :: idisp(10) LOGICAL :: abort(4) ! .. Intrinsic Functions .. INTRINSIC abs, max ! .. Executable Statements .. WRITE (nout,*) WRITE (nout,*) WRITE (nout,*) 'Example 2' WRITE (nout,*) ! Skip heading in data file READ (nin,*) READ (nin,*) ! Input N, the order of matrix A, and NZ, the number of non-zero ! elements of A. READ (nin,*) n, nz licn = 4*nz lirn = 2*nz ALLOCATE (a(licn),w(n),work1(n),x(n),icn(licn),ikeep(5*n),irn(lirn), & iw(8*n),iwork(n)) ! Input the elements of A, along with row and column information. READ (nin,*) (a(i),irn(i),icn(i),i=1,nz) ! First compute the norm of A. anorm = zero DO i = 1, n sum = zero DO j = 1, nz IF (icn(j)==i) sum = sum + abs(a(j)) END DO anorm = max(anorm,sum) END DO WRITE (nout,99999) 'Computed norm of A =', anorm ! Next estimate the norm of inverse(A). We do not form the ! inverse explicitly. ! Factorise A into L*U using F01BRF. u = tenth lblock = .TRUE. grow = .TRUE. abort(1) = .TRUE. abort(2) = .TRUE. abort(3) = .FALSE. abort(4) = .TRUE. ! ifail: behaviour on error exit ! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft ifail = 0 CALL f01brf(n,nz,a,licn,irn,lirn,icn,u,ikeep,iw,w,lblock,grow,abort, & idisp,ifail) icase = 0 LOOP: DO ifail = 0 CALL f04ycf(icase,n,x,estnrm,work1,iwork,ifail) IF (icase/=0) THEN ! Return X := inv(A)*X or X = inv(A)'*X, depending on the ! value of ICASE, by solving A*Y = X or A'*Y = X, ! overwriting Y on X. CALL f04axf(n,a,licn,icn,ikeep,x,w,icase,idisp,resid) ! Continue until icase is returned as 0. ELSE WRITE (nout,99999) 'Estimated norm of inverse(A) =', estnrm cond = anorm*estnrm WRITE (nout,99998) 'Estimated condition number of A =', cond EXIT LOOP END IF END DO LOOP 99999 FORMAT (1X,A,F8.4) 99998 FORMAT (1X,A,F5.1) END SUBROUTINE ex2 END PROGRAM f04ycfe