c======================================================================= c Test Driver: This is a test drive for using the Quadratic Assignment c Test Problem generators GEN1(subroutine QGE0L1) and GEN2 c (subroutine QGE0L2). The maximum size is limited to 100. You c may change the constants in the driver and the generators to c suit your needs (but keep them the same). c======================================================================= Program Driver Integer MaxN Parameter(MaxN = 100) Integer n, MaxA, MaxB, type, pmin(MaxN) Real A(MaxN, MaxN), B(MaxN, MaxN), rate, delta, fmin Print*, 'Please enter the value of n' Read(*,*) n Print*, 'Please enter the values of MaxA, MaxB, rate, delta' Read(*,*) MaxA, MaxB, rate, delta Call QGE0L1( A, B, n, MaxA, MaxB, rate, delta, pmin, fmin, 1 ) Call U0RMW( A, MaxN, n, 6 ) Call U0RMW( B, MaxN, n, 6 ) Call U0IVW( pmin, n, 6 ) Print*, 'fmin = ', fmin Call QGE0L2( A, B, n, MaxA, MaxB, fmin, pmin, 1, 0 ) Call U0RMW( A, MaxN, n, 6 ) Call U0RMW( B, MaxN, n, 6 ) Call U0IVW( pmin, n, 6 ) Print*, 'fmin = ', fmin Stop End c====================================================================== c Date: July 11, 1991 c Subroutine: QGE0L1 c Purpose: generate quadratic assignment test problems with known c optimal solution. Details of the algorithm see Li and c Pardalos, COAP, 1992 c Input: n, integer, the dimension of the matrix to be generated c MaxA, integer, the maximum entry in matrix A c MaxB, integer, the maximum entry in matrix B c rate, real, the value used to compute the cutoff c delta, real, scale for flow update c type, integer, the type of qap matrix to be generated c type = 1, symmetric qap c type = 2, asymmetric qap c type = 3, grid qap c Output: A, real, flow matrix c B, real, distance matrix c pmin, integer, optimal permutation c fmin, real, optimal objective function value c c modification: 7/23. Make A symmetric if type = 1, 3 c 7/27, demodify the previous one, since noninteger values c result in wrong answer, by Yong Li c 7/27, Yong Li, the parameter Maxn in calling Qu0Fun is c undefined, should be IA c c====================================================================== SUBROUTINE QGE0L1(A,B,n,MaxA,MaxB,rate,delta,pmin,fmin, type) Integer IA Parameter( IA = 100 ) REAL A(IA, n), B(IA,n), rate, delta, fmin, work(IA*IA) INTEGER ia,n,MaxA,MaxB,pmin(n),sign(IA),type,tree(IA, 2) REAL cutoff, mswtri, QU0Fun, Amin INTEGER i, j, row, col REAL RNUNF EXTERNAL RNSET, RNUNF c c create constant matrix A and a sym. matrix B, both with pos entries c CALL U0GCIM( A, IA, n, MaxA ) IF ( type .EQ. 1 ) THEN CALL U0SRIM( B, IA, n, MaxB, Work ) ELSE IF ( type .EQ. 2 ) THEN CALL U0GRIM( B, IA, n, MaxB, Work ) ELSE CALL U0Fact( n, row, col ) CALL UQ0Grd( B, IA, row, col ) END IF c c update flow matrix using min/max spanning tree of the distance matrix c CALL RNSET(0) j = abs( RNUNF()*(delta-1) ) + 1 CALL U0MST( B, IA, n, tree, -1, work ) DO 20 i = 1, n-1 A(tree(i,1), tree(i,2)) = A(tree(i,1), tree(i,2)) + j 20 CONTINUE CALL U0MST( B, IA, n, tree, 1, work ) DO 30 i = 1, n-1 A(tree(i,1), tree(i,2)) = A(tree(i,1), tree(i,2)) - j 30 CONTINUE c c prepare the sign of the triangle c sign(1) = 1 sign(2) = 1 sign(3) = -1 c c update the flow matrix A using signed triangles, spanning trees c CALL QG0CUT( B, IA, n, rate, sign, cutoff, mswtri, work ) call qg0utr( A, B, IA, n, cutoff, mswtri, delta, sign ) c c make the smallest offdiagonal entry of A zero c Amin = A(1,2) DO 32 i = 1, n DO 32 J = 1, n IF ( (i.NE.j).ANd.(A(i,j).lt.Amin) ) THEN Amin = a(i,j) END IF 32 CONTINUE DO 35 i = 1, n DO 35 J = 1, n IF (i.NE.j) THEN A(i,j) = A(i,j) - Amin END IF 35 CONTINUE c c generate optimal permutation randomly and store in sign c permute the flow matrix by pmin, and make pmin the inverse of it c call U0RPer( n, sign ) CALL U0PerM(A, ia, n, sign, work) CALL U0InvP(n, sign, pmin) c c compute the objective function value c fmin = QU0Fun(A, B, IA, n, pmin) RETURN END c====================================================================== c Date: Aug. 10, 1991 c Subroutine: QGE0L2 c Purpose: generate quadratic assignment test problems with known c optimal solution, see Li and Pardalos, COAP, 92 c Input: n, integer, the dimension of the matrix to be generated c MaxA, integer, the maximum entry in matrix A c MaxB, integer, the maximum entry in matrix B c dtype, integer, the type of dist matrix to be generated c dtype = 1, symmetric qap c dtype = 2, asymmetric qap c dtype = 3, grid qap c ftype, integer, the type of flow matrix to be generated c ftype = 0, flow matrix with known min and glb c ftype = 1, flow matrix with bad glb and known max c Output: A, B, real n*n flow and dist matrix c p, the optimal permutation ( min, or max ) c val, real, optimal value ( min or max ) c====================================================================== SUBROUTINE QGE0L2(A,B,n,MaxA,MaxB,val,p,dtype,ftype) Integer MaxN Parameter( MaxN = 100 ) REAL A(MaxN, MaxN), B(MaxN,MaxN), val INTEGER n, MaxA, MaxB, p(n), dtype, ftype INTEGER i, j, total, index, row, col + , pos(MaxN*MaxN), invpos(MaxN*MaxN) EXTERNAL SVRGP, SVRGN, RNUN REAL QU0Fun, work(MaxN*MaxN), delta(MaxN*MaxN) c c randomly generate the distance matrix B, and constant matrix A c CALL U0GCIM( A, MaxN, n, REAL(MaxA) ) IF ( dtype .EQ. 1 ) THEN CALL U0SRIM( B, MaxN, n, MaxB, Work ) ELSE IF ( dtype .EQ. 2 ) THEN CALL U0GRIM( B, MaxN, n, MaxB, Work ) ELSE CALL U0Fact( n, row, col ) CALL UQ0Grd( B, MaxN, row, col ) END IF c c store offdiag entries of B in vector work c total = 1 DO 10 i = 1, n DO 10 j = 1, n IF ( i.NE.j ) THEN work( total ) = B(i,j) total = total + 1 END IF 10 CONTINUE total = total - 1 c c if type = 1, reverse sign of x so it is sorted descendingly c IF ( ftype .EQ. 1 ) THEN DO 15 i = 1, total work(i) = -work(i) 15 CONTINUE END IF CALL U0SetI( invpos, total ) CALL SVRGP( total, work, delta, invpos ) CALL U0InvP( total, invpos, pos ) c c generate a set of n*(n-1) numbers with integer value in (1, MaxA) c CALL RNUN( total, work ) DO 30 i = 1, total work(i) = MOD( INT( work(i)*MaxA ), MaxA ) + 1 30 CONTINUE CALL SVRGN( total, work, delta ) DO 100 i = 1, n DO 100 j = 1, n IF ( i.NE.j ) THEN index = (i-1)*(n-1) + j IF ( i .LE. j ) index = index-1 A(i,j) = A(i,j) - delta( pos(index) ) END IF 100 CONTINUE c c generate optimal permutation randomly and store in sign c permute the flow matrix by pmin, and make pmin the inverse of it c CALL U0SetI( p, n ) c c compute the objective function value c val = QU0Fun(A, B, MaxN, n, p) RETURN END c======================================================================= c Subroutine: GA0Ran c Purpose: to generate a random number c======================================================================= REAL FUNCTION GA0RAN() REAL RNUNF EXTERNAL RNSET, RNUNF INTEGER seed, base SAVE /seed/ DATA seed, base/2*1234567/ c c the initial random number is set by clock c IF ( seed .EQ. base ) THEN CALL RNSET( 0 ) ELSE CALL RNSET( seed ) END IF GA0Ran = RNUNF() seed = INT( base*GA0Ran ) RETURN END c====================================================================== c date: June 26, 1991 c Subroutine: U0ZDIA c Purpose: zero the diagonal of a given matrix c c Input: a, real, n*n matrix c n, integer, the dimension of the matrix to be generated c c Output: a, real, the updated matrix with zeros on diagonal c====================================================================== SUBROUTINE U0ZDIA( a, IA, n ) REAL a(IA, n) INTEGER ia, n, i DO 10 i = 1, n a(i,i) = 0.0 10 CONTINUE RETURN END c====================================================================== c Date: June 26, 1991 c Subroutine: QG0CUT c Purpose: find the cutoff value and the maximum signed weight over c over all triangles c c Input: B, real, the distance matrix of qap c n, integer, the dimension of the matrix to be generated c ia, integer, the dimension of the matrix A when declared c rate, real, the rate specifying the cutoff value c Sign, intege vector (len=3) storing the signs of triangl c c Output: cutoff, real, the cutoff value c mswtri, real, the maximum signed weight of triangles c====================================================================== SUBROUTINE QG0CUT( b, ia, n, rate, sign, cutoff, mswtri, work ) real macmin, macmax parameter ( macmax = 1.0e+20, macmin = 1.0e-20 ) REAL b(ia, n), rate, cutoff, mswtri, work(n,2), temp INTEGER ia, n, i, j, k, sign(3) real rnunf external rnset, rnunf, svrgn c c generate n triangles associated with a pair of randomly generated c facilities, compute their signed weights, and cutoff c call rnset( 0 ) i = mod(int( rnunf()*n ), n) + 1 j = mod(i*I, n) + 1 c c make sure i and j are distinct c if (I .eq. J) then J = mod(J+1, n) + 1 end if c c compute the n signed weights of the triangles (i,j,k),k=1,n c note: if k = i, j, we set the signed weight to infinity c do 10 k = 1, n work(k, 1) = b(k,i)*sign(1) + b(k,j)*sign(2) + b(i,j)*sign(3) 10 continue work(i,1) = macmax work(j,1) = macmax c c sort the signed weights and compute cutoff c call svrgn( n, work(1,1), work(1,2) ) i = int( rate*N ) IF ( i .EQ. 0 ) THEN cutoff = work( 1, 2 ) ELSE IF ( i .GT. n ) THEN cutoff = work( n, 2 ) ELSE cutoff = work( i, 2 ) END IF c c compute mswtri c mswtri = 0 DO 50 K = 1, n DO 40 j = 1, n DO 30 i = 1, n temp = b(k,i)*sign(1) + b(k,j)*sign(2) + b(i,j)*sign(3) c c make sure i,j,k are distinct c if (((k.ne.I).and.(k.ne.j).and.(i.ne.j)).and. + (mswtri.lt.temp)) then mswtri = temp end if 30 CONTINUE 40 CONTINUE 50 CONTINUE RETURN END c====================================================================== c Date: June 26, 1991 c Subroutine: QG0UTR c Purpose: update the flow matrix using signed triangle c ( note the updated amounts are kept to be integers ) c ( also the diagonal entries are unchanged, zero ) c c Input: a, real, the flow matrix of qap c B, real, the distance matrix of qap c ia, integer, the dimension of the matrix A when declared c n, integer, the dimension of the matrix to be generated c cutoff, real, the rate specifying the cutoff value c delta, real, the weight for update c Sign, integer, sign vector of length 3 for triangls c c Output: a, real, the updated flow matrix c====================================================================== SUBROUTINE QG0UTR( a, b, ia, n, cutoff, mswtri, delta, sign ) REAL a(IA, n), b(ia, n), cutoff, mswtri, delta, temp, range INTEGER ia, n, i, j, k, sign(3) range = abs(mswtri) + 1 DO 50 K = 1, n DO 40 j = 1, n DO 30 i = 1, n temp = b(k,i)*sign(1) + b(k,j)*sign(2) + b(i,j)*sign(3) c c make sure i,j,k are distinct, and the signed weight under cutoff c if (((k.ne.I).and.(k.ne.j).and.(i.ne.j)).and. + (temp.lt.cutoff)) then c c make the update in flow to be of integer value c temp = int( (cutoff-temp)*delta/range ) IF ( temp .NE. 0.0 ) THEN a(k,i) = a(k,i) + temp*sign(1) a(k,j) = a(k,j) + temp*sign(2) a(i,j) = a(i,j) + temp*sign(3) END IF c PRINT*, 'QG0UTR: 5' end if 30 CONTINUE 40 CONTINUE 50 CONTINUE c c make sure that the entries of the flow matrix is non negative c temp = a(1,n) DO 80 j = 1, n DO 60 i = 1, n if ((i .ne. j).and.(temp .gt. a(i,j))) then temp = a(i,j) end if 60 CONTINUE 80 CONTINUE c c make the smallest entries in the flow matrix to be zero c if ( temp .ne. 0 ) then DO 95 j = 1, n DO 90 i = 1, n if (i .ne. j) then a(i,j) = a(i,j) - temp end if 90 CONTINUE 95 CONTINUE end if RETURN END c===================================================================== c Date: July 1, 1991 c c Surtoutine: QU0FUN c Purpose: objective function value c c Input: A, real, the flow matrix c B, real, the distance matrix c IA, integer, the leading dimension of A and B c n, integer, the size of the matrix A and B c p, the permutation c c Output: the objective function value of the qap with respect to P c===================================================================== FUNCTION QU0FUN(A, B, IA, N, p) integer IA, n, p(n) REAL QU0Fun, A(IA, n), B(IA, n), objval REAL f, d, epsil PARAMETER( epsil = 1.0E-6 ) objval = 0.0 do 30 j = 1, n do 20 i = 1, n f = A(i,j) d = B(P(i), P(j)) IF ( (abs(f).GE.epsil).AND.(abs(d).GE.epsil) ) THEN objval = objval + f*d END IF 20 continue 30 continue QU0Fun = objval RETURN END c======================================================================= c RNSET: to initialize the random number generator U0RANU c======================================================================= SUBROUTINE RNSET (ISEED) INTEGER ISEED REAL r, U0RANU r = U0RANU(ISEED,0) RETURN END c======================================================================= C RNUN: to generate a set of random numbers uniform in (0,1) c======================================================================= SUBROUTINE RNUN (NR, R) INTEGER NR REAL R(*), RNUNF INTEGER I Do 10 i = 1, NR R(i) = RNUNF() 10 Continue Return End c======================================================================= c Generate a random number. Note, if the random number generate U0RANU c was not set before this subroutine, then the seed 1234567 is used. c======================================================================= REAL FUNCTION RNUNF() Real U0RANU Integer seed Data seed/1234567/ RNUNF = U0RANU(seed, 1) RETURN END c======================================================================= c SVRGN: to sort a real vector RA in ascendingly and store in RB c======================================================================= SUBROUTINE SVRGN (N, RA, RB) Integer n, i, j, idx Real RA(*), RB(*), temp Do 10 i = 1, n RB(i) = RA(i) 10 Continue Do 200 i = 1, n idx = i Do 150 j = i+1, n If (RB(idx) .GT. RB(j)) idx = j 150 Continue temp = RB(i) RB(i) = RB(idx) RB(idx) = temp 200 Continue Return End c======================================================================= c Name: SVRGP c Purpose: to sort a real vector A of size n ascendingly c Input: RA - the vector to be sorted c n - the length of A c Output: RB - sorted vector c IPERM-permutation giving the order c======================================================================= SUBROUTINE SVRGP (N, RA, RB, IPERM) INTEGER N, IPERM(*), i, j, idx, itemp REAL RA(*), RB(*), temp Do 10 i = 1, n IPERM(i) = i RB(i) = RA(i) 10 Continue Do 200 i = 1, n idx = i Do 150 j = i+1, n If (RB(idx) .GT. RB(j)) idx = j 150 Continue temp = RB(i) RB(i) = RB(idx) RB(idx) = temp itemp = IPERM(i) IPERM(i) = IPERM(idx) IPERM(idx) = itemp 200 Continue Return End c====================================================================== c Date: July 21, 1991 c Subroutine: U0FACT c Purpose: factorize an integer into product of its two largest c factors c Input: n, integer, positive c Output: p, q, integer such that n = p*q c====================================================================== SUBROUTINE U0FACT( n, p, q ) INTEGER n, p, q, i p = 1 q = INT( SQRT(float(n)) ) DO 10 i = 1, q IF ( MOD(n, i).EQ.0 ) p = i 10 CONTINUE q = n/p RETURN END c====================================================================== c Date: June 25, 1991 c Subroutine: U0GCIM c Purpose: to generate an integer matrix with constant entries c and zero diagonal c Input: n, integer, the dimension of the matrix to be generated c ia, integer, the dimension of the matrix A when declared c k, real, constant c Output: A, real, flow matrix to be generated c====================================================================== SUBROUTINE U0GCIM( A, ia, n, k ) REAL A( ia, n ), k INTEGER ia, n, i, j DO 20 j = 1, n DO 10 i = 1, n A(i,j) = k 10 CONTINUE 20 CONTINUE c c zero the diagonal c call u0zdia( A, ia, n ) RETURN END c====================================================================== c Date: June 25, 1991 c Subroutine: U0GRIM c Purpose: to generate a random integer matrix with positive entries c Input: n, integer, the dimension of the matrix to be generated c ia, integer, the dimension of the matrix A when declared c MaxA, integer, maximal allowable entries in A c Work, real, working vector of size at least n*n c Output: A, real, flow matrix to be generated c====================================================================== SUBROUTINE U0GRIM( A, ia, n, MaxA, Work ) REAL A( ia, n ), Work( n*n ) INTEGER ia, n, MaxA, i, j EXTERNAL RNSET, RNUN c c generate pseudorandom number from 0-1 uniform distribution c and use them to yield integers between 0 and MaxA c CALL RNSET(0) CALL RNUN( n*n, Work ) DO 20 j = 1, n DO 10 i = 1, n A(i,j) = MOD( INT(Work((i-1)*n+j)*MaxA), MaxA ) + 1 10 CONTINUE 20 CONTINUE c c zero the diagonal c call u0zdia( A, ia, n ) RETURN END c===================================================================== c Date: July 1, 1991 c Surtoutine: U0INVP c Purpose: inverse a permutation c Input: P, integer, the permutation to be inverted c n, integer, the size of the matrix c Output: invp, integer, the inverse of p c===================================================================== subroutine U0INVP( n, p, invp ) integer n, p(n), invp(n), i do 10 i = 1, n invp(p(i)) = i 10 continue return end c====================================================================== c Date: July 11, 1991 c Subroutine: U0MST c Purpose: compute the minimum or maximum spanning tree c Input : n, the size of the graph, the vertices are 1, ..., n c A, the distance graph of a weighted graph, if there is no c edge from i to j, A(i,j) should be set to pos. inf. c Sign, the flag for minimum or maximum spanning tree c Sign = -1 minimum c work, integer working space of size 2*n c Work(i,1) to indicate what vertices are in the tree c Work(i,2) to store the verteices closest or farthest c from vertex i c Output: T, a 2*IA array, storing the edge set of the spanning tree c Sign = 1 maximum c Temp : minmax, next: the variable used to identity the next edge to c be added to the tree. c count: keep the number of edges added into the spanning tree c====================================================================== SUBROUTINE U0MST( A, IA, n, Tree, sign, work ) REAL A(IA, n) INTEGER IA, n, Tree(IA, 2), sign, work(n, 2) INTEGER EMPTY, FULL, count, minmax, next, i, j Parameter ( EMPTY = 0, FULL = 1 ) c c check sign and reset count, inv, and closet( start from vertex 1 ) c IF ( (sign.NE.1).AND.(sign.NE.-1) ) THEN PRINT*, 'U0MST: sign must be 1 or -1' RETURN END IF count = 0 work(1,1) = EMPTY DO 10 i = 2, n work(i,1) = FULL work(i,2) = 1 10 CONTINUE i = 2 DO WHILE ( i .LE. N ) c c find an edge to be added to the spanning tree c minmax = A(i, work(i,2)) next = i c c while there are more vertex in inv, compute the edge c with mininum or maximum length c DO 20 j = i, n IF (work(j,1).NE.EMPTY) THEN IF ( sign*(A(j, work(j,2)) - minmax) .GT. 0 ) THEN minmax = A(j, work(j,2)) next = j END IF END IF 20 CONTINUE c c add the edge (next, mate(next)) to the tree, and delete next from inv c don't forget to update count c count = count + 1 tree(count,1) = next tree(count,2) = work(next, 2) work(next, 1) = EMPTY c c modify the array mate(v) for i in inv c DO 100 j = 1, n IF (work(j,1).NE.EMPTY) THEN IF ( sign*( A(j,next) - A(j,work(j,2))) .GT. 0 ) THEN work(j,2) = next END IF END IF 100 CONTINUE i = 1 DO WHILE ( (i .LE. n) .AND. (work(i,1).EQ.EMPTY) ) i = i + 1 END DO END DO RETURN END c===================================================================== c Date: June 27, 1991 c Surtoutine: U0PERM c Purpose: permute a matrix A with a permutation P to get c P*A*trans(P) c Input: A, real, matrix of size n*n c ia, integer, the declared dimension of A c n, integer, the size of the matrix c p, integer, permutation c work, real, working space of size n*n c Output: A, real, the permuted matrix c===================================================================== subroutine U0PERM( A, ia, n, p, work ) integer ia, n, p(n), i, j REAL A(IA, n), work(IA, n) DO 20 j = 1, n DO 10 i = 1, n work(p(i),j) = A(i, j) 10 CONTINUE 20 CONTINUE DO 40 j = 1, n DO 30 i = 1, n A(i,p(j)) = work(i, j) 30 CONTINUE 40 CONTINUE RETURN END c===================================================================== c Date: June 27, 1991 c c Surtoutine: U0RPER c Purpose: to generate a random permutation of size n c c Input: n, integer, the size of the permutation ( input ) c Output: p, integer, the permutation ( output ) c===================================================================== subroutine U0RPER( n, p ) integer n, p(n), i, itemp, nm1, m real RNUNF, GA0RAN external RNSET, RNUNF, GA0RAN c c initialize variables, set permutation to the identity permutation c m = n nm1 = n-1 DO 5 i=1,n p(i) = i 5 CONTINUE c c randomly permute p c DO 10 i=1,nm1 j = 1 + MOD( INT(GA0RAN()*m), n) itemp = p(m) p(m) = p(j) p(j) = itemp m = m-1 10 CONTINUE RETURN END c====================================================================== c Date: Aug. 10, 1991 c Subroutine: U0SETI c Purpose: set the integer vector to be the identity permutation c Input: n, integer, the dimension of the matrix to be generated c Output: p, the identity permutation of length n c====================================================================== SUBROUTINE U0SetI( p, n ) INTEGER n, p(n), i DO 10 i = 1, n p(i) = i 10 CONTINUE RETURN END c====================================================================== c Date: July 8, 1991 c Subroutine: U0SRIM c Purpose: generate a random symmetric pos matrix with integer entry c Input: n, integer, the dimension of the matrix to be generated c ia, integer, the dimension of the matrix A when declared c MaxA, integer, maximal allowable entries in A c Work, real, working vector of size at least n*n c Output: A, real, flow matrix to be generated c====================================================================== SUBROUTINE U0SRIM( A, ia, n, MaxA, Work ) REAL A( ia, n ), Work( n*n ) INTEGER ia, n, MaxA, i, j EXTERNAL RNSET, RNUN c c generate pseudorandom number from 0-1 uniform distribution c and use them to yield integers between 0 and MaxA c CALL RNSET(0) CALL RNUN( n*n, Work ) DO 20 j = 1, n DO 10 i = 1, n A(i,j) = MOD( INT( ( Work((i-1)*n+j) + + Work((j-1)*n+i) )*MaxA/2.0 ), MaxA ) + 1 10 CONTINUE 20 CONTINUE c c zero the diagonal c call u0zdia( A, ia, n ) RETURN END c====================================================================== c Date: July 21, 1991 c Subroutine: UQ0GRD c Purpose: generate the grid distance matrix c Input: IA, integer, the dimension of the matrix A when declared c row, integer, the number of rows of the grid c col, integer, the number of cols of the grid c Output: A, n*n real distance for the grid c====================================================================== SUBROUTINE UQ0GRD( A, IA, row, col ) REAL A(IA, 1) INTEGER IA, n, row, col, numcol, i, j, ir, ic, jr, jc IF ( col .EQ. 0 ) THEN PRINT*, 'QU0Grd: error, col = 0' RETURN END IF n = row*col numcol = col DO 10 i = 1, n DO 10 j = 1, n IF ( i.EQ.j ) THEN A(i,j) = 0.0 ELSE ir = (i-1)/numcol + 1 ic = mod( i-1, numcol ) + 1 jr = (j-1)/numcol + 1 jc = mod( j-1, numcol ) + 1 A(i,j) = float( ABS(ir-jr) + ABS(ic-jc) ) END IF 10 CONTINUE RETURN END c====================================================================== c Date: 3/16/1992 c Name: U0RANU c Purpose: to generate a random real number of uniform distribution c in (0,1) c Input: seed - the seed to generate the random number, if set = 0 or c m = 0; ignore the value of seed otherwise. c SOURCE: Fortran 77 For Engineers and Scientists, by c L.R. Nyhoff & S. Leestma, Macmillan, NY. 1988. pp. 369 c modified by Yong Li c====================================================================== REAL FUNCTION U0RANU( seed, set ) Integer seed,set,m,const1 Real const2 Parameter (const1 = 2147483647, const2 = .4656613E-9) Save Data m /0/ c If (m.EQ.0) m = seed If ((set.EQ.0).OR.(m.EQ.0)) m = seed m = m * 65539 If (m .LT. 0) m = (m + 1) + const1 U0RANU = m * const2 Return End c====================================================================== c Date: July 8, 1991 c Subroutine: U0RMW c Purpose: write the data of a real matrix in the following form c to an output file specified by a unit number c a11 a12 ... a1n c ...... c an1 an2 ... ann c c Input: ia, integer, the dimension of the matrix A when declared c n, integer, the dimension of the matrix to be generated c Wfile, the file unit to which the matrix is written c A, real, an n*n matrix c c Output: NONE c====================================================================== SUBROUTINE U0RMW( A, IA, n, Wfile ) REAL A(IA, n) INTEGER IA, n, Wfile INTEGER i, j DO 20 i = 1, n WRITE(Wfile, 5) ( A(i,j), j = 1, n ) 5 FORMAT( 6(2x, f10.1) ) 20 CONTINUE RETURN END c===================================================================== c Date: July 8, 1991 c Surtoutine: U0IVW c Purpose: to write an integer vector to an output file c Input: n, integer, the size of the permutation ( input ) c Rfile, the unit number of the input file c v, integer, vector to read c Output: NONE c===================================================================== subroutine U0IVW( v, n, Wfile ) integer n, v(n), Wfile, i WRITE(Wfile, 5) ( v(i), i = 1, n ) 5 FORMAT( 10(2x, i5) ) RETURN END