← Back to branch summary

~ubuntu-branches/ubuntu/karmic/scilab/karmic

« back to all changes in this revision

Viewing changes to routines/metanet/hamil.f

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2002-03-21 16:57:43 UTC
  • Revision ID: james.westby@ubuntu.com-20020321165743-e9mv12c1tb1plztg
Tags: upstream-2.6
Import upstream version 2.6

Show diffs side-by-side

added added

removed removed

Lines of Context:
routines/metanet/hamil.f
 
1
c     algorithm 595, collected algorithms from acm.
 
2
c     algorithm appeared in acm-trans. math. software, vol.9, no. 1,
 
3
c     mar., 1983, p. 131-138.
 
4
c this program finds one or more hamiltonian circuits in a              
 
5
c directed graph of  n  vertices and  m  arcs.                          
 
6
c input : n number of vertices, ar=ls, pr=lp-1
 
7
      subroutine hamil(n,m,np1,pr,ar,s,nc,nb,
 
8
     * pc,ac,vr,vc,p,subr,rbus,tor)
 
9
      integer pr(np1),pc(np1),ar(m),ac(m),s(n),vr(n),    
 
10
     * vc(n),p(n),subr(n),rbus(n),tor(n)                  
 
11
c exact procedure to find a single hamiltonian circuit, if one exists. 
 
12
      do 1,i=1,n
 
13
         s(i)=0
 
14
 1    continue
 
15
      nc = 1                                                            
 
16
      nb = -1                                                           
 
17
      call hproc(n,pr,ar,nout,nc,nb,s,n+1,pr(n+1),pc,ac,vr, vc, 
 
18
     *p,subr,rbus,tor) 
 
19
c exact procedure to find all the hamiltonian circuits.
 
20
c      nc = -1                                                           
 
21
c       nb = -1                                                           
 
22
c       call hproc(n,pr,ar,nout,nc,nb,s,n+1,pr(n+1),pc,ac,vr,vc, 
 
23
c      * p, subr, rbus, tor)                                              
 
24
c heuristic procedure to find a single hamiltonian 
 
25
c circuit, if one exists, without performing more than 2 backtrackings.
 
26
c       nc = 1                                                            
 
27
c       nb = 2                                                            
 
28
c       call hproc(n,pr,ar,nout,nc,nb,s,n+1, pr(n+1),pc,ac,vr,vc, 
 
29
c      *p,subr,rbus,tor)                                              
 
30
c heuristic procedure to find a single hamiltonian
 
31
c circuit, if one exists, without performing more than 4 backtrackings.
 
32
c       nc = 1                                                            
 
33
c       nb = 4                                                            
 
34
c       call hproc(n,pr,ar,nout,nc,nb,s,n+1,pr(n+1),pc,ac,vr,vc, 
 
35
c      *p,subr,rbus,tor)                                              
 
36
c heuristic procedure to find  at most  2
 
37
c hamiltonian circuits, without performing more than  5  backtrackings. 
 
38
c       nc = 2                                                            
 
39
c       nb = 5                                                            
 
40
c       call hproc(n,pr,ar,-1,nc,nb,s,n+1,pr(n+1),pc,ac,vr,vc,   
 
41
c      *p,subr,rbus,tor)                                              
 
42
c       if (nc.eq.0) return
 
43
      return
 
44
      end                                                               
 
45
      subroutine hproc(n,pr,ar,kw,nc,nb,s,np1,m,pc,ac,vr,vc,
 
46
     *p,subr,rbus,tor)
 
47
c
 
48
c subroutine to find one or more hamiltonian circuits in a
 
49
c directed graph of  n  vertices ( n .gt. 1 ) represented
 
50
c by the integers  1, 2, ..., n  and  m  arcs.
 
51
c
 
52
c hproc is based on an enumerative algorithm and can be used
 
53
c either as an exact procedure or as a heuristic procedure
 
54
c (by limiting the number of backtrackings allowed).
 
55
c
 
56
c entrance to hproc is achieved by using the statement
 
