Program e04rhfe
! E04RHF Example Program Text
! Mark 28.6 Release. NAG Copyright 2022.
! 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, Intrinsic :: iso_c_binding, Only: c_null_ptr, &
c_ptr
Use nag_library, Only: e04raf, e04rff, e04rhf, e04rnf, e04rzf, e04svf, &
e04zmf, f08kbf, nag_wp, x04cbf
! .. 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 ..
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==0 .Or. ifail/=0) Then
! Final exit, solver finished.
Write (nout,99997) int(stats(1)), rinfo(1), &
sum(rinfo(2:4))/3.0_nag_wp
Flush (nout)
Exit loop
Else
! 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
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