program main !*****************************************************************************80 ! !! jacobi_rule() computes a standard Gauss-Jacobi quadrature rule. ! ! Discussion: ! ! This program computes a standard Gauss-Jacobi quadrature rule ! and writes it to a file. ! ! The user specifies: ! * the ORDER (number of points) in the rule; ! * ALPHA, the exponent of (1-x); ! * BETA, the exponent of (1+x); ! * A, the left endpoint; ! * B, the right endpoint; ! * FILENAME, fhe root name of the output files. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 28 February 2010 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) a real ( kind = rk ) alpha integer arg_num real ( kind = rk ) b real ( kind = rk ) beta character ( len = 255 ) filename integer iarg integer iargc integer ierror integer kinda integer last integer order real ( kind = rk ) r(2) character ( len = 255 ) string real ( kind = rk ), allocatable, dimension ( : ) :: w real ( kind = rk ), allocatable, dimension ( : ) :: x call timestamp ( ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'JACOBI_RULE' write ( *, '(a)' ) ' FORTRAN90 version' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Compute a Gauss-Jacobi rule for approximating' write ( *, '(a)' ) & ' Integral ( A <= x <= B ) (B-x)^alpha (x-A)^beta f(x) dx' write ( *, '(a)' ) ' of order ORDER.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) & ' The user specifies ORDER, ALPHA, BETA, A, B, and FILENAME.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' ORDER is the number of points:' write ( *, '(a)' ) ' ALPHA is the exponent of (B-X)' write ( *, '(a)' ) ' BETA is the exponent of (X-A):' write ( *, '(a)' ) ' A is the left endpoint;' write ( *, '(a)' ) ' B is the right endpoint;' write ( *, '(a)' ) ' FILENAME is used to generate 3 files:' write ( *, '(a)' ) ' filename_w.txt - the weight file' write ( *, '(a)' ) ' filename_x.txt - the abscissa file.' write ( *, '(a)' ) ' filename_r.txt - the region file.' ! ! Get the number of command line arguments. ! arg_num = iargc ( ) ! ! Get ORDER. ! if ( 1 <= arg_num ) then iarg = 1 call getarg ( iarg, string ) call s_to_i4 ( string, order, ierror, last ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Enter the rule order ORDER:' read ( *, * ) order end if ! ! Get ALPHA ! if ( 2 <= arg_num ) then iarg = 2 call getarg ( iarg, string ) call s_to_r8 ( string, alpha, ierror, last ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' ALPHA is the exponent of (B-x) in the integral:' write ( *, '(a)' ) ' Note that -1.0 < ALPHA is required.' write ( *, '(a)' ) ' Enter the value of ALPHA:' read ( *, * ) alpha end if ! ! Get BETA ! if ( 3 <= arg_num ) then iarg = 3 call getarg ( iarg, string ) call s_to_r8 ( string, beta, ierror, last ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' BETA is the exponent of (x-A) in the integral:' write ( *, '(a)' ) ' Note that -1.0 < BETA is required.' write ( *, '(a)' ) ' Enter the value of BETA:' read ( *, * ) beta end if ! ! Get A. ! if ( 4 <= arg_num ) then iarg = 4 call getarg ( iarg, string ) call s_to_r8 ( string, a, ierror, last ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Enter A, the left endpoint:' read ( *, * ) a end if ! ! Get B. ! if ( 5 <= arg_num ) then iarg = 5 call getarg ( iarg, string ) call s_to_r8 ( string, b, ierror, last ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Enter B, the right endpoint:' read ( *, * ) b end if ! ! Get FILENAME. ! if ( 6 <= arg_num ) then iarg = 6 call getarg ( iarg, filename ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) & ' Enter FILENAME, the "root name" of the quadrature files).' read ( *, '(a)' ) filename end if ! ! Input summary. ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Input summary:' write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) ' ORDER = ', order write ( *, '(a,g14.6)' ) ' ALPHA = ', alpha write ( *, '(a,g14.6)' ) ' BETA = ', beta write ( *, '(a,g14.6)' ) ' A = ', a write ( *, '(a,g14.6)' ) ' B = ', b write ( *, '(a)' ) ' OUTPUT = "' // trim ( filename ) // '".' ! ! Construct the rule. ! allocate ( w(order) ) allocate ( x(order) ) kinda = 4 call cgqf ( order, kinda, alpha, beta, a, b, x, w ) ! ! Write the rule. ! r(1) = a r(2) = b call rule_write ( order, x, w, r, filename ) ! ! Free memory. ! deallocate ( w ) deallocate ( x ) ! ! Terminate. ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'JACOBI_RULE:' write ( *, '(a)' ) ' Normal end of execution.' write ( *, '(a)' ) ' ' call timestamp ( ) stop end subroutine cdgqf ( nt, kinda, alpha, beta, t, wts ) !*****************************************************************************80 ! !! CDGQF computes a Gauss quadrature formula with default A, B and simple knots. ! ! Discussion: ! ! This routine computes all the knots and weights of a Gauss quadrature ! formula with a classical weight function with default values for A and B, ! and only simple knots. ! ! There are no moments checks and no printing is done. ! ! Use routine EIQFS to evaluate a quadrature computed by CGQFS. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 04 January 2010 ! ! Author: ! ! Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Sylvan Elhay, Jaroslav Kautsky, ! Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of ! Interpolatory Quadrature, ! ACM Transactions on Mathematical Software, ! Volume 13, Number 4, December 1987, pages 399-415. ! ! Parameters: ! ! Input, integer NT, the number of knots. ! ! Input, integer KINDA, the rule. ! 1, Legendre, (a,b) 1.0 ! 2, Chebyshev, (a,b) ((b-x)*(x-a))^(-0.5) ! 3, Gegenbauer, (a,b) ((b-x)*(x-a))^alpha ! 4, Jacobi, (a,b) (b-x)^alpha*(x-a)^beta ! 5, Generalized Laguerre, (a,inf) (x-a)^alpha*exp(-b*(x-a)) ! 6, Generalized Hermite, (-inf,inf) |x-a|^alpha*exp(-b*(x-a)^2) ! 7, Exponential, (a,b) |x-(a+b)/2.0|^alpha ! 8, Rational, (a,inf) (x-a)^alpha*(x+b)^beta ! ! Input, real ( kind = rk ) ALPHA, the value of Alpha, if needed. ! ! Input, real ( kind = rk ) BETA, the value of Beta, if needed. ! ! Output, real ( kind = rk ) T(NT), the knots. ! ! Output, real ( kind = rk ) WTS(NT), the weights. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer nt real ( kind = rk ) aj(nt) real ( kind = rk ) alpha real ( kind = rk ) beta real ( kind = rk ) bj(nt) integer kinda real ( kind = rk ) t(nt) real ( kind = rk ) wts(nt) real ( kind = rk ) zemu call parchk ( kinda, 2 * nt, alpha, beta ) ! ! Get the Jacobi matrix and zero-th moment. ! call class_matrix ( kinda, nt, alpha, beta, aj, bj, zemu ) ! ! Compute the knots and weights. ! call sgqf ( nt, aj, bj, zemu, t, wts ) return end subroutine cgqf ( nt, kinda, alpha, beta, a, b, t, wts ) !*****************************************************************************80 ! !! CGQF computes knots and weights of a Gauss quadrature formula. ! ! Discussion: ! ! The user may specify the interval (A,B). ! ! Only simple knots are produced. ! ! Use routine EIQFS to evaluate this quadrature formula. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 February 2010 ! ! Author: ! ! Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Sylvan Elhay, Jaroslav Kautsky, ! Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of ! Interpolatory Quadrature, ! ACM Transactions on Mathematical Software, ! Volume 13, Number 4, December 1987, pages 399-415. ! ! Parameters: ! ! Input, integer NT, the number of knots. ! ! Input, integer KINDA, the rule. ! 1, Legendre, (a,b) 1.0 ! 2, Chebyshev Type 1, (a,b) ((b-x)*(x-a))^-0.5) ! 3, Gegenbauer, (a,b) ((b-x)*(x-a))^alpha ! 4, Jacobi, (a,b) (b-x)^alpha*(x-a)^beta ! 5, Generalized Laguerre, (a,+oo) (x-a)^alpha*exp(-b*(x-a)) ! 6, Generalized Hermite, (-oo,+oo) |x-a|^alpha*exp(-b*(x-a)^2) ! 7, Exponential, (a,b) |x-(a+b)/2.0|^alpha ! 8, Rational, (a,+oo) (x-a)^alpha*(x+b)^beta ! 9, Chebyshev Type 2, (a,b) ((b-x)*(x-a))^(+0.5) ! ! Input, real ( kind = rk ) ALPHA, the value of Alpha, if needed. ! ! Input, real ( kind = rk ) BETA, the value of Beta, if needed. ! ! Input, real ( kind = rk ) A, B, the interval endpoints, or ! other parameters. ! ! Output, real ( kind = rk ) T(NT), the knots. ! ! Output, real ( kind = rk ) WTS(NT), the weights. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer nt real ( kind = rk ) a real ( kind = rk ) alpha real ( kind = rk ) b real ( kind = rk ) beta integer i integer kinda integer, allocatable :: mlt(:) integer, allocatable :: ndx(:) real ( kind = rk ) t(nt) real ( kind = rk ) wts(nt) ! ! Compute the Gauss quadrature formula for default values of A and B. ! call cdgqf ( nt, kinda, alpha, beta, t, wts ) ! ! Prepare to scale the quadrature formula to other weight function with ! valid A and B. ! allocate ( mlt(1:nt) ) mlt(1:nt) = 1 allocate ( ndx(1:nt) ) do i = 1, nt ndx(i) = i end do call scqf ( nt, t, mlt, wts, nt, ndx, wts, t, kinda, alpha, beta, a, b ) deallocate ( mlt ) deallocate ( ndx ) return end subroutine ch_cap ( c ) !*****************************************************************************80 ! !! CH_CAP capitalizes a single character. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 19 July 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, character C, the character to capitalize. ! implicit none character c integer itemp itemp = ichar ( c ) if ( 97 <= itemp .and. itemp <= 122 ) then c = 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: ! ! 28 July 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character C1, C2, the characters to compare. ! ! Output, logical CH_EQI, the result of the comparison. ! implicit none logical ch_eqi character c1 character c1_cap character c2 character c2_cap 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: ! ! 04 August 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! 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 class_matrix ( kinda, m, alpha, beta, aj, bj, zemu ) !*****************************************************************************80 ! !! CLASS_MATRIX computes the Jacobi matrix for a quadrature rule. ! ! Discussion: ! ! This routine computes the diagonal AJ and sub-diagonal BJ ! elements of the order M tridiagonal symmetric Jacobi matrix ! associated with the polynomials orthogonal with respect to ! the weight function specified by KIND. ! ! For weight functions 1-7, M elements are defined in BJ even ! though only M-1 are needed. For weight function 8, BJ(M) is ! set to zero. ! ! The zero-th moment of the weight function is returned in ZEMU. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 27 December 2009 ! ! Author: ! ! Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Sylvan Elhay, Jaroslav Kautsky, ! Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of ! Interpolatory Quadrature, ! ACM Transactions on Mathematical Software, ! Volume 13, Number 4, December 1987, pages 399-415. ! ! Parameters: ! ! Input, integer KINDA, the rule. ! 1, Legendre, (a,b) 1.0 ! 2, Chebyshev, (a,b) ((b-x)*(x-a))^(-0.5) ! 3, Gegenbauer, (a,b) ((b-x)*(x-a))^alpha ! 4, Jacobi, (a,b) (b-x)^alpha*(x-a)^beta ! 5, Generalized Laguerre, (a,inf) (x-a)^alpha*exp(-b*(x-a)) ! 6, Generalized Hermite, (-inf,inf) |x-a|^alpha*exp(-b*(x-a)^2) ! 7, Exponential, (a,b) |x-(a+b)/2.0|^alpha ! 8, Rational, (a,inf) (x-a)^alpha*(x+b)^beta ! ! Input, integer M, the order of the Jacobi matrix. ! ! Input, real ( kind = rk ) ALPHA, the value of Alpha, if needed. ! ! Input, real ( kind = rk ) BETA, the value of Beta, if needed. ! ! Output, real ( kind = rk ) AJ(M), BJ(M), the diagonal and subdiagonal ! of the Jacobi matrix. ! ! Output, real ( kind = rk ) ZEMU, the zero-th moment. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m real ( kind = rk ) a2b2 real ( kind = rk ) ab real ( kind = rk ) aba real ( kind = rk ) abi real ( kind = rk ) abj real ( kind = rk ) abti real ( kind = rk ) aj(m) real ( kind = rk ) alpha real ( kind = rk ) apone real ( kind = rk ) beta real ( kind = rk ) bj(m) integer i integer kinda real ( kind = rk ), parameter :: pi = 3.14159265358979323846264338327950D+00 real ( kind = rk ) temp real ( kind = rk ) temp2 real ( kind = rk ) zemu temp = epsilon ( temp ) call parchk ( kinda, 2 * m - 1, alpha, beta ) temp2 = 0.5D+00 if ( 500.0D+00 * temp < abs ( ( gamma ( temp2 ) )**2 - pi ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CLASS_MATRIX - Fatal error!' write ( *, '(a)' ) ' Gamma function does not match machine parameters.' stop end if if ( kinda == 1 ) then ab = 0.0D+00 zemu = 2.0D+00 / ( ab + 1.0D+00 ) aj(1:m) = 0.0D+00 do i = 1, m abi = i + ab * mod ( i, 2 ) abj = 2 * i + ab bj(i) = abi * abi / ( abj * abj - 1.0D+00 ) end do bj(1:m) = sqrt ( bj(1:m) ) else if ( kinda == 2 ) then zemu = pi aj(1:m) = 0.0D+00 bj(1) = sqrt ( 0.5D+00 ) bj(2:m) = 0.5D+00 else if ( kinda == 3 ) then ab = alpha * 2.0D+00 zemu = 2.0D+00**( ab + 1.0D+00 ) * gamma ( alpha + 1.0D+00 )**2 & / gamma ( ab + 2.0D+00 ) aj(1:m) = 0.0D+00 bj(1) = 1.0D+00 / ( 2.0D+00 * alpha + 3.0D+00 ) do i = 2, m bj(i) = i * ( i + ab ) / ( 4.0D+00 * ( i + alpha )**2 - 1.0D+00 ) end do bj(1:m) = sqrt ( bj(1:m) ) else if ( kinda == 4 ) then ab = alpha + beta abi = 2.0D+00 + ab zemu = 2.0D+00**( ab + 1.0D+00 ) * gamma ( alpha + 1.0D+00 ) & * gamma ( beta + 1.0D+00 ) / gamma ( abi ) aj(1) = ( beta - alpha ) / abi bj(1) = 4.0D+00 * ( 1.0 + alpha ) * ( 1.0D+00 + beta ) & / ( ( abi + 1.0D+00 ) * abi * abi ) a2b2 = beta * beta - alpha * alpha do i = 2, m abi = 2.0D+00 * i + ab aj(i) = a2b2 / ( ( abi - 2.0D+00 ) * abi ) abi = abi**2 bj(i) = 4.0D+00 * i * ( i + alpha ) * ( i + beta ) * ( i + ab ) & / ( ( abi - 1.0D+00 ) * abi ) end do bj(1:m) = sqrt ( bj(1:m) ) else if ( kinda == 5 ) then zemu = gamma ( alpha + 1.0D+00 ) do i = 1, m aj(i) = 2.0D+00 * i - 1.0D+00 + alpha bj(i) = i * ( i + alpha ) end do bj(1:m) = sqrt ( bj(1:m) ) else if ( kinda == 6 ) then zemu = gamma ( ( alpha + 1.0D+00 ) / 2.0D+00 ) aj(1:m) = 0.0D+00 do i = 1, m bj(i) = ( i + alpha * mod ( i, 2 ) ) / 2.0D+00 end do bj(1:m) = sqrt ( bj(1:m) ) else if ( kinda == 7 ) then ab = alpha zemu = 2.0D+00 / ( ab + 1.0D+00 ) aj(1:m) = 0.0D+00 do i = 1, m abi = i + ab * mod ( i, 2 ) abj = 2 * i + ab bj(i) = abi * abi / ( abj * abj - 1.0D+00 ) end do bj(1:m) = sqrt ( bj(1:m) ) else if ( kinda == 8 ) then ab = alpha + beta zemu = gamma ( alpha + 1.0D+00 ) * gamma ( - ( ab + 1.0D+00 ) ) & / gamma ( - beta ) apone = alpha + 1.0D+00 aba = ab * apone aj(1) = - apone / ( ab + 2.0D+00 ) bj(1) = - aj(1) * ( beta + 1.0D+00 ) / ( ab + 2.0D+00 ) / ( ab + 3.0D+00 ) do i = 2, m abti = ab + 2.0D+00 * i aj(i) = aba + 2.0D+00 * ( ab + i ) * ( i - 1 ) aj(i) = - aj(i) / abti / ( abti - 2.0D+00 ) end do do i = 2, m - 1 abti = ab + 2.0D+00 * i bj(i) = i * ( alpha + i ) / ( abti - 1.0D+00 ) * ( beta + i ) & / ( abti**2 ) * ( ab + i ) / ( abti + 1.0D+00 ) end do bj(m) = 0.0D+00 bj(1:m) = sqrt ( bj(1:m) ) end if return end subroutine get_unit ( iunit ) !*****************************************************************************80 ! !! GET_UNIT returns a free FORTRAN unit number. ! ! Discussion: ! ! A "free" FORTRAN unit number is an integer between 1 and 99 which ! is not currently associated with an I/O device. A free FORTRAN unit ! number is needed in order to open a file with the OPEN command. ! ! If IUNIT = 0, then no free FORTRAN unit could be found, although ! all 99 units were checked (except for units 5, 6 and 9, which ! are commonly reserved for console I/O). ! ! Otherwise, IUNIT is an integer between 1 and 99, representing a ! free FORTRAN unit. Note that GET_UNIT assumes that units 5 and 6 ! are special, and will never return those values. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 18 September 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! 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 imtqlx ( n, d, e, z ) !*****************************************************************************80 ! !! IMTQLX diagonalizes a symmetric tridiagonal matrix. ! ! Discussion: ! ! This routine is a slightly modified version of the EISPACK routine to ! perform the implicit QL algorithm on a symmetric tridiagonal matrix. ! ! The authors thank the authors of EISPACK for permission to use this ! routine. ! ! It has been modified to produce the product Q' * Z, where Z is an input ! vector and Q is the orthogonal matrix diagonalizing the input matrix. ! The changes consist (essentially) of applying the orthogonal ! transformations directly to Z as they are generated. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 27 December 2009 ! ! Author: ! ! Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Sylvan Elhay, Jaroslav Kautsky, ! Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of ! Interpolatory Quadrature, ! ACM Transactions on Mathematical Software, ! Volume 13, Number 4, December 1987, pages 399-415. ! ! Roger Martin, James Wilkinson, ! The Implicit QL Algorithm, ! Numerische Mathematik, ! Volume 12, Number 5, December 1968, pages 377-383. ! ! Parameters: ! ! Input, integer N, the order of the matrix. ! ! Input/output, real ( kind = rk ) D(N), the diagonal entries of the matrix. ! On output, the information in D has been overwritten. ! ! Input/output, real ( kind = rk ) E(N), the subdiagonal entries of the ! matrix, in entries E(1) through E(N-1). On output, the information in ! E has been overwritten. ! ! Input/output, real ( kind = rk ) Z(N). On input, a vector. On output, ! the value of Q' * Z, where Q is the matrix that diagonalizes the ! input symmetric tridiagonal matrix. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ) b real ( kind = rk ) c real ( kind = rk ) d(n) real ( kind = rk ) e(n) real ( kind = rk ) f real ( kind = rk ) g integer i integer ii integer, parameter :: itn = 30 integer j integer k integer l integer m integer mml real ( kind = rk ) p real ( kind = rk ) prec real ( kind = rk ) r real ( kind = rk ) s real ( kind = rk ) z(n) prec = epsilon ( prec ) if ( n == 1 ) then return end if e(n) = 0.0D+00 do l = 1, n j = 0 do do m = l, n if ( m == n ) then exit end if if ( abs ( e(m) ) <= prec * ( abs ( d(m) ) + abs ( d(m+1) ) ) ) then exit end if end do p = d(l) if ( m == l ) then exit end if if ( itn <= j ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'IMTQLX - Fatal error!' write ( *, '(a)' ) ' Iteration limit exceeded.' write ( *, '(a,i8)' ) ' J = ', j write ( *, '(a,i8)' ) ' L = ', l write ( *, '(a,i8)' ) ' M = ', m write ( *, '(a,i8)' ) ' N = ', n stop end if j = j + 1 g = ( d(l+1) - p ) / ( 2.0D+00 * e(l) ) r = sqrt ( g * g + 1.0D+00 ) g = d(m) - p + e(l) / ( g + sign ( r, g ) ) s = 1.0D+00 c = 1.0D+00 p = 0.0D+00 mml = m - l do ii = 1, mml i = m - ii f = s * e(i) b = c * e(i) if ( abs ( g ) <= abs ( f ) ) then c = g / f r = sqrt ( c * c + 1.0D+00 ) e(i+1) = f * r s = 1.0D+00 / r c = c * s else s = f / g r = sqrt ( s * s + 1.0D+00 ) e(i+1) = g * r c = 1.0D+00 / r s = s * c end if g = d(i+1) - p r = ( d(i) - g ) * s + 2.0D+00 * c * b p = s * r d(i+1) = g + p g = c * r - b f = z(i+1) z(i+1) = s * z(i) + c * f z(i) = c * z(i) - s * f end do d(l) = d(l) - p e(l) = g e(m) = 0.0D+00 end do end do ! ! Sorting. ! do ii = 2, n i = ii - 1 k = i p = d(i) do j = ii, n if ( d(j) < p ) then k = j p = d(j) end if end do if ( k /= i ) then d(k) = d(i) d(i) = p p = z(i) z(i) = z(k) z(k) = p end if end do return end subroutine parchk ( kinda, m, alpha, beta ) !*****************************************************************************80 ! !! PARCHK checks parameters ALPHA and BETA for classical weight functions. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 27 December 2009 ! ! Author: ! ! Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Sylvan Elhay, Jaroslav Kautsky, ! Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of ! Interpolatory Quadrature, ! ACM Transactions on Mathematical Software, ! Volume 13, Number 4, December 1987, pages 399-415. ! ! Parameters: ! ! Input, integer KINDA, the rule. ! 1, Legendre, (a,b) 1.0 ! 2, Chebyshev, (a,b) ((b-x)*(x-a))^(-0.5) ! 3, Gegenbauer, (a,b) ((b-x)*(x-a))^alpha ! 4, Jacobi, (a,b) (b-x)^alpha*(x-a)^beta ! 5, Generalized Laguerre, (a,inf) (x-a)^alpha*exp(-b*(x-a)) ! 6, Generalized Hermite, (-inf,inf) |x-a|^alpha*exp(-b*(x-a)^2) ! 7, Exponential, (a,b) |x-(a+b)/2.0|^alpha ! 8, Rational, (a,inf) (x-a)^alpha*(x+b)^beta ! ! Input, integer M, the order of the highest moment to ! be calculated. This value is only needed when kind = rk. ! ! Input, real ( kind = rk ) ALPHA, BETA, the parameters, if required ! by the value of KIND. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) alpha real ( kind = rk ) beta integer kinda integer m real ( kind = rk ) tmp if ( kinda <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'PARCHK - Fatal error!' write ( *, '(a)' ) ' KINDA <= 0.' stop end if ! ! Check ALPHA for Gegenbauer, Jacobi, Laguerre, Hermite, Exponential. ! if ( 3 <= kinda .and. alpha <= -1.0D+00 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'PARCHK - Fatal error!' write ( *, '(a)' ) ' 3 <= KINDA and ALPHA <= -1.' stop end if ! ! Check BETA for Jacobi. ! if ( kinda == 4 .and. beta <= -1.0D+00 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'PARCHK - Fatal error!' write ( *, '(a)' ) ' KINDA == 4 and BETA <= -1.0.' stop end if ! ! Check ALPHA and BETA for rational. ! if ( kinda == 8 ) then tmp = alpha + beta + m + 1.0D+00 if ( 0.0D+00 <= tmp .or. tmp <= beta ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'PARCHK - Fatal error!' write ( *, '(a)' ) ' KINDA == 8 but condition on ALPHA and BETA fails.' stop end if end if return end function r8_epsilon ( ) !*****************************************************************************80 ! !! R8_EPSILON returns the R8 roundoff unit. ! ! Discussion: ! ! The roundoff unit is a number R which is a power of 2 with the ! property that, to the precision of the computer's arithmetic, ! 1 < 1 + R ! but ! 1 = ( 1 + R / 2 ) ! ! FORTRAN90 provides the superior library routine ! ! EPSILON ( X ) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 01 September 2012 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = rk ) R8_EPSILON, the round-off unit. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) r8_epsilon r8_epsilon = 2.220446049250313D-016 return end subroutine r8mat_write ( output_filename, m, n, table ) !*****************************************************************************80 ! !! R8MAT_WRITE writes an R8MAT file. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 31 May 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) OUTPUT_FILENAME, the output file name. ! ! Input, integer M, the spatial dimension. ! ! Input, integer N, the number of points. ! ! Input, real ( kind = rk ) TABLE(M,N), the table 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 rule_write ( order, x, w, r, filename ) !*****************************************************************************80 ! !! RULE_WRITE writes a quadrature rule to a file. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 18 February 2010 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ORDER, the order of the rule. ! ! Input, real ( kind = rk ) X(ORDER), the abscissas. ! ! Input, real ( kind = rk ) W(ORDER), the weights. ! ! Input, real ( kind = rk ) R(2), defines the region. ! ! Input, character ( len = * ) FILENAME, specifies the output. ! 'filename_w.txt', 'filename_x.txt', 'filename_r.txt' defining weights, ! abscissas, and region. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer order character ( len = * ) filename character ( len = 255 ) filename_r character ( len = 255 ) filename_w character ( len = 255 ) filename_x real ( kind = rk ) r(2) real ( kind = rk ) w(order) real ( kind = rk ) x(order) filename_w = trim ( filename ) // '_w.txt' filename_x = trim ( filename ) // '_x.txt' filename_r = trim ( filename ) // '_r.txt' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Creating quadrature files.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' "Root" file name is "' // trim ( filename ) // '".' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Weight file will be "' // trim ( filename_w ) // '".' write ( *, '(a)' ) ' Abscissa file will be "' // trim ( filename_x ) // '".' write ( *, '(a)' ) ' Region file will be "' // trim ( filename_r ) // '".' call r8mat_write ( filename_w, 1, order, w ) call r8mat_write ( filename_x, 1, order, x ) call r8mat_write ( filename_r, 1, 2, r ) 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: ! ! 28 June 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! 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. ! ! Output, integer IERROR, an error flag. ! 0, no error. ! 1, an error occurred. ! ! Output, 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_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: ! ! 07 September 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! 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. ! ! Output, 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. ! ! Output, 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 ) character c logical ch_eqi 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 scqf ( nt, t, mlt, wts, nwts, ndx, swts, st, kinda, alpha, beta, a, & b ) !*****************************************************************************80 ! !! SCQF scales a quadrature formula to a nonstandard interval. ! ! Discussion: ! ! The arrays WTS and SWTS may coincide. ! ! The arrays T and ST may coincide. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 27 December 2009 ! ! Author: ! ! Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Sylvan Elhay, Jaroslav Kautsky, ! Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of ! Interpolatory Quadrature, ! ACM Transactions on Mathematical Software, ! Volume 13, Number 4, December 1987, pages 399-415. ! ! Parameters: ! ! Input, integer NT, the number of knots. ! ! Input, real ( kind = rk ) T(NT), the original knots. ! ! Input, integer MLT(NT), the multiplicity of the knots. ! ! Input, real ( kind = rk ) WTS(NWTS), the weights. ! ! Input, integer NWTS, the number of weights. ! ! Input, integer NDX(NT), used to index the array WTS. ! For more details see the comments in CAWIQ. ! ! Output, real ( kind = rk ) SWTS(NWTS), the scaled weights. ! ! Output, real ( kind = rk ) ST(NT), the scaled knots. ! ! Input, integer KINDA, the rule. ! 1, Legendre, (a,b) 1.0 ! 2, Chebyshev Type 1, (a,b) ((b-x)*(x-a))^(-0.5) ! 3, Gegenbauer, (a,b) ((b-x)*(x-a))^alpha ! 4, Jacobi, (a,b) (b-x)^alpha*(x-a)^beta ! 5, Generalized Laguerre, (a,+oo) (x-a)^alpha*exp(-b*(x-a)) ! 6, Generalized Hermite, (-oo,+oo) |x-a|^alpha*exp(-b*(x-a)^2) ! 7, Exponential, (a,b) |x-(a+b)/2.0|^alpha ! 8, Rational, (a,+oo) (x-a)^alpha*(x+b)^beta ! 9, Chebyshev Type 2, (a,b) ((b-x)*(x-a))^(+0.5) ! ! Input, real ( kind = rk ) ALPHA, the value of Alpha, if needed. ! ! Input, real ( kind = rk ) BETA, the value of Beta, if needed. ! ! Input, real ( kind = rk ) A, B, the interval endpoints. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer nt integer nwts real ( kind = rk ) a real ( kind = rk ) al real ( kind = rk ) alpha real ( kind = rk ) b real ( kind = rk ) be real ( kind = rk ) beta integer i integer k integer kinda integer l integer mlt(nt) integer ndx(nt) real ( kind = rk ) p real ( kind = rk ) shft real ( kind = rk ) slp real ( kind = rk ) st(nt) real ( kind = rk ) swts(nwts) real ( kind = rk ) t(nt) real ( kind = rk ) temp real ( kind = rk ) tmp real ( kind = rk ) wts(nwts) temp = epsilon ( temp ) call parchk ( kinda, 1, alpha, beta ) if ( kinda == 1 ) then al = 0.0D+00 be = 0.0D+00 if ( abs ( b - a ) <= temp ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SCQF - Fatal error!' write ( *, '(a)' ) ' |B - A| too small.' stop end if shft = ( a + b ) / 2.0D+00 slp = ( b - a ) / 2.0D+00 else if ( kinda == 2 ) then al = -0.5D+00 be = -0.5D+00 if ( abs ( b - a ) <= temp ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SCQF - Fatal error!' write ( *, '(a)' ) ' |B - A| too small.' stop end if shft = ( a + b ) / 2.0D+00 slp = ( b - a ) / 2.0D+00 else if ( kinda == 3 ) then al = alpha be = alpha if ( abs ( b - a ) <= temp ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SCQF - Fatal error!' write ( *, '(a)' ) ' |B - A| too small.' stop end if shft = ( a + b ) / 2.0D+00 slp = ( b - a ) / 2.0D+00 else if ( kinda == 4 ) then al = alpha be = beta if ( abs ( b - a ) <= temp ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SCQF - Fatal error!' write ( *, '(a)' ) ' |B - A| too small.' stop end if shft = ( a + b ) / 2.0D+00 slp = ( b - a ) / 2.0D+00 else if ( kinda == 5 ) then if ( b <= 0.0D+00 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SCQF - Fatal error!' write ( *, '(a)' ) ' B <= 0' stop end if shft = a slp = 1.0D+00 / b al = alpha be = 0.0D+00 else if ( kinda == 6 ) then if ( b <= 0.0D+00 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SCQF - Fatal error!' write ( *, '(a)' ) ' B <= 0.' stop end if shft = a slp = 1.0D+00 / sqrt ( b ) al = alpha be = 0.0D+00 else if ( kinda == 7 ) then al = alpha be = 0.0D+00 if ( abs ( b - a ) <= temp ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SCQF - Fatal error!' write ( *, '(a)' ) ' |B - A| too small.' stop end if shft = ( a + b ) / 2.0D+00 slp = ( b - a ) / 2.0D+00 else if ( kinda == 8 ) then if ( a + b <= 0.0D+00 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SCQF - Fatal error!' write ( *, '(a)' ) ' A + B <= 0.' stop end if shft = a slp = a + b al = alpha be = beta else if ( kinda == 9 ) then al = 0.5D+00 be = 0.5D+00 if ( abs ( b - a ) <= temp ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SCQF - Fatal error!' write ( *, '(a)' ) ' |B - A| too small.' stop end if shft = ( a + b ) / 2.0D+00 slp = ( b - a ) / 2.0D+00 end if p = slp**( al + be + 1.0D+00 ) do k = 1, nt st(k) = shft + slp * t(k) l = abs ( ndx(k) ) if ( l /= 0 ) then tmp = p do i = l, l + mlt(k) - 1 swts(i) = wts(i) * tmp tmp = tmp * slp end do end if end do return end subroutine sgqf ( nt, aj, bj, zemu, t, wts ) !*****************************************************************************80 ! !! SGQF computes knots and weights of a Gauss Quadrature formula. ! ! Discussion: ! ! This routine computes all the knots and weights of a Gauss quadrature ! formula with simple knots from the Jacobi matrix and the zero-th ! moment of the weight function, using the Golub-Welsch technique. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 04 January 2010 ! ! Author: ! ! Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Sylvan Elhay, Jaroslav Kautsky, ! Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of ! Interpolatory Quadrature, ! ACM Transactions on Mathematical Software, ! Volume 13, Number 4, December 1987, pages 399-415. ! ! Parameters: ! ! Input, integer NT, the number of knots. ! ! Input, real ( kind = rk ) AJ(NT), the diagonal of the Jacobi matrix. ! ! Input/output, real ( kind = rk ) BJ(NT), the subdiagonal of the Jacobi ! matrix, in entries 1 through NT-1. On output, BJ has been overwritten. ! ! Input, real ( kind = rk ) ZEMU, the zero-th moment of the weight function. ! ! Output, real ( kind = rk ) T(NT), the knots. ! ! Output, real ( kind = rk ) WTS(NT), the weights. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer nt real ( kind = rk ) aj(nt) real ( kind = rk ) bj(nt) real ( kind = rk ) t(nt) real ( kind = rk ) wts(nt) real ( kind = rk ) zemu ! ! Exit if the zero-th moment is not positive. ! if ( zemu <= 0.0D+00 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SGQF - Fatal error!' write ( *, '(a)' ) ' ZEMU <= 0.' stop end if ! ! Set up vectors for IMTQLX. ! t(1:nt) = aj(1:nt) wts(1) = sqrt ( zemu ) wts(2:nt) = 0.0D+00 ! ! Diagonalize the Jacobi matrix. ! call imtqlx ( nt, t, bj, wts ) wts(1:nt) = wts(1:nt)**2 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 ! 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 ( *, '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm ) return end