subroutine rqfnc(n1,n2,p,a1,y,a2,r,rhs,d1,d2,u,beta,eps,wn1,wn2,wp,nit,info) integer n1,n2,p,info,nit(3) double precision a1(p,n1),a2(p,n2),y(n1),r(n2),rhs(p),d1(n1),d2(n2),u(n1) double precision wn1(n1,9),wn2(n2,6),wp(p,p+3) double precision one,beta,eps parameter(one = 1.0d0) call lpfnc(n1,n2,p,a1,y,a2,r,rhs,d1,d2,u,beta,eps,wn1(1,1),wn2(1,1),wn1(1,2), wp(1,1),wn1(1,3),wn2(1,2),wn1(1,4),wn1(1,5),wn2(1,3),wn1(1,6), wp(1,2),wn1(1,7),wn2(1,4),wn1(1,8),wn1(1,9),wn2(1,5),wn2(1,6), wp(1,3),wp(1,4),nit,info) return end # This is a revised form of the primal-dual log barrier form of the # interior point LP solver based on Lustig, Marsten and Shanno ORSA J Opt 1992. # It generalizes rqfnb and lpfnb to accomodate inequality constraints # It is a projected Newton primal-dual logarithmic barrier method which uses # the predictor-corrector approach of Mehrotra for the mu steps. # For the sake of brevity we will call it a Frisch-Newton algorithm. # Problem: # min c1'x1 + c2'x2 s.t. A1x1 + A2x2 = b, 0<=x1<=u, 0<=x2 # # Denote dx1,dx2,dy,dw,ds,dz1,dz2 as the steps for x1,x2,y,w,s,z1,z2 # subroutine lpfnc(n1,n2,p,a1,c1,a2,c2,b,d1,d2,u,beta,eps,x1,x2,s, y,z1,z2,w,dx1,dx2,ds,dy,dz1,dz2,dw,dr1,dr2,r2, rhs,ada,nit,info) integer n1,p,i,info,nit(3),maxit double precision a1(p,n1),a2(p,n2),c1(n1),c2(n2),b(p) double precision zero,one,mone,big,ddot,dmax1,dmin1,dxdz1,dxdz2,dsdw double precision deltap,deltad,beta,eps,mu,gap,g double precision x1(n1),x2(n2),u(n1),s(n1),y(p),z1(n1),z2(n2),w(n1) double precision d1(n1),d2(n2),rhs(p),ada(p,p) double precision dx1(n1),dx2(n2),ds(n1),dy(p),dz1(n1),dz2(n2),dw(n1) double precision dr1(n1),dr2(n2),r2(n2) parameter(zero = 0.0d0) parameter(one = 1.0d0) parameter(mone = -1.0d0) parameter(big = 1.0d+20) parameter(maxit = 500) # Initialization: We try to follow the notation of LMS # On input we require: # # c1 = n1-vector of marginal costs (-y in the rq problem) # a1 = p by n1 matrix of equality constraints (x' in rq) # c2 = n2-vector of marginal costs (-r in the rq problem) # a2 = p by n2 matrix of inequality constraints (R' in rq) # b = p-vector of rhs ((1-tau)x'e in rq) # u = upper bound vector ( 1_p in rq) # beta = barrier parameter, LMS recommend .99995 # eps = convergence tolerance, LMS recommend 10d-8 # # the integer vector nit returns iteration counts # the integer info contains an error code from the Cholesky in stepy # info = 0 is fine # info < 0 invalid argument to dposv # info > 0 singular matrix nit(1)=0 nit(2)=0 nit(3)=n1 # Start at the OLS estimate for the dual vector y call dgemv('N',p,n1,one,a1,p,c1,1,zero,y,1) do i=1,n1 d1(i)=one do i=1,n2{ d2(i)=zero z2(i)=one } call stepy2(n1,n2,p,a1,d1,a2,d2,y,ada,info) if(info != 0) return # put current residual vector r1 in s (temporarily) call dcopy(n1,c1,1,s,1) call dgemv('T',p,n1,mone,a1,p,y,1,one,s,1) # Initialize remaining variables do i=1,n1{ if(dabs(s(i)) < eps){ z1(i)=dmax1(s(i),zero)+eps w(i)=dmax1(-s(i),zero)+eps } else { z1(i)=dmax1(s(i),zero) w(i)=dmax1(-s(i),zero) } s(i)=u(i)-x1(i) } gap = ddot(n1,z1,1,x1,1)+ddot(n2,z2,1,x2,1)+ddot(n1,w,1,s,1) # # Main Loop # while(gap > eps && nit(1)