NAG Library Manual, Mark 28.4
Interfaces:  FL   CL   CPP   AD
```    Program f04ydfe

!     F04YDF Example Program Text

!     Mark 28.4 Release. NAG Copyright 2022.

!     .. Use Statements ..
Use nag_library, Only: nag_wp
!     .. 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
!       .. 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,*) 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

!       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

!       The NAG name equivalent of dgetrf is f07adf
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

!           The NAG name equivalent of dgetrs is f07aef
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

!           The NAG name equivalent of dgetrs is f07aef
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
!       .. 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 ..
nout = 6
nin = 5
Write (nout,'(//1X,A/)') 'Example 2'

!       Skip heading in data file

!       Input N, the order of the matrix A, and NZ the number of nonzero
!       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.

!       Compute 1-norm of A
nrma = zero
Do i = 1, n
asum = zero
Do j = 1, nz
If (icn(j)==i) Then
asum = asum + abs(a(j))
End If
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 an 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
```