• Arts
• Business
• Computers
• Games
• Health
• Home
• Kids and Teens
• News
• Recreation
• Reference
• Regional
• Science
• Shopping
• Society
• Sports

Linear Programming

# To unbundle, sh this file echo fnc.r 1>&2 cat >fnc.r m){if(mm>n){m=2*m;break}} else return } } return end #timing function real function udt(t) real s,dtime,udt,u(2) s=dtime(u) udt=u(1) return end #This is my second attempt to code the primal-dual log barrier form of the #interior point LP solver of 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. #The primary difference between this code and the previous version fna.r is #that we assume a feasible starting point so p,d,b feasibility gaps = 0. #Note also that all variables are assumed to have upper bounds of one. #Problem: # min c'x s.t. Ax=b, 0 eps) { nit(1)=nit(1)+1 do i =1,n{ d(i) = one/(z(i)/x(i) + w(i)/s(i)) ds(i)=z(i)-w(i) dx(i)=d(i)*ds(i) } call dgemv('N',p,n,one,a,p,dx,1,zero,dy,1)#rhs call dcopy(p,dy,1,rhs,1)#save rhs call stepy(n,p,a,d,dy,ada) call dgemv('T',p,n,one,a,p,dy,1,mone,ds,1) 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)*(dx(i)/s(i) - one) dxdz(i)=dx(i)*dz(i) dsdw(i)=ds(i)*dw(i) if(dx(i)1) mu=(g/dfloat(n))*(g/acomp)**2 #else mu=acomp/(dfloat(n)**2) do i=1,n{ dz(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,dz,1,one,dy,1)#new rhs call dpotrs('U',p,1,ada,p,dy,p,info) call daxpy(p,mone,dy,1,rhs,1)#rhs=ddy call dgemv('T',p,n,one,a,p,rhs,1,zero,dw,1)#dw=A'ddy deltap=big deltad=big do i=1,n{ dx(i)=dx(i)-dz(i)-d(i)*dw(i) ds(i)=-dx(i) dz(i)=mu/x(i) - z(i)*dx(i)/x(i) - z(i) - dxdz(i)/x(i) dw(i)=mu/s(i) - w(i)*ds(i)/s(i) - w(i) - dsdw(i)/s(i) if(dx(i)&2 cat >glob.r -band){ #inside the band points s(nxt)=0 nxt=nxt+1 } else{ #below the band points call daxpy(p,one,a(1,nxt),1,glob,1) call dswap(p,a(1,nxt),1,a(1,lxt),1) blo=bl

Contact Information