Program f04ydfe
! F04YDF Example Program Text
! Mark 27.1 Release. NAG Copyright 2020.
! .. 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,'(///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
! 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 ..
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 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.
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) 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