program bisect
c **************************************************************************
c GP.F
c This program finds upper and lower bounds on the optimal solution of a
c Graph Bisection problem (i.e. m1 = m2 = n/2). Rendl-Wolkowicz upper
c bounds are found using a Bundle-Trust algorithm to perform the
c minimization, and a block Lanczos algorithm to calculate the eigenvalues.
c Lower bounds are found by doing a sort, and then applying the Kernighan-
c Lin heuristic to improve the resulting bisection.
c
c Other files required: subs.for
c
c Input file: The name of this file is provided by the user. It contains
c a sparse representation of the adjacency matrix A.
c
c When running the program, the following screen input is required:
c filename
c - The file must contain a sparse representation of the adjacency
c matrix.
c maxit
c - The maximum number of BT iterations allowed; suggestion: 30.
c set up d automatically or read from file?
c - Normally set up automatically; enter 0 for this.
c scaling factor for d
c - An initial choice for d is needed to start the BT iterations.
c This scaling factor is used to adjust the initial choice provided
c by the program. Often a good choice is 0.8.
c estimate of final upper bound
c - If 0 is entered, the program averages the initial upper bound
c and the initial feasible solution cost.
c
c At the end of the run, enter 1 if you want to save the optimal d (you
c will be asked for a filename) or some other integer to stop.
c
c Storage of the Adjacency Matrix
c This matrix is sparse, symmetric (so only the upper triangle is
c stored) and has zeros on the diagonal. The input file contains:
c n nzl nzs
c icol(1:nzl)
c irow(i) a(i) i= 1,nzs
c
c n = # of nodes
c nzl = # of last row containing non-zero above main diagonal
c nzs = total # of nonzeros above main diagonal
c icol(i) = for each row i, # non-zeros after main diagonal
c irow(i) = column number for the ith non-zero
c a(i) = weight (adjacency matrix entry) for the ith non-zero
c
c Variables and Data Structures
c n = number of nodes
c maxcom = For BT: maximum number of calls to COMPFG
c maxit = maximum number of BT iterations allowed
c iprint = controls level of printout from BT
c ifail = gives info about the result of BT
c iwork() = work space for BT
c icol(i) = For each row i, # non-zeros after main diagonal
c irow(i) = Column for the ith non-zero
c ibest() = The first subset for the best feasible partition found
c jbest() = The second subset for the best feasible partition found
c ipkeep() = Temporary storage of feasible partition
c iqkeep() = Temporary storage of feasible partition
c ad() = diagonal elements of adjacency matrix (now assumed to be 0)
c a(i) = ith non-zero weight (adjacency matrix entry)
c suma = sum of the entries in the adjacency matrix
c d = the diagonal perturbation of the A matrix
c f = the function value: f(d) = largest eigenvalue of
c Vn' * (A + diag(d)) * Vn
c g() = gradient of f
c fm = estimate of the final upper bound
c eps = stopping criterion for BT
c work() = work space for BT
c worklz() = work space for DNLASO (Lanczos algorithm)
c xconst = constant for special representation of Vn, ie -1/(n+sqrt(n))
c yconst = constant for special representation of Vn, ie -1/sqrt(n)
c vec(i,j) = the eigenvectors found by Lanczos
c bstkl = cost of best feasible solution after Kernighan-Lin
c bstev = cost of best feasible solution before Kernighan-Lin
c nfig = number of decimal digits of accuracy required from Lanczos
c nblock = number of vectors in each Lanczos block
c nval = number of eigenvalues desired
c nperm = number of eigenvalues successfully found by Lanczos
c matmul = total number of matrix accesses in BT
c infile = name of input file
c istop = tells BT method to abort if we get stuck in Lanczos
c
c last changed: July 1994 Franz Rendl
c **************************************************************************
integer nmax, reset, lwork, ntot
c-------- needs to be adapted for larger/smaller problems -(1)-
parameter (nmax = 5000, ntot = 30000, reset = 10)
parameter (nmval= 6, lwork = 100000)
integer n, i, maxcom, maxit, iprint, ifail, iwork( reset)
integer icol( nmax), irow( ntot), ibest( nmax/2)
integer ipkeep( nmax/2), iqkeep( nmax/2), jbest( nmax/2)
double precision ad( nmax), a( ntot), suma,
* d( nmax), f, g(nmax), fm, eps, work( lwork)
double precision xconst, yconst, vec(nmax, nmval), bstkl, bstev
real dummy( 5), xbst
c
c work space for dnlaso
double precision worklz(800+21*nmax)
c
character*20, infile
external compfg
common / graph1 / a, ad, suma
common / graph2 / nzl, icol, irow
common / const / xconst, yconst
c
common / iconst / matmul, nfig, nblock, nval, nperm
common / stop / istop, idummy
c
common / lzwork / worklz
common / wkspc / work
common / eij / vec
common / globl / bstev, bstkl
common / globl2 / ibest, jbest
common / keep / ipkeep, iqkeep
common / eigvls / dummy, xbst
c If run a series of tests, might choose to initialise random
c number generator each time (for repeatibility)
common /random/iurand
c
c Initialise Lanczos work array
do 5 i = 1, 800+21*nmax
worklz( i) = 0.0d0
5 continue
c bstev and bstkl contain current champions of best solutions
c with respect to eigenvectors and kernighan-lin
bstev = -1.0d10
bstkl = bstev
c *** read description of graph
print *,'filename:'
read '(a20)', infile
open( 4, file = infile, status = 'old')
read(4,*) n, nzl, nzs
c check maximum dimensions!
if (n .gt. nmax .or. nzs .gt. ntot) then
write(*,*)' Maximum dimensions exceeded!'
write(*,*)' Max number of nodes = ',nmax,' Actual = ',n
write(*,*)' Max number of edges = ',ntot,' Actual = ',nzs
stop
endif
read(4,*) (icol( i), i=1,nzl)
c check that the sum of the non-zeros per column equals the
c total number of non-zeros, nzs
isum = 0
do 6 i = 1,nzl
isum = isum + icol(i)
6 continue
if (isum .ne. nzs) then
write(*,*)' Data Error!'
write(*,*)' Number of non-zeros specified as ', nzs
write(*,*)' Sum of non-zeros per column is ',isum
stop
endif
read(4,*) (irow( i), a( i), i=1,nzs)
c read(4,*) (ad( i), i=1,n)
c main diag of graph is assumed to be zero, i.e. ad( i) = 0.0
close(4)
print *,'graph was read successully. It has'
print 7,n,nzs
7 format(2x,i6,' nodes and ',i7,' edges.')
c Some parameters are hard-coded
eps = 0.1
nfig = 5
nblock = 3
nval = 2
maxit = 15
c print *,'maxit ='
c read *,maxit
c Set up constants for dnlaso
c Set up constants for bt and Vn
matmul = 0
iprint = 2
maxcom = maxit + 20
yconst = sqrt( float( n))
xconst = - 1.0d0 / (float( n) + yconst)
yconst = - 1.0d0 / yconst
c Compute initial diagonal perturbation and sum of entries (suma)
do 14 i = 1,n
14 d( i) = 0.0d0
suma = 0.0
llast = 0
do 16 j = 1, nzl
if ( icol( j) .eq. 0) goto 16
lfirst = llast + 1
llast = llast + icol( j)
do 15 l = lfirst, llast
i = irow( l)
d( i) = d( i) + a( l)
d( j) = d( j) + a( l)
suma = suma + a( l)
15 continue
16 continue
suma = 2.0d0 * suma
do 17 i = 1,n
d( i) = d( i) - suma / float( n)
17 continue
c *** this is a good choice in general, on special problems it might
c be either chosen close to one, a trivial choice is 0
factr = .8
c for very sparse graphs factr should be close to one
if ( float(nzs) / float( n) .le. 3.) factr = .95
if (float(nzs) / float( n) .gt. 5.) factr = .6
if (float(nzs) / float(n) .gt. 10.) factr = .4
do 19 i = 1,n
19 d( i) = d( i) * factr
c Calculate an initial upper bound
call compfg( n, d, f, g)
print 8, f
8 format(5x,'initial upper bound = ', f20.4)
nval = 1
nblock = 2
c
c *******************
c Use first two eigenvectors to find an initial feasible bisection
i = 1
j = 2
print *,'use first two eigenvectors to get bisection:'
call feas(n, nzs, i, j )
c pass best solution on to BT for printout
xbst = bstkl
c *******************
c estimate of final bound
fm = 0.7 * f + 0.3 *bstkl
print 9, fm
9 format(3x,'initial guess on fm=', f20.3)
c Now perform BT iterations to improve the upper bound
call bt( n, d, f, g, compfg, fm, eps, maxcom, maxit, reset,
* iprint, ifail, iwork, work, lwork)
if (ifail .eq. 4) then
write(*,*)' BT Code does not have enough workspace'
stop
endif
c Recalculate final upper bound
nfig = 7
nval = 3
nblock = nval + 1
print 11, nfig
11 format(1x,'final check: 3 largest eigs, accuracy=', i3)
call compfg( n, d, f, g)
print 10,f
10 format(1x,'final upper bound:',f12.3)
c *******************
c Calculate a final feasible solution. Use pairs of eigenvectors.
c Note that Lanczos has only calculated nval eigenvectors. The number
c reliable eigenvectors found, nperm, may be less than nval, but we can
c use even the unreliable estimates to help find a good feasible solution.
print *,'run around circles using first eigenvectors:'
do 30 i = 1, nval - 1
do 30 j = i+1, nval
30 call feas( n, nzs, i, j)
c Now let's see what we have got:
c First round down bound to nearest integer
ibnd = f
print 802, ibnd
802 format(//2x,'** final bound (rounded down)=',i10)
print 801, bstkl
801 format(2x,'** best bisection=',f10.2)
print 803, (ibnd- bstkl)/ float(bstkl)
803 format(2x,'** final gap=' , f5.3)
c *******************
print *,'total number of matrix accesses:', matmul
c print *,'1 = save optimal d and best bisection, else stop:'
c read *,i
c this can be used to store the final result
c the two sidesof the partition are in ibest and jbest
i = 0
if (i .ne. 1) stop
print *,'filename:'
read '(a20)', infile
open( 4, file = infile, status = 'unknown')
write( 4, 90) (d(i),i=1,n)
write( 4, 91) (ibest( i), i = 1, n/2)
write( 4, 91) (jbest( i), i = 1, n/2)
91 format( 20i4)
90 format(5(1x, f12.5))
close(4)
end
c----------------------- feas.for -------------
subroutine feas(n, nzs, col1, col2)
c *** take eigenvectors in columns col1 and col2 and do:
c 1) run around circle and find best bisection, store it in *keep
c 2) take best bisection found in 1) and apply Ker-Lin
c
c We use a different representation of the matrix here, because it saves
c time in the Kernighan-Lin routine. We pretend that the matrix is not
c symmetric and store all the non-zero elements in each row.
c This new sparse matrix representation is as follows:
c nnz = number of non-zero entries
c weight - array containing the nonzero entries in row order
c iadj - indicates which column of the matrix the corresponding
c entry of weight is in
c ifirst - array of pointers into weight
c ifirst(i) indicates first weight entry for row i
c ilast - array of pointers into weight
c ilast(i) indicates last weight entry for row i
c part - This array indicates for each node which side of the
c partition it is on, i.e. ip or iq. Specifically,
c part(i) = +j if node i is in the jth position of the
c ip side of the partition and part(i) = -j if node i
c is in the jth position of the iq side.
c ----------- needs to be adapted ---------- (2) ----
parameter (nmax = 5000, ntot = 30000, nmval = 6)
parameter (nn = nmax/2, not0 = ntot*2)
integer col1, col2
integer ip(nn), iq(nn), part(nmax), icol(nmax), irow(ntot),
* ifirst(nmax), ilast(nmax), ibest(nn), jbest(nn),
* iadj(not0), kp(nn), kq(nn),ideg(nmax),
* ipkeep(nn), iqkeep(nn), iadr(nmax)
double precision sum1, sum2, x, vec(nmax,nmval), z(nmax),
* y(nmax), a(ntot), ad( nmax), suma, t, w(nmax),
* wk1(nn), wk2(nn), wk3(nn), gain(nn), sum,
* weight(not0), tcost, bstev, bstkl, cost, best
logical iwk4( nn), iwk5( nn)
common / eij / vec
common / graph1 / a, ad, suma
common / graph2 / nzl, icol, irow
common / keep / ipkeep, iqkeep
common / globl / bstev, bstkl
common / globl2 / ibest, jbest
common / wkspc / weight
c
c Calculate the new matrix representation
call trafo(n, nzl, nzs, icol, irow, a, ideg, ifirst, ilast,
* iadj, weight, iadr)
best = 0.0d0
n2 = n/2
c form z = Vn * vec(.,col1) (unproject)
c and y = Vn * vec(.,col2) (for the second eigenvector)
c first sum the elements
sum1 = 0.0d0
sum2 = 0.0d0
do 10 i = 1,n-1
sum1 = sum1 + vec(i,col1)
sum2 = sum2 + vec(i,col2)
10 continue
t = n
z(1) = -1/sqrt(t) * sum1
y(1) = -1/sqrt(t) * sum2
x = -1/(t + sqrt(t))
do 20 i = 1,n-1
z(i+1) = x * sum1 + vec(i,col1)
y(i+1) = x * sum2 + vec(i,col2)
20 continue
c now multiply by A to get linearization
i = 1
call matvec(nmax, n, i, nzl, icol, irow, a, ad, z, w)
do 25 i = 1,n
25 z( i) = w( i)
i = 1
call matvec( nmax,n, i, nzl, icol, irow, a, ad, y, w)
do 26 i = 1,n
26 y( i) = w( i)
c sum all weights and divide by 2
sum = 0.0
do 65 i = 1, ilast(n)
sum = sum + weight( i)
65 continue
sum = sum/2.0
c Now iterate, over all values of gamma
do 30 i = 1,n,5
gamma = atan(-z(i)/y(i))
do 40 j = 1, n
w(j) = cos(gamma) * z(j) + sin(gamma) * y(j)
40 continue
c find cost of bisection for w, use part and iadr as integer workspace
call fndcst(n, nzs, part, iadr, w, cost,
* ifirst, ilast, weight, iadj, ip, iq)
if ( cost. gt. best) then
best = cost
do 35 j = 1, n2
ipkeep( j) = ip( j)
iqkeep( j) = iq( j)
35 continue
endif
c update best solution found by eigenvectors
if ( cost .gt. bstev) then
bstev = cost
print 820, ' new champion:', bstev
endif
30 continue
c try Kernighan-Lin
c set up part
print 820, 'now try Ker-Lin:'
do 50 i = 1, n2
part(ipkeep(i)) = i
part(iqkeep(i)) = -i
50 continue
call partit(n2, nn, not0, .false., ipkeep, iqkeep, kp, kq,
* tcost, wk1, wk2, wk3, iwk4, iwk5,
* ifirst, ilast, weight, iadj, part, gain)
c output the original cost, the improved cost, and the difference
c update the best kerlin-solution
if (sum - tcost .gt. bstkl) then
bstkl = sum - tcost
do 60 i = 1, n2
ibest( i) = ipkeep( i)
60 jbest( i) = iqkeep( i)
print 820, ' after Ker-Lin:',bstkl
820 format(1x,a20,f10.2, 3x, 3f10.2)
endif
return
end
subroutine fndcst(n, nzs, part, index, w, cost,
* ifirst, ilast, weight, iadj, ip, iq)
integer ip(n/2), iq(n/2), part(n),
* index(n), ifirst(n), ilast(n), iadj(nzs*2)
double precision sum, cost1, w( n), weight(nzs*2), cost
c
c sort the unprojected vector, but remember the original positions
c of the elements
call sort(n,w,index)
c now round w into a feasible solution and form part
n2 = n/2
do 30 i = 1,n2
w(index(i)) = 0
part(index(i)) = -i
30 continue
do 40 i = n2+1, n
w(index(i)) = 1
part(index(i)) = i-n2
40 continue
c what is the cost of this partition?
c first set up ip and iq
j = 0
k = 0
do 50 i = 1,n
if (w(i) .eq. 1) then
j = j + 1
ip(j) = i
else
k = k + 1
iq(k) = i
endif
50 continue
c now find cost
c sum all weights and divide by 2
sum = 0.0
do 65 i = 1, ilast(n)
sum = sum + weight( i)
65 continue
sum = sum/2.0
cost1 = 0
do 60 i = 1,n2
do 70 j = ifirst(ip(i)),ilast(ip(i))
if (w(iadj(j)) .eq. 0) cost1 = cost1 + weight(j)
70 continue
60 continue
cost = sum - cost1
return
end
subroutine trafo( n, nzl, nzs, icol, irow, a, ideg, ifirst,
* ilast, iadj, weight, iadr)
integer n, nzl, nzs, icol( 1), irow( 1), ideg( 1)
integer ifirst( 1), ilast( 1), iadj( 1), iadr( 1)
double precision a( 1), weight( 1)
c set up ideg
do 10 i = 1,nzl
10 ideg(i) = icol(i)
do 20 i = nzl + 1, n
20 ideg(i) = 0
do 30 l = 1,nzs
i=irow(l)
30 ideg(i) = ideg(i) + 1
ifirst(1) = 1
ilast(1) = ideg(1)
do 40 i = 2, n
ifirst(i) = ifirst(i-1) + ideg(i-1)
40 ilast(i) = ilast(i-1) + ideg(i)
do 50 i = 1, n
50 iadr(i) = ifirst(i)
llast = 0
do 70 j = 1,nzl
if(icol(j) .eq. 0) goto 70
lfirst = llast + 1
llast = llast + icol(j)
do 60 l = lfirst, llast
i = irow(l)
iadj(iadr(i)) = j
weight(iadr(i)) = a(l)
iadr(i) = iadr(i) + 1
iadj(iadr(j)) = i
weight(iadr(j)) = a(l)
iadr(j) = iadr(j) + 1
60 continue
70 continue
return
end
subroutine sort(n,a,ipoint)
c heapsort - sort into nondecreasing order
c taken from "Combinatorial Heuristic Algorithms with FORTRAN" - H.T. Lau
integer ipoint(n)
double precision a(n),atemp
c
do 10 i = 1,n
ipoint(i) = i
10 continue
j1 = n
j2 = n/2
j3 = j2
atemp = a(j2)
jpont = ipoint(j2)
c
20 j4 = j2 + j2
if (j4 .le. j1) then
if (j4 .lt. j1) then
if (a(j4+1) .ge. a(j4)) j4 = j4 + 1
endif
if (atemp .lt. a(j4)) then
a(j2) = a(j4)
ipoint(j2) = ipoint(j4)
j2 = j4
goto 20
endif
endif
c
a(j2) = atemp
ipoint(j2) = jpont
if (j3 .gt. 1) then
j3 = j3 - 1
atemp = a(j3)
jpont = ipoint(j3)
j2 = j3
goto 20
endif
c
if (j1 .ge. 2) then
atemp = a(j1)
jpont = ipoint(j1)
a(j1) = a(1)
ipoint(j1) = ipoint(1)
j1 = j1 - 1
j2 = j3
goto 20
endif
return
end
subroutine partit( n, nn, not0, init, ip, iq, kp, kq,
* tcost, wk1, wk2, wk3, iwk4, iwk5,
* ifirst, ilast, weight, iadj, part, gain)
c the Kernighan-Lin graph bisection heuristic - takes any feasible
c bisection and iteratively improves it
integer ip( n), iq( n), kp( n), kq( n)
integer ifirst( nn*2), ilast( nn*2), iadj( not0), part( nn*2)
double precision weight( not0), wk1( n), wk2( n), wk3( n),
* gain( n), tcost, tmax, small, tot1
logical iwk4( n), iwk5( n), init
if ( init) then
c Set up an initial partition
c Put the first half of the nodes into ip and the second half into iq
do 10 i = 1,n
ip( i) = i
iq( i) = i+n
part( i) = i
part( n+i) = -i
10 continue
endif
c Initialise work arrays
20 do 30 i = 1,n
iwk4( i) = .true.
iwk5( i) = .true.
30 continue
tcost = 0.0
c External cost of a node defined to be total weight of cut edges
c from that node; similarly Internal cost is total weight of UNcut
c edges from that node. Let D(i) be the difference between the external
c cost and the internal cost for node i.
c
c Step 1: Calculate D(i) for every node i. Store in wk1 (for
c nodes in subset ip) and wk2 (for nodes in iq).
c
c Consider all nodes in subset ip
do 70 i = 1,n
tot1 = 0.0d0
tot2 = 0.0
do 80 j = ifirst( ip( i)), ilast( ip( i))
if (part( iadj( j)) .gt. 0) then
c tot2 is total weight uncut edges from node ip(i)
tot2 = tot2 + weight( j)
else
c tot1 is total weight cut edges from node ip(i)
tot1 = tot1 + weight( j)
endif
80 continue
c tcost is total weight of cut edges
tcost = tcost + tot1
wk1( i) = tot1 - tot2
70 continue
small = -2.0 * tcost
c Consider all nodes in subset iq
do 100 i = 1,n
tot1 = 0.0d0
tot2 = 0.0
do 110 j = ifirst( iq( i)), ilast( iq( i))
if (part( iadj( j)) .gt. 0) then
c tot1 is total weight cut edges from node iq(i)
tot1 = tot1 + weight( j)
else
c tot2 is total weight uncut edges from node iq(i)
tot2 = tot2 + weight( j)
endif
110 continue
wk2( i) = tot1 - tot2
100 continue
c If two nodes a and b on opposite sides of the partition are
c interchanged, the gain (reduction in cost) is
c D(a) + D(b) - 2 * the weight on the edge from a to b
c
c Step 2: Choose a and b so that the gain is maximised. Then
c remove a and b from consideration and repeat the process - do
c this until there are no nodes left to consider.
c
c Do until no nodes left to consider...
c Slight modification here: we consider swaps of at most iswap nodes,
c because this is much faster
iswap = 20
do 140 i = 1,iswap
c Initialise the maximum gain
tmax = -99999.9
c Calculate gain for nodes in ip
do 120 j = 1,n
c Is this node still available?
if ( iwk4( j)) then
c Initialise the gain array
do 115 k = 1,n
gain( k) = 0
115 continue
c
c Consider edges from node ip( j)
do 117 l = ifirst( ip( j)), ilast( ip( j))
c icol is the node that this edge connects to
icol = iadj( l)
c index gives the position of the node in the partition
index = part( icol)
c Is this node still available?
if (iwk5( abs(index)) ) then
c Is it on the opposite side of the partition?
if ( index .lt. 0)
c contribution to gain = -2 * weight on edge
* gain( -index) = -2.0 * weight( l)
endif
117 continue
c Calculate the other components of gain
do 118 k = 1,n
c Is this node still available?
if ( iwk5( k)) then
c contribution is D(a) + D(b)
gain( k) = gain( k) + wk1( j) + wk2( k)
c
if (gain( k) .gt. tmax) then
tmax = gain( k)
c ia and ib are the nodes to be interchanged (i.e. best so far)
ia = ip( j)
ib = iq( k)
ind1 = j
ind2 = k
endif
endif
118 continue
endif
120 continue
c Store max gain in wk3, and pair to be swapped in kp and kq
wk3( i ) = tmax
kp( i) = ia
kq( i) = ib
c Remove from further consideration
iwk4( ind1) = .false.
iwk5( ind2) = .false.
c Now recalculate the D values for the partitions with nodes a and b removed
c D(x)'= D(x) +2*(weight on edge from x to a) -2*(weight on edge from x to b)
c D(y)'= D(y) +2*(weight on edge from y to b) -2*(weight on edge from y to a)
c where x and a are in ip subset and y and b are in iq subset
c
c First consider edges from node a
do 130 j = ifirst( ia),ilast( ia)
icol = iadj( j)
index = part( icol)
c If on the ip side of the partition
if (index .gt. 0) then
c If allowed to consider
if (iwk4( index) )
c Add 2*(weight on edge from x to a) to D(x)
* wk1( index) = wk1( index) + 2.0*weight( j)
else
c If allowed to consider
if (iwk5( -index))
c Add -2*(weight on edge from y to a) to D(y)
* wk2( -index) = wk2( -index) - 2.0*weight( j)
endif
130 continue
c Now consider edges from node b
do 135 j = ifirst( ib),ilast( ib)
icol = iadj( j)
index = part( icol)
c If on the ip side of the partition
if (index .gt. 0) then
c If allowed to consider
if (iwk4( index))
c Add -2*(weight on edge from x to b) to D(x)
* wk1( index) = wk1( index) - 2.0*weight( j)
else
c If allowed to consider
if (iwk5( -index))
c Add 2*(weight on edge from y to b) to D(y)
* wk2( -index) = wk2( -index) + 2.0*weight( j)
endif
135 continue
140 continue
c Now that we consider swaps of at most 30 nodes, we need to set kp and
c kq for the remaining n-30 nodes
i = iswap + 1
j = i
do 150 k = 1,n
if (iwk4( k)) then
kp( i) = ip( k)
i = i + 1
endif
if (iwk5( k)) then
kq( j) = iq( k)
j = j + 1
endif
150 continue
c If we swap k nodes, the total gain is gain(1) + ... + gain(k)
c Choose k to maximise the partial sum gain(1) + ... + gain(k)
tmax = small
tot1 = 0.0d0
c Now considering swaps of at most 30 nodes
do 160 i = 1,iswap
tot1 = tot1 + wk3( i)
if ( tot1 .gt. tmax) then
tmax = tot1
k = i
endif
160 continue
c If the total gain is positive, do the swap. If it is (near) zero,
c then a local optimum has been reached and Ker-Lin terminates
if (tmax .gt. 5.0d-05) then
do 170 i = 1,k
c Update partition - do the swap
part( kq( i)) = i
part( kp( i)) = -i
ip( i) = kq( i)
170 iq( i) = kp( i)
k1 = k + 1
do 180 i = k1, n
part( kp( i)) = i
part( kq( i)) = -i
ip( i) = kp( i)
180 iq( i) = kq( i)
c Continue iterating
goto 20
endif
return
end
c -------------- compfg.for --------------
subroutine compfg( n, d, f, g)
c ---------- needs to be adapted ----- (3) ----
parameter (nmax = 5000, ntot = 30000, nmval = 6)
double precision f, d( nmax), g( nmax)
c d = diagonal shift (A - diag(d))
c f = upper bound on graph bisection (f = (1/4)*lmax(VtA(d)V) + suma/4)
c d = subgradient
c local variables
double precision suma, ad( nmax), a( ntot), sum
c variables for dnlaso
double precision val(nmval,4), worklz(800+21*nmax)
double precision qs( nmax, 45), vec( nmax, nmval)
integer ind( nmval)
c constants defining the projection matrix Vn
double precision xconst, yconst
c for temporary storage of eigenvalues
real eigval(5), xbst
external op, iovect
common / store / qs
common / graph1 / a, ad, suma
common / const / xconst, yconst
common / iconst / matmul, nfig, nblock, nval, nperm
common / stop / istop, mulply
common / lzwork / worklz
common / eigvls / eigval, xbst
common / eij / vec
istop = 0
c constants for dnlaso
nmvec = nmax
maxop = n/2
maxj = 42
nperm = 0
c do 10 i=1,n
c10 work( i) = 1.0d0
c set current perturbation of main diag (ad( i) = -d( i))
do 20 i=1,n
20 ad( i) = -d( i)
n1 = n - 1
22 continue
call dnlaso(op, iovect, n1, nval, nfig, nperm, nmval, val,
1 nmvec, vec, nblock, maxop, maxj, worklz, ind, ierr)
do 23 i=1,nperm
23 eigval( i) = val( i, 1)
do 24 i=nperm+1,5
24 eigval( i) = 0.0
mulply = ind( 1)
if ( nperm .le. 0) then
c set istop to 1 to abort BT
if ( nblock .eq. nmval ) istop = 1
nblock = nblock + 1
print 99,' **blocksize increased**:', nblock
matmul = matmul + ind( 1)
if (istop .eq. 0) goto 22
endif
matmul = matmul + ind( 1)
99 format(1x,a25,3i4)
c the upper bound f can now be computed
f = n * val( 1, 1) * 0.25 + suma * 0.25
c calculate the gradient
c first we need the sum of the components of the eigenvector
sum = 0.0d0
do 30 i = 1,n-1
30 sum = sum + vec( i, 1)
g( 1) = - yconst * yconst * sum * sum * n * 0.25
do 40 i = 1, n-1
g( i+1) = xconst * sum + vec( i, 1)
40 g( i+1) = - g( i+1) * g( i+1) * n * 0.25
c make sure that sum(g) = 0
sum = 0.0d0
do 50 i=1,n
50 sum = sum + g( i)
do 60 i=1,n
60 g( i) = g( i) - sum / float( n)
c store eigenvector for next iteration
do 70 i = 1,n1
70 worklz( i) = vec( i, 1)
return
end
subroutine op( n, m, p, q)
c computes q = ap where a is adjacency matrix ( in condensed form0
c p,q are n by m matrices
c needed for dnlaso
c constants related to the graph: nmax, ntot
c ---------- needs to be adapted -------- (4) --
parameter (nmax = 5000, ntot = 30000, nmval = 6)
double precision p( n, m), q( n, m)
c arrays defining the graph
double precision ad( nmax), a( ntot)
integer icol( nmax), irow( ntot)
c constants defining transformation matrix vn and work arrays
double precision xconst, yconst, sum, vec( nmax, nmval)
double precision v2( nmax, nmval), sump( 10), suma
common / graph1 / a, ad, suma
common / graph2 / nzl, icol, irow
common / const / xconst, yconst
c we have to compute (Vnt * A * Vn) * p
c first comes vec := Vn * p
c we need the column sum of the entries in p
np1 = n+1
do 20 j = 1,m
sum = 0.0
do 10 i = 1,n
10 sum = sum + p( i,j)
20 sump( j) = sum
c vec(1) = yconst * sump, vec(i+1) = xconst * sump + p( i)
do 30 j=1,m
do 30 i= 1,n
30 vec( i+1, j) = xconst * sump( j) + p( i, j)
do 40 j = 1,m
40 vec( 1, j) = yconst * sump( j)
c now we form v2 = A * vec
call matvec( nmax,np1, m, nzl, icol, irow, a, ad, vec, v2)
c finally we need Vnt *v2
c we need the column sums of v2 but without element in first row
do 60 j = 1,m
sum = 0.0
do 50 i = 2,np1
50 sum = sum + v2( i, j)
60 sump( j) = sum
c q( i) = yconst * v2( 1) + xconst * sum + v2( i+1)
do 80 j = 1,m
sum = yconst * v2( 1, j) + xconst * sump( j)
do 70 i = 1, n
70 q( i, j) = sum + v2( i+1, j)
80 continue
return
end
subroutine matvec( nmax, n, p, nzl, icol, irow, a, ad, w, u)
c ----------- needs to be adapted --- (5) ---
parameter ( ntot = 30000)
c sparse matrix - vektor multiplikation: u = Aw
c u und w sind n x p matrizen
integer nmax, n, p, nzl, icol( nmax), irow( ntot)
double precision a( ntot), ad( nmax), w( nmax, p), u( nmax, p)
integer i, j, k, l, lfirst, llast
do 10 i = 1,n
do 10 k = 1,p
u( i, k) = ad( i) * w( i, k)
10 continue
c columnwise
llast = 0
do 30 j = 1,nzl
if (icol( j) .eq. 0) goto 30
lfirst = llast + 1
llast = llast + icol( j)
do 20 l = lfirst, llast
i = irow( l)
do 20 k = 1, p
u( i, k) = u( i, k) + a( l) * w( j, k)
u( j, k) = u( j, k) + a( l) * w( i, k)
20 continue
30 continue
return
end
SUBROUTINE IOVECT(N,M,Q,J,K)
C
C THIS SUBROUTINE STORES AND RECALLS VECTORS. IT STORES THE VECTORS
C IN CORE ALTHOUGH MOST REAL PROBLEMS WOULD USE A DISK.
c -------- needs to be adapted -(6)-----
parameter (nmax = 5000)
DOUBLE PRECISION Q(N,M),QS( nmax,45)
COMMON /STORE/QS
IF(K.EQ.1)GO TO 30
DO 20 L=1,M
L1=J-M+L
DO 10 I=1,N
QS(I,L1)=Q(I,L)
10 CONTINUE
20 CONTINUE
RETURN
C
30 DO 50 L=1,M
L1=J-M+L
DO 40 I=1,N
Q(I,L)=QS(I,L1)
40 CONTINUE
50 CONTINUE
RETURN
END