subroutine rqfnb(n,p,a,y,rhs,d,u,beta,eps,wn,wp,nit,info) integer n,p,info,nit(3) double precision a(p,n),y(n),rhs(p),d(n),u(n),wn(n,9),wp(p,p+3) double precision one,beta,eps parameter( one = 1.0d0) call lpfnb(n,p,a,y,rhs,d,u,beta,eps,wn(1,1),wn(1,2), wp(1,1),wn(1,3),wn(1,4),wn(1,5),wn(1,6), wp(1,2),wn(1,7),wn(1,8),wn(1,9),wp(1,3),wp(1,4),nit,info) return end # This is a revised form of my primal-dual log barrier form of the # interior point LP solver based on Lustig, Marsten and Shanno ORSA J Opt 1992. # 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 c'x s.t. Ax=b, 0<=x<=u # # Denote dx,dy,dw,ds,dz as the steps for the respective variables x,y,w,s,z # subroutine lpfnb(n,p,a,c,b,d,u,beta,eps,x,s,y,z,w, dx,ds,dy,dz,dw,dr,rhs,ada,nit,info) integer n,p,pp,i,info,nit(3),maxit double precision a(p,n),c(n),b(p) double precision zero,one,mone,big,ddot,dmax1,dmin1,dxdz,dsdw double precision deltap,deltad,beta,eps,mu,gap,g double precision x(n),u(n),s(n),y(p),z(n),w(n),d(n),rhs(p),ada(p,p) double precision dx(n),ds(n),dy(p),dz(n),dw(n),dr(n) parameter( zero = 0.0d0) parameter( one = 1.0d0) parameter( mone = -1.0d0) parameter( big = 1.0d+20) parameter( maxit = 500) # Initialization: We follow the notation of LMS # On input we require: # # c = n-vector of marginal costs (-y in the rq problem) # a = p by n matrix of linear constraints (x' in rq) # b = p-vector of rhs ((1-tau)x'e in rq) # u = upper bound vector ( e 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)=n pp=p*p # Start at the OLS estimate for the dual vector y call dgemv('N',p,n,one,a,p,c,1,zero,y,1) do i=1,n d(i)=one call stepy(n,p,a,d,y,ada,info) if(info != 0) return # put current residual vector in s (temporarily) call dcopy(n,c,1,s,1) call dgemv('T',p,n,mone,a,p,y,1,one,s,1) # Initialize remaining variables do i=1,n{ if(dabs(s(i)) eps && nit(1) A'dy - ds deltap=big deltad=big do i=1,n{ dx(i)=d(i)*ds(i) ds(i)=-dx(i) dz(i)=-z(i)*(dx(i)/x(i) + one) dw(i)=-w(i)*(ds(i)/s(i) + one) if(dx(i)<0)deltap=dmin1(deltap,-x(i)/dx(i)) if(ds(i)<0)deltap=dmin1(deltap,-s(i)/ds(i)) if(dz(i)<0)deltad=dmin1(deltad,-z(i)/dz(i)) if(dw(i)<0)deltad=dmin1(deltad,-w(i)/dw(i)) } deltap=dmin1(beta*deltap,one) deltad=dmin1(beta*deltad,one) if(min(deltap,deltad) < one){ nit(2)=nit(2)+1 # Update mu mu = ddot(n,x,1,z,1)+ddot(n,s,1,w,1) g = mu + deltap*ddot(n,dx,1,z,1)+ deltad*ddot(n,dz,1,x,1) + deltap*deltad*ddot(n,dz,1,dx,1)+ deltap*ddot(n,ds,1,w,1)+ deltad*ddot(n,dw,1,s,1) + deltap*deltad*ddot(n,ds,1,dw,1) mu = mu * ((g/mu)**3) /dfloat(2*n) # Compute modified step do i=1,n{ dr(i)=d(i)*(mu*(1/s(i)-1/x(i))+ dx(i)*dz(i)/x(i)-ds(i)*dw(i)/s(i)) } call dswap(p,rhs,1,dy,1) call dgemv('N',p,n,one,a,p,dr,1,one,dy,1)# new rhs call dpotrs('U',p,1,ada,p,dy,p,info)# backsolve for dy call dgemv('T',p,n,one,a,p,dy,1,zero,u,1)#ds=A'ddy deltap=big deltad=big do i=1,n{ dxdz = dx(i)*dz(i) dsdw = ds(i)*dw(i) dx(i)= d(i)*(u(i)-z(i)+w(i))-dr(i) ds(i)= -dx(i) dz(i)= -z(i)+(mu - z(i)*dx(i) - dxdz)/x(i) dw(i)= -w(i)+(mu - w(i)*ds(i) - dsdw)/s(i) if(dx(i)<0)deltap=dmin1(deltap,-x(i)/dx(i)) if(ds(i)<0)deltap=dmin1(deltap,-s(i)/ds(i)) if(dz(i)<0)deltad=dmin1(deltad,-z(i)/dz(i)) if(dw(i)<0)deltad=dmin1(deltad,-w(i)/dw(i)) } deltap=dmin1(beta*deltap,one) deltad=dmin1(beta*deltad,one) } call daxpy(n,deltap,dx,1,x,1) call daxpy(n,deltap,ds,1,s,1) call daxpy(p,deltad,dy,1,y,1) call daxpy(n,deltad,dz,1,z,1) call daxpy(n,deltad,dw,1,w,1) gap = ddot(n,z,1,x,1)+ddot(n,w,1,s,1) } # return residuals in the vector x call daxpy(n,mone,w,1,z,1) call dswap(n,z,1,x,1) return end subroutine stepy(n,p,a,d,b,ada,info) integer n,p,pp,i,info double precision a(p,n),b(p),d(n),ada(p,p),zero parameter( zero = 0.0d0) # Solve the linear system ada'x=b by Choleski -- d is diagonal # Note that a isn't altered, and on output ada returns the upper # triangle Choleski factor, which can be reused, eg with blas dtrtrs pp=p*p do j=1,p do k=1,p ada(j,k)=zero do i=1,n call dsyr('U',p,d(i),a(1,i),1,ada,p) call dposv('U',p,1,ada,p,b,p,info) return end