57
c     call hproc(n,pr,ar,kw,nc,nb,s,n+1,pr(n+1),pc,ac,vr,vc,
 
58
c    *        p,subr,rbus,tor)
 
59
c
 
60
c the values of the first six parameters must be defined
 
61
c by the user prior to calling hproc. hproc needs  2  arrays ( pr
 
62
c and  pc ) of length  n + 1 ,  2  arrays ( ar  and  ac )
 
63
c of length  m  and  7  arrays ( s ,  vr ,  vc ,  subr ,
 
64
c rbus  and  tor ) of length  n . these arrays must be
 
65
c dimensioned by the user in the calling program.
 
66
c
 
67
c hproc calls  5  subroutines: pathp, fupd, bupd, iupd, rarc.
 
68
c meaning of the input parameters:
 
69
c n     = number of vertices.
 
70
c pr(i) = sum of the out-degrees of vertices  1, ..., i-1
 
71
c         ( pr(1) = 0 ,  pr(n+1) = m ).
 
72
c ar    = adjacency list. the elements from  ar(pr(i)+1)  to
 
73
c         ar(pr(i+1))  are a record containing,in any order,
 
74
c         all the vertices  j  such that arc  (i,j)  exists.
 
75
c         the graph should not contain arcs starting and
 
76
c         ending at the same vertex.
 
77
c kw    = unit tape on which to write the hamiltonian cir-
 
78
c         cuits found, according to format  20i5 .  kw = - 1
 
79
c         if no writing is desired. the circuits are written
 
80
c         as ordered sequences of  n  vertices.
 
81
c
 
82
c meaning of the input-output parameters:
 
83
c nc(input)  = upper bound on the number of hamiltonian
 
84
c              circuits to be found ( nc = - 1 if all the
 
85
c              hamiltonian circuits are to be found).
 
86
c nc(output) = number of hamiltonian circuits found.
 
87
c nb(input)  = - 1  if hproc must be executed as an exact
 
88
c              procedure.
 
89
c            = upper bound on the number of backtrackings if
 
90
c              hproc must be executed as a heuristic procedure.
 
91
c nb(output) = number of backtrackings performed. when hproc
 
92
c              has been executed as a heuristic procedure,
 
93
c              if  nb(output) .lt. nb(input)  then the
 
94
c              result obtained is exact.
 
95
c
 
96
c meaning of the output parameter:
 
97
c s(i) = i-th  vertex in the last hamiltonian circuit found.
 
98
c
 
99
c on return of hproc  n, pr  and  kw  are unchanged, while in
 
100
c ar  the order of the elements within each record may be
 
101
c altered.
 
102
c
 
103
c meaning of the work arrays:
 
104
c pc(i)   = sum of the in-degrees of vertices  1, ..., i-1
 
105
c           ( pc(1) = 0 ).
 
106
c ac      = adjacency list (backward). the elements from
 
107
c           ac(pc(i)+1)  to  ac(pc(i+1))  contain, in any
 
108
c           order, all the vertices  j  such that arc  (j,i)
 
109
c           exists.
 
110
c when an arc is removed from the graph at the  k-th  level
 
111
c of the branch-decision tree, the corresponding elements
 
112
c ar(q)  and  ac(t)  are set to  - (k*(n+1) + ar(q))  and
 
113
c to  - (k*(n+1) + ac(t)) , respectively.
 
114
c vr(i)   = current out-degree of vertex  i .
 
115
c vc(i)   = current in-degree of vertex  i .
 
116
c subr(i) = - (k*(n+1) + j)  if arc  (i,j)  was implied at
 
117
c           the  k-th  level of the branch-decision tree.
 
118
c         = 0  otherwise.
 
119
c rbus(i) = - j  if arc  (j,i)  is currently implied.
 
120
c         = 0  otherwise.
 
121
c tor(k)  = q*(m+1) + t  if the arc going from  s(k)  to the
 
122
c           root, corresponding to  ar(q)  and to  ac(t),
 
123
c           was removed from the graph at the  k-th  level
 
124
c           of the branch-decision tree.
 
