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