program main !*****************************************************************************80 ! !! MAIN is the main program for MONTE_CARLO_RULE. ! ! Discussion: ! ! MONTE_CARLO_RULE generates N points in the M-dimensional unit hypercube, ! and writes out files so that the data can be regarded as a quadrature rule. ! ! Usage: ! ! monte_carlo_rule m n seed ! ! where ! ! * M, the spatial dimension, ! * N, the number of points to generate, ! * SEED, the seed, a positive integer. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 23 December 2011 ! ! 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 character ( len = 3 ) m_string integer ( kind = 4 ) n character ( len = 6 ) n_string real ( kind = 8 ), allocatable, dimension ( :, : ) :: r character ( len = 255 ) r_filename integer ( kind = 4 ) s character ( len = 9 ) s_string integer ( kind = 4 ) seed character ( len = 255 ) string real ( kind = 8 ), allocatable, dimension ( : ) :: w character ( len = 255 ) w_filename real ( kind = 8 ), allocatable, dimension ( :, : ) :: x character ( len = 255 ) x_filename call timestamp ( ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'MONTE_CARLO_RULE' write ( *, '(a)' ) ' FORTRAN90 version' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Compute the abscissas and weights of a quadrature rule' write ( *, '(a)' ) ' that is simply a Monte Carlo sampling.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' The program requests input values from the user:' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' * M, the spatial dimension,' write ( *, '(a)' ) ' * N, the number of points to generate,' write ( *, '(a)' ) ' * SEED, a positive integer.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Output from the program includes' write ( *, '(a)' ) ' a set of 3 files that define the quadrature rule.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' (1) "mc_m?_n?_s?_r.txt", the ranges;' write ( *, '(a)' ) ' (2) "mc_m?_n?_s?_w.txt", the weights;' write ( *, '(a)' ) ' (3) "mc_m?_n?_s?_x.txt", the abscissas.' ! ! 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, S ! if ( 3 <= arg_num ) then iarg = 3 call getarg ( iarg, string ) call s_to_i4 ( string, s, ierror, last ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Enter the seed S (1 or greater)' read ( *, * ) s end if write ( *, '(a,i12)' ) ' Seed S = ', s ! ! Allocate the arrays. ! allocate ( r(1:m,1:2) ) allocate ( w(1:n) ) allocate ( x(1:m,1:n) ) ! ! Compute the data. ! r(1:m,1) = 0.0D+00 r(1:m,2) = +1.0D+00 w(1:n) = 1.0D+00 / real ( n, kind = 8 ) seed = s call r8mat_uniform_01 ( m, n, seed, x ) ! ! Construct appropriate file names. ! write ( n_string, '(i6)' ) n write ( m_string, '(i3)' ) m write ( s_string, '(i9)' ) s r_filename = 'mc_d' // m_string // '_n' // n_string & // '_s' // s_string // '_r.txt' w_filename = 'mc_d' // m_string // '_n' // n_string & // '_s' // s_string // '_w.txt' x_filename = 'mc_d' // m_string // '_n' // n_string & // '_s' // s_string // '_x.txt' call s_blank_delete ( r_filename ) call s_blank_delete ( w_filename ) call s_blank_delete ( x_filename ) ! ! Write the rule to files. ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Creating R file = "' // trim ( r_filename ) // '".' call r8mat_write ( r_filename, m, 2, r ) write ( *, '(a)' ) ' Creating W file = "' // trim ( w_filename ) // '".' call r8mat_write ( w_filename, 1, n, w ) write ( *, '(a)' ) ' Creating X file = "' // trim ( x_filename ) // '".' call r8mat_write ( x_filename, m, n, x ) ! ! Free memory. ! deallocate ( r ) deallocate ( w ) deallocate ( x ) ! ! Terminate. ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'MONTE_CARLO_RULE' write ( *, '(a)' ) ' Normal end of execution.' write ( *, '(a)' ) ' ' call timestamp ( ) stop 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 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 subroutine r8mat_uniform_01 ( m, n, seed, r ) !*****************************************************************************80 ! !! R8MAT_UNIFORM_01 returns a unit pseudorandom R8MAT. ! ! Discussion: ! ! An R8MAT is an array of R8's. ! ! For now, the input quantity SEED is an integer variable. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 31 May 2007 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Paul Bratley, Bennett Fox, Linus Schrage, ! A Guide to Simulation, ! Second Edition, ! Springer, 1987, ! ISBN: 0387964673, ! LC: QA76.9.C65.B73. ! ! Bennett Fox, ! Algorithm 647: ! Implementation and Relative Efficiency of Quasirandom ! Sequence Generators, ! ACM Transactions on Mathematical Software, ! Volume 12, Number 4, December 1986, pages 362-376. ! ! Pierre L'Ecuyer, ! Random Number Generation, ! in Handbook of Simulation, ! edited by Jerry Banks, ! Wiley, 1998, ! ISBN: 0471134031, ! LC: T57.62.H37. ! ! Peter Lewis, Allen Goodman, James Miller, ! A Pseudo-Random Number Generator for the System/360, ! IBM Systems Journal, ! Volume 8, Number 2, 1969, pages 136-143. ! ! Parameters: ! ! Input, integer ( kind = 4 ) M, N, the number of rows and columns ! in the array. ! ! Input/output, integer ( kind = 4 ) SEED, the "seed" value, which ! should NOT be 0. On output, SEED has been updated. ! ! Output, real ( kind = 8 ) R(M,N), the array of pseudorandom values. ! implicit none integer ( kind = 4 ) m integer ( kind = 4 ) n integer ( kind = 4 ) i integer ( kind = 4 ), parameter :: i4_huge = 2147483647 integer ( kind = 4 ) j integer ( kind = 4 ) k integer ( kind = 4 ) seed real ( kind = 8 ) r(m,n) if ( seed == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8MAT_UNIFORM_01 - Fatal error!' write ( *, '(a)' ) ' Input value of SEED = 0.' stop end if do j = 1, n do i = 1, m k = seed / 127773 seed = 16807 * ( seed - k * 127773 ) - k * 2836 if ( seed < 0 ) then seed = seed + i4_huge end if r(i,j) = real ( seed, kind = 8 ) * 4.656612875D-10 end do end do return end subroutine r8mat_write ( output_filename, m, n, table ) !*****************************************************************************80 ! !! R8MAT_WRITE writes an R8MAT file. ! ! 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 greater precision in the output file, try: ! ! '(', m, 'g', 24, '.', 16, ')' ! write ( string, '(a1,i8,a1,i8,a1,i8,a1)' ) '(', m, 'g', 16, '.', 8, ')' ! ! Write the data. ! do j = 1, n write ( output_unit, string ) table(1:m,j) end do ! ! Close the file. ! close ( unit = output_unit ) return end subroutine s_blank_delete ( s ) !*****************************************************************************80 ! !! S_BLANK_DELETE removes blanks from a string, left justifying the remainder. ! ! Discussion: ! ! All TAB characters are also removed. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 26 July 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, character ( len = * ) S, the string to be transformed. ! implicit none character ch integer ( kind = 4 ) get integer ( kind = 4 ) put character ( len = * ) s integer ( kind = 4 ) s_length character, parameter :: tab = achar ( 9 ) put = 0 s_length = len_trim ( s ) do get = 1, s_length ch = s(get:get) if ( ch /= ' ' .and. ch /= tab ) then put = put + 1 s(put:put) = ch end if end do s(put+1:s_length) = ' ' 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 software is released under the GNU LGPL license. ! ! Modified: ! ! 12 November 2007 ! ! 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 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 == ' ' ) 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 == ' ' ) 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.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