125
c         = 0  otherwise.
 
126
c p(i)    = pointer for the forward step. the next arc
 
127
c           starting from  i  to be considered in the
 
128
c           branch-decision tree is  (i,ar(pr(i)+p(i)).
 
129
c
 
130
c meaning of the main work simple variables:
 
131
c jr  = root. the hamiltonian circuits are determined as
 
132
c       paths starting and ending at  jr .
 
133
c k   = current level of the branch-decision tree.
 
134
c m   = number of arcs.
 
135
c mp1 = m + 1 (used for packing  tor ).
 
136
c np1 = n + 1 (used for packing  ar ,  ac  and  subr )
 
137
c
 
138
      integer pr(np1), pc(np1), ar(m), ac(m), s(n), vr(n), vc(n), p(n),
 
139
     * subr(n), rbus(n), tor(n)
 
140
c
 
141
c s t e p   0   (initialize).
 
142
c
 
143
      nco = nc
 
144
      nc = 0
 
145
      nbo = nb
 
146
      nb = 0
 
147
      do 10 i=1,n
 
148
        vc(i) = 0
 
149
        subr(i) = 0
 
150
        rbus(i) = 0
 
151
        p(i) = 1
 
152
   10 continue
 
153
      do 30 i=1,n
 
154
        j1 = pr(i) + 1
 
155
        j2 = pr(i+1)
 
156
        vr(i) = j2 - j1 + 1
 
157
        if (vr(i).eq.0) return
 
158
        do 20 j=j1,j2
 
159
          ja = ar(j)
 
160
          vc(ja) = vc(ja) + 1
 
161
   20   continue
 
162
   30 continue
 
163
      pc(1) = 0
 
164
      do 40 i=1,n
 
165
        if (vc(i).eq.0) return
 
166
        pc(i+1) = pc(i) + vc(i)
 
167
        vc(i) = 0
 
168
   40 continue
 
169
      do 60 i=1,n
 
170
        j1 = pr(i) + 1
 
171
        j2 = pr(i+1)
 
172
        do 50 j=j1,j2
 
173
          jj = ar(j)
 
174
          vc(jj) = vc(jj) + 1
 
175
          ja = pc(jj) + vc(jj)
 
176
          ac(ja) = i
 
177
   50   continue
 
178
   60 continue
 
179
      mp1 = m + 1
 
180
c select as root  jr  the vertex  i  with maximum  vc(i)
 
181
c (break ties by choosing  i  with minimum  vr(i) ).
 
182
      maxe = vc(1)
 
183
      minu = vr(1)
 
184
      jr = 1
 
185
      do 100 i=2,n
 
186
        if (vc(i)-maxe) 100, 70, 80
 
187
   70   if (vr(i).ge.minu) go to 100
 
188
        go to 90
 
189
   80   maxe = vc(i)
 
190
   90   minu = vr(i)
 
191
        jr = i
 
192
  100 continue
 
193
      k1 = -np1
 
194
      k = 1
 
195
      s(1) = jr
 
196
c
 
197
c s t e p   1   (search for implied arcs).
 
198
c
 
199
  110 do 120 j=1,n
 
200
        if (vr(j).eq.1) go to 130
 
201
        if (vc(j).eq.1) go to 170
 
202
  120 continue
 
203
c no further arc can be implied.
 
204
      go to 220
 
205
c arc  (j,jl)  is implied because  vr(j) = 1 .
 
206
  130 l1 = pr(j) + 1
 
207
      l2 = pr(j+1)
 
208
      do 140 l=l1,l2
 
209
        if (ar(l).gt.0) go to 150
 
210
  140 continue
 
211
  150 jl = ar(l)
 
212
c find the starting vertex  it1  and the ending vertex  it2
 
213
c of the largest path of implied arcs containing  (j,jl) .
 
214
      call pathp(j, jl, subr, rbus, ar, pr, s, n, np, it1, it2, k, jr,
 
215
     * m, np1)
 
216
      if (np.eq.0) go to 160
 
217
      if (np.eq.(-1)) go to 340
 
218
c subroutine pathp found a hamiltonian circuit.
 
219
      k = k + 1
 
220
      go to 320
 
221
  160 subr(j) = k1 - jl
 
222
      rbus(jl) = j
 
223
c remove from the graph all arcs terminating at  jl .
 
224
      call iupd(j, jl, l, ac, ar, pc, pr, vc, vr, k1, n, m, np1)
 
225
      if (j.eq.0) go to 340
 
226
      go to 210
 
227
c arc  (jl,j)  is implied because  vc(j) = 1 .
 
228
  170 l1 = pc(j) + 1
 
229
      l2 = pc(j+1)
 
230
      do 180 l=l1,l2
 
231
        if (ac(l).gt.0) go to 190
 
232
  180 continue
 
233
  190 jl = ac(l)
 
234
c find the starting vertex  it1  and the ending vertex  it2
 
235
c of the largest path of implied arcs containing  (jl,j) .
 
236
      call pathp(jl, j, subr, rbus, ar, pr, s, n, np, it1, it2, k, jr,
 
237
     * m, np1)
 
238
      if (np.eq.0) go to 200
 
239
      if (np.eq.(-1)) go to 340
 
240
c subroutine pathp found a hamiltonian circuit.
 
241
      i = s(k)
 
242
      k = k + 1
 
243
      go to 320
 
244
  200 subr(jl) = k1 - j
 
245
      rbus(j) = jl
 
246
c remove from the graph all arcs emanating from  jl .
 
247
      call iupd(j, jl, l, ar, ac, pr, pc, vr, vc, k1, n, m, np1)
 
248
      if (j.eq.0) go to 340
 
249
c if arc  (it2,it1)  is in the graph, remove it.
 
250
  210 call rarc(it2, it1, ar, ac, pr, pc, vr, vc, k1, jj, ll, n, m, np1)
 
251
      if (jj.eq.(-1)) go to 340
 
252
      go to 110
 
253
c
 
254
c s t e p   2   (add implied arcs to  s ).
 
255
c
 
256
  220 i = s(k)
 
257
      if (subr(i).eq.0) go to 230
 
258
      jsubr = -subr(i) + subr(i)/np1*np1
 
259
      if (jsubr.eq.jr) go to 340
 
260
      k = k + 1
 
261
      s(k) = jsubr
 
262
      if (k.ne.n) go to 220
 
263
      if (subr(jsubr).lt.0) go to 320
 
264
      go to 340
 
265
c
 
266
c s t e p   3   (branch).
 
267
c
 
268
  230 l1 = pr(i) + p(i)
 
269
      l2 = pr(i+1)
 
270
      if (l1.gt.l2) go to 340
 
271
c find the next arc  (i,jl)  to be added to  s .
 
272
      dens = n**3
 
273
      j1 = 0
 
274
      j2 = 0
 
275
      do 310 j=l1,l2
 
276
        jl = ar(j)
 
277
        if (jl.lt.0) go to 310
 
278
        if (vr(jl).gt.0) go to 260
 
279
        if (subr(jl).eq.0) go to 310
 
280
        if (jl.eq.jr) go to 310
 
281
        iend = jl
 
282
  240   iend = -subr(iend) + subr(iend)/np1*np1
 
283
        if (subr(iend).ne.0) go to 240
 
284
        if (vc(jl).lt.vr(iend)) go to 250
 
285
        score = vr(iend)*n + vc(jl)
 
286
        go to 280
 
287
  250   score = vc(jl)*n + vr(iend)
 
288
        go to 280
 
289
  260   if (vc(jl).lt.vr(jl)) go to 270
 
290
        score = vr(jl)*n + vc(jl)
 
291
        go to 280
 
292
  270   score = vc(jl)*n + vr(jl)
 
293
  280   if (dens.le.score) go to 290
 
294
        dens = score
 
295
        ipi = j
 
296
  290   if (j1.eq.0) go to 300
 
297
        if (j2.eq.0) j2 = j
 
298
        go to 310
 
299
  300   j1 = j
 
300
  310 continue
 
301
      if (j1.eq.0) go to 340
 
302
      jl = ar(ipi)
 
303
      ar(ipi) = ar(j1)
 
304
      ar(j1) = jl
 
305
      if (j2.eq.0) j2 = pr(i+1) + 1
 
306
      p(i) = j2 - pr(i)
 
307
      k = k + 1
 
308
      s(k) = jl
 
309
      k1 = -k*np1
 
310
c remove from the graph all arcs emanating from  i .
 
311
      call fupd(ar, ac, pr, pc, vr, vc, i, k1, n, m, np1)
 
312
c remove from the graph all arcs terminating at jl .
 
313
      call fupd(ac, ar, pc, pr, vc, vr, jl, k1, n, m, np1)
 
314
      tor(k) = 0
 
315
c if arc  (jl,jr)  is in the graph, remove it.
 
316
      call rarc(jl, jr, ar, ac, pr, pc, vr, vc, k1, jj, ll, n, m, np1)
 
317
      if (jj.eq.0) go to 110
 
318
      if (jj.eq.(-1)) go to 340
 
319
      tor(k) = jj*mp1 + ll
 
320
      go to 110
 
321
c
 
322
c s t e p   4   (hamiltonian circuit found).
 
323
c
 
324
  320 nc = nc + 1
 
325
      if (kw.eq.(-1)) go to 330
 
326
c ww s(kj),kj=1,n)
 
