subroutine get_seed ( seed ) !*****************************************************************************80 ! !! GET_SEED returns a seed for the random number generator. ! ! Discussion: ! ! The seed depends on the current time, and ought to be (slightly) ! different every millisecond. Once the seed is obtained, a random ! number generator should be called a few times to further process ! the seed. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 02 August 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer ( kind = 4 ) SEED, a pseudorandom seed value. ! implicit none integer ( kind = 4 ), parameter :: i4_huge = 2147483647 integer ( kind = 4 ) seed real ( kind = 8 ) temp character ( len = 10 ) time character ( len = 8 ) today integer ( kind = 4 ) values(8) character ( len = 5 ) zone call date_and_time ( today, time, zone, values ) temp = 0.0D+00 temp = temp + real ( values(2) - 1, kind = 8 ) / 11.0D+00 temp = temp + real ( values(3) - 1, kind = 8 ) / 30.0D+00 temp = temp + real ( values(5), kind = 8 ) / 23.0D+00 temp = temp + real ( values(6), kind = 8 ) / 59.0D+00 temp = temp + real ( values(7), kind = 8 ) / 59.0D+00 temp = temp + real ( values(8), kind = 8 ) / 999.0D+00 temp = temp / 6.0D+00 do while ( temp <= 0.0D+00 ) temp = temp + 1.0D+00 end do do while ( 1.0D+00 < temp ) temp = temp - 1.0D+00 end do seed = int ( real ( i4_huge, kind = 8 ) * temp ) ! ! Never use a seed of 0 or maximum integer. ! if ( seed == 0 ) then seed = 1 end if if ( seed == i4_huge ) then seed = seed - 1 end if return end subroutine get_unit ( iunit ) !*****************************************************************************80 ! !! GET_UNIT returns a free FORTRAN unit number. ! ! Discussion: ! ! A "free" FORTRAN unit number is an integer 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 an integer 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 GNU LGPL license. ! ! Modified: ! ! 26 October 2008 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer ( kind = 4 ) IUNIT, the free unit number. ! implicit none integer ( kind = 4 ) i integer ( kind = 4 ) ios integer ( kind = 4 ) 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 function i4_uniform ( a, b, seed ) !*****************************************************************************80 ! !! I4_UNIFORM returns a scaled pseudorandom I4. ! ! Discussion: ! ! An I4 is an integer ( kind = 4 ) value. ! ! The pseudorandom number will be scaled to be uniformly distributed ! between A and B. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 12 November 2006 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Paul Bratley, Bennett Fox, Linus Schrage, ! A Guide to Simulation, ! Springer Verlag, pages 201-202, 1983. ! ! Pierre L'Ecuyer, ! Random Number Generation, ! in Handbook of Simulation, ! edited by Jerry Banks, ! Wiley Interscience, page 95, 1998. ! ! Bennett Fox, ! Algorithm 647: ! Implementation and Relative Efficiency of Quasirandom ! Sequence Generators, ! ACM Transactions on Mathematical Software, ! Volume 12, Number 4, pages 362-376, 1986. ! ! Peter Lewis, Allen Goodman, James Miller ! A Pseudo-Random Number Generator for the System/360, ! IBM Systems Journal, ! Volume 8, pages 136-143, 1969. ! ! Parameters: ! ! Input, integer ( kind = 4 ) A, B, the limits of the interval. ! ! Input/output, integer ( kind = 4 ) SEED, the "seed" value, which ! should NOT be 0. On output, SEED has been updated. ! ! Output, integer ( kind = 4 ) I4_UNIFORM, a number between A and B. ! implicit none integer ( kind = 4 ) a integer ( kind = 4 ) b integer ( kind = 4 ) i4_uniform integer ( kind = 4 ) k real ( kind = 4 ) r integer ( kind = 4 ) seed integer ( kind = 4 ) value if ( seed == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4_UNIFORM - Fatal error!' write ( *, '(a)' ) ' Input value of SEED = 0.' stop end if k = seed / 127773 seed = 16807 * ( seed - k * 127773 ) - k * 2836 if ( seed < 0 ) then seed = seed + 2147483647 end if r = real ( seed, kind = 4 ) * 4.656612875E-10 ! ! Scale R to lie between A-0.5 and B+0.5. ! r = ( 1.0E+00 - r ) * ( real ( min ( a, b ), kind = 4 ) - 0.5E+00 ) & + r * ( real ( max ( a, b ), kind = 4 ) + 0.5E+00 ) ! ! Use rounding to convert R to an integer between A and B. ! value = nint ( r, kind = 4 ) value = max ( value, min ( a, b ) ) value = min ( value, max ( a, b ) ) i4_uniform = value return end subroutine latin_edge ( dim_num, point_num, seed, x ) !*****************************************************************************80 ! !! LATIN_EDGE returns edge points in a Latin square. ! ! Discussion: ! ! In each spatial dimension, there will be exactly one ! point with the coordinate value ! ! ( 0, 1, 2, ..., point_num-1 ) / ( point_num - 1 ) ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 10 September 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) DIM_NUM, the spatial dimension. ! ! Input, integer ( kind = 4 ) POINT_NUM, the number of points, which should ! be at least 2! ! ! Input/output, integer ( kind = 4 ) SEED, a seed for the random ! number generator, ! needed if the portable UNIFORM routine is being used. ! ! Output, real ( kind = 8 ) X(DIM_NUM,POINT_NUM), the points. ! implicit none integer ( kind = 4 ) dim_num integer ( kind = 4 ) point_num integer ( kind = 4 ) :: base = 1 integer ( kind = 4 ) i integer ( kind = 4 ) j integer ( kind = 4 ) perm(point_num) integer ( kind = 4 ) seed real ( kind = 8 ) x(dim_num,point_num) if ( point_num == 1 ) then x(1:dim_num,1) = 0.5D+00 else do i = 1, dim_num call perm_uniform ( point_num, base, seed, perm ) do j = 1, point_num x(i,j) = real ( perm(j) - 1, kind = 8 ) & / real ( point_num - 1, kind = 8 ) end do end do end if return end subroutine perm_uniform ( n, base, seed, p ) !*****************************************************************************80 ! !! PERM_UNIFORM selects a random permutation of N objects. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 18 November 2008 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Albert Nijenhuis, Herbert Wilf, ! Combinatorial Algorithms, ! Academic Press, 1978, second edition, ! ISBN 0-12-519260-6. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of objects to be permuted. ! ! Input, integer ( kind = 4 ) BASE, is 0 for a 0-based permutation and 1 for ! a 1-based permutation. ! ! Input/output, integer ( kind = 4 ) SEED, a seed for the random ! number generator. ! ! Output, integer ( kind = 4 ) P(N), the permutation. P(I) is the "new" ! location of the object originally at I. ! implicit none integer ( kind = 4 ) n integer ( kind = 4 ) base integer ( kind = 4 ) i integer ( kind = 4 ) i4_uniform integer ( kind = 4 ) j integer ( kind = 4 ) k integer ( kind = 4 ) p(n) integer ( kind = 4 ) seed do i = 1, n p(i) = ( i - 1 ) + base end do do i = 1, n j = i4_uniform ( i, n, seed ) k = p(i) p(i) = p(j) p(j) = k end do return end subroutine r8mat_write ( output_filename, m, n, table ) !*****************************************************************************80 ! !! R8MAT_WRITE writes an R8MAT file. ! ! Discussion: ! ! An R8MAT is an array of R8 values. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 31 May 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) OUTPUT_FILENAME, the output file name. ! ! Input, integer ( kind = 4 ) M, the spatial dimension. ! ! Input, integer ( kind = 4 ) N, the number of points. ! ! Input, real ( kind = 8 ) TABLE(M,N), the table data. ! implicit none integer ( kind = 4 ) m integer ( kind = 4 ) n integer ( kind = 4 ) j character ( len = * ) output_filename integer ( kind = 4 ) output_status integer ( kind = 4 ) output_unit character ( len = 30 ) string real ( kind = 8 ) table(m,n) ! ! Open the file. ! call get_unit ( output_unit ) open ( unit = output_unit, file = output_filename, & status = 'replace', iostat = output_status ) if ( output_status /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8MAT_WRITE - Fatal error!' write ( *, '(a,i8)' ) ' Could not open the output file "' // & trim ( output_filename ) // '" on unit ', output_unit output_unit = -1 stop end if ! ! Create a format string. ! ! For less precision in the output file, try: ! ! '(', m, 'g', 14, '.', 6, ')' ! if ( 0 < m .and. 0 < n ) then write ( string, '(a1,i8,a1,i8,a1,i8,a1)' ) '(', m, 'g', 24, '.', 16, ')' ! ! Write the data. ! do j = 1, n write ( output_unit, string ) table(1:m,j) end do end if ! ! Close the file. ! close ( unit = output_unit ) return end subroutine random_initialize ( seed_input ) !*****************************************************************************80 ! !! RANDOM_INITIALIZE initializes the FORTRAN90 random number seed. ! ! Discussion: ! ! The G95 compiler is currently generating internal compiler errors ! when it tries to compile this routine. Just one more exasperation ! on the mountain of complications because of the ragged interface with ! the nonstandard random number generator standard! ! ! If you don't initialize the FORTRAN90 random number generator ! routine RANDOM_NUMBER, which is used by calls like ! ! call random_number ( harvest = x ) ! ! then its behavior is not specified. That may be OK for you. But ! you may want to be able to force the same sequence to be generated ! each time, or to force a different sequence. ! ! To control the sequence of random numbers, you need to set the seed. ! In FORTRAN90, this is done by calling the RANDOM+SEED routine. ! You can call it with no arguments, in fact. But if you call ! it with no arguments: ! ! call random_seed ( ) ! ! then its behavior (or more particularly, the behavior of RANDOM_NUMBER) ! is still not specified. You might hope that the system will go out ! and pick a nice random seed for you, but there's no guarantee. ! ! ! For example, on the DEC ALPHA, if you compile a program that calls ! RANDOM_NUMBER, then every time you run it, you get the same sequence ! of "random" values. If you compile a program that calls RANDOM_SEED ! with no arguments, and then calls RANDOM_NUMBER, you still get the ! same sequence each time you run the program. ! ! In order to actually try to scramble up the random number generator ! a bit, this routine goes through the tedious process of getting the ! size of the random number seed, making up values based on the current ! time, and setting the random number seed. ! ! Unfortunately, the RANDOM_SEED routine has a very elastic definition. ! It does not use a single scalar integer SEED. Instead, it communicates ! with the user through an integer vector whose size is not specified. ! You actually have to "ask" the routine to tell you the size of this ! vector. Then you can fill up the vector with values to be used to ! influence the seeding of the random number routine. The details of ! how the seed affects the sequence are also unspecified, but the only ! thing we are reasonably confident about is that specifying the same ! seed should result in the same sequence, and specifying different ! seeds should result in different sequences! ! ! I assume this is far more than you wanted to know. (It's certainly ! more than I wanted to know!) ! ! The argument SEED is an input quantity only, so it is legal ! to type ! ! call random_initialize ( 0 ) ! ! or ! ! call random_initialize ( 18867 ) ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 21 February 2006 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) SEED_INPUT, the user's "suggestion" for a seed. ! However, if the input value is 0, the routine will come up with ! its own "suggestion", based on the system clock. ! implicit none integer ( kind = 4 ) count integer ( kind = 4 ) count_max integer ( kind = 4 ) count_rate logical, parameter :: debug = .false. integer ( kind = 4 ) i integer ( kind = 4 ) seed integer ( kind = 4 ) seed_input integer ( kind = 4 ), allocatable :: seed_vector(:) integer ( kind = 4 ) seed_size real ( kind = 4 ) t integer ( kind = 4 ), parameter :: warm_up = 100 seed = seed_input ! ! Initialize the random seed routine. ! call random_seed ( ) ! ! Determine the size of the random number seed vector. ! call random_seed ( size = seed_size ) ! ! Allocate a vector of the right size to be used as a random seed. ! allocate ( seed_vector(seed_size) ) ! ! If the user supplied a SEED value, use that. ! ! Otherwise, use the system clock value to make up a value that is ! likely to change based on when this routine is called. ! if ( seed /= 0 ) then if ( debug ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'RANDOM_INITIALIZE' write ( *, '(a,i20)' ) ' Initialize RANDOM_NUMBER, user SEED = ', seed end if else call system_clock ( count, count_rate, count_max ) seed = count if ( debug ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'RANDOM_INITIALIZE' write ( *, '(a,i20)' ) ' Initialize RANDOM_NUMBER, arbitrary SEED = ', & seed end if end if ! ! Set the seed vector. We don't know the significance of the ! individual entries of the internal seed vector, so we'll just set ! all entries to SEED. ! seed_vector(1:seed_size) = seed ! ! Now call RANDOM_SEED, and tell it to use this seed vector. ! call random_seed ( put = seed_vector(1:seed_size) ) ! ! Free up the seed space. ! deallocate ( seed_vector ) ! ! Call the random number routine a bunch of times just to "warm it up". ! do i = 1, warm_up call random_number ( harvest = t ) end do return end subroutine s_to_i4 ( s, value, ierror, length ) !*****************************************************************************80 ! !! S_TO_I4 reads an integer value from a string. ! ! Discussion: ! ! Instead of ICHAR, we now use the IACHAR function, which ! guarantees the ASCII collating sequence. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 12 January 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) S, a string to be examined. ! ! Output, integer ( kind = 4 ) VALUE, the integer value read from the string. ! If the string is blank, then VALUE will be returned 0. ! ! Output, integer ( kind = 4 ) IERROR, an error flag. ! 0, no error. ! 1, an error occurred. ! ! Output, integer ( kind = 4 ) LENGTH, the number of characters ! of S used to make the integer. ! implicit none character c integer ( kind = 4 ) i integer ( kind = 4 ) ierror integer ( kind = 4 ) isgn integer ( kind = 4 ) length character ( len = * ) s integer ( kind = 4 ) state character :: TAB = achar ( 9 ) integer ( kind = 4 ) value value = 0 ierror = 0 length = 0 state = 0 isgn = 1 do i = 1, len_trim ( s ) c = s(i:i) ! ! STATE = 0, haven't read anything. ! if ( state == 0 ) then if ( c == ' ' .or. c == TAB ) then else if ( c == '-' ) then state = 1 isgn = -1 else if ( c == '+' ) then state = 1 isgn = +1 else if ( lle ( '0', c ) .and. lle ( c, '9' ) ) then state = 2 value = iachar ( c ) - iachar ( '0' ) else ierror = 1 return end if ! ! STATE = 1, have read the sign, expecting digits or spaces. ! else if ( state == 1 ) then if ( c == ' ' .or. c == TAB ) then else if ( lle ( '0', c ) .and. lle ( c, '9' ) ) then state = 2 value = iachar ( c ) - iachar ( '0' ) else ierror = 1 return end if ! ! STATE = 2, have read at least one digit, expecting more. ! else if ( state == 2 ) then if ( lle ( '0', c ) .and. lle ( c, '9' ) ) then value = 10 * value + iachar ( c ) - iachar ( '0' ) else value = isgn * value ierror = 0 length = i - 1 return end if end if end do ! ! If we read all the characters in the string, see if we're OK. ! if ( state == 2 ) then value = isgn * value ierror = 0 length = len_trim ( s ) else value = 0 ierror = 1 length = 0 end if return 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 GNU LGPL license. ! ! Modified: ! ! 18 May 2013 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! None ! implicit none character ( len = 8 ) ampm integer ( kind = 4 ) d integer ( kind = 4 ) h integer ( kind = 4 ) m integer ( kind = 4 ) mm character ( len = 9 ), parameter, dimension(12) :: month = (/ & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' /) integer ( kind = 4 ) n integer ( kind = 4 ) s integer ( kind = 4 ) values(8) integer ( kind = 4 ) 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,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