C Search n*n non-negative integer matrix elements that give, when
C permuted in all possible ways (n^2)!/(n!)^2 different determinants.
C The maximum occurring matrix element shall be minimized.
C Currently anything beyond n=3 is out of reach.
C
C Author:
C Hugo Pfoertner http://www.pfoertner.org/
C
C Program tested and run on SGI Altix 3000 (16 processors)
C http://www.sgi.com/products/servers/altix/index.html
C using the Intel Fortran 8 compiler for Linux
C http://www.intel.com/software/products/compilers/clin/
C
C Thanks to MTU Aero Engines http://www.mtu.de/
C for making these powerful resources available
C
C
C Version history:
C Sep 10 2004 Sorting only applied to output
C Sep 9 2004 adjust scaling to current largest matrix element
C Sep 8 2004 draw matrix elements from
C individual triangular distributions
C Aug 25 2004 Exit immediately when first collision is found
C Aug 19 2005 Initial version
C
implicit integer (a-z)
parameter ( n=3, m=5000000, mxm=n*n )
real t0, x(0:mxm), rantri
dimension v(-m:m), p(9,10080), seed(2)
C Determinants calculated using 64bit floating point arithmetic
C gives larger range unless 64 bit integers are available.
C On Intel Itanium 2 processors also a speed advantage is expected.
doubleprecision i1,i2,i3,i4,i5,i6,i7,i8,i9,
& d, d2, d3, b(mxm), bmax(mxm)
external rantri
C Statement functions to compute determinant
D2(I1,I2,I3,I4)=I1*I4-I2*I3
D3(I1,I2,I3,I4,I5,I6,I7,I8,I9)=
& I1*D2(I5,I6,I8,I9)
& -I2*D2(I4,I6,I7,I9)
& +I3*D2(I4,I5,I7,I8)
C...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
C This matrix produces all 10080 determinants, used as prototype
C to determine a set of permutations leading to the maximum number
C of different determinants.
data bmax / 0.0D0, 3.0D0, 11.0D0, 46.0D0, 72.0D0, 103.0D0,
& 109.0D0, 113.0D0, 117.0D0 /
C...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
C
nn = n * n
C generate array of permutation indices by recording the permutations
C that lead to the first occurrence in the hit count list for
C determinant values produced by permuting the "bmax" matrix elements
do 200 j = 1, 9
200 b(j) = dble(j)
do 205 j = -m, m
205 v(j) = 0
kper = 0
C Start of skipped loop over all 9! permutations
210 continue
det = nint(d3 (bmax(1),bmax(2),bmax(3),bmax(4),bmax(5),bmax(6),
& bmax(7),bmax(8),bmax(9)))
v(det) = v(det) + 1
C permutation leading to the first hit is recorded
if ( v(det) .eq. 1 ) then
kper = kper + 1
do 220 j = 1, nn
p(j,kper) = nint(b(j))
220 continue
endif
C permute index array and prototype matrix simultaneously
C lpg is available at http://www.randomwalk.de/sequences/lpg.txt
call lpg ( nn, b, nexb )
call lpg ( nn, bmax, next )
C unless reverse order is found, skip back
if ( next .ne. 0 ) goto 210
C Result 10080 is expected, print for appeasement of user
write (*,*) ' Permutations:', kper
C
C Preset max det and diff # of dets
dmax = 0
mulmax = 0
C
C Ask user for starting value of largest matrix element and
C maximum number of search loop.
C
write (*,*) ' npmax, rep'
read (*,*) npmax, rep
C Get time in seconds since midnight from system. This is used to
C initialize the seed for the random number generator.
C Any other initialization can be used, taking into account that
C parallel jobs should be started with different seeds.
C Luxury version: Use "thread-safe" RNGs for parallel jobs.
T0 = SECNDS ( 0.0 )
C The Ecuyer generator used in Intel Fortran needs two seeds
seed(1) = NINT ( 100.0 * T0 )
seed(2) = rep
C Fortan 90 library function
call random_seed ( put=seed(1:2) )
C
C Preset global target
limopt = 150
C
C Main iteration loop
C
do 100 l = 1, rep
C
C Adjust target every 10000 search steps. It is assumed that a
C file named "globopt.txt" is kept in the directory where the
C executable is started. This file also serves to adjust the
C search limit dynamically in cases for jobs running in parallel
C on multiple CPUs
C
if ( mod(l,10000) .eq. 1 ) then
open (unit=21,file='globopt.txt',status='old',
& iostat=ioret)
if ( ioret .eq. 0 ) then
read (21,*) limopt
if ( npmax .gt. limopt-1 )
& write (*,*) 'NPmax adjusted to:', limopt-1
npmax = limopt - 1
endif
close ( unit=21 )
endif
C Clear occurrence counts for determinant values
do 5 i = -m, m
v(i) = 0
5 continue
C Trace maximum determinant
dmax = 0
C Generate random vector from uniform (0,1) distribution
C using the standard Fortran 90 function
C Any suitable other RNG may be used here, provided its repetition
C period is sufficiently large to avoid duplication of searches
C after long running time. Intel Fortran 8.0 provides a uniform RNG
C with > 10^18 period length.
C See:
C
call random_number ( harvest=x(0:10) )
C
C Generate matrix elements by drawing from triangular
C distributions that were derived from previously computed
C successful designs drawn from uniform distributions.
C see: http://www.randomwalk.de/scimath/m33distr.pdf
C
C Set largest element cyclically to npmax,npmax-1,-2,..-5
b(1) = dble (npmax - mod(l,5))
C set smallest element to 0 in 70% of all cases
if ( x(1) .gt. 0.7 ) then
x(9) = rantri ( x(9), 0.0, 0.1, 0.01 )
else
x(9) = 0.0
endif
x(2) = rantri ( x(2), 0.8, 0.995, 0.98 )
x(3) = rantri ( x(3), 0.72, 0.98, 0.92 )
x(4) = rantri ( x(4), 0.55, 0.96, 0.87 )
x(5) = rantri ( x(5), 0.35, 0.95, 0.8 )
x(6) = rantri ( x(6), 0.15, 0.85, 0.49 )
C x(8) is drawn from 2 uniform distributions:
C 90% (0..0.2), 10%(0.2..0.4) decision based on additional
C random value in x(0)
x(8) = x(8) * 0.2
if ( x(0) .gt. 0.9 ) x(8) = x(8) + 0.2
C x(7) is drawn from a symmetric triangular distribution
C that is dynamically inserted between x(6) and x(8)
C Sort such that x(6)>x(8)
if ( x(6) .lt. x(8) ) then
t0 = x(6)
x(6) = x(8)
x(8) = t0
endif
C write (*,*) x(8), x(6)
x(7) = rantri ( x(7), x(8), x(6), 0.5*(x(8)+x(6)) )
C
C set matrix elements
C Transform random numbers from (0,1) to range 0..b(1)
do 20 i = 2, nn
b(i) = dble ( int( x(i) * (b(1)+1.0D0) ))
20 continue
C
C loop over the permutations covering all possible arrangements that
C produce a full set of distinct determinant values
do 26 j = 1, kper
C
det = nint(d3( b(p(1,j)),b(p(2,j)),b(p(3,j)),b(p(4,j)),b(p(5,j)),
& b(p(6,j)),b(p(7,j)),b(p(8,j)),b(p(9,j))))
C
C Determine maximum of determinant
if ( det .gt. dmax ) then
dmax = det
endif
C Increment hit counts
v(det) = v(det) + 1
C Check for "early termination" of loop over permutations
C Terminate permutation loop when double hit is found
if ( abs(v(det)) .gt. 1 ) then
goto 100
endif
C
C End of permutation loop
26 continue
C
C Count hits
muldet = 0
C The following code is superfluous when "early termination" is active
do 40 k = -dmax, dmax
if ( v(k) .ne. 0 ) muldet = muldet + 1
40 continue
C Report good results (only effective when early termination is
C not activated
if ( muldet .ge. 10076 ) then
mulmax = muldet
do 31 k = 1, nn
31 bmax(k) = b(k)
C
C Sort output. Any suitable sorting routine can be used,
C e.g. SORT from NAPACK (adjusted for double prec.)
C or DSORT from SLATEC
C See e.g. http://gams.nist.gov/serve.cgi/Class/N6a2b/
C
call vsorti ( NN, BMAX )
C
write (*,1000) dmax, mulmax, (nint(bmax(k)),k=nn,1,-1)
1000 format ( i9, i11, 9I4 )
C Get global minimum of npmax and look for improvement
open (unit=21,file='globopt.txt',status='old',
& iostat=ioret)
C Skip global update if open fails
if ( ioret .ne. 0 ) goto 32
read (21,*) limopt
if ( mulmax .eq. 10080 ) then
btop = 0
do 33 j = 1, nn
33 btop = max(nint(bmax(j)),btop)
if ( btop .ge. limopt ) goto 34
rewind 21
write (21,1021) btop, (nint(bmax(k)),k=1,nn)
1021 format ( i4, 9i4 )
34 continue
C
close ( unit=21 )
C Reduce search target
npmax = btop - 1
C
else
C
close ( unit=21 )
if ( npmax .eq. limopt-1 ) goto 100
C If global target is more ambitious than local target then adjust
npmax = limopt - 1
endif
C
mulmax = 0
goto 100
C
32 continue
C if global opt could not be read, restart with reduced target
npmax = npmax - 1
mulmax = 0
goto 100
endif
C
C Target for jumps terminating the processing of the
C current set of matrix elements
C
100 continue
C
C Output of final result (also found in globopt.txt)
write (*,1000) npmax+1, mulmax, (nint(bmax(k)),k=1,nn)
end
C...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
real function rantri ( ran01, a, b, m )
C create a random number from a triangular distribution
C with minimum a, maximum b and 3rd vertex at m
C ran01 is a uniform random variable calculated externally.
C
C Author: Hugo Pfoertner http://www.pfoertner.org/
C
C Change history:
C Sep 8, 2004 include ran01 in parameter list (i.e. the uniform random
C numbers have to be created outside rantri)
C Aug 30, 2004 Initial version
C
real ran01, a, b, m, dab, dam, dmb
C
C...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
dab = b - a
dam = m - a
dmb = b - m
C
C determine if ran01 is in rising or falling part of the pdf
if ( ran01 .lt. dam/dab ) then
rantri = a + sqrt ( dab * dam * ran01 )
else
rantri = b - sqrt ( dab * dmb * (1.0 - ran01 ) )
endif
CTest write (*,*) a,b,m,rantri
C End of function rantri
end