Program e04rnfe

!     E04RNF Example Program Text
!     Mark 26.1 Release. NAG Copyright 2016.

!     Compute E-optimal experiment design via semidefinite programming,
!     this can be done as follows
!       max {lambda_min(A) | A = sum x_i*v_i*v_i^T, x_i>=0, sum x_i = 1}
!     where v_i are given vectors.

!     Use nag_library

!     .. Use Statements ..
      Use, Intrinsic                   :: iso_c_binding, Only: c_null_ptr,     &
                                          c_ptr
      Use nag_library, Only: e04raf, e04rff, e04rhf, e04rjf, e04rnf, e04rzf,   &
                             e04svf, e04zmf, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Real (Kind=nag_wp), Parameter    :: big = 1E+20_nag_wp
      Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
      Type (c_ptr)                     :: h
      Real (Kind=nag_wp)               :: tol
      Integer                          :: dima, i, idblk, idlc, idx, ifail,    &
                                          inform, j, k, m, nblk, nnzasum,      &
                                          nnzb, nnzc, nnzu, nnzua, nnzuc,      &
                                          nvar, p
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: a(:), b(:), bl(:), bu(:), c(:),      &
                                          v(:,:), x(:)
      Real (Kind=nag_wp)               :: rdummy(1), rinfo(32), stats(32)
      Integer, Allocatable             :: blksizea(:), icola(:), icolb(:),     &
                                          idxc(:), irowa(:), irowb(:), nnza(:)
      Integer                          :: idummy(1)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: repeat
!     .. Executable Statements ..
      Continue

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

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

!     Read in the number of vectors and their size.
      Read (nin,*) m
      Read (nin,*) p

      Allocate (v(p,m))

!     Read in the vectors v_j.
      Do j = 1, m
        Read (nin,*)(v(i,j),i=1,p)
      End Do

!     Initialize handle.
      h = c_null_ptr

!     Variables of the problem will be x_1, ..., x_m (weights of the vectors)
!     and t (artificial variable for minimum eigenvalue).
      nvar = m + 1

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

!     Add the objective function to the handle: max t.
      nnzc = 1
      Allocate (idxc(nnzc),c(nnzc))
      idxc(:) = (/m+1/)
      c(:) = (/1._nag_wp/)

      ifail = 0
      Call e04rff(h,nnzc,idxc,c,0,idummy,idummy,rdummy,ifail)

      Allocate (bl(nvar),bu(nvar))
      bl(1:m) = 0.0_nag_wp
      bl(m+1) = -big
      bu(1:m+1) = big

!     Add simple bounds on variables, x_i>=0.
      ifail = 0
      Call e04rhf(h,nvar,bl,bu,ifail)

      nnzb = m
      Allocate (irowb(nnzb),icolb(nnzb),b(nnzb))
      irowb(:) = 1
      icolb(:) = (/(j,j=1,m)/)
      b(:) = 1.0_nag_wp

!     Add the linear constraint: sum x_i = 1.
      idlc = 0
      ifail = 0
      Call e04rjf(h,1,(/1.0_nag_wp/),(/1.0_nag_wp/),nnzb,irowb,icolb,b,idlc,   &
        ifail)

!     Generate matrix constraint as:
!       sum_{i=1}^m x_i*v_i*v_i^T - t*I >=0

      nblk = 1
      dima = p

!     Total number of nonzeros
      nnzasum = p + m*(p+1)*p/2

      Allocate (nnza(nvar+1),irowa(nnzasum),icola(nnzasum),a(nnzasum),x(nvar))
!     A_0 is empty
      nnza(1) = 0
!     A_1, A_2, ..., A_m are v_i*v_i^T
      nnza(2:m+1) = (p+1)*p/2
      idx = 0
      Do k = 1, m
        Do i = 1, p
          Do j = i, p
            idx = idx + 1
            irowa(idx) = i
            icola(idx) = j
            a(idx) = v(i,k)*v(j,k)
          End Do
        End Do
      End Do
!     A_{m+1} is the -identity
      nnza(m+2) = p
      Do i = 1, p
        idx = idx + 1
        irowa(idx) = i
        icola(idx) = 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,'Task = Maximize',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)

!     Print results
      Write (nout,*)
      tol = 0.00001_nag_wp
      Write (nout,*) '  Weight        Row of design matrix'
      Write (nout,*) repeat('-',13+p*8)
      Do j = 1, m
        If (x(j)>tol) Then
          Write (*,99999) x(j), v(1:p,j)
        End If
      End Do
      Write (nout,99998) 'only those rows with a weight > ', tol, ' are shown'

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

99999 Format (1X,F7.2,5X,10(1X,F7.2))
99998 Format (1X,A,E8.1,A)
    End Program e04rnfe