Program e04svfe

!     E04SVF Example Program Text
!     Mark 27.1 Release. NAG Copyright 2020.

!     Compute Lovasz theta number of the given graph G on the input
!     via semidefinite programming as
!       min {lambda_max(H) | H is nv x nv symmetric matrix where
!                            h_ij=1 if ij is not an edge or if i==j}

!     .. Use Statements ..
      Use, Intrinsic                   :: iso_c_binding, Only: c_null_ptr,     &
                                          c_ptr
      Use nag_library, Only: e04raf, e04rff, e04rnf, e04rzf, e04svf, e04zmf,   &
                             e04znf, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
      Type (c_ptr)                     :: h
      Real (Kind=nag_wp)               :: rvalue
      Integer                          :: dima, i, idblk, idx, ifail, inform,  &
                                          ivalue, j, maxe, nblk, ne, nnzasum,  &
                                          nnzu, nnzua, nnzuc, nv, nvar, optype
      Character (40)                   :: cvalue
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: a(:), x(:)
      Real (Kind=nag_wp)               :: rdummy(1), rinfo(32), stats(32)
      Integer, Allocatable             :: blksizea(:), icola(:), irowa(:),     &
                                          nnza(:), va(:), vb(:)
      Integer                          :: idummy(1)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: trim
!     .. Executable Statements ..
      Continue

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

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

!     Read in the number of vertices and edges of the graph.
      Read (nin,*) nv
      Read (nin,*) ne

      Allocate (va(ne),vb(ne))

!     Read in edges of the graph, accept only 1<=i<j<=nv.
      maxe = ne
      ne = 0
      Do idx = 1, maxe
        Read (nin,*) i, j
        If (i>=1 .And. i<j .And. j<=nv) Then
          ne = ne + 1
          va(ne) = i
          vb(ne) = j
        End If
      End Do

!     Initialize handle.
      h = c_null_ptr

!     Number of variables (same as edges in the graph plus one).
      nvar = ne + 1

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

!     Add the objective function to the handle.
      ifail = 0
      Call e04rff(h,1,(/1/),(/1.0_nag_wp/),0,idummy,idummy,rdummy,ifail)

!     Generate matrix constraint as:
!       sum_{ij is edge in G} x_ij*E_ij + t*I - J >=0
!     where J is the all-ones matrix.

!     Just one matrix inequality.
      nblk = 1
      dima = nv

!     Total number of nonzeros
      nnzasum = ne + nv + (nv+1)*nv/2

      Allocate (nnza(nvar+1),irowa(nnzasum),icola(nnzasum),a(nnzasum),x(nvar))
!     A_0 is all ones matrix
      nnza(1) = (nv+1)*nv/2
      idx = 0
      Do i = 1, nv
        Do j = i, nv
          idx = idx + 1
          irowa(idx) = i
          icola(idx) = j
          a(idx) = 1.0_nag_wp
        End Do
      End Do
!     A_1 is the identity
      nnza(2) = nv
      Do i = 1, nv
        idx = idx + 1
        irowa(idx) = i
        icola(idx) = i
        a(idx) = 1.0_nag_wp
      End Do
!     A_2, A_3, ..., A_{ne+1} match the E_ij matrices
      nnza(3:ne+2) = 1
      Do i = 1, ne
        idx = idx + 1
        irowa(idx) = va(i)
        icola(idx) = vb(i)
        a(idx) = 1.0_nag_wp
      End Do

!     Add the constraint to the problem formulation.
      Allocate (blksizea(nblk))
      blksizea(:) = (/dima/)

      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,'Initial X = Automatic',ifail)

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

      If (ifail==0 .Or. ifail==50) Then
!       Retrieve some of the settings
        ifail = 0
        Call e04znf(h,'Hessian Density',ivalue,rvalue,cvalue,optype,ifail)
        Write (nout,*) 'The solver chose to use ', trim(cvalue), ' hessian'
        ifail = 0
        Call e04znf(h,'Linesearch Mode',ivalue,rvalue,cvalue,optype,ifail)
        Write (nout,*) 'and ', trim(cvalue), ' as linesearch.'

        Write (nout,Fmt=99999) 'Lovasz theta number of the given graph is',    &
          rinfo(1)
      End If

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

99999 Format (1X,A,1X,F7.2)
    End Program e04svfe