program root_main ! !***************************************************************************** ! !! ROOT_MAIN controls the interactive root finding program. ! ! ! Discussion: ! ! The primeval version of ROOT was supplied by: ! ! Professor James Gentilesco, ! Computer Science Department, ! University of Pittsburgh at Johnstown. ! ! You must supply the following: ! ! a function f, and possibly its first and second derivatives. ! see below for the form of the formulas. ! ! you will have to supply at least one, and possibly as many ! as three starting points for the method. ! some methods require two points, and that the function be ! positive at one point and negative at the other. ! ! you will also be required to choose a solution method. ! ! Methods: ! ! #1 Bisection ! 2 starting points. ! 1 formula, F(X). ! Requires initial bracketing interval. ! given xl and xr on which f changes sign, set xm=0.5*(xl+xr). ! for next pair of points, take xm, and either xl or xr, whichever ! has the opposite function sign to xm. ! ! #2 Regula falsi. ! 2 starting points. ! 1 formula, F(X). ! Requires initial bracketing interval. ! given xl and xr, set xm=( f(xl)*xr-f(xr)*xl ) / (f(xl)-f(xr)) ! keep xm, and whichever of xl and xr has f of opposite sign. ! ! #3 Modified regula falsi. ! 2 starting points. ! 1 formula, F(X). ! Requires initial bracketing interval. ! ! #4 Secant method. ! 2 starting points. ! 1 formula, F(X). ! given x(k-1) and x(k), ! set top=(f(x(k-1))*x(k)-f(x(k))*x(k-1)) ! set bot = f(x(k-1)-f(x(k)) ! and set x(k+1)=top/bot. on next step, use x(k) and x(k+1). ! (compare regula falsi, where you use the same formula to get ! x(k+1) but you then might save x(k-1) instead of x(k)). ! ! #5 Muller's method, restricted to real values. ! 3 starting points. ! 1 formula, F(X). ! ! #6 Muller's method, allowing complex values. ! 3 starting points. ! 1 formula, F(X). ! ! #7 Newton-Raphson ! 1 starting point. ! 2 formulas, F(X), F'(X). ! given x(k), set x(k+1)=x(k)-f(x(k))/f'(x(k)) ! ! #8 Steffenson's method. ! 1 starting point. ! 1 formula, F(X). ! given x(k), set xtemp=x(k)+f(x(k)), ! approximate fprime=(f(xtemp)-f(x(k)))/f(x(k)) ! set x(k+1)=x(k)-f(x(k))/fprime ! ! #9 Laguerre's method for polynomials. ! 2 starting points. ! 3 formulas, F(X), F'(X), F''(X). ! f must be a polynomial, and you must supply the degree of f. ! ! #10 Fixed point iteration. ! 1 starting point. ! 1 formula, F(X). ! for fixed point iteration, we try to solve the equation f(x)=0 ! by setting g(x)=x-f(x), and searching for a point x* where ! g(x*)=x*, or g(x*)=x*-f(x*)=x* or f(x*)=0. note that ! for any function f(x), there are many other ways of setting ! up a fixed point iteration. very often, a poor formulation ! will cause the iteration to diverge. ! ! #11 Halley's method. ! 1 starting point; ! 3 formulas, F(X), F'(X), F''(X). ! given x(k), set top=f'(x(k))*f(x(k)), ! set bot=f'(x(k))*f'(x(k))-0.5*f''(x(k))*f(x(k)) ! set x(k+1)=x(k)-top/bot ! ! #12 Brent's method ! 2 starting points; ! 1 formula F(X). ! Requires initial bracketing interval. ! implicit none ! integer, parameter :: maxfix = 80 integer, parameter :: maxmet = 12 integer, parameter :: maxrpn = 320 ! character com character ( len = 80 ) command real fatol integer ierror character ( len = maxfix ) iform0 character ( len = maxfix ) iform1 character ( len = maxfix ) iform2 character ( len = maxfix ) iformg integer ifrm integer iloc character ( len = maxfix ) infix integer irpn integer itmax integer lchar integer method integer mform integer mstart character ( len = 35 ) names(maxmet) character ( len = 10 ) namvar integer nform(maxmet) integer nstart(maxmet) logical s_eqi real value real x1 real x2 real x3 real xatol real xrtol ! common /rpncom/ irpn(320) ! call timestamp ( ) ! ! initialize ! namvar = ' ' ifrm = 0 com = 'I' call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, namvar, value ) com = 'A' namvar = 'X' call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, namvar, value ) iform0 = ' ' iform1 = ' ' iform2 = ' ' iformg = ' ' ! ! Set method data ! nform(1) = 1 nform(2) = 1 nform(3) = 1 nform(4) = 1 nform(5) = 1 nform(6) = 1 nform(7) = 2 nform(8) = 1 nform(9) = 3 nform(10) = 1 nform(11) = 3 nform(12) = 1 names(1) = 'Bisection' names(2) = 'Regula Falsi' names(3) = 'Modified Regula Falsi' names(4) = 'Secant method' names(5) = 'Muller''s method (real values)' names(6) = 'Muller''s method (complex values)' names(7) = 'Newton''s method' names(8) = 'Steffenson''s method' names(9) = 'Laguerre''s method for polynomials' names(10) = 'Fixed point iteration' names(11) = 'Halley''s method' names(12) = 'Brent''s method' nstart(1) = 2 nstart(2) = 2 nstart(3) = 2 nstart(4) = 2 nstart(5) = 3 nstart(6) = 3 nstart(7) = 1 nstart(8) = 1 nstart(9) = 2 nstart(10) = 1 nstart(11) = 1 nstart(12) = 2 fatol = 0.000001E+00 x1 = 1.0E+00 x2 = 2.0E+00 x3 = 3.0E+00 itmax = 25 method = 1 mform = 1 mstart = 2 namvar = ' ' xatol = 0.000001E+00 xrtol = 0.000001E+00 ! ! Say hello. ! write ( *, * ) write ( *, * ) 'ROOT' write ( *, * ) ' An interactive program for testing simple' write ( *, * ) ' root finding algorithms.' do write ( *, * ) ' ' write ( *, * ) 'Command? (Type H for help)' read ( *, '(a)' ) command call s_blank_delete ( command ) ! ! F ! F = function ! FX ! FX = function ! F(X) ! F(X) = function ! if ( s_eqi ( command, 'F' ) .or. & s_eqi ( command(1:2), 'F=' ) .or. & s_eqi ( command, 'FX' ) .or. & s_eqi ( command(1:3), 'FX=' ) .or. & s_eqi ( command, 'F(X)' ) .or. & s_eqi ( command(1:5), 'F(X)=' ) ) then iloc = index ( command, '=' ) if ( iloc == 0 ) then write ( *, * ) 'Enter the formula for F(X)' read ( *, '(a)' ) iform0 else iform0 = command(iloc+1:) end if infix = iform0 ifrm = 1 com = 'F' call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, namvar, value ) ! ! F' = derivative ! else if ( s_eqi ( command, 'F''' ) .or. & s_eqi ( command(1:3), 'F''=' ) ) then iloc = index ( command, '=' ) if ( iloc == 0 ) then write ( *, * ) 'Enter the formula for F''(X)' read ( *, '(a)' ) iform1 else iform1 = command(iloc+1:) end if infix = iform1 ifrm = 2 com = 'F' call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, namvar, value ) ! ! F" = second derivative ! else if ( s_eqi ( command(1:2), 'F"' ) ) then iloc = index ( command, '=' ) if ( iloc == 0 ) then write ( *, * ) 'Enter the formula for F"(X)' read ( *, '(a)' ) iform2 else iform2 = command(iloc+1:) end if infix = iform2 ifrm = 3 com = 'F' call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, namvar, value ) ! ! F'' = second derivative ! else if ( s_eqi ( command(1:3), 'F''''' ) ) then iloc = index ( command, '=' ) if ( iloc == 0 ) then write ( *, * ) 'Enter the formula for F"(X)' read ( *, '(a)' ) iform2 else iform2 = command(iloc+1:) end if infix = iform2 ifrm = 3 com = 'F' call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, namvar, value ) ! ! FATOL = Function absolute tolerance. ! else if ( s_eqi ( command(1:5), 'FATOL' ) ) then if ( command(6:6) /= '=' ) then write ( *, * ) 'Enter FATOL' read ( *, '(a)' ) command else command(1:5) = ' ' end if call s_to_r ( command, fatol, ierror, lchar ) write ( *, * ) ' ' write ( *, * ) 'Function absolute tolerance FATOL set to ', fatol ! ! GO ! else if ( s_eqi ( command, 'GO' ) ) then call solve ( fatol, itmax, method, mform, x1, x2, x3, xatol, xrtol ) ! ! G = fixed point function. ! else if ( s_eqi ( command(1:1), 'G' ) ) then iloc = index ( command, '=' ) if ( iloc == 0 ) then write ( *, * ) 'Enter the formula for G(X)' read ( *, '(a)' ) iformg else iformg = command(iloc+1:) end if infix = iformg ifrm = 4 com = 'F' call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, namvar, value ) ! ! HELP ! else if ( s_eqi ( command(1:1), 'H' ) ) then call help ! ! ITMAX = Iteration maximum. ! else if ( s_eqi ( command(1:5), 'ITMAX' ) ) then if ( command(6:6) /= '=' ) then write ( *, * ) 'Enter iteration maximum' read ( *, '(a)' ) command else command(1:6) = ' ' end if call s_to_i ( command, itmax, ierror, lchar ) write ( *, * ) ' ' write ( *, * ) 'Iteration maximum ITMAX set to ', itmax ! ! P ! else if ( s_eqi ( command(1:1), 'P' ) ) then call prlist ( fatol, iform0, iform1, iform2, iformg, itmax, & method, mstart, names, maxmet, x1, x2, x3, xatol, xrtol ) ! ! M ! else if ( s_eqi ( command(1:1), 'M' ) ) then call method_choice ( method, mform, mstart, names, nform, maxmet, & nstart, ierror ) ! ! Q = QUIT ! else if ( s_eqi ( command(1:1), 'Q' ) ) then exit ! ! X1 = First starting point. ! else if ( s_eqi ( command(1:2), 'X1' ) ) then if ( command(3:3) /= '=' ) then write ( *, * ) 'Enter X1' read ( *, '(a)' ) command else command(1:3) = ' ' end if call s_to_r ( command, x1, ierror, lchar ) write ( *, * ) ' ' write ( *, * ) 'First starting point X1 set to ', x1 ! ! X2 = Second starting point. ! else if ( s_eqi ( command(1:2), 'X2' ) ) then if ( command(3:3) /= '=' ) then write ( *, * ) 'Enter X2' read ( *, '(a)' ) command else command(1:3) = ' ' end if call s_to_r ( command, x2, ierror, lchar ) write ( *, * ) ' ' write ( *, * ) 'Second starting point X2 set to ', x2 ! ! X3 = Third starting point. ! else if ( s_eqi ( command(1:2), 'X3' ) ) then if ( command(3:3) /= '=' ) then write ( *, * ) 'Enter X3' read ( *, '(a)' ) command else command(1:3) = ' ' end if call s_to_r ( command, x3, ierror, lchar ) write ( *, * ) ' ' write ( *, * ) 'Third starting point X3 set to ', x3 ! ! XATOL = X absolute tolerance. ! else if ( s_eqi ( command(1:5), 'XATOL' ) ) then if ( command(6:6) /= '=' ) then write ( *, * ) 'Enter XATOL' read ( *, '(a)' ) command else command(1:5) = ' ' end if call s_to_r ( command, xatol, ierror, lchar ) write ( *, * ) ' ' write ( *, * ) 'X absolute tolerance XATOL set to ', xatol ! ! XRTOL = X relative tolerance. ! else if ( s_eqi ( command(1:5), 'XRTOL' ) ) then if ( command(6:6) /= '=' ) then write ( *, * ) 'Enter XRTOL' read ( *, '(a)' ) command else command(1:5) = ' ' end if call s_to_r ( command, xrtol, ierror, lchar ) write ( *, * ) ' ' write ( *, * ) 'X relative tolerance XRTOL set to ', xrtol else write ( *, * ) 'Unrecognized command!' end if end do write ( *, * ) ' ' write ( *, * ) 'ROOT:' write ( *, * ) ' Normal end of execution.' stop end subroutine bisection ( fatol, itmax, x1, x2, xatol, xrtol ) ! !******************************************************************************* ! !! BISECTION carries out the bisection algorithm. ! ! ! Modified: ! ! 19 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real FATOL, the absolute error tolerance for F(X). ! ! Input, integer ITMAX, the maximum number of steps allowed. ! ! Input, real X1, X2, two points where the function differs in sign. ! ! Input, real XATOL, absolute error tolerance for the root. ! implicit none ! real fatol real fxl real fxm real fxr integer iterate integer itmax real x_ave real x1 real x2 real xatol real xl real xm real xr real xrtol ! write ( *, * ) ' ' write ( *, * ) ' ' write ( *, * ) 'The bisection method' write ( *, * ) ' ' write ( *, * ) ' ' xl = x1 xr = x2 call f ( xl, fxl, 0 ) call f ( xr, fxr, 0 ) ! ! Check for change of sign. ! if ( fxl * fxr > 0.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'WARNING:' write ( *, * ) ' ' write ( *, * ) ' This method requires a change of sign ' write ( *, * ) ' interval, but your function does NOT change ' write ( *, * ) ' sign over the given interval, so this method ' write ( *, * ) ' cannot be used.' write ( *, * ) ' ' write ( *, * ) 'SUGGESTION:' write ( *, * ) ' ' write ( *, * ) ' Try a different method;' write ( *, * ) ' Try changing X1 or X2;' write ( *, * ) ' Perhaps your function F(X) is wrong.' return end if write ( *, * ) ' Step Left Middle Right' iterate = 0 write ( *, * ) ' ' write ( *, '(i4,'' X'',g16.8,16x,f16.8)' ) iterate, xl, xr write ( *, '(4x,''FX'',g16.8,16x,g16.8)' ) fxl, fxr ! ! Take another step of the iteration. ! do iterate = iterate + 1 x_ave = ( abs ( xr ) + abs ( xl ) ) / 2.0E+00 if ( iterate > itmax ) then write ( *, * ) ' ' write ( *, * ) ' Maximum number of steps taken.' exit end if if ( abs ( xr - xl ) <= xatol ) then write ( *, * ) ' ' write ( *, * ) ' Absolute convergence of the X interval.' exit end if if ( abs ( xr - xl ) <= xrtol * x_ave ) then write ( *, * ) ' ' write ( *, * ) ' Absolute convergence of the X interval.' exit end if if ( abs ( fxl ) <= fatol ) then write ( *, * ) ' ' write ( *, * ) ' Absolute convergence of |F(X)|.' exit end if if ( abs ( fxr ) <= fatol ) then write ( *, * ) ' ' write ( *, * ) ' Absolute convergence of |F(X)|.' exit end if ! ! Compute the next iterate. ! xm = ( xl + xr ) / 2.0E+00 call f ( xm, fxm, 0 ) write ( *, * ) ' ' write ( *, '(i4,'' X'',3f16.8)' ) iterate, xl, xm, xr write ( *, '(4x,''FX'',3g16.8)' ) fxl, fxm, fxr ! ! Prepare for the next step. ! if ( fxl * fxm > 0.0E+00 ) then xl = xm fxl = fxm else xr = xm fxr = fxm end if end do return end subroutine brent ( fatol, itmax, x1, x2, xatol, xrtol ) ! !******************************************************************************* ! !! BRENT implements the Brent bisection-based zero finder. ! ! ! Reference: ! ! Richard Brent, ! Algorithms for Minimization without Derivatives, ! Prentice Hall, 1973. ! ! Modified: ! ! 20 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real FATOL, an absolute error tolerance for the function ! value of the root. If an approximate root X satisfies ! ABS ( F ( X ) ) <= FATOL, then X will be accepted as the ! root and the iteration will be terminated. ! ! Input, integer ITMAX, the maximum number of steps allowed. ! ! Input, real X1, X2, two points at which the function differs in sign. ! ! Input, real XATOL, XRTOL, absolute and relative error tolerances ! for the root. ! implicit none ! real d real e real fatol real fxa real fxb real fxc integer iterate integer itmax real p real q real r real s real x_ave real x_int real x1 real x2 real xa real xb real xc real xm real xatol real xrtol ! ! Initialization. ! write ( *, * ) ' ' write ( *, * ) 'Brent''s Method' write ( *, * ) ' ' write ( *, * ) 'Step XA XB F(XA) F(XB)' write ( *, * ) ' ' iterate = 0 xa = x1 xb = x2 call f ( xa, fxa, 0 ) call f ( xb, fxb, 0 ) ! ! Check for change of sign. ! if ( fxa * fxb > 0.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'WARNING:' write ( *, * ) ' ' write ( *, * ) ' This method requires a change of sign ' write ( *, * ) ' interval, but your function does NOT change ' write ( *, * ) ' sign over the given interval, so this method ' write ( *, * ) ' cannot be used.' write ( *, * ) ' ' write ( *, * ) 'SUGGESTION:' write ( *, * ) ' ' write ( *, * ) ' Try a different method;' write ( *, * ) ' Try changing X1 or X2;' write ( *, * ) ' Perhaps your function F(X) is wrong.' return end if xc = xa fxc = fxa d = xb - xa e = d do write ( *, '(i4,2x,2g16.8,2g14.6)' ) iterate, xb, xc, fxb, fxc iterate = iterate + 1 if ( iterate > itmax ) then write ( *, * ) ' ' write ( *, * ) ' Maximum number of steps taken.' exit end if ! ! If necessary, swap data so that |F(XB)| <= |F(XC)|. ! if ( abs ( fxc ) < abs ( fxb ) ) then call r_swap ( xb, xc ) call r_swap ( fxb, fxc ) end if ! ! XM is the halfwidth of the current change-of-sign interval. ! x_int = xc - xb x_ave = ( abs ( xc ) + abs ( xb ) ) / 2.0E+00 xm = 0.5E+00 * ( xc - xb ) if ( abs ( x_int ) <= xatol ) then write ( *, * ) ' ' write ( *, * ) ' Absolute convergence of the X interval.' exit end if if ( abs ( x_int ) <= xrtol * x_ave ) then write ( *, * ) ' ' write ( *, * ) ' Relative convergence of the X interval.' exit end if if ( abs ( fxb ) <= fatol ) then write ( *, * ) ' ' write ( *, * ) ' Absolute convergence of |F(X)|.' exit end if ! ! See if a bisection is forced. ! if ( abs ( e ) < xatol .or. abs ( fxa ) <= abs ( fxb ) ) then d = xm e = d else s = fxb / fxa ! ! Linear interpolation. ! if ( xa == xc ) then p = 2.0E+00 * xm * s q = 1.0E+00 - s ! ! Inverse quadratic interpolation. ! else q = fxa / fxc r = fxb / fxc p = s * ( 2.0E+00 * xm * q * ( q - r ) - ( xb - xa ) * ( r - 1.0E+00 ) ) q = ( q - 1.0E+00 ) * ( r - 1.0E+00 ) * ( s - 1.0E+00 ) end if if ( p > 0.0E+00 ) then q = - q else p = - p end if s = e e = d if ( 2.0E+00 * p >= 3.0E+00 * xm * q - abs ( xatol * q ) .or. & p >= abs ( 0.5E+00 * s * q ) ) then d = xm e = d else d = p / q end if end if ! ! Save in XA, FXA the previous values of XB, FXB. ! xa = xb fxa = fxb ! ! Compute the next iterate. ! if ( abs ( d ) > xatol ) then xb = xb + d else if ( xm > 0.0E+00 ) then xb = xb + xatol else if ( xm <= 0.0E+00 ) then xb = xb - xatol end if call f ( xb, fxb, 0 ) ! ! If the new FXB has the same sign as FXC, then replace XC by XA. ! if ( ( fxb > 0.0E+00 .and. fxc > 0.0E+00 ) .or. & ( fxb < 0.0E+00 .and. fxc < 0.0E+00 ) ) then xc = xa fxc = fxa d = xb - xa e = d end if end do return end subroutine c_swap ( x, y ) ! !******************************************************************************* ! !! C_SWAP swaps two complex values. ! ! ! Modified: ! ! 26 July 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, complex X, Y. On output, the values of X and ! Y have been interchanged. ! implicit none ! complex x complex y complex z ! z = x x = y y = z return end subroutine ch_cap ( c ) ! !******************************************************************************* ! !! CH_CAP capitalizes a single character. ! ! ! 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 ) ! !******************************************************************************* ! !! CH_EQI is a case insensitive comparison of two characters for equality. ! ! ! Examples: ! ! CH_EQI ( 'A', 'a' ) is .TRUE. ! ! 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 function ch_is_alpha ( c ) ! !******************************************************************************* ! !! CH_IS_ALPHA returns TRUE if C is an alphabetic character. ! ! ! Modified: ! ! 05 August 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character C, a character to check. ! ! Output, logical C_IS_ALPHA is TRUE if C is an alphabetic character. ! implicit none ! character c logical ch_is_alpha ! if ( ( lle ( 'a', c ) .and. lle ( c, 'z' ) ) .or. & ( lle ( 'A', c ) .and. lle ( c, 'Z' ) ) ) then ch_is_alpha = .true. else ch_is_alpha = .false. end if return end subroutine ch_to_digit ( c, digit ) ! !******************************************************************************* ! !! CH_TO_DIGIT returns the integer value of a base 10 digit. ! ! ! Example: ! ! C DIGIT ! --- ----- ! '0' 0 ! '1' 1 ! ... ... ! '9' 9 ! ' ' 0 ! 'X' -1 ! ! 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 integer 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 comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, namvar, value ) ! !******************************************************************************* ! !! COMRPN can translate formulas you type in and evaluate them. ! ! ! This means that you can design interactive FORTRAN programs ! which can input their formulas at run time, rather than ! recompiling every time you want a new formula. ! ! If you declare vectors or matrices, you may reference entries ! such as X(3) or A(3,2) or even X(2+3). ! ! ! Formulas are made from operators, constants, punctuation, ! variables, and functions. ! ! ! The list of legal operators includes: ! ! + - * / ** ^ = \ ! ! * means multiplication, and is standard matrix multiplication if ! both arguments are (conformable) matrices. If one argument is a ! scalar, it multiplies all the entries of the other argument. ! If both arguments are row or column vectors, the dot product ! is taken. ! ! Thus, ! ! x is a scalar equal to 2, ! y is a scalar equal to 3, ! u is a vector equal to (1, 2, 3), ! v is a vector equal to (1, 0, 2), ! ! then ! ! 2*3 results in 6. ! ! x*y results in 6. ! ! u*v results in 7. ! ! x*v would result in the vector (2, 0, 4). ! ! ! / means division, as in A/B, but for matrices the only allowed ! form requires that B be a scalar, in which case each element of ! A is divided by B. ! ! ! + and - are allowed for pairs of scalars, vectors or matrices of ! the same order, and also for a square matrix and a scalar, ! in which case A+2 is interpreted as A+2*I. ! ! ! ** or ^ means exponentiation. In the scalar case, any ! nonnegative scalar can be taken to any power, positive, zero, or ! negative. A negative scalar may only be taken to an integer ! power. ! ! ** is also supported for square matrices, which can be taken ! to any nonnegative integer power. Nonsingular matrices may ! also be taken to negative integer powers. ! ! ! = is used to assign values. For vectors or matrices, this may ! also be used to assign a single entry, as in A(3,2)=7. ! ! ! / is standard scalar division. It is used in a formula like ! ! X/Y ! ! or ! ! (A*B)/(X+1) ! ! Normally, X would be a scalar quantity. However, as NONSTANDARD ! shorthand, you may 'divide' a matrix by a scalar, in which case ! each entry of the matrix is divided by the scalar, e.g. ! ! A/2 = (1/2)*A. ! ! MATALG will NOT allow you to use the "/" operator to represent ! multiplication by an inverse matrix. Thus, if A is a matrix, ! the statement ! ! B/A ! ! is illegal. However, the statement ! ! INV(A)*B ! ! will compute the inverse of A and multiply by B, and the ! statement ! ! A \ B ! ! will compute the LU factorization of A, and use that to ! compute the inverse of A times B. ! ! ! Constants: ! ! ! Real and integer constants may appear in your formulas, as well ! as the symbol 'PI' and the machine unit roundoff number 'EPS' ! which is the power of two such that 1+EPS>1 but 1+(EPS/2)=1. ! ! ! Punctuation: ! ! ! Blanks may be used anywhere. ! ! Commas are used to separate arguments in a function, such as ! MAX(X,Y). ! ! Parentheses are used: ! ! to group quantities: (X+Y)*Z ! ! to reference an entry of a vector or matrix: A(2,2) ! ! to mark the argument of a function: SIN(X), MAX(X,Y) ! ! ! Functions and operators: ! ! ! S, S1, S2: Arguments which may only be scalar. ! V : Arguments may only be a scalar or vector. ! M : Arguments may only be scalar or square matrix. ! MV : Arguments may only be scalar, vector, or square matrix. ! * : Arguments which may be scalar, vector, or matrix. ! I : Arguments which may only be positive integers. ! ! ABS(*) Absolute value of *. ! ! ACOS(S) The arc cosine of S. ! -1 <= S <=1. ! ! ASIN(S) The arc sine of S. ! -1 <= S <=1. ! ! ATAN(S) The arc tangent of S. ! ! ATAN2(S1,S2) Arc tangent of (S1/S2) ! Correctly computes ATAN2(1.0,0.0)=PI/2. ! ! COS(S) The cosine of S, with S measured in radians. ! ! COSD(S) The cosine of S, with S measured in degrees. ! ! COSH(S) Hyperbolic cosine of S. ! ! DET(M) Determinant of square matrix M. ! ! DIAG(*) The diagonal entries of *, stored in a vector. ! ! E The base of the natural logarithm system. ! E=2.71828... ! You may not change the value of E. ! ! EPS The machine epsilon, or unit roundoff number. ! EPS is the power of 2 such that ! ! 1 < 1+EPS and 1=1+(EPS/2). ! ! You may not change the value of EPS. ! ! EVAL(M) Real and imaginary parts of eigenvalues of matrix ! M, stored in a matrix of N rows and 2 columns, ! with the real parts in column 1, and the ! imaginary parts in column 2. ! ! EVEC(M) Eigenvectors of square matrix M, stored in a ! square matrix of same size as M. ! ! (Eigenvectors corresponding to a complex pair: ! If eigenvalues j and j+1 are a complex pair, ! then the eigenvector for eigenvalue j is ! column j + i*column j+1, and the eigenvector ! for eigenvalue j+1 is column j - i*column j+1.) ! ! EXP(S) Exponential of S. ! ! GCF(I,J) Greatest common factor of two integers. ! ! HILBERT(I) The Hilbert matrix of order I. ! ! HILBINV(I) The inverse of the Hilbert matrix of order I. ! ! HOUSE(V) The Householder elementary reflector matrix H, ! defined as ! ! H = I - 2*(V * TRANSPOSE(V)) / (TRANS(V)*V) ! ! ID(I) The matrix identity of order I. ! ID(3) is the 3 by 3 identity, for instance. ! ! INT(*) Truncates real values to their integer part. ! INT(1.0) = INT(1.1) = INT(1.9) = 1. ! ! INV(M) The inverse matrix of the square matrix M. ! If M is singular, there will be no inverse. ! ! IVAL(M) Imaginary parts of eigenvalues of square matrix M ! stored in a column vector. ! ! LCM(I,J) Least common multiple of two integers. ! ! LN(S) The natural logarithm of S. ! S must be greater than 0. ! ! LOG(S) The natural logarithm of S. ! S must be greater than 0. ! ! LOG10(S) The logarithm base 10 of S. ! S must be greater than 0. ! ! LOG2(S) The logarithm base 2 of S. ! S must be greater than 0. ! ! MAX(S1,S2) The maximum of S1 and S2. ! ! MIN(S1,S2) Minimum of S1 and S2. ! ! NEG(*) Changes sign of *. ! ! NINT(*) Nearest integer value to entries of *. ! ! NORM0(*) Maximum or infinity norm of a vector or matrix. ! NORM0(S) = ABS(S) ! NORM0(V) returns the maximum of the absolute ! values of the entries of V. ! NORM0(M) returns the maximum of the sum of the ! absolute values of the entries in each row. ! ! NORM1(*) The L1-norm of a vector or matrix. ! NORM1(S) = ABS(S). ! NORM1(V) returns the sum of the absolute values ! of the entries of V. ! NORM0(M) returns the maximum of the sum of the ! absolute values of the entries in each column. ! ! NORM2(MV) L2-norm, Euclidean norm or root-mean-square norm ! of a vector or a square matrix M. NORM2(V) returns ! the square root of the sum of the squares of the ! entries of V. NORM2(M) returns the maximum magnitude ! of the eigenvalues of Transpose(M)*M. ! ! NORMF(*) The Frobenius norm. NORMF(M) returns the square ! root of the sum of the squares of all the entries ! of a matrix M. NORMF may also be applied to a ! vector, giving the same results as NORM2(V). ! ! NORMS(MV) The spectral "norm". NORMS(M) returns the ! value of the maximum of the absolute values of ! all the eigenvalues of the square matrix M. ! ! Note that the spectral norm is NOT a true ! norm, although it is frequently used as though ! it were. ! ! PI 3.14159265358979... ! You may not change the value of PI. ! ! POLY(V,M) Polynomial evaluation. V contains the ! coefficients of the polynomial, with V(1) the ! coefficient of M**(N-1) and V(N) the constant ! term. ! ! M is the scalar or square matrix argument of the ! polynomial. ! ! RAN(*) Fills * with random numbers between 0 and 1. ! ! ROUND(*) Replaces all "small" entries of * by 0. ! Here "small" means less than RTOL in absolute ! value. ! ! RTOL The tolerance for the ROUND command, initially set ! to EPS. To change the value of RTOL, simply ! reassign it with a statement like "RTOL=0.1E-6". ! ! RVAL(M) Real parts of eigenvalues of square matrix M, ! stored in a column vector. ! ! SIN(S) The sine of S, with S measured in radians. ! ! SIND(S) The sine of S, with S measured in degrees. ! ! SINH(S) Hyperbolic sine of S. ! ! SQRT(S) Square root of S. S must be non-negative. ! ! STEP(S) Heavyside step function. ! ! STEP=0 if S < 0, ! STEP=1 if 0 <= S. ! ! TAN(S) Tangent of S. ! ! TAND(S) Tangent of S measured in degrees. ! ! TANH(S) Hyperbolic tangent of S. ! ! TRACE(*) Sum of diagonal elements of *. ! ! TRANS(*) Transpose of matrix or vector. ! ! A=TRANS(A) is a legal formula if A is square. ! ! ZEROS(I) Matrix zero of order I. ! ZEROS(4) is the 4 by 4 zero matrix. ! ! I! Factorial of I. I should be an integer from 0 to ! 25. 0!=1, 1!=1, 2!=1*2=2, 3!=1*2*3=6, and so on. ! For large I, the exact value is not returned. ! !******************************************************************************* ! ! To add a new function to COMRPN: ! ! Update the information in the IOPSYM, IPRSYM, and SYMBOL ! arrays. ! ! Insert text into FUNVAL or FUNSCL which evaluates the ! function, and determines its dimensions. ! !******************************************************************************* ! ! Parameters: ! ! COM Input, character COM. ! ! The command to be executed, which determines what other ! input must be supplied. ! ! ! COM='A' Add (=declare) a variable: ! ! Causes the program to add the variable named in NAMVAR ! to its list of symbols. From now on, it may be used in ! formulas, and values may be assigned to it. If you ! declare a name which already exists, then if that name ! was already declared by the program, your command is ! ignored. Otherwise, your new definition overrides your ! old one. ! ! Input for the "A" command includes NAMVAR, NROW, ! and NCOL. ! ! ! COM='D' delete a variable: ! ! This command can be used to free up memory. You ! must also give the name of the variable to be deleted. ! The variable must be one you have already defined. ! ! Input for the "D" command includes NAMVAR. ! ! ! COM='E' evaluate a formula: ! ! Using current values of variables, evaluate the formula. ! The formula number is given by IFRM. The value of the ! formula is returned in VALUE, and its dimensions in NROW ! and NCOL. ! ! Input to the "E" command includes IFRM. ! Output from the "E" command includes VALUE. ! ! ! COM='F', Formula is to be compiled: ! ! The user assigns an index, IFRM, to the formula, so that ! it can be referred to later. The formula is stored in ! INFIX on input. The formula may refer to variables which ! have not yet been declared. ! ! Output from the "F" command includes IFRM. ! ! ! COM='G', Formula is to be compiled, all variables ! in the formula have already been declared: ! ! The same as 'F' except that the formula may not refer to ! any variables which have not been declared. ! ! Output from the "G" command includes IFRM. ! ! ! COM='I' Initialize: ! ! This must be the first call, to initialize COMRPN's internal ! data. Also, if old data is to be flushed out, this command ! can be used. ! ! ! COM='R' Read value of symbol: ! ! The program returns the current value of the variable ! NAMVAR. VALUE will contain the value, and NROW and NCOL ! the dimensions. ! ! ! COM='V' assign value to symbol: ! ! This command assigns the value of the numbers in VALUE ! to the variable NAMVAR, with NROW and NCOL specifying ! the dimension of the variable. ! ! IERROR Output, integer IERROR. ! ! IERROR is an error flag. ! ! 0 Means no error occurred. ! ! 1 means some kind of error. These errors are usually ! signaled by a printed message. They can include the ! following problems: ! ! On adding a variable, if the name was already in use, ! either by you or the program, the program returns ! with this error warning. ! ! On setting a variable, if the name supplied does ! not represent any recognized variable, the program ! ignores the command and returns. ! ! On compiling a formula, the program will return with ! an error warning if any of the following occur: ! ! The formula does not compile. An unknown variable, ! missing parentheses, garbled information or mistyping ! may be responsible. For that matter, the compiler ! may fail to compile perfectly good, but complicated ! formulas. ! ! There is not enough space in IRPN to store the RPN code. ! In this case, you must either re-initialize, or increase ! the storage available in IRPN and signal this with ! an increased value of MAXRPN. ! ! On evaluating a formula, an error can occur if there ! is no formula corresponding to the given index IFRM. ! ! Certain arithmetic errors, such as division ! by zero, may be caught and flagged by the compiler. ! ! IFRM Input/output, integer IFRM. ! ! COMRPN can store more than one formula internally, at ! one time. COMRPN refers to the formulas by an index ! number. When a formula is entered, with the "F" or "G" ! command, it is assigned an index number, whose value is ! returned to the user. Then, whenever the formula is to ! be evaluated, the user must input the value of IFRM, so ! that COMRPN knows which formula to evaluate. In many ! cases, only one formula is being used, so that IFRM is ! usually just 1. ! ! For the "E" command, IFRM is the index of the formula to ! evaluate. ! ! For the "F" or "G" commands, the index assigned to the ! formula. ! ! INFIX character ( len = * ) INFIX, for the F or G commands, ! contains the user's formula to be compiled and evaluated. ! ! IRPN integer IRPN(MAXRPN), for the F, G or E commands, ! a vector used to store the compiled versions of the input ! formulas. ! ! MAXRPN Input, integer MAXRPN. ! ! The maximum length of IRPN. 80 should be enough. ! ! NAMVAR Input, character ( len = * ) NAMVAR, used for the A, D, L, R and ! V commands, containing the name of the variable to be ! added, deleted, listed, read or valued. ! ! NAMVAR should not be a blank, nor should it be the name ! 'VOID'. ! ! VALUE real VALUE(1,1), used to pass values in and out ! of the program. The actual size of the object in VALUE ! is (NROW,NCOL). Thus to pass a scalar, set VALUE(1,1) ! to the value and set NROW=NCOL=1. To pass a 2 by 3 ! matrix, set the values of VALUE and set NROW=2, NCOL=3. ! ! When the program is asked to return a value in VALUE, it ! sets NROW and NCOL as well. ! ! For the V command, VALUE is the input value to be ! assigned to a variable. ! ! For the E command, VALUE is the output value of the formula. ! ! For the R command, VALUE is the output value of the variable. ! !******************************************************************************* ! ! Internal variables: ! ! IOPSYM Internal, integer IOPSYM(MAXSYM), contains, for each ! symbol, the number of operands. In particular, if ! a symbol represents a constant, IOPSYM(I) is 0. ! If a symbol represents a unary operator such as ABS, ! IOPSYM(I) is 1. IOPSYM may be as large as 3. ! ! IPRSYM Internal, integer IPRSYM(MAXSYM), contains, for each ! symbol, the operator "priority" of that symbol. This ! is simply to help COMRPN deal with ambiguous formulas like ! A*B+C. ! ! IPVAL Internal, integer IPVAL(MAXSYM), contains, for each ! symbol, the address in VALSYM where the value of the ! symbol is stored, if it is a scalar. If the symbol ! represents a vector or matrix, IPVAL(IARG) points to ! the location of the first entry. ! ! MAXSYM Internal, integer MAXSYM, controls the maximum number ! of symbols allowed. This includes permanent symbols ! defined by the program, as well as symbols the user ! declares while running the program. ! !******************************************************************************* ! implicit none ! integer, parameter :: maxsym = 150 integer, parameter :: maxval = 2500 integer, parameter :: maxlen = 20 ! integer maxrpn ! character com integer ierror integer, save :: ifinis integer, save :: ifree integer ifreesv integer ifrm integer ihi integer ilo integer imin integer implic integer, save :: indx1 integer, save :: indx2 integer, save :: ineg character ( len = * ) infix integer, dimension ( 80 ), save :: intsym integer, dimension ( maxsym ), save :: iopsym integer, dimension ( maxsym ), save :: iprsym integer, dimension ( maxsym ), save :: ipval integer, save :: irad integer irpn(maxrpn) integer, save :: irtol integer, dimension ( maxsym ), save :: istack integer maxfix integer maxfrm character ( len = * ) namvar character ( len = maxlen ) namvr integer, save :: nints integer nrpn integer, save :: nsym integer, save :: nsymp integer, save :: nsyms integer nuvar logical s_eqi logical s_paren_check character ( len = maxlen ), dimension ( maxsym ), save :: symbol real value real, dimension ( maxval ), save :: valsym ! ierror = 0 implic = 0 namvr = namvar maxfix = len ( infix ) if ( maxfix > 0 ) then maxfrm = ( maxrpn / maxfix ) else maxfrm = 0 end if call s_blank_delete ( namvr ) if ( s_eqi ( com, 'E' ) .or. s_eqi ( com, 'G' ) ) then if ( ifrm <= 0 .or. ifrm > maxfrm ) then write ( *, * ) ' ' write ( *, * ) 'COMRPN - Error!' write ( *, * ) ' Illegal formula index= ', ifrm ierror = 1 return end if end if ! ! COM=I initialize. ! if ( s_eqi ( com, 'I' ) ) then call inicom ( ifinis, ifree, indx1, indx2, ineg, infix, intsym, iopsym, & iprsym, ipval, irad, irpn, irtol, istack, maxrpn, maxsym, & maxval, nints, nsym, nsymp, nsyms, symbol, valsym, value ) return ! ! COM=A, Add variable. ! else if (s_eqi(com,'a') ) then call symadd ( ierror, ifree, iopsym, iprsym, ipval, & maxsym, maxval, namvr, nsym, nsymp, symbol, valsym ) if ( ierror == 0 ) then write ( *, * ) ' ' write ( *, * ) 'COMRPN - Note:' write ( *, * ) ' Added the variable ', namvr end if return ! ! COM=V, Set variable value. ! COM=R, Get variable value. ! else if ( s_eqi ( com, 'r' ) .or. s_eqi ( com, 'v' ) ) then call symbol_value ( com, ierror, ifree, iopsym, iprsym, ipval, maxsym, & maxval, namvr, nsym, nsymp, symbol, valsym, value ) return ! ! COM=F, compile formula. ! else if ( s_eqi ( com, 'f' ) ) then nuvar = 1 ! ! Remove all blanks from the formula. ! call s_blank_delete ( infix ) ! ! Reject the formula if it has no information in it. ! if ( len_trim ( infix ) <= 0 ) then ierror = 1 write ( *, * ) ' ' write ( *, * ) 'COMRPN - Fatal error!' write ( *, * ) ' The formula has zero length.' return end if ! ! Do a simple check on parentheses. ! if ( .not. s_paren_check ( infix ) ) then ierror = 1 write ( *, * ) ' ' write ( *, * ) 'COMRPN - Fatal error!' write ( *, * ) ' The formula did not pass the parentheses checks!' return end if ! ! Convert INFIX, the string of characters, into INTSYM, an array of ! integers which are indices of tokens. ! call tokens ( ierror, ifinis, ifree, implic, indx1, indx2, ineg, infix, & intsym, iopsym, iprsym, ipval, maxsym, maxval, nints, nsym, nsymp, & nuvar, symbol, valsym ) if ( ierror /= 0 ) then write ( *, * ) ' ' write ( *, * ) 'COMRPN - Fatal error!' write ( *, * ) ' Could not convert formula into tokens.' if ( implic /= 0 ) then write ( *, * ) ' Cancelling the variable ' // trim ( symbol(implic) ) symbol(implic) = 'VOID' end if return end if ! ! Convert INTSYM, which is an infix formula, into IRPN, which ! is an RPN formula. ! imin = (ifrm-1) * 80 + 1 call rpnset ( ierror, ifinis, imin, intsym, iopsym, iprsym, irpn, & istack, maxrpn, maxsym, nints, nrpn, symbol ) if ( ierror /= 0 ) then write ( *, * ) ' ' write ( *, * ) 'COMRPN - Fatal error!' write ( *, * ) ' Error! Could not compile formula.' if ( implic /= 0 ) then write ( *, * ) ' Cancelling the variable ' // trim ( symbol(implic) ) symbol(implic) = 'VOID' end if return end if ! ! Check that operators and arguments are balanced. ! ihi = imin + nrpn call rpnchk ( ierror, ihi, ilo, iopsym, irpn, maxrpn, maxsym ) if ( ierror /= 0 .or. ilo /= imin + 1 ) then write ( *, * ) ' ' write ( *, * ) 'COMRPN - Error!' write ( *, * ) ' Parentheses mismatch.' write ( *, * ) ' ILO = ', ilo write ( *, * ) ' IMIN = ', imin ierror = 1 if ( implic /= 0 ) then write ( *, * ) ' Cancelling the variable ' // trim ( symbol(implic) ) symbol(implic) = 'VOID' end if return end if ! ! For implicit definition via equality, evaluate formula to ! get dimensions. ! if ( implic /= 0 ) then go to 10 end if else if ( s_eqi ( com, 'e' ) ) then go to 10 else write ( *, * ) ' ' write ( *, * ) 'COMRPN - Fatal error!' write ( *, * ) 'Error! Unknown command = ' // com ierror = 1 end if return ! ! COM=E, Evaluate formula. ! 10 continue imin = (ifrm-1) * 80 + 1 nrpn = 80 ifreesv = ifree call rpnval ( ierror, ifree, imin, iopsym, iprsym, ipval, irad, irpn, & irtol, istack, maxrpn, maxsym, maxval, nrpn, nsym, & nsyms, symbol, valsym, value ) ifree = ifreesv if ( ierror /= 0 ) then value = 0.0E+00 write ( *, * ) ' ' write ( *, * ) 'COMRPN - Fatal error!' write ( *, * ) ' The formula could not be evaluated!' if ( implic /= 0 ) then write ( *, * ) ' Cancelling the variable ' // trim ( symbol(implic) ) symbol(implic) = 'VOID' end if return end if return end subroutine digit_to_ch ( digit, c ) ! !******************************************************************************* ! !! DIGIT_TO_CH returns the character representation of a decimal digit. ! ! ! Example: ! ! DIGIT C ! ----- --- ! 0 '0' ! 1 '1' ! ... ... ! 9 '9' ! 17 '*' ! ! Modified: ! ! 04 August 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer DIGIT, the digit value between 0 and 9. ! ! Output, character C, the corresponding character, or '*' if DIGIT ! was illegal. ! implicit none ! character c integer digit ! if ( 0 <= digit .and. digit <= 9 ) then c = char ( digit + 48 ) else c = '*' end if return end subroutine f ( x, fx, indx ) ! !******************************************************************************* ! !! F is a function evaluation interface routine. ! ! ! Discussion: ! ! F calls an interactive compiler to evaluate the function. ! ! Modified: ! ! 18 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real X, the point at which the formula is to be evaluated. ! ! Output, real FX, the value of the formula. ! ! Input, integer INDX, the index of the formula, normally ! 0 for the formula F(X), ! 1 for the formula F'(X), ! 2 for the formula F''(X), ! 3 for a fixed point formula. ! implicit none ! real fx integer ierror integer indx character ( len = 80 ) infix integer irpn integer maxrpn real value real x ! common /rpncom/ irpn(320) ! maxrpn = 320 value = x call comrpn ( 'V', ierror, indx+1, infix, irpn, maxrpn, 'X', value ) call comrpn ( 'E', ierror, indx+1, infix, irpn, maxrpn, ' ', value ) fx = value return end subroutine fc ( x, fx, indx ) ! !******************************************************************************* ! !! FC is a complex function evaluation interface routine. ! ! ! Discussion: ! ! This is just a dummy routine. We can't handle complex data now. ! ! Modified: ! ! 18 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real X, the point at which the formula is to be evaluated. ! ! Output, real FX, the value of the formula. ! ! Input, integer INDX, the index of the formula, normally ! 0 for the formula F(X), ! 1 for the formula F'(X), ! 2 for the formula F''(X), ! 3 for a fixed point formula. ! implicit none ! complex fx integer indx complex x ! return end subroutine falsi ( fatol, itmax, x1, x2, xatol, xrtol ) ! !******************************************************************************* ! !! FALSI carries out the Regula Falsi method. ! ! ! Modified: ! ! 18 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real FATOL, the absolute error tolerance for F(X). ! ! Input, integer ITMAX, the maximum number of steps allowed. ! ! Input, real X1, X2, two points where the function differs in sign. ! ! Input, real XATOL, absolute error tolerance for the root. ! implicit none ! real fatol real fxl real fxm real fxmold real fxr integer iterate integer itmax real wate1 real wate2 real x_ave real x1 real x2 real xatol real xl real xm real xr real xrtol ! ! The starting points ! xl = x1 xr = x2 call f ( xl, fxl, 0 ) call f ( xr, fxr, 0 ) ! ! Check for change of sign. ! if ( fxl * fxr > 0.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'WARNING:' write ( *, * ) ' ' write ( *, * ) ' This method requires a change of sign ' write ( *, * ) ' interval, but your function does NOT change ' write ( *, * ) ' sign over the given interval, so this method ' write ( *, * ) ' cannot be used.' write ( *, * ) ' ' write ( *, * ) 'SUGGESTION:' write ( *, * ) ' ' write ( *, * ) ' Try a different method;' write ( *, * ) ' Try changing X1 or X2;' write ( *, * ) ' Perhaps your function F(X) is wrong.' return end if write ( *, * ) ' ' write ( *, * ) ' ' write ( *, * ) 'Regula Falsi' write ( *, * ) ' ' write ( *, * ) ' ' write ( *, * ) 'Step Left Middle Right' iterate = 0 write ( *, * ) ' ' write ( *, '(1x,i4,'' x'',g16.8,16x,g16.8)' ) iterate, xl, xr write ( *, '(1x,4x,''fx'',g16.8,16x,g16.8)' ) fxl, fxr ! ! Initialize weights ! wate1 = fxl wate2 = fxr if ( abs ( fxr ) >= abs ( fxl ) ) then xm = xl fxm = fxl else xm = xr fxm = fxr end if do x_ave = ( abs ( xr ) + abs ( xl ) ) / 2.0E+00 if ( abs ( xr - xl ) <= xatol ) then write ( *, * ) ' ' write ( *, * ) ' Absolute convergence of the X interval.' return end if if ( abs ( xr - xl ) <= xrtol * x_ave ) then write ( *, * ) ' ' write ( *, * ) ' Relative convergence of the X interval.' return end if if ( abs ( fxl ) <= fatol ) then write ( *, * ) ' ' write ( *, * ) ' Absolute convergence of |F(X)|.' exit end if if ( abs ( fxr ) <= fatol ) then write ( *, * ) ' ' write ( *, * ) ' Absolute convergence of |F(X)|.' exit end if iterate = iterate + 1 if ( iterate > itmax ) then write ( *, * ) ' ' write ( *, * ) ' Maximum number of steps taken.' return end if ! ! Compute the next iterate. ! fxmold = fxm if ( wate1 - wate2 == 0.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'The method has broken down.' write ( *, * ) 'A zero divisor was computed.' return end if xm = ( wate1 * xr - wate2 * xl ) / ( wate1 - wate2 ) call f ( xm, fxm, 0 ) write ( *, * ) ' ' write ( *, '(1x,i4,'' x'',3g16.8)' ) iterate, xl, xm, xr write ( *, '(1x,4x,''fx'',3g16.8)' ) fxl, fxm, fxr ! ! Set up the next interval. ! if ( fxl * fxm >= 0.0E+00 ) then xl = xm fxl = fxm wate1 = fxm else xr = xm fxr = fxm wate2 = fxm end if end do return end subroutine falsi2 ( fatol, itmax, x1, x2, xatol, xrtol ) ! !******************************************************************************* ! !! FALSI2 carries out the modified Regula Falsi method. ! ! ! Modified: ! ! 18 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real FATOL, the absolute error tolerance for F(X). ! ! Input, integer ITMAX, the maximum number of steps allowed. ! ! Input, real X1, X2, two points where the function differs in sign. ! ! Input, real XATOL, absolute error tolerance for the root. ! implicit none ! real fatol real fxl real fxm real fxmold real fxr integer iterate integer itmax real wate1 real wate2 real x_ave real x1 real x2 real xatol real xl real xm real xr real xrtol ! ! The starting points ! xl = x1 xr = x2 call f ( xl, fxl, 0 ) call f ( xr, fxr, 0 ) ! ! Check for change of sign. ! if ( fxl * fxr > 0.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'WARNING:' write ( *, * ) ' ' write ( *, * ) ' This method requires a change of sign ' write ( *, * ) ' interval, but your function does NOT change ' write ( *, * ) ' sign over the given interval, so this method ' write ( *, * ) ' cannot be used.' write ( *, * ) ' ' write ( *, * ) 'SUGGESTION:' write ( *, * ) ' ' write ( *, * ) ' Try a different method;' write ( *, * ) ' Try changing X1 or X2;' write ( *, * ) ' Perhaps your function F(X) is wrong.' return end if write ( *, * ) ' ' write ( *, * ) ' ' write ( *, * ) 'Modified Regula Falsi' write ( *, * ) ' ' write ( *, * ) ' ' write ( *, * ) 'Step Left Middle Right' iterate = 0 write ( *, * ) ' ' write ( *, '(1x,i4,'' x'',g16.8,16x,g16.8)' ) iterate, xl, xr write ( *, '(1x,4x,''fx'',g16.8,16x,g16.8)' ) fxl, fxr ! ! Initialize weights ! wate1 = fxl wate2 = fxr if ( abs ( fxr ) >= abs ( fxl ) ) then xm = xl fxm = fxl else xm = xr fxm = fxr end if do if ( abs ( xr - xl ) <= xatol ) then write ( *, * ) ' ' write ( *, * ) ' Absolute convergence of the X interval.' return end if x_ave = ( abs ( xr ) + abs ( xl ) ) / 2.0E+00 if ( abs ( xr - xl ) <= xrtol * x_ave ) then write ( *, * ) ' ' write ( *, * ) ' Relative convergence of the X interval.' return end if if ( abs ( fxl ) <= fatol ) then write ( *, * ) ' ' write ( *, * ) ' Absolute convergence of |F(X)|.' exit end if if ( abs ( fxr ) <= fatol ) then write ( *, * ) ' ' write ( *, * ) ' Absolute convergence of |F(X)|.' exit end if iterate = iterate + 1 if ( iterate > itmax ) then write ( *, * ) ' ' write ( *, * ) ' Maximum number of steps taken.' return end if ! ! Compute the next iterate. ! fxmold = fxm if ( wate1 - wate2 == 0.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'The method has broken down.' write ( *, * ) 'A zero divisor was computed.' return end if xm = ( wate1 * xr - wate2 * xl ) / ( wate1 - wate2 ) call f ( xm, fxm, 0 ) write ( *, * ) ' ' write ( *, '(1x,i4,'' x'',3g16.8)' ) iterate, xl, xm, xr write ( *, '(1x,4x,''fx'',3g16.8)' ) fxl, fxm, fxr ! ! Set up the next interval. ! if ( fxl * fxm >= 0.0E+00 ) then xl = xm fxl = fxm wate1 = fxm if ( fxm * fxmold > 0.0E+00 ) then wate2 = 0.5E+00 * wate2 end if else xr = xm fxr = fxm wate2 = fxm if ( fxm * fxmold > 0.0E+00 ) then wate1 = 0.5E+00 * wate1 end if end if end do return end subroutine fixed ( fatol, itmax, x1, xatol, xrtol ) ! !******************************************************************************* ! !! FIXED implements the fixed point method for a nonlinear equation. ! ! ! THIS NEEDS TO BE REVISED!!! ! ! Modified: ! ! 18 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real FATOL, the absolute error tolerance for F(X). ! ! Input, integer ITMAX, the maximum number of steps allowed. ! ! Input, real X1, a point where the iteration begins. ! ! Input, real XATOL, absolute error tolerance for the root. ! implicit none ! real fatol real fx real gx integer iterate integer itmax real x_ave real x_inc real x1 real xatol real xnew real xold real xrtol ! ! Get the starting point. ! xnew = x1 iterate = 0 call f ( xnew, fx, 0 ) gx = xnew - fx write ( *, * ) ' ' write ( *, * ) ' ' write ( *, * ) 'Fixed Point Iteration' write ( *, * ) ' ' write ( *, * ) ' ' write ( *, * ) 'Iteration x fx gx' write ( *, * ) ' ' write ( *, '(1x,i6,3f20.10)' ) iterate, xnew, fx, gx do iterate = iterate + 1 if ( iterate > itmax ) then write ( *, * ) ' ' write ( *, * ) ' Maximum number of steps taken.' return end if ! ! Compute the next iterate. ! call f ( xnew, gx, 3 ) xold = xnew xnew = gx x_inc = xnew - xold call f ( xnew, fx, 0 ) write ( *, '(1x,i6,3f20.10)' ) iterate, xnew, fx, x_inc ! ! Check for acceptance of function value or interval ! x_ave = ( abs ( xnew ) + abs ( xold ) ) / 2.0E+00 if ( abs ( x_inc ) <= xatol ) then write ( *, * ) ' ' write ( *, * ) ' Absolute convergence of the X increment.' return end if if ( abs ( x_inc ) <= xrtol * x_ave ) then write ( *, * ) ' ' write ( *, * ) ' Relative convergence of the X increment.' return end if if ( abs ( fx ) <= fatol ) then write ( *, * ) ' ' write ( *, * ) ' Absolute convergence of |F(X)|.' exit end if end do return end subroutine funscl ( angle_to_radian, arg1, arg2, ierror, result, rtol, sym ) ! !******************************************************************************* ! !! FUNSCL evaluates a scalar function of one or more scalar arguments. ! ! ! Modified: ! ! 15 April 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ARG2, ARG2, the values of the arguments of the function. ! ! Output, integer IERROR, error flag. ! 0, no error occurred. ! nonzero, an error occurred. ! ! Output, real RESULT, the value of the scalar function. ! ! Input, real RTOL, the rounding tolerance. ! ! Input, character ( len = * ) SYM, the symbolic name of the ! function to be evaluated. ! implicit none ! real, parameter :: degrad = 3.14159265358979 / 180.0E+00 ! real angle_to_radian real arg1 real arg2 integer i_lcm integer iarg integer iarg1 integer iarg2 integer ierror integer i_gcd real result real rtol logical s_eqi real sarg character ( len = * ) sym real temp ! ierror = 0 result = 0.0E+00 ! ! Addition. ! if ( sym == '+' ) then result = arg1 + arg2 ! ! Subtraction. ! else if ( sym == '-' ) then result = arg1 - arg2 ! ! Division. ! else if ( sym == '/' ) then if ( arg2 == 0.0E+00 ) then ierror = 1 write ( *, * ) ' ' write ( *, * ) 'FUNSCL - Error!' write ( *, * ) ' Attempt to divide by 0!' return end if result = arg1 / arg2 ! ! Multiplication. ! else if ( sym == '*' ) then result = arg1 * arg2 ! ! Exponentiation. ! else if ( sym == '**' .or. sym == '^' ) then if ( arg1 == 0.0E+00 .and. arg2 == 0.0 ) then write ( *, * ) ' ' write ( *, * ) 'FUNSCL - Error!' write ( *, * ) ' Attempt to compute 0**0 !' ierror = 1 else if ( arg1 < 0.0E+00 .and. real(int(arg2)) /= arg2 ) then write ( *, * ) ' ' write ( *, * ) 'FUNSCL - Error!' write ( *, * ) ' Illegal exponentiation:' write ( *, * ) ' ', arg1, ' ** ', arg2 ierror = 1 else sarg = 1.0E+00 if ( arg1 < 0.0E+00 ) then iarg = int ( arg2 ) if ( 2 * ( iarg / 2 ) /= iarg ) then sarg = - 1.0E+00 end if end if result = sarg * abs ( arg1 )**arg2 end if ! ! Assignment. ! else if ( sym == '=' ) then result = arg2 ! ! Absolute value. ! else if ( s_eqi ( sym, 'ABS' ) ) then result = abs ( arg1 ) ! ! Arc Cosine. ! else if ( s_eqi ( sym, 'ACOS' ) ) then if ( arg1 < -1.0E+00 .or. arg1 > 1.0 ) then write ( *, * ) ' ' write ( *, * ) 'FUNSCL - Error!' write ( *, * ) ' Illegal inverse cosine of ', arg1 ierror = 1 return end if result = acos ( arg1 ) / angle_to_radian ! ! Arc Sine. ! else if ( s_eqi ( sym, 'ASIN' ) ) then if ( arg1 < - 1.0E+00 .or. arg1 > 1.0 ) then write ( *, * ) ' ' write ( *, * ) 'FUNSCL - Error!' write ( *, * ) ' Illegal inverse sine of ', arg1 ierror = 1 end if result = asin ( arg1 ) / angle_to_radian ! ! Arc Tangent. ! else if ( s_eqi ( sym, 'ATAN' ) ) then result = atan ( arg1 ) / angle_to_radian ! ! Arc Tangent of ratio. ! else if ( s_eqi ( sym, 'ATAN2' ) ) then if ( arg1 == 0.0E+00 .and. arg2 == 0.0 ) then result = 0.0E+00 else result = atan2 ( arg1, arg2 ) / angle_to_radian end if ! ! Cosine. ! else if ( s_eqi ( sym, 'COS' ) ) then result = cos ( angle_to_radian * arg1 ) ! ! Cotangent of angle. ! else if ( s_eqi ( sym, 'COT' ) ) then result = cos ( angle_to_radian * arg1 ) / sin ( angle_to_radian * arg1 ) ! ! Hyperbolic cosine. ! else if ( s_eqi ( sym, 'COSH' ) ) then result = cosh ( angle_to_radian * arg1 ) ! ! Cosecant of angle. ! else if ( s_eqi ( sym, 'CSC' ) ) then result = 1.0E+00 / sin ( angle_to_radian * arg1 ) ! ! Determinant (of a scalar). ! else if ( s_eqi ( sym, 'DET' ) ) then result = arg1 ! ! Diagonal (of a scalar). ! else if ( s_eqi ( sym, 'DIAG' ) ) then result = arg1 ! ! Exponential. ! else if ( s_eqi ( sym, 'EXP' ) ) then result = exp ( arg1 ) ! ! Greatest common factor. ! else if ( s_eqi ( sym, 'GCD' ) ) then iarg1 = nint ( arg1 ) iarg2 = nint ( arg2 ) result = real ( i_gcd ( iarg1, iarg2 ) ) ! ! Hilbert. ! else if ( s_eqi ( sym, 'HILBERT' ) ) then result = 1.0E+00 ! ! Hilbert inverse. ! else if ( s_eqi ( sym, 'HILBINV' ) ) then result = 1.0E+00 ! ! Identity. ! else if ( s_eqi ( sym, 'ID' ) ) then result = 1.0E+00 ! ! Integer value. ! else if ( s_eqi ( sym, 'INT' ) ) then result = real ( int ( arg1 ) ) ! ! Inverse (of a scalar). ! else if ( s_eqi ( sym, 'INV' ) ) then if ( arg1 == 0.0E+00 ) then ierror = 1 write ( *, * ) ' ' write ( *, * ) 'FUNSCL - Error!' write ( *, * ) ' Attempt to compute 1/0!' return end if result = 1.0E+00 / arg1 ! ! Least common multiple. ! else if ( s_eqi ( sym, 'LCM' ) ) then iarg1 = nint ( arg1 ) iarg2 = nint ( arg2 ) result = real ( i_lcm ( iarg1, iarg2 ) ) ! ! Natural logarithm. ! else if ( s_eqi ( sym, 'ALOG' ) .or. s_eqi ( sym, 'LN' ) .or. & s_eqi ( sym, 'LOG' ) ) then if ( arg1 <= 0.0E+00 ) then ierror = 1 write ( *, * ) ' ' write ( *, * ) 'FUNSCL - Error!' write ( *, * ) ' Illegal LOG of ', arg1 return end if result = log ( arg1 ) ! ! Logarithm base 10. ! else if ( s_eqi ( sym, 'ALOG10' ) .or. s_eqi ( sym, 'LOG10' ) ) then if ( arg1 <= 0.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'FUNSCL - Error!' write ( *, * ) ' Illegal LOG10 of ', arg1 ierror = 1 return end if result = log10 ( arg1 ) ! ! Logarithm base 2. ! else if ( s_eqi ( sym, 'LOG2' ) ) then if ( arg1 <= 0.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'FUNSCL - Error!' write ( *, * ) ' Illegal LOG2 of ', arg1 ierror = 1 return end if result = log ( arg1 ) / log ( 2.0E+00 ) ! ! Maximum of two values. ! else if ( s_eqi ( sym, 'MAX' ) ) then if ( arg1 > arg2 ) then result = arg1 else result = arg2 end if ! ! Minimum of two values. ! else if ( s_eqi ( sym, 'MIN' ) ) then if ( arg1 < arg2 ) then result = arg1 else result = arg2 end if ! ! Negation ! else if ( s_eqi ( sym, 'NEG' ) ) then result = - arg1 ! ! Nearest integer value ! else if ( s_eqi ( sym, 'NINT' ) ) then result = anint ( arg1 ) ! ! NORM0, NORM1, NORM2, NORME or NORMF of a scalar. ! else if ( s_eqi ( sym, 'NORM0' ) .or. s_eqi ( sym, 'NORM1' ) .or. & s_eqi ( sym, 'NORM2' ) .or. s_eqi ( sym, 'NORME' ) .or. & s_eqi ( sym, 'NORMF' ) ) then result = abs ( arg1 ) ! ! POLY (scalar coefficient array, means constant polynomial). ! else if ( s_eqi ( sym, 'POLY' ) ) then result = arg1 ! ! Random value. ! else if ( s_eqi ( sym, 'RAN' ) ) then call random_number ( result ) ! ! Rounding. ! else if ( s_eqi ( sym, 'ROUND' ) ) then if ( abs ( arg1 ) < rtol ) then result = 0.0E+00 else result = arg1 end if ! ! Secant. ! else if ( s_eqi ( sym, 'SEC' ) ) then result = 1.0E+00 / cos ( angle_to_radian * arg1 ) ! ! Sine. ! else if ( s_eqi ( sym, 'SIN' ) ) then result = sin ( angle_to_radian * arg1 ) ! ! Hyperbolic sine. ! else if ( s_eqi ( sym, 'SINH' ) ) then result = sinh ( angle_to_radian * arg1 ) ! ! Square root. ! else if ( s_eqi ( sym, 'SQRT' ) ) then if ( arg1 < 0.0E+00 ) then ierror = 1 write ( *, * ) ' ' write ( *, * ) 'FUNSCL - Error!' write ( *, * ) ' Illegal SQRT of ', arg1 return end if result = sqrt ( arg1 ) ! ! Step function. ! else if ( s_eqi ( sym, 'STEP' ) ) then if ( arg1 < 0.0E+00 ) then result = 0.0E+00 else result = 1.0E+00 end if ! ! Tangent. ! else if ( s_eqi ( sym, 'TAN' ) ) then result = tan ( angle_to_radian * arg1 ) ! ! Hyperbolic tangent. ! else if ( s_eqi ( sym, 'TANH' ) ) then result = tanh ( angle_to_radian * arg1 ) ! ! Trace (of a scalar) ! else if ( s_eqi ( sym, 'TRACE' ) ) then result = arg1 ! ! Transpose (of a scalar) ! else if ( s_eqi ( sym, 'TRANS' ) ) then result = arg1 ! ! Zero ! else if ( s_eqi ( sym, 'ZEROS' ) ) then result = 0.0E+00 ! ! Factorial function. ! else if ( sym == '!' ) then if ( arg1 < 0.0E+00 .or. arg1 > 25.0 ) then write ( *, * ) ' ' write ( *, * ) 'FUNSCL - Error!' write ( *, * ) ' Argument out of range for factorial!' ierror = 1 return end if temp = arg1 result = 1.0E+00 do while ( temp > 1.0E+00 ) result = result * temp temp = temp - 1.0E+00 end do ! ! Unknown function. ! else result = 0.0E+00 ierror = 1 write ( *, * ) ' ' write ( *, * ) 'FUNSCL - Error!' write ( *, * ) ' Unrecognized function name:', trim ( sym ) end if return end subroutine funval ( iarg1, iarg2, ierror, ifree, iopsym, iprsym, & ipset, ipval, irad, irtol, itemp, maxsym, maxval, nsyms, & sym, symbol, valsym ) ! !******************************************************************************* ! !! FUNVAL evaluates a scalar, vector or matrix valued function. ! ! ! Modified: ! ! 17 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer IARG1, IARG2, are the indices ! of the arguments to the current function. ! ! Output, integer IERROR. ! ! IERROR is an error flag. If it is zero on return, ! then no error was detected, and the formula was ! evaluated successfully. If it is nonzero on return, ! then an error was found, and the evaluation of the ! formula could not be carried out. ! ! Input, integer IFREE, the address of the next free ! memory location in VALSYM. ! ! Input, integer IOPSYM(MAXSYM), contains, for each ! symbol, the number of operands. In particular, if ! a symbol represents a constant, IOPSYM(I) is 0. ! If a symbol represents a unary operator such as ABS, ! IOPSYM(I) is 1. IOPSYM may be as large as 3. ! ! Input, integer IPRSYM(MAXSYM), the priorities of the functions. ! ! Output, integer IPSET. If the function was an 'INDX1' or ! 'INDX2', then IPSET is the relative offset of the value. ! ! Input, integer IPVAL(MAXSYM), contains, for each symbol, ! the address in VALSYM where the value of the symbol ! is stored, if it is a scalar. If the symbol represents ! a vector or matrix, IPVAL(IARG) points to location of ! the first entry. ! ! Input, integer ITEMP, the index of the constant. ! ! Input, integer MAXIW, the amount of space in IWORK, which ! should be at least MAXROW. ! ! Input, integer MAXSYM, the maximum number of symbols allowed. ! ! Input, integer MAXVAL, the maximum number of values allowed ! in VALSYM. ! ! Input, integer NSYMS, the number of declared symbols. ! ! Input, integer NUMDIM(2,MAXSYM). ! For each symbol I, NUMDIM(1,I) is the number of rows in its ! value, and NUMDIM(2,I) is the number of columns. ! ! Input, character ( len = MAXLEN ) SYM, the symbolic name of the function. ! ! Input, character ( len = MAXLEN ) SYMBOL(MAXSYM), the symbolic names ! of the internally defined constants and functions. ! ! Input, real VALSYM(MAXVAL), contains the values of all the ! symbolic variables. ! implicit none ! integer, parameter :: maxlen = 20 integer maxsym integer maxval ! real angle_to_radian real arg1 real arg2 character ( len = 3 ) ctemp integer i integer iarg1 integer iarg2 integer ierror integer ifree integer index1 integer index4 integer indx integer info integer iopsym(maxsym) integer iprsym(maxsym) integer ipset integer ipval(maxsym) integer irad integer irtol integer itemp integer jseed integer nsyms real result real rtol logical s_eqi character ( len = maxlen ) sym character ( len = maxlen ) symbol(maxsym) real valsym(maxval) ! ipset = 0 ierror = 0 ! ! If the operator is assignment, then the user may be trying to ! reset the values of certain reserved variables, including ! PI and EPS. ! if ( sym == '=' ) then if ( s_eqi ( symbol(iarg1), 'E' ) .or. & s_eqi ( symbol(iarg1), 'EPS' ) .or. & s_eqi ( symbol(iarg1), 'PI' ) ) then ierror = 1 write ( *, * ) ' ' write ( *, * ) 'FUNVAL - Error!' write ( *, * ) ' You may not change the value of ', trim ( symbol(iarg1) ) return end if end if index1 = ipval(iarg1) if ( nsyms >= maxsym ) then ierror = 1 write ( *, * ) ' ' write ( *, * ) 'FUNVAL - Error!' write ( *, * ) ' Not enough free memory left!' write ( *, * ) ' The KILL or INIT commands may help!' return end if nsyms = nsyms + 1 iopsym(nsyms) = 0 iprsym(nsyms) = 10 itemp = itemp + 1 call i_to_s_zero ( itemp, ctemp ) symbol(nsyms) = 'STK000' symbol(nsyms)(4:6) = ctemp ipval(nsyms) = ifree index4 = ipval(nsyms) arg1 = valsym(index1) if ( iarg2 == 0 ) then arg2 = 0.0E+00 else arg2 = valsym(ipval(iarg2)) end if angle_to_radian = valsym(irad) rtol = valsym(irtol) call funscl ( angle_to_radian, arg1, arg2, ierror, result, rtol, sym ) valsym(index4) = result if ( sym == '=' ) then arg1 = arg2 valsym(index1) = arg2 end if return end subroutine halley ( fatol, itmax, x1, xatol, xrtol ) ! !******************************************************************************* ! !! HALLEY implements Halley's method for a nonlinear equation. ! ! ! Modified: ! ! 18 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real FATOL, the absolute error tolerance for F(X). ! ! Input, integer ITMAX, the maximum number of steps allowed. ! ! Input, real X1, a point where the iteration begins. ! ! Input, real XATOL, absolute error tolerance for the root. ! implicit none ! real bot real fatol real fprime real ftwo real fx integer iterate integer itmax real x_ave real x_inc real x1 real xatol real xnew real xold real xrtol ! write ( *, * ) ' ' write ( *, * ) ' ' write ( *, * ) 'Halley''s method' write ( *, * ) ' ' write ( *, * ) ' ' write ( *, * ) 'Iteration X FX' write ( *, * ) ' ' iterate = 0 xnew = x1 call f ( xnew, fx, 0 ) write ( *, '(i6,f20.10,f20.10)' ) iterate, xnew, fx if ( abs ( fx ) <= fatol ) then write ( *, * ) ' ' write ( *, * ) ' Absolute convergence of |F(X)|.' return end if do iterate = iterate + 1 if ( iterate > itmax ) then write ( *, * ) ' ' write ( *, * ) ' Maximum number of steps taken.' return end if xold = xnew call f ( xold, fx, 0 ) call f ( xold, fprime, 1 ) call f ( xold, ftwo, 2 ) bot = fprime**2 - 0.5E+00 * ftwo * fx if ( bot == 0.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'The method has broken down.' write ( *, * ) 'A zero divisor was computed.' exit end if if ( fprime == 0.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'The method has broken down.' write ( *, * ) 'F''(X) = 0.' exit end if ! ! Compute the next iterate. ! x_inc = - ( fx * fprime ) / bot xnew = xold + x_inc call f ( xnew, fx, 0 ) write ( *, '(i6,f20.10,f20.10)' ) iterate, xnew, fx x_ave = ( abs ( xnew ) + abs ( xold ) ) / 2.0E+00 if ( abs ( x_inc ) <= xatol ) then write ( *, * ) ' ' write ( *, * ) ' Absolute convergence of the X increment.' exit end if if ( abs ( x_inc ) <= xrtol * x_ave ) then write ( *, * ) ' ' write ( *, * ) ' Relative convergence of the X increment.' exit end if if ( abs ( fx ) <= fatol ) then write ( *, * ) ' ' write ( *, * ) ' Absolute convergence of |F(X)|.' exit end if end do return end subroutine help ! !******************************************************************************* ! !! HELP prints out a list of commands. ! ! ! Modified: ! ! 19 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! None ! implicit none ! write ( *, * ) ' ' write ( *, * ) 'To SET something:' write ( *, * ) ' ' write ( *, * ) ' F = set F(X).' write ( *, * ) ' F'' = set F''(X).' write ( *, * ) ' F" = set F"(X).' write ( *, * ) ' FATOL = function absolute tolerance;' write ( *, * ) ' G = set G(X), the fixed point function;' write ( *, * ) ' ITMAX = maximum number of iterations;' write ( *, * ) ' M set the method' write ( *, * ) ' X1 = First starting point;' write ( *, * ) ' X2 = Second starting point;' write ( *, * ) ' X3 = Third starting point.' write ( *, * ) ' XATOL = X absolute tolerance;' write ( *, * ) ' XRTOL = X relative tolerance;' write ( *, * ) ' ' write ( *, * ) 'To DO something:' write ( *, * ) ' ' write ( *, * ) ' GO Go for the root;' write ( *, * ) ' H Help! (print this list);' write ( *, * ) ' P Print the current values;' write ( *, * ) ' Q Quit.' return end function i_gcd ( i, j ) ! !******************************************************************************* ! !! I_GCD finds the greatest common divisor of I and J. ! ! ! Modified: ! ! 03 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer I, J, two numbers whose greatest common divisor ! is desired. ! ! Output, integer I_GCD, the greatest common divisor of I and J. ! ! Note that only the absolute values of I and J are ! considered, so that the result is always nonnegative. ! ! If I or J is 0, I_GCD is returned as max ( 1, abs ( I ), abs ( J ) ). ! ! If I and J have no common factor, I_GCD is returned as 1. ! ! Otherwise, using the Euclidean algorithm, I_GCD is the ! largest common factor of I and J. ! implicit none ! integer i integer i_gcd integer ip integer iq integer ir integer j ! i_gcd = 1 ! ! Return immediately if either I or J is zero. ! if ( i == 0 ) then i_gcd = max ( 1, abs ( j ) ) return else if ( j == 0 ) then i_gcd = max ( 1, abs ( i ) ) return end if ! ! Set IP to the larger of I and J, IQ to the smaller. ! This way, we can alter IP and IQ as we go. ! ip = max ( abs ( i ), abs ( j ) ) iq = min ( abs ( i ), abs ( j ) ) ! ! Carry out the Euclidean algorithm. ! do ir = mod ( ip, iq ) if ( ir == 0 ) then exit end if ip = iq iq = ir end do i_gcd = iq return end function i_lcm ( i, j ) ! !******************************************************************************* ! !! I_LCM computes the least common multiple of two integers. ! ! ! Definition: ! ! The least common multiple may be defined as ! ! LCM(I,J) = ABS( I * J ) / GCF(I,J) ! ! where GCF(I,J) is the greatest common factor of I and J. ! ! Modified: ! ! 01 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer I, J, the integers whose I_LCM is desired. ! ! Output, integer I_LCM, the least common multiple of I and J. ! I_LCM is never negative. I_LCM is 0 if either I or J is zero. ! implicit none ! integer i integer i_gcd integer j integer i_lcm ! i_lcm = abs ( i * ( j / i_gcd ( i, j ) ) ) return end subroutine i_to_s_zero ( intval, s ) ! !******************************************************************************* ! !! I_TO_S_ZERO converts an integer to a string, with zero padding. ! ! ! Examples: ! ! Assume that S is 6 characters long: ! ! INTVAL S ! ! 1 000001 ! -1 -00001 ! 0 000000 ! 1952 001952 ! 123456 123456 ! 1234567 ****** <-- Not enough room! ! ! Modified: ! ! 04 August 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer INTVAL, an integer to be converted. ! ! Output, character ( len = * ) S, the representation of the integer. ! The integer will be right justified, and zero padded. ! If there is not enough space, the string will be filled with stars. ! implicit none ! character c integer i integer idig integer ihi integer ilo integer intval integer ipos integer ival character ( len = * ) s ! s = ' ' ilo = 1 ihi = len ( s ) if ( ihi <= 0 ) then return end if ! ! Make a copy of the integer. ! ival = intval ! ! Handle the negative sign. ! if ( ival < 0 ) then if ( ihi <= 1 ) then s(1:1) = '*' return end if ival = - ival s(1:1) = '-' ilo = 2 end if ! ! Working from right to left, strip off the digits of the integer ! and place them into S(ILO:IHI). ! ipos = ihi do while ( ival /= 0 .or. ipos == ihi ) idig = mod ( ival, 10 ) ival = ival / 10 if ( ipos < ilo ) then do i = 1, ihi s(i:i) = '*' end do return end if call digit_to_ch ( idig, c ) s(ipos:ipos) = c ipos = ipos - 1 end do ! ! Fill the empties with zeroes. ! do i = ilo, ipos s(i:i) = '0' end do return end subroutine inicom ( ifinis, ifree, indx1, indx2, ineg, infix, intsym, iopsym, & iprsym, ipval, irad, irpn, irtol, istack, maxrpn, maxsym, & maxval, nints, nsym, nsymp, nsyms, symbol, valsym, value ) ! !******************************************************************************* ! !! INICOM initializes data for COMRPN. ! ! ! Modified: ! ! 29 October 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IFINIS, the index in SYMBOL of the symbol ! '$', meaning the end of the formula. ! ! Output, integer IFREE, the index of the first free address in VALSYM. ! ! Output, integer INDX1, the index of the symbol 'INDEX1'. ! ! Output, integer INDX2, the index of the symbol 'INDEX2'. ! ! Output, integer INEG, the index of the symbol 'NEG'. ! ! Output, character ( len = * ) INFIX, space for an infix formula. ! ! Output, integer INTSYM(80), a set of integers which ! are the indices of symbols, representing an infix formula. ! ! Output, integer IPVAL(MAXSYM), contains, for each symbol, ! the address in VALSYM where the value of the symbol ! is stored, if it is a scalar. If the symbol represents ! a vector or matrix, IPVAL(IARG) points to location of ! the first entry. ! ! Output, integer IRAD, the index of the symbol 'ANGLE_TO_RADIAN'. ! ! Output, integer IRPN(MAXRPN), used to store the compiled ! versions of user formulas. ! ! Workspace, integer ISTACK(MAXSYM), workspace for interpreting ! the formula. ! ! Input, integer MAXRPN, specifies the length of IRPN. ! ! Input, integer MAXSYM, the maximum number of symbols allowed. ! ! Input, integer MAXVAL, the maximum number of values allowed in VALSYM. ! ! Output, integer NINTS, the number of useful entries in INTSYM. ! ! Output, integer NSYM, the total number of symbols, including temporaries. ! ! Output, integer NSYMP, the number of permanent symbols. ! ! Output, integer NYMS, the number of declared symbols. ! ! Output, integer NUMDIM(2,MAXSYM). ! For each symbol I, NUMDIM(1,I) is the number of rows in its ! value, and NUMDIM(2,I) is the number of columns. ! ! Output, character ( len = MAXLEN ) SYMBOL(MAXSYM), the symbolic names ! of the internally defined constants and functions. ! ! Output, real VALSYM(MAXVAL), contains the values of all the ! symbolic variables. ! ! Workspace, real VALUE, space used to hold ! the value of the variable to be saved. ! implicit none ! integer, parameter :: maxlen = 20 integer maxrpn integer maxsym integer maxval ! integer i integer ifinis integer ifree integer indx1 integer indx2 integer ineg character ( len = * ) infix integer intsym(80) integer iopsym(maxsym) integer iprsym(maxsym) integer ipval(maxsym) integer irad integer irpn(maxrpn) integer irtol integer istack(maxsym) integer j integer nints integer nsym integer nsymp integer nsyms logical s_eqi character ( len = maxlen ) symbol(maxsym) real temp real valsym(maxval) real value ! nsym = 0 nsym = nsym + 1 symbol(nsym) = 'ABS' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'ACOS' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'ALOG' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'ALOG10' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'ANGLE_TO_RADIAN' iopsym(nsym) = 0 iprsym(nsym) = 10 valsym(nsym) = 1.0E+00 nsym = nsym + 1 symbol(nsym) = 'ASIN' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'ATAN' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'ATAN2' iopsym(nsym) = 2 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'COS' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'COSH' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'COT' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'CSC' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'DET' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'DIAG' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'E' iopsym(nsym) = 0 iprsym(nsym) = 10 valsym(nsym) = exp ( 1.0E+00 ) nsym = nsym + 1 symbol(nsym) = 'EPS' iopsym(nsym) = 0 iprsym(nsym) = 10 valsym(nsym) = epsilon ( 1.0E+00 ) nsym = nsym + 1 symbol(nsym) = 'EVAL' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'EVEC' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'EXP' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'GCD' iopsym(nsym) = 2 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'HILBERT' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'HILBINV' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'HOUSE' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'ID' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'INT' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'INV' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'IVAL' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'KRYLOV' iopsym(nsym) = 3 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'LCM' iopsym(nsym) = 2 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'LN' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'LOG' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'LOG10' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'LOG2' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'MAX' iopsym(nsym) = 2 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'MIN' iopsym(nsym) = 2 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'NEG' iopsym(nsym) = 1 iprsym(nsym) = 5 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'NINT' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'NORM0' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'NORM1' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'NORM2' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'NORME' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'NORMF' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'NORMS' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'PI' iopsym(nsym) = 0 iprsym(nsym) = 10 valsym(nsym) = 4.0E+00 * atan2 ( 1.0, 1.0 ) nsym = nsym + 1 symbol(nsym) = 'POLY' iopsym(nsym) = 2 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'RAN' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'RCOND' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'ROUND' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'RTOL' iopsym(nsym) = 0 iprsym(nsym) = 10 valsym(nsym) = epsilon ( 1.0E+00 ) nsym = nsym + 1 symbol(nsym) = 'RVAL' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'SEC' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'SEED' iopsym(nsym) = 0 iprsym(nsym) = 10 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'SIN' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'SINH' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'SQRT' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'STEP' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'TAN' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'TANH' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'TRACE' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'TRANS' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'ZEROS' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = '**' iopsym(nsym) = 2 iprsym(nsym) = 8 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = '^' iopsym(nsym) = 2 iprsym(nsym) = 8 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = '+' iopsym(nsym) = 2 iprsym(nsym) = 5 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = '-' iopsym(nsym) = 2 iprsym(nsym) = 5 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = '/' iopsym(nsym) = 2 iprsym(nsym) = 7 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = '*' iopsym(nsym) = 2 iprsym(nsym) = 6 valsym(nsym) = 0.0E+00 ! ! Take care of fact that ninny UNIX FORTRAN compilers won't let us ! explicitly use a backslash character! ! nsym = nsym + 1 symbol(nsym) = char ( 92 ) iopsym(nsym) = 2 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = '(' iopsym(nsym) = -1 iprsym(nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = ')' iopsym(nsym) = -1 iprsym(nsym) = 2 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = '=' iopsym(nsym) = 2 iprsym(nsym) = 3 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = ',' iopsym(nsym) = -1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = '!' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'INDEX1' iopsym(nsym) = 2 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'INDEX2' iopsym(nsym) = 3 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = '$' iopsym(nsym) = -1 iprsym(nsym) = 0 valsym(nsym) = 0.0E+00 nsymp = nsym nsyms = nsymp ifree = nsymp + 1 infix = ' ' intsym(1:80) = 0 irpn(1:maxrpn) = 0 istack(1:maxsym) = 0 nints = 0 do i = 1, nsymp ipval(i) = i end do do i = 1, nsymp if ( s_eqi ( symbol(i), 'ANGLE_TO_RADIAN' ) ) then irad = i else if ( s_eqi ( symbol(i), 'NEG' ) ) then ineg = i else if ( symbol(i) == '$' ) then ifinis = i else if ( s_eqi ( symbol(i), 'INDEX1' ) ) then indx1 = i else if ( s_eqi ( symbol(i), 'INDEX2' ) ) then indx2 = i else if ( s_eqi ( symbol(i), 'RTOL' ) ) then irtol = i end if end do call random_seed value = 0.0E+00 return end subroutine method_choice ( method, mform, mstart, names, nform, maxmet, & nstart, ierror ) ! !******************************************************************************* ! !! METHOD_CHOICE allows the user to choose the method. ! ! ! Modified: ! ! 19 November 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer METHOD, the index of the method chosen by the user. ! ! Output, integer MFORM, the number of formulas (f(x), f'(x), ...) ! needed for this method. ! ! Output, integer MSTART, the number of starting points needed for ! the method chosen. ! ! Input, character ( len = 35 ) NAMES(MAXMET), the names of the methods. ! ! Input, integer NFORM(MAXMET), the number of formulas (f(x), f'(x), ...) ! needed for each method. ! ! Input, integer MAXMET, the number of methods available. ! ! Input, integer NSTART(MAXMET), the number of starting points needed ! for each method. ! ! Output, integer IERROR, error flag. ! implicit none ! integer maxmet ! integer i integer ierror integer ios integer method integer mform integer mstart character ( len = 35 ) names(maxmet) integer nform(maxmet) integer nstart(maxmet) ! ierror = 0 write ( *, * ) ' ' write ( *, * ) 'The available methods are:' write ( *, * ) ' ' do i = 1, maxmet write ( *, '(i4,2x,a)' ) i, names(i) end do write ( *, * ) ' ' write ( *, * ) 'Enter the number of the method.' read ( *, *, iostat = ios ) method if ( ios /= 0 ) then ierror = 1 write ( *, * ) ' ' write ( *, * ) 'METHOD_CHOICE - Warning!' write ( *, * ) ' Your input was not acceptable.' return end if if ( method < 1 .or. method > maxmet ) then write ( *, * ) ' ' write ( *, * ) 'METHOD_CHOICE - Warning!' write ( *, * ) ' The method number you chose was unacceptable.' write ( *, * ) ' The default method was chosen.' method = 1 end if mform = nform(method) mstart = nstart(method) return end subroutine laguerre ( fatol, itmax, x1, xatol, xrtol ) ! !******************************************************************************* ! !! LAGUERRE implements Laguerre's algorithm. ! ! ! Modified: ! ! 18 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real FATOL, the absolute error tolerance for F(X). ! ! Input, integer ITMAX, the maximum number of steps allowed. ! ! Input, real X1, two points where the function differs in sign. ! ! Input, real XATOL, absolute error tolerance for the root. ! implicit none ! real degree real denom real f0 real f1 real f2 real fatol real fx real hx real hxn integer iterate integer itmax real x_ave real x_inc real x1 real xatol real xnew real xold real xrtol ! iterate = 0 ! ! DEGREE is the polynomial degree, the highest exponent on X. ! write ( *, * ) ' ' write ( *, * ) 'Laguerre''s method for polynomials.' write ( *, * ) 'Enter the degree of the polynomial:' read ( *, * ) degree xnew = x1 xold = x1 call f ( xnew, fx, 0 ) write ( *, * ) ' ' write ( *, * ) ' ' write ( *, * ) 'Laguerre''s method for polynomials' write ( *, * ) ' ' write ( *, * ) ' ' write ( *, * ) 'Iteration x fx' write ( *, * ) ' ' write ( *, '(1x,i6,f20.10,f20.10)' ) iterate, xnew, fx do iterate = iterate + 1 if ( iterate > itmax ) then write ( *, * ) ' ' write ( *, * ) ' Maximum number of steps taken.' return end if call f ( xnew, f0, 0 ) call f ( xnew, f1, 1 ) call f ( xnew, f2, 2 ) hx = ( degree - 1.0E+00 ) * & ( ( degree - 1.0E+00 ) * f1**2 - degree * f0 * f2 ) ! ! Check for attempt to take negative square root. ! if ( hx < 0.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'The algorithm has broken down.' write ( *, * ) 'Square root of a negative value.' return end if ! ! Determine the sign of hx in denominator of iterative formula. ! if ( f1 >= 0.0E+00 ) then hxn = sqrt ( hx ) else hxn = - sqrt ( hx ) end if denom = f1 + hxn ! ! Check for division by zero ! if ( denom == 0.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'The algorithm breaks down.' write ( *, * ) 'A zero divisor was computed.' return end if x_inc = - ( degree * f0 ) / denom xold = xnew xnew = xold + x_inc call f ( xnew, fx, 0 ) write ( *, '(1x,i6,f20.10,f20.10)' ) iterate, xnew, fx x_ave = ( abs ( xnew ) + abs ( xold ) ) / 2.0E+00 if ( abs ( x_inc ) < xatol * ( 1.0E+00 + x_ave ) ) then write ( *, * ) ' ' write ( *, * ) 'Convergence of the X increment.' return end if if ( abs ( fx ) <= fatol ) then write ( *, * ) ' ' write ( *, * ) ' Absolute convergence of |F(X)|.' return end if end do return end subroutine muller_c ( fatol, itmax, x1, x2, x3, xatol, xrtol ) ! !******************************************************************************* ! !! MULLER_C carries out Muller's method for complex values. ! ! ! Discussion: ! ! This can't work until your compiler can handle complex data. ! ! Modified: ! ! 18 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real FATOL, the absolute error tolerance for F(X). ! ! Input, integer ITMAX, the maximum number of steps allowed. ! ! Input, complex X1, X2, X3, three points to start the iteration. ! ! Input, real XATOL, absolute error tolerance for the root. ! implicit none ! double precision a double precision b double precision c double precision discrm real fatol complex fminus complex fplus complex fxmid complex fxnew complex fxold integer iterate integer itmax real x_ave real x_inc complex x1 complex x2 complex x3 real xatol complex xlast complex xmid complex xminus complex xnew complex xold complex xplus real xrtol ! xnew = x1 xmid = x2 xold = x3 call fc ( xnew, fxnew, 0 ) call fc ( xmid, fxmid, 0 ) call fc ( xold, fxold, 0 ) write ( *, * ) ' ' write ( *, * ) ' ' write ( *, * ) 'Muller''s method (complex root version)' write ( *, * ) ' ' write ( *, * ) ' ' write ( *, * ) 'Iteration x fx disc' write ( *, * ) ' ' iterate = 0 write ( *, '(1x,i6,f20.10,f20.10,f20.10)' ) iterate, xnew, fxnew if ( abs ( fxnew ) < fatol ) then write ( *, * ) ' ' write ( *, * ) '|F(X)| is below the tolerance.' return end if do if ( abs ( fxnew ) >= abs ( fxmid ) ) then call c_swap ( xnew, xmid ) call c_swap ( fxnew, fxmid ) end if xlast = xnew iterate = iterate + 1 if ( iterate > itmax ) then write ( *, * ) ' ' write ( *, * ) ' Maximum number of steps taken.' exit end if a = dble ( ( xmid - xnew ) * ( fxold - fxnew ) & - ( xold - xnew ) * ( fxmid - fxnew ) ) b = dble ( ( xold - xnew )**2 * ( fxmid - fxnew ) & - ( xmid - xnew )**2 * ( fxold - fxnew ) ) c = dble ( ( xold - xnew ) * ( xmid - xnew ) * ( xold - xmid ) * fxnew ) xold = xmid xmid = xnew ! ! Apply the quadratic formula to get roots XPLUS and XMINUS. ! discrm = b**2 - 4.0E+00 * a * c if ( a == 0.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'The algorithm has broken down.' write ( *, * ) 'The quadratic coefficient A is zero.' exit end if xplus = xnew + sngl ( ( - b + sqrt ( discrm ) ) / ( 2.0E+00 * a ) ) call fc ( xplus, fplus, 0 ) xminus = xnew + sngl ( ( - b - sqrt ( discrm ) ) / ( 2.0E+00 * a ) ) call fc ( xminus, fminus, 0 ) ! ! Take whichever of the two quadratic roots is closest to a root of the function. ! if ( abs ( fminus ) < abs ( fplus ) ) then xnew = xminus else xnew = xplus end if fxold = fxmid fxmid = fxnew call fc ( xnew, fxnew, 0 ) write ( *, '(1x,i6,f20.10,f20.10,f20.10)' ) iterate, xnew, fxnew, discrm ! ! Check for convergence. ! x_ave = ( abs ( xnew ) + abs ( xmid ) + abs ( xold ) ) / 3.0E+00 x_inc = xnew - xmid if ( abs ( x_inc ) <= xatol ) then write ( *, * ) ' ' write ( *, * ) ' Absolute convergence of the X increment.' exit end if if ( abs ( x_inc ) <= xrtol * x_ave ) then write ( *, * ) ' ' write ( *, * ) ' Relative convergence of the X increment.' exit end if if ( abs ( fxnew ) <= fatol ) then write ( *, * ) ' ' write ( *, * ) ' Absolute convergence of |F(X)|.' exit end if end do return end subroutine muller_r ( fatol, itmax, x1, x2, x3, xatol, xrtol ) ! !******************************************************************************* ! !! MULLER_R carries out Muller's method restricted to real roots. ! ! ! Modified: ! ! 18 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real FATOL, the absolute error tolerance for F(X). ! ! Input, integer ITMAX, the maximum number of steps allowed. ! ! Input, real X1, X2, X3, three points to start the iteration. ! ! Input, real XATOL, absolute error tolerance for the root. ! implicit none ! double precision a double precision b double precision c double precision discrm real fatol real fminus real fplus real fxmid real fxnew real fxold integer iterate integer itmax real x_ave real x_inc real x1 real x2 real x3 real xatol real xlast real xmid real xminus real xnew real xold real xplus real xrtol ! xnew = x1 xmid = x2 xold = x3 call f ( xnew, fxnew, 0 ) call f ( xmid, fxmid, 0 ) call f ( xold, fxold, 0 ) write ( *, * ) ' ' write ( *, * ) ' ' write ( *, * ) 'Muller''s method (real root version)' write ( *, * ) ' ' write ( *, * ) 'Note that the existence of a nearby complex root' write ( *, * ) 'is suggested when the value of the discriminant' write ( *, * ) 'is negative.' write ( *, * ) ' ' write ( *, * ) ' ' write ( *, * ) 'Iteration x fx disc' write ( *, * ) ' ' iterate = 0 discrm = 0.0E+00 write ( *, '(1x,i6,f20.10,f20.10,f20.10)' ) iterate, xnew, fxnew if ( abs ( fxnew ) < fatol ) then write ( *, * ) ' ' write ( *, * ) '|F(X)| is below the tolerance.' return end if do if ( abs ( fxnew ) >= abs ( fxmid ) ) then call r_swap ( xnew, xmid ) call r_swap ( fxnew, fxmid ) end if xlast = xnew iterate = iterate + 1 if ( iterate > itmax ) then write ( *, * ) ' ' write ( *, * ) ' Maximum number of steps taken.' exit end if a = dble ( ( xmid - xnew ) * ( fxold - fxnew ) & - ( xold - xnew ) * ( fxmid - fxnew ) ) b = dble ( ( xold - xnew )**2 * ( fxmid - fxnew ) & - ( xmid - xnew )**2 * ( fxold - fxnew ) ) c = dble ( ( xold - xnew ) * ( xmid - xnew ) * ( xold - xmid ) * fxnew ) xold = xmid xmid = xnew ! ! Apply the quadratic formula to get roots XPLUS and XMINUS. ! discrm = b**2 - 4.0E+00 * a * c if ( discrm <= 0.0E+00 ) then discrm = 0.0E+00 end if if ( a == 0.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'The algorithm has broken down.' write ( *, * ) 'The quadratic coefficient A is zero.' exit end if xplus = xnew + sngl ( ( - b + sqrt ( discrm ) ) / ( 2.0E+00 * a ) ) call f ( xplus, fplus, 0 ) xminus = xnew + sngl ( ( - b - sqrt ( discrm ) ) / ( 2.0E+00 * a ) ) call f ( xminus, fminus, 0 ) ! ! Take whichever of the two quadratic roots is closest to a root of ! the function. ! if ( abs ( fminus ) < abs ( fplus ) ) then xnew = xminus else xnew = xplus end if fxold = fxmid fxmid = fxnew call f ( xnew, fxnew, 0 ) write ( *, '(1x,i6,f20.10,f20.10,f20.10)' ) iterate, xnew, fxnew, discrm ! ! Check for convergence. ! x_ave = ( abs ( xnew ) + abs ( xmid ) + abs ( xold ) ) / 3.0E+00 x_inc = xnew - xmid if ( abs ( x_inc ) <= xatol ) then write ( *, * ) ' ' write ( *, * ) ' Absolute convergence of the X increment.' exit end if if ( abs ( x_inc ) <= xrtol * x_ave ) then write ( *, * ) ' ' write ( *, * ) ' Relative convergence of the X increment.' exit end if if ( abs ( fxnew ) <= fatol ) then write ( *, * ) ' ' write ( *, * ) ' Absolute convergence of |F(X)|.' exit end if end do return end subroutine newton ( fatol, itmax, x1, xatol, xrtol ) ! !******************************************************************************* ! !! NEWTON carries out Newton's method. ! ! ! Modified: ! ! 20 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real FATOL, the absolute error tolerance for F(X). ! ! Input, integer ITMAX, the maximum number of steps allowed. ! ! Input, real X1, a point where the iteration begins. ! ! Input, real XATOL, absolute error tolerance for the root. ! implicit none ! integer, parameter :: MAX_ROOT = 20 ! real fatol real fprime real fx integer i integer iroot character isay integer iterat integer itmax real roots(MAX_ROOT) real x_ave real x_inc real x1 real xatol real xnew real xold real xrtol ! roots(1:MAX_ROOT) = 0.0E+00 iroot = 0 10 continue iroot = iroot + 1 if ( iroot > MAX_ROOT ) then return end if write ( *, * ) ' ' write ( *, * ) ' ' write ( *, * ) ' Newton''s method' write ( *, * ) ' ' write ( *, * ) ' ' write ( *, * ) 'Iteration X FX' write ( *, * ) ' ' iterat = 0 xnew = x1 call f ( xnew, fx, 0 ) write ( *, '(i6,f20.10,f20.10)' ) iterat, xnew, fx if ( abs ( fx ) <= fatol ) then write ( *, * ) ' ' write ( *, * ) '|F(X)| is below the tolerance.' return end if do iterat = iterat + 1 if ( iterat > itmax ) then write ( *, * ) ' ' write ( *, * ) ' Maximum number of steps taken.' return end if xold = xnew call f ( xold, fx, 0 ) call f ( xold, fprime, 1 ) do i = 1, iroot - 1 if ( xold - roots(i) == 0.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'The algorithm has broken down.' write ( *, * ) 'A zero divisor was computed.' exit end if fprime = fprime - fx / ( xold - roots(i) ) end do if ( fprime == 0.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'The algorithm has broken down.' write ( *, * ) 'A zero divisor was computed.' exit end if x_inc = - fx / fprime xnew = xold + x_inc call f ( xnew, fx, 0 ) write ( *, '(i6,f20.10,f20.10)' ) iterat, xnew, fx ! ! Check for acceptance of function value or interval ! x_ave = abs ( xnew ) + abs ( xold ) if ( abs ( x_inc ) <= xatol ) then write ( *, * ) ' ' write ( *, * ) ' Absolute convergence of the X increment.' exit end if if ( abs ( x_inc ) <= xrtol * x_ave ) then write ( *, * ) ' ' write ( *, * ) ' Relative convergence of the X increment.' exit end if if ( abs ( fx ) <= fatol ) then write ( *, * ) ' ' write ( *, * ) ' Absolute convergence of |F(X)|.' exit end if end do return end subroutine prlist ( fatol, iform0, iform1, iform2, iformg, itmax, & method, mstart, names, maxmet, x1, x2, x3, xatol, xrtol ) ! !******************************************************************************* ! !! PRLIST prints out the value of the algorithm parameters. ! ! ! Modified: ! ! 19 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! implicit none ! integer maxmet ! real fatol character ( len = * ) iform0 character ( len = * ) iform1 character ( len = * ) iform2 character ( len = * ) iformg integer itmax integer method integer mstart character ( len = 35 ) names(maxmet) real x1 real x2 real x3 real xatol real xrtol ! write ( *, * ) ' ' write ( *, * ) 'Current parameters:' write ( *, * ) ' ' write ( *, * ) 'FATOL Function absolute tolerance = ', fatol write ( *, * ) 'ITMAX Maximum number of iterations =', itmax write ( *, * ) 'METHOD ', names(method) write ( *, * ) 'F(X) = ', trim ( iform0 ) write ( *, * ) 'F''(X) = ', trim ( iform1 ) write ( *, * ) 'F"(X) = ', trim ( iform2 ) write ( *, * ) 'G(X) = ', trim ( iformg ) write ( *, * ) 'Number of starting values needed = ', mstart write ( *, * ) 'X1 Starting point #1 = ', x1 write ( *, * ) 'X2 Starting point #2 = ', x2 write ( *, * ) 'X3 Starting point #3 = ', x3 write ( *, * ) 'XATOL X absolute tolerance = ', xatol write ( *, * ) 'XRTOL X relative tolerance = ', xrtol return end subroutine r_swap ( x, y ) ! !******************************************************************************* ! !! R_SWAP switches two real values. ! ! ! Modified: ! ! 30 November 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, real X, Y. On output, the values of X and ! Y have been interchanged. ! implicit none ! real x real y real z ! z = x x = y y = z return end subroutine rpnchk ( ierror, ihi, ilo, iopsym, irpn, maxrpn, maxsym ) ! !******************************************************************************* ! !! RPNCHK examines an RPN formula, looking for a complete RPN expression. ! ! ! Discussion: ! ! The routine starts at location IHI, and finds the position ILO ! such that IRPN(ILO)...IRPN(IHI) represents a single argument. ! ! Modified: ! ! 24 October 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IERROR, 0, no error; 1 an error. ! ! Input, integer IHI, the location in IRPN where the search begins. ! ! Output, integer ILO, the location in IRPN such that IRPN(ILO) through ! IRPN(IHI) represents a single argument, or IHI+1 if no such location ! was found. ! ! Input, integer IOPSYM(MAXSYM), contains, for each symbol, the number ! of operands. In particular, if a symbol represents a constant, ! IOPSYM(I) is 0. If a symbol represents a unary operator such as ABS, ! IOPSYM(I) is 1. IOPSYM may be as large as 3. ! ! Workspace, integer IRPN(MAXRPN), used to store the compiled ! versions of user formulas. ! ! Input, integer MAXRPN, specifies the length of IRPN. ! ! Input, integer MAXSYM, the maximum number of symbols allowed. ! implicit none ! integer maxrpn integer maxsym ! integer ierror integer ihi integer ilo integer iopsym(maxsym) integer irpn(maxrpn) integer isum ! isum = 0 ierror = 0 ilo = ihi + 1 do ilo = ilo - 1 if ( ilo <= 0 ) then write ( *, * ) ' ' write ( *, * ) 'RPNCHK - Error!' write ( *, * ) ' Reached beginning of formula without matchup.' ierror = 1 return end if if ( iopsym(irpn(ilo)) < 0 .and. isum == 0 ) then cycle end if isum = isum + 1 - iopsym(irpn(ilo)) if ( isum == 1 ) then exit end if end do return end subroutine rpnset ( ierror, ifinis, imin, intsym, iopsym, iprsym, irpn, & istack, maxrpn, maxsym, nints, nrpn, symbol ) ! !******************************************************************************* ! !! RPNSET converts the infix formula into an RPN formula. ! ! ! Modified: ! ! 25 October 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IERROR, error flag. ! 0, no error occurred. ! nonzero, an error occurred. ! ! Input, integer IFINIS, the index in SYMBOL of the symbol ! '$', meaning the end of the formula. ! ! Input, integer IMIN, an offset for accessing IRPN. ! ! Input, integer INTSYM(80), a set of integers which ! are the indices of symbols, representing an infix formula. ! ! Input, integer IOPSYM(MAXSYM), contains, for each ! symbol, the number of operands. In particular, if ! a symbol represents a constant, IOPSYM(I) is 0. ! If a symbol represents a unary operator such as ABS, ! IOPSYM(I) is 1. IOPSYM may be as large as 3. ! ! Input, integer IPRSYM(MAXSYM), the priority of each symbol. ! ! Output, integer IRPN(MAXRPN), the RPN version of the infix formula. ! ! Workspace, integer ISTACK(MAXSYM). ! ! Input, integer MAXRPN, the maximum number of RPN symbols allowed. ! ! Input, integer MAXSYM, the maximum number of symbols allowed. ! ! Input, integer NINTS, the number of useful entries in INTSYM. ! ! Output, integer NRPN, the number of useful entries in IRPN. ! ! Input, character ( len = MAXLEN ) SYMBOL(MAXSYM), the symbolic names ! of the internally defined constants and functions. ! implicit none ! integer, parameter :: maxlen = 20 integer maxrpn integer maxsym ! integer ierror integer ifinis integer imin integer intsym(80) integer iopsym(maxsym) integer iprsym(maxsym) integer iread integer irpn(maxrpn) integer istack(maxsym) integer istak integer isym integer jsym integer nints integer nrpn character ( len = maxlen ) sym character ( len = maxlen ) sym2 character ( len = maxlen ) symbol(maxsym) ! ierror = 0 nrpn = 0 ! ! An error if the formula seems to have nothing in it. ! if ( nints <= 0 ) then write ( *, * ) ' ' write ( *, * ) 'RPNSET - Error!' write ( *, * ) ' This formula seems to be blank!' ierror = 1 return end if iread = 0 istak = 0 ! ! Read symbol number IREAD from INTSYM. ! do iread = iread + 1 if ( iread > nints ) then ierror = 1 write ( *, * ) ' ' write ( *, * ) 'RPNSET - Error!' write ( *, * ) ' This formula does not make sense!' return end if isym = intsym(iread) if ( isym == 0 ) then cycle end if sym = symbol(isym) if ( sym == ',' ) then cycle end if ! ! If the symbol is "$", then it's time to pop the stack. ! if ( sym == '$' ) then do if ( istak <= 0 ) then nrpn = nrpn + 1 irpn(nrpn+imin) = ifinis if ( iread < nints ) then write ( *, * ) ' ' write ( *, * ) 'RPNSET - Error!' write ( *, * ) ' Some of the formula is left over!' ierror = 1 end if return end if nrpn = nrpn + 1 irpn(nrpn+imin) = istack(istak) istak = istak - 1 end do end if ! ! A left parenthesis goes onto the stack. ! if ( sym == '(' ) then istak = istak + 1 istack(istak) = isym cycle end if ! ! A variable or constant goes immediately into IRPN. ! if ( iopsym(isym) == 0 ) then nrpn = nrpn + 1 irpn(nrpn+imin) = isym cycle end if ! ! Put operators onto the stack. ! do if ( istak <= 0 ) then istak = istak + 1 istack(istak) = isym exit end if jsym = istack(istak) if ( iprsym(jsym) <= iprsym(isym) ) then jsym = istack(istak) sym2 = symbol(jsym) if ( sym == ')' .and. sym2 == '(' ) then istak = istak - 1 else istak = istak + 1 istack(istak) = isym end if exit end if nrpn = nrpn + 1 irpn(nrpn+imin) = istack(istak) istak = istak - 1 end do end do return end subroutine rpnval ( ierror, ifree, imin, iopsym, iprsym, ipval, & irad, irpn, irtol, istack, maxrpn, maxsym, maxval, & nrpn, nsym, nsyms, symbol, valsym, value ) ! !******************************************************************************* ! !! RPNVAL evaluates the symbolic functions in an RPN formula. ! ! ! Discussion: ! ! RPNVAL determines the number of arguments, gathering their values, ! and calling FUNVAL to evaluate the given function. ! ! Modified: ! ! 01 November 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real EPS, the value of the machine precision. ! ! Output, integer IERROR. ! 0, no error occurred. ! nonzero, an error occurred. ! ! Input, integer IFREE, the index of the next free location in VALSYM. ! ! Input, integer IMIN, an offset for indexing IRPN. ! ! Input, integer IOPSYM(MAXSYM), contains, for each ! symbol, the number of operands. In particular, if ! a symbol represents a constant, IOPSYM(I) is 0. ! If a symbol represents a unary operator such as ABS, ! IOPSYM(I) is 1. IOPSYM may be as large as 3. ! ! Input, integer IPRSYM(MAXSYM), the relative priorities ! of the functions. ! ! Input, integer IPVAL(MAXSYM), contains, for each symbol, ! the address in VALSYM where the value of the symbol ! is stored, if it is a scalar. If the symbol represents ! a vector or matrix, IPVAL(IARG) points to location of ! the first entry. ! ! Workspace, integer IRPN(MAXRPN), used to store the compiled ! versions of user formulas. ! ! Workspace, integer ISTACK(MAXSYM), workspace for interpreting ! the formula. ! ! Input, integer MAXRPN, specifies the length of IRPN. ! ! Input, integer MAXSYM, the maximum number of symbols allowed. ! ! Input, integer MAXVAL, the maximum number of values allowed in VALSYM. ! ! Output, integer NCOL, the number of columns in the result. ! ! Output, integer NROW, the number of rows in the result. ! ! Input, integer NRPN, the number of useful entries in IRPN. ! ! Input, integer NSYM, the total number of symbols, including temporaries. ! ! Input, integer NSYMS, the number of declared symbols. ! ! Input, integer NUMDIM(2,MAXSYM). ! For each symbol I, NUMDIM(1,I) is the number of rows in its ! value, and NUMDIM(2,I) is the number of columns. ! ! Input, character ( len = * ) SYMBOL(MAXSYM), the symbolic names ! of the internally defined constants and functions. ! ! Input, real VALSYM(MAXVAL), contains the values of all the ! symbolic variables. ! ! Workspace, real VALUE, space used to hold ! the value of the variable to be saved. ! implicit none ! integer, parameter :: maxlen = 20 integer maxrpn integer maxsym ! character ( len = 3 ) ctemp integer iarg1 integer iarg2 integer ierror integer ifree integer imin integer index1 integer index4 integer iopsym(maxsym) integer iprsym(maxsym) integer ipset integer ipval(maxsym) integer irad integer iread integer irpn(maxrpn) integer irtol integer istack(maxsym) integer istak integer isym integer isymo integer itemp integer maxval integer nrpn integer nsym integer nsyms character ( len = maxlen ) sym character ( len = maxlen ) symbol(maxsym) real valsym(maxval) real value ! value = 0.0E+00 if ( nrpn <= 1 ) then return end if istak = 0 isymo = 0 nsyms = nsym itemp = 0 iread=0 do iread = iread + 1 if ( iread > nrpn ) then ierror = 1 write ( *, * ) ' ' write ( *, * ) 'RPNVAL - Error!' write ( *, * ) ' The formula ended unexpectedly!' exit end if isym = irpn(iread+imin) if ( 0 < isym .and. isym <= maxsym ) then sym = symbol(isym) else sym = '$' end if if ( sym == '$' .or. iread > nrpn ) then index4 = ipval(isymo) value = valsym(index4) exit end if ! ! Constants and variables go into stack. ! if ( iopsym(isym) == 0 ) then if ( isym <= nsym ) then istak = istak+1 istack(istak) = isym isymo = isym else if ( nsyms >= maxsym ) then write ( *, * ) ' ' write ( *, * ) 'RPNVAL - Error!' write ( *, * ) ' There is no memory left for more symbols!' write ( *, * ) ' Perhaps the KILL or INIT command would help.' ierror = 1 exit end if nsyms = nsyms+1 isymo = nsyms istak = istak+1 istack(istak) = nsyms iprsym(nsyms) = 10 iopsym(nsyms) = 0 itemp = itemp+1 ctemp = ' ' call i_to_s_zero ( itemp, ctemp ) symbol(nsyms) = 'STK000' symbol(nsyms)(4:6) = ctemp ipval(nsyms) = ifree ifree = ifree + 1 index1 = ipval(isym) index4 = ipval(nsyms) valsym(index4) = valsym(index1) end if cycle end if ! ! Pull off arguments. ! iarg1 = istack(istak) iarg2 = 0 if ( iopsym(isym) == 2 ) then iarg2 = istack(istak) istak = istak-1 iarg1 = istack(istak) else if (iopsym(isym) == 3 ) then istak = istak-1 iarg2 = istack(istak) istak = istak-1 iarg1 = istack(istak) end if sym = symbol(isym) call funval ( iarg1, iarg2, ierror, ifree, iopsym, iprsym, & ipset, ipval, irad, irtol, itemp, maxsym, maxval, nsyms, & sym, symbol, valsym ) isymo = nsyms if ( ierror /= 0 ) then write ( *, * ) ' ' write ( *, * ) 'RPNVAL - Warning!' write ( *, * ) ' Evaluation of this formula is abandoned.' value = 0.0E+00 exit end if if ( ipset == 0 ) then ipval(nsyms) = ifree ifree = ifree + 1 else ipval(nsyms) = ipset ipset = 0 end if istack(istak) = nsyms end do return end subroutine s_blank_delete ( s ) ! !******************************************************************************* ! !! S_BLANK_DELETE removes blanks from a string, left justifying the remainder. ! ! ! Comment: ! ! All TAB characters are also removed. ! ! Modified: ! ! 26 July 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, character ( len = * ) S, the string to be transformed. ! implicit none ! character c integer iget integer iput character ( len = * ) s character ( len = 1 ), parameter :: TAB = char ( 9 ) ! iput = 0 do iget = 1, len ( s ) c = s(iget:iget) if ( c /= ' ' .and. c /= TAB ) then iput = iput + 1 s(iput:iput) = c end if end do s(iput+1:) = ' ' return end subroutine s_cap ( string ) ! !******************************************************************************* ! !! S_CAP replaces any lowercase letters by uppercase ones in a string. ! ! ! Modified: ! ! 16 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, character ( len = * ) STRING, the string to be transformed. ! implicit none ! character c integer i integer nchar character ( len = * ) string ! nchar = len ( string ) do i = 1, nchar c = string(i:i) call ch_cap ( c ) string(i:i) = c end do return end function s_eqi ( strng1, strng2 ) ! !******************************************************************************* ! !! S_EQI is a case insensitive comparison of two strings for equality. ! ! ! Examples: ! ! S_EQI ( 'Anjana', 'ANJANA' ) is .TRUE. ! ! Modified: ! ! 14 April 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) STRNG1, STRNG2, the strings to compare. ! ! Output, logical S_EQI, the result of the comparison. ! implicit none ! integer i integer len1 integer len2 integer lenc logical s_eqi character s1 character s2 character ( len = * ) strng1 character ( len = * ) strng2 ! len1 = len ( strng1 ) len2 = len ( strng2 ) lenc = min ( len1, len2 ) s_eqi = .false. do i = 1, lenc s1 = strng1(i:i) s2 = strng2(i:i) call ch_cap ( s1 ) call ch_cap ( s2 ) if ( s1 /= s2 ) then return end if end do do i = lenc + 1, len1 if ( strng1(i:i) /= ' ' ) then return end if end do do i = lenc + 1, len2 if ( strng2(i:i) /= ' ' ) then return end if end do s_eqi = .true. return end function s_is_alpha ( s ) ! !******************************************************************************* ! !! S_IS_ALPHA returns .TRUE. if the string contains only alphabetic characters. ! ! ! Discussion: ! ! Here, alphabetic characters are 'A' through 'Z' and 'a' through 'z'. ! ! Modified: ! ! 14 April 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) S, the string to be checked. ! ! Output, logical S_IS_ALPHA, .TRUE. if the string contains only ! alphabetic characters, .FALSE. otherwise. ! implicit none ! logical ch_is_alpha integer i character ( len = * ) s logical s_is_alpha ! s_is_alpha = .false. do i = 1, len ( s ) if ( .not. ch_is_alpha ( s(i:i) ) ) then return end if end do s_is_alpha = .true. return end function s_paren_check ( s ) ! !******************************************************************************* ! !! S_PAREN_CHECK checks the parentheses in a string. ! ! ! Discussion: ! ! Blanks are removed from the string, and then the following checks ! are made: ! ! 1) as we read the string left to right, there must never be more ! right parentheses than left ones; ! 2) there must be an equal number of left and right parentheses; ! 3) there must be no occurrences of the abutting packages '...)(...'. ! 4) there must be no occurrences of the empty package '()'. ! ! Modified: ! ! 20 November 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) S, the string to check. ! ! Output, logical S_PAREN_CHECK is TRUE if the string passed the checks. ! implicit none ! integer i integer isum character ( len = * ) s character ( len = 256 ) s_copy integer s_len logical s_paren_check ! s_copy = s call s_blank_delete ( s_copy) s_len = len_trim ( s_copy ) ! ! 1) Letting '(' = +1 and ')' = -1, check that the running parentheses sum ! is always nonnegative. ! isum = 0 do i = 1, s_len if ( s_copy(i:i) == '(' ) then isum = isum + 1 end if if ( s_copy(i:i) == ')' ) then isum = isum - 1 if ( isum < 0 ) then s_paren_check = .false. return end if end if end do ! ! 2) Check that the final parentheses sum is zero. ! if ( isum /= 0 ) then s_paren_check = .false. return end if ! ! 3) Check that there are no ")(" pairs. ! do i = 2, s_len if ( s_copy(i-1:i) == ')(' ) then s_paren_check = .false. return end if end do ! ! 4) Check that there are no "()" pairs. ! do i = 2, s_len if ( s_copy(i-1:i) == '()' ) then s_paren_check = .false. return end if end do ! ! The checks were passed. ! s_paren_check = .true. return end subroutine s_to_i ( s, ival, ierror, last ) ! !******************************************************************************* ! !! S_TO_I reads an integer value from a string. ! ! ! 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 LAST, the last character of S used to make IVAL. ! implicit none ! character c integer i integer ierror integer isgn integer istate integer ival integer last 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 last = 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 last = len_trim ( s ) else ierror = 1 last = 0 end if return end subroutine s_to_r ( s, r, ierror, lchar ) ! !******************************************************************************* ! !! S_TO_R reads a real number from a string. ! ! ! Discussion: ! ! This 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 real number. ! ! Legal input is: ! ! 1 blanks, ! 2 '+' or '-' sign, ! 2.5 spaces ! 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. ! ! Examples: ! ! S R ! ! '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) ! ! Modified: ! ! 12 February 2001 ! ! 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 R, the real value that was 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 LCHAR, the number of characters read from ! the string to form the number, including any terminating ! characters such as a trailing comma or blanks. ! implicit none ! logical ch_eqi character c integer ierror integer ihave integer isgn integer iterm integer jbot integer jsgn integer jtop integer lchar integer nchar integer ndig real r real rbot real rexp real rtop character ( len = * ) s character, parameter :: TAB = char ( 9 ) ! nchar = len_trim ( s ) ierror = 0 r = 0.0E+00 lchar = - 1 isgn = 1 rtop = 0.0E+00 rbot = 1.0E+00 jsgn = 1 jtop = 0 jbot = 1 ihave = 1 iterm = 0 do lchar = lchar + 1 c = s(lchar+1:lchar+1) ! ! Blank or TAB character. ! if ( c == ' ' .or. c == TAB ) then if ( ihave == 2 ) then else if ( ihave == 6 .or. ihave == 7 ) then iterm = 1 else if ( ihave > 1 ) then ihave = 11 end if ! ! Comma. ! else if ( c == ',' .or. c == ';' ) then if ( ihave /= 1 ) then iterm = 1 ihave = 12 lchar = lchar + 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 ( ihave >= 6 .and. ihave <= 8 ) then ihave = 9 else iterm = 1 end if ! ! 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. lge ( c, '0' ) .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.0E+00 * rtop + real ( ndig ) else if ( ihave == 5 ) then rtop = 10.0E+00 * rtop + real ( ndig ) rbot = 10.0E+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 .or. lchar+1 >= nchar ) 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 LCHAR is equal to NCHAR. ! if ( iterm /= 1 .and. lchar+1 == nchar ) then lchar = 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 return end if ! ! Number seems OK. Form it. ! if ( jtop == 0 ) then rexp = 1.0E+00 else if ( jbot == 1 ) then rexp = 10.0E+00**( jsgn * jtop ) else rexp = jsgn * jtop rexp = rexp / jbot rexp = 10.0E+00**rexp end if end if r = isgn * rexp * rtop / rbot return end subroutine secant ( fatol, itmax, x1, x2, xatol, xrtol ) ! !******************************************************************************* ! !! SECANT carries out the secant algorithm. ! ! ! Modified: ! ! 18 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real FATOL, the absolute error tolerance for F(X). ! ! Input, integer ITMAX, the maximum number of steps allowed. ! ! Input, real X1, X2, two points at which the iteration begins. ! ! Input, real XATOL, absolute error tolerance for the root. ! implicit none ! real fatol real fx real fxnew real fxold integer iterate integer itmax real x_ave real x_inc real x1 real x2 real xatol real x real xnew real xold real xrtol ! write ( *, * ) ' ' write ( *, * ) ' ' write ( *, * ) 'The secant method:' write ( *, * ) ' ' write ( *, * ) ' ' write ( *, * ) 'Step X FX' write ( *, * ) ' ' iterate = -1 xold = x1 call f ( xold, fxold, 0 ) write ( *, '(i4,2g16.8)' ) iterate, xold, fxold if ( abs ( fxold ) <= fatol ) then write ( *, * ) ' ' write ( *, * ) '|F(X)| is below the tolerance.' return end if iterate = 0 x = x2 call f ( x, fx, 0 ) write ( *, '(i4,2g16.8)' ) iterate, x, fx if ( abs ( fx ) <= fatol ) then write ( *, * ) ' ' write ( *, * ) ' Absolute convergence of |F(X)|.' return end if do iterate = iterate + 1 if ( iterate > itmax ) then write ( *, * ) ' ' write ( *, * ) ' Maximum number of steps taken.' exit end if x_inc = x - xold x_ave = ( abs ( x ) + abs ( xold ) ) / 2.0E+00 if ( abs ( x_inc ) <= xatol ) then write ( *, * ) ' ' write ( *, * ) ' Absolute convergence of the X increment.' exit end if if ( abs ( x_inc ) <= xrtol * x_ave ) then write ( *, * ) ' ' write ( *, * ) ' Relative convergence of the X increment.' exit end if if ( abs ( fx ) <= fatol ) then write ( *, * ) ' ' write ( *, * ) ' Absolute convergence of |F(X)|.' exit end if ! ! Compute the next iterate. ! xnew = ( fxold * x - fx * xold ) / ( fxold - fx ) call f ( xnew, fxnew, 0 ) write ( *, '(i4,2g16.8)' ) iterate, xnew, fxnew ! ! Prepare for the next step. ! xold = x fxold = fx x = xnew fx = fxnew end do return end subroutine solve ( fatol, itmax, method, mform, x1, x2, x3, xatol, xrtol ) ! !******************************************************************************* ! !! SOLVE carries out the iteration procedure. ! ! ! Modified: ! ! 18 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real FATOL, the absolute error tolerance for F(X). ! ! Input, integer ITMAX, the maximum number of steps allowed. ! ! Input, real X1, X2, X3, three starting points. ! ! Input, real XATOL, XRTOL, the absolute and relative error ! tolerances for the root. ! implicit none ! real fatol integer itmax integer method integer mform real x1 real x2 real x3 real xatol real xrtol ! ! Call the requested method ! if ( method == 1 ) then call bisection ( fatol, itmax, x1, x2, xatol, xrtol ) else if ( method == 2 ) then call falsi ( fatol, itmax, x1, x2, xatol, xrtol ) else if ( method == 3 ) then call falsi2 ( fatol, itmax, x1, x2, xatol, xrtol ) else if ( method == 4 ) then call secant ( fatol, itmax, x1, x2, xatol, xrtol ) else if ( method == 5 ) then call muller_r ( fatol, itmax, x1, x2, x3, xatol, xrtol ) else if ( method == 6 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SOLVE - Warning!' write ( *, '(a)' ) ' MULLER_C is not operational.' ! call muller_c ( fatol, itmax, x1, x2, x3, xatol, xrtol ) else if ( method == 7 ) then call newton ( fatol, itmax, x1, xatol, xrtol ) else if ( method == 8 ) then call steffenson ( fatol, itmax, x1, xatol, xrtol ) else if ( method == 9 ) then call laguerre ( fatol, itmax, x1, xatol, xrtol ) else if ( method == 10 ) then call fixed ( fatol, itmax, x1, xatol, xrtol ) else if ( method == 11 ) then call halley ( fatol, itmax, x1, xatol, xrtol ) else if ( method == 12 ) then call brent ( fatol, itmax, x1, x2, xatol, xrtol ) else write ( *, * ) ' ' write ( *, * ) 'Unrecognized method!' end if return end subroutine steffenson ( fatol, itmax, x1, xatol, xrtol ) ! !******************************************************************************* ! !! STEFFENSON carries out Steffenson's method. ! ! ! Modified: ! ! 20 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real FATOL, the absolute error tolerance for F(X). ! ! Input, integer ITMAX, the maximum number of steps allowed. ! ! Input, real X1, a point where the iteration begins. ! ! Input, real XATOL, absolute error tolerance for the root. ! implicit none ! real fatol real fx real fxtemp integer iterate integer itmax real x_ave real x_inc real x1 real xatol real xnew real xold real xsum real xrtol real xtemp ! write ( *, * ) ' ' write ( *, * ) ' ' write ( *, * ) 'Steffenson''s method' write ( *, * ) ' ' write ( *, * ) ' ' write ( *, * ) 'Iteration X FX' write ( *, * ) ' ' iterate = 0 xnew = x1 call f ( xnew, fx, 0 ) write ( *, '(1x,i6,f20.10,f20.10)' ) iterate, xnew, fx if ( abs ( fx ) <= fatol ) then write ( *, * ) ' ' write ( *, * ) '|F(X)| is below the tolerance.' return end if do iterate = iterate + 1 if ( iterate > itmax ) then write ( *, * ) ' ' write ( *, * ) ' Maximum number of steps taken.' exit end if xold = xnew call f ( xold, fx, 0 ) xtemp = xold + fx call f ( xtemp, fxtemp, 0 ) if ( fxtemp == fx ) then write ( *, * ) ' ' write ( *, * ) 'The algorithm has broken down.' write ( *, *) 'A zero divisor was computed.' return end if x_inc = - fx**2 / ( fxtemp - fx ) xnew = xold + x_inc call f ( xnew, fx, 0 ) write ( *, '(1x,i6,f20.10,f20.10)' ) iterate, xnew, fx ! ! Check for convergence. ! x_ave = ( abs ( xnew ) + abs ( xold ) ) / 2.0E+00 x_inc = abs ( xnew - xold ) if ( abs ( x_inc ) <= xatol ) then write ( *, * ) ' ' write ( *, * ) ' Absolute convergence of the X increment.' exit end if if ( abs ( x_inc ) <= xrtol * x_ave ) then write ( *, * ) ' ' write ( *, * ) ' Relative convergence of the X increment.' exit end if if ( abs ( fx ) <= fatol ) then write ( *, * ) ' ' write ( *, * ) ' Absolute convergence of |F(X)|.' exit end if end do return end subroutine symadd ( ierror, ifree, iopsym, iprsym, ipval, & maxsym, maxval, namvr, nsym, nsymp, symbol, valsym ) ! !*********************************************************************** ! !! SYMADD adds a symbol name to the list of symbolic names. ! ! ! Modified: ! ! 01 November 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IERROR, error flag. ! 0, no error occurred, variable was added to list. ! 1, error occurred, variable was not added to list. ! ! Input/output, integer IFREE, the index of the next free ! memory location in VALSYM. ! ! Input, integer IOPSYM(MAXSYM), contains, for each ! symbol, the number of operands. In particular, if ! a symbol represents a constant, IOPSYM(I) is 0. ! If a symbol represents a unary operator such as ABS, ! IOPSYM(I) is 1. IOPSYM may be as large as 3. ! ! Input, integer IPRSYM(MAXSYM), the relative priorities ! of the functions. ! ! Input, integer IPVAL(MAXSYM), contains, for each symbol, ! the address in VALSYM where the value of the symbol ! is stored, if it is a scalar. If the symbol represents ! a vector or matrix, IPVAL(IARG) points to location of ! the first entry. ! ! Input, integer MAXSYM, the maximum number of symbols allowed. ! ! Input, integer MAXVAL, the maximum number of values allowed ! in VALSYM. ! ! Input, character ( len = MAXLEN ) NAMVR, the name of the variable. ! ! Input, integer NSYM, the total number of symbols, including temporaries. ! ! Input, integer NSYMP, the number of permanent symbols. ! ! Input, integer NUMDIM(2,MAXSYM). ! For each symbol I, NUMDIM(1,I) is the number of rows in its ! value, and NUMDIM(2,I) is the number of columns. ! ! Input, character ( len = MAXLEN ) SYMBOL(MAXSYM), the symbolic names ! of the internally defined constants and functions. ! ! Input, real VALSYM(MAXVAL), contains the values of all the ! symbolic variables. ! implicit none ! integer, parameter :: maxlen = 20 integer maxsym integer maxval ! integer i integer ierror integer ifree integer indx integer iopsym(maxsym) integer iprsym(maxsym) integer ipval(maxsym) integer lennam character ( len = maxlen ) namvr integer nsym integer nsymp logical s_eqi character ( len = maxlen ) symbol(maxsym) real valsym(maxval) ! ! Get the length of the variable name. ! ierror = 0 lennam = len_trim ( namvr ) if ( lennam <= 0 ) then write ( *, * ) ' ' write ( *, * ) 'SYMADD - Error!' write ( *, * ) ' This variable name has zero length!' ierror = 1 return end if ! ! Check if the name is already in use. ! do i = 1, nsymp if ( s_eqi ( namvr, symbol(i) ) ) then write ( *, * ) ' ' write ( *, * ) 'SYMADD - Error!' write ( *, * ) ' The name ', trim ( namvr ), ' is already in use.' ierror = 1 return end if end do ! ! Check for user overriding earlier use of same name. ! do i = nsymp+1, nsym if ( s_eqi ( namvr, symbol(i) ) ) then symbol(i) = 'VOID' end if end do ! ! Is there enough space in SYMBOL for another symbol name? ! if ( nsym >= maxsym ) then write ( *, * ) ' ' write ( *, * ) 'SYMADD - Error!' write ( *, * ) ' There is not enough memory for ', trim ( namvr ) write ( *, * ) ' Perhaps the KILL or INIT commands would help.' ierror = 1 return end if ! ! Is there enough space in VALSYM for the symbol's value? ! if ( ifree > maxval ) then ierror = 1 write ( *, * ) ' ' write ( *, * ) 'SYMADD - Error!' write ( *, * ) ' There is not enough memory for ', trim ( namvr ) write ( *, * ) ' Perhaps the KILL or INIT commands would help.' return end if ! ! Insert information about the symbol into various tables. ! nsym = nsym + 1 symbol(nsym) = namvr iopsym(nsym) = 0 iprsym(nsym) = 10 ipval(nsym) = ifree ifree = ifree + 1 indx = ipval(nsym) valsym(indx) = 0.0E+00 return end subroutine symbol_value ( com, ierror, ifree, iopsym, iprsym, ipval, maxsym, & maxval, namvar, nsym, nsymp, symbol, valsym, value ) ! !******************************************************************************* ! !! SYMBOL_VALUE sets, evaluates, or deletes a variable. ! ! ! Modified: ! ! 14 March 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character COM, determines what task is to be done: ! 'R', "reads" the variable, returning its value and dimensions. ! 'V', "sets" the value and dimensions of the variable. ! 'D', "deletes" the variable from memory. ! ! Output, integer ierror, error flag. ! 0, no error occurred. ! nonzero, an error occurred. ! ! Input/output, integer IFREE, the index of the next free ! memory location in VALSYM. ! ! Input, integer IOPSYM(MAXSYM), contains, for each ! symbol, the number of operands. In particular, if ! a symbol represents a constant, IOPSYM(i) is 0. ! If a symbol represents a unary operator such as abs, ! IOPSYM(i) is 1. IOPSYM may be as large as 3. ! ! Input, integer IPRSYM(MAXSYM), the function priorities. ! ! Input, integer IPVAL(MAXSYM), contains, for each symbol, ! the address in VALSYM where the value of the symbol ! is stored, if it is a scalar. If the symbol represents ! a vector or matrix, IPVAL(IARG) points to location of ! the first entry. ! ! Input, integer MAXSYM, the maximum number of symbols allowed. ! ! Input, integer MAXVAL, the maximum number of values allowed. ! ! Input, character ( len = maxlen ) namvar, the name of the variable. ! ! Input, integer NSYM, the total number of symbols, including temporaries. ! ! Input, integer NSYMP, the number of permanent symbols. ! ! Input, character ( len = maxlen ) SYMBOL(MAXSYM), symbolic names. ! ! Input, real VALSYM(MAXVAL), the values of the symbolic variables. ! ! Input/output, real VALUE, holds the value of the variable. ! implicit none ! integer, parameter :: maxlen = 20 integer maxsym integer maxval ! character com integer i integer ierror integer ifree integer indx integer iopsym(maxsym) integer iprsym(maxsym) integer ipval(maxsym) integer match character ( len = maxlen ) namvar integer nsym integer nsymp logical s_eqi character ( len = maxlen ) symbol(maxsym) real valsym(maxval) real value ! ierror = 0 if ( len_trim ( namvar ) <= 0 ) then ierror = 1 write ( *, * ) ' ' write ( *, * ) 'SYMBOL_VALUE - Error!' write ( *, * ) ' Bad variable name ' // trim ( namvar ) return end if ! ! Search for a match between namvar and any predefined symbol. ! match = 0 do i = 1, nsym if ( s_eqi ( namvar, symbol(i) ) ) then match = i exit end if end do ! ! If no such variable seems to exist, but we have been asked ! to DELETE or READ such a variable, then exit. ! if ( match == 0 ) then if ( s_eqi ( com, 'D' ) .or. s_eqi ( com, 'R' ) ) then write ( *, * ) ' ' write ( *, * ) 'SYMBOL_VALUE - Error!' write ( *, * ) ' No variable named ' // namvar if ( s_eqi ( com, 'R' ) ) then ierror = 1 return end if end if ! ! Add the name of the new variable to the internal list. ! call symadd ( ierror, ifree, iopsym, iprsym, ipval, & maxsym, maxval, namvar, nsym, nsymp, symbol, valsym ) if ( ierror /= 0 ) then return end if match = nsym ! ! The matched symbol has index MATCH. ! else ! ! If the symbol is not the name of a variable, but ! rather the name of an operator, then this is an error. ! if ( iopsym(match) /= 0 ) then ierror = 1 write ( *, * ) ' ' write ( *, * ) 'SYMBOL_VALUE - Error!' write ( *, * ) ' You are trying to assign a value to' write ( *, * ) ' the name of a function or operator:' write ( *, * ) symbol(match) return end if ! ! Some variables may not be revalued. ! if ( s_eqi ( com, 'V' ) ) then if ( s_eqi ( symbol(match), 'E' ) ) ierror = 1 if ( s_eqi ( symbol(match), 'EPS' ) ) ierror = 1 if ( s_eqi ( symbol(match), 'PI' ) ) ierror = 1 if ( symbol(match)(1:1) == '"' ) ierror = 1 end if ! ! Some variables may not be deleted. ! if ( s_eqi ( com, 'D' ) ) then if ( s_eqi ( symbol(match), 'E' ) ) ierror = 1 if ( s_eqi ( symbol(match), 'EPS' ) ) ierror = 1 if ( s_eqi ( symbol(match), 'PI' ) ) ierror = 1 if ( s_eqi ( symbol(match), 'RTOL' ) ) ierror = 1 if ( s_eqi ( symbol(match), 'SEED' ) ) ierror = 1 if ( symbol(match)(1:1) == '"' ) ierror = 1 end if if ( ierror /= 0 ) then ierror = 1 write ( *, * ) ' ' write ( *, * ) 'SYMBOL_VALUE - Error!' write ( *, * ) ' Try to change special variable:' write ( *, * ) symbol(match) return end if end if indx = ipval(match) if ( s_eqi ( com, 'V' ) ) then valsym(indx) = value else if ( s_eqi ( com, 'R' ) ) then value = valsym(indx) end if return end subroutine timestamp ( ) ! !******************************************************************************* ! !! TIMESTAMP prints the current YMDHMS date as a time stamp. ! ! ! Example: ! ! May 31 2001 9:45:54.872 AM ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! None ! implicit none ! character ( len = 8 ) ampm integer d character ( len = 8 ) date 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 character ( len = 10 ) time integer values(8) integer y character ( len = 5 ) zone ! call date_and_time ( date, time, zone, values ) y = values(1) m = values(2) d = values(3) h = values(5) n = values(6) s = values(7) mm = values(8) if ( h < 12 ) then ampm = 'AM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Noon' else ampm = 'PM' end if else h = h - 12 if ( h < 12 ) then ampm = 'PM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Midnight' else ampm = 'AM' end if end if end if write ( *, '(a,1x,i2,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & trim ( month(m) ), d, y, h, ':', n, ':', s, '.', mm, trim ( ampm ) return end subroutine tokens ( ierror, ifinis, ifree, implic, indx1, indx2, ineg, infix, & intsym, iopsym, iprsym, ipval, maxsym, maxval, nints, nsym, nsymp, & nuvar, symbol, valsym ) ! !******************************************************************************* ! !! TOKENS parses a character string into recognized symbols and constants. ! ! ! Discussion: ! ! The routine produces the output array of integers INTSYM, containing ! whose entries stand for recognized symbols, variable names or ! constants. ! ! Modified: ! ! 02 November 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IERROR, error flag. ! 0, no error occurred. ! nonzero, an error occurred. ! ! Input, integer IFINIS, the index in SYMBOL of the symbol ! '$', meaning the end of the formula. ! ! Input, integer IFREE, the address of the next free ! memory location in VALSYM. ! ! Output, integer IMPLIC, implicit variable flag. ! 0, the assignment variable was not implicit. ! nonzero, the assignment variable was implicit, and has been ! added to the symbol table, and assigned the symbol index of IMPLIC. ! ! Input, integer INDX1, the index of the symbol 'INDEX1'. ! ! Input, integer INDX2, the index of the symbol 'INDEX2'. ! ! Input, integer INEG, the index of the symbol 'NEG'. ! ! Input, character ( len = * ) INFIX, a sequence of characters representing ! an infix formula. ! ! Output, integer INTSYM(80), a sequence of integers which are ! the indices of the symbols used in INFIX. ! ! Input, integer IOPSYM(MAXSYM), contains, for each ! symbol, the number of operands. In particular, if ! a symbol represents a constant, IOPSYM(I) is 0. ! If a symbol represents a unary operator such as ABS, ! IOPSYM(I) is 1. IOPSYM may be as large as 3. ! ! Input, integer IPRSYM(MAXSYM), the relative priorities ! of the functions. ! ! Input, integer IPVAL(MAXSYM), contains, for each symbol, ! the address in VALSYM where the value of the symbol ! is stored, if it is a scalar. If the symbol represents ! a vector or matrix, IPVAL(IARG) points to location of ! the first entry. ! ! Input, integer MAXSYM, the maximum number of symbols allowed. ! ! Input, integer MAXVAL, the maximum number of values allowed in VALSYM. ! ! Output, integer NINTS, the number of useful entries in INTSYM. ! ! Input, integer NSYM, the total number of symbols, including temporaries. ! ! Input, integer NSYMP, the number of permanent symbols. ! ! Input, integer NUVAR. ! If NUVAR is 0, then the formula may refer to variables which ! have not yet been declared. ! Otherwise, all variables on the right hand side of the "=" ! sign must have been previously declared. ! ! Input, character ( len = MAXLEN ) SYMBOL(MAXSYM), the symbolic names ! of the internally defined constants and functions. ! ! Input, real VALSYM(MAXVAL), contains the values of all the ! symbolic variables. ! implicit none ! integer, parameter :: maxlen = 20 integer maxsym ! character ctemp integer i integer iback integer icom1 integer icom2 integer ierror integer ifinis integer ifree integer ilpr integer implic integer indx1 integer indx2 integer ineg character ( len = * ) infix integer intsym(80) integer iopsym(maxsym) integer ipoint integer ipos integer iprsym(maxsym) integer ipval(maxsym) integer irpr integer j integer jback integer lchar integer lenfix integer lenmat integer lens integer loc integer loc1 integer loc2 integer match integer matcho integer maxval character ( len = maxlen ) namvr integer nints logical notnum integer nsym integer nsymp integer nuvar real rval logical s_eqi logical s_is_alpha character sym character ( len = maxlen ) sym1 character ( len = maxlen ) sym2 character ( len = maxlen ) sym3 character ( len = maxlen ) symbol(maxsym) real valsym(maxval) ! implic = 0 ierror = 0 nints = 0 match = 0 ipos = 1 lenfix = len_trim ( infix ) ! ! 1: Does this formula begin with a variable name followed ! by an equals sign? ! sym = '=' loc1 = index ( infix, sym ) if ( loc1 <= 1 ) then go to 20 end if ! ! The formula begins with a variable name followed by an equals sign. ! ! Pick off the name of the variable, but first make sure ! you stop before any left hand parenthesis, in case the ! formula is something like "alfred(7)=5". ! sym = '(' loc2 = index ( infix, sym ) if ( loc2 <= 0 ) then loc = loc1 else loc = min ( loc1, loc2 ) end if if ( loc <=1 .or. loc > maxlen+1 ) then go to 20 end if loc = loc - 1 namvr = infix(1:loc) ! ! Check to see whether the name of the assignment variable has ! already been defined. ! do i = 1, nsym if ( s_eqi ( namvr, symbol(i) ) ) then go to 20 end if end do ! ! The assignment variable has not previously been defined. ! Prepare to add it to the table of symbols. ! ipos = loc + 1 write ( *, * ) ' ' write ( *, * ) 'TOKENS:' write ( *, * ) ' Implicit definition of ', trim ( namvr ) call symadd ( ierror, ifree, iopsym, iprsym, ipval, maxsym, maxval, namvr, & nsym, nsymp, symbol, valsym ) if ( ierror /= 0 ) then return end if nints = nints + 1 intsym(nints) = nsym match = nsym implic = nsym ! ! Process the next chunk of the formula. ! 20 continue ! ! If we have reached the end of the input string, we're done. ! Append an "End-of-formula" marker "$" and return. ! if ( ipos > lenfix ) then nints = nints + 1 intsym(nints) = ifinis return end if ! ! We're not done yet. Remember the last match in MATCHO, and ! look for the next match. Go through all the symbols, and ! try to find the longest match possible. ! matcho = match match = 0 lenmat = 0 do i = 1, nsym lens = len_trim ( symbol(i) ) if ( lens > lenmat .and. ipos + lens - 1 <= lenfix ) then if ( s_eqi ( infix(ipos:ipos+lens-1), symbol(i) ) ) then lenmat = lens match = i end if end if end do if ( match == 0 ) then go to 110 end if ! ! We found a match. ! ! But watch out! The user might be implicitly defining a name, ! but because the first part matches an earlier symbol, you don't ! realize it! For example, in a formula involving SINGLE, you ! might think you see "SIN". So make sure that the match you've ! got is not followed by an alphabetic character. ! if ( iopsym(match) /= 0 ) then go to 40 end if if ( ipos+lenmat-1 >= lenfix ) then go to 40 end if ctemp = infix(ipos+lenmat:ipos+lenmat) if ( s_is_alpha ( ctemp ) ) then go to 140 end if ! ! We have matched an old symbol. ! 40 continue sym1 = symbol(match) sym2 = ' ' if ( matcho /= 0 ) then sym2 = symbol(matcho) end if ! ! Check for unary minus or plus. ! if ( sym1 == '-' ) then if ( sym2 == '**' .or. sym2 == ',' .or. & sym2 == '(' .or. sym2 == '=' .or. & ipos == 1 ) then sym = infix(ipos+1:ipos+1) if ( lge ( sym, '0' ) .and. lle ( sym, '9' ) ) then go to 110 end if if ( sym == '.' ) then go to 110 end if end if end if if ( sym1 == '-' ) then if ( ipos == 1 .or. sym2 == '(' .or. sym2 == ',' .or. sym2 == '=' ) then nints = nints + 1 intsym(nints) = ineg ipos = ipos + 1 go to 20 end if end if if ( sym1 == '+' ) then if ( ipos == 1 .or. sym2 == '(' .or. sym2 == '=' ) then ipos = ipos + 1 go to 20 end if end if ! ! If we've hit a right parenthesis, then we assume we're dealing with: ! ! a vector index: A(3) ! a matrix index: A(1,2) ! a two place operator: MAX(A,B), or MIN or POLY ! a three place operator: KRYLOV(A,X,N) ! ! We need to find the matching left parenthesis, count the ! arguments by assuming that commas separate them, and read ! off the name preceding the left parenthesis. ! if ( sym1 == ')' ) then icom1 = 0 icom2 = 0 ilpr = 0 irpr = 1 do iback = 1, nints i = nints + 1 - iback if ( intsym(i) <= 0 ) then go to 90 end if sym3 = symbol ( intsym(i) ) if ( sym3 == ',' .and. irpr-ilpr == 1 ) then if ( icom1 == 0 ) then icom1 = i else if ( icom2 == 0 ) then icom2 = i else ierror = 1 write ( *, * ) ' ' write ( *, * ) 'TOKENS - Error!' write ( *, * ) ' Too many commas in formula.' return end if else if ( sym3 == '(' ) then ilpr = ilpr + 1 else if ( sym3 == ')' ) then irpr = irpr + 1 end if ! ! If ILPR=IRPR, we've reached the matching left parenthesis. ! if ( ilpr == irpr ) then if ( i <= 1 ) then go to 100 end if if ( intsym(i-1) <= 0 ) then go to 100 end if sym3 = symbol ( intsym(i-1) ) if ( iopsym(intsym(i-1)) /= 0 .and. & .not. s_eqi ( sym3, 'MAX' ) .and. & .not. s_eqi ( sym3, 'MIN' ) .and. & .not. s_eqi ( sym3, 'ATAN2' ) .and. & .not. s_eqi ( sym3, 'POLY' ) .and. & .not. s_eqi ( sym3, 'KRYLOV' ) ) then go to 100 end if ! ! Copy ) into INTSYM. ! nints = nints + 1 intsym(nints) = match ! ! INDX1 is our name for a vector reference. ! if ( icom1 == 0 ) then match = indx1 ! ! INDX2 is our name for a matrix reference. ! else if ( icom2 == 0 ) then intsym(icom1) = match do jback = 1, nints-icom1 j = nints + 1 - jback intsym(j+1) = intsym(j) end do nints = nints + 1 intsym(icom1+1) = match - 1 if ( iopsym(intsym(i-1)) /= 0 ) then ipos = ipos + 1 go to 20 end if match = indx2 ! ! Insert )( to divide a trio of arguments. ! else intsym(icom1) = match do jback = 1, nints-icom1 j = nints + 1 - jback intsym(j+1) = intsym(j) end do intsym(icom1+1) = match - 1 nints = nints + 1 intsym(icom2) = match do jback = 1, nints-icom2 j = nints + 1 - jback intsym(j+1) = intsym(j) end do intsym(icom2+1) = match - 1 nints = nints + 1 if ( iopsym(intsym(i-1)) /= 0 ) then ipos = ipos + 1 go to 20 end if match = indx2 end if ipos = ipos + 1 nints = nints + 1 intsym(nints) = match go to 20 end if 90 continue end do end if 100 continue lens = len_trim ( symbol(match) ) nints = nints + 1 intsym(nints) = match ipos = ipos + lens go to 20 ! ! Check for a constant. ! 110 continue call s_to_r ( infix(ipos:), rval, ierror, lchar ) if ( lchar <= 0 ) then go to 140 end if if ( ierror /= 0 ) then write ( *, * ) ' ' write ( *, * ) 'TOKENS - Error!' write ( *, * ) ' Illegal real number.' return end if ! ! Look for a constant already defined with the same value. ! do i = nsymp+1, nsym if ( iprsym(i) == 10 ) then if ( symbol(i)(1:1) == '"' ) then ipoint = ipval(i) if ( valsym(ipoint) == rval ) then match = i go to 130 end if end if end if end do ! ! This is a new constant. Make up a name for it. ! if ( real ( int ( rval ) ) == rval ) then write ( namvr, '(''"'',i19)' ) int ( rval ) else write ( namvr, '(''"'', g19.10)' ) rval end if call s_blank_delete ( namvr ) ! ! Add the constant to the symbol list. ! call symadd ( ierror, ifree, iopsym, iprsym, ipval, maxsym, maxval, namvr, & nsym, nsymp, symbol, valsym ) if ( ierror /= 0 ) then return end if ipoint = ipval(nsym) valsym(ipoint) = rval match = nsym 130 continue ! ! Advance LCHAR characters, except that CHRCTR will read in a ! comma as the end of the number, and we have to back up in ! such a case. ! if ( infix(ipos+lchar-1:ipos+lchar-1) /= ',' ) then ipos = ipos + lchar else ipos = ipos + lchar - 1 end if nints = nints + 1 intsym(nints) = match go to 20 ! ! Consider implicit variable declaration on right hand side. ! 140 continue if ( nuvar == 0 ) then ierror = 1 write ( *, * ) ' ' write ( *, * ) 'TOKENS - Error!' write ( *, * ) ' Undeclared variable in formula.' return end if namvr = ' ' do i = 1, maxlen sym = infix(ipos-1+i:ipos-1+i) notnum = llt ( sym, '0' ) .or. lgt ( sym, '0' ) if ( ( .not. s_is_alpha(sym) ) .and. notnum ) then exit end if namvr(i:i) = sym end do lchar = len_trim ( namvr ) if ( lchar <= 0 ) then ierror = 1 write ( *, * ) ' ' write ( *, * ) 'TOKENS - Error!' write ( *, '(a)' ) ' The string starting at location ', ipos, & ' is unreadable.' return end if ipos = ipos + lchar call symadd ( ierror, ifree, iopsym, iprsym, ipval, maxsym, maxval, namvr, & nsym, nsymp, symbol, valsym ) nints = nints + 1 intsym(nints) = nsym match = nsym write ( *, * ) ' ' write ( *, * ) 'TOKENS:' write ( *, * ) ' The undeclared variable ', trim ( namvr ), & ' will be assumed to be a scalar.' go to 20 end