program main !*****************************************************************************80 ! !! MAIN is the main program for LATIN_CENTER_DATASET. ! ! Discussion: ! ! LATIN_CENTER_DATASET generates and saves a Latin Center Square dataset. ! ! A "Latin Center Square" dataset is created by dividing each side ! of the unit hypercube into N subintervals, and then choosing ! a set of N subcubes. The subcubes are selected in such a way ! that if we project the subcubes onto any coordinate direction, ! there is exactly one subcube in each subinterval. ! ! That's the "Latin Square" part of the dataset. ! ! Then we choose the center of each subcube as a "representative". ! This set of points constitutes the "Latin Center Square" dataset, ! a terminology I made up to distinguish it from the "Latin Edge Square" ! and "Latin Random Square" datasets. ! ! Usage: ! ! latin_center_dataset m n seed ! ! where ! ! * M, the spatial dimension, ! * N, the number of points to generate, ! * SEED, the seed, a positive integer. ! ! creates an M by N dataset and writes it to the ! file "latin_center_M_N.txt". ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 13 December 2009 ! ! Author: ! ! John Burkardt ! implicit none integer ( kind = 4 ) arg_num integer ( kind = 4 ) iarg integer ( kind = 4 ) iargc integer ( kind = 4 ) ierror integer ( kind = 4 ) ios integer ( kind = 4 ) last integer ( kind = 4 ) m integer ( kind = 4 ) n character ( len = 255 ) output_filename real ( kind = 8 ), allocatable, dimension ( :, : ) :: r integer ( kind = 4 ) seed character ( len = 255 ) string call timestamp ( ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'LATIN_CENTER_DATASET' write ( *, '(a)' ) ' FORTRAN90 version' write ( *, '(a)' ) ' Generate a Latin Center Square dataset.' ! ! Get the number of command line arguments. ! arg_num = iargc ( ) ! ! Get the spatial dimension M. ! if ( 1 <= arg_num ) then iarg = 1 call getarg ( iarg, string ) call s_to_i4 ( string, m, ierror, last ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Enter the spatial dimension M (1 or greater)' read ( *, * ) m end if write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) ' Spatial dimension M = ', m ! ! Get the number of points N. ! if ( 2 <= arg_num ) then iarg = 2 call getarg ( iarg, string ) call s_to_i4 ( string, n, ierror, last ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Enter the number of points N (1 or greater)' read ( *, * ) n end if write ( *, '(a,i8)' ) ' Number of points N = ', n ! ! Get the seed, SEED ! if ( 3 <= arg_num ) then iarg = 3 call getarg ( iarg, string ) call s_to_i4 ( string, seed, ierror, last ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Enter the seed SEED (1 or greater)' read ( *, * ) seed end if write ( *, '(a,i12)' ) ' Seed SEED = ', seed if ( seed == 0 ) then call get_seed ( seed ) write ( *, '(a,i12)' ) ' Randomized SEED = ', seed end if ! ! Compute the data ! allocate ( r(1:m,1:n) ) call latin_center ( m, n, seed, r ) ! ! Write it to a file. ! write ( output_filename, '(a,i2.2,a,i5.5,a)' ) & 'latin_center_', m, '_', n, '.txt' call r8mat_write ( output_filename, m, n, r ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) & ' The Latin Center Square data was written to the file "' & // trim ( output_filename ) // '".' ! ! Free memory. ! deallocate ( r ) ! ! Terminate. ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'LATIN_CENTER_DATASET' write ( *, '(a)' ) ' Normal end of execution.' write ( *, '(a)' ) ' ' call timestamp ( ) stop end 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 ), parameter :: i4_huge = 2147483647 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 + i4_huge 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_center ( dim_num, point_num, seed, x ) !*****************************************************************************80 ! !! LATIN_CENTER returns points in a Latin Center Square. ! ! Discussion: ! ! In each spatial dimension, there will be exactly one ! point with the coordinate value ! ! ( 1, 3, 5, ..., 2*point_num-1 ) / ( 2 * point_num ) ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 20 March 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) DIM_NUM, the spatial dimension. ! ! Input, integer ( kind = 4 ) POINT_NUM, the number of points. ! ! Input/output, integer ( kind = 4 ) SEED, a seed for the random ! number generator. ! ! 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) do i = 1, dim_num call perm_uniform ( point_num, base, seed, perm ) do j = 1, point_num x(i,j) = real ( 2 * perm(j) - 1, kind = 8 ) & / real ( 2 * point_num, kind = 8 ) end do end do 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 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