Program e04rffe

!     E04RFF Example Program Text

!     Compute the nearest correlation matrix in Frobenius norm
!     using nonlinear semidefinite programming. By default,
!     preserve the nonzero structure of the input matrix
!     (preserve_structure = .True.).

!     Mark 27.3 Release. NAG Copyright 2021.

!     .. Use Statements ..
      Use, Intrinsic                   :: iso_c_binding, Only: c_null_ptr,     &
                                          c_ptr
      Use nag_library, Only: e04raf, e04rff, e04rnf, e04rzf, e04svf, e04zmf,   &
                             nag_wp, x04caf
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: nin = 5, nout = 6
      Logical, Parameter               :: preserve_structure = .True.
!     .. Local Scalars ..
      Type (c_ptr)                     :: h
      Integer                          :: dima, i, idblk, idx, ifail, inform,  &
                                          j, n, nblk, nnzasum, nnzc, nnzh,     &
                                          nnzu, nnzua, nnzuc, nvar
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: a(:), g(:,:), hmat(:), x(:)
      Real (Kind=nag_wp)               :: rdummy(1), rinfo(32), stats(32)
      Integer, Allocatable             :: blksizea(:), icola(:), icolh(:),     &
                                          irowa(:), irowh(:), nnza(:)
      Integer                          :: idummy(1)
!     .. Executable Statements ..
      Continue

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

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

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

      Allocate (g(n,n))

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

!     Symmetrize G:  G = (G + G')/2
      Do j = 2, n
        Do i = 1, j - 1
          g(i,j) = (g(i,j)+g(j,i))/2.0_nag_wp
          g(j,i) = g(i,j)
        End Do
      End Do

!     Initialize handle.
      h = c_null_ptr

!     There are as many variables as nonzeros above the main diagonal in
!     the input matrix. The variables are corrections of these elements.
      nvar = 0
      Do j = 2, n
        Do i = 1, j - 1
          If (.Not. preserve_structure .Or. g(i,j)/=0.0_nag_wp) Then
            nvar = nvar + 1
          End If
        End Do
      End Do
      Allocate (x(nvar))

!     Initialize an empty problem handle with NVAR variables.
      ifail = 0
      Call e04raf(h,nvar,ifail)

!     Set up the objective - minimize Frobenius norm of the corrections.
!     Our variables are stored as a vector thus, just minimize
!     sum of squares of the corrections --> H is identity matrix, c = 0.
      nnzc = 0
      nnzh = nvar
      Allocate (irowh(nnzh),icolh(nnzh),hmat(nnzh))
      Do i = 1, nvar
        irowh(i) = i
        icolh(i) = i
        hmat(i) = 1.0_nag_wp
      End Do

!     Add the quadratic objective to the handle.
      ifail = 0
      Call e04rff(h,nnzc,idummy,rdummy,nnzh,irowh,icolh,hmat,ifail)

!     Construct linear matrix inequality to request that
!     matrix G with corrections X is positive semidefinite.
!     (Don't forget the sign at A_0!)

!     How many nonzeros do we need? Full triangle for A_0 and
!     one nonzero element for each A_i.
      nnzasum = n*(n+1)/2 + nvar

      Allocate (nnza(nvar+1),irowa(nnzasum),icola(nnzasum),a(nnzasum))
      nnza(1) = n*(n+1)/2
      nnza(2:nvar+1) = 1

!     Copy G to A_0, only upper triangle with different sign (because -A_0)
!     and set the diagonal to 1.0 as that's what we want independently
!     of what was in G.
      idx = 1
      Do j = 1, n
        Do i = 1, j - 1
          irowa(idx) = i
          icola(idx) = j
          a(idx) = -g(i,j)
          idx = idx + 1
        End Do
!       Unit diagonal.
        irowa(idx) = j
        icola(idx) = j
        a(idx) = -1.0_nag_wp
        idx = idx + 1
      End Do
!     A_i has just one nonzero - it binds x_i with its position as
!     a correction.
      Do j = 2, n
        Do i = 1, j - 1
          If (.Not. preserve_structure .Or. g(i,j)/=0.0_nag_wp) Then
            irowa(idx) = i
            icola(idx) = j
            a(idx) = 1.0_nag_wp
            idx = idx + 1
          End If
        End Do
      End Do

!     Just one matrix inequality of the dimension of the original matrix.
      nblk = 1
      Allocate (blksizea(nblk))
      dima = n
      blksizea(:) = (/dima/)

!     Add the constraint to the problem formulation.
      idblk = 0
      ifail = 0
      Call e04rnf(h,nvar,dima,nnza,nnzasum,irowa,icola,a,nblk,blksizea,idblk,  &
        ifail)

!     Set optional arguments of the solver.
      ifail = 0
      Call e04zmf(h,'Print Options = No',ifail)
      ifail = 0
      Call e04zmf(h,'Initial X = Automatic',ifail)

!     Pass the handle to the solver, we are not interested in
!     Lagrangian multipliers.
      nnzu = 0
      nnzuc = 0
      nnzua = 0
      ifail = 0
      Call e04svf(h,nvar,x,nnzu,rdummy,nnzuc,rdummy,nnzua,rdummy,rinfo,stats,  &
        inform,ifail)

!     Destroy the handle.
      ifail = 0
      Call e04rzf(h,ifail)

!     Form the new nearest correlation matrix as the sum
!     of G and the correction X.
      idx = 1
      Do j = 1, n
        Do i = 1, j - 1
          If (.Not. preserve_structure .Or. g(i,j)/=0.0_nag_wp) Then
            g(i,j) = g(i,j) + x(idx)
            idx = idx + 1
          End If
        End Do
        g(j,j) = 1.0_nag_wp
      End Do

!     Print the matrix.
      ifail = 0
      Call x04caf('Upper','N',n,n,g,n,'Nearest Correlation Matrix',ifail)

    End Program e04rffe