C Maximization of determinant of n*n integer matrix with elements
C 1..n^2
C
C Hugo Pfoertner http://www.pfoertner.org/
C Version history:
C 23.09.03 Search by exchange steps
C 19.09.03 5*5 matrix included
C 17.09.03 Initial version
      parameter ( n=7, nn=n*n )
      doubleprecision d, db, dglob, ai, aj
C...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
      doubleprecision A(n,N), b(n,N), det(2), work(n), as(nn), bs(nn)
      real x(nn)
      integer ip(nn), ipvt(nn)
      integer genrand_int32
      equivalence (a,as)
      equivalence (b,bs)
      external genrand_int32
C
C Get starting time of computation from system
      T0 = SECNDS ( 0.0 )
C Initial values for seed generation
      IRAN = NINT ( 100.0 * T0 )
      IF ( MOD ( IRAN, 2 ) .EQ. 0 ) IRAN = IRAN + 1
      j = init_genrand(iran)
      dglob = 0
500   continue
c Generate random vector using the Mersenne Twister
      do 10 i = 1, nn
      x(i) = float ( genrand_int32 () )
10    continue
c Permutation vector by sorting (SLATEC)
      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)
      call dgefa ( B, n, n, ipvt, info )
C Determinant (Linpack)
      call dgedi ( B, n, n, ipvt, det, work, 10 )
      d = det(1) * 10.0d0**det(2)
      l = 0
C Perform 10 sweeps of interchange attempts
      do 100 it = 1, 15
      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
      call dgefa ( b, n, n, ipvt, info )
      call dgedi ( b, n, n, ipvt, det, work, 10 )
      db = det(1) * 10.0d0**det(2)
      if ( abs(db) .gt. abs(d) ) then
        l = l + 1
        lt = it
c        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
100   continue
C Update best value and print
      if ( abs(d) .gE. dglob ) then
        dglob = abs(d)
        write (*,*) ' d=', d, ' l, lt:', l, lt
        write (*,1000) (nint(as(k)),k=1,nn)
1000    format ( (25i3),:,/,(25i3) )
      endif
C Skip back to generation of next random matrix
C (infinite loop)
      goto 500
      end