327
  330 if (nc.eq.nco) go to 430
 
328
      k = k - 1
 
329
c
 
330
c s t e p   5   (backtrack).
 
331
c
 
332
  340 if (k.le.1) go to 430
 
333
      ja = s(k)
 
334
      p(ja) = 1
 
335
      ja = s(k-1)
 
336
      if (subr(ja).eq.0) go to 350
 
337
c backtracking for an implied arc.
 
338
      k = k - 1
 
339
      go to 340
 
340
  350 if (nb.eq.nbo) go to 430
 
341
      nb = nb + 1
 
342
      k1 = -k*np1
 
343
      k2 = -(k+1)*np1
 
344
      i = s(k-1)
 
345
c backtracking for the arcs implied at level  k .
 
346
      iff = 0
 
347
      do 360 j=1,n
 
348
        if (subr(j).gt.k1) go to 360
 
349
        if (subr(j).lt.k2) go to 360
 
350
        ja = k1 - subr(j)
 
351
        rbus(ja) = 0
 
352
        subr(j) = 0
 
353
        iff = 1
 
354
  360 continue
 
355
      if (iff.eq.1) go to 370
 
356
c no arc was implied at level  k .
 
357
      call bupd(ar, ac, pr, pc, vr, vc, i, k1, k2, n, m, np1)
 
