Program f04ydfe ! f04ydf Example Program Text ! Mark 24 Release. NAG Copyright 2012. ! .. Implicit None Statement .. Implicit None ! .. Parameters .. Integer, Parameter :: nout = 6 ! .. Executable Statements .. Write (nout,*) 'F04YDF Example Program Results' Call ex1 Call ex2 Contains Subroutine ex1 ! .. Use Statements .. Use nag_library, Only: dgetrf, dgetrs, f04ydf, f06raf, nag_wp ! .. Implicit None Statement .. Implicit None ! .. Parameters .. Integer, Parameter :: nin = 5, nout = 6 ! .. Local Scalars .. Real (Kind=nag_wp) :: cond, nrma, nrminv Integer :: i, ifail, irevcm, lda, ldx, ldy, & m, n, seed, t ! .. Local Arrays .. Real (Kind=nag_wp), Allocatable :: a(:,:), work(:), x(:,:), y(:,:) Integer, Allocatable :: ipiv(:), iwork(:) ! .. Executable Statements .. Write (nout,*) 'Example 1' Write (nout,*) ! Skip heading in data file Read (nin,'(///A)') Read (nin,*) m, n, t lda = m ldx = n ldy = m Allocate (a(lda,n),x(ldx,t),y(ldy,t),work(m*t),iwork(2*n+5*t+20), & ipiv(n)) ! Read A from data file Read (nin,*)(a(i,1:n),i=1,m) ! Compute 1-norm of A nrma = f06raf('1',m,n,a,lda,work) Write (nout,99999) 'The norm of A is: ', nrma ! Estimate the norm of A^(-1) without explicitly forming A^(-1) ! Perform an LU factorization so that A=LU where L and U are lower ! and upper triangular. ifail = 0 Call dgetrf(m,n,a,lda,ipiv,ifail) seed = 354 irevcm = 0 loop: Do Call f04ydf(irevcm,m,n,x,ldx,y,ldy,nrminv,t,seed,work,iwork,ifail) If (irevcm==0) Then Exit loop Else If (irevcm==1) Then ! Compute y = inv(A)*x Call dgetrs('N',n,t,a,lda,ipiv,x,ldx,ifail) ! x was overwritten by dgetrs, so set y=x y(1:n,1:t) = x(1:n,1:t) Else ! Compute x = transpose(inv(A))*y Call dgetrs('T',n,t,a,lda,ipiv,y,ldy,ifail) ! y was overwritten by dgetrs so set x=y x(1:n,1:t) = y(1:n,1:t) End If End Do loop Write (nout,99999) 'The estimated norm of inverse(A) is: ', nrminv ! Compute and print the estimated condition number cond = nrminv*nrma Write (nout,*) Write (nout,99999) 'Estimated condition number of A: ', cond Write (nout,*) 99999 Format (1X,A,F6.2) End Subroutine ex1 Subroutine ex2 ! .. Use Statements .. Use nag_library, Only: f01brf, f04axf, f04ydf, nag_wp ! .. Implicit None Statement .. Implicit None ! .. Parameters .. Real (Kind=nag_wp), Parameter :: zero = 0.0_nag_wp ! .. Local Scalars .. Real (Kind=nag_wp) :: asum, cond, nrma, nrminv, pivot, & resid Integer :: i, ifail, irevcm, j, ldx, ldy, & licn, lirn, n, nin, nout, nz, & seed, t Logical :: grow, lblock ! .. Local Arrays .. Real (Kind=nag_wp), Allocatable :: a(:), w(:), work(:), x(:,:), y(:,:) Integer, Allocatable :: icn(:), ikeep(:), irn(:), iw(:), & iwork(:) Integer :: idisp(10) Logical :: abort(4) ! .. Intrinsic Procedures .. Intrinsic :: abs, max ! .. Executable Statements .. Continue nout = 6 nin = 5 Write (nout,'(//1X,A/)') 'Example 2' ! Skip heading in data file Read (nin,'(//A)') ! Input N, the order of the matrix A, and NZ the number of non-zero ! elements of A, together with t the norm estimation parameter. Read (nin,*) n, nz, t licn = 4*nz lirn = 2*nz ldx = n ldy = n ! Allocate the required memory Allocate (a(licn),icn(licn),irn(lirn),x(ldx,t),y(ldy,t),work(n*t), & ikeep(5*n),iwork(2*n+5*t+20),iw(8*n),w(n)) ! Input the elements of A, along with row and column information. Read (nin,*)(a(i),irn(i),icn(i),i=1,nz) ! Compute 1-norm of A nrma = zero Do i = 1, n asum = zero Do j = 1, nz If (icn(j)==i) asum = asum + abs(a(j)) End Do nrma = max(nrma,asum) End Do Write (nout,99999) 'The norm of A is: ', nrma ! Perform an LU factorization so that A=LU where L and U are lower ! and upper triangular, using F01BRF pivot = 0.1_nag_wp grow = .True. lblock = .True. abort(1) = .True. abort(2) = .True. abort(3) = .False. abort(4) = .True. ifail = 0 Call f01brf(n,nz,a,licn,irn,lirn,icn,pivot,ikeep,iw,w,lblock,grow, & abort,idisp,ifail) ! Compute and estimate of the 1-norm of inv(A) seed = 412 irevcm = 0 loop: Do Call f04ydf(irevcm,n,n,x,ldx,y,ldy,nrminv,t,seed,work,iwork,ifail) If (irevcm==0) Then Exit loop Else If (irevcm==1) Then ! Compute y = inv(A)*x Do i = 1, t Call f04axf(n,a,licn,icn,ikeep,x(1,i),w,irevcm,idisp,resid) End Do ! x was overwritten by f04axf, so set y=x y(1:n,1:t) = x(1:n,1:t) Else ! Compute x = transpose(inv(A))*y Do i = 1, t Call f04axf(n,a,licn,icn,ikeep,y(1,i),w,irevcm,idisp,resid) End Do ! y was overwritten by f04axf so set x=y x(1:n,1:t) = y(1:n,1:t) End If End Do loop Write (nout,99999) 'The estimated norm of inverse(A) is: ', nrminv ! Compute and print the estimated condition number cond = nrminv*nrma Write (nout,*) Write (nout,99999) 'Estimated condition number of A: ', cond Write (nout,*) 99999 Format (1X,A,F6.2) End Subroutine ex2 End Program f04ydfe