Program f01mdfe

!     F01MDFE Example Program Text

!     Mark 27.1 Release. NAG Copyright 2020.

!     .. 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
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
Else If (uplo=='L') Then
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