Program f04ycfe ! F04YCF Example Program Text ! Mark 24 Release. NAG Copyright 2012. ! .. 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, dgetrf, dgetrs, f04ycf, nag_wp ! .. Parameters .. Real (Kind=nag_wp), Parameter :: zero = 0.0E+0_nag_wp Integer, Parameter :: nin = 5 ! .. Local Scalars .. Real (Kind=nag_wp) :: anorm, cond, estnrm Integer :: i, icase, ifail, info, j, lda, n ! .. Local Arrays .. Real (Kind=nag_wp), Allocatable :: a(:,:), p(:), work(:), x(:) Integer, Allocatable :: ipiv(:), iwork(:) ! .. Intrinsic Procedures .. 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),ipiv(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 ! LU Factorize A ! The NAG name equivalent of dgetrf is f07adf Call dgetrf(n,n,a,lda,ipiv,info) 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(A)*X ! The NAG name equivalent of dgetrs is f07aef Call dgetrs('N',n,1,a,lda,ipiv,x,n,info) Else If (icase==2) Then ! Return the vector inv(A)'*X Call dgetrs('T',n,1,a,lda,ipiv,x,n,info) 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 ! .. 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 Procedures .. 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