program main !*****************************************************************************80 ! !! tsp_anneal_test() tests tsp_anneal(). ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 20 June 2026 ! ! Author: ! ! John Burkardt. ! implicit none integer, parameter :: rk8 = kind ( 1.0D+00 ) integer, parameter :: n = 42 real ( kind = rk8 ) cost integer i integer i2 integer j integer order(n) character ( len = 9 ) :: prefix = 'dantzig42' real ( kind = rk8 ) x(n,2) real ( kind = rk8 ), allocatable :: x2(:,:) call timestamp ( ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'tsp_anneal_test():' write ( *, '(a)' ) ' Fortran90 version' write ( *, '(a)' ) ' tsp_anneal() uses simulated annealing ' write ( *, '(a)' ) ' to solve the traveling salesman problem (TSP).' ! ! Assign coordinates to the cities. ! x = transpose ( reshape ( (/ & 170.0, 85.0, & 166.0, 88.0, & 133.0, 73.0, & 140.0, 70.0, & 142.0, 55.0, & 126.0, 53.0, & 125.0, 60.0, & 119.0, 68.0, & 117.0, 74.0, & 99.0, 83.0, & 73.0, 79.0, & 72.0, 91.0, & 37.0, 94.0, & 6.0, 106.0, & 3.0, 97.0, & 21.0, 82.0, & 33.0, 67.0, & 4.0, 66.0, & 3.0, 42.0, & 27.0, 33.0, & 52.0, 41.0, & 57.0, 59.0, & 58.0, 66.0, & 88.0, 65.0, & 99.0, 67.0, & 95.0, 55.0, & 89.0, 55.0, & 83.0, 38.0, & 85.0, 25.0, & 104.0, 35.0, & 112.0, 37.0, & 112.0, 24.0, & 113.0, 13.0, & 125.0, 30.0, & 135.0, 32.0, & 147.0, 18.0, & 147.0, 36.0, & 154.0, 45.0, & 157.0, 54.0, & 158.0, 61.0, & 172.0, 82.0, & 174.0, 87.0 /), (/ 2, 42 /) ) ) ! ! Apply simulated annealing. ! call tsp_anneal ( n, x, cost, order ) ! ! Report the result. ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' The final path:' write ( *, '(a)' ) ' ' write ( *, '(a,t13,a,t23,a)') 'City','X','Y' do i = 1, n j = order(i) write ( *, '(i4,2f10.4)' ) j, x(j,1), x(j,2) end do write ( *, '(a)' ) '' write ( *, '(a,g14.6)' ) ' Cost of path = ', cost ! ! Create the circuit data with transposition, reordering and wrap-around. ! allocate ( x2 ( 1 : n, 1 : 2 ) ) do i = 1, n i2 = order(i) x2(i,1:2) = x(i2,1:2) end do ! ! Send the data to the plotter. ! call tsp_display ( prefix, n, x2 ) deallocate ( x2 ) ! ! Terminate. ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'tsp_anneal_test():' write ( *, '(a)' ) ' Normal end of execution.' write ( *, '(a)' ) '' call timestamp ( ) stop end subroutine timestamp ( ) !*****************************************************************************80 ! !! timestamp() prints the current YMDHMS date as a time stamp. ! ! Example: ! ! 31 May 2001 9:45:54.872 AM ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 15 August 2021 ! ! Author: ! ! John Burkardt ! implicit none character ( len = 8 ) ampm integer d integer h integer m integer mm character ( len = 9 ), parameter, dimension(12) :: month = (/ & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' /) integer n integer s integer values(8) integer y call date_and_time ( values = values ) y = values(1) m = values(2) d = values(3) h = values(5) n = values(6) s = values(7) mm = values(8) if ( h < 12 ) then ampm = 'AM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Noon' else ampm = 'PM' end if else h = h - 12 if ( h < 12 ) then ampm = 'PM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Midnight' else ampm = 'AM' end if end if end if write ( *, '(i2.2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm ) return end subroutine tsp_display ( prefix, n, x2 ) !*****************************************************************************80 ! !! tsp_display() plots cities and a tour connecting them. ! ! Discussion: ! ! The code automatically includes a final stop at the first city. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 27 July 2026 ! ! Author: ! ! John Burkardt ! ! Input: ! ! character ( len = * ) prefix: a string that defines the name of ! the problem. It is also used to create file names needed for ! input to gnuplot. ! ! integer n: the number of cities. ! ! real ( kind = rk8 ) x(n,2): the x and y coordinates of the cities. ! implicit none integer, parameter :: rk8 = kind ( 1.0D+00 ) integer n character ( len = 80 ) command_filename integer command_unit character ( len = 80 ) data_filename integer data_unit integer i character ( len = * ) prefix real ( kind = rk8 ) x2(n,2) ! ! Create the data file. ! call get_unit ( data_unit ) data_filename = prefix // '_data.txt' open ( unit = data_unit, file = data_filename, status = 'replace' ) do i = 1, n write ( data_unit, '(2x,g14.6,2x,g14.6)' ) x2(i,1), x2(i,2) end do i = 1 write ( data_unit, '(2x,g14.6,2x,g14.6)' ) x2(i,1), x2(i,2) close ( unit = data_unit ) write ( *, '(a)' ) ' Created data file "' // trim ( data_filename ) // '".' ! ! Create command file. ! call get_unit ( command_unit ) command_filename = prefix // '_commands.txt' open ( unit = command_unit, file = command_filename, status = 'replace' ) write ( command_unit, '(a)' ) '# ' // trim ( command_filename ) write ( command_unit, '(a)' ) '#' write ( command_unit, '(a)' ) '# Usage:' write ( command_unit, '(a)' ) '# gnuplot < ' // trim ( command_filename ) write ( command_unit, '(a)' ) '#' write ( command_unit, '(a)' ) 'set term png' write ( command_unit, '(a)' ) 'set output "' // prefix // '.png"' write ( command_unit, '(a)' ) 'set xlabel "<-- X -->"' write ( command_unit, '(a)' ) 'set ylabel "<-- Y -->"' write ( command_unit, '(a)' ) 'set title "' // prefix // '"' write ( command_unit, '(a)' ) 'set grid' write ( command_unit, '(a)' ) 'unset key' write ( command_unit, '(a)' ) 'set style data lines' write ( command_unit, '(a)' ) 'plot "' // trim ( data_filename ) // & '" using 1:2 lw 3 linecolor rgb "blue", \' write ( command_unit, '(a)' ) ' "' // trim ( data_filename ) // & '" using 1:2 with points pointtype 7 pointsize 3 linecolor rgb "green"' write ( command_unit, '(a)' ) 'quit' close ( unit = command_unit ) write ( *, '(a)' ) & ' Created command file "' // trim ( command_filename ) // '".' return end subroutine get_unit ( iunit ) !*****************************************************************************80 ! !! get_unit() returns a free Fortran unit number. ! ! Discussion: ! ! A "free" Fortran unit number is a value between 1 and 99 which ! is not currently associated with an I/O device. A free Fortran unit ! number is needed in order to open a file with the OPEN command. ! ! If IUNIT = 0, then no free Fortran unit could be found, although ! all 99 units were checked (except for units 5, 6 and 9, which ! are commonly reserved for console I/O). ! ! Otherwise, IUNIT is a value between 1 and 99, representing a ! free Fortran unit. Note that GET_UNIT assumes that units 5 and 6 ! are special, and will never return those values. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 26 October 2008 ! ! Author: ! ! John Burkardt ! ! Output: ! ! integer IUNIT, the free unit number. ! implicit none integer i integer ios integer iunit logical lopen iunit = 0 do i = 1, 99 if ( i /= 5 .and. i /= 6 .and. i /= 9 ) then inquire ( unit = i, opened = lopen, iostat = ios ) if ( ios == 0 ) then if ( .not. lopen ) then iunit = i return end if end if end if end do return end