! E04PTF Example Program Text
! Mark 27.0 Release. NAG Copyright 2019.
Module e04ptfe_mod
! .. Use Statements ..
Use nag_library, Only: nag_wp
! .. Implicit None Statement ..
Implicit None
! .. Accessibility Statements ..
Private
Public :: monit
Contains
Subroutine monit(handle,rinfo,stats,iuser,ruser,cpuser,inform)
! Monitoring function can be used to monitor the progress
! or, for example, to implement bespoke stopping criteria
! .. Use Statements ..
Use iso_c_binding, Only: c_ptr
! .. Scalar Arguments ..
Type (c_ptr), Intent (In) :: cpuser, handle
Integer, Intent (Inout) :: inform
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (In) :: rinfo(100), stats(100)
Real (Kind=nag_wp), Intent (Inout) :: ruser(*)
Integer, Intent (Inout) :: iuser(*)
! .. Local Scalars ..
Real (Kind=nag_wp) :: tol
Integer :: nout, tol_reached
! .. Executable Statements ..
nout = iuser(1)
tol_reached = iuser(2)
tol = ruser(1)
! If x is close to the solution, print a message
If (rinfo(15)<tol .And. rinfo(16)<tol .And. rinfo(17)<tol .And. &
rinfo(18)<tol) Then
If (tol_reached==0) Then
Write (nout,*)
Write (nout,99999) &
'monit() reports good approximate solution (tol =', tol, ')'
iuser(2) = 1
End If
End If
Return
99999 Format (5X,A,Es9.2,A)
End Subroutine monit
End Module e04ptfe_mod
Program e04ptfe
! .. Use Statements ..
Use e04ptfe_mod, Only: monit
Use iso_c_binding, Only: c_null_ptr, c_ptr
Use nag_library, Only: dsyevd, e04ptf, e04raf, e04rbf, e04ref, e04rhf, &
e04rjf, e04rzf, e04zmf, nag_wp
! .. Implicit None Statement ..
Implicit None
! .. Parameters ..
Integer, Parameter :: nin = 5, nout = 6, nqc = 1
! .. Local Scalars ..
Type (c_ptr) :: cpuser, handle
Real (Kind=nag_wp) :: r1
Integer :: i, idgroup, idlc, idxa, ifail, j, &
liwork, lwork, m, n, na, nnza, &
nnzp0, nnzp1, nnzu, nnzuc, nu, nv, &
nvarc1, nvarc2, rptr, x_idx
Logical :: verbose_output
Character (8) :: ctype1, ctype2
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: a(:), bla(:), bua(:), c(:), f0(:,:), &
f1(:,:), lambda0(:), lambda1(:), &
p0(:), p1(:), q0(:), q1(:), u(:), &
uc(:), work(:), x(:), xl(:), xu(:)
Real (Kind=nag_wp) :: rinfo(100), ruser(1), stats(100)
Integer, Allocatable :: icola(:), icolp0(:), icolp1(:), &
idxc1(:), idxc2(:), irowa(:), &
irowp0(:), irowp1(:), iwork(:)
Integer :: iuser(2)
! .. Intrinsic Procedures ..
Intrinsic :: max, sqrt
! .. Executable Statements ..
Write (nout,*) 'E04PTF Example Program Results'
! Skip Header in data file
Read (nin,*)
! Read dimensions of the problem
Read (nin,*) n, nnzp0, nnzp1
! Initialize size of linear constraints in SOCP
m = nqc
na = n + nqc + 1
nnza = nqc + n
! Initialize size of cone constraints
nvarc1 = 2
nvarc2 = 2
! Allocate memory to read data
lwork = max(2*n*n+6*n+1,120+9*n)
liwork = 5*n + 3
Allocate (irowp0(nnzp0),icolp0(nnzp0),p0(nnzp0),irowp1(nnzp1), &
icolp1(nnzp1),p1(nnzp1),q0(n),q1(n),f0(n,n),f1(n,n),lambda0(n), &
lambda1(n),work(lwork),iwork(liwork))
! Read problem data
Read (nin,*) irowp0(1:nnzp0)
Read (nin,*) icolp0(1:nnzp0)
Read (nin,*) p0(1:nnzp0)
Read (nin,*) irowp1(1:nnzp1)
Read (nin,*) icolp1(1:nnzp1)
Read (nin,*) p1(1:nnzp1)
Read (nin,*) q0(1:n)
Read (nin,*) q1(1:n)
Read (nin,*) r1
! Store full P0 and P1 in F0 and F1
f0(:,:) = 0.0_nag_wp
Do i = 1, nnzp0
f0(irowp0(i),icolp0(i)) = p0(i)
End Do
f1(:,:) = 0.0_nag_wp
Do i = 1, nnzp1
f1(irowp1(i),icolp1(i)) = p1(i)
End Do
! Factorize P0 and P1 via eigenvalue decomposition
Call dsyevd(job='V',uplo='L',n=n,a=f0,lda=n,w=lambda0,work=work, &
lwork=lwork,iwork=iwork,liwork=liwork,info=ifail)
Call dsyevd(job='V',uplo='L',n=n,a=f1,lda=n,w=lambda1,work=work, &
lwork=lwork,iwork=iwork,liwork=liwork,info=ifail)
! Fomulate F0 and F1 in P0 = F0'*F0, P1 = F1'*F1
nu = 0
nv = 0
Do i = 1, n
If (lambda0(i)>0) Then
f0(:,i) = f0(:,i)*sqrt(lambda0(i))
m = m + 1
nu = nu + 1
nnza = nnza + n
End If
If (lambda1(i)>0) Then
f1(:,i) = f1(:,i)*sqrt(lambda1(i))
m = m + 1
nv = nv + 1
nnza = nnza + n
End If
End Do
nnza = nnza + nu + nv
na = na + nu + nv
nvarc1 = nvarc1 + nu
nvarc2 = nvarc2 + nv
! Add two fixed variable for two rotated quadratic cones
na = na + 2
m = m + 2
nnza = nnza + 2
! Compute size of multipliers
nnzu = 2*na + 2*m
nnzuc = nvarc1 + nvarc2
! Allocate memory to build SOCP
Allocate (icola(nnza),irowa(nnza),a(nnza),bla(m),bua(m),xl(na),xu(na), &
c(na),x(na),u(nnzu),uc(nnzuc),idxc1(nvarc1),idxc2(nvarc2))
! Build objective function parameter c
! [x, t1, u, v, y1, y2, t0]
c(:) = 0.0_nag_wp
c(1:n) = q0(1:n)
c(na) = 1.0_nag_wp
! Build linear constraints parameter A
idxa = 0
rptr = 0
! q1 in First row
rptr = rptr + 1
irowa(1:n) = rptr
icola(1:n) = (/(j,j=1,n)/)
a(1:n) = q1(1:n)
idxa = idxa + n
! F0 in F0*x-u=0 row
Do i = 1, n
If (lambda0(i)>0) Then
rptr = rptr + 1
irowa(idxa+1:idxa+n) = rptr
icola(idxa+1:idxa+n) = (/(j,j=1,n)/)
a(idxa+1:idxa+n) = f0(1:n,i)
idxa = idxa + n
End If
End Do
! F1 in F1*x-v=0 row
Do i = 1, n
If (lambda1(i)>0) Then
rptr = rptr + 1
irowa(idxa+1:idxa+n) = rptr
icola(idxa+1:idxa+n) = (/(j,j=1,n)/)
a(idxa+1:idxa+n) = f1(1:n,i)
idxa = idxa + n
End If
End Do
! Rest of A, a diagonal matrix
irowa(idxa+1:idxa+m) = (/(j,j=1,m)/)
icola(idxa+1:idxa+m) = (/(j+n,j=1,m)/)
a(idxa+1:idxa+m) = 1.0_nag_wp
a(idxa+2:idxa+1+nu+nv) = -1.0_nag_wp
! RHS in linear constraints
bla(:) = 0.0_nag_wp
bua(:) = 0.0_nag_wp
bla(1) = -r1
bua(1) = -r1
bla(m-nqc:m) = 1.0_nag_wp
bua(m-nqc:m) = 1.0_nag_wp
! Box constraints, all variables are free
xl(:) = -1E+20_nag_wp
xu(:) = 1E+20_nag_wp
! Cone constraints
! First cone
idxc1(1) = na
idxc1(2) = n + 1 + nu + nv + 1
idxc1(3:nvarc1) = (/(j,j=n+2,n+1+nu)/)
ctype1 = 'RQUAD'
! Second cone
idxc2(1) = n + 1
idxc2(2) = n + 1 + nu + nv + 2
idxc2(3:nvarc2) = (/(j,j=n+nu+2,n+1+nu+nv)/)
ctype2 = 'RQUAD'
! Create the problem handle
ifail = 0
Call e04raf(handle,na,ifail)
! Set objective function
ifail = 0
Call e04ref(handle,na,c,ifail)
! Set box constraints
ifail = 0
Call e04rhf(handle,na,xl,xu,ifail)
! Set linear constraints
ifail = 0
idlc = 0
Call e04rjf(handle,m,bla,bua,nnza,irowa,icola,a,idlc,ifail)
! Set first cone constraint
ifail = 0
idgroup = 0
Call e04rbf(handle,ctype1,nvarc1,idxc1,idgroup,ifail)
! Set second cone constraint
ifail = 0
idgroup = 0
Call e04rbf(handle,ctype2,nvarc2,idxc2,idgroup,ifail)
! Turn on monitoring
ifail = 0
Call e04zmf(handle,'SOCP Monitor Frequency = 1',ifail)
! Set this to .True. to cause e04ptf to produce intermediate
! progress output
verbose_output = .False.
If (verbose_output) Then
! Require printing of primal and dual solutions at the end of the solve
ifail = 0
Call e04zmf(handle,'Print Solution = YES',ifail)
Else
! Turn off printing of intermediate progress output
ifail = 0
Call e04zmf(handle,'Print Level = 1',ifail)
End If
! Call SOCP interior point solver
cpuser = c_null_ptr
iuser(1) = nout
iuser(2) = 0
ruser(1) = 1.0E-07_nag_wp
ifail = -1
Call e04ptf(handle,na,x,nnzu,u,nnzuc,uc,rinfo,stats,monit,iuser,ruser, &
cpuser,ifail)
! Print solution if optimal or suboptimal solution found
If (ifail==0 .Or. ifail==50) Then
Write (nout,99999) 'Optimal X:'
Write (nout,99997) 'x_idx', ' Value '
Do x_idx = 1, n
Write (nout,99998) x_idx, x(x_idx)
End Do
End If
! Free the handle memory
ifail = 0
Call e04rzf(handle,ifail)
99999 Format (1X,A)
99998 Format (2X,I5,3X,Es11.3e2)
99997 Format (2X,A5,3X,A12)
End Program e04ptfe