358
      call bupd(ac, ar, pc, pr, vc, vr, s(k), k1, k2, n, m, np1)
 
359
      if (tor(k).eq.0) go to 420
 
360
      j1 = tor(k)/mp1
 
361
      j2 = tor(k) - j1*mp1
 
362
      ar(j1) = jr
 
363
      ja = s(k)
 
364
      vr(ja) = vr(ja) + 1
 
365
      ac(j2) = s(k)
 
366
      vc(jr) = vc(jr) + 1
 
367
      go to 420
 
368
c at least one arc was implied at level  k .
 
369
  370 do 410 j=1,n
 
370
        l1 = pr(j) + 1
 
371
        l2 = pr(j+1)
 
372
        do 400 l=l1,l2
 
373
          jl = ar(l)
 
374
          if (jl.gt.k1) go to 400
 
375
          if (jl.lt.k2) go to 400
 
376
          jl = k1 - jl
 
377
          ar(l) = jl
 
378
          vr(j) = vr(j) + 1
 
379
          ll1 = pc(jl) + 1
 
380
          ll2 = pc(jl+1)
 
381
          do 380 ll=ll1,ll2
 
382
            if (k1-ac(ll).eq.j) go to 390
 
383
  380     continue
 
384
  390     ac(ll) = j
 
