This is the mail archive of the
gcc-bugs@gcc.gnu.org
mailing list for the GCC project.
Bug report (1/3)
- To: gcc-bugs at gcc dot gnu dot org
- Subject: Bug report (1/3)
- From: Jean-Philippe Chancelier <Jean-Philippe dot Chancelier at cereve dot enpc dot fr>
- Date: Fri, 3 Sep 1999 09:23:10 +0200
To: gcc-bugs@gcc.gnu.org
Subject: Bug report
FCC: ~/mail/outmail
Mime-Version: 1.0 (generated by tm-edit 7.106)
Content-Type: text/plain; charset=US-ASCII
Bug report for
--------------------------------------------
g77 version 2.95.1 19990816 (release) (from FSF-g77 version 0.5.25 19990816 (release))
On
deedee% uname -a
SunOS deedee 5.6 Generic sun4u sparc SUNW,Ultra-1
--------------------------------------------
deedee% g77 -v --save-temps -O visitor.f
g77 version 2.95.1 19990816 (release) (from FSF-g77 version 0.5.25 19990816 (release))
Driving: g77 -v -save-temps -O visitor.f -lg2c -lm
Reading specs from /usr/local/lib/gcc-lib/sparc-sun-solaris2.6/2.95.1/specs
gcc version 2.95.1 19990816 (release)
/usr/local/lib/gcc-lib/sparc-sun-solaris2.6/2.95.1/f771 visitor.f -quiet -dumpbase visitor.f -O -version -fversion -o visitor.s
GNU F77 version 2.95.1 19990816 (release) (sparc-sun-solaris2.6) compiled by GNU C version 2.95.1 19990816 (release).
GNU Fortran Front End version 0.5.25 19990816 (release)
visitor.f: In subroutine `setupo':
visitor.f:2066: Internal compiler error in `final_scan_insn', at final.c:2920
Please submit a full bug report.
See <URL:http://www.gnu.org/software/gcc/faq.html#bugreport> for instructions.
deedee%
--------------------------------------------
--save-temps only save a .s for fortran files
but the Fortran file is self contained so I add it to my mail
-----------------------------visitor.f-----------------------------------
c algorithm 750, collected algorithms from acm.
c this work published in transactions on mathematical software,
c vol. 21, no. 4, december, 1995, p. 410--415.
c *************************************************************************
subroutine visitor(n,a,nstac,iroad,x,fstar)
integer a(n,n),x(nstac)
c real alpha,avson,spars
c integer active,err,i,inf,lb0,lbc,lopt,maxnd,n,nass,nexp,nprobq,
c + ordx,zeur,zstar
integer fstar(n),iroad(n),active,err,ordx,zeur,zstar
do 11,i=1,n
iroad(i)=0
11 continue
ordx = nstac
inf = 99999999
alpha = -1.
maxnd = -1
zeur = -1
do 1 i=1,nstac
x(i)=0.
1 continue
do 2 i=1,n
do 2 j=1,n
k=i+(j-1)*n
x(k)=a(i,j)
2 continue
c
call cdt(n,ordx,x,maxnd,inf,alpha,zeur,zstar,fstar,lb0,lbc,nexp,
+ nprobq,nass,active,lopt,spars,avson,err)
c
if (err.eq.0) go to 10
c wwfmt='('' solution not optimal '')'
10 continue
iroad(1)=1
do 20 i=2,n
iroad(i)=fstar(iroad(i-1))
20 continue
return
end
c *************************************************************************
subroutine cdt(n,ordx,x,maxnd,inf,alpha,zeur,zstar,fstar,lb0,lbc,
+ nexp,nprobq,nass,active,lopt,spars,avson,err)
c i/o parameter --------------------------------------------------------
c input parameters
c alpha = if .gt. 0 define the artificial upper bound ub=alpha*lb0
c if .le. 0 no effect
c n = number of nodes
c inf = very large number
c maxnd = if greater or equal to zero it is the maximum number of
c node to be explored on the branch decision-tree
c ordx = dimension of vector x
c zeur = if greater then zero it must define a valid upper bound
c the cost matrix is stored (column by column) in the first n**2 elem.
c of the vector x. the cost matrix it is not available in output.
c
c output value
c
c active = number of active problems in queue when the program stops
c (not yet and not currently examined)
c avson = average number of son nodes for the explored problems
c err = .gt. zero if an error occurred
c = -1 if no solution with value less than the artificial
c upper bound was found.
c in this case zstar is the artificial upper bound value
c fstar(i)= successor of node i in the current best solution
c lb0 = value of the optimal assignment at the root node
c lbc = value of the current best lb, i.e. value of the lower
c bound of the last problem extracted from the queue;
c lbc is a valid lower bound for the instance
c nass = number of assignment solved
c nexp = number of explored problems
c nprobq = number of problems stored in the queue
c spars = sparsity of the reduced matrix
c zeur = minimum value between the value of zeur in input and the
c value of the approximate solution computed using the
c patching heuristic
c zstar = value of the current best solution
c
c internal storage of vectors and matrices in vector x
c
c x(mm1) = the queue v(ordv) (dinamically allocated)
c x(mm2) = first element of matrices mv,mf (dinamically allocated)
c x(mm3) = vector ic(*)
c x(mm4) = vectors c(*) (cra(*) in ctcs,
c x(mm5) = ica(*) in ctcs)
c x(mm6) = iva(n)
c x(mm7) = ivb(n)
c x(mm8) = cr(n+1)
c x(mm9) = primal solution of the map associated with the current
c branched node f(n)
c x(mm10) = dual var. of the columns in the solution of ap, dualv(n)
c x(mm11) = primal solution of the map at the current node fc(n)
c x(mm12) = current dual variables of map, vd(n)
c x(mm13) = the heuristic solution of atsp
c x(mm14)..x(mm26) = temporary vectors
c
c fmvf = pointer to the first free element to store columns of
c matrices mv,mf
c x(fmvf) = pointer to the next free element for matrices mv,mf
c
c local variables
c .. scalar arguments ..
c real alpha,avson,spars
c integer active,err,inf,lb0,lbc,lopt,maxnd,n,nass,nexp,nprobq,ordx,
c + zeur,zstar
integer active,err,ordx,zeur,zstar
c .. array arguments ..
integer fstar(n),x(ordx)
c .. local scalars ..
c integer fmvf,i,ibig,idelta,ienlrg,imp,inf2,ipp,isalva,
c + isalvi,ivai,k,maxdix,maxica,mm1,mm10,mm11,mm12,mm13,mm14,
c + mm15,mm16,mm17,mm18,mm19,mm2,mm20,mm21,mm22,mm23,mm24,
c + mm25,mm26,mm3,mm4,mm5,mm6,mm7,mm8,mm9,nc,nci,ncodal,ngen,
c + nlsp,nnodi,nnodin,nodoim,nprob,nprobv,nq,offv,orcr,orcra,
c + ordsp,ordv,pnuovo,psalvo,punta,puntb,puntlv,r,sc1,totass,
c + ur,vimpa,vimpb,z,zarf,zc,zeri
integer fmvf,offv,orcr,orcra,ordsp,ordv,pnuovo,psalvo,punta,
+puntb,puntlv,r,sc1,totass,ur,vimpa,vimpb,z,zarf,zc,zeri
logical artif,eur
c subroutines ..
c actpro,agmhp,calcud,calcur,cercsb,clearq,contci,copyx,creams,ctcs,
c enlarg,enlini,errors,exque,genson,inque,inquer,karp
c ..
inf2 = float(inf)/2.
nq = n**2
err = 0
c
mm26 = ordx - n + 1
mm25 = mm26 - n
mm24 = mm25 - n
mm23 = mm24 - n
mm22 = mm23 - n
mm21 = mm22 - n
mm20 = mm21 - n
mm19 = mm20 - n
mm18 = mm19 - n
mm17 = mm18 - n
mm16 = mm17 - n
mm15 = mm16 - n
mm14 = mm15 - n
mm13 = mm14 - n
mm12 = mm13 - n
mm11 = mm12 - n
mm10 = mm11 - n
mm9 = mm10 - n
mm8 = mm9 - n - 1
mm7 = mm8 - n
mm6 = mm7 - n
if (mm6.gt.nq) then
mm5 = nq + (mm6-nq)*2./3.
mm4 = nq + (mm6-nq)/3.
mm3 = nq + 1
maxica = mm6 - mm5 + 1
orcra = mm5 - mm4 + 1
maxdix = mm4 - mm3 + 1
orcr = n + 1
c
nprob = 1
psalvo = 0
ncodal = 0
puntlv = 1
vimpa = 0
vimpb = 0
sc1 = -1
ngen = 1
fmvf = 0
nass = 1
c statistics information
nexp = 1
nprobq = 1
lopt = 1
avson = 0.
spars = 0.
active = 0
c solve the initial assignment problem
call ctcs(n,z,1,x(1),x(mm9),x(mm13),x(mm14),x(mm15),x(mm16),
+ x(mm17),x(mm4),x(mm5),x(mm8),x(mm3),x(mm26),x(mm10),
+ x(mm20),x(mm21),x(mm23),x(mm22),x(mm24),maxdix,
+ maxica,x(mm25),orcr,orcra,inf2,err)
if (err.eq.0) then
lb0 = z
lbc = z
c
c totass stores the value of the ap solution at the root node
totass = z
call contci(x(mm9),nci,n,x(mm14))
call enlini(n,x(mm8),x(mm3),x(mm9),x(mm26),x(mm10),
+ x(mm14),x(mm15),x(mm16),x(mm17),x(1),x(mm4),
+ x(mm5),nc,inf2)
zstar = totass
if (nc.eq.1) then
c optimal solution found at the root node by ctcs or enlini
call copyx(x(mm9),fstar,n)
zstar = totass
return
else
c
call copyx(x(mm9),x(mm13),n)
call karp(n,x(1),x(mm13),x(mm15),x(mm16),x(mm17),
+ x(mm18),totass,zstar,inf2)
c check if the heuristic solution is an optimum
if (zeur.gt.0 .and. zeur.lt.zstar) then
eur = .true.
zstar = zeur
else
eur = .false.
call copyx(x(mm13),fstar,n)
end if
if (zstar.eq.totass) return
avson = nc
artif = .false.
if (alpha.gt.0.0) then
if (totass*alpha.lt.zstar) then
zstar = totass*alpha
zarf = zstar
artif = .true.
end if
end if
c
idelta = zstar - lb0
call creams(idelta,x(1),n,x(mm8),x(1),mm6,mm4,mm3,
+ x(mm26),x(mm10),zeri,nlsp,err)
if (err.eq.0) then
ordsp = nlsp
spars = nlsp*100./float(n* (n-1))
ienlrg = 0
if (float(zeri).gt.2.5*float(n)) ienlrg = 1
c
c choose the subtour for the branch phase
c
call cercsb(x(mm9),x(mm14),x(mm10),n,nnodi,isalva,
+ inf2)
nnodin = nnodi
isalvi = isalva
c
c define the working arrays to store the queue:
c scalar and vectorial information
mm1 = 1
mm2 = mm3
ordv = mm2 - mm1
c
c mm1 = first element containing scalar informations
c mm2 = first position occupied by matrices mf,mv
c
c insert the root node in queue and initialize the vectors
c for the branching phase
puntlv = 1
call inquer(isalvi,nnodin,psalvo,ngen,n,totass,
+ vimpa,vimpb,x(mm9),x(mm10),puntlv,
+ punta,puntb,ordv,x(mm6),x(mm7),x(mm1),
+ pnuovo,nprob,sc1,offv,err)
if (err.eq.0) then
psalvo = pnuovo
active = 0
c ww(6,fmt='('' root node: zstar='',i10,'' lb0='',i10)') zstar,lb0
else
call errors(err,4)
return
end if
else
call errors(err,3)
return
end if
end if
else
call errors(err,2)
return
end if
else
call errors(err,1)
return
end if
c psalvo is the pointer to the actual problem
100 do 200 i = 1,nprob
c
c generate the descending nodes of psalvo
c
call genson(n,x(mm6),x(mm7),x(mm10),nprob,i,vimpa,vimpb,
+ nodoim,x(mm8),ordsp,x(mm3),x(mm4),inf2)
call copyx(x(mm9),x(mm11),n)
call copyx(x(mm10),x(mm12),n)
zc = totass
ivai = x(mm6+i-1)
c
c solve the new map problem
c
call calcur(n,x(mm8),ordsp,x(mm3),x(mm4),x(mm11),ivai,x(mm12),
+ ur,inf2)
call agmhp(n,ivai,x(mm8),ordsp,x(mm3),x(mm4),x(mm11),x(mm12),
+ ur,zc,x(mm14),x(mm15),x(mm16),x(mm17),x(mm18),
+ zstar,imp,inf2)
if (nass.eq.maxnd) then
call errors(err,-1)
return
else if (imp.ne.1) then
if (zc.lt.zstar) then
nass = nass + 1
c
c count the subtours
c
call contci(x(mm11),nc,n,x(mm15))
if (ienlrg.ne.0) then
call calcud(n,ordsp,x(mm3),x(mm8),x(mm4),x(mm12),
+ x(mm11),x(mm14))
call enlarg(n,x(mm8),ordsp,x(mm3),x(mm4),x(mm11),
+ x(mm14),x(mm12),x(mm15),x(mm16),
+ x(mm17),x(mm18),nc,inf2)
end if
if (nc.eq.1) then
c
c found a new feasible solution
c
zstar = zc
c
c clear the queue of the active problems
c
call clearq(zstar,puntb,fmvf,ordx,x(1),ncodal,
+ ordv,x(mm1),active)
call copyx(x(mm11),fstar,n)
c
c wwfmt='('' n.exp='',i8,'' n.prob.q='',i8,'' n.ass='',i8, '' active='
c ',i8,'' av.son='',f8.2)' nexp,nprobq,nass,active,avson
c
c wwfmt='('' zstar='',i10,'' lbc='',i10,'' nprobq='',i8,'' active='',i8)'
c zstar,lbc,nprobq,active
c statistics
lopt = 2
ipp = psalvo
105 ipp = x(mm1+ipp+2)
lopt = lopt + 1
if (ipp.gt.1) go to 105
c
c if the cost zc of the solution of the current problem p is ugual to
c the lb of the problem father of p, say pf, no other son of pf
c have to be generated
c
ibig = x(mm1+psalvo+1)
if (zc.le.ibig) then
nprobv = i
go to 300
end if
else
call cercsb(x(mm11),x(mm14),x(mm12),n,nnodi,
+ isalva,inf2)
ngen = i
c
c insert the problem in queue
c
if (fmvf.ne.0) then
k = fmvf
fmvf = x(fmvf)
else
mm2 = mm2 - 2*n
k = mm2
end if
ordv = mm2 - mm1
if (mm1+puntlv+offv+nnodi.lt.mm2) then
call inque(isalva,nnodi,sc1,psalvo,ngen,
+ x(mm12),x(mm11),n,zc,x(k),x(k+n),k,
+ vimpa,vimpb,ordv,x(mm1),ncodal,
+ puntlv,punta,puntb,offv,inf2)
nprobq = nprobq + 1
active = active + 1
r = int(nprobq/1000.)*1000
if (r.eq.nprobq) continue
c wwfmt='('' zstar='',i10,'' lbc='',i10,'' nprobq='',i8, '' active='
c ',i8)' zstar,lbc,nprobq,active
else
call errors(err,5)
return
end if
end if
end if
end if
c
200 continue
c
nprobv = nprob
300 if (ncodal.ne.0) then
c
c extract the problem with lowest bound
call exque(ordv,x(mm1),ordx,x(1),fmvf,x(mm10),x(mm9),
+ pnuovo,ngen,punta,nprob,x(mm6),x(mm7),ncodal,n,sc1,
+ totass,offv)
avson = (avson*nexp+nprob)/float(nexp+1)
nexp = nexp + 1
active = active - 1
if (lbc.lt.totass) then
lbc = totass
c wwfmt='('' zstar='',i10,'' lbc='',i10,'' nprobq='',
c i8,'' active='',i8)' zstar,lbc,nprobq,active
end if
lbc = totass
c
c update sets of included and excluded arcs
c
call actpro(psalvo,n,nprobv,pnuovo,ngen,vimpa,vimpb,ordv,
+ x(mm1),x(mm8),ordsp,x(mm3),x(mm4),offv,inf2)
psalvo = pnuovo
go to 100
end if
c
c
c the tree search is terminated
c
if (artif .and. zstar.eq.zarf) call errors(err,6)
if (eur) then
c wwfmt='('' zeur is the optimal solution value'')'
c wwfmt='('' fstar() do not contain the optimal tour'')'
end if
return
end
c *************************************************************************
subroutine contci(f,nc,n,flag)
integer f(n),flag(n)
c .. local scalars i,k
nc = 0
do 100 i = 1,n
flag(i) = 0
100 continue
do 200 i = 1,n
if (flag(i).le.0) then
k = i
nc = nc + 1
120 flag(k) = 1
k = f(k)
if (k.ne.i) go to 120
end if
200 continue
return
end
c *************************************************************************
subroutine loadfv(f,v,n,vroot,froot)
integer f(n),froot(n),v(n),vroot(n)
do 100 i = 1,n
froot(i) = f(i)
vroot(i) = v(i)
100 continue
return
end
c *************************************************************************
subroutine backfv(vroot,froot,n,vd,f)
integer f(n),froot(n),vd(n),vroot(n)
do 100 i = 1,n
f(i) = froot(i)
vd(i) = vroot(i)
100 continue
return
end
c *************************************************************************
subroutine inquer(isalva,nnodi,psalvo,ngen,n,totass,vimpa,vimpb,f,
+ vd,puntlv,punta,puntb,ordv,iva,ivb,v,p2p,nprob,
+ sc1,offv,ierr)
integer offv,ordv,p2p,psalvo,punta,
+ puntb,puntlv,sc1,totass,vimpa,vimpb
integer f(n),iva(n),ivb(n),v(ordv),vd(n)
integer i,i32,ksalva,nm1,nodob,nodopa,plvp
c ..
i32 = 32000
offv = 7
if (puntlv+offv+nnodi.le.ordv) then
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) = 1
v(puntlv+6) = vimpa*i32 + vimpb
v(puntlv+7) = sc1
plvp = puntlv + offv
i = 1
ksalva = isalva
nodopa = isalva
50 nodob = f(nodopa)
v(plvp+i) = nodopa*i32 + nodob
iva(i) = nodopa
ivb(i) = nodob
i = i + 1
nodopa = nodob
if (nodopa.ne.ksalva) go to 50
punta = puntlv
puntb = puntlv + 1
puntlv = puntlv + offv + nnodi + 1
p2p = punta
nprob = nnodi
sc1 = 0
return
end if
ierr = 2
return
end
c *************************************************************************
subroutine copyx(f,fc,n)
integer f(n),fc(n)
do 100 i = 1,n
fc(i) = f(i)
100 continue
return
end
c *************************************************************************
subroutine ctcs(n,z,ks,a,f,r,vs,fs,uold,sur,cra,ica,cr,ic,u,v,fb,
+ pi,c,dm,ll,maxdix,maxica,rx,orcr,orcra,inf,ierr)
c solution of the linear min-sum assignment problem through the
c algorithm presented in the paper:
c g. carpaneto and p. toth, "primal-dual algorithms for the
c assignment problem", discrete applied mathematics 18, 137-153,1987
integer ierr,inf,ks,maxdix,maxica,n,orcr,orcra,z
integer a(n,n),c(n),cr(orcr),cra(orcra),dm(n),f(n),fb(n),fs(n),
+ ic(maxdix),ica(maxica),ll(n),pi(n),r(n),rx(n),sur(n),u(n),
+ uold(n),v(n),vs(n)
c .. local scalars ..
integer i,iflag,j,k,kfeas,kopt,krip,kth,m,m1,max,maxdim,
+ maxfea,maxrip,nfeas,nkk,nur,zz
c subroutines apmmix,asmixm,feaso,indus3,opto,setupo
c .. data statements ..
data maxrip/5/,maxfea/1/
c ..
maxdim = maxdix - n
krip = 0
nfeas = 0
do 100 i = 1,n
f(i) = 0
fb(i) = 0
cra(i) = 0
a(i,i) = inf
100 continue
nkk = 0
cr(1) = -1
m1 = -1
call indus3(n,a,f,m,u,v,fb,pi,max,inf)
z = 0
zz = 0
do 200 i = 1,n
z = z + v(i)
zz = zz + u(i)
200 continue
if (m.ne.n) then
if (ks.eq.1) then
c
if (float(m)/float(n).gt.0.6) then
ierr = 0
call setupo(n,a,u,v,m,ic,cr,kth,maxdim,iflag,inf)
if (iflag.ne.1) then
if (cr(1).ge.0) then
do 202 i = 1,n
vs(i) = v(i)
fs(i) = f(i)
uold(i) = u(i)
202 continue
204 call asmixm(n,a,ic,cr,f,fb,u,v,z,pi,r,c,cra,
+ ica,sur,nur,dm,ll,rx,inf)
if (ierr.ne.0) return
if (nur.gt.0) then
if (nfeas.ne.maxfea) then
nfeas = nfeas + 1
call feaso(n,a,kfeas,kth,uold,vs,
+ nfeas,cra,ica,nkk,sur,nur,
+ maxica)
if (kfeas.eq.1) go to 204
end if
do 206 j = 1,n
fb(j) = 0
206 continue
do 208 i = 1,n
j = fs(i)
if (j.gt.0) fb(j) = i
u(i) = uold(i)
v(i) = vs(i)
f(i) = fs(i)
208 continue
else
call opto(n,a,u,v,kopt,cra,ica,nkk,f,fb,
+ uold,maxica,ierr)
if (ierr.eq.0) then
if (kopt.eq.1) return
if (kopt.ne.m1) then
if (krip.ne.maxrip) then
krip = krip + 1
go to 204
end if
end if
end if
end if
end if
end if
end if
end if
call apmmix(n,a,f,z,fb,u,v,pi,r,c,dm,ll,inf)
return
end if
z = 0
do 300 k = 1,n
z = z + u(k) + v(k)
300 continue
return
end
c *************************************************************************
subroutine opto(n,a,u,v,kopt,cra,ica,nkk,f,fb,uold,maxica,ierr)
c check the feasibility of the dual solution , and hence the
c optimality of the primal solution found by the sparse matrix
c procedure.
integer ierr,kopt,maxica,n,nkk
integer a(n,n),cra(*),f(n),fb(n),ica(maxica),u(n),uold(n),v(n)
integer i,ia,j,jc,k,kk,kkn,min,neg
c ..
kopt = 0
kk = nkk
neg = 0
do 100 i = 1,n
if (u(i).ne.uold(i)) then
min = 0
k = i
do 20 j = 1,n
ia = a(i,j) - u(i) - v(j)
if (ia.lt.0) then
if (ia.lt.0) neg = neg + 1
kk = kk + 1
if (kk.le.maxica) then
kkn = kk + n
2 if (cra(k).eq.0) then
cra(k) = kkn
cra(kkn) = 0
ica(kk) = j
else
k = cra(k)
go to 2
end if
end if
if (ia.lt.min) min = ia
end if
20 continue
u(i) = u(i) + min
uold(i) = u(i)
if (min.ne.0) then
jc = f(i)
f(i) = 0
fb(jc) = 0
kopt = kopt - 1
end if
end if
100 continue
nkk = kk
if (nkk.le.maxica) then
if (kopt.lt.0) then
kopt = 0
return
else
kopt = 1
return
end if
end if
ierr = 1
kopt = -1
return
end
c *************************************************************************
subroutine asmixm(n,a,ic,cr,f,fb,dualu,dualv,z,pre,uv,sr,cra,ica,
+ sur,nur,dm,ll,r,inf)
c shortest augmenting path and hungarian method :
c version for sparse matrices
integer inf,n,nur,z
integer a(n,n),cr(*),cra(*),dm(n),dualu(n),dualv(n),f(n),fb(n),
+ ic(*),ica(*),ll(n),pre(n),r(n),sr(n),sur(n),uv(n)
integer aa,d,i,i1,i2,i2p1,ii,iix,ik,ind,index,indexv,iv,j,j1,k,k1,
+ k2,kdir,kmn,knuv,li,nuv,u,w
c ..
nur = 0
do 500 u = 1,n
if (f(u).gt.0) go to 500
do 50 i = 1,n
ll(i) = 0
dm(i) = inf
50 continue
nuv = n
j1 = 0
knuv = 0
k1 = cr(u)
k2 = cr(u+1) - 1
do 100 ik = k1,k2
i = ic(ik)
dm(i) = a(u,i) - dualu(u) - dualv(i)
pre(i) = u
j1 = j1 + 1
sr(j1) = i
100 continue
k = u
150 if (cra(k).eq.0) then
r(1) = u
else
k = cra(k)
kmn = k - n
i = ica(kmn)
dm(i) = a(u,i) - dualu(u) - dualv(i)
pre(i) = u
j1 = j1 + 1
sr(j1) = i
go to 150
end if
200 d = inf
index = 0
i2 = 1
i1 = i2
c
if (nuv.lt.j1) then
kdir = 1
if (knuv.ne.1) then
iv = 0
do 210 i = 1,n
if (ll(i).ne.1) then
iv = iv + 1
uv(iv) = i
end if
210 continue
knuv = 1
end if
do 220 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
index = i
if (d.eq.0) go to 300
end if
i2 = i2 + 1
r(i2) = iv
end if
220 continue
else
kdir = 0
c
do 240 li = 1,j1
i = sr(li)
if (dm(i).le.d) then
if (ll(i).ne.1) then
if (dm(i).ne.d) then
index = 0
i2 = i1
end if
d = dm(i)
if (fb(i).le.0) then
index = i
if (d.eq.0) go to 300
end if
i2 = i2 + 1
r(i2) = i
end if
end if
240 continue
end if
if (d.eq.inf) then
c
nur = nur + 1
sur(nur) = u
go to 500
else if (index.le.0) then
i1 = i1 + 1
i2p1 = i2 + i1
c
do 280 iix = i1,i2
ii = i2p1 - iix
if (kdir.eq.1) then
indexv = r(ii)
index = uv(indexv)
uv(indexv) = uv(nuv)
else
index = r(ii)
end if
ll(index) = 1
nuv = nuv - 1
w = fb(index)
k1 = cr(w)
k2 = cr(w+1) - 1
c
do 250 ik = k1,k2
i = ic(ik)
if (ll(i).ne.1) then
aa = d + a(w,i) - dualu(w) - dualv(i)
if (dm(i).gt.aa) then
if (dm(i).ge.inf) then
j1 = j1 + 1
sr(j1) = i
end if
dm(i) = aa
pre(i) = w
end if
end if
250 continue
k = w
260 if (cra(k).ne.0) then
k = cra(k)
kmn = k - n
i = ica(kmn)
if (ll(i).ne.1) then
aa = d + a(w,i) - dualu(w) - dualv(i)
if (dm(i).gt.aa) then
if (dm(i).ge.inf) then
j1 = j1 + 1
sr(j1) = i
end if
dm(i) = aa
pre(i) = w
end if
end if
go to 260
end if
280 continue
go to 200
end if
c
300 do 350 j = 1,n
if (dm(j).lt.d) then
dualv(j) = dualv(j) + dm(j) - d
i = fb(j)
dualu(i) = dualu(i) - dm(j) + d
end if
350 continue
dualu(u) = dualu(u) + d
c
400 w = pre(index)
fb(index) = w
ind = f(w)
f(w) = index
if (w.ne.u) then
index = ind
go to 400
end if
500 continue
if (nur.le.0) then
c
z = 0
do 550 i = 1,n
z = z + dualu(i) + dualv(i)
550 continue
return
end if
z = -1
return
end
c *************************************************************************
subroutine feaso(n,a,kfeas,kth,us,vs,nfeas,cra,ica,nkk,sur,nur,
+ maxica)
c insert new entries in rows sur (l) (l = 1,nur) in order to find
c a feasible assignment in the sparse cost matrix.
integer kfeas,kth,maxica,n,nfeas,nkk,nur
integer a(n,n),cra(*),ica(maxica),sur(n),us(n),vs(n)
integer coef,i,ia,ii,iukth,iukthn,j,k,kk,kkn,kth0,kthn
data coef/3/
c ..
kk = nkk
kth0 = kth
if (kth0.eq.0) kth0 = 1
do 200 ii = 1,nur
i = sur(ii)
cc = 0.
50 cc = cc + float(coef*nfeas)
kthn = cc*float(kth0)
iukthn = kthn + us(i)
iukth = kth + us(i)
k = i
do 100 j = 1,n
ia = a(i,j) - vs(j)
if (ia.le.iukthn) then
if (ia.gt.iukth) then
kk = kk + 1
if (kk.gt.maxica) go to 300
kkn = kk + n
cra(k) = kkn
cra(kkn) = 0
ica(kk) = j
end if
end if
100 continue
if (cra(i).eq.0) go to 50
200 continue
nkk = kk
kfeas = 1
return
300 kfeas = -1
return
end
c
c *************************************************************************
subroutine indus3(n,a,f,m,u,v,fb,fu,max,inf)
c search for an initial dual solution ( u(i),v(j) ) and an initial
c partial primal solution ( f(i),fb(j) ) of the ap problem
integer inf,m,max,n
integer a(n,n),f(n),fb(n),fu(n),u(n),v(n)
integer i,ia,ii,imin,j,j1,jj,jmin,l,maxl,min
c ..
m = 0
do 100 j = 1,n
u(j) = 0
min = inf
do 50 i = 1,n
ia = a(i,j)
if (ia.le.min) then
if (ia.ge.min) then
if (f(i).ne.0) go to 50
end if
min = ia
imin = i
end if
50 continue
v(j) = min
if (f(imin).eq.0) then
m = m + 1
fb(j) = imin
f(imin) = j
fu(imin) = j + 1
end if
100 continue
max = 0
do 400 i = 1,n
if (f(i).ne.0) go to 400
min = inf
maxl = 0
do 150 j = 1,n
l = a(i,j) - v(j)
if (l.gt.max) maxl = l
if (l.le.min) then
if (l.ge.min) then
if (fb(j).ne.0) go to 150
if (fb(jmin).eq.0) go to 150
end if
min = l
jmin = j
end if
150 continue
u(i) = min
if (maxl-min.gt.max) max = maxl - min
j = jmin
if (fb(j).eq.0) go to 300
do 200 j = jmin,n
if (a(i,j)-v(j).le.min) then
ii = fb(j)
j1 = fu(ii)
if (j1.le.n) then
do 155 jj = j1,n
if (fb(jj).le.0) then
if (a(ii,jj)-u(ii).eq.v(jj)) go to 250
end if
155 continue
fu(ii) = n + 1
end if
end if
200 continue
go to 400
250 f(ii) = jj
fb(jj) = ii
fu(ii) = jj + 1
300 m = m + 1
fb(j) = i
f(i) = j
fu(i) = j + 1
400 continue
return
end
c *************************************************************************
subroutine apmmix(n,a,f,z,fb,dualu,dualv,pre,uv,r,dm,dp,inf)
c shortest augmenting path and hungarian method for complete matrices