module iolog save integer :: iu=20,lo=21 integer :: logu real :: time end module iolog module nodbb save integer :: nvare,listl,listsa,noint,idepth,nsolen integer :: ipart(1:26000) real :: xincva,xival,revbnd real :: dfpart(1:26000),xinc(1:55000),pu(1:26000),pd(1:26000),pg(1:26000) character:: qfix*8 type nudo_bb real :: xiobnd,vbnd integer :: ivid1,ivid2 end type nudo_bb type (nudo_bb) :: nudo(1:12500) end module nodbb module defpa save integer :: nrmax=26000,ncolma=55000,nemax=4050000 integer :: ntmax=500000,maxnod=12500,namax=299500 integer :: itsinv=9999999 integer :: iuno=1,idos=2 real :: epmcd =2.2204e-16 real :: ztolze=1.1921e-7 real :: ztolpv=1.4901e-8 real :: ztcost=2.4222e-5 real :: toleta=1.1921e-7 real :: plinfy=1.e20,cero=0 real :: uno=1,dos=2 end module defpa module parlp save integer :: jh(1:26000),la(1:55001),ia(1:299500),le(1:4050001),ie(1:4050000) integer :: mreg(1:26000),hreg(1:26000),vreg(1:26000) integer :: kinbas(1:55000),kinben(1:55000),hrtype(1:26000) integer :: inibb,invfrq,nrow,ncol,itdmax,itpmax integer :: nelem,neta,nleta,nueta,ninf,minomax,nopt integer :: logite,nrwm1 real :: a(1:299500),x(1:26000),y(1:26000),b(1:26000),yen(1:26000) real :: xci(1:55000),xlb(1:55000),xub(1:55000),xcs(1:55000),e(1:4050000) real :: suminf,rcost,rpivot,itrfrq,itcnt character:: qname(55000)*8,qnampo*8,qstat*6 end module parlp PROGRAM Bbmi ! !****************************************************************** ! ! BBMI resuelve problemas de programación lineal, !C entera y programación mixta lineal-entera. !C !C !C E S T R U C T U R A d e S U B R U T I N A S !C !C ------------------------------------------------------- !C | Nivel0 | Nivel1 | Nivel2 | Nivel3 | Nivel4 | Nivel5 | !C ------------------------------------------------------- !C | | | | | | | !C | | !C | / INPUT | !C | | | !C | | | !C | | / / WRETA | !C | | | INVERT < | !C | | | \ FTRAN | !C | | | | !C | | NORMAL < FORMC | !C | | | BTRAN | !C | | | PRICE | !C | | | FTRAN | !C | | | CHUZR | !C | | \ WRETA | !C | | | !C | | WRILPR | !C | | | !C | | / | BTRAN | !C | | | TESTX < | !C | | | | WRILPR | !C | | | | !C | | | | FTRAN | !C | | | PENLTS < | !C | | | | BRANCH | !C | | | | !C | | | / | !C | | | | BKTRAK < FTRAN | !C | | | | | !C | | | | | !C | BBMI < | | / / WRETA | !C | | | | | INVERT < | !C | | | | | \ FTRAN | !C | | | | | | !C | | | CYCLE < NORMAL < FORMC | !C | | | | | BTRAN | !C | | BANDB < | | PRICE | !C | | | | | FTRAN | !C | | | | | CHUZR | !C | | | | \ WRETA | !C | | | | | !C | | | | | BTRAN | !C | | | | TESTX < | !C | | | \ | WRILPR | !C | | | | !C | | | DCHUZC < FTRAN | !C | | | | !C | | | CYCLE | !C | | | INVERT | !C | | | WRETA | !C | | \ DCHUZR | !C | | | !C | \ WRILPR | !C | | !C | | | | | | | !C ------------------------------------------------------- !C !C !C****************************************************************** ! use iolog; use defpa; use nodbb; use parlp implicit none call input ! Entrada de datos time=cputim() call normal ! Se resuelve programa lineal original call wrilpr (iuno) if (nvare<=ncol) then call bandb ! Problema con variables entereas call wrilpr (3) ! Se resuelve con branch-and-bound endif write (lo,9000) cputim()-time 9000 format(/'Tiempo total de CPU:',f9.2,' segundos') contains subroutine bandb ! ! Rutina que dirige el proceso de optimización del programa ! entero o mixto mediante un procedimiento Branch-and-Bound ! con relajaciones lineales. !****************************************************************** ! use iolog; use defpa; use nodbb; use parlp ! implicit none save integer :: irowpd,idir,iptypd,icolit,iax,jcolpd,itdual character*8 :: qupdo,qplus ! if (itcnt>=itrfrq) return listl=0; listsa=0; idepth=0; itdual=0; call testx ! Se comprueba si la solución inicial es entera if (qstat=='ENTERA') return inibb = 1 ! Indicador de estar en proceso de optimi. entera bucle1: do ! ! Nudo analizándose no rechazado; calcular penalizaciones por separar una u otra variable. ! Las nuevas ramas del árbol se generan en BRANCH a la que llama PENLTS. ! Si en PENLTS se rechaza el nudo (IDIR=0), buscar en BKTRAK el siguiente a analizar; ! si no se desecha, utilizar el método dual del simplex para reoptimizar. ! call penlts (irowpd,idir,iptypd) if (idir==0) then if (logite==1) write (logu,550) call bucle_bktrak; if (listl==0) return cycle bucle1 endif icolit = jh(irowpd) if (icolit<=nrow) then icolit=icolit+ncol-nrow-1 else icolit=icolit-nrow endif qupdo = ' ' if (idir==1) then iax=ipart(irowpd) write (qupdo,10) 'X<',iax else iax=ipart(irowpd)+1 write (qupdo,10) 'X>',iax endif dual_simplex: do itdual=1,itdmax ! *** Ciclo iterativo dual del simplex. call dchuzc (jcolpd,irowpd,iptypd) ! Variable no básica jcolpd a entrar en la base itcnt =itcnt+1 itsinv=itsinv+1 if (itdual==itdmax)jcolpd=0 if (jcolpd==0) then idepth=idepth+1 if (logite==1) then write (logu,600) icolit,qfix,idepth, & qupdo,listsa,noint,int(itcnt),' D','no factible' ! else ! if (itdual==itdmax) write (21,*) 'D-inf' endif call bucle_bktrak; if (listl==0) return cycle bucle1 endif if (nelem>nemax.or.itsinv>=invfrq) then ! ¿Hay que reinvertir la base? call invert itsinv=0 else call wreta (irowpd) endif if (itcnt>=itrfrq) then write (lo,'("*** Nº iteraciones sobrepasado; el programa se para")') call wrilpr (3) write (lo,9000) cputim()-time stop endif call dchuzr (irowpd,iptypd) ! Escoger variable básica irowpd fuera de límites ! que saldrá de la base if (irowpd==0) exit dual_simplex enddo dual_simplex idepth = idepth+1 if (x(1)<=xincva+ztolze) then qplus = 'ABANDONO' if (logite==1) write (logu,700) icolit,qfix,idepth,qupdo, & listsa,noint,int(itcnt),' D',x(1)*minomax,qplus call bucle_bktrak; if (listl==0) return cycle bucle1 endif qplus = ' ' if (logite==1) write (logu,700) icolit,qfix,idepth,qupdo, & listsa,noint,int(itcnt),' D',x(1)*minomax,qplus call testx if (qstat=='ENTERA') then call bucle_bktrak; if (listl==0) return endif enddo bucle1 10 format(a,i4) 550 format('----- Nudo desechado en PENLTS -----') 600 format(i6,a8,i5,7x,a,i7,4x,i7,8x,i4,a,5x,a) 700 format(i6,a8,i5,7x,a,i7,4x,i7,8x,i4,a,1pd18.7,1x,a) 9000 format(/'Tiempo total de CPU:',f9.2,' segundos') end subroutine bandb subroutine bucle_bktrak ! ! Esta rutina lleva a cabo la resolución de la relajación ! lineal del nudo determinado en BKTRAK. Se utiliza ! el método primal del simplex. !****************************************************************** ! use iolog; use nodbb; use parlp implicit none integer :: iax,inco,icolit,iabs character*8 :: qupdo,qplus do call bktrak if (listl==0) return ! Lista de nudos a analizar vacia; se termina call normal ! Si no esta vacía, optimizar nudo mediante simplex primal. inco = nudo(listl)%ivid1 icolit= iabs(inco) ! Se imprimen resultados intermedios si ha lugar qupdo = ' ' if (inco>0) then iax=int(xub(icolit)) write (qupdo,10) 'X<',iax else iax=int(xlb(icolit)) write (qupdo,10) 'X>',iax endif if (icolit<=nrow) then icolit=icolit+ncol-nrow-1 else icolit=icolit-nrow endif idepth=idepth+1 if (qstat=='FEAS ') then if (x(1)<=xincva+ztolze) then ! Función obj. peor que actual qplus = 'ABANDONO' if (logite==1) write (logu,20) icolit,qfix,idepth,qupdo, & listsa,noint,int(itcnt),' P',x(1)*minomax,qplus,cputim() cycle endif qplus = ' ' if (logite==1) write (logu,20) icolit,qfix,idepth,qupdo, & listsa,noint,int(itcnt),' P',x(1)*minomax,qplus,cputim() else if (logite==1) write (logu,30) icolit,qfix,idepth,qupdo, & listsa,noint,int(itcnt),' P',cputim() cycle endif call testx ! Comprobar si solución conseguida es entera factible. if (qstat/='ENTERA') return ! No; retornar enddo ! Si; continuar con otro nudo 10 format(a,i4) 20 format(i6,a8,i5,7x,a,i7,4x,i7,8x,i4,a,1pd18.7,1x,a,f10.2) 30 format(i6,a8,i5,7x,a,i7,4x,i7,8x,i4,a,5x,'no factible',f10.2) end subroutine bucle_bktrak subroutine bktrak ! ! Se escoge un nudo del árbol enumerativo siguiendo la ! regla LIFO (Last In First Out). !****************************************************************** ! use iolog; use defpa; use nodbb; use parlp implicit none integer :: ntemp3,inco,i,j,ir,icolit,iax real :: dp,de,dy character*8 :: qupdo ntemp3=0; qfix=' ' do if (listl==0) return ! Se termina cuando la lista de nudos vacía if (nudo(listl)%xiobnd>xincva+ztolze.and.nudo(listl)%xiobnd/=-1.e50) then inco = nudo(listl)%ivid1 if (inco>=0) then dp =xlb(inco) xlb(inco)=xub(inco)+uno xub(inco)=nudo(listl)%vbnd if (kinbas(inco)<=0) then kinbas(inco)=0 ntemp3 =1 endif else inco =-inco dp =xub(inco) xub(inco)=xlb(inco)-uno xlb(inco)=nudo(listl)%vbnd if (kinbas(inco)<=0) then kinbas(inco)=-1 ntemp3 = 1 endif endif nudo(listl)%ivid1 =-nudo(listl)%ivid1 idepth = nudo(listl)%ivid2 nudo(listl)%vbnd =dp nudo(listl)%xiobnd=-1.e50 listsa =listsa-1 if (ntemp3==0) return ! -1 -1 ! Calcular la solución que determina nudo LISTL: x = B b - B Nx ! N forall (i=1:nrow) y(i)=b(i) do j = 1,ncol if (kinbas(j)==-1) then de=xub(j) elseif (kinbas(j)==0) then de=xlb(j) else cycle endif do i=la(j),la(j+1)-1 ir =ia(i) y(ir)=y(ir)-a(i)*de end do end do call ftran (1) forall (i=1:nrow) x(i)=y(i) return else ! Nudo no interesa; adaptar límites de las variables y escoger otro. inco = nudo(listl)%ivid1 icolit = iabs(inco) if (nudo(listl)%xiobnd/=-1.e50) then listsa=listsa-1 inco =nudo(listl)%ivid1 qupdo =' ' if (inco>0) then iax=int(xub(icolit))+1 write (qupdo,10) 'X>',iax else iax=int(xlb(icolit))-1 write (qupdo,10) 'X<',iax endif if (icolit<=nrow) then icolit=icolit+ncol-nrow-1 else icolit=icolit-nrow endif if (logite==1) write (logu,550) icolit,qfix,nudo(listl)%ivid2+1, & qupdo,listsa,nudo(listl)%xiobnd*minomax,cputim() endif if (inco>=0) then if (kinbas(inco)==-1) then ntemp3=1 dp =xub(inco)-xlb(inco) dy =nudo(listl)%vbnd-xub(inco) if (dpmaxnod) stop 'El número de nudos excede el establecido; incrementar "MAXNOD"' if (idir==-1) nudo(listl)%vbnd=xlb(icol) if (idir== 1) nudo(listl)%vbnd=xub(icol) nudo(listl)%ivid1 =idir*icol nudo(listl)%ivid2 =idepth nudo(listl)%xiobnd=xival end subroutine branch subroutine btran ! ! "Backward transformation". t -1 ! Se calculan los multiplicadores simplex pi=c B realizando ! B ! transformación inversa del vector c ; es decir ! B ! t ! pi=(((c E )E )...E ) ! B k) k-1 1 ! !******* Descripción de vectores afectados ***** ! Y := Vector de costes c ! B ! ! LE, IE, E := Vectores que definen en forma dispersa las ! matrices elementales eta. ! !****************************************************************** ! use defpa; use parlp implicit none integer :: i,j,ipiv,kk,ll real :: dsum if (neta<=0) return do i = neta,1,-1 ll =le(i) kk =le(i+1)-1 ipiv=ie(ll) dsum=cero if (kk>ll) then do j=ll+1,kk dsum=dsum+e(j)*y(ie(j)) end do endif y(ipiv)=(y(ipiv)-dsum)/e(ll) end do end subroutine btran subroutine chuzr(jcolp,irowp,npivot,iptype,theta,ivout) ! ! Esta rutina, una vez efectuada la transformación directa ! -1 ! Y = B a ! jcolp ! del vector columna de la matriz de condiciones A, a , ! jcolp ! que entrará en la base, selecciona la fila IROWP sobre la ! que pivotará dicha columna en esta iteración del método ! simplex; selecciona de esta forma la variable básica X ! que ha de salir de la base. irowp ! !******* Descripción de parámetros ***** ! KINBAS := Estado en el que se encuentran los componentes del ! del vector X de variables de decisión: ! =0, variable en su límite inferior; ! =-1, variable en su límite superior; ! =-2, variable libre: especificada FR; ! =>0, variable básica, el número indica ! sobre qué vector fila pivota. ! ! HRTYPE := Se calcula en FORMC de tal forma que: ! =-2 si X(j) < XLB(j)-ZTOLZE; ! = 0 si X(j) es factible y ! = 2 si X(j) > XUB(j)+ZTOLZE. ! ! XLB, XUB := Límites inferior y superior de las variables. ! ! JH := Vector que indica la columna con la que pivota ! cada vector fila. ! ! X := Valor que en esta iteración toma el vector de ! variables de decisión. ! -1 ! Y := El ya mencionado B a ! jcolp ! !********************************************************************* ! use iolog; use defpa; use nodbb use parlp ! implicit none integer :: jcolp,is,i,j,inco,jtype,jhit1,jhit2 integer :: irowp,ivout,npivot,iptype real :: ynorm,stepmn,stepmx,theta1,theta2,pertbn,tol real :: pivot,pivabs,res,pthet1,pthet2,theta real :: pivmx2,pivmx1,sum logical :: hlow1,hlow2,hitlow,move,unbndd ! jh(nrwm1) = jcolp ! ! Sabiendo que se está maximizando una función objetivo, ! comprobar si la variable a entrar en la base se incrementa ! desde su límite inferior, RCOST>0, o decrementa desde ! su superior, RCOST<0. ! if (rcost>cero) then x(nrwm1) = xlb(jcolp) if (kinbas(jcolp)==-2) x(nrwm1)=cero y(nrwm1) = -uno is = 1 else x(nrwm1) = xub(jcolp) if (kinbas(jcolp)==-2) x(nrwm1)=cero y(nrwm1) = uno is = -1 forall (i=1:nrow) y(i)=-y(i) endif sum=cero do i=2,nrow sum=sum+abs(y(i)) end do ynorm =sum/(nrow-1) stepmn=ztolpv/(uno+ynorm); stepmx=1.e50/(uno+ynorm) ! Antes 1.e11 theta1=stepmx theta2=cero pertbn=1.1*ztolze tol =ztolpv*ynorm ! ! PERTBN debe ser mayor que ZTOLZE. ! ! Primera pasada. Se perturba RES, la distancia a cada límite o ! cota, de tal forma que THETA1 sea ligeramente mayor que el ! verdadero paso a dar y THETA2 ligeramente inferior. En casos ! de degeneración, esta estrategia pemite cierta libertad en una ! segunda pasada. ! do j=2,nrwm1 inco =jh(j) pivot =y(j) pivabs=abs(pivot) if (pivabs>tol) then jtype=hrtype(j) if (pivot>cero) then ! ! La posible variable a entrar en la base decrecería. ! Chequear un THETA1 más pequeño si se satisface su límite ! inferior. ! if (jtype>=0) then res=x(j)-xlb(inco)+pertbn if (theta1*pivot>res) theta1=res/pivot ! ! Comprobar si existe un THETA2 mayor si se viola su límite ! superior. ! if (jtype>0) then res=x(j)-xub(inco)-pertbn if (theta2*pivotres) theta1=res/pivabs ! ! Comprobar si existe un THETA2 mayor si se viola su límite ! inferior. ! if (jtype<0) then res=xlb(inco)-x(j)-pertbn if (theta2*pivabstol) then jtype=hrtype(j) if (pivot>cero) then ! ! La posible variable a entrar en la base decrecería. ! Chequear un THETA1 más pequeño si se satisface su límite ! inferior. ! if (jtype>=0) then res = x(j)-xlb(inco) if (pthet1*pivot>=res.and.pivmx10) then res = x(j)-xub(inco) if (pthet2*pivot<=res.and.pivmx2=res.and.pivmx1theta2) then theta =theta2 irowp =jhit2 hitlow=hlow2 endif ! move =theta>stepmn unbndd=theta>=stepmx.or.irowp==0 if (.not.unbndd) then if (hitlow) then; unbndd=xlb(jh(irowp))<-plinfy else; unbndd=xub(jh(irowp))>plinfy endif endif if (.not.move) theta=cero if (unbndd) stop 'Problema no acotado' ! ! Se adapta el vector solución al final de una pivotacióón del ! método simplex. ! ! x =x -THETA*y ; 1<=i<=m. ! j j i ! i i ! if (irowp==nrwm1) then ! ! No ha habido pivotación; una variable no básica va de uno ! de sus límites al otro. Cambiar su estatus: si = -1, ! hacerlo 0; si = 0, hacerlo -1. ! kinbas(jcolp)=-(kinbas(jcolp)+1) ivout =jcolp rpivot =cero npivot =0 else npivot =1 ! Pivotación normal. rpivot =y(irowp)*is ivout =jh(irowp) kinbas(jcolp)=irowp if (hitlow) then kinbas(ivout)= 0 iptype=0 else kinbas(ivout)= -1 iptype=-1 endif jh(irowp)=jcolp endif if (move) forall (i=1:nrwm1) x(i)=x(i)-y(i)*theta ! Adaptar el vector x. x(irowp)=x(nrwm1) if (rcostztolze.and.(kinbas(j)==0.or.kinbas(j)==-1)) then forall (i=1:nrow) y(i)=cero; forall (i=la(j):la(j+1)-1) y(ia(i))=a(i) ! Unpack en Y() call ftran (iuno) select case (kinbas(j)) case(0) ! La candidata no básica j en su límite inferior if (y(irowp)>ztolpv) then de = y(1)/y(irowp) if (deztolze.and.(kinbas(j)==0.or.kinbas(j)==-1)) then forall (i=1:nrow) y(i)=cero; forall (i=la(j):la(j+1)-1) y(ia(i))=a(i) ! Unpack en Y() call ftran (iuno) select case (kinbas(j)) case (-1) ! La candidata no básica j en su límite superior if (y(irowp)>ztolpv) then de = y(1)/y(irowp) if (de>dp) then jcolp=j; dp=de endif endif case (0) ! La candidata no básica j en su límite inferior if (y(irowp)<-ztolpv) then de = y(1)/y(irowp) if (de>dp) then jcolp=j; dp=de endif endif end select endif end do if (jcolp==0) return ! ! -1 ! Se calcula B N y el paso ! jcolp ! forall (i=1:nrow) y(i)=cero; forall (i=la(jcolp):la(jcolp+1)-1) y(ia(i))=a(i) ! Unpack en Y() call ftran (iuno) theta=(x(irowp)-xlb(jh(irowp)))/y(irowp) endif forall (i=1:nrow) x(i)=x(i)-y(i)*theta if (kinbas(jcolp)==-1) then ! Se recalcula el valor de la nueva variable básica jcolp x(irowp)=xub(jcolp)+theta else x(irowp)=xlb(jcolp)+theta endif kinbas(jcolp) =irowp ! La variable jcolp pasa o ocupar la posición de irowp kinbas(jh(irowp))=iptype ! La antigua irowp para a su límite sup. o inf. jh(irowp) =jcolp ! write(21,'(a,e12.5,i4)')'I-d, t=',theta,jcolp end subroutine dchuzc subroutine dchuzr(irowp,iptype) ! ! Se determina para una iteración del método simplex dual la ! variable x que ha de salir de la base. ! irowp ! ! Si el problema es factible en el primal se acaba; se ha ! llegado al óptimo haciéndose IROWP=0. ! !****************************************************************** ! use defpa; use nodbb; use parlp implicit none integer :: irowp,i,inco,iptype,ninfdual real :: dp,de,suminf ! ! Se escoge la fila con mayor no factibilidad. ! irowp =0 qstat ='FEAS ' dp =-plinfy ninfdual=0 suminf =0 do i = 2,nrow inco=jh(i) if (x(i)dp) then dp =de iptype=0 irowp =i endif elseif (x(i)>xub(inco)+ztolze) then ! ! La variable básica que pivota en la fila I esta por encima de ! su límite superior ! qstat ='INFEAS' ninfdual=ninfdual+1 de =x(i)-xub(inco) suminf =suminf+de if (de>dp) then dp =de iptype=-1 irowp = i endif endif end do end subroutine dchuzr subroutine formc ! t -1 ! Se forma en Y el vector de coeficientes c de la fun. objetivo para la expresión -c B a . ! B B j ! ! En la Phase II es Y(1)=1 y todos los demás componentes cero. ! En Phase I, Y(1)=0 añadiéndose una variable artificial Y(i)=1 si x u (recordemos que el problema es de maximización), ! i i ! por lo que los costes se incluyen se incluyen en c ! B. ! !******* Descripción de vectores afectados ***** ! XLB, XUB := Límites inferior y superior de las variables ! del problema. ! ! JH := Vector que indica la columna con la que pivota ! cada vector fila. ! ! -1 ! X := Vector solución actual =B b. ! ! Y := Vector de coste a formar. ! ! QSTAT := Estado de la solución = 'INFEAS' si no es factible; ! = 'FEAS' si es factible. ! !****************************************************************** ! use defpa; use parlp implicit none integer :: i,icol qstat='INFEAS'; suminf=cero; y(1)=cero; ninf=0 do i=2,nrow icol=jh(i) if (x(i)>xlb(icol)-ztolze) then if (x(i)nle) return do ik=nfe,nle ll=le(ik); kk=le(ik+1)-1; ipiv=ie(ll); y(ipiv)=y(ipiv)/e(ll); yip=y(ipiv); if (kk>ll) forall (j=ll+1:kk) y(ie(j))=y(ie(j))-e(j)*yip end do end subroutine ftran subroutine input ! ! Entrada de los datos del problema a resolver. ! !******* Descripción de vectores afectados ***** ! ! B := Término de la derecha del problema. ! ! KINBAS := Estado en el que se encuentran los componentes ! del vector X de variables de decisión: ! = 0, variable en su límite inferior; ! =-1, variable en su límite superior; ! =-2, variable libre: especificada FR; ! =>0, variable básica, el número indica ! el vector fila sobre el que pivota. ! ! QNAME := Vector CHARACTER donde se guardan los nombres ! atribuidos en el fichero MPS a las variables ! de decisión y a las condiciones. ! ! XLB, XUB := Límites inferior y superior de las variables. ! ! JH := Vector que indica la columna con la que pivota ! cada vector fila. ! ! LA, IA, A := Vectores donde se guarda toda la información ! relativa a la matriz de los coeficientes ! de las condiciones. ! ! QNAMPO := Variable CHARACTER con el nombre del caso a ! resolver. !****************************************************************** ! use iolog; use defpa; use nodbb; use parlp implicit none integer :: initbd,ind,l,j,i real :: vtemp1,vtemp2,dabs character :: idf*20,char1*80,qn(3)*8,qcs*8 character*1 :: k1,k2,k3,k4 minomax= -1; inibb = 0; nrow =1; itcnt = 0; nvare = 0 invfrq =100; itrfrq =2.5e8; logite=0; qcs = ' ' a(1) =uno; b(1) = cero; ia(1) =1; la(1) = 1 jh(1) = 1; kinbas(1)= 1; xincva=-1.e75 nsolen = 0; initbd = 0; l =0; itdmax=50; itpmax=500 write (*,'(A)') 'Fichero de datos del problema?' read (*,'(A)') idf open (iu,file=idf) open (lo,file=idf(:index(idf,' ')-1)//'.pbb') do ! Leer las especificaciones del proceso de optimización. read (iu,'(A)') char1 if (index(char1,'NAME')==1) exit if (index(char1,'MAXIMIZAR')/=0) then minomax = 1 cycle endif ind = index(char1,'FRECUENCIA DE REINVERSION') if (ind/=0) then read (char1(ind+25:),'(BN,I50)') invfrq cycle endif ind = index(char1,'LIMITE DE ITERACIONES') if (ind/=0) then read (char1(ind+21:),'(BN,f50.0)') itrfrq cycle endif ind = index(char1,'ITERACIONES DUAL-SIMPLEX') if (ind/=0) then read (char1(ind+24:),'(BN,i10)') itdmax cycle endif ind = index(char1,'COTA INICIAL F.O.') if (ind/=0) then read (char1(ind+17:),'(BN,I50)') initbd cycle endif ind = index(char1,'PROCESO ITERATIVO') if (ind/=0) then read (char1(ind+17:),'(BN,I50)') logu if (logu/=0) then open(logu,file=idf//'.LOG') else logu = lo endif logite = 1 cycle endif ind = index(char1,'TOLERANCIA DE CERO') if (ind/=0) then read (char1(ind+18:),'(bn,d50.0)') ztolze cycle endif ind = index(char1,'TOLERANCIA DE COSTES') if (ind/=0) then read (char1(ind+20:),'(bn,d50.0)') ztcost cycle endif ind = index(char1,'TOLERANCIA DE PIVOTE') if (ind/=0) then read (char1(ind+20:),'(bn,d50.0)') ztolpv cycle endif enddo if (initbd/=0) xincva = dble(initbd)*minomax backspace iu Principal: do ! Leer el fichero de datos en formato MPS. read (iu,8000) k1,k2,k3,k4,qn(1),qn(2),vtemp1,qn(3),vtemp2 if (k1=='E'.and.k2=='N'.and.k3=='D'.and.k4=='A') exit Principal if (k1/=" ") then select case (k1//k2//k3) case ("NAM") qnampo = qn(2) case ("ROW") l = 1 case ("COL") l = 2 case ("RHS") l = 3 case ("BOU") l = 4 case ("RAN") l = 5 case ("INT") nvare = ncol+1 case ("* ") case default write (lo,'('' *** Clave '',A,'' no reconocida'')') k1//k2//k3//k4 stop end select cycle Principal endif select case (l) case (1) ! ROWS: nueva fila nrow = nrow+1 if (nrow>=nrmax) then write (lo,8554) nrmax stop endif qname(nrow) = qn(1) select case (k2//k3) ! Determinar tipo de fila case ("L") xlb(nrow) = cero; xci(nrow) = cero xub(nrow) = plinfy; xcs(nrow) = plinfy a(nrow) = uno case ("E") xlb(nrow) = cero; xci(nrow) = cero xub(nrow) = cero; xcs(nrow) = cero a(nrow) = uno case ("G") xlb(nrow) = cero; xci(nrow) = cero xub(nrow) = plinfy; xcs(nrow) = plinfy a(nrow) = -uno case ("N") nrow = nrow-1 qname(1) = qn(1) ncol = nrow cycle Principal end select b(nrow) =0.0; ia(nrow)=nrow; kinbas(nrow)=nrow la(nrow)=nrow; jh(nrow)=nrow nelem =nrow; ncol =nrow case (2) ! COLUMNS: coeficientes de las condiciones j = 2 if (abs(vtemp1)ncolma) then write (lo,8555) ncolma stop endif ncol = ncol+1 qcs = qn(1) qname(ncol) = qcs la(ncol) = nelem+1 xlb(ncol) = cero; xci(ncol) = cero xub(ncol) = plinfy; xcs(ncol) = plinfy kinbas(ncol) = 0 endif bucle_columns: do ! Comprobar si la fila existe. do i = 1,nrow if (qn(j)==qname(i)) then nelem = nelem+1 if (nelem>=namax) then write (lo,8550) namax stop endif if (qn(j)==qname(1).and.minomax==1) vtemp1 = -vtemp1 ia(nelem) = i a(nelem) = vtemp1 la(ncol+1) = nelem+1 if (j==3.or.abs(vtemp2)ztolze) then xub(i)=vtemp1 xcs(i)=vtemp1 a(i) =-1.0 else xub(i)=abs(vtemp1) xcs(i)=abs(vtemp1) endif else xub(i)=abs(vtemp1) xcs(i)=abs(vtemp1) endif if (j==3) cycle Principal j =3 vtemp1=vtemp2 if (qn(j)==' ') cycle Principal cycle bucle_ranges endif end do write (lo,8450) qn(j),qn(1) stop enddo bucle_ranges end select enddo Principal if (nvare==0) nvare = ncol+1 ! Fin de INPUT; se imprimen algunas estadísticas. write (lo,8650) qnampo write (lo,8700) nrow,nrow-1,ncol-nrow,nrow-1,nelem-nrow,ncol-nvare+1 & ,dble(nelem-nrow)*100./(nrow*(ncol-nrow)) nrwm1 = nrow+1 close(iu) if (itdmax== 50.and.nrow*2>itdmax) itdmax=nrow*2+10 if (itpmax==500.and.nrow*2>itdmax) itdmax=ncol+100 8000 format(4a1,a,2x,a,2x,f12.0,3x,a,2x,f12.0) 8250 format('*** Variable ',a,' previamente definida') 8300 format('*** La condición ',a,' en la línea COLUMN ',a,' no existe.') 8350 format('*** La condición ',a,' en la línea RHS ',a,' no existe.') 8400 format('¡Cuidado! La variable ',a,' de la línea BOUNDS ',a, & ' no existe o los coeficientes de la misma son nulos') 8450 format('*** La condición ',a,' en la línea RANGES ',a,' no existe.') 8550 format('*** El número de coeficientes de la matriz A',' excede de', & i7,'; el programa se para.') 8554 format('*** El número de condiciones del problema',' excede de',i7, & '; el programa se para.') 8555 format('*** El número de variables del problema',' excede de',i7, & '; el programa se para.') 8650 format('Problema ',a) 8700 format(//'*** Estadísticas del Problema'/4x,i7,' Fila(s)'/4x,i7, & ' Restricción(es)'/4x,i7,' Variable(s) de decisión'/4x,i7, & ' Variable(s) de ','holgura/artificiales'/4x,i7, & ' Elementos no cero'/4x,i7,' Variable(s) entera(s)'/4x, & 'Densidad de la matriz de coeficientes A:',f8.3,'%'/) end subroutine input subroutine invert ! ! Se refactoriza la inversa de la matriz básica B en la ! forma LU. ! !******* Descripción de vectores afectados ***** ! ! B := Término de la derecha del problema. ! ! KINBAS := Estado en el que se encuentran los componentes ! del vector X de variables de decisión: ! =0, variable en su límite inferior; ! =-1, variable en su límite superior; ! =-2, variable libre: especificada FR; ! =>0, variable básica, el número indica ! el vector fila sobre el que pivota. ! ! XLB, XUB := Límites inferior y superior de las variables. ! ! JH := Vector que indica la columna con la que pivota ! cada vector fila. ! ! LA, IA, A := Vectores donde se guarda toda la información ! relativa a la matriz de los coeficientes ! de las condiciones. ! ! -1 -1 -1 ! X := vector B b - B Nx = B (b-Nx ) ! ! LE, IE, E := Vectores donde se guarda toda la información ! relativa a la matriz de las ! transformaciones elementales eta. ! !****************************************************************** ! use iolog; use defpa; use parlp implicit none integer :: nlelem,nuelem,nabove,lr1,kr1,lr4,kr4,i,k integer :: lr3,ir,nbnonz,j,iv,ll,kk,ircnt,irp,nvrem integer :: irowp,ii,iir,imcnt,jj integer :: nslck,nbelow,nelast,ntlast,lr,jk integer :: ircmin,incr,nf,iaux,idif,nofd,nstr,inc real :: ztrel,ep,ztoly,ymax,yval,dmax1,ztmin,de real :: x2max,eemax,ermax,v,w ztrel = 0.05 10 continue neta =0; nleta =0; nueta =0; nelem=0 nlelem=0; nuelem=0; nabove=0; lr4 =nrwm1 le(1) =1; lr1 =1; kr1 =0; kr4 =nrow do i=1,nrow ! Se ponen las variables de holgura if (jh(i)<=nrow) then ! y artificiales en la parte 4, lr4 =lr4-1 ! el resto en la parte 1 mreg(lr4)=jh(i) vreg(lr4)=jh(i) else kr1 =kr1+1 vreg(kr1)=jh(i) endif hreg(i) = -1 jh(i) = 0 end do lr3=lr4 do i=lr4,kr4 ir =mreg(i) hreg(ir) =0 jh(ir) =ir kinbas(ir)=ir end do nbnonz = kr4-lr4+1 if (kr1/=0) then j=lr1 ! Recuperar vectores de debajo del pico do ! y contar filas. iv =vreg(j) ll =la(iv) kk =la(iv+1)-1 ircnt=0 do i=ll,kk nbnonz = nbnonz+1 ir = ia(i) if (hreg(ir)<0) then ircnt =ircnt+1 hreg(ir)=hreg(ir)-1 irp =ir endif end do select case (ircnt) case(0) write (lo,8000) kinbas(iv)=0 vreg(j) =vreg(kr1) kr1 =kr1-1 if (j>kr1) exit cycle case (1) vreg(j) =vreg(kr1) kr1 =kr1-1 lr3 =lr3-1 vreg(lr3) =iv mreg(lr3) =irp hreg(irp) =0 jh(irp) =iv kinbas(iv)=irp if (j>kr1) exit cycle case default if (j>=kr1) exit j = j+1 end select enddo ! ! Recuperar el resto de los vectores de encima y debajo del ! pico y establecer contador de mérito de las columnas. ! 310 nvrem = 0 if (kr1/=0) then j = lr1 do1: do iv =vreg(j) ll =la(iv) kk =la(iv+1)-1 ircnt=0 do i=ll,kk ir = ia(i) if (hreg(ir)==-2) then nabove=nabove+1 ! Pivota por encima del pico (parte de L). irowp =ir forall (k=1:nrow) y(k)=cero; forall (k=la(iv):la(iv+1)-1) y(ia(k))=a(k) ! Unpack en Y() call wreta (irowp) nleta =neta jh(ir) =iv kinbas(iv)=ir vreg(j) =vreg(kr1) kr1 =kr1-1 nvrem =nvrem+1 hreg(ir) =iv do ii=ll,kk ! Cambia contador de filas. iir=ia(ii) if (hreg(iir)<0) hreg(iir)=hreg(iir)+1 end do if (j>kr1) exit do1 cycle do1 elseif (hreg(ir)<0) then ircnt=ircnt+1 irp =ir endif end do select case (ircnt) case (0) write (lo,8000) kinbas(iv) = 0 vreg(j) = vreg(kr1) nvrem = nvrem+1 kr1 = kr1-1 if (j>kr1) exit do1 cycle do1 case (1) ! Pone vector debajo del pico. vreg(j) = vreg(kr1) nvrem = nvrem+1 kr1 = kr1-1 lr3 = lr3-1 vreg(lr3) = iv mreg(lr3) = irp hreg(irp) = 0 jh(irp) = iv kinbas(iv) = irp do ii = ll,kk ! Cambia contador de filas. iir = ia(ii) if (hreg(iir)<0) hreg(iir) = hreg(iir)+1 end do if (j>kr1) exit do1 cycle do1 case default if (j>=kr1) exit do1 j = j+1 end select enddo do1 if (nvrem>0) go to 310 if (kr1/=0) then do j = lr1,kr1 ! Obtener contador de meritos. iv = vreg(j) ll = la(iv) kk = la(iv+1)-1 imcnt = 0 do i = ll,kk ir = ia(i) if (hreg(ir)<0) imcnt=imcnt-hreg(ir)-1 end do mreg(j) = imcnt end do inc=1 ! Ordenar columnas segun mérito usando el algo. SHELL. do inc=2*inc+1 if (inc>kr1) exit enddo do inc=inc/2 do i=inc+1,kr1 v=mreg(i) w=vreg(i) j=i do if (mreg(i-inc)<=v) exit mreg(j)=mreg(j-inc) vreg(j)=vreg(j-inc) j=j-inc if (j<=inc) exit enddo mreg(j)=v vreg(j)=w enddo if (inc<=1) exit enddo endif endif endif nslck =0 ! Sacar etas debajo del pico (parte de U). nbelow =0 nelast =nemax ntlast =ntmax le(ntlast+1)=nelast+1 lr = lr3 if (lr3>=lr4) lr=lr4 if (lr<=kr4) then jk = kr4+1 do jj = lr,kr4 jk = jk-1 iv = vreg(jk) i = mreg(jk) nbelow = nbelow+1 if (iv<=nrow) nslck = nslck+1 ll = la(iv) kk = la(iv+1)-1 if (kk>ll.or.abs(a(ll)-1.0)>ztolze) then nueta = nueta+1 do j = ll,kk ir = ia(j) if (ir==i) then ep = dble(a(j)) else ie(nelast)=ir e(nelast) =a(j) nelast =nelast-1 nuelem =nuelem+1 endif end do ie(nelast) = i e(nelast) = ep le(ntlast) = nelast nelast = nelast-1 ntlast = ntlast-1 nuelem = nuelem+1 endif end do endif if (kr1/=0) then do j = lr1,kr1 ! Descomponer en forma LU en pico. iv = vreg(j) forall (i=1:nrow) y(i)=cero; forall (i=la(iv):la(iv+1)-1) y(ia(i))=a(i) ! Unpack en Y() call ftran (idos) ztoly = ztolpv 2070 continue ymax =cero irowp =0 ircmin=-999999999 do i = 1,nrow yval = abs(y(i)) if (yval>ztoly.and.hreg(i)<0) then ymax = max(yval,ymax) if (hreg(i)>ircmin) then ircmin=hreg(i) irowp =i endif endif end do if (irowp==0) then write (lo,8000) kinbas(iv) = 0 cycle endif ztmin = ztrel*ymax if (abs(y(irowp))toleta) then if (hreg(i)<0) then nelem =nelem+1 ! Elementos eta de L. ie(nelem)=i e(nelem) =y(i) else ie(nelast)=i ! Elementos eta de U. e(nelast) =y(i) nelast =nelast-1 nuelem =nuelem+1 endif endif end do jh(irowp) =iv kinbas(iv)=irowp nueta =nueta+1 ie(nelast)=irowp if (j==kr1) then e(nelast)=y(irowp) else e(nelast) =uno neta =neta+1 le(neta+1)=nelem+1 endif nuelem =nuelem+1 le(ntlast)=nelast nelast =nelast-1 ntlast =ntlast-1 do i=1,nrow ! Adaptar contadores de filas. if (abs(y(i))>toleta.and.hreg(i)<0) then hreg(i) = hreg(i)-incr if (hreg(i)>=0) hreg(i)=-1 endif end do hreg(irowp) = 0 end do endif ! Mezclar etas L y U. nleta =neta neta =nleta+nueta nlelem=nelem nelem =nlelem+nuelem if (nelem>=nemax) then write (lo,8005) stop endif if (nuelem>0) then nf =nemax-nuelem+1 ! Desplazar vectores IE y E de los elementos de U. incr=0 do i=nf,nemax incr =incr+1 iaux =nlelem+incr ie(iaux)=ie(i) e(iaux) =e(i) end do idif=nemax-nlelem-nuelem ! Redefinir vector LE de los elementos eta. nf =ntmax-nueta+1 incr=0 do i=nf,ntmax incr =incr+1 le(nleta+incr)=le(i)-idif end do le(neta+1) = nelem+1 endif do i=1,nrow ! Insertar holgura en columnas borradas. if (jh(i)/=0) cycle jh(i)=i irowp=i forall (j=1:nrow) y(j)=cero; forall (j=la(i):la(i+1)-1) y(ia(j))=a(j) ! Unpack en Y() call ftran (iuno) call wreta (irowp) end do ! ! -1 -1 -1 ! Adaptar x: x = B b - B Nx = B (b-Nx ) ! N N ! forall (i=1:nrow) y(i)=b(i) do j=1,ncol select case (kinbas(j)) case (-2); de=0 case (-1); de=xub(j) case (0); de=xlb(j) case default; cycle end select do i=la(j),la(j+1)-1 ir =ia(i) y(ir)=y(ir)-a(i)*de end do end do call ftran (iuno) forall (i=1:nrow) x(i)=y(i) nofd=nelem-neta ! Se imprimen estadísticas de inversión de la base. nstr=nrow-nslck if (inibb==0) write (lo,8500) nbnonz,nstr,nabove,nbelow,nlelem, & nleta,nuelem,nueta,nofd,neta ! ! -1 ! Calcular en y: b-B (B(b-Nx )-Nx ) y chequear error relativo ! N N ! forall (i=1:nrow) y(i)=b(i) x2max = cero do j=1,ncol select case (kinbas(j)) case (-2); de=0 case (-1); de=xub(j) case (0); de=xlb(j) case default; de=x(kinbas(j)) end select do i=la(j),la(j+1)-1 ir =ia(i) y(ir)=y(ir)-a(i)*de end do end do call ftran (iuno) eemax=maxval(abs(y(1:nrow))) eemax=max(eemax,cero) forall (i=1:nrow) x(i)=x(i)+y(i) x2max=maxval(abs(x(1:nrow))) x2max=max(x2max,cero) if (x2max/=cero) ermax=eemax/x2max if (inibb==0) write (lo,'(4X,''Error relativo en x:'',D14.6)') ermax if (ermax0.1) go to 8800 ztrel = 0.5 go to 10 8800 continue if (inibb==0) write (lo,8900) 8000 format('*** Matriz singular') 8005 format('*** Espacio para vectores ETA insuficiente; ', & 'el programa se para'//' incrementar E(.) y NEMAX') 8500 format(/'*** Estadísticas de INVERT'/4x,i7, & ' elementos no cero en la base'/4x,i7, & ' columnas estructurales en la base'/4x,i7, & ' vectores columna antes del "bump"'/4x,i7, & ' vectores columna después del "bump"'/4x,'L:',i7, & ' elementos no cero;',i7,' vectores ETA'/4x,'U:',i7, & ' elementos no cero;',i7,' vectores ETA'/4x,'Total:',i7, & ' elementos no en la diagonal;',i5,' vectores ETA') 8900 format(//' Ite. Inf. Nopt. Sum.Inf/F.Obj', & ' Entra de-a Cost.Red Sale de-a Paso ',' Pivote El.Eta') end subroutine invert subroutine normal ! ! Se coordinan todas las subrutinas necesarias para llevar a ! cabo la resolución de un programa lineal por el método ! simplex revisado. ! !******* Descripción de vectores afectados ***** ! ! B := Termino de la derecha del problema. ! ! KINBAS := Estado en el que se encuentran los componentes ! del vector X de variables de decisión: ! = 0, variable en su límite inferior; ! =-1, variable en su límite superior; ! =-2, variable libre: especificada FR; ! =>0, variable básica, el número indica ! el vector fila sobre el que pivota. ! ! XLB, XUB := Límites inferior y superior de las variables. ! ! JH := Vector que indica la columna con la que pivota ! cada vector fila. ! ! LA, IA, A := Vectores donde se guarda toda la información ! relativa a la matriz de los coeficientes ! de las condiciones. ! ! -1 -1 -1 ! X := vector B b - B Nx = B (b-Nx ) ! ! LE, IE, E := Vectores donde se guarda toda la información ! relativa a la matriz de las ! transformaciones elementales eta. ! ! YTEMP1, YTEMP2, YTEMP3 := vectores de trabajo. ! !****************************************************************** ! use iolog; use defpa; use nodbb; use parlp implicit none save integer :: jcolp,irowp,npivot,iptype,itent,i integer :: ivout,ivoute,ivines real :: sumobj,theta character :: ast*1,move1*4,move2*4 if (itsinv>=invfrq) then call invert itsinv = 0 endif itent=0 bucle_simplex: do ! ITERACIONES MÉTODO SIMPLEX call formc ! Asignación precios call btran call price (jcolp) ! Cálculo variable jcolp a entrar en base itent =itent+1 itcnt =itcnt+1 itsinv=itsinv+1 if (jcolp==0) return ! Se ha llegado al óptimo; precios óptimos forall (i=1:nrow) y(i)=cero; forall (i=la(jcolp):la(jcolp+1)-1) y(ia(i))=a(i) ! Unpack en Y() call ftran (1) call chuzr (jcolp,irowp,npivot,iptype,theta,ivout) ! Obtener variable irowp que sale de base if (inibb==0) then ! Escritura de resultados intermedios ast = ' ' if (npivot==1) then if (rcost>0.0) then move1 = 'L->B' else move1 = 'U->B' endif if (iptype==-1) then move2 = 'B->U' else move2 = 'B->L' endif else ast = '*' if (rcost>0.0) then move1 = 'L->U' move2 = 'L->U' else move1 = 'U->L' move2 = 'U->L' endif endif if (jcolp<=nrow) then ivines = jcolp+ncol-nrow-1 else ivines = jcolp-nrow endif if (ivout<=nrow) then ivoute = ivout+ncol-nrow-1 else ivoute = ivout-nrow endif sumobj = x(1)*minomax if (ninf>0) sumobj=suminf write (lo,8000) itcnt,ninf,nopt,sumobj,ivines,ast,move1, & rcost*minomax,ivoute,ast,move2,theta,rpivot,nelem endif if (npivot==1) then ! Adaptación estructuras de datos if (nelem>=nemax) then call invert itsinv = 0 cycle else call wreta (irowp) endif endif if (itsinv>=invfrq) then call invert itsinv = 0 cycle endif if (itcnt>=itrfrq) then write (lo,'("Límite de iteraciones excedido; el programa se para.")') if (inibb==0) stop call wrilpr (3) stop endif if (inibb==1.and.itent>=itpmax) then qstat='INFEAS' ! write (21,*) 'P-inf' return endif enddo bucle_simplex 8000 format(f8.0,2i5,e15.7,i5,a,1x,a,e9.2,i5,a,a,2e9.2,i8) end subroutine normal subroutine penlts(irowp,idir,iptype) ! ! * * ! Se calculan las penalizaciones P , P y P para cada variable ! U D G ! básica no entera debiendo serlo. También se comprueba la ! existencia de separaciones forzadas para variables no ! no básicas. Como variable de separación se escoge aquella ! con la mayor penalización. La siguiente acción se lleva ! a cabo en la dirección opuesta a la de máxima penalización. ! La de máxima penalización se añade a la lista de las ! que se analizarán más adelante. El añadido de ramas al ! árbol enumerativo se hace en la subrutina BRANCH llamada ! desde aquí. ! ! /--\ ! | | ! \--/ ! / \ ! / \ X >= Z + 1 (entero de x+1) ! / \ P P ! X <= Z /P P \ ! P P / D U \ ! / \ IDIR=1 ! / IDIR=-1 \ ! /--\ /--\ ! | | | | ! \--/ \--/ ! !****************************************************************** ! use defpa; use nodbb; use parlp implicit none integer :: i,j,idir,icol,ntemp2,irowp,iptype real :: ta0j,de,dp,pen qfix = ' ' do i = 2,nrow if (dfpart(i)<=ztolze) then pu(i)=cero; pd(i)=cero; pg(i)=cero else pu(i)=plinfy; pd(i)=plinfy; pg(i)=plinfy endif end do ! ! Se analizan todas las variables (COL=J=2,NCOL) no básicas: ! KINBAS(J)=-1 o KINBAS(J)=0. ! do j = 2,ncol if ((kinbas(j)==-1.or.kinbas(j)==0).and.abs(xub(j)-xlb(j))>ztolze) then ! ! Se calcula -1 ! B N ! J ! forall (i=1:nrow) y(i)=cero; forall (i=la(j):la(j+1)-1) y(ia(i))=a(i) ! Unpack en Y() call ftran (1) ! ! Se comprueban incrementos positivos para variables no básicas ! en su límite inferior y negativos para no básicas en su superior. ! if (kinbas(j)==-1) then forall (i=1:nrow) y(i)=-y(i) endif ! ! Para las variables no básicas que no han de tomar ! valores enteros hacemos el valor TA0K=0; ! para aquellas que sí, TA0j=Y(1). ! Y(1) = A (ver referencia). ! 00 if (j=nvare) then xival = x(1)-y(1) if (xival<=xincva+ztolze) then idir = 2*kinbas(j)+1 icol = j call branch (icol,idir) if (idir==-1) xlb(icol)=xub(j) if (idir== 1) xub(icol)=xlb(j) if (icol<=nrow) then icol = icol+ncol-nrow-1 else icol = icol-nrow endif if (qfix==' ') then write (qfix,'(A,I4,A)') "FIJ",icol," " else write (qfix,'(A,I4,A)') "FIJ",icol,"*" endif cycle endif endif ! do i = 2,nrow if (jh(i)>=nvare.and.dfpart(i)>ztolze) then ! ! Calcular la penalización P que determinara el acotar la ! U ! variable básica x(i), que toma un valor no entero, ! superiormente y consecuentemente una no básica x saldrá ! de su límite. j ! if (y(i)<-ztolpv) then de = max(y(1)*(dfpart(i)-1)/y(i),ta0j) pu(i) = min(de,pu(i)) endif ! ! Calcular la penalización P que determinara el acotar la ! D ! variable básica x(i), que toma un valor no entero, ! inferiormente y consecuentemente una no básica x saldrá ! de su límite. j ! if (y(i)>ztolpv) then de =max(y(1)*dfpart(i)/y(i),ta0j) pd(i)=min(de,pd(i)) endif ! ! Calcular la penalización P asociada a la introducción de un ! corte de Gomory. G ! ! Primero, analizar las variables no básicas x que no han ! de tomar valores enteros. j ! if (jztolpv) then de =dfpart(i)*y(1)/y(i) pg(i)=min(de,pg(i)) elseif (y(i)<-ztolpv) then de = (dfpart(i)-1)*y(1)/y(i) pg(i) = min(de,pg(i)) endif else ! ! Segundo, analizar las variables no básicas x que sí han de ! tomar valores enteros. j ! dp=abs(y(i))-floor(abs(y(i))) if (dp>ztolze.and.dp<1-ztolze) then if (y(i)<-ztolpv) dp=1-dp if (dp>dfpart(i)) then de=(1-dfpart(i))*y(1)/(1-dp) else de=dfpart(i)*y(1)/dp endif pg(i)=min(de,pg(i)) endif endif endif enddo endif end do ! ! Calcular la penalización P máxima y comprobar si se pueden ! G G ! abandonar las ramas que partirían del nudo en que nos ! encontramos y que resultarían de imponer nuevos límites a ! variables que debiendo ser enteras no cumplen este requisito ! de momento. ! pen = cero do i = 2,nrow if (jh(i)>=nvare) pen=max(pen,pg(i)) end do xival=x(1)-pen if (xival<=xincva+ztolze) then idir=0 return endif ! ! Comprobar si las penalizaciones P o P por imponer nuevos ! U D ! límites a la variable x(i) hacen que las soluciones ! continuas obtenibles son peores que la "titular" ! entera hasta este momento. ! ntemp2 = 0 do i = 2,nrow if (jh(i)>=nvare.and.x(1)-max(pd(i),pu(i))<=xincva+ztolze) then if (pu(i)<=pd(i)) then ! ! Asi es: se añade a la lista de nudos a analizar ! posteriormente el que reulta de bifurcar desde el nudo en ! que nos encontramos en la dirección resultante de acotar ! la variable x(i) inferiormente. Seguiremos inmediatamente ! analizando la rama resultante de acotar x(i) superiormente. ! xival =x(1)-pd(i) idir =-1 irowp =i icol =jh(irowp) ntemp2=1 call branch (icol,idir) xlb(icol)=ipart(i)+1 else ! ! Asi es: se añade a la lista de nudos a analizar ! posteriormente el que resulta de bifurcar desde el nudo en ! que nos encontramos en la dirección resultante de acotar ! la variable x(i) superiormente. Seguiremos inmediatamente ! analizando la rama resultante de acotar x(i) inferiormente. ! xival =x(1)-pu(i) idir =1 irowp =i icol =jh(irowp) ntemp2=1 call branch (icol,idir) xub(icol)=ipart(i) endif endif end do ! ! Las penalizaciones P o P por imponer nuevos límites a la ! U D ! variable x(i) no hacen presagiar que las soluciones ! continuas obtenibles son peores que la "titular" entera ! hasta este momento; se elige con respecto a qué variable ! bifurcar y en qué dirección. ! if (ntemp2==0) then pen =0.0 irowp=0 ! ! Se determina como variable de bifurcación aquella que tenga ! una mayor penalización P o P esperando así agotar ! D U ! pronto el árbol que parta de aquí. ! do i = 2,nrow if (jh(i)>=nvare) then if (pu(i)<=pd(i)) then if (pd(i)>pen) then pen =pd(i) irowp =i idir =-1 revbnd=dble(ipart(i)+1) endif elseif (pu(i)>pen) then pen =pu(i) irowp =i idir =1 revbnd=dble(ipart(i)) endif endif end do ! ! Problemas... Todas las penalizaciones P y P de las variables ! son cero: degeneración. D U ! Elegir cualquier variable básica no entera como variable ! de bifurcación; anadir a la lista de nudos a analizar ! la rama resultante de acotar superiormente la variable x(i) ! y analizar a continuación la derivada de acotar ! inferiormente la variable. ! if (irowp==0) then do irowp=2,nrow if (jh(irowp)>=nvare.and.dfpart(irowp)>ztolze) exit end do pen =pu(irowp) idir =1 revbnd=ipart(irowp) endif icol =jh(irowp) xival=x(1)-pen call branch (icol,idir) if (idir==-1) xlb(icol)=revbnd if (idir== 1) xub(icol)=revbnd endif if (idir==-1) iptype=0 if (idir== 1) iptype=-1 end subroutine penlts subroutine price(jcolp) ! ! Se asignan precios a las columnas no básicas y se determina ! la columna no básica JCOLP a entrar en la base. ! ! t -1 t t -1 t ! Z = c B b + (c -c B N)x ; en este caso c = 0. ! B N B N N ! ! El coste reducido de la columna no básica a es ! j ! t -1 ! -c B a . ! B j ! ! Se utiliza el criterio de Dantzig de coger la columna con un ! coste reducido que augura la mejora máxima de la función objetivo ! !******* Descripción de vectores afectados ***** ! ! KINBAS := Estado en el que se encuentran los componentes ! del vector X de variables de decisión: ! =0, variable en su límite inferior; ! =-1, variable en su límite superior; ! =-2, variable libre: especificada FR; ! =>0, variable básica, el número indica ! el vector fila sobre el que pivota. ! ! XLB, XUB := Límites inferior y superior de las variables. ! ! LA, IA, A := Vectores donde se guarda toda la información ! relativa a la matriz de los coeficientes ! de las condiciones. ! ! Y := El vector columna a . ! j !****************************************************************** ! use defpa; use nodbb; use parlp implicit none integer :: i,jcolp,j,jcol1,jcol2 real :: cmin,cmax,dsum cmin = plinfy cmax = -plinfy nopt = 0 ! ! Cálculo de los costes reducidos de todas las variables no ! básicas. ! ! t -1 ! r = c - c B a . ! j j B j ! ! en este caso c =0. ! j ! ! La función onjetivo se podra mejorar (este es un problema de ! maximización), si ! ! r > 0 para x =l ---> KINBAS(J)=0 ! j j j ! O ! r < 0 FOR x =u ---> KINBAS(J)=-1 ! j j j ! ! jcolp=0; jcol1=0; jcol2=0 do j = 2,ncol if (abs(xub(j)-xlb(j))>ztolze.and.kinbas(j)<=0) then dsum=cero do i=la(j),la(j+1)-1 dsum=dsum-a(i)*y(ia(i)) end do select case (kinbas(j)) case (-2) ! Variable libre; poner inmediatamente en la base jcolp=j rcost=dsum if (abs(dsum)>ztcost) nopt=nopt+1 case (0) if (dsum>ztcost) then nopt=nopt+1 if (dsum>cmax) then cmax = dsum; jcol1 = j ! if (inibb==1) then ! rcost=dsum; jcolp=j; return ! Implementaría regla de Bland en opt. entera ! endif endif endif case (-1) if (dsum<-ztcost) then nopt=nopt+1 if (dsum-ztcost) then jcolp=0; rcost=cero else jcolp=jcol2; rcost=cmin endif elseif (cmin<-ztcost) then if (cmax<=abs(cmin)) then jcolp=jcol2; rcost=cmin else jcolp=jcol1; rcost=cmax endif else jcolp=jcol1; rcost=cmax endif end subroutine price subroutine testx ! ! Se comprueba la solución óptima obtenida del programa lineal ! rechazándose el nudo no generándose, por tanto, a partir ! de el ninguna rama si: ! (1) el valor de la función objetivo de PL es peor que el ! de la mejor solución entera obtenida hasta este ! momento; ! (2) la solución PL es entera. !****************************************************************** ! use iolog; use defpa; use nodbb; use parlp implicit none integer :: i xival=x(1); noint=0 do i=2,nrow ! Parte entera y fraccionaria de cada variable básica. if (x(i)>-ztolze) then ipart(i)=floor(x(i)+ztolze) else ipart(i)=ceiling(x(i)+ztolze) endif ! dfpart(i) = dabs(x(i)-dble(ipart(i))) dfpart(i)=x(i)-ipart(i) if (jh(i)>=nvare.and.dfpart(i)>ztolze) noint=noint+1 end do if (noint/=0) return qstat = 'ENTERA' ! Si; si es mejor que la actual, guardarla if (inibb==0) then write (lo,'("*** La solución inicial del programa lineal es ENTERA ÓPTIMA ***")') return endif if (xival<=xincva+ztolze) return xincva=xival ! Nueva solución entera mejor que la actual do i=1,ncol kinben(i)=kinbas(i) select case (kinbas(i)) case (-1) xinc(i)=xub(i) if (abs(xub(i)-xlb(i))ztolpv) then nelem =nelem+1 ie(nelem)=i e(nelem) =y(i) endif end do neta =neta+1 le(neta+1)=nelem+1 end subroutine wreta subroutine wrilpr(ipa) ! ! Se escriben los resultados del problema. ! !****************************************************************** ! use iolog; use defpa; use nodbb; use parlp implicit none integer :: lode,i,j,ipa real :: ytempi,yi,bi,xubi,d1,act2,xtemp,redco real :: bact,d2 character*3 :: qll=" LI",qllb="LIB",qul=" LS",qulb="LSB",qeq=" EQ", & qeqb="EQB",qbas=" BS",qfx=" FX",qfr=" FR",qbsd="BSD",qbs character*14 :: char1,char2 select case (ipa) ! Se imprime función objetivo y estadísticas. case (1) if (qstat=='INFEAS') then write (lo,910) stop endif write (lo,10) lode = lo case (2) lode = logu write (lo,13) x(1)*minomax,cputim() if (logite==0) return write (logu,12) case (3) if (inibb==0) return lode = lo if (nsolen==0) then write (lo,16) return endif if (itcnt>=itrfrq) then write (lo,17) else write (lo,11) endif do i = 1,ncol if (kinben(i)== 0) xlb(i)=xinc(i) if (kinben(i)==-1) xub(i)=xinc(i) if (kinben(i)==-3) then xlb(i)=xinc(i) xub(i)=xinc(i) endif if (kinben(i)>0) x(kinben(i))=xinc(i) kinbas(i)=kinben(i) ! ¡OJO! con esto, sólo se hace al final de todo end do forall (i=1:nrow) y(i)=yen(i) end select write (lode,14) qnampo,int(itcnt),x(1)*minomax ! Resultados de condiciones char1 = ' Ninguno ' write (lode,20) 1,qname(1),qbas,x(1)*minomax,(-x(1)*minomax),char1,char1,1.0 do i = 2,nrow ytempi = 0.0 yi = y(i) bi = b(i) xubi = xub(i) if (a(i)>0) then d1 = bi-xubi if (xubi==plinfy) d1 = -plinfy d2 = bi if (kinbas(i)>0) then ! Básica ytempi=x(kinbas(i)) bact =bi-ytempi qbs =qbas if (ytempi0) qbs=qbsd bact=bi d1 =bi d2 =bi endif write (char1,'(G14.8)') d1 if (d1<=(-plinfy)+ztolze) char1 = ' Ninguno ' write (char2,'(G14.8)') d2 if (d2>=plinfy-ztolze) char2 = ' Ninguno ' write (lode,20) i,qname(i),qbs,bact,ytempi,char1,char2,yi*minomax end do write (lode,30) ! Resultados de las variables: COLUMNS. do i = nrow+1,ncol redco = cero select case (kinbas(i)) case (-2) ! Variable libre xtemp = cero qbs = qfr case (-1) ! Variable en su límite superior xtemp = xub(i) if (abs(xub(i)-xcs(i))= 1 xtemp=x(kinbas(i)) qbs =qbas end select if (abs(xub(i)-xlb(i))= plinfy-ztolze) char2=' Ninguno ' act2 = 0.0 do j = la(i),la(i+1)-1 if (ia(j)==1) act2=-a(j)*minomax end do write (lode,40) i-nrow,qname(i),qbs,xtemp,act2,char1,char2,redco*minomax end do if (logite==1.and.ipa/=3) write (logu,150) 10 format(/21x,'--- SOLUCIÓN INICIAL PROGRAMA LINEAL ---'/25x,32('-')/) 11 format(/26x,'--- SOLUCIÓN ENTERA OPTIMA ---'/30x,22('-')/) 17 format(/26x,'--- SOLUCIÓN ENTERA MEJOR ---'/30x,21('-')/) 16 format(/31x,'--- NO SE HA ENCONTRADO SOLUCIÓN ENTERA ---'/34x,35('-')/) 12 format(/34x,'--- SOLUCION ENTERA ---'/38x,15('-')) 13 format(/'* Nueva solución entera; z(PE)=',f17.5,'; Tiempo:',f11.3,' seg.') 14 format(4x,'Nombre del problema: ',a/4x,'No. de iteraciones: ',i11/ & 4x,'Valor de la función objetivo: ',g30.15///'*** FILAS'// & ' No. ..Fila.. en ....Valor.... ', & '...Holgura... .Lí.Inferior. ','.Lí.Superior. Val.Dual.'/) 20 format(i4,1x,a,a,2g14.8,2a14,g10.4) 30 format(//'*** COLUMNAS'//' No. .Columna en ','....Valor.... Coste en F.O. .Lí.',& 'Inferior. .Lí.Superior. Cos.Red.'/) 40 format(i4,1x,a,a,2g14.8,2a14,g9.3) 150 format(//' Variable',28x,'Nudos', & ' Variables Valor '/ & 'Separación Nivel Dirección en Lista ', & 'Ent. No Ent. Iteración Func. Obj.') 910 format('El Problema No es Factible') end subroutine wrilpr real function cputim() use iolog implicit none real*4 :: cp(2),etime cputim=etime(cp) end function cputim end PROGRAM Bbmi