385
          vc(jl) = vc(jl) + 1
 
386
  400   continue
 
387
  410 continue
 
388
  420 k = k - 1
 
389
      go to 230
 
390
c
 
391
c re-store the original vector  ar .
 
392
c
 
393
  430 do 440 j=1,m
 
394
        if (ar(j).gt.0) go to 440
 
395
        ar(j) = -ar(j) + ar(j)/np1*np1
 
396
  440 continue
 
397
      return
 
398
      end
 
399
      subroutine pathp(i,j,subr,rbus,ar,pr,s,n,np,i1,i2,k,    
 
400
     * jr, m, np1)
 
401
c subroutine to find the starting vertex  i1  and the ending
 
402
c vertex  i2  of the largest path of implied arcs containing
 
403
c arc  (i,j) .
 
404
c meaning of the output parameter  np :
 
405
c np =  0  if the path contains  l .lt. n  vertices.
 
406
c    =  1  if the path contains  n  vertices and arc (i2,i1)
 
407
c          exists (the hamiltonian circuit is stored in  s )
 
408
c    = -1  if the path contains  n  vertices but arc (i2,i1)
 
409
c          does not exist.
 
410
      integer subr(n), rbus(n), ar(m), pr(np1), s(n)
 
411
      np = 0
 
412
      l = 1
 
413
      i1 = i
 
414
   10 if (rbus(i1).eq.0) go to 20
 
415
      i1 = rbus(i1)
 
416
      l = l + 1
 
417
      go to 10
 
418
   20 i2 = j
 
419
      l = l + 1
 
420
   30 if (subr(i2).eq.0) go to 40
 
421
      i2 = -subr(i2) + subr(i2)/np1*np1
 
422
      l = l + 1
 
423
      go to 30
 
424
   40 continue
 
425
      if (l.lt.n) return
 
426
c the path contains  n  vertices.
 
427
      k1 = -k*np1
 
428
      l1 = pr(i2) + 1
 
429
      l2 = pr(i2+1)
 
430
      do 60 l=l1,l2
 
431
        if (ar(l).lt.0) go to 50
 
432
        if (ar(l).eq.i1) go to 70
 
433
        go to 60
 
434
   50   if (k1-ar(l).eq.i1) go to 70
 
435
   60 continue
 
436
c no hamiltonian circuit can be determined.
 
437
      np = -1
 
438
      return
 
439
c a hamiltonian circuit exists. store it in  s .
 
440
   70 np = 1
 
441
      rbus(j) = i
 
442
      rbus(i1) = i2
 
443
      s(n) = rbus(jr)
 
444
      l = n - 1
 
445
   80 if (l.eq.k) go to 90
 
446
      ja = s(l+1)
 
447
      s(l) = rbus(ja)
 
448
      l = l - 1
 
449
      go to 80
 
450
   90 rbus(i1) = 0
 
451
      rbus(j) = 0
 
452
      return
 
453
      end
 
454
      subroutine fupd(a1, a2, p1, p2, v1, v2, i1, k1, n, m, np1)
 
455
c forward step updating
 
456
      integer a1(m), a2(m), p1(np1), p2(np1), v1(n), v2(n)
 
457
      j1 = p1(i1) + 1
 
458
      j2 = p1(i1+1)
 
459
      do 30 j=j1,j2
 
460
        if (a1(j).lt.0) go to 30
 
461
        ia = a1(j)
 
462
        l1 = p2(ia) + 1
 
463
        l2 = p2(ia+1)
 
464
        do 10 l=l1,l2
 
465
          if (a2(l).eq.i1) go to 20
 
466
   10   continue
 
467
   20   v2(ia) = v2(ia) - 1
 
468
        a2(l) = k1 - a2(l)
 
469
        a1(j) = k1 - ia
 
470
   30 continue
 
471
      v1(i1) = 0
 
472
      return
 
473
      end
 
474
      subroutine bupd(a1, a2, p1, p2, v1, v2, ii, k1, k2, n, m, np1)
 
475
c backtracking step updating
 
