Bug report (2/3)
Jean-Philippe Chancelier
Jean-Philippe.Chancelier@cereve.enpc.fr
Fri Sep 3 00:23:00 GMT 1999
integer inf,n,z
integer a(n,n),dm(n),dp(n),dualu(n),dualv(n),f(n),fb(n),pre(n),
+ r(n),uv(n)
integer d,i,i1,i2,i2p1,ii,iix,ind,index,indexv,iv,nuv,u,vgl,w
c ..
do 400 u = 1,n
if (f(u).gt.0) go to 400
do 50 i = 1,n
pre(i) = u
uv(i) = i
dp(i) = inf
dm(i) = a(u,i) - dualu(u) - dualv(i)
50 continue
nuv = n
i2 = 0
r(1) = u
dp(u) = 0
100 d = inf
index = 0
i1 = i2
do 150 iv = 1,nuv
i = uv(iv)
if (dm(i).le.d) then
if (dm(i).ne.d) then
index = 0
i2 = i1
end if
d = dm(i)
if (fb(i).le.0) then
if (d.eq.0) go to 250
index = i
end if
i2 = i2 + 1
r(i2) = iv
end if
150 continue
if (index.gt.0) go to 300
i1 = i1 + 1
i2p1 = i2 + i1
do 200 iix = i1,i2
ii = i2p1 - iix
indexv = r(ii)
index = uv(indexv)
uv(indexv) = uv(nuv)
nuv = nuv - 1
w = fb(index)
dp(w) = d
do 160 iv = 1,nuv
i = uv(iv)
vgl = d + a(w,i) - dualu(w) - dualv(i)
if (dm(i).gt.vgl) then
dm(i) = vgl
pre(i) = w
end if
160 continue
200 continue
go to 100
c
250 index = i
300 w = pre(index)
fb(index) = w
ind = f(w)
f(w) = index
if (w.eq.u) then
c
do 320 i = 1,n
if (dp(i).ne.inf) dualu(i) = dualu(i) + d - dp(i)
if (dm(i).lt.d) dualv(i) = dualv(i) + dm(i) - d
320 continue
else
index = ind
go to 300
end if
400 continue
c
z = 0
do 500 i = 1,n
z = z + dualu(i) + dualv(i)
500 continue
return
end
c *************************************************************************
subroutine agmhp(n,r,cr,ordsp,ic,c,f,v,ur,z,p,q,hp,dm,fb,zs,imp,
+ inf)
c agmheap search for an augmenting path starting at row r, using
c johnson algorithm, implemented via a heap-queue
integer imp,inf,n,ordsp,r,ur,z,zs
integer c(ordsp),cr(n+1),dm(n),f(n),fb(n),hp(n),ic(ordsp),p(n),
+ q(n),v(n)
integer delta,dj,dp1,du,fb1,fbj,h,hp1,hp2,hp3,i,ifin,ind,iniz,j,k,
+ k2,nhp,w
c ..
delta = zs - z + ur
imp = 0
f(r) = 0
do 100 i = 1,n
fb(i) = 0
100 continue
do 200 i = 1,n
q(i) = 0
p(i) = 0
j = f(i)
if (j.gt.0) fb(j) = i
dm(i) = inf
200 continue
nhp = 0
p(r) = 0
du = 0
w = r
iniz = cr(w)
ifin = cr(w+1) - 1
if (iniz.le.ifin) go to 500
go to 700
300 w = fb(i)
c
iniz = cr(w)
ifin = cr(w+1) - 1
if (iniz.gt.ifin) go to 700
c
do 400 h = iniz,ifin
j = ic(h)
if (j.eq.i) then
du = du - c(h)
go to 500
end if
400 continue
c
500 do 600 h = iniz,ifin
j = ic(h)
dj = du + c(h) - v(j)
if (w.ne.r) dj = dj + v(i)
if (dm(j).gt.dj) then
dm(j) = dj
p(j) = w
fbj = fb(j)
if (q(j).eq.0) then
c
nhp = nhp + 1
q(j) = nhp
end if
c
k = q(j)
520 k2 = k/2
if (k2.gt.0) then
hp2 = hp(k2)
if (dj.le.dm(hp2)) then
if (dj.ge.dm(hp2)) then
if (fbj.ne.0 .or. fb(hp2).eq.0) go to 540
end if
hp(k) = hp2
q(hp2) = k
k = k2
go to 520
end if
end if
540 hp(k) = j
q(j) = k
end if
600 continue
c
700 i = hp(1)
du = dm(i)
if (du.gt.delta) go to 1000
if (fb(i).eq.0) then
750 w = p(i)
fb(i) = w
ind = f(w)
f(w) = i
if (w.ne.r) then
c
i = ind
go to 750
else
c
do 760 i = 1,n
if (dm(i).lt.du) v(i) = v(i) + dm(i) - du
760 continue
z = z + du - ur
return
end if
else
q(i) = 0
nhp = nhp - 1
c
if (nhp.lt.0) then
c
imp = 1
return
else if (nhp.eq.0) then
go to 300
else
c
hp1 = hp(nhp+1)
dp1 = dm(hp1)
fb1 = fb(hp1)
k = 1
end if
end if
800 k2 = k*2
if (k2.lt.nhp) then
hp2 = hp(k2)
hp3 = hp(k2+1)
if (dm(hp2).ge.dm(hp3)) then
if (dm(hp2).ne.dm(hp3) .or. fb(hp2).ne.0) then
hp2 = hp3
k2 = k2 + 1
end if
end if
else if (k2.ne.nhp) then
go to 900
end if
hp2 = hp(k2)
if (dp1.le.dm(hp2)) then
if (dp1.lt.dm(hp2)) go to 900
if (fb1.eq.0 .or. fb(hp2).ne.0) go to 900
end if
hp(k) = hp2
q(hp2) = k
k = k2
go to 800
900 hp(k) = hp1
q(hp1) = k
go to 300
1000 z = inf
return
end
c *************************************************************************
subroutine enlini(n,cr,ic,f,u,v,s,p,flag,cont,a,cra,ica,nc,inf)
c try to connect subtours of the solution of the assignment at
c the root node, with zero reduced cost arcs
integer inf,n,nc
integer a(n,n),cont(n),cr(n+1),cra(*),f(n),flag(n),ic(*),ica(*),
+ p(n),s(n),u(n),v(n)
integer card,corid,flaga,fsalva,i,iall,ii,ij,ijn,ik,iolds,ip1,ip2,
+ istart,j,k,kc,kk,l,lm,lm2,m1,nc2,ninf2,pj,si,ult
c ..
ninf2 = -float(inf)/2.
m1 = -1
kc = 0
iall = 0
do 100 i = 1,n
j = f(i)
s(i) = j
p(j) = i
flag(i) = 0
cont(i) = 0
100 continue
do 200 istart = 1,n
if (flag(istart).ne.0) go to 200
kc = kc + 1
k = istart
card = 0
150 flag(k) = -kc
k = s(k)
if (v(k).gt.ninf2) card = card + 1
if (k.ne.istart) go to 150
if (card.ne.n) then
j = 2*kc
do 160 i = 1,j,2
if (cont(i).le.card) then
if (cont(i).ne.0) then
ip2 = i + 2
do 152 l = ip2,j
lm = j + ip2 - l
lm2 = lm - 2
cont(lm) = cont(lm2)
152 continue
end if
cont(i) = card
ip1 = i + 1
cont(ip1) = istart
go to 200
end if
160 continue
else
nc = 1
return
end if
200 continue
nc = kc
nc2 = nc*2
do 600 ii = 2,nc2,2
istart = cont(ii)
if (flag(istart).gt.0) go to 600
fsalva = iabs(flag(istart))
flaga = 0
i = istart
250 j = f(i)
if (v(j).le.ninf2) go to 550
if (cr(1).ne.m1) then
c
ip1 = i + 1
ult = cr(ip1) - 1
ik = i
k = cr(i)
j = ic(k)
go to 400
end if
300 k = 1
ult = n
350 j = k
c
400 if (i.eq.j) go to 500
450 corid = a(i,j) - u(i) - v(j)
if (corid.le.0) then
if (iabs(flag(j)).ne.fsalva) then
if (flaga.ne.0 .or. flag(j).le.0) then
pj = p(j)
if (v(j).ge.ninf2) then
si = f(i)
if (a(pj,si)-u(pj).le.v(si)) then
c
flaga = 1
iall = iall + 1
nc = nc - 1
si = f(i)
pj = p(j)
f(i) = j
p(j) = i
f(pj) = si
p(si) = pj
iolds = s(i)
kk = j
if (flag(j).gt.0) s(i) = pj
if (flag(j).le.0) s(i) = j
s(pj) = iolds
452 flag(kk) = fsalva
kk = f(kk)
if (kk.ne.si) go to 452
if (cr(1).eq.m1) go to 300
k = cr(i)
j = ic(k)
go to 400
end if
end if
end if
end if
end if
500 k = k + 1
if (cr(1).ne.m1) then
if (k.le.ult) then
j = ic(k)
go to 400
else if (cra(ik).ne.0) then
ij = cra(ik)
ik = ij
ijn = ij - n
j = ica(ijn)
go to 450
end if
else if (k.le.ult) then
go to 350
end if
550 flag(i) = iabs(flag(i))
i = s(i)
if (i.ne.istart) go to 250
600 continue
return
end
c *************************************************************************
subroutine enlarg(n,cr,ordsp,ic,c,f,u,vd,s,p,flag,cont,nc,inf)
c try to connect subtours of the current modified assignment
c problem, with zero reduced cost arcs
integer inf,n,nc,ordsp
integer c(ordsp),cont(n),cr(n+1),f(n),flag(n),ic(ordsp),p(n),s(n),
+ u(n),vd(n)
integer card,corid,flaga,fsalva,i,iall,ii,incr,iolds,ip1,ip2,ipj,
+ istart,j,jj,k,kc,kk,l,lm,lm2,m1,meta,minf2,nc2,pj,si,ult,
+ ult2
c ..
minf2 = -float(inf)/2.
kc = 0
m1 = -1
iall = 0
do 100 i = 1,n
j = f(i)
s(i) = j
p(j) = i
flag(i) = 0
cont(i) = 0
100 continue
do 200 istart = 1,n
if (flag(istart).ne.0) go to 200
kc = kc + 1
k = istart
card = 0
150 flag(k) = -kc
k = s(k)
if (vd(k).gt.minf2) card = card + 1
if (k.ne.istart) go to 150
if (card.ne.n) then
j = 2*kc
do 160 i = 1,j,2
if (cont(i).le.card) then
if (cont(i).ne.0) then
ip2 = i + 2
do 152 l = ip2,j
lm = j + ip2 - l
lm2 = lm - 2
cont(lm) = cont(lm2)
152 continue
end if
cont(i) = card
ip1 = i + 1
cont(ip1) = istart
go to 200
end if
160 continue
else
nc = 1
return
end if
200 continue
nc = kc
nc2 = nc*2
do 400 ii = 2,nc2,2
istart = cont(ii)
if (flag(istart).gt.0) go to 400
fsalva = iabs(flag(istart))
flaga = 0
i = istart
250 j = f(i)
if (vd(j).le.minf2) go to 350
ip1 = i + 1
ult = cr(ip1) - 1
k = cr(i)
300 j = ic(k)
corid = c(k) - u(i) - vd(j)
if (corid.le.0) then
if (iabs(flag(j)).ne.fsalva) then
if (flaga.ne.0 .or. flag(j).le.0) then
pj = p(j)
if (vd(j).ge.minf2) then
si = f(i)
ipj = pj + 1
ult2 = cr(ipj) - 1
kk = cr(pj)
meta = float(ic(ult2)-ic(kk))/2.
if (si.gt.meta) then
ipj = pj + 1
kk = cr(ipj) - 1
ult2 = cr(pj)
incr = -1
else
incr = 1
end if
302 jj = ic(kk)
if (jj.eq.si) then
if (c(kk)-u(pj).le.vd(si)) then
flaga = 1
iall = iall + 1
nc = nc - 1
si = f(i)
pj = p(j)
f(i) = j
p(j) = i
f(pj) = si
p(si) = pj
iolds = s(i)
kk = j
if (flag(j).gt.0) s(i) = pj
if (flag(j).le.0) s(i) = j
s(pj) = iolds
304 flag(kk) = fsalva
kk = f(kk)
if (kk.ne.si) go to 304
k = cr(i)
go to 300
end if
else if ((jj.le.si) .or. (incr.ne.1)) then
if ((jj.ge.si) .or. (incr.ne.m1)) then
if (kk.ne.ult2) then
kk = kk + incr
go to 302
end if
end if
end if
end if
end if
end if
end if
k = k + 1
if (k.le.ult) go to 300
350 flag(i) = iabs(flag(i))
i = s(i)
if (i.ne.istart) go to 250
400 continue
return
end
c *************************************************************************
subroutine karp(n,a,f,p,flag,na,nd,totass,ub,inf)
c heuristic solution of the atsp using the patching algorithm
integer inf,n,totass,ub
integer a(n,n),f(n),flag(n),na(n),nd(n),p(n)
integer cont,costo,fsalva,i,im1,istart,j,k,kk,kp,min,mn1,mn2,mn3,
+ nc,nim1,pmn1,rr,s2,temp,z
c ..
do 100 i = 1,n
flag(i) = 0
100 continue
cont = 0
nc = 0
do 200 i = 1,n
if (flag(i).le.0) then
k = i
nc = nc + 1
120 flag(k) = nc
cont = cont + 1
kk = f(k)
p(kk) = k
k = kk
if (k.ne.i) go to 120
na(nc) = cont
nd(nc) = i
cont = 0
end if
200 continue
if (nc.eq.1) return
c
c sort the subtours
c
ub = totass
z = nc
300 min = n
do 400 i = 1,z
if (na(i).lt.min) then
min = na(i)
k = i
end if
400 continue
na(k) = na(z)
na(z) = min
temp = nd(k)
nd(k) = nd(z)
nd(z) = temp
z = z - 1
if (z.ge.2) go to 300
rr = 2
500 istart = nd(rr)
i = istart
im1 = rr - 1
nim1 = nd(im1)
fsalva = flag(nim1)
mn2 = inf
600 do 700 k = 1,n
if (flag(k).eq.fsalva) then
j = f(i)
kp = p(k)
costo = a(i,k) + a(kp,j) - a(kp,k) - a(i,j)
if (costo.lt.mn2) then
mn1 = k
mn2 = costo
mn3 = i
if (costo.eq.0) go to 800
end if
end if
700 continue
i = f(i)
if (i.ne.istart) go to 600
800 s2 = f(mn3)
pmn1 = p(mn1)
f(mn3) = mn1
p(mn1) = mn3
f(pmn1) = s2
p(s2) = pmn1
ub = ub + mn2
c
c update flag
c
fsalva = flag(istart)
kk = mn1
900 flag(kk) = fsalva
kk = f(kk)
if (kk.ne.s2) go to 900
nc = nc - 1
if (nc.eq.1) return
rr = rr + 1
go to 500
end
c *************************************************************************
subroutine cercsb(f,flag,vd,n,cmin,isalva,inf)
c chose the subtour for the branch phase
integer cmin,inf,isalva,n
integer f(n),flag(n),vd(n)
integer card,i,istart,k,kc,minf2
c ..
minf2 = -float(inf)/2.
cmin = inf
kc = 0
do 100 i = 1,n
flag(i) = 0
100 continue
do 200 istart = 1,n
if (flag(istart).ne.0) go to 200
kc = kc + 1
card = 0
k = istart
150 flag(k) = kc
k = f(k)
if (vd(k).ge.minf2) card = card + 1
if (k.ne.istart) go to 150
if (card.lt.cmin) then
cmin = card
isalva = istart
end if
200 continue
return
end
c *************************************************************************
subroutine clearq(zstar,puntb,fmvf,ordx,x,ncodal,ordv,v,active)
c make not active nodes with lb greater than zstar
integer active,fmvf,ncodal,ordv,ordx,puntb,zstar
integer v(ordv),x(ordx)
integer j
c ..
100 if (ncodal.eq.0) return
if (v(puntb+1).lt.zstar) return
j = v(puntb+4)
x(j) = fmvf
fmvf = j
puntb = v(puntb)
ncodal = ncodal - 1
active = active - 1
go to 100
end
c *************************************************************************
subroutine inque(isalva,nnodi,sc1,psalvo,ngen,vd,f,n,totass,mv,mf,
+ mm2,vimpa,vimpb,ordv,v,ncodal,puntlv,punta,puntb,
+ ioffv,inf)
c psalvo = pointer to the father
c nnodi = number of nodes in the subtour
integer inf,ioffv,isalva,mm2,n,ncodal,ngen,nnodi,ordv,psalvo,
+ punta,puntb,puntlv,sc1,totass,vimpa,vimpb
integer f(n),mf(n),mv(n),v(ordv),vd(n)
integer i,i32,ksalva,lb,lb1,minf2,nm1,
+ nodob,nodopa,pcorr,plvp,puntol,sc2,sc3
c ..
i32 = 32000
c
c insert the problem in the queue in this order
c (from the third location):
c lb,
c psalvo,
c the generation number (among the sons of the same father) and the
c number of not imposed arcs in the subtour choosen for the branch
c phase (packed),
c the pointer to the column of matrices mv,mf,
c the additional arc to exclude from solution
c (start node + end node, packed),
c the score sc1 ((n.of excluded arcs)*10+number of imposed arcs),
c the not imposed arcs in the subtour choosen for the
c branch phase (packed)
c
minf2 = -float(inf)/2.
nm1 = ngen - 1
sc1 = sc1 + (nm1)*10 + 1
v(puntlv+2) = totass
v(puntlv+3) = psalvo
v(puntlv+4) = ngen*i32 + nnodi
v(puntlv+5) = mm2
do 100 i = 1,n
mv(i) = vd(i)
mf(i) = f(i)
100 continue
v(puntlv+6) = vimpa*i32 + vimpb
v(puntlv+7) = sc1
plvp = puntlv + ioffv
i = 1
ksalva = isalva
nodopa = isalva
200 nodob = f(nodopa)
if (vd(nodob).ge.minf2) then
v(plvp+i) = nodopa*i32 + nodob
i = i + 1
end if
nodopa = nodob
if (nodopa.ne.ksalva) go to 200
c
c puntlv: pointer at the first empty position of v()
c punta : pointer at the active problem with best lb
c puntb : pointer at the active problem with worst lb
c
ncodal = ncodal + 1
if (ncodal.eq.1) then
c
c insert the first problem
c
punta = puntlv
puntb = puntlv + 1
else
sc2 = v(punta+ioffv)
sc3 = v(puntb+ioffv-1)
lb = totass
if (lb.lt.v(punta+2) .or. (lb.eq.v(punta+2).and.(sc1.ge.sc2)))
+ then
c
v(puntlv) = punta
v(punta+1) = puntlv + 1
punta = puntlv
else if (lb.gt.v(puntb+1) .or. (lb.eq.v(puntb+1).and.
+ (sc1.le.sc3))) then
c
v(puntlv+1) = puntb
v(puntb-1) = puntlv
puntb = puntlv + 1
else
c
c find the first active problem with cost greater than lb
c
pcorr = v(punta)
220 sc2 = v(pcorr+ioffv)
lb1 = v(pcorr+2)
if ((lb.lt.lb1) .or. ((lb.eq.lb1).and. (sc1.ge.sc2))) then
c
c insert the problem in queue
c
puntol = v(pcorr+1) - 1
v(puntol) = puntlv
v(puntlv) = pcorr
v(pcorr+1) = puntlv + 1
v(puntlv+1) = puntol + 1
else
pcorr = v(pcorr)
go to 220
end if
end if
end if
puntlv = puntlv + ioffv + nnodi + 1
return
end
c *************************************************************************
subroutine exque(ordv,v,ordx,x,fmvf,vd,f,p2,ngen,punta,nprob,
+ iva,ivb,ncodal,n,sc1,totass,ioffv)
integer fmvf,ioffv,n,ncodal,ngen,nprob,ordv,ordx,p2,punta,sc1,
+ totass
integer f(n),iva(n),ivb(n),v(ordv),vd(n),x(ordx)
c .. local scalars r32,i,i32,ix,j,jp
c ..
i32 = 32000
r32 = 32000.
p2 = punta
ngen = float(v(punta+4))/r32
nprob = v(punta+4) - ngen*i32
totass = v(punta+2)
jp = v(punta+5)
j = 1
do 100 i = 1,nprob
ix = float(v(punta+ioffv+i))/r32
iva(j) = ix
ivb(j) = v(punta+ioffv+i) - ix*i32
j = j + 1
100 continue
c
c load the new vectors vd and f
c
do 200 i = 0,n - 1
vd(i+1) = x(jp+i)
f(i+1) = x(jp+n+i)
200 continue
c
x(jp) = fmvf
fmvf = jp
sc1 = v(punta+ioffv)
punta = v(punta)
ncodal = ncodal - 1
return
end
c *************************************************************************
subroutine actpro(patt,n,nprobv,pnuovo,ngen,vimpa,vimpb,ordv,v,cr,
+ ordsp,ic,c,ioffv,inf)
integer ordsp,ordv,patt,pnuovo,vimpa,vimpb
integer c(ordsp),cr(n+1),ic(ordsp),v(ordv)
integer p1,pold,pp
c subroutines ..modmat
m1 = -1
p1 = 1
i32 = 32000
r32 = 32000.
c update the cost matrix from the actual one(associated with patt)
c to the new one(associated with pnuovo)
if (v(pnuovo+3).ne.patt) then
c the problem patt is not the father of pnuovo
c
c "mark" the street from pnuovo to the root node
pp = pnuovo
50 if (pp.eq.0) then
ngenol = nprobv
c find a marked problem
60 nodoa = float(v(patt + ioffv + ngenol))/r32
nodob = v(patt + ioffv + ngenol) - nodoa*i32
nodoc = vimpa
nodod = vimpb
c update the matrix from the one associated with the son of patt
c to the one associated with patt
call modmat(cr,ordsp,ic,c,nodoa,nodob,m1,inf)
if (nodoc.ne.0) call modmat(cr,ordsp,ic,c,nodoc,nodod,m1,
+ inf)
if (v(patt+5).le.0) then
c the actual matrix (associated with the problem patt)
c is the one associated with common predecessor of pnew
c and the first patt
c
c pold : common predecessor problem
pp = pnuovo
pold = patt
c
70 patt = v(pp+3)
ngenpp = float(v(pp+4))/r32
c prohibit arcs
nodoa = float(v(patt + ioffv + ngenpp))/r32
nodob = v(patt + ioffv + ngenpp) - nodoa*i32
nodoc = float(v(pp+6))/r32
nodod = v(pp+6) - nodoc*i32
call modmat(cr,ordsp,ic,c,nodoa,nodob,p1,inf)
if (nodoc.ne.0) call modmat(cr,ordsp,ic,c,nodoc,nodod,
+ p1,inf)
if (patt.eq.pold) then
c
c update the marked problems
c
pp = pnuovo
75 v(pp+5) = -v(pp+5)
pp = v(pp+3)
if (pp.ne.0) go to 75
else
pp = patt
go to 70
end if
else
ngenol = float(v(patt+4))/r32
vimpa = float(v(patt+6))/r32
vimpb = v(patt+6) - vimpa*i32
patt = v(patt+3)
go to 60
end if
else
c
v(pp+5) = -v(pp+5)
pp = v(pp+3)
go to 50
end if
else
c
c the problem patt is the father of pnuovo
c
if (ngen.eq.nprobv) return
nodoa = float(v(patt + ioffv + nprobv))/r32
nodob = v(patt + ioffv + nprobv) - nodoa*i32
c
c update the cost of the prohibited arc
c
call modmat(cr,ordsp,ic,c,nodoa,nodob,m1,inf)
nodoc = vimpa
nodod = vimpb
c update the cost of the ''implicit'' prohibited arc
if (nodoc.ne.0) call modmat(cr,ordsp,ic,c,nodoc,nodod,m1,
+ inf)
nodoa = float(v(patt + ioffv + ngen))/r32
nodob = v(patt + ioffv + ngen) - nodoa*i32
c prohibit the arc nodoa,nodob
call modmat(cr,ordsp,ic,c,nodoa,nodob,p1,inf)
nodoc = float(v(pnuovo+6))/r32
nodod = v(pnuovo+6) - nodoc*i32
c prohibit the ''implicit''arc
if (nodoc.ne.0) call modmat(cr,ordsp,ic,c,nodoc,nodod,p1,
+ inf)
return
end if
return
end
c *************************************************************************
subroutine genson(n,iva,ivb,vd,nprob,i,vimpa,vimpb,nodoim,cr,
+ ordsp,ic,c,inf)
integer i,inf,n,nodoim,nprob,ordsp,vimpa,vimpb
integer c(ordsp),cr(n+1),ic(ordsp),iva(n),ivb(n),vd(n)
integer i2,ibm1,im1,m1,p1
c subroutines modmat
m1 = -1
p1 = 1
vimpa = 0
vimpb = 0
c prohibit i-th arc
call modmat(cr,ordsp,ic,c,iva(i),ivb(i),p1,inf)
if (i.ne.1) then
im1 = i - 1
ibm1 = ivb(im1)
i2 = ibm1
call modmat(cr,ordsp,ic,c,iva(im1),i2,m1,inf)
vd(ibm1) = vd(ibm1) - inf
c don't prohibit 2 times the same arc
if (nodoim.ne.ivb(i)) then
i2 = nodoim
call modmat(cr,ordsp,ic,c,iva(i),i2,p1,inf)
vimpa = iva(i)
vimpb = nodoim
else
vimpa = 0
vimpb = 0
end if
else
nodoim = ivb(nprob)
return
end if
c
c if necessary adjust the 'implicit' prohibited arc
c
i2 = nodoim
if (i.ge.3) call modmat(cr,ordsp,ic,c,iva(im1),i2,m1,inf)
return
end
c *************************************************************************
subroutine calcur(n,cr,ordsp,ic,c,f,ivai,vd,ur,inf)
c compute the value of the dual variable associated with row ivai
integer inf,ivai,n,ordsp,ur
integer c(ordsp),cr(n+1),f(n),ic(ordsp),vd(n)
integer jk,kl,pr,ul,vr
c ..
jk = f(ivai)
vr = vd(jk)
pr = cr(ivai)
ul = cr(ivai+1) - 1
do 100 kl = pr,ul
if (ic(kl).eq.jk) then
ur = c(kl) - vr - inf
go to 200
end if
100 continue
200 return
end
c *************************************************************************
subroutine calcud(n,ordsp,ic,cr,c,vc,fc,uc)
c compute the current dual variables u
integer n,ordsp
integer c(ordsp),cr(n+1),fc(n),ic(ordsp),uc(n),vc(n)
integer i,ifs,ii,ils,jk
c ..
do 100 i = 1,n
jk = fc(i)
ifs = cr(i)
ils = cr(i+1) - 1
do 50 ii = ifs,ils
if (ic(ii).eq.jk) uc(i) = c(ii) - vc(jk)
50 continue
100 continue
return
end
c *************************************************************************
subroutine modmat(cr,ordsp,ic,c,nodoa,nodob,flag,inf)
c if flag=1 exclude arc (nodoa,nodob)
c if flag=-1 remove exclusion constraint form arc (nodoa,nodob)
integer flag,inf,nodoa,nodob,ordsp
integer c(ordsp),cr(*),ic(ordsp)
integer j,k,pr,ul
c ..
pr = cr(nodoa)
ul = cr(nodoa+1) - 1
do 100 j = pr,ul
k = ic(j)
if (k.eq.nodob) then
c(j) = c(j) + flag*inf
go to 200
end if
100 continue
200 return
end
c *************************************************************************
subroutine creams(gap,a,n,cr,x,mm6,mm4,mm3,u,v,zeri,nlsp,err)
c reduction procedure
c store the sparse matrix in x() and defines mm3,mm4
c vector ic() is stored from mm3 and vector c() from mm4
integer err,gap,mm3,mm4,mm6,n,nlsp,zeri
integer a(n,n),cr(n+1),u(n),v(n),x(*)
integer gap2,i,ia,iu,j,k,middle
c ..
zeri = 0
err = 0
c
mm4 = mm6
middle = (mm4+mm3)/2.
k = middle
do 100 i = n,1,-1
cr(i+1) = mm4
iu = u(i)
gap2 = gap + iu
do 50 j = n,1,-1
if (i.ne.j) then
ia = a(i,j) - v(j)
a(i,j) = ia - iu
if (ia.le.gap2) then
if (ia.eq.iu) zeri = zeri + 1
mm4 = mm4 - 1
if (mm4.le.middle) go to 500
x(k) = j
x(mm4) = ia - u(i)
k = k - 1
end if
end if
50 continue
100 continue
More information about the Gcc-bugs
mailing list