NAG Library Manual, Mark 27.3
Interfaces:  FL   CL   CPP   AD 

NAG FL Interface Introduction
Example description
    Program f04ydfe

!     F04YDF Example Program Text

!     Mark 27.3 Release. NAG Copyright 2021.

!     .. 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