Program f01mdfe

!     F01MDFE Example Program Text

!     Mark 28.3 Release. NAG Copyright 2022.

!     .. Use Statements ..
      Use nag_library, Only: dsyev, f01mdf, f01mef, f06rcf, nag_wp, x02ajf,    &
                             x04caf
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: delta, res
      Integer                          :: i, ifail, j, lda, lwork, n
      Character (1)                    :: uplo
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: a(:,:), d(:,:), e(:,:), eigs(:),     &
                                          offdiag(:), orig_a(:,:), work(:)
      Integer, Allocatable             :: ipiv(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: sqrt
!     .. Executable Statements ..

      Write (nout,*) 'F01MDF Example Program Results'
!     Skip heading in data file
      Read (nin,*)
      Read (nin,*) n, uplo
      lda = n
      lwork = 66*n
      Allocate (a(lda,n),offdiag(n),ipiv(n),work(lwork),orig_a(n,n),d(n,n),    &
        eigs(n),e(n,n))

!     Read A from data file
      If (uplo=='U') Then
        Read (nin,*)(a(i,i:n),i=1,n)
      Else If (uplo=='L') Then
        Read (nin,*)(a(i,1:i),i=1,n)
      End If

!     Store the original matrix A for use later
      Do j = 1, n
        Do i = j, n
          If (uplo=='L') Then
            orig_a(i,j) = a(i,j)
            orig_a(j,i) = a(i,j)
          Else
            orig_a(i,j) = a(j,i)
            orig_a(j,i) = a(j,i)
          End If
        End Do
      End Do

!     Display A
      Write (nout,*)
      Flush (nout)
      Call x04caf('G','Nonunit',n,n,orig_a,lda,'Original Matrix A',ifail)
      Write (nout,*)

!     Calculate a delta
      res = f06rcf('I',uplo,n,orig_a,lda,work)
      delta = sqrt(x02ajf())*res
!     Print delta.
      Write (nout,99999) 'Delta:', delta
      Flush (nout)

!     ifail: behaviour on error exit
!             =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
      ifail = 0

!     Compute the modified Cholesky factorization of A
      Call f01mdf(uplo,n,a,lda,offdiag,ipiv,delta,ifail)

!     Construct the matrix D and display it
      d(1:n,1:n) = 0.0_nag_wp
      Do i = 1, n
        d(i,i) = a(i,i)
      End Do

      If (uplo=='L') Then
        Do i = 1, n - 1
          d(i+1,i) = offdiag(i)
          d(i,i+1) = offdiag(i)
        End Do
      Else
        Do i = 2, n
          d(i-1,i) = offdiag(i)
          d(i,i-1) = offdiag(i)
        End Do
      End If

      Write (nout,*)
      Flush (nout)
      Call x04caf('G','Nonunit',n,n,d,lda,'Block Diagonal Matrix D',ifail)

      ifail = 0
!     Compute the eigenvalues of D and print them
      Call dsyev('N',uplo,n,d,n,eigs,work,lwork,ifail)
      Write (nout,*)
      Flush (nout)
      Call x04caf('General',' ',1,n,eigs,1,'Eigenvalues of D',ifail)

!     Compute the perturbed matrix A + E
      Call f01mef(uplo,n,a,lda,offdiag,ipiv,ifail)

!     Construct the matrix E and display it
      If (uplo=='L') Then
        Do j = 1, n
          e(j,j) = orig_a(j,j) - a(j,j)
          Do i = j + 1, n
            e(i,j) = orig_a(i,j) - a(i,j)
            e(j,i) = e(i,j)
          End Do
        End Do
      Else
        Do j = 1, n
          e(j,j) = orig_a(j,j) - a(j,j)
          Do i = 1, j - 1
            e(i,j) = orig_a(i,j) - a(i,j)
            e(j,i) = e(i,j)
          End Do
        End Do
      End If

      Write (nout,*)
      Flush (nout)
      Call x04caf('G','Nonunit',n,n,e,lda,'Perturbation Matrix E',ifail)
      Write (nout,*)

!     Calculate the Frobenius norm of E and print it
      res = f06rcf('F',uplo,n,e,n,work)
      Write (nout,99998) 'Frobenius Norm of E:', res
      Flush (nout)

99999 Format (1X,A,E12.4)
99998 Format (1X,A,F8.4)

    End Program f01mdfe