subroutine bezier_patch_evaluate ( node_num, node_xyz, rectangle_num, & rectangle_node, patch, point_num, point_uv, point_xyz ) !*****************************************************************************80 ! !! bezier_patch_evaluate() evaluates a Bezier patch. ! ! Discussion: ! ! Given a Bezier surface defined as a collection of patches, ! it is desired to evaluate the surface, in a single patch, ! for a variety of (U,V) parameter values. ! ! The user specifies the patch index, and the (U,V) coordinates of ! several points. The routine returns the (X,Y,Z) coordinates of the ! points. The special parameter values (0,0), (1,0), (0,1) and (1,1) should ! return the coordinates of the SW, SE, NW and NE corners of the patch. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 September 2021 ! ! Author: ! ! John Burkardt ! ! Input: ! ! integer NODE_NUM, the number of nodes. ! ! integer NODE_XYZ(3,NODE_NUM), the node coordinates. ! ! integer RECTANGLE_NUM, the number of rectangles. ! ! integer RECTANGLE_NODE(16,RECTANGLE_NUM), the nodes that make up ! each rectangle. ! ! integer PATCH, the index of the Bezier patch. ! ! integer POINT_NUM, the number of points at which evaluation ! is desired. ! ! real ( kind = rk ) POINT_UV(2,POINT_NUM), the (U,V) parameter ! coordinates of the points. ! ! Output: ! ! real ( kind = rk ) POINT_XYZ(3,POINT_NUM), the (X,Y,Z) ! physical coordinates of the points. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer node_num integer point_num integer rectangle_num real ( kind = rk ) node_xyz(3,node_num) integer patch real ( kind = rk ) patch_xmat(4,4) real ( kind = rk ) patch_ymat(4,4) real ( kind = rk ) patch_zmat(4,4) integer point real ( kind = rk ) point_uv(2,point_num) real ( kind = rk ) point_xyz(3,point_num) integer rectangle_node(16,rectangle_num) real ( kind = rk ) u real ( kind = rk ) uvec(4) real ( kind = rk ) v real ( kind = rk ) vvec(4) patch_xmat(1:4,1:4) = reshape ( node_xyz(1,rectangle_node(1:16,patch)), & (/ 4, 4 /) ) patch_ymat(1:4,1:4) = reshape ( node_xyz(2,rectangle_node(1:16,patch)), & (/ 4, 4 /) ) patch_zmat(1:4,1:4) = reshape ( node_xyz(3,rectangle_node(1:16,patch)), & (/ 4, 4 /) ) do point = 1, point_num u = point_uv(1,point) v = point_uv(2,point) uvec(1:4) = (/ ( 1.0D+00 - u )**3, & 3.0D+00 * u * ( 1.0D+00 - u )**2, & 3.0D+00 * u**2 * ( 1.0D+00 - u ), & u**3 /) vvec(1:4) = (/ ( 1.0D+00 - v )**3, & 3.0D+00 * v * ( 1.0D+00 - v )**2, & 3.0D+00 * v**2 * ( 1.0D+00 - v ), & v**3 /) point_xyz(1,point) = dot_product ( uvec(1:4), & matmul ( patch_xmat(1:4,1:4), vvec(1:4) ) ) point_xyz(2,point) = dot_product ( uvec(1:4), & matmul ( patch_ymat(1:4,1:4), vvec(1:4) ) ) point_xyz(3,point) = dot_product ( uvec(1:4), & matmul ( patch_zmat(1:4,1:4), vvec(1:4) ) ) end do return end subroutine bezier_surface_neighbors ( rectangle_num, & rectangle_node, rectangle_neighbor ) !*****************************************************************************80 ! !! bezier_surface_neighbors() determines Bezier rectangle neighbors. ! ! Discussion: ! ! A (bicubic) Bezier surface is constructed of patches. Each ! patch is an (X,Y,Z) image of the unit rectangle in (U,V) space. ! Two patches may be said to be "neighbors" if their images share ! a side. ! ! It is perfectly possible for each Bezier patch to have NO neighbors. ! It is perfectly possible for each Bezier patch to have many neighbors ! sharing a single side. ! ! However, in the most common case, and the one handled here, we ! may assume that the (X,Y,Z) patches fit together in a way similar ! to the way quilt patches form a surface, so that most rectangles ! have four neighbors, one on each side, while boundary rectangles ! might have fewer neighbors. ! ! This routine creates a data structure recording the neighbor information. ! ! The primary amount of work occurs in sorting a list of ! 4 * RECTANGLE_NUM data items representing the sides of each patch. ! ! The nodes of a single rectangle are assumed to be numbered as follows: ! ! V ! | 13-14-15-16 ! | | | | | ! | 9-10-11-12 ! | | | | | ! | 5--6--7--8 ! | | | | | ! | 1--2--3--4 ! | ! +--------------->U ! ! We assume that a neighbor patch will agree in its (X,Y,Z) values ! for all four nodes along a single side (although the nodes might ! be listed in reverse order.) ! ! We assume that there is at most one neighbor patch on each side. ! (However, even in a simple case like the Utah teapot, this is not ! true, because at the top and bottom, many patches have one side ! that degenerates to a single point, and several such patches ! meet at that point!) ! ! We choose to number the sides of each patch as follows: ! ! ! V SIDE 3 ! | ! | 13-14-15-16 ! | S | | | | ! | I 9-10-11-12 ! | D | | | | SIDE 2 ! | E 5--6--7--8 ! | | | | | ! | 4 1--2--3--4 ! | SIDE 1 ! +--------------->U ! ! And these indices for the sides correspond to the first index ! in the RECTANGLE_NEIGHBOR array. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 September 2021 ! ! Author: ! ! John Burkardt ! ! Input: ! ! integer RECTANGLE_NUM, the number of rectangles. ! ! integer RECTANGLE_NODE(16,RECTANGLE_NUM), the nodes that make up ! each rectangle. ! ! Output: ! ! integer RECTANGLE_NEIGHBOR(4,RECTANGLE_NUM), the indices of ! the rectanges that are direct neighbors of a given rectangle. ! RECTANGLE_NEIGHBOR(1,I) is the index of the rectangle which touches ! side 1, defined by nodes 1 and 4, and so on. RECTANGLE_NEIGHBOR(1,I) ! is negative if there is no neighbor on that side. ! implicit none integer rectangle_num integer i1 integer i2 integer i3 integer i4 integer irow integer rectangle integer rectangle1 integer rectangle2 integer row(4*rectangle_num,4) integer rectangle_node(16,rectangle_num) integer rectangle_neighbor(4,rectangle_num) integer side1 integer side2 ! ! Step 1. ! From the list of vertices for rectangle T, ! construct records that describe the four side segments, by listing ! the first and last nodes of each segment. ! To make matching easier later, we sort each pair of nodes. ! do rectangle = 1, rectangle_num i1 = rectangle_node(1,rectangle) i2 = rectangle_node(4,rectangle) i3 = rectangle_node(13,rectangle) i4 = rectangle_node(16,rectangle) if ( i1 < i2 ) then row(4*(rectangle-1)+1,1:4) = (/ i1, i2, 1, rectangle /) else row(4*(rectangle-1)+1,1:4) = (/ i2, i1, 1, rectangle /) end if if ( i2 < i4 ) then row(4*(rectangle-1)+2,1:4) = (/ i2, i4, 2, rectangle /) else row(4*(rectangle-1)+2,1:4) = (/ i4, i2, 2, rectangle /) end if if ( i4 < i3 ) then row(4*(rectangle-1)+3,1:4) = (/ i4, i3, 3, rectangle /) else row(4*(rectangle-1)+3,1:4) = (/ i3, i4, 3, rectangle /) end if if ( i3 < i1 ) then row(4*(rectangle-1)+4,1:4) = (/ i3, i1, 4, rectangle /) else row(4*(rectangle-1)+4,1:4) = (/ i1, i3, 4, rectangle /) end if end do ! ! Step 2. Perform an ascending dictionary sort on the neighbor relations. ! We only intend to sort on columns 1 and 2; the routine we call here ! sorts on columns 1 through 3 but that won't hurt us. ! ! What we need is to find cases where two rectangles share an edge. ! Say they share an edge defined by the nodes I and J. Then there are ! two rows of ROW that start out ( I, J, ?, ? ). By sorting ROW, ! we make sure that these two rows occur consecutively. That will ! make it easy to notice that the rectangles are neighbors. ! call i4row_sort_a ( 4*rectangle_num, 4, row ) ! ! Step 3. Neighboring rectangles show up as consecutive rows with ! identical first two entries. Whenever you spot this happening, ! make the appropriate entries in RECTANGLE_NEIGHBOR. ! rectangle_neighbor(1:4,1:rectangle_num) = -1 irow = 1 do if ( 4 * rectangle_num <= irow ) then exit end if if ( row(irow,1) /= row(irow+1,1) .or. row(irow,2) /= row(irow+1,2) ) then irow = irow + 1 cycle end if side1 = row(irow,3) rectangle1 = row(irow,4) side2 = row(irow+1,3) rectangle2 = row(irow+1,4) rectangle_neighbor(side1,rectangle1) = rectangle2 rectangle_neighbor(side2,rectangle2) = rectangle1 irow = irow + 2 end do return end subroutine bezier_surface_node_print ( node_num, node_xyz ) !*****************************************************************************80 ! !! bezier_surface_node_print() prints nodes. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 September 2021 ! ! Author: ! ! John Burkardt ! ! Input: ! ! integer NODE_NUM, the number of nodes. ! ! real ( kind = rk ) NODE_XYZ(3,NODE_NUM), the coordinates of the ! nodes. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: dim_num = 3 integer node_num real ( kind = rk ) node_xyz(dim_num,node_num) call r8mat_transpose_print ( dim_num, node_num, node_xyz, & ' Bezier Surface Nodes:' ) return end subroutine bezier_surface_node_read ( node_file_name, node_num, node_xyz ) !*****************************************************************************80 ! !! bezier_surface_node_read() reads nodes from a node file. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 September 2021 ! ! Author: ! ! John Burkardt ! ! Input: ! ! character ( len = * ) NODE_FILE_NAME, the name of the node file. ! ! integer NODE_NUM, the number of nodes. ! ! Output: ! ! real ( kind = rk ) NODE_XYZ(3,NODE_NUM), the coordinates of the ! nodes. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: dim_num = 3 integer node_num character ( len = * ) node_file_name real ( kind = rk ) node_xyz(dim_num,node_num) call r8mat_data_read ( node_file_name, dim_num, node_num, node_xyz ) return end subroutine bezier_surface_node_size ( node_file_name, node_num ) !*****************************************************************************80 ! !! bezier_surface_node_size() counts nodes in a node file. ! ! Discussion: ! ! This version of the routine simply counts the number of lines ! in the file (although it ignores comment lines beginning with "#"). ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 September 2021 ! ! Author: ! ! John Burkardt ! ! Input: ! ! character ( len = * ) NODE_FILE_NAME, the name of the node file. ! ! Output: ! ! integer NODE_NUM, the number of nodes. ! implicit none character ( len = * ) node_file_name integer node_num call file_row_count ( node_file_name, node_num ) return end subroutine bezier_surface_node_write ( node_file_name, node_num, node_xyz ) !*****************************************************************************80 ! !! bezier_surface_node_write() writes nodes to a node file. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 September 2021 ! ! Author: ! ! John Burkardt ! ! Input: ! ! character ( len = * ) NODE_FILE_NAME, the name of the node file. ! ! integer NODE_NUM, the number of nodes. ! ! real ( kind = rk ) NODE_XYZ(3,NODE_NUM), the coordinates of the ! nodes. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: dim_num = 3 integer node_num character ( len = * ) node_file_name real ( kind = rk ) node_xyz(dim_num,node_num) call r8mat_write ( node_file_name, dim_num, node_num, node_xyz ) return end subroutine bezier_surface_rectangle_print ( rectangle_num, rectangle_node ) !*****************************************************************************80 ! !! bezier_surface_rectangle_print() prints rectangles. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 September 2021 ! ! Author: ! ! John Burkardt ! ! Input: ! ! integer RECTANGLE_NUM, the number of rectangles. ! ! integer RECTANGLE_NODE(16,RECTANGLE_NUM), the nodes ! that make up each rectangle. ! implicit none integer, parameter :: dim_num = 16 integer rectangle_num integer rectangle_node(dim_num,rectangle_num) call i4mat_transpose_print ( dim_num, rectangle_num, rectangle_node, & ' Bezier Rectangles:' ) return end subroutine bezier_surface_rectangle_read ( rectangle_file_name, & rectangle_num, rectangle_node ) !*****************************************************************************80 ! !! bezier_surface_rectangle_read() reads rectangles from a rectangle file. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 September 2021 ! ! Author: ! ! John Burkardt ! ! Input: ! ! character ( len = * ) RECTANGLE_FILE_NAME, the name of the ! rectangle file. ! ! integer RECTANGLE_NUM, the number of rectangles. ! ! Output: ! ! integer RECTANGLE_NODE(16,RECTANGLE_NUM), ! the nodes that make up each rectangle. ! implicit none integer, parameter :: dim_num = 16 integer rectangle_num character ( len = * ) rectangle_file_name integer rectangle_node(dim_num,rectangle_num) call i4mat_data_read ( rectangle_file_name, dim_num, rectangle_num, & rectangle_node ) return end subroutine bezier_surface_rectangle_size ( rectangle_file_name, rectangle_num ) !*****************************************************************************80 ! !! bezier_surface_rectangle_size() counts rectangles in a rectangle file. ! ! Discussion: ! ! This version of the routine simply counts the number of lines ! in the file (although it ignores comment lines beginning with "#"). ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 September 2021 ! ! Author: ! ! John Burkardt ! ! Input: ! ! character ( len = * ) RECTANGLE_FILE_NAME, the name of the ! rectangle file. ! ! Output: ! ! integer NODE_NUM, the number of rectangles. ! implicit none character ( len = * ) rectangle_file_name integer rectangle_num call file_row_count ( rectangle_file_name, rectangle_num ) return end subroutine bezier_surface_rectangle_write ( rectangle_file_name, & rectangle_num, rectangle_node ) !*****************************************************************************80 ! !! bezier_surface_rectangle_write() writes rectangle to a rectangle file. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 September 2021 ! ! Author: ! ! John Burkardt ! ! Input: ! ! character ( len = * ) RECTANGLE_FILE_NAME, the name of the ! rectangle file. ! ! integer RECTANGLE_NUM, the number of rectangle. ! ! integer RECTANGLE_NODE(16,RECTANGLE_NUM), the nodes that ! make up each rectangle. ! implicit none integer, parameter :: dim_num = 16 integer rectangle_num character ( len = * ) rectangle_file_name integer rectangle_node(dim_num,rectangle_num) call i4mat_write ( rectangle_file_name, dim_num, rectangle_num, & rectangle_node ) return end subroutine ch_cap ( ch ) !*****************************************************************************80 ! !! ch_cap() capitalizes a single character. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 September 2021 ! ! Author: ! ! John Burkardt ! ! Input: ! ! character CH, the character to capitalize. ! ! Output: ! ! character CH, the capitalized character. ! implicit none character ch integer itemp itemp = ichar ( ch ) if ( 97 <= itemp .and. itemp <= 122 ) then ch = char ( itemp - 32 ) end if return end function ch_eqi ( c1, c2 ) !*****************************************************************************80 ! !! ch_eqi() is a case insensitive comparison of two characters for equality. ! ! Example: ! ! CH_EQI ( 'A', 'a' ) is TRUE. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 September 2021 ! ! Author: ! ! John Burkardt ! ! Input: ! ! character C1, C2, the characters to compare. ! ! Output: ! ! logical CH_EQI, the result of the comparison. ! implicit none character c1 character c1_cap character c2 character c2_cap logical ch_eqi c1_cap = c1 c2_cap = c2 call ch_cap ( c1_cap ) call ch_cap ( c2_cap ) if ( c1_cap == c2_cap ) then ch_eqi = .true. else ch_eqi = .false. end if return end subroutine ch_to_digit ( c, digit ) !*****************************************************************************80 ! !! ch_to_digit() returns the value of a base 10 digit. ! ! Example: ! ! C DIGIT ! --- ----- ! '0' 0 ! '1' 1 ! ... ... ! '9' 9 ! ' ' 0 ! 'X' -1 ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 September 2021 ! ! Author: ! ! John Burkardt ! ! Input: ! ! character C, the decimal digit, '0' through '9' or blank ! are legal. ! ! Output: ! ! integer DIGIT, the corresponding value. If C was ! 'illegal', then DIGIT is -1. ! implicit none character c integer digit if ( lge ( c, '0' ) .and. lle ( c, '9' ) ) then digit = ichar ( c ) - 48 else if ( c == ' ' ) then digit = 0 else digit = -1 end if return end subroutine file_column_count ( input_filename, column_num ) !*****************************************************************************80 ! !! file_column_count() counts the number of columns in the first line of a file. ! ! Discussion: ! ! The file is assumed to be a simple text file. ! ! Most lines of the file is presumed to consist of COLUMN_NUM words, ! separated by spaces. There may also be some blank lines, and some ! comment lines, ! which have a "#" in column 1. ! ! The routine tries to find the first non-comment non-blank line and ! counts the number of words in that line. ! ! If all lines are blanks or comments, it goes back and tries to analyze ! a comment line. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 September 2021 ! ! Author: ! ! John Burkardt ! ! Input: ! ! character ( len = * ) INPUT_FILENAME, the name of the file. ! ! Output: ! ! integer COLUMN_NUM, the number of columns in the file. ! implicit none integer column_num logical got_one character ( len = * ) input_filename integer input_unit integer ios character ( len = 255 ) line ! ! Open the file. ! call get_unit ( input_unit ) open ( unit = input_unit, file = input_filename, status = 'old', & form = 'formatted', access = 'sequential', iostat = ios ) if ( ios /= 0 ) then column_num = -1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FILE_COLUMN_COUNT - Fatal error!' write ( *, '(a)' ) ' Could not open the file:' write ( *, '(a)' ) ' ' // trim ( input_filename ) return end if ! ! Read one line, but skip blank lines and comment lines. ! got_one = .false. do read ( input_unit, '(a)', iostat = ios ) line if ( ios /= 0 ) then exit end if if ( len_trim ( line ) == 0 ) then cycle end if if ( line(1:1) == '#' ) then cycle end if got_one = .true. exit end do if ( .not. got_one ) then rewind ( input_unit ) do read ( input_unit, '(a)', iostat = ios ) line if ( ios /= 0 ) then exit end if if ( len_trim ( line ) == 0 ) then cycle end if got_one = .true. exit end do end if close ( unit = input_unit ) if ( .not. got_one ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FILE_COLUMN_COUNT - Warning!' write ( *, '(a)' ) ' The file does not seem to contain any data.' column_num = -1 return end if call s_word_count ( line, column_num ) return end subroutine file_row_count ( input_filename, row_num ) !*****************************************************************************80 ! !! file_row_count() counts the number of row records in a file. ! ! Discussion: ! ! It does not count lines that are blank, or that begin with a ! comment symbol '#'. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 September 2021 ! ! Author: ! ! John Burkardt ! ! Input: ! ! character ( len = * ) INPUT_FILENAME, the name of the input file. ! ! Output: ! ! integer ROW_NUM, the number of rows found. ! implicit none integer bad_num integer comment_num integer ierror character ( len = * ) input_filename integer input_unit integer ios character ( len = 255 ) line integer record_num integer row_num call get_unit ( input_unit ) open ( unit = input_unit, file = input_filename, status = 'old', & iostat = ios ) if ( ios /= 0 ) then row_num = -1; ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FILE_ROW_COUNT - Fatal error!' write ( *, '(a)' ) ' Could not open the input file: ' // & trim ( input_filename ) stop end if comment_num = 0 row_num = 0 record_num = 0 bad_num = 0 do read ( input_unit, '(a)', iostat = ios ) line if ( ios /= 0 ) then ierror = record_num exit end if record_num = record_num + 1 if ( line(1:1) == '#' ) then comment_num = comment_num + 1 cycle end if if ( len_trim ( line ) == 0 ) then comment_num = comment_num + 1 cycle end if row_num = row_num + 1 end do close ( unit = input_unit ) 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 MIT license. ! ! Modified: ! ! 02 September 2021 ! ! 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 subroutine i4mat_data_read ( input_filename, m, n, table ) !*****************************************************************************80 ! !! i4mat_data_read() reads data from an I4MAT file. ! ! Discussion: ! ! An I4MAT is an array of I4's. ! ! The file may contain more than N points, but this routine ! will return after reading N points. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 September 2021 ! ! Author: ! ! John Burkardt ! ! Input: ! ! character ( len = * ) INPUT_FILENAME, the name of the input file. ! ! integer M, the spatial dimension. ! ! integer N, the number of points. ! ! Output: ! ! integer TABLE(M,N), the data. ! implicit none integer m integer n integer ierror character ( len = * ) input_filename integer input_status integer input_unit integer j character ( len = 255 ) line integer table(m,n) integer x(m) ierror = 0 call get_unit ( input_unit ) open ( unit = input_unit, file = input_filename, status = 'old', & iostat = input_status ) if ( input_status /= 0 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4MAT_DATA_READ - Fatal error!' write ( *, '(a,i8)' ) ' Could not open the input file "' // & trim ( input_filename ) // '" on unit ', input_unit stop end if j = 0 do while ( j < n ) read ( input_unit, '(a)', iostat = input_status ) line if ( input_status /= 0 ) then ierror = 2 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4MAT_DATA_READ - Fatal error!' write ( *, '(a)' ) ' Error while reading lines of data.' write ( *, '(a,i8)' ) ' Number of values expected per line M = ', m write ( *, '(a,i8)' ) ' Number of data lines read, J = ', j write ( *, '(a,i8)' ) ' Number of data lines needed, N = ', n stop end if if ( line(1:1) == '#' .or. len_trim ( line ) == 0 ) then cycle end if call s_to_i4vec ( line, m, x, ierror ) if ( ierror /= 0 ) then cycle end if j = j + 1 table(1:m,j) = x(1:m) end do close ( unit = input_unit ) return end subroutine i4mat_header_read ( input_filename, m, n ) !*****************************************************************************80 ! !! i4mat_header_read() reads the header from an I4MAT. ! ! Discussion: ! ! An I4MAT is an array of I4's. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 September 2021 ! ! Author: ! ! John Burkardt ! ! Input: ! ! character ( len = * ) INPUT_FILENAME, the name of the input file. ! ! Output: ! ! integer M, spatial dimension. ! ! integer N, the number of points. ! implicit none character ( len = * ) input_filename integer m integer n call file_column_count ( input_filename, m ) if ( m <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4MAT_HEADER_READ - Fatal error!' write ( *, '(a)' ) ' There was some kind of I/O problem while trying' write ( *, '(a)' ) ' to count the number of data columns in' write ( *, '(a)' ) ' the file "' // trim ( input_filename ) // '".' stop end if call file_row_count ( input_filename, n ) if ( n <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4MAT_HEADER_READ - Fatal error!' write ( *, '(a)' ) ' There was some kind of I/O problem while trying' write ( *, '(a)' ) ' to count the number of data rows in' write ( *, '(a)' ) ' the file "' // trim ( input_filename ) // '".' stop end if return end subroutine i4mat_print ( m, n, a, title ) !*****************************************************************************80 ! !! i4mat_print() prints an I4MAT. ! ! Discussion: ! ! An I4MAT is a rectangular array of integer values. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 September 2021 ! ! Author: ! ! John Burkardt ! ! Input: ! ! integer M, the number of rows in A. ! ! integer N, the number of columns in A. ! ! integer A(M,N), the matrix to be printed. ! ! character ( len = * ) TITLE, a title to be printed first. ! TITLE may be blank. ! implicit none integer m integer n integer a(m,n) integer ihi integer ilo integer jhi integer jlo character ( len = * ) title ilo = 1 ihi = m jlo = 1 jhi = n call i4mat_print_some ( m, n, a, ilo, jlo, ihi, jhi, title ) return end subroutine i4mat_print_some ( m, n, a, ilo, jlo, ihi, jhi, title ) !*****************************************************************************80 ! !! i4mat_print_some() prints some of an I4MAT. ! ! Discussion: ! ! An I4MAT is a rectangular array of integer values. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 September 2021 ! ! Author: ! ! John Burkardt ! ! Input: ! ! integer M, N, the number of rows and columns. ! ! integer A(M,N), an M by N matrix to be printed. ! ! integer ILO, JLO, the first row and column to print. ! ! integer IHI, JHI, the last row and column to print. ! ! character ( len = * ) TITLE, an optional title. ! implicit none integer, parameter :: incx = 10 integer m integer n integer a(m,n) character ( len = 8 ) ctemp(incx) integer i integer i2hi integer i2lo integer ihi integer ilo integer inc integer j integer j2 integer j2hi integer j2lo integer jhi integer jlo character ( len = * ) title write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) do j2lo = max ( jlo, 1 ), min ( jhi, n ), incx j2hi = j2lo + incx - 1 j2hi = min ( j2hi, n ) j2hi = min ( j2hi, jhi ) inc = j2hi + 1 - j2lo write ( *, '(a)' ) ' ' do j = j2lo, j2hi j2 = j + 1 - j2lo write ( ctemp(j2), '(i8)' ) j end do write ( *, '('' Col '',10a8)' ) ctemp(1:inc) write ( *, '(a)' ) ' Row' write ( *, '(a)' ) ' ' i2lo = max ( ilo, 1 ) i2hi = min ( ihi, m ) do i = i2lo, i2hi do j2 = 1, inc j = j2lo - 1 + j2 write ( ctemp(j2), '(i8)' ) a(i,j) end do write ( *, '(i5,1x,10a8)' ) i, ( ctemp(j), j = 1, inc ) end do end do return end subroutine i4mat_transpose_print ( m, n, a, title ) !*****************************************************************************80 ! !! i4mat_transpose_print() prints an I4MAT, transposed. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 September 2021 ! ! Author: ! ! John Burkardt ! ! Input: ! ! integer M, N, the number of rows and columns. ! ! integer A(M,N), an M by N matrix to be printed. ! ! character ( len = * ) TITLE, an optional title. ! implicit none integer m integer n integer a(m,n) character ( len = * ) title call i4mat_transpose_print_some ( m, n, a, 1, 1, m, n, title ) return end subroutine i4mat_transpose_print_some ( m, n, a, ilo, jlo, ihi, jhi, title ) !*****************************************************************************80 ! !! i4mat_transpose_print_some() prints some of the transpose of an I4MAT. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 September 2021 ! ! Author: ! ! John Burkardt ! ! Input: ! ! integer M, N, the number of rows and columns. ! ! integer A(M,N), an M by N matrix to be printed. ! ! integer ILO, JLO, the first row and column to print. ! ! integer IHI, JHI, the last row and column to print. ! ! character ( len = * ) TITLE, an optional title. ! implicit none integer, parameter :: incx = 10 integer m integer n integer a(m,n) character ( len = 7 ) ctemp(incx) integer i integer i2 integer i2hi integer i2lo integer ihi integer ilo integer inc integer j integer j2hi integer j2lo integer jhi integer jlo character ( len = * ) title write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) do i2lo = max ( ilo, 1 ), min ( ihi, m ), incx i2hi = i2lo + incx - 1 i2hi = min ( i2hi, m ) i2hi = min ( i2hi, ihi ) inc = i2hi + 1 - i2lo write ( *, '(a)' ) ' ' do i = i2lo, i2hi i2 = i + 1 - i2lo write ( ctemp(i2), '(i7)') i end do write ( *, '('' Row '',10a7)' ) ctemp(1:inc) write ( *, '(a)' ) ' Col' write ( *, '(a)' ) ' ' j2lo = max ( jlo, 1 ) j2hi = min ( jhi, n ) do j = j2lo, j2hi do i2 = 1, inc i = i2lo - 1 + i2 write ( ctemp(i2), '(i7)' ) a(i,j) end do write ( *, '(i5,1x,10a7)' ) j, ( ctemp(i), i = 1, inc ) end do end do return end subroutine i4mat_write ( output_filename, m, n, table ) !*****************************************************************************80 ! !! i4mat_write() writes an I4MAT file. ! ! Discussion: ! ! An I4MAT is an array of I4's. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 September 2021 ! ! Author: ! ! John Burkardt ! ! Input: ! ! character ( len = * ) OUTPUT_FILENAME, the output file name. ! ! integer M, the spatial dimension. ! ! integer N, the number of points. ! ! integer TABLE(M,N), the data. ! implicit none integer m integer n integer j character ( len = * ) output_filename integer output_status integer output_unit character ( len = 30 ) string integer 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)' ) 'I4MAT_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. ! if ( 0 < m .and. 0 < n ) then write ( string, '(a1,i8,a4)' ) '(', m, 'i10)' ! ! 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 i4row_compare ( m, n, a, i, j, isgn ) !*****************************************************************************80 ! !! i4row_compare() compares two rows of an I4ROW. ! ! Discussion: ! ! An I4ROW is an M by N array of integer values, regarded ! as an array of M rows of length N. ! ! Example: ! ! Input: ! ! M = 3, N = 4, I = 2, J = 3 ! ! A = ( ! 1 2 3 4 ! 5 6 7 8 ! 9 10 11 12 ) ! ! Output: ! ! ISGN = -1 ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 September 2021 ! ! Author: ! ! John Burkardt ! ! Input: ! ! integer M, N, the number of rows and columns. ! ! integer A(M,N), an array of M rows of vectors of length N. ! ! integer I, J, the rows to be compared. ! I and J must be between 1 and M. ! ! Output: ! ! integer ISGN, the results of the comparison: ! -1, row I < row J, ! 0, row I = row J, ! +1, row J < row I. ! implicit none integer m integer n integer a(m,n) integer i integer isgn integer j integer k ! ! Check that I and J are legal. ! if ( i < 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4ROW_COMPARE - Fatal error!' write ( *, '(a)' ) ' Row index I is less than 1.' write ( *, '(a,i8)' ) ' I = ', i stop else if ( m < i ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4ROW_COMPARE - Fatal error!' write ( *, '(a)' ) ' Row index I is out of bounds.' write ( *, '(a,i8)' ) ' I = ', i write ( *, '(a,i8)' ) ' Maximum legal value is M = ', m stop end if if ( j < 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4ROW_COMPARE - Fatal error!' write ( *, '(a)' ) ' Row index J is less than 1.' write ( *, '(a,i8)' ) ' J = ', j stop else if ( m < j ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4ROW_COMPARE - Fatal error!' write ( *, '(a)' ) ' Row index J is out of bounds.' write ( *, '(a,i8)' ) ' J = ', j write ( *, '(a,i8)' ) ' Maximum legal value is M = ', m stop end if isgn = 0 if ( i == j ) then return end if k = 1 do while ( k <= n ) if ( a(i,k) < a(j,k) ) then isgn = -1 return else if ( a(j,k) < a(i,k) ) then isgn = +1 return end if k = k + 1 end do return end subroutine i4row_sort_a ( m, n, a ) !*****************************************************************************80 ! !! i4row_sort_a() ascending sorts the rows of an I4ROW. ! ! Discussion: ! ! An I4ROW is an M by N array of integer values, regarded ! as an array of M rows of length N. ! ! In lexicographic order, the statement "X < Y", applied to two ! vectors X and Y of length M, means that there is some index I, with ! 1 <= I <= M, with the property that ! ! X(J) = Y(J) for J < I, ! and ! X(I) < Y(I). ! ! In other words, X is less than Y if, at the first index where they ! differ, the X value is less than the Y value. ! ! Example: ! ! Input: ! ! M = 5, N = 3 ! ! A = ! 3 2 1 ! 2 4 3 ! 3 1 8 ! 2 4 2 ! 1 9 9 ! ! Output: ! ! A = ! 1 9 9 ! 2 4 2 ! 2 4 3 ! 3 1 8 ! 3 2 1 ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 September 2021 ! ! Author: ! ! John Burkardt ! ! Input: ! ! integer M, the number of rows of A. ! ! integer N, the number of columns of A. ! ! integer A(M,N), the array of M rows of N-vectors. ! ! Output: ! ! integer A(M,N), the rows of A have been sorted in ascending ! lexicographic order. ! implicit none integer m integer n integer a(m,n) integer i integer indx integer isgn integer j if ( m <= 1 ) then return end if if ( n <= 0 ) then return end if ! ! Initialize. ! i = 0 indx = 0 isgn = 0 j = 0 ! ! Call the external heap sorter. ! do call sort_heap_external ( m, indx, i, j, isgn ) ! ! Interchange the I and J objects. ! if ( 0 < indx ) then call i4row_swap ( m, n, a, i, j ) ! ! Compare the I and J objects. ! else if ( indx < 0 ) then call i4row_compare ( m, n, a, i, j, isgn ) else if ( indx == 0 ) then exit end if end do return end subroutine i4row_swap ( m, n, a, i1, i2 ) !*****************************************************************************80 ! !! i4row_swap() swaps two rows of an I4ROW. ! ! Discussion: ! ! An I4ROW is an M by N array of integer values, regarded ! as an array of M rows of length N. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 September 2021 ! ! Author: ! ! John Burkardt ! ! Input: ! ! integer M, N, the number of rows and columns. ! ! integer A(M,N), an array of data. ! ! integer I1, I2, the two rows to swap. ! ! Output: ! ! integer A(M,N), an array after row swapping. ! implicit none integer m integer n integer a(m,n) integer i1 integer i2 integer row(n) ! ! Check. ! if ( i1 < 1 .or. m < i1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4ROW_SWAP - Fatal error!' write ( *, '(a)' ) ' I1 is out of range.' stop end if if ( i2 < 1 .or. m < i2 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4ROW_SWAP - Fatal error!' write ( *, '(a)' ) ' I2 is out of range.' stop end if if ( i1 == i2 ) then return end if row(1:n) = a(i1,1:n) a(i1,1:n) = a(i2,1:n) a(i2,1:n) = row(1:n) return end subroutine r8mat_data_read ( input_filename, m, n, table ) !*****************************************************************************80 ! !! r8mat_data_read() reads data from an R8MAT file. ! ! Discussion: ! ! An R8MAT is an array of R8 values. ! ! Discussion: ! ! The file may contain more than N points, but this routine will ! return after reading N of them. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 September 2021 ! ! Author: ! ! John Burkardt ! ! Input: ! ! character ( len = * ) INPUT_FILENAME, the name of the input file. ! ! integer M, the spatial dimension. ! ! integer N, the number of points. ! ! Output: ! ! real ( kind = rk ) TABLE(M,N), the data. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n integer ierror character ( len = * ) input_filename integer input_status integer input_unit integer j character ( len = 255 ) line real ( kind = rk ) table(m,n) real ( kind = rk ) x(m) ierror = 0 call get_unit ( input_unit ) open ( unit = input_unit, file = input_filename, status = 'old', & iostat = input_status ) if ( input_status /= 0 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8MAT_DATA_READ - Fatal error!' write ( *, '(a,i8)' ) ' Could not open the input file "' // & trim ( input_filename ) // '" on unit ', input_unit stop end if j = 0 do while ( j < n ) read ( input_unit, '(a)', iostat = input_status ) line if ( input_status /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8MAT_DATA_READ - Fatal error!' write ( *, '(a)' ) ' Error while reading lines of data.' write ( *, '(a,i8)' ) ' Number of values expected per line M = ', m write ( *, '(a,i8)' ) ' Number of data lines read, J = ', j write ( *, '(a,i8)' ) ' Number of data lines needed, N = ', n stop end if if ( line(1:1) == '#' .or. len_trim ( line ) == 0 ) then cycle end if call s_to_r8vec ( line, m, x, ierror ) if ( ierror /= 0 ) then cycle end if j = j + 1 table(1:m,j) = x(1:m) end do close ( unit = input_unit ) return end subroutine r8mat_header_read ( input_filename, m, n ) !*****************************************************************************80 ! !! r8mat_header_read() reads the header from an R8MAT file. ! ! Discussion: ! ! An R8MAT is an array of R8 values. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 September 2021 ! ! Author: ! ! John Burkardt ! ! Input: ! ! character ( len = * ) INPUT_FILENAME, the name of the input file. ! ! Output: ! ! integer M, spatial dimension. ! ! integer N, the number of points. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) character ( len = * ) input_filename integer m integer n call file_column_count ( input_filename, m ) if ( m <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8MAT_HEADER_READ - Fatal error!' write ( *, '(a)' ) ' There was some kind of I/O problem while trying' write ( *, '(a)' ) ' to count the number of data columns in' write ( *, '(a)' ) ' the file "' // trim ( input_filename ) // '".' stop end if call file_row_count ( input_filename, n ) if ( n <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8MAT_HEADER_READ - Fatal error!' write ( *, '(a)' ) ' There was some kind of I/O problem while trying' write ( *, '(a)' ) ' to count the number of data rows in' write ( *, '(a)' ) ' the file "' // trim ( input_filename ) // '".' stop end if return end subroutine r8mat_print ( m, n, a, title ) !*****************************************************************************80 ! !! r8mat_print() prints an R8MAT. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 September 2021 ! ! Author: ! ! John Burkardt ! ! Input: ! ! integer M, the number of rows in A. ! ! integer N, the number of columns in A. ! ! real ( kind = rk ) A(M,N), the matrix. ! ! character ( len = * ) TITLE, a title to be printed. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a(m,n) character ( len = * ) title call r8mat_print_some ( m, n, a, 1, 1, m, n, title ) return end subroutine r8mat_print_some ( m, n, a, ilo, jlo, ihi, jhi, title ) !*****************************************************************************80 ! !! r8mat_print_some() prints some of an R8MAT. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 September 2021 ! ! Author: ! ! John Burkardt ! ! Input: ! ! integer M, N, the number of rows and columns. ! ! real ( kind = rk ) A(M,N), an M by N matrix to be printed. ! ! integer ILO, JLO, the first row and column to print. ! ! integer IHI, JHI, the last row and column to print. ! ! character ( len = * ) TITLE, an optional title. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: incx = 5 integer m integer n real ( kind = rk ) a(m,n) character ( len = 14 ) ctemp(incx) integer i integer i2hi integer i2lo integer ihi integer ilo integer inc integer j integer j2 integer j2hi integer j2lo integer jhi integer jlo character ( len = * ) title if ( 0 < len_trim ( title ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) end if do j2lo = max ( jlo, 1 ), min ( jhi, n ), incx j2hi = j2lo + incx - 1 j2hi = min ( j2hi, n ) j2hi = min ( j2hi, jhi ) inc = j2hi + 1 - j2lo write ( *, '(a)' ) ' ' do j = j2lo, j2hi j2 = j + 1 - j2lo write ( ctemp(j2), '(i7,7x)') j end do write ( *, '('' Col '',5a14)' ) ctemp(1:inc) write ( *, '(a)' ) ' Row' write ( *, '(a)' ) ' ' i2lo = max ( ilo, 1 ) i2hi = min ( ihi, m ) do i = i2lo, i2hi do j2 = 1, inc j = j2lo - 1 + j2 if ( a(i,j) == real ( int ( a(i,j) ), kind = rk ) ) then write ( ctemp(j2), '(f8.0,6x)' ) a(i,j) else write ( ctemp(j2), '(g14.6)' ) a(i,j) end if end do write ( *, '(i5,1x,5a14)' ) i, ( ctemp(j), j = 1, inc ) end do end do return end subroutine r8mat_transpose_print ( m, n, a, title ) !*****************************************************************************80 ! !! r8mat_transpose_print() prints an R8MAT, transposed. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 September 2021 ! ! Author: ! ! John Burkardt ! ! Input: ! ! integer M, N, the number of rows and columns. ! ! real ( kind = rk ) A(M,N), an M by N matrix to be printed. ! ! character ( len = * ) TITLE, an optional title. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a(m,n) character ( len = * ) title call r8mat_transpose_print_some ( m, n, a, 1, 1, m, n, title ) return end subroutine r8mat_transpose_print_some ( m, n, a, ilo, jlo, ihi, jhi, title ) !*****************************************************************************80 ! !! r8mat_transpose_print_some() prints some of an R8MAT, transposed. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 September 2021 ! ! Author: ! ! John Burkardt ! ! Input: ! ! integer M, N, the number of rows and columns. ! ! real ( kind = rk ) A(M,N), an M by N matrix to be printed. ! ! integer ILO, JLO, the first row and column to print. ! ! integer IHI, JHI, the last row and column to print. ! ! character ( len = * ) TITLE, an optional title. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: incx = 5 integer m integer n real ( kind = rk ) a(m,n) character ( len = 14 ) ctemp(incx) integer i integer i2 integer i2hi integer i2lo integer ihi integer ilo integer inc integer j integer j2hi integer j2lo integer jhi integer jlo character ( len = * ) title write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) do i2lo = max ( ilo, 1 ), min ( ihi, m ), incx i2hi = i2lo + incx - 1 i2hi = min ( i2hi, m ) i2hi = min ( i2hi, ihi ) inc = i2hi + 1 - i2lo write ( *, '(a)' ) ' ' do i = i2lo, i2hi i2 = i + 1 - i2lo write ( ctemp(i2), '(i7,7x)') i end do write ( *, '('' Row '',5a14)' ) ctemp(1:inc) write ( *, '(a)' ) ' Col' j2lo = max ( jlo, 1 ) j2hi = min ( jhi, n ) do j = j2lo, j2hi do i2 = 1, inc i = i2lo - 1 + i2 write ( ctemp(i2), '(g14.6)' ) a(i,j) end do write ( *, '(i5,1x,5a14)' ) j, ( ctemp(i), i = 1, inc ) end do 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 MIT license. ! ! Modified: ! ! 02 September 2021 ! ! Author: ! ! John Burkardt ! ! Input: ! ! character ( len = * ) OUTPUT_FILENAME, the output file name. ! ! integer M, the spatial dimension. ! ! integer N, the number of points. ! ! real ( kind = rk ) TABLE(M,N), the data. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n integer j character ( len = * ) output_filename integer output_status integer output_unit character ( len = 30 ) string real ( kind = rk ) 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_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 MIT license. ! ! Modified: ! ! 02 September 2021 ! ! Author: ! ! John Burkardt ! ! Input: ! ! character ( len = * ) S, the string to be transformed. ! ! Output: ! ! character ( len = * ) S, the transformed string. ! implicit none character c integer get integer put integer nchar character ( len = * ) s character, parameter :: TAB = char ( 9 ) put = 0 nchar = len_trim ( s ) do get = 1, nchar c = s(get:get) if ( c /= ' ' .and. c /= TAB ) then put = put + 1 s(put:put) = c end if end do s(put+1:nchar) = ' ' return end subroutine s_to_i4 ( s, ival, ierror, length ) !*****************************************************************************80 ! !! s_to_i4() reads an I4 from a string. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 September 2021 ! ! Author: ! ! John Burkardt ! ! Input: ! ! character ( len = * ) S, a string to be examined. ! ! Output: ! ! integer IVAL, the integer value read from the string. ! If the string is blank, then IVAL will be returned 0. ! ! integer IERROR, an error flag. ! 0, no error. ! 1, an error occurred. ! ! integer LENGTH, the number of characters of S used to make IVAL. ! implicit none character c integer i integer ierror integer isgn integer istate integer ival integer length character ( len = * ) s ierror = 0 istate = 0 isgn = 1 ival = 0 do i = 1, len_trim ( s ) c = s(i:i) ! ! Haven't read anything. ! if ( istate == 0 ) then if ( c == ' ' ) then else if ( c == '-' ) then istate = 1 isgn = -1 else if ( c == '+' ) then istate = 1 isgn = + 1 else if ( lle ( '0', c ) .and. lle ( c, '9' ) ) then istate = 2 ival = ichar ( c ) - ichar ( '0' ) else ierror = 1 return end if ! ! Have read the sign, expecting digits. ! else if ( istate == 1 ) then if ( c == ' ' ) then else if ( lle ( '0', c ) .and. lle ( c, '9' ) ) then istate = 2 ival = ichar ( c ) - ichar ( '0' ) else ierror = 1 return end if ! ! Have read at least one digit, expecting more. ! else if ( istate == 2 ) then if ( lle ( '0', c ) .and. lle ( c, '9' ) ) then ival = 10 * ival + ichar ( c ) - ichar ( '0' ) else ival = isgn * ival 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 ( istate == 2 ) then ival = isgn * ival length = len_trim ( s ) else ierror = 1 length = 0 end if return end subroutine s_to_i4vec ( s, n, ivec, ierror ) !*****************************************************************************80 ! !! s_to_i4vec() reads an I4VEC from a string. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 September 2021 ! ! Author: ! ! John Burkardt ! ! Input: ! ! character ( len = * ) S, the string to be read. ! ! integer N, the number of values expected. ! ! Output: ! ! integer IVEC(N), the values read from the string. ! ! integer IERROR, error flag. ! 0, no errors occurred. ! -K, could not read data for entries -K through N. ! implicit none integer n integer i integer ierror integer ilo integer ivec(n) integer length character ( len = * ) s i = 0 ierror = 0 ilo = 1 do while ( i < n ) i = i + 1 call s_to_i4 ( s(ilo:), ivec(i), ierror, length ) if ( ierror /= 0 ) then ierror = -i exit end if ilo = ilo + length end do return end subroutine s_to_r8 ( s, dval, ierror, length ) !*****************************************************************************80 ! !! s_to_r8() reads an R8 from a string. ! ! Discussion: ! ! The routine will read as many characters as possible until it reaches ! the end of the string, or encounters a character which cannot be ! part of the number. ! ! Legal input is: ! ! 1 blanks, ! 2 '+' or '-' sign, ! 2.5 blanks ! 3 integer part, ! 4 decimal point, ! 5 fraction part, ! 6 'E' or 'e' or 'D' or 'd', exponent marker, ! 7 exponent sign, ! 8 exponent integer part, ! 9 exponent decimal point, ! 10 exponent fraction part, ! 11 blanks, ! 12 final comma or semicolon, ! ! with most quantities optional. ! ! Example: ! ! S DVAL ! ! '1' 1.0 ! ' 1 ' 1.0 ! '1A' 1.0 ! '12,34,56' 12.0 ! ' 34 7' 34.0 ! '-1E2ABCD' -100.0 ! '-1X2ABCD' -1.0 ! ' 2E-1' 0.2 ! '23.45' 23.45 ! '-4.2E+2' -420.0 ! '17d2' 1700.0 ! '-14e-2' -0.14 ! 'e2' 100.0 ! '-12.73e-9.23' -12.73 * 10.0^(-9.23) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 September 2021 ! ! Author: ! ! John Burkardt ! ! Input: ! ! character ( len = * ) S, the string containing the ! data to be read. Reading will begin at position 1 and ! terminate at the end of the string, or when no more ! characters can be read to form a legal real. Blanks, ! commas, or other nonnumeric data will, in particular, ! cause the conversion to halt. ! ! Output: ! ! real ( kind = rk ) DVAL, the value read from the string. ! ! integer IERROR, error flag. ! 0, no errors occurred. ! 1, 2, 6 or 7, the input number was garbled. The ! value of IERROR is the last type of input successfully ! read. For instance, 1 means initial blanks, 2 means ! a plus or minus sign, and so on. ! ! integer LENGTH, the number of characters read ! to form the number, including any terminating ! characters such as a trailing comma or blanks. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) logical ch_eqi character c real ( kind = rk ) dval integer ierror integer ihave integer isgn integer iterm integer jbot integer jsgn integer jtop integer length integer nchar integer ndig real ( kind = rk ) rbot real ( kind = rk ) rexp real ( kind = rk ) rtop character ( len = * ) s nchar = len_trim ( s ) ierror = 0 dval = 0.0D+00 length = -1 isgn = 1 rtop = 0 rbot = 1 jsgn = 1 jtop = 0 jbot = 1 ihave = 1 iterm = 0 do length = length + 1 if ( nchar < length+1 ) then exit end if c = s(length+1:length+1) ! ! Blank character. ! if ( c == ' ' ) then if ( ihave == 2 ) then else if ( ihave == 6 .or. ihave == 7 ) then iterm = 1 else if ( 1 < ihave ) then ihave = 11 end if ! ! Comma. ! else if ( c == ',' .or. c == ';' ) then if ( ihave /= 1 ) then iterm = 1 ihave = 12 length = length + 1 end if ! ! Minus sign. ! else if ( c == '-' ) then if ( ihave == 1 ) then ihave = 2 isgn = -1 else if ( ihave == 6 ) then ihave = 7 jsgn = -1 else iterm = 1 end if ! ! Plus sign. ! else if ( c == '+' ) then if ( ihave == 1 ) then ihave = 2 else if ( ihave == 6 ) then ihave = 7 else iterm = 1 end if ! ! Decimal point. ! else if ( c == '.' ) then if ( ihave < 4 ) then ihave = 4 else if ( 6 <= ihave .and. ihave <= 8 ) then ihave = 9 else iterm = 1 end if ! ! Scientific notation exponent marker. ! else if ( ch_eqi ( c, 'E' ) .or. ch_eqi ( c, 'D' ) ) then if ( ihave < 6 ) then ihave = 6 else iterm = 1 end if ! ! Digit. ! else if ( ihave < 11 .and. lle ( '0', c ) .and. lle ( c, '9' ) ) then if ( ihave <= 2 ) then ihave = 3 else if ( ihave == 4 ) then ihave = 5 else if ( ihave == 6 .or. ihave == 7 ) then ihave = 8 else if ( ihave == 9 ) then ihave = 10 end if call ch_to_digit ( c, ndig ) if ( ihave == 3 ) then rtop = 10.0D+00 * rtop + real ( ndig, kind = rk ) else if ( ihave == 5 ) then rtop = 10.0D+00 * rtop + real ( ndig, kind = rk ) rbot = 10.0D+00 * rbot else if ( ihave == 8 ) then jtop = 10 * jtop + ndig else if ( ihave == 10 ) then jtop = 10 * jtop + ndig jbot = 10 * jbot end if ! ! Anything else is regarded as a terminator. ! else iterm = 1 end if ! ! If we haven't seen a terminator, and we haven't examined the ! entire string, go get the next character. ! if ( iterm == 1 ) then exit end if end do ! ! If we haven't seen a terminator, and we have examined the ! entire string, then we're done, and LENGTH is equal to NCHAR. ! if ( iterm /= 1 .and. length+1 == nchar ) then length = nchar end if ! ! Number seems to have terminated. Have we got a legal number? ! Not if we terminated in states 1, 2, 6 or 7! ! if ( ihave == 1 .or. ihave == 2 .or. ihave == 6 .or. ihave == 7 ) then ierror = ihave write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'S_TO_r8 - Serious error!' write ( *, '(a)' ) ' Illegal or nonnumeric input:' write ( *, '(a)' ) ' ' // trim ( s ) return end if ! ! Number seems OK. Form it. ! if ( jtop == 0 ) then rexp = 1.0D+00 else if ( jbot == 1 ) then rexp = 10.0D+00 ** ( jsgn * jtop ) else rexp = 10.0D+00 ** ( real ( jsgn * jtop, kind = rk ) & / real ( jbot, kind = rk ) ) end if end if dval = real ( isgn, kind = rk ) * rexp * rtop / rbot return end subroutine s_to_r8vec ( s, n, rvec, ierror ) !*****************************************************************************80 ! !! s_to_r8vec() reads an R8VEC from a string. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 September 2021 ! ! Author: ! ! John Burkardt ! ! Input: ! ! character ( len = * ) S, the string to be read. ! ! integer N, the number of values expected. ! ! Output: ! ! real ( kind = rk ) RVEC(N), the values read from the string. ! ! integer IERROR, error flag. ! 0, no errors occurred. ! -K, could not read data for entries -K through N. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n integer i integer ierror integer ilo integer lchar real ( kind = rk ) rvec(n) character ( len = * ) s i = 0 ierror = 0 ilo = 1 do while ( i < n ) i = i + 1 call s_to_r8 ( s(ilo:), rvec(i), ierror, lchar ) if ( ierror /= 0 ) then ierror = -i exit end if ilo = ilo + lchar end do return end subroutine s_word_count ( s, nword ) !*****************************************************************************80 ! !! s_word_count() counts the number of "words" in a string. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 September 2021 ! ! Author: ! ! John Burkardt ! ! Input: ! ! character ( len = * ) S, the string to be examined. ! ! Output: ! ! integer NWORD, the number of "words" in the string. ! Words are presumed to be separated by one or more blanks. ! implicit none logical blank integer i integer lens integer nword character ( len = * ) s nword = 0 lens = len ( s ) if ( lens <= 0 ) then return end if blank = .true. do i = 1, lens if ( s(i:i) == ' ' ) then blank = .true. else if ( blank ) then nword = nword + 1 blank = .false. end if end do return end subroutine sort_heap_external ( n, indx, i, j, isgn ) !*****************************************************************************80 ! !! SORT_HEAP_EXTERNAL externally sorts a list of items into ascending order. ! ! Discussion: ! ! The actual list of data is not passed to the routine. Hence this ! routine may be used to sort integers, reals, numbers, names, ! dates, shoe sizes, and so on. After each call, the routine asks ! the user to compare or interchange two items, until a special ! return value signals that the sorting is completed. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 September 2021 ! ! Author: ! ! FORTRAN77 original version by Albert Nijenhuis, Herbert Wilf. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Albert Nijenhuis, Herbert Wilf, ! Combinatorial Algorithms, ! Academic Press, 1978, second edition, ! ISBN 0-12-519260-6. ! ! Input: ! ! integer N, the number of items to be sorted. ! ! integer INDX, the main communication signal. ! The user must set INDX to 0 before the first call. ! ! integer ISGN, results of comparison of elements I and J. ! (Used only when the previous call returned INDX less than 0). ! ISGN <= 0 means I is less than or equal to J; ! 0 <= ISGN means I is greater than or equal to J. ! ! Output: ! ! integer INDX, the main communication signal. ! On return, if INDX is ! # greater than 0, ! * interchange items I and J; ! * call again. ! # less than 0, ! * compare items I and J; ! * set ISGN = -1 if I < J, ISGN = +1 if J < I; ! * call again. ! # equal to 0, the sorting is done. ! ! integer I, J, the indices of two items. ! On return with INDX positive, elements I and J should be interchanged. ! On return with INDX negative, elements I and J should be compared, and ! the result reported in ISGN on the next call. ! implicit none integer i integer, save :: i_save = 0 integer indx integer isgn integer j integer, save :: j_save = 0 integer, save :: k = 0 integer, save :: k1 = 0 integer n integer, save :: n1 = 0 ! ! INDX = 0: This is the first call. ! if ( indx == 0 ) then i_save = 0 j_save = 0 k = n / 2 k1 = k n1 = n ! ! INDX < 0: The user is returning the results of a comparison. ! else if ( indx < 0 ) then if ( indx == -2 ) then if ( isgn < 0 ) then i_save = i_save + 1 end if j_save = k1 k1 = i_save indx = -1 i = i_save j = j_save return end if if ( 0 < isgn ) then indx = 2 i = i_save j = j_save return end if if ( k <= 1 ) then if ( n1 == 1 ) then i_save = 0 j_save = 0 indx = 0 else i_save = n1 n1 = n1 - 1 j_save = 1 indx = 1 end if i = i_save j = j_save return end if k = k - 1 k1 = k ! ! 0 < INDX, the user was asked to make an interchange. ! else if ( indx == 1 ) then k1 = k end if do i_save = 2 * k1 if ( i_save == n1 ) then j_save = k1 k1 = i_save indx = -1 i = i_save j = j_save return else if ( i_save <= n1 ) then j_save = i_save + 1 indx = -2 i = i_save j = j_save return end if if ( k <= 1 ) then exit end if k = k - 1 k1 = k end do if ( n1 == 1 ) then i_save = 0 j_save = 0 indx = 0 i = i_save j = j_save else i_save = n1 n1 = n1 - 1 j_save = 1 indx = 1 i = i_save j = j_save end if return end subroutine timestamp ( ) !*****************************************************************************80 ! !! timestamp() prints the current YMDHMS date as a time stamp. ! ! Example: ! ! May 31 2001 9:45:54.872 AM ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 September 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 ( *, '(a,1x,i2,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & trim ( month(m) ), d, y, h, ':', n, ':', s, '.', mm, trim ( ampm ) return end