476
      integer a1(m), a2(m), p1(np1), p2(np1), v1(n), v2(n)
 
477
      l1 = p1(ii) + 1
 
478
      l2 = p1(ii+1)
 
479
      do 30 l=l1,l2
 
480
        if (a1(l).gt.k1) go to 30
 
481
        if (a1(l).lt.k2) go to 30
 
482
        ia = k1 - a1(l)
 
483
        a1(l) = ia
 
484
        v1(ii) = v1(ii) + 1
 
485
        ll1 = p2(ia) + 1
 
486
        ll2 = p2(ia+1)
 
487
        do 10 ll=ll1,ll2
 
488
          if (k1-a2(ll).eq.ii) go to 20
 
489
   10   continue
 
490
   20   a2(ll) = ii
 
491
        v2(ia) = v2(ia) + 1
 
492
   30 continue
 
493
      return
 
494
      end
 
495
      subroutine iupd(ia, ib, l, a1, a2, p1, p2, v1, v2, k1, n, m, np1)
 
496
      integer a1(m), a2(m), p1(np1), p2(np1), v1(n), v2(n)
 
497
c updating for implied arc
 
498
      m1 = p1(ib) + 1
 
499
      m2 = p1(ib+1)
 
500
      do 40 mm=m1,m2
 
501
        iarc = a1(mm)
 
502
        if (iarc.lt.0) go to 40
 
503
        if (v2(iarc).ne.1) go to 10
 
504
        if (iarc.ne.ia) go to 50
 
505
        jj = l
 
506
        go to 30
 
507
   10   j1 = p2(iarc) + 1
 
508
        j2 = p2(iarc+1)
 
509
        do 20 jj=j1,j2
 
510
          if (a2(jj).eq.ib) go to 30
 
511
   20   continue
 
512
   30   a2(jj) = k1 - a2(jj)
 
513
        v2(iarc) = v2(iarc) - 1
 
514
        a1(mm) = k1 - iarc
 
515
        v1(ib) = v1(ib) - 1
 
516
   40 continue
 
517
      return
 
518
   50 ia = 0
 
519
      return
 
520
      end
 
521
      subroutine rarc(ia, ib, ar, ac, pr, pc, vr, vc, k1, jj, ll, n, m,
 
522
     * np1)
 
523
c subroutine to remove arc  (ia,ib)  from the graph.
 
524
c meaning of the output parameters  jj  and  ll :
 
525
c jj =  location of the element of  ar  corresponding to the removed arc.
 
526
c    =  0  if arc  (ia,ib)  is not in the graph.
 
527
c    = -1  if, after the removal of arc  (ia,ib) , the graph
 
528
c          would admit no hamiltonian circuit.
 
529
c ll =  location of the element of  ac  corresponding to the
 
530
c       removed arc (defined only if  jj .gt. 0 ).
 
531
      integer ar(m), ac(m), pr(np1), pc(np1), vr(n), vc(n)
 
532
      j1 = pr(ia) + 1
 
533
      j2 = pr(ia+1)
 
534
      do 20 jj=j1,j2
 
535
        if (ar(jj).lt.0) go to 20
 
536
        if (ar(jj).ne.ib) go to 20
 
537
        l1 = pc(ib) + 1
 
538
        l2 = pc(ib+1)
 
539
        do 10 ll=l1,l2
 
540
          if (ac(ll).eq.ia) go to 30
 
541
   10   continue
 
542
   20 continue
 
543
c arc  (ia,ib)  is not in the graph.
 
544
      jj = 0
 
545
      return
 
546
   30 if (vr(ia).eq.1) go to 40
 
547
      if (vc(ib).eq.1) go to 40
 
548
      ar(jj) = k1 - ib
 
549
      vr(ia) = vr(ia) - 1
 
550
      ac(ll) = k1 - ia
 
551
      vc(ib) = vc(ib) - 1
 
552
      return
 
553
c arc  (ia,ib)  cannot be removed from the graph.
 
554
   40 jj = -1
 
555
      return
 
556
      end
 
557