Program e04rhfe

!     E04RHF Example Program Text
!     Mark 26 Release. NAG Copyright 2016.

!     Matrix completion problem (rank minimization) solved approximately
!     by SDP via nuclear norm minimization formulated as:
!        min    trace(X1) + trace(X2)
!        s.t.   [ X1, Y; Y', X2 ] >=0
!               0 <= Y_ij <= 1

!     .. Use Statements ..
      Use nag_library, Only: e04raf, e04rff, e04rhf, e04rnf, e04rzf, e04svf,   &
                             e04zmf, f08kbf, nag_wp, x04cbf
      Use, Intrinsic                   :: iso_c_binding, Only: c_null_ptr,     &
                                          c_ptr
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Real (Kind=nag_wp), Parameter    :: stol = 1E-5_nag_wp
      Real (Kind=nag_wp), Parameter    :: time_limit = 120.0_nag_wp
      Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
      Type (c_ptr)                     :: h
      Integer                          :: dima, i, idblk, idx, idxobj, idxx,   &
                                          ifail, info, inform, j, lwork, m, n, &
                                          nblk, nnz, nnzasum, nnzc, nnzh,      &
                                          nnzu, nnzua, nnzuc, nvar, rank
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: a(:), bl(:), bu(:), c(:), s(:),      &
                                          work(:), x(:), y(:,:)
      Real (Kind=nag_wp)               :: rdummy(1), rinfo(32), stats(32)
      Integer, Allocatable             :: blksizea(:), icola(:), idxc(:),      &
                                          irowa(:), nnza(:)
      Integer                          :: idummy(1)
      Character (1)                    :: cdummy(1)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: int, max, min, sum
!     .. Executable Statements ..
      Continue

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

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

!     Read in the problem size and allocate space for the input data.
      Read (nin,*) m, n
      Allocate (y(m,n))

!     Read in the matrix Y.
      Read (nin,*)(y(i,1:n),i=1,m)

!     Count the number of specified elements (i.e., nonnegative)
      nnz = 0
      Do i = 1, m
        Do j = 1, n
          If (y(i,j)>=0.0_nag_wp) Then
            nnz = nnz + 1
          End If
        End Do
      End Do

!     Initialize handle.
      h = c_null_ptr

!     There are as many variables as missing entries in the Y matrix
!     plus two full symmetric matrices m x m and n x n.
      nvar = m*(m+1)/2 + n*(n+1)/2 + m*n - nnz
      Allocate (x(nvar),bl(nvar),bu(nvar))

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

!     Create bounds for the missing entries in Y matrix to be between 0 and 1
      bl(:) = -1E+20_nag_wp
      bu(:) = 1E+20_nag_wp
      bl(m*(m+1)/2+n*(n+1)/2+1:nvar) = 0.0_nag_wp
      bu(m*(m+1)/2+n*(n+1)/2+1:nvar) = 1.0_nag_wp
      ifail = 0
      Call e04rhf(h,nvar,bl,bu,ifail)

!     Allocate space for the objective - minimize trace of the matrix
!     constraint. There is no quadratic part in the objective.
      nnzc = m + n
      nnzh = 0
      Allocate (idxc(nnzc),c(nnzc))

!     Construct linear matrix inequality to request that
!     [ X1, Y; Y', X2] is positive semidefinite.

!     How many nonzeros do we need? As many as number of variables
!     and the number of specified elements together.
      nnzasum = m*(m+1)/2 + n*(n+1)/2 + m*n

      Allocate (nnza(nvar+1),irowa(nnzasum),icola(nnzasum),a(nnzasum))
      nnza(1) = nnz
      nnza(2:nvar+1) = 1

!     Copy Y to the upper block of A_0 with the different sign
!     (because of the sign at A_0!)
!     (upper triangle)
      idx = 1
      Do i = 1, m
        Do j = 1, n
          If (y(i,j)>=0.0_nag_wp) Then
            irowa(idx) = i
            icola(idx) = m + j
            a(idx) = -y(i,j)
            idx = idx + 1
          End If
        End Do
      End Do
!     One matrix for each variable, A_i has just one nonzero - it binds
!     x_i with its position in the matrix constraint. Set also the objective.
!     1,1 - block, X1 matrix (mxm)
      idxobj = 1
      idxx = 1
      Do i = 1, m
!       the next element is diagonal ==> part of the objective as a trace()
        idxc(idxobj) = idxx
        c(idxobj) = 1.0_nag_wp
        idxobj = idxobj + 1
        Do j = i, m
          irowa(idx) = i
          icola(idx) = j
          a(idx) = 1.0_nag_wp
          idx = idx + 1
          idxx = idxx + 1
        End Do
      End Do
!     2,2 - block, X2 matrix (nxn)
      Do i = 1, n
!       the next element is diagonal ==> part of the objective as a trace()
        idxc(idxobj) = idxx
        c(idxobj) = 1.0_nag_wp
        idxobj = idxobj + 1
        Do j = i, n
          irowa(idx) = m + i
          icola(idx) = m + j
          a(idx) = 1.0_nag_wp
          idx = idx + 1
          idxx = idxx + 1
        End Do
      End Do
!     1,2 - block, missing element in Y we are after
      Do i = 1, m
        Do j = 1, n
          If (y(i,j)<0.0_nag_wp) Then
            irowa(idx) = i
            icola(idx) = m + j
            a(idx) = 1.0_nag_wp
            idx = idx + 1
          End If
        End Do
      End Do

!     Add the sparse linear objective to the handle.
      ifail = 0
      Call e04rff(h,nnzc,idxc,c,nnzh,idummy,idummy,rdummy,ifail)

!     Just one matrix inequality of the dimension of the extended matrix.
      nblk = 1
      Allocate (blksizea(nblk))
      dima = m + 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.
!     Completely turn off printing, allow timing and
!     turn on the monitor mode to stop every iteration.
      ifail = 0
      Call e04zmf(h,'Print File = -1',ifail)
      ifail = 0
      Call e04zmf(h,'Stats Time = Yes',ifail)
      ifail = 0
      Call e04zmf(h,'Monitor Frequency = 1',ifail)
      ifail = 0
      Call e04zmf(h,'Initial X = Automatic',ifail)
      ifail = 0
      Call e04zmf(h,'Dimacs = Check',ifail)

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

        If (inform==1) Then
!         Monitor stop
          Write (nout,99998) int(stats(1)), rinfo(1),                          &
            sum(rinfo(2:4))/3.0_nag_wp
          Flush (nout)

!         Check time limit and possibly stop the solver.
          If (stats(8)>time_limit) Then
            inform = 0
          End If
        Else
!         Final exit, solver finished.
          Write (nout,99997) int(stats(1)), rinfo(1),                          &
            sum(rinfo(2:4))/3.0_nag_wp
          Flush (nout)
          Exit loop
        End If

      End Do loop

      If (ifail==0 .Or. ifail==50) Then
!       Successful run, fill the missing elements in the matrix Y.
        idx = m*(m+1)/2 + n*(n+1)/2 + 1
        Do i = 1, m
          Do j = 1, n
            If (y(i,j)<0.0_nag_wp) Then
              y(i,j) = x(idx)
              idx = idx + 1
            End If
          End Do
        End Do

!       Print the matrix.
        ifail = 0
        Call x04cbf('General','N',m,n,y,m,'F7.1','Completed Matrix','Integer', &
          cdummy,'Integer',cdummy,80,0,ifail)

!       Compute rank of the matrix via SVD, use the fact that the order
!       of the singular values is descending.
        lwork = 20*max(m,n)
        Allocate (s(min(m,n)),work(lwork))
        Call f08kbf('No','No',m,n,y,m,s,rdummy,1,rdummy,1,work,lwork,info)
        If (info==0) Then
lp_rank:  Do rank = 1, min(m,n)
            If (s(rank)<=stol) Then
              Exit lp_rank
            End If
          End Do lp_rank
          Write (nout,99999) 'Rank is', rank - 1
99999     Format (1X,A,I20)
        End If
      Else If (ifail==20) Then
        Write (nout,*) 'The given time limit was reached, run aborted.'
      End If

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

99998 Format (1X,'Monitor at iteration  ',I2,': objective ',F7.2,              &
        ', avg.error ',Es9.2e2)
99997 Format (1X,'Finished at iteration ',I2,': objective ',F7.2,              &
        ', avg.error ',Es9.2e2)
    End Program e04rhfe