NAG Library Manual, Mark 30
Interfaces:  FL   CL   CPP   AD 

NAG FL Interface Introduction
Example description
    Program e04rpfe

!     E04RPF Example Program Text

!     Read a 'generic' LMI/BMI SDP problem from the input file,
!     formulate the problem via a handle and pass it to the solver.

!     Mark 30.0 Release. NAG Copyright 2024.

!     .. Use Statements ..
      Use, Intrinsic                   :: iso_c_binding, Only: c_ptr
      Use nag_library, Only: e04raf, e04ref, e04rff, e04rjf, e04rnf, e04rpf,   &
                             e04ryf, e04rzf, e04svf, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
      Type (c_ptr)                     :: handle
      Integer                          :: blkidx, dimaq, idblk, idlc, idx,     &
                                          idxend, ifail, inform, midx, nblk,   &
                                          nclin, nnzasum, nnzb, nnzc, nnzh,    &
                                          nnzqsum, nnzu, nnzua, nnzuc, nq,     &
                                          nvar
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: a(:), b(:), bl(:), bu(:), cvec(:),   &
                                          h(:), q(:), x(:)
      Real (Kind=nag_wp)               :: rdummy(1), rinfo(32), stats(32)
      Integer, Allocatable             :: icola(:), icolb(:), icolh(:),        &
                                          icolq(:), idxc(:), irowa(:),         &
                                          irowb(:), irowh(:), irowq(:),        &
                                          nnza(:), nnzq(:), qi(:), qj(:)
      Integer                          :: idummy(1)
!     .. Executable Statements ..

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

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

!     Read the problem size.
      Read (nin,*) nvar
      Read (nin,*) nnzh
      Read (nin,*) nclin, nnzb
      Read (nin,*) nblk

!     Initialize handle to an empty problem.
      ifail = 0
      Call e04raf(handle,nvar,ifail)

!     Read the linear part of the objective function.
      Allocate (cvec(nvar))
      Read (nin,*) cvec(1:nvar)

      If (nnzh==0) Then
!       Add the linear objective function to the problem formulation.
        ifail = 0
        Call e04ref(handle,nvar,cvec,ifail)
        Deallocate (cvec)

      Else
!       The linear part of the objective was read in as dense, E04RFF needs
!       the sparse format.
        nnzc = nvar
        Allocate (idxc(nnzc))
        Do idx = 1, nvar
          idxc(idx) = idx
        End Do

!       Read nonzeros for H (quadratic part of the objective) if present.
        Allocate (irowh(nnzh),icolh(nnzh),h(nnzh))
        Do idx = 1, nnzh
          Read (nin,*) h(idx), irowh(idx), icolh(idx)
        End Do

!       Add the quadratic objective function to the problem formulation.
        ifail = 0
        Call e04rff(handle,nnzc,idxc,cvec,nnzh,irowh,icolh,h,ifail)
        Deallocate (idxc,cvec,irowh,icolh,h)
      End If

!     Read a block of linear constraints and its bounds if present.
      If (nclin>0 .And. nnzb>0) Then
        Allocate (irowb(nnzb),icolb(nnzb),b(nnzb),bl(nclin),bu(nclin))
        Do idx = 1, nnzb
          Read (nin,*) b(idx), irowb(idx), icolb(idx)
        End Do
        Read (nin,*) bl(1:nclin)
        Read (nin,*) bu(1:nclin)

!       Add the block of linear constraints.
        idlc = 0
        ifail = 0
        Call e04rjf(handle,nclin,bl,bu,nnzb,irowb,icolb,b,idlc,ifail)
        Deallocate (irowb,icolb,b,bl,bu)
      End If

!     Read all matrix inequalities.
      Do blkidx = 1, nblk
        Read (nin,*) dimaq
        Read (nin,*) nnzasum, nnzqsum
        idblk = 0

        If (nnzasum>0) Then
!         Read a linear matrix inequality composed of (NVAR+1) matrices.
          Allocate (nnza(nvar+1),irowa(nnzasum),icola(nnzasum),a(nnzasum))
          idx = 1
          Do midx = 1, nvar + 1
!           Read matrix A_{midx-1}.
            Read (nin,*) nnza(midx)
            idxend = idx + nnza(midx) - 1
            Do idx = idx, idxend
              Read (nin,*) a(idx), irowa(idx), icola(idx)
            End Do
          End Do

!         Add the linear matrix inequality to the problem formulation.
          idblk = 0
          ifail = 0
          Call e04rnf(handle,nvar,dimaq,nnza,nnzasum,irowa,icola,a,1,idummy,   &
            idblk,ifail)
          Deallocate (nnza,irowa,icola,a)
        End If

        If (nnzqsum>0) Then
!         Read bilinear part of the matrix inequality composed of NQ matrices.
          Read (nin,*) nq
          Allocate (qi(nq),qj(nq),nnzq(nq),irowq(nnzqsum),icolq(nnzqsum),      &
            q(nnzqsum))
          idx = 1
          Do midx = 1, nq
!           Read matrix Q_ij where i=QI(midx), j=QJ(midx).
            Read (nin,*) qi(midx), qj(midx)
            Read (nin,*) nnzq(midx)
            idxend = idx + nnzq(midx) - 1
            Do idx = idx, idxend
              Read (nin,*) q(idx), irowq(idx), icolq(idx)
            End Do
          End Do

!         Expand the existing linear matrix inequality with the bilinear terms
!         or (if linear part was not present) create a new matrix inequality.
          ifail = 0
          Call e04rpf(handle,nq,qi,qj,dimaq,nnzq,nnzqsum,irowq,icolq,q,idblk,  &
            ifail)
          Deallocate (qi,qj,nnzq,irowq,icolq,q)
        End If

      End Do

!     Problem was successfully decoded.
      Write (nout,*) 'SDP problem was read, passing it to the solver.'
      Write (nout,*)
      Flush (nout)

!     Print overview of the handle.
      ifail = 0
      Call e04ryf(handle,nout,'Overview,Matrix Constraints',ifail)

!     Allocate memory for the solver.
      Allocate (x(nvar))

!     Call the solver, ignore Lagrangian multipliers.
      nnzu = 0
      nnzuc = 0
      nnzua = 0
      inform = 0
      x(:) = 0.0_nag_wp

      ifail = 0
      Call e04svf(handle,nvar,x,nnzu,rdummy,nnzuc,rdummy,nnzua,rdummy,rinfo,   &
        stats,inform,ifail)

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

    End Program e04rpfe