C Maximization of determinant of n*n integer matrix with elements C 1..n^2; C attempts to find next terms in C http://www.research.att.com/~njas/sequences/A085000 C Author: Hugo Pfoertner http://www.pfoertner.org/ C C URL of this file: C http://www.randomwalk.de/sequences/a085000.txt C C Recent record result found with this program C C 2007 Oct 28 C d7:= det mat ( C (46, 42, 15, 2, 27, 24, 18), C ( 9, 48, 36, 30, 7, 14, 31), C (39, 11, 44, 34, 13, 29, 5), C (26, 22, 17, 41, 47, 1, 21), C (20, 8, 40, 6, 33, 23, 45), C ( 4, 28, 19, 25, 38, 49, 12), C (32, 16, 3, 37, 10, 35, 43) ); C d7 := 762150368499 C C The following results were produced by a more sophisticated search method C developed by Hermann Jurksch C 2008 Jan 6, private communication C d8 := det mat ( C ( 1, 34, 55, 19, 59, 26, 45, 20), C (23, 2, 18, 61, 35, 57, 39, 25), C (21, 56, 3, 27, 47, 41, 12, 53), C (50, 42, 31, 4, 17, 60, 43, 13), C (24, 33, 64, 40, 5, 38, 8, 48), C (49, 6, 30, 14, 36, 16, 46, 63), C (29, 54, 22, 52, 11, 7, 58, 28), C (62, 32, 37, 44, 51, 15, 9, 10) ); C d8 := 440943507851753 C C 2007 Oct 1 found by Hugo Pfoertner and Hermann Jurksch C d9:= det mat ( C (68, 7, 12, 62, 73, 26, 29, 58, 34), C (67, 37, 43, 10, 3, 61, 33, 78, 36), C (30, 20, 79, 53, 49, 71, 40, 25, 2), C (56, 50, 8, 27, 42, 60, 81, 4, 41), C (23, 14, 54, 63, 11, 18, 72, 44, 70), C ( 1, 38, 32, 21, 65, 66, 22, 48, 76), C (45, 74, 31, 80, 17, 46, 5, 24, 47), C (15, 77, 35, 39, 51, 16, 59, 69, 9), C (64, 52, 75, 13, 57, 6, 28, 19, 55) ) C d9 := 3.46254605664224e+17 C C 2008 Jan 5, found by Herman Jurksch C d10:= det mat ( C ( 1, 2, 61, 84, 81, 82, 39, 54, 41, 60), C (53, 57, 3, 65, 94, 20, 91, 22, 66, 33), C (46, 63, 47, 4, 45, 78, 83, 28, 13, 98), C (79, 42, 49, 71, 5, 95, 51, 10, 77, 26), C (17, 75, 87, 58, 30, 6, 38, 27, 86, 80), C (68, 93, 76, 50, 85, 56, 7, 37, 14, 19), C (100,16, 31, 35, 62, 34, 8, 64, 67, 88), C (21, 72, 29, 9, 48, 73, 43, 97, 89, 25), C (52, 70, 23, 96, 11, 36, 55, 92, 12, 59), C (69, 15, 99, 32, 44, 24, 90, 74, 40, 18) ) C d10 = -3.56944784623E+020 C C Version history: (result updates see above) C 2008 Jan 22 Improved results for all n>=7 included as comments C 2007 Aug 21 Including the information received from Hermann Jurksch: C Improved result for n=8 and unpublished upper bound for C determinant value C 2007 Jun 1 C preprocessor-like alternatives for target system dependencies C 2007 May 29 Counter for random restarts C 04.12.06 Generic program for arbitrary n, include parameter statement C from file a85kdim.for and create name of output file C a085000_nn.txt internally (nn replaced by actual dimension) C 21.11.06 Append good results to external result file, C use default target when searching for maximal determinant C 26.07.06 Search for target value of determinant C 12.08.04 Use F90 random generator RANDOM_NUMBER C 23.09.03 Search by exchange steps C 19.09.03 5*5 matrix included C 17.09.03 Initial version C doubleprecision d, db, dglob, ai, aj, dtargt, good C Include external file with problem dimension and cut-off value C for result inclusion into permanent output file C Examples for contents of file a85kdim.for: C parameter ( n=8, good=4.40876D14 ) C parameter ( n=10, good=3.5690D20 ) C etc. include 'a85kdim.for' parameter ( nn=n*n ) C...+....1....+....2....+....3....+....4....+....5....+....6....+....7.. doubleprecision A(n,N), b(n,N), det(2), work(n), as(nn), bs(nn) & , dmue real x(nn) integer ip(nn), ipvt(nn), seed(2) character datstr*8, timstr*10, outfil*14 C Alternative RNG: Mersenne twister CTwister integer genrand_int32 C equivalence (a,as) equivalence (b,bs) CTwister external genrand_int32 C C Build name of permanent result file if ( n .lt. 10 ) then outfil = 'a085000_x.txt' write (outfil(9:9),'(I1)') n else outfil = 'a085000_xx.txt' write (outfil(9:10),'(I2.2)') n endif C C Get starting time of computation from system T0 = SECNDS ( 0.0 ) C Initial values for seed generation SEED(1) = NINT ( 100.0 * T0 ) SEED(2) = 1 C Twister IF ( MOD ( IRAN, 2 ) .EQ. 0 ) IRAN = IRAN + 1 C Twister j = init_genrand(iran) CALL RANDOM_SEED ( PUT=SEED(1:2) ) C C Target replaced by Hermann Jurksch's upper bound C (private communication, 2007 Aug 19) C lambda := sqrt(n^2*(n+1)*(n^2+1)/12) C mu := (n^2+1)/2 C D_max := n * my * lambda^(n-1) C dmue = 0.5D0 * dble(nn+1) C use variable dtargt temporarily as Lambda dtargt = sqrt ( dble(nn) * dble(n+1) * dmue / 6.0D0 ) dtargt = dble(n) * dmue * dtargt**(n-1) C dglob = 2.0D0*dtargt C Ask for target value (activate if required) C write (*,*) ' Target value' C read (*,*) dtargt write (*,*) ' Target value:', dtargt C C Counter for number of restarts nrstrt = 0 500 continue C C increment counter nrstrt = nrstrt + 1 c Generate random vector using the Mersenne Twister C or equivalent RNG call RANDOM_NUMBER ( HARVEST=X(1:NN) ) CTwister do 10 i = 1, nn CTwister x(i) = float ( genrand_int32 () ) CTwister 10 continue c Permutation vector by sorting (SLATEC) C Source available via C http://gams.nist.gov/serve.cgi/Module/SLATEC/SPSORT/11543/ call spsort ( x, nn, ip, -1, ier ) c set matrix elements do 20 i = 1, nn as(ip(i)) = dble(i) bs(ip(i)) = dble(i) 20 continue C Determinant of new random matrix C Gauss elimination (Linpack) C#ifdef LINPC C old linpack C call dgefa ( B, n, n, ipvt, info ) C#else C Lapack call dgetrf ( n, n, B, n, ipvt, info ) C#endif C Determinant (Linpack) call dgedi ( B, n, n, ipvt, det, work, 10 ) d = det(1) * 10.0d0**det(2) C C Counter for successful interchanges l = 0 C Perform sweeps of interchange attempts do 100 it = 1, 50*n ltset = 0 do 30 i = 1, nn-1 do 40 j = i+1, nn C Copy to work matrix B, because Gauss elimination C overwrites matrix. DO 41 k = 1, nn bs(k) = as(k) 41 continue c Interchange matrix elements i and j in 1-dim stored bS(j) = aS(i) ai = as(i) bS(i) = aS(j) aj = as(j) C b (bs) are destroyed by dgefa or dgetrf C #ifdef LINPC C old linpack C call dgefa ( b, n, n, ipvt, info ) C#else C Lapack call dgetrf ( n, n, B, n, ipvt, info ) C#endif call dgedi ( b, n, n, ipvt, det, work, 10 ) db = det(1) * 10.0d0**det(2) C better approximation to target? if ( abs(abs(db)-dtargt) .lt. abs(abs(d)-dtargt) ) then l = l + 1 lt = it ltset = 1 CDiagnosis write (*,*) l, i, j, db C perform successful interchange in original matrix aS(j) = ai aS(i) = aj d = db endif 40 continue 30 continue C Stop when last sweep went through without update if ( ltset .eq. 0 ) goto 110 100 continue 110 continue C Update best value and print if ( abs(abs(d)-dtargt) .LE. abs(dglob-dtargt) ) then dglob = abs(d) write (*,1002) d, l, lt, nrstrt 1002 format ( ' d=', 1P, D25.16, & ' l=', I6, ' lt=', I4, ' nrstrt=', I11 ) if ( n .le. 10 ) then write (*,1000) (nint(as(k)),k=1,nn) C...+....1....+....2....+....3....+....4....+....5....+....6....+....7.. 1000 format ( (1X,25i3),:,/,(1X,25i3),:,/,(1X,25i3),:, & /,(1X,25i3),:,/,(1X,25i3) ) else write (*,1010) (nint(as(k)),k=1,nn) C...+....1....+....2....+....3....+....4....+....5....+....6....+....7.. 1010 format ( (1X,19i4),:,/,(1X,19i4),:,/,(1X,19i4),:, & /,(1X,19i4),:,/,(1X,19i4),/,(1X,19i4),:,/,(1X,19i4),:, & /,(1X,19i4),:,/,(1X,19i4),/,(1X,19i4),:,/,(1X,19i4) ) endif C if ( dglob .gt. good ) then C Append good results to global logfile open ( unit=10, file=outfil, form='formatted', & status='old', position='append', iostat=ios ) if ( ios .eq. 0 ) then call date_and_time ( date=datstr, time=timstr ) write (10,1001) datstr(1:4), datstr(5:6), datstr(7:8), & timstr(1:2), timstr(3:4), timstr(5:6) 1001 format ( /, 1X, A, '/', A, '/', A, 2X, A, ':', A, ':', A ) write (10,*) ' d=', d if ( n .lt. 10 ) then write (10,1000) (nint(as(k)),k=1,nn) else write (10,1010) (nint(as(k)),k=1,nn) endif endif close ( unit=10 ) endif endif C Skip back to generation of next random matrix C (infinite loop) goto 500 end