Program g02aafe

!     G02AAF Example Program Text

!     Mark 28.6 Release. NAG Copyright 2022.

!     .. Use Statements ..
      Use nag_library, Only: f06rcf, g02aaf, nag_wp, x04caf
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: dist, errtol, nrmgrd
      Integer                          :: feval, i, ifail, iter, j, lda, ldg,  &
                                          ldx, maxit, maxits, n
      Character (1)                    :: uplo
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: a(:,:), g(:,:), work(:), x(:,:)
      Integer, Allocatable             :: iwork(:)
!     .. Executable Statements ..

      Write (nout,*) 'G02AAF Example Program Results'
      Write (nout,*)
      Flush (nout)

!     Skip heading in data file
      Read (nin,*)

!     Read in the problem size
      Read (nin,*) n

      lda = n
      ldg = n
      ldx = n
      Allocate (a(lda,n),g(ldg,n),iwork(n),work(n),x(ldx,n))

!     Read in the matrix G
      Read (nin,*)(g(i,1:n),i=1,n)

      uplo = 'L'
!     Copy G into A to compute a bound on the distance to the NCM
!     We assume G is symmetric and copy only the lower triangle
      Do j = 1, n
        Do i = j, n
          a(i,j) = g(i,j)
        End Do
      End Do

      Call compute_bound(dist,uplo,n,g,ldg,a,lda,work,iwork,ifail)
      Write (nout,*)
      Write (nout,99999) 'Upper bound on the distance to the NCM:', dist
      Write (nout,*)

!     A was overwritten by compute_bound. In order to compute the actual
!     distance to the nearest correlation matrix, store another copy of G in A
      Do j = 1, n
        Do i = j, n
          a(i,j) = g(i,j)
        End Do
      End Do

!     Use the defaults for ERRTOL, MAXITS and MAXIT
      errtol = 0.0E0_nag_wp
      maxits = 0
      maxit = 0

!     Calculate nearest correlation matrix
      ifail = 0
      Call g02aaf(g,ldg,n,errtol,maxits,maxit,x,ldx,iter,feval,nrmgrd,ifail)

!     Compute distance between the original input matrix and the nearest
!     correlation matrix
      Do j = 1, n
        Do i = j, n
          a(i,j) = a(i,j) - x(i,j)
        End Do
      End Do

      dist = f06rcf('F',uplo,n,a,lda,work)

!     Display results
      ifail = 0
      Call x04caf('General',' ',n,n,x,ldx,'Nearest Correlation Matrix',ifail)
      Write (nout,*)
      Write (nout,99999) 'Distance between G and X:', dist
      Write (nout,99998) 'Number of Newton steps taken:', iter
      Write (nout,99997) 'Number of function evaluations:', feval

99999 Format (1X,A,F15.3)
99998 Format (1X,A,I11)
99997 Format (1X,A,I9)

    Contains
      Subroutine compute_bound(bound,uplo,n,g,ldg,a,lda,offdiag,ipiv,ifail)
!       Compute a bound on the distance between G and the nearest
!       correlation matrix using the modified Cholesky factorization.

!       .. Use Statements ..
        Use nag_library, Only: f01mdf, f01mef
!       .. Implicit None Statement ..
        Implicit None
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: bound
        Integer, Intent (Out)          :: ifail
        Integer, Intent (In)           :: lda, ldg, n
        Character (1), Intent (In)     :: uplo
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Inout) :: a(lda,n)
        Real (Kind=nag_wp), Intent (In) :: g(ldg,n)
        Real (Kind=nag_wp), Intent (Out) :: offdiag(n)
        Integer, Intent (Out)          :: ipiv(n)
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: delta
        Integer                        :: i, j
!       .. Local Arrays ..
        Real (Kind=nag_wp)             :: work(1)
!       .. Intrinsic Procedures ..
        Intrinsic                      :: sqrt
!       .. Executable Statements ..

        delta = 1.0E-7_nag_wp
!       Compute the modified Cholesky factorization of G
        ifail = 0
        Call f01mdf(uplo,n,a,lda,offdiag,ipiv,delta,ifail)

!       Compute the perturbed matrix G + E and store in A
        ifail = 0
        Call f01mef(uplo,n,a,lda,offdiag,ipiv,ifail)

!       Overwrite A with D^(-1/2)(G + E)D^(-1/2) - G
!       where D = diag(G + E)
        Do j = 1, n
          Do i = j + 1, n
            a(i,j) = a(i,j)/sqrt(a(i,i)*a(j,j)) - g(i,j)
          End Do
          a(j,j) = 1.0E0_nag_wp - g(j,j)
        End Do

!       Bound is given by ||D^(-1/2)(G + E)D^(-1/2) - G||
        bound = f06rcf('F',uplo,n,a,lda,work)

      End Subroutine compute_bound
    End Program g02aafe