function i4_log_10 ( i ) !*****************************************************************************80 ! !! I4_LOG_10 returns the integer part of the logarithm base 10 of an I4. ! ! Discussion: ! ! I4_LOG_10 ( I ) + 1 is the number of decimal digits in I. ! ! An I4 is an integer value. ! ! Example: ! ! I I4_LOG_10 ! ----- -------- ! 0 0 ! 1 0 ! 2 0 ! 9 0 ! 10 1 ! 11 1 ! 99 1 ! 100 2 ! 101 2 ! 999 2 ! 1000 3 ! 1001 3 ! 9999 3 ! 10000 4 ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 08 June 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer I, the number whose logarithm base 10 ! is desired. ! ! Output, integer I4_LOG_10, the integer part of the ! logarithm base 10 of the absolute value of X. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer i integer i_abs integer i4_log_10 integer ten_pow if ( i == 0 ) then i4_log_10 = 0 else i4_log_10 = 0 ten_pow = 10 i_abs = abs ( i ) do while ( ten_pow <= i_abs ) i4_log_10 = i4_log_10 + 1 ten_pow = ten_pow * 10 end do end if return end function i4_uniform_ab ( a, b ) !*****************************************************************************80 ! !! i4_uniform_ab returns a scaled pseudorandom I4 between A and B. ! ! Discussion: ! ! An I4 is an integer value. ! ! The pseudorandom number will be scaled to be uniformly distributed ! between A and B. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 October 2012 ! ! 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. ! ! Input: ! ! integer A, B, the limits of the interval. ! ! Output: ! ! integer I4_UNIFORM_AB, a number between A and B. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer a integer b integer i4_uniform_ab real r integer value call random_number ( harvest = r ) ! ! Scale R to lie between A-0.5 and B+0.5. ! r = ( 1.0E+00 - r ) * ( real ( min ( a, b ) ) - 0.5E+00 ) & + r * ( real ( max ( a, b ) ) + 0.5E+00 ) ! ! Use rounding to convert R to an integer between A and B. ! value = nint ( r ) value = max ( value, min ( a, b ) ) value = min ( value, max ( a, b ) ) i4_uniform_ab = value return end subroutine i4vec_print ( n, a, title ) !*****************************************************************************80 ! !! I4VEC_PRINT prints an I4VEC. ! ! Discussion: ! ! An I4VEC is a vector of I4's. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 May 2010 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of components of the vector. ! ! Input, integer A(N), the vector to be printed. ! ! Input, character ( len = * ) TITLE, a title. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n integer a(n) integer i character ( len = * ) title if ( 0 < len_trim ( title ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) end if write ( *, '(a)' ) ' ' do i = 1, n write ( *, '(2x,i8,a,2x,i12)' ) i, ':', a(i) end do return end function perm1_check ( n, p ) !*****************************************************************************80 ! !! PERM1_CHECK checks a 1-based permutation. ! ! Discussion: ! ! The routine verifies that each of the integers from 1 to ! to N occurs among the N entries of the permutation. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 24 May 2015 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of entries. ! ! Input, integer P(N), the array to check. ! ! Output, integer PERM1_CHECK: ! 0, P is a 1-based permutation. ! 1, P is not a 1-based permutation. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n integer ierror integer location integer p(n) integer perm1_check integer value ierror = 0 do value = 1, n ierror = 1 do location = 1, n if ( p(location) == value ) then ierror = 0 exit end if end do if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'PERM1_CHECK - Fatal error!' write ( *, '(a,i4)' ) ' Permutation is missing value ', value exit end if end do perm1_check = ierror return end subroutine r8col_compare ( m, n, a, i, j, value ) !*****************************************************************************80 ! !! R8COL_COMPARE compares columns in an R8COL. ! ! Discussion: ! ! An R8COL is an M by N array of R8's, regarded as an array of N columns, ! each of length M. ! ! Example: ! ! Input: ! ! M = 3, N = 4, I = 2, J = 4 ! ! A = ( ! 1. 2. 3. 4. ! 5. 6. 7. 8. ! 9. 10. 11. 12. ) ! ! Output: ! ! VALUE = -1 ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 05 December 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns. ! ! Input, real ( kind = rk ) A(M,N), the M by N array. ! ! Input, integer I, J, the columns to be compared. ! I and J must be between 1 and N. ! ! Output, integer VALUE, the results of the comparison: ! -1, column I < column J, ! 0, column I = column J, ! +1, column J < column I. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a(m,n) integer i integer j integer k integer value ! ! Check. ! if ( i < 1 .or. n < i ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8COL_COMPARE - Fatal error!' write ( *, '(a)' ) ' Column index I is out of bounds.' write ( *, '(a,i8)' ) ' I = ', i stop 1 end if if ( j < 1 .or. n < j ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8COL_COMPARE - Fatal error!' write ( *, '(a)' ) ' Column index J is out of bounds.' write ( *, '(a,i8)' ) ' J = ', j stop 1 end if value = 0 if ( i == j ) then return end if k = 1 do while ( k <= m ) if ( a(k,i) < a(k,j) ) then value = -1 return else if ( a(k,j) < a(k,i) ) then value = +1 return end if k = k + 1 end do return end subroutine r8col_duplicates ( m, n, n_unique, a ) !*****************************************************************************80 ! !! R8COL_DUPLICATES generates an R8COL with some duplicate columns. ! ! Discussion: ! ! An R8COL is an M by N array of R8's, regarded as an array of N columns, ! each of length M. ! ! This routine generates a random R8COL with a specified number of ! duplicate columns. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 21 July 2010 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, the number of rows in each column of A. ! ! Input, integer N, the number of columns in A. ! ! Input, integer N_UNIQUE, the number of unique columns in A. ! 1 <= N_UNIQUE <= N. ! ! Output, real ( kind = rk ) A(M,N), the array. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a(m,n) integer i4_uniform_ab integer j1 integer j2 integer n_unique real ( kind = rk ) temp(m) if ( n_unique < 1 .or. n < n_unique ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8COL_DUPLICATES - Fatal error!' write ( *, '(a)' ) ' 1 <= N_UNIQUE <= N is required.' stop 1 end if call random_number ( harvest = a(1:m,1:n_unique) ) ! ! Randomly copy unique columns. ! do j1 = n_unique + 1, n j2 = i4_uniform_ab ( 1, n_unique ) a(1:m,j1) = a(1:m,j2) end do ! ! Permute the columns. ! do j1 = 1, n j2 = i4_uniform_ab ( j1, n ) temp(1:m) = a(1:m,j1) a(1:m,j1) = a(1:m,j2) a(1:m,j2) = temp(1:m) end do return end subroutine r8col_find ( m, n, a, x, col ) !*****************************************************************************80 ! !! R8COL_FIND seeks a column value in an R8COL. ! ! Discussion: ! ! An R8COL is an M by N array of R8's, regarded as an array of N columns, ! each of length M. ! ! Example: ! ! Input: ! ! M = 3, ! N = 4, ! ! A = ( ! 1. 2. 3. 4. ! 5. 6. 7. 8. ! 9. 10. 11. 12. ) ! ! x = ( 3., ! 7., ! 11. ) ! ! Output: ! ! COL = 3 ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 05 December 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns. ! ! Input, real ( kind = rk ) A(M,N), a table of numbers, regarded as ! N columns of vectors of length M. ! ! Input, real ( kind = rk ) X(M), a vector to be matched with a column of A. ! ! Output, integer COL, the index of the first column of A ! which exactly matches every entry of X, or -1 if no match ! could be found. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a(m,n) integer col integer i integer j real ( kind = rk ) x(m) col = -1 do j = 1, n col = j do i = 1, m if ( x(i) /= a(i,j) ) then col = -1 exit end if end do if ( col /= -1 ) then return end if end do return end subroutine r8col_first_index ( m, n, a, tol, first_index ) !*****************************************************************************80 ! !! R8COL_FIRST_INDEX indexes the first occurrence of values in an R8COL. ! ! Discussion: ! ! An R8COL is an M by N array of R8's, regarded as an array of N columns, ! each of length M. ! ! For element A(1:M,J) of the matrix, FIRST_INDEX(J) is the index in A of ! the first column whose entries are equal to A(1:M,J). ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 24 November 2008 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns of A. ! The length of an "element" of A, and the number of "elements". ! ! Input, real ( kind = rk ) A(M,N), the array. ! ! Input, real ( kind = rk ) TOL, a tolerance for equality. ! ! Output, integer FIRST_INDEX(N), the first occurrence index. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a(m,n) integer first_index(n) integer j1 integer j2 real ( kind = rk ) tol first_index(1:n) = -1 do j1 = 1, n if ( first_index(j1) == -1 ) then first_index(j1) = j1 do j2 = j1 + 1, n if ( maxval ( abs ( a(1:m,j1) - a(1:m,j2) ) ) <= tol ) then first_index(j2) = j1 end if end do end if end do return end subroutine r8col_flip ( m, n, a ) !*****************************************************************************80 ! !! R8COL_FLIP flips the entries in each column of an R8COL. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 05 May 2017 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns. ! ! Input/output, real ( kind = rk ) A(M,N), the array to be flipped. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a(m,n) integer i integer ihi integer j real ( kind = rk ) t ihi = m / 2 do j = 1, n do i = 1, ihi t = a(i,j) a(i,j) = a(m+1-i,j) a(m+1-j,j) = t end do end do return end subroutine r8col_indicator ( m, n, table ) !*****************************************************************************80 ! !! R8COL_INDICATOR sets up an "indicator" R8COL. ! ! Discussion: ! ! An R8COL is an M by N array of R8's, regarded as an array of N columns, ! each of length M. ! ! The value of each entry suggests its location, as in: ! ! 11 12 13 14 ! 21 22 23 24 ! 31 32 33 34 ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 28 May 2008 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, the number of rows of the matrix. ! M must be positive. ! ! Input, integer N, the number of columns of the matrix. ! N must be positive. ! ! Output, real ( kind = rk ) TABLE(M,N), the table. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n integer fac integer i integer i4_log_10 integer j real ( kind = rk ) table(m,n) fac = 10 ** ( i4_log_10 ( n ) + 1 ) do i = 1, m do j = 1, n table(i,j) = real ( fac * i + j, kind = rk ) end do end do return end subroutine r8col_insert ( n_max, m, n, a, x, col ) !*****************************************************************************80 ! !! R8COL_INSERT inserts a column into an R8COL. ! ! Discussion: ! ! An R8COL is an M by N array of R8's, regarded as an array of N columns, ! each of length M. ! ! Example: ! ! Input: ! ! N_MAX = 10, ! M = 3, ! N = 4, ! ! A = ( ! 1. 2. 3. 4. ! 5. 6. 7. 8. ! 9. 10. 11. 12. ) ! ! X = ( 3., 4., 18. ) ! ! Output: ! ! N = 5, ! ! A = ( ! 1. 2. 3. 3. 4. ! 5. 6. 4. 7. 8. ! 9. 10. 18. 11. 12. ) ! ! COL = 3 ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 05 December 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N_MAX, the maximum number of columns in A. ! ! Input, integer M, the number of rows. ! ! Input/output, integer N, the number of columns. ! If the new column is inserted into the table, then the output ! value of N will be increased by 1. ! ! Input/output, real ( kind = rk ) A(M,N_MAX), a table of numbers, regarded ! as an array of columns. The columns must have been sorted ! lexicographically. ! ! Input, real ( kind = rk ) X(M), a vector of data which will be inserted ! into the table if it does not already occur. ! ! Output, integer COL. ! I, X was inserted into column I. ! -I, column I was already equal to X. ! 0, N = N_MAX. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n_max real ( kind = rk ) a(m,n_max) integer col integer high integer isgn integer j integer low integer mid integer n real ( kind = rk ) x(m) ! ! Refuse to work if N_MAX <= N. ! if ( n_max <= n ) then col = 0 return end if ! ! Stick X temporarily in column N+1, just so it's easy to use R8COL_COMPARE. ! a(1:m,n+1) = x(1:m) ! ! Do a binary search. ! low = 1 high = n do if ( high < low ) then col = low exit end if mid = ( low + high ) / 2 call r8col_compare ( m, n + 1, a, mid, n + 1, isgn ) if ( isgn == 0 ) then col = -mid return else if ( isgn == -1 ) then low = mid + 1 else if ( isgn == +1 ) then high = mid - 1 end if end do ! ! Shift part of the table up to make room. ! do j = n, col, -1 a(1:m,j+1) = a(1:m,j) end do ! ! Insert the new column. ! a(1:m,col) = x(1:m) n = n + 1 return end subroutine r8col_max ( m, n, a, amax ) !*****************************************************************************80 ! !! R8COL_MAX returns the maximums in an R8COL. ! ! Discussion: ! ! An R8COL is an M by N array of R8's, regarded as an array of N columns, ! each of length M. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 15 September 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns. ! ! Input, real ( kind = rk ) A(M,N), the array to be examined. ! ! Output, real ( kind = rk ) AMAX(N), the maximums of the columns. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a(m,n) real ( kind = rk ) amax(n) integer j do j = 1, n amax(j) = maxval ( a(1:m,j) ) end do return end subroutine r8col_max_index ( m, n, a, imax ) !*****************************************************************************80 ! !! R8COL_MAX_INDEX returns the indices of column maximums in an R8COL. ! ! Discussion: ! ! An R8COL is an M by N array of R8's, regarded as an array of N columns, ! each of length M. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 15 September 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns. ! ! Input, real ( kind = rk ) A(M,N), the array to be examined. ! ! Output, integer IMAX(N); IMAX(I) is the row of A in which ! the maximum for column I occurs. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a(m,n) real ( kind = rk ) amax integer i integer imax(n) integer j do j = 1, n imax(j) = 1 amax = a(1,j) do i = 2, m if ( amax < a(i,j) ) then imax(j) = i amax = a(i,j) end if end do end do return end subroutine r8col_max_one ( m, n, a ) !*****************************************************************************80 ! !! R8COL_MAX_ONE rescales an R8COL so each column maximum is 1. ! ! Discussion: ! ! An R8COL is an M by N array of R8's, regarded as an array of N columns, ! each of length M. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 08 May 2010 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns. ! ! Input/output, real ( kind = rk ) A(M,N), the array to be rescaled. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a(m,n) integer i integer i_big integer j do j = 1, n i_big = 1 do i = 2, m if ( abs ( a(i_big,j) ) < abs ( a(i,j) ) ) then i_big = i end if end do if ( a(i_big,j) /= 0.0D+00 ) then a(1:m,j) = a(1:m,j) / a(i_big,j) end if end do return end subroutine r8col_mean ( m, n, a, mean ) !*****************************************************************************80 ! !! R8COL_MEAN returns the column means of an R8COL. ! ! Discussion: ! ! An R8COL is an M by N array of R8's, regarded as an array of N columns, ! each of length M. ! ! Example: ! ! A = ! 1 2 3 ! 2 6 7 ! ! MEAN = ! 1.5 4.0 5.0 ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 05 December 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns. ! ! Input, real ( kind = rk ) A(M,N), the array to be examined. ! ! Output, real ( kind = rk ) MEAN(N), the means, or averages, of the columns. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a(m,n) integer j real ( kind = rk ) mean(n) do j = 1, n mean(j) = sum ( a(1:m,j) ) end do mean(1:n) = mean(1:n) / real ( m, kind = rk ) return end subroutine r8col_min ( m, n, a, amin ) !*****************************************************************************80 ! !! R8COL_MIN returns the column minimums of an R8COL. ! ! Discussion: ! ! An R8COL is an M by N array of R8's, regarded as an array of N columns, ! each of length M. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 05 December 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns. ! ! Input, real ( kind = rk ) A(M,N), the array to be examined. ! ! Output, real ( kind = rk ) AMIN(N), the minimums of the columns. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a(m,n) real ( kind = rk ) amin(n) integer j do j = 1, n amin(j) = minval ( a(1:m,j) ) end do return end subroutine r8col_min_index ( m, n, a, imin ) !*****************************************************************************80 ! !! R8COL_MIN_INDEX returns the indices of column minimums in an R8COL. ! ! Discussion: ! ! An R8COL is an M by N array of R8's, regarded as an array of N columns, ! each of length M. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 15 September 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns. ! ! Input, real ( kind = rk ) A(M,N), the array to be examined. ! ! Output, integer IMIN(N); IMIN(I) is the row of A in which ! the minimum for column I occurs. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a(m,n) real ( kind = rk ) amin integer i integer imin(n) integer j do j = 1, n imin(j) = 1 amin = a(1,j) do i = 2, m if ( a(i,j) < amin ) then imin(j) = i amin = a(i,j) end if end do end do return end subroutine r8col_normalize_li ( m, n, a ) !*****************************************************************************80 ! !! R8COL_NORMALIZE_LI normalizes an R8COL with the column infinity norm. ! ! Discussion: ! ! Each column is scaled so that the entry of maximum norm has the value 1. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 08 February 2012 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns. ! ! Input/output, real ( kind = rk ) A(M,N), the array to be normalized. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a(m,n) real ( kind = rk ) c integer i integer j do j = 1, n c = a(1,j) do i = 2, m if ( abs ( c ) < abs ( a(i,j) ) ) then c = a(i,j) end if end do if ( c /= 0.0D+00 ) then a(1:m,j) = a(1:m,j) / c end if end do return end subroutine r8col_part_quick_a ( m, n, a, l, r ) !*****************************************************************************80 ! !! R8COL_PART_QUICK_A reorders the columns of an R8COL. ! ! Discussion: ! ! An R8COL is an M by N array of R8's, regarded as an array of N columns, ! each of length M. ! ! The routine reorders the columns of A. Using A(1:M,1) as a ! key, all entries of A that are less than or equal to the key will ! precede the key, which precedes all entries that are greater than the key. ! ! Example: ! ! Input: ! ! M = 2, N = 8 ! A = ( 2 8 6 0 10 10 0 5 ! 4 8 2 2 6 0 6 8 ) ! ! Output: ! ! L = 2, R = 4 ! ! A = ( 0 0 2 8 6 10 10 5 ! 2 6 4 8 2 6 0 8 ) ! ---- ------------- ! LEFT KEY RIGHT ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 05 December 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, the row dimension of A, and the length of ! a column. ! ! Input, integer N, the column dimension of A. ! ! Input/output, real ( kind = rk ) A(M,N). On input, the array to be checked. ! On output, A has been reordered as described above. ! ! Output, integer L, R, the indices of A that define the three ! segments. Let KEY = the input value of A(1:M,1). Then ! I <= L A(1:M,I) < KEY; ! L < I < R A(1:M,I) = KEY; ! R <= I KEY < A(1:M,I). ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a(m,n) integer j integer k real ( kind = rk ) key(m) integer l integer r logical ( kind = 4 ) r8vec_eq logical ( kind = 4 ) r8vec_gt logical ( kind = 4 ) r8vec_lt if ( n < 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8COL_PART_QUICK_A - Fatal error!' write ( *, '(a)' ) ' N < 1.' return end if if ( n == 1 ) then l = 0 r = 2 return end if key(1:m) = a(1:m,1) k = 1 ! ! The elements of unknown size have indices between L+1 and R-1. ! l = 1 r = n + 1 do j = 2, n if ( r8vec_gt ( m, a(1:m,l+1), key(1:m) ) ) then r = r - 1 call r8vec_swap ( m, a(1:m,r), a(1:m,l+1) ) else if ( r8vec_eq ( m, a(1:m,l+1), key(1:m) ) ) then k = k + 1 call r8vec_swap ( m, a(1:m,k), a(1:m,l+1) ) l = l + 1 else if ( r8vec_lt ( m, a(1:m,l+1), key(1:m) ) ) then l = l + 1 end if end do ! ! Shift small elements to the left. ! do j = 1, l - k a(1:m,j) = a(1:m,j+k) end do ! ! Shift KEY elements to center. ! do j = l - k + 1, l a(1:m,j) = key(1:m) end do ! ! Update L. ! l = l - k return end subroutine r8col_permute ( m, n, p, a ) !*****************************************************************************80 ! !! R8COL_PERMUTE permutes an R8COL in place. ! ! Discussion: ! ! An R8COL is an M by N array of R8's, regarded as an array of N columns, ! each of length M. ! ! The same logic can be used to permute an array of objects of any ! arithmetic type, or an array of objects of any complexity. The only ! temporary storage required is enough to store a single object. The number ! of data movements made is N + the number of cycles of order 2 or more, ! which is never more than N + N/2. ! ! Example: ! ! Input: ! ! M = 2 ! N = 5 ! P = ( 2, 4, 5, 1, 3 ) ! A = ( 1.0, 2.0, 3.0, 4.0, 5.0 ) ! (11.0, 22.0, 33.0, 44.0, 55.0 ) ! ! Output: ! ! A = ( 2.0, 4.0, 5.0, 1.0, 3.0 ) ! ( 22.0, 44.0, 55.0, 11.0, 33.0 ). ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 09 December 2006 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, the dimension of objects. ! ! Input, integer N, the number of objects. ! ! Input, integer P(N), the permutation. P(I) = J means ! that the I-th element of the output array should be the J-th ! element of the input array. ! ! Input/output, real ( kind = rk ) A(M,N), the array to be permuted. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a(m,n) real ( kind = rk ) a_temp(m) integer ierror integer iget integer iput integer istart integer p(n) integer perm1_check ierror = perm1_check ( n, p ) if ( ierror .ne. 0 ) then write ( *, '(a)' ) '' write ( *, '(a)' ) 'R8COL_PERMUTE - Fatal error!' write ( *, '(a)' ) ' PERM1_CHECK returned error.' stop 1 end if ! ! Search for the next element of the permutation that has not been used. ! do istart = 1, n if ( p(istart) < 0 ) then cycle else if ( p(istart) == istart ) then p(istart) = - p(istart) cycle else a_temp(1:m) = a(1:m,istart) iget = istart ! ! Copy the new value into the vacated entry. ! do iput = iget iget = p(iget) p(iput) = - p(iput) if ( iget < 1 .or. n < iget ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8COL_PERMUTE - Fatal error!' write ( *, '(a)' ) ' A permutation index is out of range.' write ( *, '(a,i8,a,i8)' ) ' P(', iput, ') = ', iget stop 1 end if if ( iget == istart ) then a(1:m,iput) = a_temp(1:m) exit end if a(1:m,iput) = a(1:m,iget) end do end if end do ! ! Restore the signs of the entries. ! p(1:n) = - p(1:n) return end subroutine r8col_print ( m, n, a, title ) !*****************************************************************************80 ! !! R8COL_PRINT prints an R8COL. ! ! Discussion: ! ! An R8COL is an M by N array of R8's, regarded as an array of N columns, ! each of length M. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 12 September 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, the number of rows in A. ! ! Input, integer N, the number of columns in A. ! ! Input, real ( kind = rk ) A(M,N), the matrix. ! ! Input, character ( len = * ) TITLE, a title. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a(m,n) character ( len = * ) title call r8col_print_some ( m, n, a, 1, 1, m, n, title ) return end subroutine r8col_print_some ( m, n, a, ilo, jlo, ihi, jhi, title ) !*****************************************************************************80 ! !! R8COL_PRINT_SOME prints some of an R8COL. ! ! Discussion: ! ! An R8COL is an M by N array of R8's, regarded as an array of N columns, ! each of length M. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 10 September 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns. ! ! Input, real ( kind = rk ) A(M,N), an M by N matrix to be printed. ! ! Input, integer ILO, JLO, the first row and column to print. ! ! Input, integer IHI, JHI, the last row and column to print. ! ! Input, character ( len = * ) TITLE, a 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 write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) if ( m <= 0 .or. n <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' (None)' return 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), '(i8,6x)' ) 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,a,5a14)' ) i, ':', ( ctemp(j), j = 1, inc ) end do end do return end subroutine r8col_reverse ( m, n, a ) !*****************************************************************************80 ! !! R8COL_REVERSE reverses the order of columns in an R8COL. ! ! Discussion: ! ! To reverse the columns is to start with something like ! ! 11 12 13 14 15 ! 21 22 23 24 25 ! 31 32 33 34 35 ! 41 42 43 44 45 ! 51 52 53 54 55 ! ! and return ! ! 15 14 13 12 11 ! 25 24 23 22 21 ! 35 34 33 32 31 ! 45 44 43 42 41 ! 55 54 53 52 51 ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 06 May 2013 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns. ! ! Input/output, real ( kind = rk ) A(M,N), the matrix. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a(m,n) integer j integer jhi real ( kind = rk ) t(m) jhi = n / 2 do j = 1, jhi t(1:m) = a(1:m,j) a(1:m,j) = a(1:m,n+1-j) a(1:m,n+1-j) = t(1:m) end do return end subroutine r8col_separation ( m, n, a, d_min, d_max ) !*****************************************************************************80 ! !! R8COL_SEPARATION returns the "separation" of an R8COL. ! ! Discussion: ! ! D_MIN is the minimum distance between two columns, ! D_MAX is the maximum distance between two columns. ! ! The distances are measured using the Loo norm. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 24 February 2014 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns ! in the array. If N < 2, it does not make sense to call this routine. ! ! Input, real ( kind = rk ) A(M,N), the array whose variances are desired. ! ! Output, real ( kind = rk ) D_MIN, D_MAX, the minimum and maximum distances. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a(m,n) real ( kind = rk ) d real ( kind = rk ) d_max real ( kind = rk ) d_min integer j1 integer j2 d_min = huge ( d_min ) d_max = 0.0D+00 do j1 = 1, n do j2 = j1 + 1, n d = maxval ( abs ( a(1:m,j1) - a(1:m,j2) ) ) d_min = min ( d_min, d ) d_max = max ( d_max, d ) end do end do return end subroutine r8col_sort_heap_a ( m, n, a ) !*****************************************************************************80 ! !! R8COL_SORT_HEAP_A ascending heapsorts an R8COL. ! ! Discussion: ! ! An R8COL is an M by N array of R8's, regarded as an array of N columns, ! each of length M. ! ! In lexicographic order, the statement "X < Y", applied to two real ! 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, the first time they differ, X is smaller. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 05 December 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns. ! ! Input/output, real ( kind = rk ) A(M,N). ! On input, the array of N columns of M-vectors. ! On output, the columns of A have been sorted in lexicographic order. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a(m,n) integer i integer indx integer isgn integer j1 integer j2 real ( kind = rk ) t(m) if ( m <= 0 ) then return end if if ( n <= 1 ) then return end if ! ! Initialize. ! indx = 0 isgn = 0 j1 = 0 j2 = 0 ! ! Call the external heap sorter. ! do call sort_heap_external ( n, indx, j1, j2, isgn ) ! ! Interchange columns J1 and J2. ! if ( 0 < indx ) then t(1:m) = a(1:m,j1) a(1:m,j1) = a(1:m,j2) a(1:m,j2) = t(1:m) ! ! Compare columns J1 and J2. ! else if ( indx < 0 ) then isgn = 0 do i = 1, m if ( a(i,j1) < a(i,j2) ) then isgn = -1 exit else if ( a(i,j2) < a(i,j1) ) then isgn = +1 exit end if end do ! ! The columns are sorted. ! else if ( indx == 0 ) then exit end if end do return end subroutine r8col_sort_heap_index_a ( m, n, a, indx ) !*****************************************************************************80 ! !! R8COL_SORT_HEAP_INDEX_A does an indexed heap ascending sort of an R8COL. ! ! Discussion: ! ! An R8COL is an M by N array of R8's, regarded as an array of N columns, ! each of length M. ! ! The sorting is not actually carried out. Rather an index array is ! created which defines the sorting. This array may be used to sort ! or index the array, or to sort or index related arrays keyed on the ! original array. ! ! A(*,J1) < A(*,J2) if the first nonzero entry of A(*,J1)-A(*,J2) ! is negative. ! ! Once the index array is computed, the sorting can be carried out ! "implicitly: ! ! A(*,INDX(*)) is sorted, ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 27 October 2008 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, the number of rows in each column of A. ! ! Input, integer N, the number of columns in A. ! ! Input, real ( kind = rk ) A(M,N), the array. ! ! Output, integer INDX(N), the sort index. The I-th element ! of the sorted array is column INDX(I). ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a(m,n) real ( kind = rk ) column(m) integer i integer indx(n) integer indxt integer ir integer isgn integer j integer l if ( n < 1 ) then return end if do i = 1, n indx(i) = i end do if ( n == 1 ) then return end if l = ( n / 2 ) + 1 ir = n do if ( 1 < l ) then l = l - 1 indxt = indx(l) column(1:m) = a(1:m,indxt) else indxt = indx(ir) column(1:m) = a(1:m,indxt) indx(ir) = indx(1) ir = ir - 1 if ( ir == 1 ) then indx(1) = indxt exit end if end if i = l j = l + l do while ( j <= ir ) if ( j < ir ) then call r8vec_compare ( m, a(1:m,indx(j)), a(1:m,indx(j+1)), isgn ) if ( isgn < 0 ) then j = j + 1 end if end if call r8vec_compare ( m, column, a(1:m,indx(j)), isgn ) if ( isgn < 0 ) then indx(i) = indx(j) i = j j = j + j else j = ir + 1 end if end do indx(i) = indxt end do return end subroutine r8col_sort_quick_a ( m, n, a ) !*****************************************************************************80 ! !! R8COL_SORT_QUICK_A ascending quick sorts an R8COL. ! ! Discussion: ! ! An R8COL is an M by N array of R8's, regarded as an array of N columns, ! each of length M. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 21 May 2012 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, the row order of A, and the length of ! a column. ! ! Input, integer N, the number of columns of A. ! ! Input/output, real ( kind = rk ) A(M,N). ! On input, the array to be sorted. ! On output, the array has been sorted. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: level_max = 30 integer m integer n real ( kind = rk ) a(m,n) integer base integer l_segment integer level integer n_segment integer rsave(level_max) integer r_segment if ( m <= 0 ) then return end if if ( n < 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8COL_SORT_QUICK_A - Fatal error!' write ( *, '(a)' ) ' N < 1.' write ( *, '(a,i8)' ) ' N = ', n stop 1 end if if ( n == 1 ) then return end if level = 1 rsave(level) = n + 1 base = 1 n_segment = n do ! ! Partition the segment. ! call r8col_part_quick_a ( m, n_segment, a(1:m,base:base+n_segment-1), & l_segment, r_segment ) ! ! If the left segment has more than one element, we need to partition it. ! if ( 1 < l_segment ) then if ( level_max < level ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8COL_SORT_QUICK_A - Fatal error!' write ( *, '(a,i8)' ) ' Exceeding recursion maximum of ', level_max stop 1 end if level = level + 1 n_segment = l_segment rsave(level) = r_segment + base - 1 ! ! The left segment and the middle segment are sorted. ! Must the right segment be partitioned? ! else if ( r_segment < n_segment ) then n_segment = n_segment + 1 - r_segment base = base + r_segment - 1 ! ! Otherwise, we back up a level if there is an earlier one. ! else do if ( level <= 1 ) then return end if base = rsave(level) n_segment = rsave(level-1) - rsave(level) level = level - 1 if ( 0 < n_segment ) then exit end if end do end if end do return end subroutine r8col_sorted_tol_undex ( m, n, a, unique_num, tol, undx, xdnu ) !*****************************************************************************80 ! !! R8COL_SORTED_TOL_UNDEX indexes tolerably unique entries in a sorted R8COL. ! ! Discussion: ! ! An R8COL is an M by N array of R8's, regarded as an array of N columns, ! each of length M. ! ! The goal of this routine is to determine a vector UNDX, ! which points, to the tolerably unique elements of A, in sorted order, ! and a vector XDNU, which identifies, for each entry of A, the index of ! the unique sorted element of A. ! ! This is all done with index vectors, so that the elements of ! A are never moved. ! ! Assuming A is already sorted, we examine the entries of A in order, ! noting the unique entries, creating the entries of XDNU and ! UNDX as we go. ! ! Once this process has been completed, the vector A could be ! replaced by a compressed vector XU, containing the unique entries ! of A in sorted order, using the formula ! ! XU(*) = A(UNDX(*)). ! ! We could then, if we wished, reconstruct the entire vector A, or ! any element of it, by index, as follows: ! ! A(I) = XU(XDNU(I)). ! ! We could then replace A by the combination of XU and XDNU. ! ! Later, when we need the I-th entry of A, we can locate it as ! the XDNU(I)-th entry of XU. ! ! Here is an example of a vector A, the unique sort and ! inverse unique sort vectors and the compressed unique sorted vector. ! ! I A XU Undx Xdnu ! ----+------+------+-----+-----+ ! 1 | 11.0 | 11.0 1 1 ! 2 | 11.0 | 22.0 5 1 ! 3 | 11.0 | 33.0 8 1 ! 4 | 11.0 | 55.0 9 1 ! 5 | 22.0 | 2 ! 6 | 22.0 | 2 ! 7 | 22.0 | 2 ! 8 | 33.0 | 3 ! 9 | 55.0 | 4 ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 17 July 2010 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, the dimension of the data values. ! ! Input, integer N, the number of data values. ! ! Input, real ( kind = rk ) A(M,N), the data values. ! ! Input, integer UNIQUE_NUM, the number of unique values ! in A. This value is only required for languages in which the size of ! UNDX must be known in advance. ! ! Input, real ( kind = rk ) TOL, a tolerance for equality. ! ! Output, integer UNDX(UNIQUE_NUM), the UNDX vector. ! ! Output, integer XDNU(N), the XDNU vector. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n integer unique_num real ( kind = rk ) a(m,n) real ( kind = rk ) diff integer i integer i2 integer j integer k real ( kind = rk ) tol integer undx(unique_num) logical ( kind = 4 ) unique integer xdnu(n) ! ! Consider entry I = 1. ! It is unique, so set the number of unique items to K. ! Set the K-th unique item to I. ! Set the representative of item I to the K-th unique item. ! i = 1 k = 1 undx(k) = i xdnu(i) = k ! ! Consider entry I. ! ! If it is unique, increase the unique count K, set the ! K-th unique item to I, and set the representative of I to K. ! ! If it is not unique, set the representative of item I to a ! previously determined unique item that is close to it. ! do i = 2, n unique = .true. do j = 1, k i2 = undx(j) diff = maxval ( abs ( a(1:m,i) - a(1:m,i2) ) ) if ( diff <= tol ) then unique = .false. xdnu(i) = j exit end if end do if ( unique ) then k = k + 1 undx(k) = i xdnu(i) = k end if end do return end subroutine r8col_sorted_tol_unique ( m, n, a, tol, unique_num ) !*****************************************************************************80 ! !! R8COL_SORTED_TOL_UNIQUE keeps tolerably unique elements in a sorted R8COL. ! ! Discussion: ! ! An R8COL is an M by N array of R8's, regarded as an array of N columns, ! each of length M. ! ! The columns of the array may be ascending or descending sorted. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 July 2010 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns. ! ! Input/output, real ( kind = rk ) A(M,N). ! On input, the sorted array of N columns of M-vectors. ! On output, a sorted array of columns of M-vectors. ! ! Input, real ( kind = rk ) TOL, a tolerance for equality. ! ! Output, integer UNIQUE_NUM, the number of unique columns. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a(m,n) real ( kind = rk ) diff integer i integer j real ( kind = rk ) tol logical ( kind = 4 ) unique integer unique_num if ( n <= 0 ) then unique_num = 0 return end if unique_num = 1 do i = 2, n unique = .true. do j = 1, unique_num diff = maxval ( abs ( a(1:m,j) - a(1:m,i) ) ) if ( diff <= tol ) then unique = .false. exit end if end do if ( unique ) then unique_num = unique_num + 1 a(1:m,unique_num) = a(1:m,i) end if end do return end subroutine r8col_sorted_tol_unique_count ( m, n, a, tol, unique_num ) !*****************************************************************************80 ! !! R8COL_SORTED_TOL_UNIQUE_COUNT: tolerably unique elements in a sorted R8COL. ! ! Discussion: ! ! An R8COL is an M by N array of R8's, regarded as an array of N columns, ! each of length M. ! ! The columns of the array may be ascending or descending sorted. ! ! If the tolerance is large enough, then the concept of uniqueness ! can become ambiguous. If we have a tolerance of 1.5, then in the ! list ( 1, 2, 3, 4, 5, 6, 7, 8, 9 ) is it fair to say we have only ! one unique entry? That would be because 1 may be regarded as unique, ! and then 2 is too close to 1 to be unique, and 3 is too close to 2 to ! be unique and so on. ! ! This seems wrongheaded. So I prefer the idea that an item is not ! unique under a tolerance only if it is close to something that IS unique. ! Thus, the unique items are guaranteed to cover the space if we include ! a disk of radius TOL around each one. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 19 July 2010 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns. ! ! Input, real ( kind = rk ) A(M,N), a sorted array, containing ! N columns of data. ! ! Input, real ( kind = rk ) TOL, a tolerance for equality. ! ! Output, integer UNIQUE_NUM, the number of unique columns. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a(m,n) real ( kind = rk ) diff integer i integer i2 integer j integer k real ( kind = rk ) tol integer undx(n) logical ( kind = 4 ) unique integer unique_num ! ! Consider entry I = 1. ! It is unique, so set the number of unique items to K. ! Set the K-th unique item to I. ! i = 1 k = 1 undx(k) = i ! ! Consider entry I. ! ! If it is unique, increase the unique count K and set the ! K-th unique item to I. ! do i = 2, n unique = .true. do j = 1, k i2 = undx(j) diff = maxval ( abs ( a(1:m,i) - a(1:m,i2) ) ) if ( diff <= tol ) then unique = .false. exit end if end do if ( unique ) then k = k + 1 undx(k) = i end if end do unique_num = k return end subroutine r8col_sorted_undex ( m, n, a, unique_num, undx, xdnu ) !*****************************************************************************80 ! !! R8COL_SORTED_UNDEX returns unique sorted indexes for a sorted R8COL. ! ! Discussion: ! ! An R8COL is an M by N array of R8's, regarded as an array of N columns, ! each of length M. ! ! The goal of this routine is to determine a vector UNDX, ! which points, to the unique elements of A, in sorted order, ! and a vector XDNU, which identifies, for each entry of A, the index of ! the unique sorted element of A. ! ! This is all done with index vectors, so that the elements of ! A are never moved. ! ! Assuming A is already sorted, we examine the entries of A in order, ! noting the unique entries, creating the entries of XDNU and ! UNDX as we go. ! ! Once this process has been completed, the vector A could be ! replaced by a compressed vector XU, containing the unique entries ! of A in sorted order, using the formula ! ! XU(*) = A(UNDX(*)). ! ! We could then, if we wished, reconstruct the entire vector A, or ! any element of it, by index, as follows: ! ! A(I) = XU(XDNU(I)). ! ! We could then replace A by the combination of XU and XDNU. ! ! Later, when we need the I-th entry of A, we can locate it as ! the XDNU(I)-th entry of XU. ! ! Here is an example of a vector A, the sort and inverse sort ! index vectors, and the unique sort and inverse unique sort vectors ! and the compressed unique sorted vector. ! ! I A XU Undx Xdnu ! ----+------+------+-----+-----+ ! 1 | 11.0 | 11.0 1 1 ! 2 | 11.0 | 22.0 5 1 ! 3 | 11.0 | 33.0 8 1 ! 4 | 11.0 | 55.0 9 1 ! 5 | 22.0 | 2 ! 6 | 22.0 | 2 ! 7 | 22.0 | 2 ! 8 | 33.0 | 3 ! 9 | 55.0 | 4 ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 17 July 2010 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, the dimension of the data values. ! ! Input, integer N, the number of data values. ! ! Input, real ( kind = rk ) AL(M,N), the data values. ! ! Input, integer UNIQUE_NUM, the number of unique values ! in A. This value is only required for languages in which the size of ! UNDX must be known in advance. ! ! Output, integer UNDX(UNIQUE_NUM), the UNDX vector. ! ! Output, integer XDNU(N), the XDNU vector. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n integer unique_num real ( kind = rk ) a(m,n) integer i integer j integer undx(unique_num) integer xdnu(n) ! ! Walk through the sorted array. ! i = 1 j = 1 undx(j) = i xdnu(i) = j do i = 2, n if ( any ( a(1:m,i) /= a(1:m,j) ) ) then j = j + 1 undx(j) = i end if xdnu(i) = j end do return end subroutine r8col_sorted_unique ( m, n, a, unique_num ) !*****************************************************************************80 ! !! R8COL_SORTED_UNIQUE keeps unique elements in a sorted R8COL. ! ! Discussion: ! ! An R8COL is an M by N array of R8's, regarded as an array of N columns, ! each of length M. ! ! The columns of the array may be ascending or descending sorted. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 July 2010 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns. ! ! Input/output, real ( kind = rk ) A(M,N). ! On input, the sorted array of N columns of M-vectors. ! On output, a sorted array of columns of M-vectors. ! ! Output, integer UNIQUE_NUM, the number of unique columns. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a(m,n) integer j1 integer j2 integer unique_num if ( n <= 0 ) then unique_num = 0 return end if j1 = 1 do j2 = 2, n if ( any ( a(1:m,j1) /= a(1:m,j2) ) ) then j1 = j1 + 1 a(1:m,j1) = a(1:m,j2) end if end do unique_num = j1 return end subroutine r8col_sorted_unique_count ( m, n, a, unique_num ) !*****************************************************************************80 ! !! R8COL_SORTED_UNIQUE_COUNT counts unique elements in a sorted R8COL. ! ! Discussion: ! ! An R8COL is an M by N array of R8's, regarded as an array of N columns, ! each of length M. ! ! The columns of the array may be ascending or descending sorted. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 July 2010 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns. ! ! Input, real ( kind = rk ) A(M,N), a sorted array, containing ! N columns of data. ! ! Output, integer UNIQUE_NUM, the number of unique columns. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a(m,n) integer j1 integer j2 integer unique_num unique_num = 0 if ( n <= 0 ) then return end if unique_num = 1 j1 = 1 do j2 = 2, n if ( any ( a(1:m,j1) /= a(1:m,j2) ) ) then unique_num = unique_num + 1 j1 = j2 end if end do return end subroutine r8col_sortr_a ( m, n, a, key ) !*****************************************************************************80 ! !! R8COL_SORTR_A ascending sorts one column of an R8COL, adjusting all columns. ! ! Discussion: ! ! An R8COL is an M by N array of R8's, regarded as an array of N columns, ! each of length M. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 05 December 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns. ! ! Input/output, real ( kind = rk ) A(M,N). ! On input, an unsorted M by N array. ! On output, rows of the array have been shifted in such ! a way that column KEY of the array is in nondecreasing order. ! ! Input, integer KEY, the column in which the "key" value ! is stored. On output, column KEY of the array will be ! in nondecreasing order. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a(m,n) integer i integer indx integer isgn integer j integer key real ( kind = rk ) temp(n) if ( m <= 0 ) then return end if if ( key < 1 .or. n < key ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8COL_SORTR_A - Fatal error!' write ( *, '(a)' ) ' The value of KEY is not a legal column index.' write ( *, '(a,i8)' ) ' KEY = ', key write ( *, '(a,i8)' ) ' N = ', n stop 1 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 temp(1:n) = a(i,1:n) a(i,1:n) = a(j,1:n) a(j,1:n) = temp(1:n) ! ! Compare the I and J objects. ! else if ( indx < 0 ) then if ( a(i,key) < a(j,key) ) then isgn = -1 else isgn = +1 end if else if ( indx == 0 ) then exit end if end do return end subroutine r8col_sum ( m, n, a, colsum ) !*****************************************************************************80 ! !! R8COL_SUM sums the columns of an R8COL. ! ! Discussion: ! ! An R8COL is an M by N array of R8's, regarded as an array of N columns, ! each of length M. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 05 December 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns. ! ! Input, real ( kind = rk ) A(M,N), the array to be examined. ! ! Output, real ( kind = rk ) COLSUM(N), the sums of the columns. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a(m,n) real ( kind = rk ) colsum(n) integer j do j = 1, n colsum(j) = sum ( a(1:m,j) ) end do return end subroutine r8col_swap ( m, n, a, j1, j2 ) !*****************************************************************************80 ! !! R8COL_SWAP swaps columns I and J of an R8COL. ! ! Discussion: ! ! An R8COL is an M by N array of R8's, regarded as an array of N columns, ! each of length M. ! ! Example: ! ! Input: ! ! M = 3, N = 4, J1 = 2, J2 = 4 ! ! A = ( ! 1. 2. 3. 4. ! 5. 6. 7. 8. ! 9. 10. 11. 12. ) ! ! Output: ! ! A = ( ! 1. 4. 3. 2. ! 5. 8. 7. 6. ! 9. 12. 11. 10. ) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 05 December 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns. ! ! Input/output, real ( kind = rk ) A(M,N), the M by N array. ! ! Input, integer J1, J2, the columns to be swapped. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a(m,n) real ( kind = rk ) col(m) integer j1 integer j2 if ( j1 < 1 .or. n < j1 .or. j2 < 1 .or. n < j2 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8COL_SWAP - Fatal error!' write ( *, '(a)' ) ' J1 or J2 is out of bounds.' write ( *, '(a,i8)' ) ' J1 = ', j1 write ( *, '(a,i8)' ) ' J2 = ', j2 write ( *, '(a,i8)' ) ' NCOL = ', n stop 1 end if if ( j1 == j2 ) then return end if col(1:m) = a(1:m,j1) a(1:m,j1) = a(1:m,j2) a(1:m,j2) = col(1:m) return end subroutine r8col_to_r8vec ( m, n, a, x ) !*****************************************************************************80 ! !! R8COL_TO_R8VEC converts an R8COL to an R8VEC. ! ! Discussion: ! ! An R8COL is an M by N array of R8's, regarded as an array of N columns, ! each of length M. ! ! An R8VEC is a vector of R8's. ! ! Example: ! ! M = 3, N = 4 ! ! A = ! 11 12 13 14 ! 21 22 23 24 ! 31 32 33 34 ! ! X = ( 11, 21, 31, 12, 22, 32, 13, 23, 33, 14, 24, 34 ) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 05 December 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns. ! ! Input, real ( kind = rk ) A(M,N), the array. ! ! Output, real ( kind = rk ) X(M*N), a vector containing the N columns of A. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a(m,n) integer j integer k real ( kind = rk ) x(m*n) k = 1 do j = 1, n x(k:k+m-1) = a(1:m,j) k = k + m end do return end subroutine r8col_tol_undex ( m, n, a, unique_num, tol, undx, xdnu ) !*****************************************************************************80 ! !! R8COL_TOL_UNDEX indexes tolerably unique entries of an R8COL. ! ! Discussion: ! ! An R8COL is an M by N array of R8's, regarded as an array of N columns, ! each of length M. ! ! The goal of this routine is to determine a vector UNDX, ! which points, to the unique elements of A, in sorted order, ! and a vector XDNU, which identifies, for each entry of A, the index of ! the unique sorted element of A. ! ! This is all done with index vectors, so that the elements of ! A are never moved. ! ! The first step of the algorithm requires the indexed sorting ! of A, which creates arrays INDX and XDNI. (If all the entries ! of A are unique, then these arrays are the same as UNDX and XDNU.) ! ! We then use INDX to examine the entries of A in sorted order, ! noting the unique entries, creating the entries of XDNU and ! UNDX as we go. ! ! Once this process has been completed, the object X could be ! replaced by a compressed object XU, containing the unique entries ! of X in sorted order, using the formula ! ! XU(*) = A(UNDX(*)). ! ! We could then, if we wished, reconstruct the entire vector A, or ! any element of it, by index, as follows: ! ! A(I) = XU(XDNU(I)). ! ! We could then replace A by the combination of XU and XDNU. ! ! Later, when we need the I-th entry of A, we can locate it as ! the XDNU(I)-th entry of XU. ! ! Here is an example of a vector A, the sort and inverse sort ! index vectors, and the unique sort and inverse unique sort vectors ! and the compressed unique sorted vector. ! ! I A Indx Xdni XU Undx Xdnu ! ----+-----+-----+-----+--------+-----+-----+ ! 1 | 11. 1 1 | 11. 1 1 ! 2 | 22. 3 5 | 22. 2 2 ! 3 | 11. 6 2 | 33. 4 1 ! 4 | 33. 9 8 | 55. 5 3 ! 5 | 55. 2 9 | 4 ! 6 | 11. 7 3 | 1 ! 7 | 22. 8 6 | 2 ! 8 | 22. 4 7 | 2 ! 9 | 11. 5 4 | 1 ! ! INDX(2) = 3 means that sorted item(2) is A(3). ! XDNI(2) = 5 means that A(2) is sorted item(5). ! ! UNDX(3) = 4 means that unique sorted item(3) is at A(4). ! XDNU(8) = 2 means that A(8) is at unique sorted item(2). ! ! XU(XDNU(I))) = A(I). ! XU(I) = A(UNDX(I)). ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 17 July 2010 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, the dimension of the data values. ! ! Input, integer N, the number of data values. ! ! Input, real ( kind = rk ) A(M,N), the data values. ! ! Input, integer UNIQUE_NUM, the number of unique values ! in A. This value is only required for languages in which the size of ! UNDX must be known in advance. ! ! Input, real ( kind = rk ) TOL, a tolerance for equality. ! ! Output, integer UNDX(UNIQUE_NUM), the UNDX vector. ! ! Output, integer XDNU(N), the XDNU vector. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n integer unique_num real ( kind = rk ) a(m,n) real ( kind = rk ) diff integer i integer indx(n) integer j integer k real ( kind = rk ) tol integer undx(unique_num) logical ( kind = 4 ) unique integer xdnu(n) ! ! Implicitly sort the array. ! call r8col_sort_heap_index_a ( m, n, a, indx ) ! ! Consider entry I = 1. ! It is unique, so set the number of unique items to K. ! Set the K-th unique item to I. ! Set the representative of item I to the K-th unique item. ! i = 1 k = 1 undx(k) = indx(i) xdnu(indx(i)) = k ! ! Consider entry I. ! ! If it is unique, increase the unique count K, set the ! K-th unique item to I, and set the representative of I to K. ! ! If it is not unique, set the representative of item I to a ! previously determined unique item that is close to it. ! do i = 2, n unique = .true. do j = 1, k diff = maxval ( abs ( a(1:m,indx(i)) - a(1:m,undx(j)) ) ) if ( diff <= tol ) then unique = .false. xdnu(indx(i)) = j exit end if end do if ( unique ) then k = k + 1 undx(k) = indx(i) xdnu(indx(i)) = k end if end do return end subroutine r8col_tol_unique_count ( m, n, a, tol, unique_num ) !*****************************************************************************80 ! !! R8COL_TOL_UNIQUE_COUNT counts tolerably unique entries in an R8COL. ! ! Discussion: ! ! An R8COL is an M by N array of R8's, regarded as an array of N columns, ! each of length M. ! ! If the tolerance is large enough, then the concept of uniqueness ! can become ambiguous. If we have a tolerance of 1.5, then in the ! list ( 1, 2, 3, 4, 5, 6, 7, 8, 9 ) is it fair to say we have only ! one unique entry? That would be because 1 may be regarded as unique, ! and then 2 is too close to 1 to be unique, and 3 is too close to 2 to ! be unique and so on. ! ! This seems wrongheaded. So I prefer the idea that an item is not ! unique under a tolerance only if it is close to something that IS unique. ! Thus, the unique items are guaranteed to cover the space if we include ! a disk of radius TOL around each one. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 19 July 2010 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, the number of rows. ! ! Input, integer N, the number of columns. ! ! Input, real ( kind = rk ) A(M,N), the array of N columns of data. ! ! Input, real ( kind = rk ) TOL, a nonnegative tolerance for equality. ! ! Output, integer UNIQUE_NUM, the number of unique columns. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a(m,n) real ( kind = rk ) diff integer i integer indx(n) integer j integer k real ( kind = rk ) tol integer undx(n) logical ( kind = 4 ) unique integer unique_num ! ! Implicitly sort the array. ! call r8col_sort_heap_index_a ( m, n, a, indx ) ! ! Consider entry I = 1. ! It is unique, so set the number of unique items to K. ! Set the K-th unique item to I. ! Set the representative of item I to the K-th unique item. ! i = 1 k = 1 undx(k) = indx(i) ! ! Consider entry I. ! ! If it is unique, increase the unique count K, set the ! K-th unique item to I, and set the representative of I to K. ! ! If it is not unique, set the representative of item I to a ! previously determined unique item that is close to it. ! do i = 2, n unique = .true. do j = 1, k diff = maxval ( abs ( a(1:m,indx(i)) - a(1:m,undx(j)) ) ) if ( diff <= tol ) then unique = .false. exit end if end do if ( unique ) then k = k + 1 undx(k) = indx(i) end if end do unique_num = k return end subroutine r8col_tol_unique_index ( m, n, a, tol, unique_index ) !*****************************************************************************80 ! !! R8COL_TOL_UNIQUE_INDEX indexes tolerably unique entries in an R8COL. ! ! Discussion: ! ! An R8COL is an M by N array of R8's, regarded as an array of N columns, ! each of length M. ! ! For element A(1:M,J) of the matrix, UNIQUE_INDEX(J) is the uniqueness index ! of A(1:M,J). That is, if A_UNIQUE contains the unique elements of A, ! gathered in order, then ! ! A_UNIQUE ( 1:M, UNIQUE_INDEX(J) ) = A(1:M,J) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 17 July 2010 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns of A. ! ! Input, real ( kind = rk ) A(M,N), the array. ! ! Input, real ( kind = rk ) TOL, a tolerance for equality. ! ! Output, integer UNIQUE_INDEX(N), the first occurrence index. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a(m,n) real ( kind = rk ) diff integer j1 integer j2 real ( kind = rk ) tol integer unique_index(n) integer unique_num unique_index(1:n) = -1 unique_num = 0 do j1 = 1, n if ( unique_index(j1) == -1 ) then unique_num = unique_num + 1 unique_index(j1) = unique_num do j2 = j1 + 1, n diff = maxval ( abs ( a(1:m,j1) - a(1:m,j2) ) ) if ( diff <= tol ) then unique_index(j2) = unique_num end if end do end if end do return end subroutine r8col_transpose_print ( m, n, a, title ) !*****************************************************************************80 ! !! R8COL_TRANSPOSE_PRINT prints an R8COL, transposed. ! ! Discussion: ! ! An R8COL is an M by N array of R8's, regarded as an array of N columns, ! each of length M. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 14 June 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns. ! ! Input, real ( kind = rk ) A(M,N), an M by N matrix to be printed. ! ! Input, character ( len = * ) TITLE, a title. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a(m,n) character ( len = * ) title call r8col_transpose_print_some ( m, n, a, 1, 1, m, n, title ) return end subroutine r8col_transpose_print_some ( m, n, a, ilo, jlo, ihi, jhi, title ) !*****************************************************************************80 ! !! R8COL_TRANSPOSE_PRINT_SOME prints some of an R8COL, transposed. ! ! Discussion: ! ! An R8COL is an M by N array of R8's, regarded as an array of N columns, ! each of length M. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 10 September 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns. ! ! Input, real ( kind = rk ) A(M,N), an M by N matrix to be printed. ! ! Input, integer ILO, JLO, the first row and column to print. ! ! Input, integer IHI, JHI, the last row and column to print. ! ! Input, character ( len = * ) TITLE, a 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 ) if ( m <= 0 .or. n <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' (None)' return end if 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), '(i8,6x)' ) i end do write ( *, '('' Row '',5a14)' ) 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), '(g14.6)' ) a(i,j) end do write ( *, '(i5,a,5a14)' ) j, ':', ( ctemp(i), i = 1, inc ) end do end do return end subroutine r8col_undex ( m, n, a, unique_num, undx, xdnu ) !*****************************************************************************80 ! !! R8COL_UNDEX returns unique sorted indexes for an R8COL. ! ! Discussion: ! ! An R8COL is an M by N array of R8's, regarded as an array of N columns, ! each of length M. ! ! The goal of this routine is to determine a vector UNDX, ! which points, to the unique elements of A, in sorted order, ! and a vector XDNU, which identifies, for each entry of A, the index of ! the unique sorted element of A. ! ! This is all done with index vectors, so that the elements of ! A are never moved. ! ! The first step of the algorithm requires the indexed sorting ! of A, which creates arrays INDX and XDNI. (If all the entries ! of A are unique, then these arrays are the same as UNDX and XDNU.) ! ! We then use INDX to examine the entries of A in sorted order, ! noting the unique entries, creating the entries of XDNU and ! UNDX as we go. ! ! Once this process has been completed, the object X could be ! replaced by a compressed object XU, containing the unique entries ! of X in sorted order, using the formula ! ! XU(*) = A(UNDX(*)). ! ! We could then, if we wished, reconstruct the entire vector A, or ! any element of it, by index, as follows: ! ! A(I) = XU(XDNU(I)). ! ! We could then replace A by the combination of XU and XDNU. ! ! Later, when we need the I-th entry of A, we can locate it as ! the XDNU(I)-th entry of XU. ! ! Here is an example of a vector A, the sort and inverse sort ! index vectors, and the unique sort and inverse unique sort vectors ! and the compressed unique sorted vector. ! ! I A Indx Xdni XU Undx Xdnu ! ----+-----+-----+-----+--------+-----+-----+ ! 1 | 11. 1 1 | 11. 1 1 ! 2 | 22. 3 5 | 22. 2 2 ! 3 | 11. 6 2 | 33. 4 1 ! 4 | 33. 9 8 | 55. 5 3 ! 5 | 55. 2 9 | 4 ! 6 | 11. 7 3 | 1 ! 7 | 22. 8 6 | 2 ! 8 | 22. 4 7 | 2 ! 9 | 11. 5 4 | 1 ! ! INDX(2) = 3 means that sorted item(2) is A(3). ! XDNI(2) = 5 means that A(2) is sorted item(5). ! ! UNDX(3) = 4 means that unique sorted item(3) is at A(4). ! XDNU(8) = 2 means that A(8) is at unique sorted item(2). ! ! XU(XDNU(I))) = A(I). ! XU(I) = A(UNDX(I)). ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 17 July 2010 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, the dimension of the data values. ! ! Input, integer N, the number of data values. ! ! Input, real ( kind = rk ) A(M,N), the data values. ! ! Input, integer UNIQUE_NUM, the number of unique values ! in A. This value is only required for languages in which the size of ! UNDX must be known in advance. ! ! Output, integer UNDX(UNIQUE_NUM), the UNDX vector. ! ! Output, integer XDNU(N), the XDNU vector. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n integer unique_num real ( kind = rk ) a(m,n) real ( kind = rk ) diff integer i integer indx(n) integer j integer undx(unique_num) integer xdnu(n) ! ! Implicitly sort the array. ! call r8col_sort_heap_index_a ( m, n, a, indx ) ! ! Walk through the implicitly sorted array. ! i = 1 j = 1 undx(j) = indx(i) xdnu(indx(i)) = j do i = 2, n diff = maxval ( abs ( a(1:m,indx(i)) - a(1:m,undx(j)) ) ) if ( 0.0D+00 < diff ) then j = j + 1 undx(j) = indx(i) end if xdnu(indx(i)) = j end do return end subroutine r8col_uniform_ab ( m, n, a, b, r ) !*****************************************************************************80 ! !! R8COL_UNIFORM_AB returns a scaled pseudorandom R8COL. ! ! Discussion: ! ! An R8COL is an array of R8 values, regarded as a set of column vectors. ! ! A <= R(I,J) <= B. ! ! Licensing: ! ! This code is distributed under the MIT 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 M, N, the number of rows and columns ! in the array. ! ! Input, real ( kind = rk ) A, B, the lower and upper limits. ! ! Output, real ( kind = rk ) R(M,N), the array of pseudorandom values. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a real ( kind = rk ) b real ( kind = rk ) r(m,n) call random_number ( harvest = r(1:m,1:n) ) r(1:m,1:n) = a + ( b - a ) * r(1:m,1:n) return end subroutine r8col_uniform_abvec ( m, n, a, b, r ) !*****************************************************************************80 ! !! R8COL_UNIFORM_ABVEC fills an R8COL with scaled pseudorandom numbers. ! ! Discussion: ! ! An R8COL is an array of R8 values, regarded as a set of column vectors. ! ! The user specifies a minimum and maximum value for each row. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 19 December 2011 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Paul Bratley, Bennett Fox, Linus Schrage, ! A Guide to Simulation, ! Springer Verlag, pages 201-202, 1983. ! ! 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 M, N, the number of rows and columns in ! the array. ! ! Input, real ( kind = rk ) A(M), B(M), the lower and upper limits. ! ! Output, real ( kind = rk ) R(M,N), the array of pseudorandom values. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a(m) real ( kind = rk ) b(m) integer i real ( kind = rk ) r(m,n) call random_number ( harvest = r(1:m,1:n) ) do i = 1, m r(i,1:n) = a(i) + ( b(i) - a(i) ) * r(i,1:n) end do return end subroutine r8col_unique_count ( m, n, a, unique_num ) !*****************************************************************************80 ! !! R8COL_UNIQUE_COUNT counts the unique columns in an unsorted R8COL. ! ! Discussion: ! ! An R8COL is an M by N array of R8's, regarded as an array of N columns, ! each of length M. ! ! Because the array is unsorted, this algorithm is O(N^2). ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 17 July 2010 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, the number of rows. ! ! Input, integer N, the number of columns. ! ! Input, real ( kind = rk ) A(M,N), the array of N columns of data. ! ! Input, real ( kind = rk ) TOL, a nonnegative tolerance for equality. ! ! Output, integer UNIQUE_NUM, the number of unique columns. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a(m,n) real ( kind = rk ) diff integer j1 integer j2 logical ( kind = 4 ) unique(n) integer unique_num unique_num = 0 do j1 = 1, n unique_num = unique_num + 1 unique(j1) = .true. do j2 = 1, j1 - 1 if ( unique(j2) ) then diff = maxval ( abs ( a(1:m,j1) - a(1:m,j2) ) ) if ( diff == 0.0D+00 ) then unique_num = unique_num - 1 unique(j1) = .false. exit end if end if end do end do return end subroutine r8col_unique_index ( m, n, a, unique_index ) !*****************************************************************************80 ! !! R8COL_UNIQUE_INDEX indexes the unique occurrence of values in an R8COL. ! ! Discussion: ! ! An R8COL is an M by N array of R8's, regarded as an array of N columns, ! each of length M. ! ! For element A(1:M,J) of the matrix, UNIQUE_INDEX(J) is the uniqueness index ! of A(1:M,J). That is, if A_UNIQUE contains the unique elements of A, ! gathered in order, then ! ! A_UNIQUE ( 1:M, UNIQUE_INDEX(J) ) = A(1:M,J) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 17 July 2010 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns of A. ! The length of an "element" of A, and the number of "elements". ! ! Input, real ( kind = rk ) A(M,N), the array. ! ! Output, integer UNIQUE_INDEX(N), the first occurrence index. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a(m,n) real ( kind = rk ) diff integer j1 integer j2 integer unique_index(n) integer unique_num unique_index(1:n) = -1 unique_num = 0 do j1 = 1, n if ( unique_index(j1) == -1 ) then unique_num = unique_num + 1 unique_index(j1) = unique_num do j2 = j1 + 1, n diff = maxval ( abs ( a(1:m,j1) - a(1:m,j2) ) ) if ( diff == 0.0D+00 ) then unique_index(j2) = unique_num end if end do end if end do return end subroutine r8col_variance ( m, n, a, variance ) !*****************************************************************************80 ! !! R8COL_VARIANCE returns the variances of an R8COL. ! ! Discussion: ! ! An R8COL is an M by N array of R8's, regarded as an array of N columns, ! each of length M. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 05 December 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns in ! the array. ! ! Input, real ( kind = rk ) A(M,N), the array whose variances are desired. ! ! Output, real ( kind = rk ) VARIANCE(N), the variances of the rows. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a(m,n) integer i integer j real ( kind = rk ) mean real ( kind = rk ) variance(n) do j = 1, n mean = sum ( a(1:m,j) ) / real ( m, kind = rk ) variance(j) = 0.0D+00 do i = 1, m variance(j) = variance(j) + ( a(i,j) - mean )**2 end do if ( 1 < m ) then variance(j) = variance(j) / real ( m - 1, kind = rk ) else variance(j) = 0.0D+00 end if end do return end subroutine r8vec_compare ( n, a1, a2, isgn ) !*****************************************************************************80 ! !! R8VEC_COMPARE compares two R8VEC's. ! ! Discussion: ! ! An R8VEC is a vector of R8's. ! ! The lexicographic ordering is used. ! ! Example: ! ! Input: ! ! A1 = ( 2.0, 6.0, 2.0 ) ! A2 = ( 2.0, 8.0, 12.0 ) ! ! Output: ! ! ISGN = -1 ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 23 February 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of entries in the vectors. ! ! Input, real ( kind = rk ) A1(N), A2(N), the vectors to be compared. ! ! Output, integer ISGN, the results of the comparison: ! -1, A1 < A2, ! 0, A1 = A2, ! +1, A1 > A2. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ) a1(n) real ( kind = rk ) a2(n) integer isgn integer k isgn = 0 k = 1 do while ( k <= n ) if ( a1(k) < a2(k) ) then isgn = -1 return else if ( a2(k) < a1(k) ) then isgn = + 1 return end if k = k + 1 end do return end function r8vec_eq ( n, a1, a2 ) !*****************************************************************************80 ! !! R8VEC_EQ is true if two R8VECs are equal. ! ! Discussion: ! ! An R8VEC is a vector of R8's. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 05 December 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of entries in the vectors. ! ! Input, real ( kind = rk ) A1(N), A2(N), two vectors to compare. ! ! Output, logical ( kind = 4 ) R8VEC_EQ, is TRUE if every pair of elements ! A1(I) and A2(I) are equal, and FALSE otherwise. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ) a1(n) real ( kind = rk ) a2(n) logical ( kind = 4 ) r8vec_eq r8vec_eq = ( all ( a1(1:n) == a2(1:n) ) ) return end function r8vec_gt ( n, a1, a2 ) !*****************************************************************************80 ! !! R8VEC_GT == ( A1 > A2 ) for R8VEC's. ! ! Discussion: ! ! An R8VEC is a vector of R8's. ! ! The comparison is lexicographic. ! ! A1 > A2 <=> A1(1) > A2(1) or ! ( A1(1) == A2(1) and A1(2) > A2(2) ) or ! ... ! ( A1(1:N-1) == A2(1:N-1) and A1(N) > A2(N) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 05 December 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the dimension of the vectors. ! ! Input, real ( kind = rk ) A1(N), A2(N), the vectors to be compared. ! ! Output, logical ( kind = 4 ) R8VEC_GT, is TRUE if and only if A1 > A2. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ) a1(n) real ( kind = rk ) a2(n) integer i logical ( kind = 4 ) r8vec_gt r8vec_gt = .false. do i = 1, n if ( a2(i) < a1(i) ) then r8vec_gt = .true. exit else if ( a1(i) < a2(i) ) then r8vec_gt = .false. exit end if end do return end function r8vec_lt ( n, a1, a2 ) !*****************************************************************************80 ! !! R8VEC_LT == ( A1 < A2 ) for R8VEC's. ! ! Discussion: ! ! An R8VEC is a vector of R8's. ! ! The comparison is lexicographic. ! ! A1 < A2 <=> A1(1) < A2(1) or ! ( A1(1) == A2(1) and A1(2) < A2(2) ) or ! ... ! ( A1(1:N-1) == A2(1:N-1) and A1(N) < A2(N) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 05 December 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the dimension of the vectors. ! ! Input, real ( kind = rk ) A1(N), A2(N), the vectors to be compared. ! ! Output, logical ( kind = 4 ) R8VEC_LT, is TRUE if and only if A1 < A2. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ) a1(n) real ( kind = rk ) a2(n) logical ( kind = 4 ) r8vec_lt integer i r8vec_lt = .false. do i = 1, n if ( a1(i) < a2(i) ) then r8vec_lt = .true. exit else if ( a2(i) < a1(i) ) then r8vec_lt = .false. exit end if end do return end subroutine r8vec_print ( n, a, title ) !*****************************************************************************80 ! !! R8VEC_PRINT prints an R8VEC. ! ! Discussion: ! ! An R8VEC is a vector of R8's. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 22 August 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of components of the vector. ! ! Input, real ( kind = rk ) A(N), the vector to be printed. ! ! Input, character ( len = * ) TITLE, a title. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ) a(n) integer i character ( len = * ) title write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) write ( *, '(a)' ) ' ' do i = 1, n write ( *, '(2x,i8,a,1x,g16.8)' ) i, ':', a(i) end do return end subroutine r8vec_swap ( n, a1, a2 ) !*****************************************************************************80 ! !! R8VEC_SWAP swaps the entries of two R8VECs. ! ! Discussion: ! ! An R8VEC is a vector of R8's. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 04 December 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of entries in the arrays. ! ! Input/output, real ( kind = rk ) A1(N), A2(N), the vectors to swap. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ) a1(n) real ( kind = rk ) a2(n) real ( kind = rk ) a3(n) a3(1:n) = a1(1:n) a1(1:n) = a2(1:n) a2(1:n) = a3(1:n) 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: ! ! 05 February 2004 ! ! Author: ! ! Original FORTRAN77 version by Albert Nijenhuis, Herbert Wilf. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Albert Nijenhuis, Herbert Wilf, ! Combinatorial Algorithms for Computers and Calculators, ! Academic Press, 1978, ! ISBN: 0-12-519260-6, ! LC: QA164.N54. ! ! Parameters: ! ! Input, integer N, the number of items to be sorted. ! ! Input/output, integer INDX, the main communication signal. ! The user must set INDX to 0 before the first call. ! Thereafter, the user should not change the value of INDX until ! the sorting is done. ! 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. ! ! Output, 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. ! ! Input, 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. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) 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: ! ! 31 May 2001 9:45:54.872 AM ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 18 May 2013 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! None ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) character ( len = 8 ) ampm integer d integer h integer m integer mm character ( len = 9 ), parameter, dimension(12) :: month = (/ & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' /) integer n integer s integer values(8) integer y call date_and_time ( values = values ) y = values(1) m = values(2) d = values(3) h = values(5) n = values(6) s = values(7) mm = values(8) if ( h < 12 ) then ampm = 'AM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Noon' else ampm = 'PM' end if else h = h - 12 if ( h < 12 ) then ampm = 'PM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Midnight' else ampm = 'AM' end if end if end if write ( *, '(i2.2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm ) return end