C
c this is version 2.22b of march 1991.
C
C WRITTEN BY: DR. JOSEPH F. GRCAR
C SANDIA NATIONAL LABORATORY
C SCIENTIFIC COMPUTING DIVISION 8233
C LIVERMORE, CALIFORNIA 94551-0969 USA
C
C (415) 294-2662
C (FTS) 234-2662
C
C///////////////////////////////////////////////////////////////////////
C
C CHANGES FROM THE PREVIOUS VERSION:
C
C 1) ALLOW NEGATIVE STEP S0 AS FLAG TO AVOID NEWTON SEARCH.
C
C///////////////////////////////////////////////////////////////////////
C
C ARGUMENTS FOR SUBROUTINE TWOPNT:
C
C ERROR output logical - error flag. if true, then a programming
C error or a severe numerical error blocks execution. error
C messages appear in the output text file (no matter the
C output level).
C
C TEXT input integer - unit number for an output file. all text
C output goes to this unit. 'level' controls the quantity
C of output, but zero and negative values for 'text' suspend
C all output (including error messages).
C
C LEVEL input integer dimensioned (2) - output level. in the tree
C of subroutine calls stemming from twopnt, 'level(1)'
C specifies the highest routines that write informative
C messages to the output text file. if 'level(1)' equals 1,
C then only twopnt writes messages; if 'level(1)' equals 2,
C then subroutines called directly by twopnt also write
C messages; and so on. if 'level(1)' is zero or negative,
C then only error messages appear. if the unit number
C 'text' is zero or negative, then not even error messages
C appear.
C
C 'level(2)' specifies the highest routines that write the
C solution to the output file. if 'level(2)' equals 1, then
C the solution appears at the start and finish of twopnt; if
C 'level(2)' equals 2, then the solution appears at the
C start of twopnt and at the finish of every subroutine
C called directly by twopnt (the solution does not reappear
C when twopnt itself finishes); and so on. 'level(2)' must
C be less than or equal to 'level(2)'.
C
C ISIZE input integer - size of the integer work space array
C 'iwork'. must be at least 3 * 'pmax'.
C
C IWORK integer dimensioned (isize) - integer work space.
C
C RSIZE real integer - size of the real work space array 'rwork'.
C must be at least 7 * 'scalrs' + (7 * 'comps' + 2) *
C 'pmax'.
C
C RWORK real dimensioned (isize) - real work space.
C
C ABOVE input/output real dimensioned (scalrs + comps * pmax) -
C upper bounds. each solution variable has its own bound.
C no solution variable exceeds its upper bound at any time.
C if mesh refinement occurs, then new points receive bounds
C interpolated from surrounding points. the user may
C overwrite these during 'update' reverse communication
C requests.
C
C ACTIVE input logical dimensioned (comps) - markers for active
C components. only active components of the solution govern
C mesh refinement. component j is active when 'active(j)'
C is true.
C
C ADAPT input logical - refinement directive. if true, then
C permit mesh refinement.
C
C BELOW input/output real dimensioned (scalrs + comps * pmax) -
C lower bounds. each solution variable has its own bound.
C no solution solution variable falls short of its lower
C bound at any time. if mesh refinement occurs, then new
C points receive bounds interpolated from surrounding
C points. the user may overwrite these during 'update'
C reverse communication requests.
C
C BINARY input integer - unit number for an output file. this unit
C receives binary output that summarizes the progress of
C twopnt. program 'show' graphs the data. zero and
C negative values for 'binary' suspend all output.
C
C BUFFER input/output real dimensioned (scalrs + comps * pmax) -
C reverse communication buffer.
C
C COMPS input integer - number of components.
C
C CONDIT input real - condition number. 'condit' passes the
C condition number of the latest jacobian matrix. zero and
C negative values mean the number is unknown.
C
C FUNCT output logical - reverse communication signal. if true,
C then (1) evaluate the residual function at the approximate
C solution in 'buffer', (2) store the result in 'buffer',
C and (3) reenter twopnt. when 'time' is also true, then
C the residual function should be the one for the
C time-dependent problem.
C
C JACOB output logical - reverse communication signal. if true,
C then (1) evaluate the jacobian matrix of the residual
C function at the approximate solution in 'buffer', (2)
C prepare to solve systems of linear equations having this
C matrix of coefficients (for example, factor the matrix),
C (3) optionally place the condition number of the matrix in
C 'condit', and (4) call twopnt again. when 'time' is also
C true, then the residual function should be the one for the
C time-dependent problem.
C
C MARK output logical dimensioned (pmax) - markers for new
C points. if 'update' is true, then 'mesh(k)' is new
C exactly when 'mark(k)' is true.
C
C MESH input/output real dimensioned (pmax) - mesh points. this
C array contains the values of the independent spatial
C variable at which the steady-state and time-dependent
C functions discretize the differential equations. the user
C chooses the initial points. the points must lie in
C increasing or decreasing order. if 'adapt' is true, then
C twopnt may add points up to the capacity of the array.
C 'points' is the current number.
C
C PASS output integer - number of the latest problem pass on the
C latest mesh. see 'passes'.
C
C PASSES input integer - number of passes per mesh. if it is
C desirable to pose several problems (perhaps of increasing
C difficulty) on each mesh, then set passes to the number of
C problems. otherwise, set passes to one.
C
C PLIMIT input integer - maximum number of points that may be
C simultaneously added to the mesh.
C
C PMAX input integer - maximum number of mesh points. pmax
C determines the size of many arrays.
C
C POINTS input/output integer - current number of mesh points. the
C components for each point lie together after the scalar
C variables in arrays 'above', 'below', 'buffer', and 'x'.
C
C REENT input/output logical - reenter signal. set false on input
C to mark the start of a new problem. set true on output
C when twopnt requests information. in this case, perform
C the task indicated by one of the reverse communication
C signals ('funct', 'jacob', 'save', 'show', 'solve',
C 'store', or 'update') and reenter twopnt.
C
C REPORT output integer - failure note. if both 'error' and
C 'succes' are false, then 'report' explains why (otherwise
C 'report' returns zero). (1) no change: the newton
C algorithm did not solve the steady-state problem, and the
C subsequent time step did not change the solution. the
C time step is probably too small. (2) no solution: the
C newton algorithm did not solve a time-dependent problem.
C the time step is probably too large. (3) no space: the
C mesh refinement criteria are not satisfied on the largest
C possible grid.
C
C SAVE output logical - reverse communication signal. if true,
C then the latest approximate solution may be used as a
C starting estimate should twopnt have to be restarted (due
C to system crash or time out, for example). thus, (1)
C optionally save 'mesh', 'points', and the solution found
C in 'buffer', and (2) reenter twopnt.
C
C SCALRS input integer - number of scalar variables. these precede
C the mesh variables in the arrays x and buffer.
C
C SHOW output logical - reverse communication signal. if true,
C then (1) write to unit 'text' the approximate solution in
C 'buffer', and (2) reenter twopnt. the user can display
C the solution with problem dependent data, such as
C component labels, that are unknown to twopnt. 'level(2)'
C determines when twopnt requests the solution be shown.
C
C SOLVE output logical - reverse communication signal. if true,
C then (1) solve a system of linear equations having the
C right hand side in 'buffer' and using the most recent
C jacobian matrix, (2) store the solution in 'buffer', and
C (3) reenter twopnt.
C
C SSABS input real - absolute bound for the steady-state problem.
C the newton algorithm terminates if each entry of the next
C (undamped) step is either absolutely small or small
C relative to the corresponding entry of the approximate
C solution.
C
C SSAGE input integer - age of the steady state jacobian. number
C of newton steps to be taken before recomputing the
C jacobian.
C
C SSREL input real - relative bound for the steady-state problem.
C the newton terminates if each entry of the next (undamped)
C step is either absolutely small or small relative to the
C corresponding entry of the approximate solution.
C
C STEPS0 input integer - number of initial time steps. if
C positive, then this many steps should be completed of the
C time-dependent problem before attempting the newton
C algorithm. if negative, newton is not attempted.
C
C STEPS1 input integer - number of smoothing time steps. if the
C newton algorithm cannot solve the steady-state problem,
C then at most this many steps should be attempted of the
C time-dependent problem.
C
C STEPS2 input integer - number of time steps to be taken before
C increasing the time step.
C
C STORE output logical - reverse communication signal. if true,
C then twopnt is about to perform or is performing a time
C integration. 'buffer' contains either the initial value
C or the solution from the last time step. the
C time-dependent function needs this to evaluate time
C derivatives. thus, (1) store the solution, and (2)
C reenter twopnt.
C
C SUCCES output logical - completion status. if true, then twopnt
C solves the boundary value problem to the extent that the
C solution in 'x' satisfies the mesh refinement criteria.
C if false, then the solution in 'x' does not satisfy the
C criteria ('report' may explain why).
C
C TDABS input real - absolute bound for the time-dependent
C problem. the newton algorithm terminates if each entry of
C the next (undamped) step is either absolutely small or
C small relative to the corresponding entry of the
C approximate solution.
C
C TDAGE input integer - age of the time dependent jacobian.
C number of newton steps to be taken before recomputing the
C jacobian during time stepping.
C
C TDEC input real - factor by which to divide the time step when
C necessary.
C
C TDREL input real - relative bound for the time-dependent
C problem. the newton algorithm terminates if each entry of
C the next (undamped) step is either absolutely small or
C small relative to the corresponding entry of the
C approximate solution.
C
C TIME output logical - reverse communication signal. if true,
C then satisfy the accompanying request for 'funct' or
C 'jacob' using the time-dependent function.
C
C TINC input real - factor by which to multiply the time step
C when the number of steps specified by 'retire' has passed.
C
C TMAX input real - maximum time step.
C
C TMIN input real - minimum time step.
C
C TOLER0 input real - absolute and relative significance level.
C the mesh refinement criteria ('toler1' and 'toler2') apply
C only to active, significant components. a component is
C insignificant when its change (maximum less minimum) is
C either absolutely small or small relative to its maximum
C magnitude.
C
C TOLER1 input real - bound upon the range of a component. the
C first mesh refinement criterion requires that each
C component's change over adjacent mesh points be no more
C than the fraction 'toler1' of the component's change
C (maximum less minimum) over all mesh points. if some
C component violates this criterion over some mesh interval,
C then twopnt adds the interval's midpoint to the mesh.
C
C TOLER2 input real - bound upon the range of a component's
C derivative. the second mesh refinement criterion requires
C that the change in each component's derivative over
C adjacent mesh intervals be no more than the fraction
C 'toler2' of the derivative's change (maximum less minimum)
C over all mesh intervals. if some component violates this
C criterion over some pair of adjacent mesh intervals, then
C twopnt adds the intervals' midpoints to the mesh.
C
C TSTEP0 input real - initial time step.
C
C TSTEP1 output real - latest size of the time step to be used in
C evaluating the function (when time steps are being taken).
C
C UPDATE output logical - reverse communication signal. if true,
C then 'mesh' contains some new points (marked by 'mark').
C twopnt chooses interpolated bounds and interpolated
C solution estimates for the new points. these may be
C inadequate. thus (1) optionally adjust 'above', 'below',
C and the approximate solution in 'buffer', and (2) reenter
C twopnt.
C
C X input/output real dimensioned (scalrs + comps * pmax) -
C solution. on input, an estimate. on output, the
C solution.
C
C///////////////////////////////////////////////////////////////////////
c
SUBROUTINE COPY (N, X, Y)
IMPLICIT COMPLEX (A - Z)
C*****precision > double
DOUBLE PRECISION X, Y
C*****END precision > double
INTEGER J, N
C*****precision > single
C REAL X, Y
C*****END precision > single
DIMENSION X(N), Y(N)
DO 0100 J = 1, N
Y(J) = X(J)
0100 CONTINUE
END
c
SUBROUTINE CPUTIM (TIMER)
C///////////////////////////////////////////////////////////////////////
C
C CPUTIM
C
C OBTAIN ELAPSED CPU JOB TIME IN SECONDS.
C
C///////////////////////////////////////////////////////////////////////
IMPLICIT COMPLEX (A - Z)
C*****precision > double
DOUBLE PRECISION TIMER
C*****END precision > double
C*****precision > single
C REAL TIMER
C*****END precision > single
C THE SYSTEM LIBRARY BASELIB CONTAINS SUBROUTINE IZM00.
C*****CRAY/CTSS (LIVERMORE ENVIRONMENT)
C EXTERNAL IZM00
C INTEGER FLAG, IZM00, VALUE
C REAL TEMP
C DIMENSION VALUE(5)
C
C FLAG = IZM00 (VALUE)
C IF(FLAG.EQ.0) THEN
C TEMP = REAL (VALUE(1)) / 1.0E6
C ELSE
C TEMP = 0.0
C ENDIF
c*****end cray/ctss (livermore environment)
c the system library cftlib contains subroutine isecond.
c*****cray/ctss (los alamos environment)
C EXTERNAL ISECOND
C INTEGER VALUE
C REAL TEMP
C
C CALL ISECOND (VALUE)
C TEMP = REAL (VALUE) / 1.0E6
C*****END CRAY/CTSS (LOS ALAMOS ENVIRONMENT)
C*****CRAY/UNICOS
C REAL SECOND, TEMP
C TEMP = SECOND ()
C*****END CRAY/UNICOS
C*****SUN/UNIX
C EXTERNAL ETIME
C REAL ETIME, LIST
C DIMENSION LIST(2)
C
C TEMP = ETIME (LIST)
C*****END SUN/UNIX
C THE VAX/VMS VERSION USES SYSTEM SERVICE CALLS.
C*****VAX/VMS
C INTEGER*2 LIST
C INTEGER*4 ADDRESS, TICS
C REAL*4 TEMP
C
C DIMENSION LIST(6)
C
C EQUIVALENCE (LIST(3), ADDRESS)
C
C DATA LIST(1), LIST(2), LIST(5), LIST(6) / 4, '0407'X, 0, 0 /
C
C ADDRESS = %LOC (TICS)
C CALL SYS$GETJPIW (, , , LIST, , , )
C TEMP = REAL (TICS) / 1.0E2
C*****END VAX/VMS
C*****PHOENICS
REAL TEMP
INTEGER ITIME(20)
external MCLOCK
CALL MCLOCK(ITIME)
TEMP=60.*(ITIME(5)+60.*ITIME(4))+ITIME(6)
1 +0.01*(ITIME(7)+0.01*ITIME(8))
C*****END PHOENICS
C*****precision > double
TIMER = DBLE (TEMP)
C*****END precision > double
C*****precision > single
C TIMER = TEMP
C*****END precision > single
END
c
SUBROUTINE TPELPS (TIMER)
IMPLICIT COMPLEX (A - Z)
C*****precision > double
DOUBLE PRECISION TEMP, TIMER
C*****END precision > double
C*****precision > single
C REAL TEMP, TIMER
C*****END precision > single
CALL CPUTIM (TEMP)
TIMER = TEMP - TIMER
END
c
SUBROUTINE EXTENT (LENGTH, STRING)
C///////////////////////////////////////////////////////////////////////
C
C EXTENT
C
C FIND THE LAST NONBLANK CHARACTER IN A STRING.
C
C///////////////////////////////////////////////////////////////////////
IMPLICIT COMPLEX (A - P, R - Z), INTEGER (Q)
CHARACTER STRING*(*)
INTEGER J, LENGTH
LENGTH = 1
DO 0100 J = 1, LEN (STRING)
IF(STRING(J : J).NE.' ') LENGTH = J
0100 CONTINUE
END
c
SUBROUTINE GRAB (ERROR, TEXT, LAST, FIRST, NUMBER)
C///////////////////////////////////////////////////////////////////////
C
C GRAB
C
C RESERVE SPACE IN A DATA SPACE.
C
C///////////////////////////////////////////////////////////////////////
IMPLICIT COMPLEX (A - Z)
CHARACTER ID*9
INTEGER FIRST, LAST, NUMBER, TEXT
LOGICAL ERROR
PARAMETER (ID = ' GRAB: ')
C/// CHECK.
ERROR = .NOT. (0.LE.LAST)
IF(ERROR) GO TO 9001
ERROR = .NOT. (0.LE.NUMBER)
IF(ERROR) GO TO 9002
C/// GRAB THE SPACE.
FIRST = LAST + 1
LAST = LAST + MAX (1, NUMBER)
C/// ERROR MESSAGES.
GO TO 99999
9001 IF(0.LT.TEXT) WRITE (TEXT, 99001) ID, LAST
GO TO 99999
9002 IF(0.LT.TEXT) WRITE (TEXT, 99002) ID, NUMBER
GO TO 99999
99001 FORMAT
+ (/1X, A9, 'ERROR. THE POINTERS FOR THE DATA SPACE ARE'
+ /10X, 'INCONSISTENT.'
+ //10X, I10, ' LAST ENTRY CURRENTLY IN USE')
99002 FORMAT
+ (/1X, A9, 'ERROR. THE NUMBER OF WORDS TO BE RESERVED IN THE'
+ /10X, 'DATA SPACE MUST NOT BE NEGATIVE.'
+ //10X, I10, ' NUMBER')
C/// EXIT.
99999 CONTINUE
END
c
SUBROUTINE LOGSTR (FORMAT, LENGTH, STRING, VALUE)
IMPLICIT COMPLEX (A - Z)
CHARACTER FORMAT*(*), STRING*(*), TEMP*80
INTEGER FIRST, J, LAST, LENGTH
LOGICAL FOUND
C*****precision > double
DOUBLE PRECISION VALUE, ZERO
C*****END precision > double
C*****precision > single
C REAL VALUE, ZERO
C*****END precision > single
C/// FORM THE WORD.
ZERO = 0.0
IF(ZERO.LT.VALUE) THEN
WRITE (TEMP, FORMAT, ERR = 0200) LOG10 (VALUE)
ELSEIF(VALUE.EQ.ZERO) THEN
TEMP = 'NULL'
ELSE
TEMP = '?'
ENDIF
C/// FIND THE FIRST AND LAST CHARACTERS.
FOUND = .FALSE.
LAST = 0
DO 0100 J = 1, LEN (TEMP)
IF(TEMP (J : J).NE.' ') THEN
LAST = J
IF(.NOT. FOUND) FIRST = J
FOUND = .TRUE.
ENDIF
0100 CONTINUE
C/// PACK THE WORD INTO THE STRING.
IF(0.LT.LENGTH .AND. LENGTH.LE.LEN (STRING)) THEN
STRING (1 : LENGTH) = ' '
IF(FOUND) THEN
IF((LAST - FIRST).LT.LENGTH) THEN
STRING (LENGTH - (LAST - FIRST) : LENGTH)
+ = TEMP (FIRST : LAST)
ELSE
GO TO 0200
ENDIF
ENDIF
ELSE
GO TO 0200
ENDIF
C/// SET THE STRING TO '!!! ...'.
GO TO 99999
0200 CONTINUE
DO 0300 J = 1, MIN (LEN (STRING), LENGTH)
STRING (J : J) = '!'
0300 CONTINUE
C/// EXIT.
99999 CONTINUE
END
c
SUBROUTINE NEWTON
+ (ERROR, TEXT, ABOVE, ACCEPT, BELOW, BUFFER, CONDIT, DECIDE,
+ EXIST, FUNCT, JACOB, LEVEL1, LEVEL2, N, REENT, REPORT, RETIRE,
+ S, SHOW, SOLVE, STEPS, STMAX, STRIAL, SUCCES, X, XTRIAL, Y,
+ YNORM, YTRIAL)
IMPLICIT COMPLEX (A - P, R - Z), INTEGER (Q)
CHARACTER COLUMN*16, HEADER*80, ID*9
EXTERNAL COPY, LOGSTR
INTEGER
+ AGE, COUNT, ENTRY, FORMED, INDEX, J, K, LEVEL1, LEVEL2, LIMIT,
+ N, REDUCE, REPORT, RETIRE, RETURN, ROUTE, STEPS, STMAX, TEXT
LOGICAL
+ ACCEPT, DECIDE, ERROR, EXIST, FORCE, FUNCT, JACOB, PRINTJ,
+ PRINTS, PRINTY, REENT, SHOW, SOLVE, SUCCES
C*****precision > double
DOUBLE PRECISION
C*****END precision > double
C*****precision > single
C REAL
C*****END precision > single
+ ABOVE, BELOW, BUFFER, COEFF, CONDIT, HALF, ONE, S, SNORM,
+ STNORM, STRIAL, TEMP, VALUE, X, XTRIAL, Y, YNORM, YTNORM,
+ YTRIAL, ZERO
PARAMETER (ID = 'NEWTON: ')
PARAMETER (LIMIT = 5)
C REPORT CODES
PARAMETER
+ (QNULL = 0, QDONE = 1, QDVRG = 2, QBNDS = 3, QSTDY = 4,
+ QSMAX = 5)
DIMENSION
+ ABOVE(N), BELOW(N), BUFFER(N), COLUMN(7), HEADER(4, 2), S(N),
+ STRIAL(N), X(N), XTRIAL(N), Y(N), YTRIAL(N)
C/// SAVE LOCAL VALUES DURING RETURNS FOR REVERSE COMMUNCIATION.
SAVE
C FORMED: SOLUTION INDEX AT WHICH A JACOBIAN WAS LAST COMPUTED.
C INDEX: NUMBER OF THE LATEST ACCEPTED SOLUTION.
C///////////////////////////////////////////////////////////////////////
C
C INITIALIZATION.
C
C///////////////////////////////////////////////////////////////////////
C/// EVERY-TIME INITIALIZATION.
C TURN OFF ALL REVERSE COMMUNICATION FLAGS.
DECIDE = .FALSE.
FUNCT = .FALSE.
JACOB = .FALSE.
SHOW = .FALSE.
SOLVE = .FALSE.
C TURN OFF ALL COMPLETION STATUS FLAGS.
ERROR = .FALSE.
REPORT = QNULL
SUCCES = .FALSE.
C PRECISION-INDEPENDENT CONSTANTS.
ONE = 1.0
HALF = 0.5
ZERO = 0.0
C/// IF THIS IS A RETURN CALL, THEN CONTINUE WHERE THE PROGRAM PAUSED.
IF(REENT) THEN
REENT = .FALSE.
GO TO ROUTE
ENDIF
C/// ONE-TIME INITIALIZATION.
FORMED = - 1
INDEX = 0
PRINTJ = .FALSE.
PRINTS = .FALSE.
PRINTY = .FALSE.
C/// CHECK THE ARGUMENTS.
C ERROR = .NOT. (1.LE.N)
IF(.NOT. (1.LE.N)) ERROR = .TRUE.
IF(ERROR) GO TO 9001
COUNT = 0
DO 1100 J = 1, N
IF(.NOT. BELOW(J).LT.ABOVE(J)) COUNT = COUNT + 1
1100 CONTINUE
C ERROR = COUNT.NE.0
IF(COUNT.NE.0) ERROR = .TRUE.
IF(ERROR) GO TO 9002
COUNT = 0
DO 1200 J = 1, N
IF(.NOT. (BELOW(J).LE.X(J) .AND. X(J).LE.ABOVE(J)))
+ COUNT = COUNT + 1
1200 CONTINUE
C ERROR = COUNT.NE.0
IF(COUNT.NE.0) ERROR = .TRUE.
IF(ERROR) GO TO 9003
C/// PRINT THE HEADER.
C 123456789_123456789_123456789_123456789_12345
C 123456 123456 123456 123456 123456
HEADER (1, 1) = ' LOG10 MAX NORM '
HEADER (2, 1) = ' --------------------------------- '
HEADER (3, 1) = 'SOLUTN FUNCTN STEP TRIAL TRIAL '
HEADER (4, 1) = 'NUMBER VALUE SIZE VALUE STEP '
C 123456789_12345
C 123456 123456
HEADER (1, 2) = ' LOG10'
HEADER (2, 2) = ' LOG10 JACOBN'
HEADER (3, 2) = 'DAMPNG CONDTN'
HEADER (4, 2) = ' COEFF NUMBER'
IF(LEVEL1 .GE. 1) THEN
IF(TEXT.GT.0) WRITE (TEXT, 10001)
+ ID, ((HEADER(J, K), K = 1, 2), J = 1, 4)
ENDIF
C/// EVALUATE THE FUNCTION AT THE INITIAL X AND STORE IN Y.
PRINTY = .TRUE.
ASSIGN 1300 TO RETURN
GO TO 9911
1300 CONTINUE
CALL TPNORM (N, YNORM, Y)
C///////////////////////////////////////////////////////////////////////
C
C BLOCK FOR ONE STEP OF NEWTON'S ALGORITHM.
C
C///////////////////////////////////////////////////////////////////////
C/// FORM THE JACOBIAN AT THE LATEST X, AND RE-EVALUATE THE FUNCTION
C/// IN CASE THE FUNCTION HAS CHANGED, UNLESS THE ALGORITHM HAS JUST
C/// BEGUN AND A JACOBIAN ALREADY EXISTS.
IF(INDEX.EQ.0 .AND. EXIST) GO TO 2300
0200 CONTINUE
FORMED = INDEX
PRINTJ = .TRUE.
ASSIGN 2100 TO RETURN
GO TO 9921
2100 CONTINUE
PRINTY = .TRUE.
ASSIGN 2200 TO RETURN
GO TO 9911
2200 CONTINUE
CALL TPNORM (N, YNORM, Y)
2300 CONTINUE
C/// APPLY THE INVERSE JACOBIAN TO Y AND STORE THE STEP DIRECTION IN S.
PRINTS = .TRUE.
ASSIGN 2400 TO RETURN
GO TO 9931
2400 CONTINUE
CALL TPNORM (N, SNORM, S)
C/// IF THE SOLUTION IS ACCEPTABLE, THEN EXIT WITH SUCCESS.
ASSIGN 2500 TO RETURN
GO TO 9941
2500 CONTINUE
ACCEPT = ACCEPT .AND. (1.LE.INDEX)
IF(.NOT. ACCEPT) GO TO 2800
IF(LEVEL1 .GE. 1) THEN
C PRINT A LINE FOR THE FINAL SOLUTION.
DO 2600 J = 1, 7
COLUMN(J) = ' '
2600 CONTINUE
WRITE (COLUMN(1), '(I6)') INDEX
CALL LOGSTR ('(F6.2)', 6, COLUMN(2), YNORM)
CALL LOGSTR ('(F6.2)', 6, COLUMN(3), SNORM)
IF(PRINTJ) CALL LOGSTR ('(F6.2)', 6, COLUMN(7), CONDIT)
IF(TEXT.GT.0) WRITE (TEXT, 10002) COLUMN
C PRINT THE FINAL SOLUTION.
IF(LEVEL1 .GE. LEVEL2 .AND. LEVEL2 .GE. 1) THEN
IF(TEXT.GT.0) WRITE (TEXT, 10003) ID
ASSIGN 2700 TO RETURN
IF(TEXT.GT.0) GO TO 9951
ELSEIF(LEVEL1 .GE. 1) THEN
IF(TEXT.GT.0) WRITE (TEXT, 10004) ID
ENDIF
ENDIF
2700 CONTINUE
REPORT = QDONE
SUCCES = .TRUE.
GO TO 99999
2800 CONTINUE
C/// FIND THE LARGEST DAMPING COEFFICIENT BETWEEN 0 AND 1 THAT KEEPS
C/// THE TRIAL X WITHIN BOUNDS. IF THE NEXT STEP WILL PLACE SOME
C/// VARIABLES AT THEIR UPPER OR LOWER BOUNDS, THEN MAKE PROVISIONS TO
C/// FORCE ONE OF THESE VARIABLES TO THE APPROPRIATE BOUNDARY.
0300 CONTINUE
COEFF = ONE
FORCE = .FALSE.
DO 3100 J = 1, N
IF(S(J).GT.ZERO) THEN
TEMP = (X(J) - BELOW(J)) / S(J)
IF(TEMP.LE.COEFF) THEN
COEFF = TEMP
ENTRY = J
FORCE = .TRUE.
VALUE = BELOW(J)
ENDIF
ELSEIF(S(J).LT.ZERO) THEN
TEMP = (X(J) - ABOVE(J)) / S(J)
IF(TEMP.LE.COEFF) THEN
COEFF = TEMP
ENTRY = J
FORCE = .TRUE.
VALUE = ABOVE(J)
ENDIF
ENDIF
3100 CONTINUE
C ERROR = COEFF.LT.ZERO
IF(COEFF.LT.ZERO) ERROR = .TRUE.
IF(ERROR) GO TO 9004
REDUCE = 0
C/// IF THE COEFFICIENT IS ZERO, THEN SOME VARIABLES LIE AT THEIR
C/// UPPER OR LOWER BOUNDS AND THE PRESENT STEP DIRECTION S WOULD
C/// CARRY THEM OUT OF BOUNDS. S MUST BE BAD. PRINT A LINE FOR THE
C/// FAILED S.
IF(COEFF.EQ.ZERO) THEN
IF(LEVEL1 .GE. 1) THEN
DO 3200 J = 1, 7
COLUMN(J) = ' '
3200 CONTINUE
IF(PRINTY) WRITE (COLUMN(1), '(I6)') INDEX
IF(PRINTY) CALL LOGSTR ('(F6.2)', 6, COLUMN(2), YNORM)
IF(PRINTS) CALL LOGSTR ('(F6.2)', 6, COLUMN(3), SNORM)
COLUMN(6) = ' ZERO'
IF(PRINTJ) CALL LOGSTR ('(F6.2)', 6, COLUMN(7), CONDIT)
IF(TEXT.GT.0) WRITE (TEXT, 10002) COLUMN
PRINTJ = .FALSE.
PRINTS = .FALSE.
PRINTY = .FALSE.
ENDIF
C/// IF THE JACOBIAN IS NOT CURRENT, THEN CONTINUE THE ALGORITHM FROM
C/// THE POINT OF EVALUATING A NEW JACOBIAN AND FORMING A NEW S.
C/// OTHERWISE, THE ALGORITHM FAILS.
IF(FORMED.LT.INDEX) GO TO 0200
REPORT = QBNDS
SUCCES = .FALSE.
IF(LEVEL1 .GE. 1) THEN
IF(TEXT.GT.0) WRITE (TEXT, 10005) ID
DO 0400 J = 1, N
IF(BELOW(J).EQ.X(J) .OR. X(J).EQ.ABOVE(J)) THEN
IF(TEXT.GT.0)
+ WRITE (TEXT, 10006) J, BELOW(J), X(J), ABOVE(J)
ENDIF
0400 CONTINUE
ENDIF
GO TO 99999
ENDIF
C/// FORM THE TRIAL X, TAKING CARE THAT ROUNDING ERRORS DO NOT PLACE IT
C/// OUT OF BOUNDS. BUT IF ONE OF THE VARIABLES SHOULD LIE AT A
C/// BOUNDARY, THEN THEN PLACE IT THERE BECAUSE ROUNDING ERRORS MAY
C/// LEAVE THE VARIABLE SHORT OF THE BOUNDARY.
0500 CONTINUE
DO 5100 J = 1, N
XTRIAL(J) = X(J) - COEFF * S(J)
5100 CONTINUE
IF(REDUCE.EQ.0) THEN
DO 5200 J = 1, N
XTRIAL(J) = MIN (XTRIAL(J), ABOVE(J))
5200 CONTINUE
DO 5300 J = 1, N
XTRIAL(J) = MAX (XTRIAL(J), BELOW(J))
5300 CONTINUE
ENDIF
IF(FORCE .AND. REDUCE.EQ.0) THEN
FORCE = .FALSE.
XTRIAL(ENTRY) = VALUE
ENDIF
C/// EVALUATE THE FUNCTION AT THE TRIAL X AND STORE IN THE TRIAL Y.
ASSIGN 5400 TO RETURN
GO TO 9961
5400 CONTINUE
CALL TPNORM (N, YTNORM, YTRIAL)
C/// APPLY THE INVERSE JACOBIAN TO THE TRIAL Y AND STORE IN THE TRIAL
C/// S.
ASSIGN 5500 TO RETURN
GO TO 9971
5500 CONTINUE
CALL TPNORM (N, STNORM, STRIAL)
C///////////////////////////////////////////////////////////////////////
C
C PRINT.
C
C///////////////////////////////////////////////////////////////////////
IF(LEVEL1 .GE. 1) THEN
DO 6100 J = 1, 7
COLUMN(J) = ' '
6100 CONTINUE
IF(PRINTY) WRITE (COLUMN(1), '(I6)') INDEX
IF(PRINTY) CALL LOGSTR ('(F6.2)', 6, COLUMN(2), YNORM)
IF(PRINTS) CALL LOGSTR ('(F6.2)', 6, COLUMN(3), SNORM)
CALL LOGSTR ('(F6.2)', 6, COLUMN(4), YTNORM)
CALL LOGSTR ('(F6.2)', 6, COLUMN(5), STNORM)
IF(COEFF.LT.ONE) CALL LOGSTR ('(F6.2)', 6, COLUMN(6), COEFF)
IF(PRINTJ) CALL LOGSTR ('(F6.2)', 6, COLUMN(7), CONDIT)
IF(TEXT.GT.0) WRITE (TEXT, 10002) COLUMN
PRINTJ = .FALSE.
PRINTS = .FALSE.
PRINTY = .FALSE.
ENDIF
C///////////////////////////////////////////////////////////////////////
C
C BLOCK FOR THE DECISION TREE.
C
C///////////////////////////////////////////////////////////////////////
C/// IF THE NORM OF THE TRIAL S IS SMALLER THAN THE NORM OF THE PRESENT
C/// S, THEN THE STEP IS SUCCESSFUL. REPLACE X, Y, AND S BY THE TRIAL
C/// VALUES.
C THIS IF-BLOCK IS IMPLEMENTED USING GO-TO'S IN ORDER TO PERMIT
C THE CODE INSIDE THE BLOCK TO USE REVERSE COMMUNICATION.
C IF(STNORM.LT.SNORM) THEN
IF(.NOT. (STNORM.LT.SNORM)) GO TO 7500
INDEX = INDEX + 1
CALL COPY (N, XTRIAL, X)
CALL COPY (N, YTRIAL, Y)
CALL COPY (N, STRIAL, S)
SNORM = STNORM
YNORM = YTNORM
PRINTS = .TRUE.
PRINTY = .TRUE.
C/// IF THE SOLUTION IS ACCEPTABLE, THEN EXIT WITH SUCCESS.
ASSIGN 7100 TO RETURN
GO TO 9941
7100 CONTINUE
ACCEPT = ACCEPT .AND. (1.LE.INDEX)
IF(.NOT. ACCEPT) GO TO 7400
IF(LEVEL1 .GE. 1) THEN
C PRINT A LINE FOR THE FINAL SOLUTION.
DO 7200 J = 1, 7
COLUMN(J) = ' '
7200 CONTINUE
WRITE (COLUMN(1), '(I6)') INDEX
CALL LOGSTR ('(F6.2)', 6, COLUMN(2), YNORM)
CALL LOGSTR ('(F6.2)', 6, COLUMN(3), SNORM)
IF(PRINTJ) CALL LOGSTR ('(F6.2)', 6, COLUMN(7), CONDIT)
IF(TEXT.GT.0) WRITE (TEXT, 10002) COLUMN
C PRINT THE FINAL SOLUTION.
IF(LEVEL1 .GE. LEVEL2 .AND. LEVEL2 .GE. 1) THEN
IF(TEXT.GT.0) WRITE (TEXT, 10003) ID
ASSIGN 7300 TO RETURN
IF(TEXT.GT.0) GO TO 9951
ELSEIF(LEVEL1 .GE. 1) THEN
IF(TEXT.GT.0) WRITE (TEXT, 10004) ID
ENDIF
ENDIF
7300 CONTINUE
REPORT = QDONE
SUCCES = .TRUE.
GO TO 99999
7400 CONTINUE
C/// OTHERWISE, IF THE LIMITING NUMBER OF STEPS HAVE BEEN TAKEN, THE
C/// ALGORITHM FAILS.
IF(INDEX.EQ.STMAX) THEN
REPORT = QSMAX
SUCCES = .FALSE.
IF(LEVEL1 .GE. 1) THEN
IF(TEXT.GT.0) WRITE (TEXT, 10007) ID
ENDIF
GO TO 99999
ENDIF
C/// OTHERWISE, CONTINUE THE ALGORITHM. IF THE JACOBIAN IS YOUNG,
C/// THEN CONTINUE AT THE POINT OF SELECTING THE NEW DAMPING
C/// COEFFICIENT. BUT IF THE JACOBIAN IS OF RETIREMENT AGE, THEN
C/// CONTINUE AT THE POINT OF EVALUATING A NEW JACOBIAN.
AGE = INDEX - MAX (0, FORMED)
IF(AGE.LT.RETIRE) GO TO 0300
GO TO 0200
C/// IF THE NORM OF THE TRIAL S IS NO SMALLER THAN THE NORM OF THE
C/// PRESENT S, THEN THE STEP IS NOT SUCCESSFUL.
C ELSE
7500 CONTINUE
C/// IF THE DAMPING COEFFICIENT HAS BEEN REDUCED ONLY A FEW TIMES, THEN
C/// DECREASE THE COEFFICIENT AND CONTINUE THE ALGORITHM AT THE POINT
C/// OF FORMING A NEW TRIAL X.
COEFF = HALF * COEFF
REDUCE = REDUCE + 1
IF(REDUCE.LE.LIMIT) GO TO 0500
C/// OTHERWISE, S MUST BE BAD. IF THE JACOBIAN IS NOT CURRENT, THEN
C/// CONTINUE THE ALGORITHM FROM THE POINT OF EVALUATING A NEW
C/// JACOBIAN AND FORMING A NEW S.
IF(FORMED.LT.INDEX) GO TO 0200
C/// OTHERWISE, THE ALGORITHM FAILS.
REPORT = QDVRG
SUCCES = .FALSE.
IF(LEVEL1 .GE. 1) THEN
IF(TEXT.GT.0) WRITE (TEXT, 10008) ID
ENDIF
GO TO 99999
C ENDIF
C///////////////////////////////////////////////////////////////////////
C
C REQUEST REVERSE COMMUNICATION.
C
C///////////////////////////////////////////////////////////////////////
C/// EVALUATE THE FUNCTION AT X AND STORE IN Y.
9911 CONTINUE
CALL COPY (N, X, BUFFER)
FUNCT = .TRUE.
REENT = .TRUE.
ASSIGN 9912 TO ROUTE
GO TO 99999
9912 CONTINUE
CALL COPY (N, BUFFER, Y)
GO TO RETURN
C/// FORM THE JACOBIAN AT THE LATEST X.
9921 CONTINUE
CALL COPY (N, X, BUFFER)
JACOB = .TRUE.
REENT = .TRUE.
ASSIGN 9922 TO ROUTE
GO TO 99999
9922 CONTINUE
GO TO RETURN
C/// APPLY THE INVERSE JACOBIAN TO Y AND STORE THE STEP DIRECTION IN S.
9931 CONTINUE
CALL COPY (N, Y, BUFFER)
SOLVE = .TRUE.
REENT = .TRUE.
ASSIGN 9932 TO ROUTE
GO TO 99999
9932 CONTINUE
CALL COPY (N, BUFFER, S)
GO TO RETURN
C/// CHECK WHETHER THE NEW S IS SUFFICIENTLY SMALL THAT THE SOLUTION IS
C/// ACCEPTABLE.
9941 CONTINUE
DECIDE = .TRUE.
REENT = .TRUE.
ASSIGN 9942 TO ROUTE
GO TO 99999
9942 CONTINUE
GO TO RETURN
C/// PRINT THE SOLUTION.
9951 CONTINUE
CALL COPY (N, X, BUFFER)
SHOW = .TRUE.
REENT = .TRUE.
ASSIGN 9952 TO ROUTE
GO TO 99999
9952 CONTINUE
GO TO RETURN
C/// EVALUATE THE FUNCTION AT THE TRIAL X AND STORE IN THE TRIAL Y.
9961 CONTINUE
CALL COPY (N, XTRIAL, BUFFER)
FUNCT = .TRUE.
REENT = .TRUE.
ASSIGN 9962 TO ROUTE
GO TO 99999
9962 CONTINUE
CALL COPY (N, BUFFER, YTRIAL)
GO TO RETURN
C/// APPLY THE INVERSE JACOBIAN TO THE TRIAL Y AND STORE THE TRIAL
C/// STEP DIRECTION IN THE TRIAL S.
9971 CONTINUE
CALL COPY (N, YTRIAL, BUFFER)
SOLVE = .TRUE.
REENT = .TRUE.
ASSIGN 9972 TO ROUTE
GO TO 99999
9972 CONTINUE
CALL COPY (N, BUFFER, STRIAL)
GO TO RETURN
C///////////////////////////////////////////////////////////////////////
C
C FORMAT STATEMENTS FOR INFORMATIVE MESSAGES.
C
C///////////////////////////////////////////////////////////////////////
10001 FORMAT
+ (/1X, A9, 'LOG OF THE DAMPED AND SIMPLIFIED NEWTON ALGORITHM.'
+ /4(/10X, A45, A15)/)
10002 FORMAT
+ (10X, A6, 6(3X, A6))
10003 FORMAT
+ (/1X, A9, 'SUCCESS. THE SOLUTION:')
10004 FORMAT
+ (/1X, A9, 'SUCCESS.')
10005 FORMAT
+ (/1X, A9, 'FAILURE. SOME VARIABLES LIE AT THEIR UPPER OR'
+ /10X, 'LOWER BOUNDS. THE NEXT STEP WOULD CARRY THEM OUT OF'
+ /10X, 'BOUNDS.'
C 123456789 123456789 123456789 123456789
+ //10X, ' VARIABLE LOWER UPPER'
+ /10X, ' NUMBER BOUND VALUE BOUND'
+ /)
10006 FORMAT
+ (10X, I9, 3(3X, E9.2))
10007 FORMAT
+ (/1X, A9, 'FAILURE. THE LIMITING NUMBER OF STEPS HAS BEEN'
+ /10X, 'TAKEN.')
10008 FORMAT
+ (/1X, A9, 'FAILURE. NO DAMPING COEFFICIENT PRODUCED A TRIAL'
+ /10X, 'SOLUTION WITH A STEP VECTOR SHORTER THAN THE PRESENT'
+ /10X, 'UNDAMPED STEP.')
C///////////////////////////////////////////////////////////////////////
C
C ERROR MESSAGES.
C
C///////////////////////////////////////////////////////////////////////
9001 IF(TEXT.GT.0)
+ WRITE (TEXT, 99001) ID, N
GO TO 99999
9002 IF(TEXT.GT.0) THEN
WRITE (TEXT, 99002) ID, N, COUNT
COUNT = 0
DO 9100 J = 1, N
IF(.NOT. BELOW(J).LT.ABOVE(J)) THEN
COUNT = COUNT + 1
IF(COUNT.LE.20) THEN
WRITE (TEXT, 98001) J, BELOW(J), ABOVE(J)
ELSEIF(COUNT.EQ.21) THEN
WRITE (TEXT, 98002)
ENDIF
ENDIF
9100 CONTINUE
ENDIF
GO TO 99999
9003 IF(TEXT.GT.0)
+ WRITE (TEXT, 99003) ID, COUNT
GO TO 99999
9004 IF(TEXT.GT.0)
+ WRITE (TEXT, 99004) ID, INDEX, COEFF
GO TO 99999
99001 FORMAT
+ (/1X, A9, 'ERROR. THE NUMBER OF VARIABLES IS NOT POSITIVE.'
+ //1X, I10, ' NUMBER')
99002 FORMAT
+ (/1X, A9, 'ERROR. THE LEFT BOUNDS UPON SOME SOLUTION VALUES'
+ /10X, 'EQUAL OR EXCEED THE RIGHT BOUNDS.'
+ //1X, I10, ' SOLUTION VARIABLES'
+ /1X, I10, ' VARIABLES WITH INCONSISENT BOUNDS'
C 123456789_ 123456789_ 123456789_
+ //1X, ' INDEX LEFT BOUND RIGHT')
98001 FORMAT
+ (1X, I10, 1P, 2X, E10.2, 2X, E10.2)
98002 FORMAT
+ (1X, ' ...')
99003 FORMAT
+ (/1X, A9, 'ERROR. SOME INITIAL SOLUTION VALUES LIE OUTSIDE'
+ /10X, 'THEIR BOUNDS.'
+ //1X, I10, ' VALUES OUT OF BOUNDS')
99004 FORMAT
+ (/1X, A9, 'ERROR. THE INITIAL DAMPING COEFFICIENT IS NEGATIVE.'
+ /10X, 'THE PROGRAM LOGIC SHOULD NOT PERMIT THIS. EITHER THE'
+ /10X, 'LOGIC IS WRONG, OR THE MACHINE ARITHMETIC IS'
+ /10X, 'INCONSISTENT.'
+ //1X, I10, ' LATEST SOLUTION VECTOR'
+ /1X, E10.2, ' COEFFICIENT')
C/// EXIT.
99999 CONTINUE
C COPY THE PROTECTED LOCAL VARIABLE INTO THE ARGUMENT LIST.
STEPS = INDEX
END
c
SUBROUTINE TPNORM (N, VALUE, X)
IMPLICIT COMPLEX (A - Z)
C*****precision > double
DOUBLE PRECISION VALUE, X
C*****END precision > double
INTEGER J, N
C*****precision > single
C REAL VALUE, X
C*****END precision > single
DIMENSION X(N)
VALUE = 0.0
DO 0100 J = 1, N
VALUE = MAX (VALUE, ABS (X(J)))
0100 CONTINUE
END
c
SUBROUTINE REFINE
+ (ERROR, TEXT, ABOVE, ACTIVE, BELOW, BUFFER, COMPS, LEVEL1,
+ LEVEL2, MARK, MESH, NEWMSH, PLIMIT, PMAX, POINTS, RATIO1,
+ RATIO2, REENT, SHOW, SUCCES, TOLER0, TOLER1, TOLER2, UPDATE,
+ VARY, VARY1, VARY2, X)
IMPLICIT COMPLEX (A - Z)
CHARACTER ID*9
EXTERNAL COPY
INTEGER
+ ACT, ACTSIG, BOUND, COMPS, COUNT, COUNT1, COUNT2, EQUAL,
+ FORMER, HIGHER, IDEAL, J, K, LEVEL1, LEVEL2, MORE, NEW, OLD,
+ PLIMIT, PMAX, PMIN, POINTS, ROUTE, STRICT, TEXT, TOTAL, VARY,
+ VARY1, VARY2
LOGICAL
+ ACTIVE, ERROR, MARK, NEWMSH, REENT, SHOW, SIGNIF, SUCCES,
+ UPDATE
C*****precision > double
DOUBLE PRECISION
+ ABOVE, BELOW, BUFFER, DIFFER, HALF, LEFT, LENGTH, LOWER,
+ MAXMAG, MEAN, MESH, ONE, RANGE, RATIO1, RATIO2, RIGHT, TEMP,
+ TEMP1, TEMP2, TOLER0, TOLER1, TOLER2, UPPER, X, ZERO
C*****END precision > double
C*****precision > single
C REAL
C + ABOVE, BELOW, BUFFER, DIFFER, HALF, LEFT, LENGTH, LOWER,
C + MAXMAG, MEAN, MESH, ONE, RANGE, RATIO1, RATIO2, RIGHT, TEMP,
C + TEMP1, TEMP2, TOLER0, TOLER1, TOLER2, UPPER, X, ZERO
C*****END precision > single
PARAMETER (ID = 'REFINE: ')
PARAMETER (PMIN = 3)
DIMENSION
+ ABOVE(COMPS, PMAX), ACTIVE(COMPS), BELOW(COMPS, PMAX),
+ BUFFER(COMPS * PMAX), MARK(PMAX),
+ MESH(PMAX), RATIO1(PMAX), RATIO2(PMAX), VARY(PMAX),
+ VARY1(PMAX), VARY2(PMAX), X(COMPS, PMAX)
C/// SAVE LOCAL VARIABLES DURING RETURNS FOR REVERSE COMMUNICATION.
SAVE
C///////////////////////////////////////////////////////////////////////
C
C RETURN FROM REVERSE COMMUNICATION.
C
C///////////////////////////////////////////////////////////////////////
C TURN OFF ALL REVERSE COMMUNICATION FLAGS.
SHOW = .FALSE.
UPDATE = .FALSE.
C TURN OFF ALL COMPLETION STATUS FLAGS.
ERROR = .FALSE.
NEWMSH = .FALSE.
SUCCES = .FALSE.
C/// CONTINUE WHERE THE PROGRAM PAUSED.
IF(REENT) THEN
REENT = .FALSE.
GO TO ROUTE
ENDIF
C///////////////////////////////////////////////////////////////////////
C
C ENTER.
C
C///////////////////////////////////////////////////////////////////////
C/// ONE-TIME INITIALIZATION.
C PRECISION-INDEPENDENT CONSTANTS.
HALF = 0.5
ONE = 1.0
ZERO = 0.0
C/// LEVEL1 PRINTING.
IF(LEVEL1 .GE. 1 .AND. TEXT.GT.0) WRITE (TEXT, 10001) ID
C/// CHECK THE ARGUMENTS.
ERROR = .NOT. (0.LE.PLIMIT)
IF(ERROR) GO TO 9001
ERROR = .NOT. (PMIN.LE.POINTS .AND. POINTS.LE.PMAX)
IF(ERROR) GO TO 9002
COUNT = 0
DO 1100 J = 1, COMPS
IF(ACTIVE(J)) COUNT = COUNT + 1
1100 CONTINUE
ERROR = .NOT. (1.LE.COUNT)
IF(ERROR) GO TO 9003
ERROR = .NOT. (ZERO.LE.TOLER0)
IF(ERROR) GO TO 9004
ERROR = .NOT. (ZERO.LE.TOLER1 .AND. TOLER1.LE.ONE
+ .AND. ZERO.LE.TOLER2 .AND. TOLER2.LE.ONE)
IF(ERROR) GO TO 9005
COUNT = 0
DO 1200 K = 1, POINTS - 1
IF(MESH(K).LT.MESH(K + 1)) COUNT = COUNT + 1
1200 CONTINUE
ERROR = .NOT. (COUNT.EQ.0 .OR. COUNT.EQ.POINTS - 1)
IF(ERROR) GO TO 9006
LENGTH = ABS (MESH(POINTS) - MESH(1))
COUNT = 0
DO 1300 K = 1, POINTS - 1
DIFFER = ABS (MESH(K) - MESH(K + 1))
TEMP = LENGTH + DIFFER
IF(TEMP.EQ.LENGTH) COUNT = COUNT + 1
1300 CONTINUE
ERROR = .NOT. (COUNT.EQ.0)
IF(ERROR) GO TO 9007
C///////////////////////////////////////////////////////////////////////
C
C AT EACH INTERVAL, COUNT THE ACTIVE, SIGNIFICANT COMPONENTS THAT
C VARY TOO GREATLY.
C
C///////////////////////////////////////////////////////////////////////
ACT = 0
ACTSIG = 0
DO 2100 K = 1, POINTS
MARK(K) = .FALSE.
RATIO1(K) = ZERO
RATIO2(K) = ZERO
VARY1(K) = 0
VARY2(K) = 0
2100 CONTINUE
C/// TOP OF THE LOOP OVER THE COMPONENTS.
DO 2600 J = 1, COMPS
IF(ACTIVE(J)) THEN
ACT = ACT + 1
C/// FIND THE RANGE AND MAXIMUM MAGNITUDE OF THE COMPONENT.
LOWER = X(J, 1)
UPPER = X(J, 1)
MAXMAG = ABS (X(J, 1))
DO 2200 K = 2, POINTS
LOWER = MIN (LOWER, X(J, K))
UPPER = MAX (UPPER, X(J, K))
MAXMAG = MAX (MAXMAG, ABS (X(J, K)))
2200 CONTINUE
RANGE = UPPER - LOWER
C/// DECIDE WHETHER THE COMPONENT IS SIGNIFICANT.
SIGNIF = ABS (RANGE).GT.TOLER0 * MAX (ONE, MAXMAG)
IF(SIGNIF) THEN
ACTSIG = ACTSIG + 1
C/// AT EACH INTERVAL, SEE WHETHER THE COMPONENT'S CHANGE EXCEEDS SOME
C/// FRACTION OF THE COMPONENT'S GLOBAL CHANGE.
DO 2300 K = 1, POINTS - 1
DIFFER = ABS (X(J, K + 1) - X(J, K))
IF(RANGE.NE.ZERO)
+ RATIO1(K) = MAX (RATIO1(K), DIFFER / RANGE)
IF(DIFFER.GT.TOLER1 * RANGE)
+ VARY1(K) = VARY1(K) + 1
2300 CONTINUE
C/// FIND THE GLOBAL CHANGE OF THE COMPONENT'S DERIVATIVE.
TEMP = (X(J, 2) - X(J, 1)) / (MESH(2) - MESH(1))
LOWER = TEMP
UPPER = TEMP
DO 2400 K = 2, POINTS - 1
TEMP = (X(J, K + 1) - X(J, K))
+ / (MESH(K + 1) - MESH(K))
LOWER = MIN (LOWER, TEMP)
UPPER = MAX (UPPER, TEMP)
2400 CONTINUE
RANGE = UPPER - LOWER
C/// AT EACH INTERIOR POINT, SEE WHETHER THE DERIVATIVE'S CHANGE
C/// EXCEEDS SOME FRACTION OF THE DERIVATIVE'S GLOBAL CHANGE.
DO 2500 K = 2, POINTS - 1
LEFT = (X(J, K) - X(J, K - 1))
+ / (MESH(K) - MESH(K - 1))
RIGHT = (X(J, K + 1) - X(J, K))
+ / (MESH(K + 1) - MESH(K))
DIFFER = ABS (LEFT - RIGHT)
IF(RANGE.NE.ZERO)
+ RATIO2(K) = MAX (RATIO2(K), DIFFER / RANGE)
IF(DIFFER.GT.TOLER2 * RANGE)
+ VARY2(K) = VARY2(K) + 1
2500 CONTINUE
C/// BOTTOM OF THE LOOP OVER THE COMPONENTS.
ENDIF
ENDIF
2600 CONTINUE
C/// COUNT THE INTERVALS IN WHICH VARIATIONS THAT ARE TOO LARGE OCCUR.
IDEAL = 0
DO 2700 K = 1, POINTS - 1
VARY(K) = VARY1(K)
IF(1.LT.K) VARY(K) = VARY(K) + VARY2(K)
IF(K.LT.POINTS - 1) VARY(K) = VARY(K) + VARY2(K + 1)
IF(0.LT.VARY(K)) IDEAL = IDEAL + 1
2700 CONTINUE
C///////////////////////////////////////////////////////////////////////
C
C SELECT THE INTERVALS TO HALVE.
C
C HALVE INTERVALS ON WHICH COMPONENTS OR DERIVATIVES VARY TOO
C GREATLY. IF IT IS NOT POSSIBLE TO HALVE ALL SUCH INTERVALS, THEN
C FIND A LOWER BOUND FOR THE NUMBER OF VARYING COMPONENTS PER
C INTERVAL SO THAT ALL INTERVALS WITH STRICTLY MORE CAN BE HALVED.
C
C///////////////////////////////////////////////////////////////////////
C/// FIND THE BOUND.
BOUND = 0
MORE = MAX (0, MIN (PMAX - POINTS, PLIMIT))
3100 CONTINUE
EQUAL = 0
STRICT = 0
DO 3200 K = 1, POINTS - 1
IF(VARY(K).EQ.BOUND) THEN
EQUAL = EQUAL + 1
ELSEIF(VARY(K).GT.BOUND) THEN
STRICT = STRICT + 1
IF(STRICT.EQ.1) THEN
HIGHER = VARY(K)
ELSE
HIGHER = MIN (HIGHER, VARY(K))
ENDIF
ENDIF
3200 CONTINUE
IF(STRICT.GT.MORE) THEN
BOUND = HIGHER
GO TO 3100
ENDIF
C/// DETERMINE HOW MANY INTERVALS TO HALVE OF THOSE WHOSE NUMBER OF
C/// VARIATIONS EXACTLY EQUAL THE BOUND.
C IF(0.LT.BOUND .AND. 0.EQ.STRICT) THEN
C EQUAL = MIN (EQUAL, MORE)
C ELSE
C EQUAL = 0
C ENDIF
IF(0.LT.BOUND) THEN
EQUAL = MIN (EQUAL, MORE - STRICT)
ELSE
EQUAL = 0
ENDIF
C/// MARK THE INTERVALS TO HALVE.
COUNT = 0
DO 3300 K = 1, POINTS - 1
IF(BOUND.LT.VARY(K)) THEN
MARK(K) = .TRUE.
ELSEIF(BOUND.EQ.VARY(K) .AND. COUNT.LT.EQUAL) THEN
COUNT = COUNT + 1
MARK(K) = .TRUE.
ELSE
MARK(K) = .FALSE.
ENDIF
3300 CONTINUE
C///////////////////////////////////////////////////////////////////////
C
C HALVE THE INTERVALS, IF ANY.
C
C///////////////////////////////////////////////////////////////////////
C/// FORM THE TOTAL NUMBER OF POINTS IN THE NEW MESH.
TOTAL = POINTS + EQUAL + STRICT
IF(0.EQ.EQUAL + STRICT) GO TO 4700
C/// TOP OF THE BLOCK TO CREATE THE NEW MESH. CHECK THE INTERVALS.
COUNT1 = 0
COUNT2 = 0
LENGTH = ABS (MESH(POINTS) - MESH(1))
DO 4100 K = 1, POINTS - 1
IF(MARK(K)) THEN
MEAN = HALF * (MESH(K) + MESH(K + 1))
TEMP = LENGTH + ABS (MESH(K + 1) - MEAN)
IF(TEMP.EQ.LENGTH) COUNT1 = COUNT1 + 1
TEMP = LENGTH + ABS (MEAN - MESH(K))
IF(TEMP.EQ.LENGTH) COUNT1 = COUNT1 + 1
IF(.NOT. ((MESH(K).LT.MEAN .AND. MEAN.LT.MESH(K + 1))
+ .OR. (MESH(K + 1).LT.MEAN .AND. MEAN.LT.MESH(K))))
+ COUNT2 = COUNT2 + 1
ENDIF
4100 CONTINUE
ERROR = .NOT. (COUNT1.EQ.0 .AND. COUNT2.EQ.0)
IF(ERROR) GO TO 9008
C/// ADD THE NEW POINTS. INTERPOLATE X AND THE BOUNDS.
NEW = TOTAL
DO 4400 OLD = POINTS, 2, - 1
MESH(NEW) = MESH(OLD)
DO 4200 J = 1, COMPS
ABOVE(J, NEW) = ABOVE(J, OLD)
BELOW(J, NEW) = BELOW(J, OLD)
X(J, NEW) = X(J, OLD)
4200 CONTINUE
NEW = NEW - 1
IF(MARK(OLD - 1)) THEN
MESH(NEW) = HALF * (MESH(OLD) + MESH(OLD - 1))
DO 4300 J = 1, COMPS
ABOVE(J, NEW)
+ = HALF * (ABOVE(J, OLD) + ABOVE(J, OLD - 1))
BELOW(J, NEW)
+ = HALF * (BELOW(J, OLD) + BELOW(J, OLD - 1))
X(J, NEW) = HALF * (X(J, OLD) + X(J, OLD - 1))
4300 CONTINUE
NEW = NEW - 1
ENDIF
4400 CONTINUE
C/// MARK THE NEW POINTS.
NEW = TOTAL
DO 4500 OLD = POINTS, 2, - 1
MARK(NEW) = .FALSE.
NEW = NEW - 1
IF(MARK(OLD - 1)) THEN
MARK(NEW) = .TRUE.
NEW = NEW - 1
ENDIF
4500 CONTINUE
MARK(NEW) = .FALSE.
C/// UPDATE THE NUMBER OF POINTS.
FORMER = POINTS
POINTS = TOTAL
C/// ALLOW THE USER TO UPDATE THE SOLUTION AND BOUNDS.
CALL COPY (COMPS * POINTS, X, BUFFER)
REENT = .TRUE.
UPDATE = .TRUE.
ASSIGN 4600 TO ROUTE
GO TO 99999
4600 CONTINUE
CALL COPY (COMPS * POINTS, BUFFER, X)
C/// BOTTOM OF THE BLOCK TO CREATE A NEW MESH.
4700 CONTINUE
C///////////////////////////////////////////////////////////////////////
C
C PRINT.
C
C///////////////////////////////////////////////////////////////////////
C/// LEVEL1 PRINTING. SHOW THE COUNTS AND TOLERANCES.
IF(LEVEL1 .GE. 1 .AND. TEXT.GT.0) THEN
WRITE (TEXT, 10002) COMPS - ACT, ACT - ACTSIG
IF(ACTSIG.EQ.0) THEN
WRITE (TEXT, 10003) ID
ELSE
TEMP1 = RATIO1(1)
DO 5100 K = 2, FORMER - 1
TEMP1 = MAX (TEMP1, RATIO1(K))
5100 CONTINUE
TEMP2 = RATIO2(2)
DO 5200 K = 3, FORMER - 1
TEMP2 = MAX (TEMP2, RATIO2(K))
5200 CONTINUE
WRITE (TEXT, 10004)
+ TEMP1, TOLER1, TEMP2, TOLER2, IDEAL, EQUAL + STRICT
IF(IDEAL.GT.0 .AND. MORE.EQ.0) WRITE (TEXT, 10005)
C/// PRINT THE MESH.
IF(IDEAL.EQ.0) THEN
WRITE (TEXT, 10006)
ELSE
WRITE (TEXT, 10007)
ENDIF
WRITE (TEXT, 10008)
OLD = 0
DO 5300 K = 1, POINTS
IF(.NOT. MARK(K)) THEN
OLD = OLD + 1
IF(1.LT.K) THEN
IF(MARK(K - 1)) THEN
WRITE (TEXT, 10009)
+ K - 1, MESH(K - 1),
+ RATIO1(OLD - 1), VARY1(OLD - 1)
ELSE
WRITE (TEXT, 10010)
+ RATIO1(OLD - 1), VARY1(OLD - 1)
ENDIF
ENDIF
IF(1.LT.K .AND. K.LT.POINTS) THEN
WRITE (TEXT, 10011)
+ K, MESH(K), RATIO2(OLD), VARY2(OLD)
ELSE
WRITE (TEXT, 10011) K, MESH(K)
ENDIF
ENDIF
5300 CONTINUE
ENDIF
ENDIF
C/// LEVEL2 PRINTING.
IF(1.LE.LEVEL2 .AND. LEVEL2.LE.LEVEL1 .AND.
+ 0.LT.EQUAL + STRICT .AND. TEXT.GT.0) THEN
WRITE (TEXT, 10012) ID
CALL COPY (COMPS * POINTS, X, BUFFER)
SHOW = .TRUE.
REENT = .TRUE.
ASSIGN 5400 TO ROUTE
GO TO 99999
ENDIF
5400 CONTINUE
C///////////////////////////////////////////////////////////////////////
C
C DEFINE COMPLETION STATUS FLAGS.
C
C///////////////////////////////////////////////////////////////////////
SUCCES = 0.EQ.IDEAL
NEWMSH = 0.LT.EQUAL + STRICT
C///////////////////////////////////////////////////////////////////////
C
C INFORMATIVE MESSAGES.
C
C///////////////////////////////////////////////////////////////////////
10001 FORMAT
+ (/1X, A9, 'REFINE THE MESH, IF NECESSARY.')
10002 FORMAT
+ (/10X, I10, ' INACTIVE COMPONENTS'
+ /10X, I10, ' INSIGNIFICANT, ACTIVE COMPONENTS')
10003 FORMAT
+ (/1X, A9, 'WARNING. THERE ARE NO ACTIVE, SIGNIFICANT'
+ /10X, 'COMPONENTS.')
10004 FORMAT
+ (/10X, F10.5, ' LARGEST CHANGE IN VALUE (LOCAL / GLOBAL)'
+ /10X, F10.5, ' DESIRED BOUND'
+ //10X, F10.5, ' LARGEST CHANGE IN DERIVATIVE (LOCAL / GLOBAL)'
+ /10X, F10.5, ' DESIRED BOUND'
+ //10X, I10, ' POTENTIAL NUMBER OF NEW MESH POINTS'
+ /10X, I10, ' ACTUAL NUMBER')
10005 FORMAT
+ (/10X, 'WARNING. MORE POINTS WOULD BE ADDED TO THE MESH WERE'
+ /10X, 'IT NOT FOR THE LIMIT ON NEW POINTS OR THE LIMIT ON'
+ /10X, 'TOTAL POINTS.')
10006 FORMAT
+ (/10X, 'NO NEW POINTS ARE NEEDED. THE MESH:')
10007 FORMAT
+ (/10X, 'THE MESH, "*" DENOTES A NEW POINT:')
10008 FORMAT
C 123 123 1 123 1
C 123456 1234567890123456 1234 1234 1234 1234
+ (/10X, ' LARGEST NUMBER OF'
+ /10X, ' VARIATION TOO LARGE'
+ /10X, ' INDEX MESH POINT RATIO RATIOS'
+ /10X, '------ ---------------- --------- ---------'
+ /10X, ' 1ST 2ND 1ST 2ND')
10009 FORMAT
+ (10X, I6, '* ', E16.9, 3X, F4.2, 8X, I4)
10010 FORMAT
+ (38X, F4.2, 8X, I4)
10011 FORMAT
+ (10X, I6, ' ', E16.9, 8X, F4.2, 8X, I4)
10012 FORMAT
+ (/1X, A9, 'THE SOLUTION ESTIMATE ON THE NEW MESH:')
C///////////////////////////////////////////////////////////////////////
C
C ERROR MESSAGES.
C
C///////////////////////////////////////////////////////////////////////
GO TO 99999
9001 IF(TEXT.GT.0) WRITE (TEXT, 99001)
+ ID, PLIMIT
GO TO 99999
9002 IF(TEXT.GT.0) WRITE (TEXT, 99002)
+ ID, POINTS, PMIN, PMAX
GO TO 99999
9003 IF(TEXT.GT.0) WRITE (TEXT, 99003)
+ ID, COUNT, COMPS
GO TO 99999
9004 IF(TEXT.GT.0) WRITE (TEXT, 99004)
+ ID, TOLER0
GO TO 99999
9005 IF(TEXT.GT.0) WRITE (TEXT, 99005)
+ ID, TOLER1, TOLER2
GO TO 99999
9006 IF(TEXT.GT.0) WRITE (TEXT, 99006)
+ ID, COUNT, POINTS - 1 - COUNT
GO TO 99999
9007 IF(TEXT.GT.0) WRITE (TEXT, 99007)
+ ID, LENGTH, COUNT
GO TO 99999
9008 IF(TEXT.GT.0) WRITE (TEXT, 99008) ID, COUNT1, COUNT2
GO TO 99999
99001 FORMAT
+ (/1X, A9, 'ERROR. THE LIMIT ON THE NUMBER OF NEW POINTS MUST'
+ /10X, 'BE POSITIVE.'
+ //10X, I10, ' LIMIT')
99002 FORMAT
+ (/1X, A9, 'ERROR. THERE ARE TOO MANY OR TOO FEW POINTS.'
+ //10X, I10, ' POINTS'
+ /10X, I10, ' MINIMUM (SET BY THIS SUBROUTINE)'
+ /10X, I10, ' MAXIMUM')
99003 FORMAT
+ (/1X, A9, 'ERROR. THERE ARE NO ACTIVE COMPONENTS.'
+ //10X, I10, ' ACTIVE COMPONENTS'
+ /10X, I10, ' TOTAL COMPONENTS')
99004 FORMAT
+ (/1X, A9, 'ERROR. THE SIGNIFICANCE LEVEL IS NEGATIVE.'
+ //10X, E10.2, ' TOLER0')
99005 FORMAT
+ (/1X, A9, 'ERROR. THE BOUNDS UPON THE RELATIVE VARIATIONS'
+ /10X, 'DO NOT LIE BETWEEN ZERO AND ONE, INCLUSIVE.'
+ //10X, E10.2, ' TOLER1'
+ /10X, E10.2, ' TOLER2')
99006 FORMAT
+ (/1X, A9, 'ERROR. THE POINTS ARE NOT ORDERED.'
+ //10X, I10, ' INTERVALS WITH ENDPOINTS ORDERED LEFT TO RIGHT'
+ /10X, I10, ' INTERVALS WITH ENDPOINTS ORDERED RIGHT TO LEFT')
99007 FORMAT
+ (/1X, A9, 'ERROR. SOME INTERVALS ARE SMALLER THAN MACHINE'
+ /10X, 'EPSILON TIMES THE TOTAL LENGTH OF ALL INTERVALS.'
+ //10X, E10.2, ' TOTAL LENGTH'
+ /10X, I10, ' SMALL INTERVALS')
99008 FORMAT
+ (/1X, A9, 'ERROR. EITHER (1) SOME NEW INTERVALS WOULD BE'
+ /10X, 'SMALLER THAN MACHINE EPSILON TIMES THE TOTAL LENGTH OF'
+ /10X, 'ALL INTERVALS, OR (2) ROUNDING ERRORS PLACE SOME NEW'
+ /10X, 'POINTS OUTSIDE THE INTERVALS THEY SHOULD HALVE.'
+ //10X, I10, ' NEW INTERVALS TOO SMALL'
+ /10X, I10, ' NEW POINTS OUT OF ORDER')
99999 CONTINUE
END
c
SUBROUTINE TIMSTP
+ (ERROR, TEXT, ABOVE, ACCEPT, BELOW, BUFFER, CONDIT, DECIDE,
+ DESIRE, FUNCT, JACOB, LEVEL1, LEVEL2, N, REENT, REPORT,
+ RETIRE, S, SHOW, SIZE1, SOLVE, STEP, STMAX, STORE, STRIAL,
+ SUCCES, TDAGE, TDEC, TIME, TINC, TMAX, TMIN, TOLER, X, XSAVE,
+ XTRIAL, Y, YNORM, YTRIAL)
IMPLICIT COMPLEX (A - P, R - Z), INTEGER (Q)
CHARACTER COLUMN*16, HEADER*80, ID*9, REMARK*80, STRING*80
EXTERNAL COPY, EXTENT, LOGSTR, NEWTON
INTEGER
+ AGE, BEGIN, COUNT, DESIRE, J, K, LENGTH, LEVEL1, LEVEL2, N,
+ NUMBER, REPORT, RETIRE, ROUTE, STEP, STMAX, TDAGE, TEXT,
+ TRIALS, TRMAX, XREPOR
LOGICAL
+ ACCEPT, AGAIN, DECIDE, ERROR, EXIST, FUNCT, JACOB, REENT,
+ REENT1, SHOW, SOLVE, STORE, SUCCES, SUCCE1, TIME
C*****precision > double
DOUBLE PRECISION
C*****END precision > double
C*****precision > single
C REAL
C*****END precision > single
+ ABOVE, BELOW, BUFFER, CHANGE, CONDIT, CSAVE, DUMMY, S, SIZE0,
+ SIZE1, STRIAL, TDEC, TINC, TMAX, TMIN, TOLER, X, XSAVE, XTRIAL,
+ Y, YNORM, YTRIAL
PARAMETER (ID = 'TIMSTP: ')
PARAMETER (TRMAX = 10)
C TIME STEP SIZE STATUS CODES
PARAMETER
+ (QBAD = 1, QDEC = 2, QGOOD = 3, QNONE = 4, QNOCH = 5, QINC = 6)
C REPORT CODES
PARAMETER
+ (QNULL = 0, QDONE = 1, QDVRG = 2, QBNDS = 3, QSTDY = 4,
+ QSMAX = 5)
DIMENSION
+ ABOVE(N), BELOW(N), BUFFER(N), COLUMN(7), HEADER(5, 2), S(N),
+ STRIAL(N), X(N), XSAVE(N), XTRIAL(N), Y(N), YTRIAL(N)
C/// SAVE LOCAL VALUES DURING RETURNS FOR REVERSE COMMUNCIATION.
SAVE
C///////////////////////////////////////////////////////////////////////
C
C PROLOGUE.
C
C///////////////////////////////////////////////////////////////////////
C/// INITIALIZE.
C TURN OFF ALL REVERSE COMMUNICATION FLAGS.
DECIDE = .FALSE.
FUNCT = .FALSE.
JACOB = .FALSE.
SHOW = .FALSE.
SOLVE = .FALSE.
STORE = .FALSE.
TIME = .FALSE.
C TURN OFF ALL COMPLETION STATUS FLAGS.
ERROR = .FALSE.
REPORT = QNULL
SUCCES = .FALSE.
C/// IF THIS IS A RETURN CALL, THEN CONTINUE WHERE THE PROGRAM PAUSED.
IF(REENT) THEN
REENT = .FALSE.
GO TO ROUTE
ENDIF
C/// CHECK THE ARGUMENTS.
ERROR = .NOT. (0.LT.DESIRE)
IF(ERROR) GO TO 9001
ERROR = .NOT. (1.0.LE.TDEC .AND. 1.0.LE.TINC)
IF(ERROR) GO TO 9002
ERROR = .NOT. (0.LT.RETIRE)
IF(ERROR) GO TO 9003
ERROR = .NOT. (0.LE.STEP)
IF(ERROR) GO TO 9004
C/// LEVEL 1 AND LEVEL 2 PRINTING.
C 123456789_123456789_123456789_123456
C 123456 123456 123456 123456
HEADER(1, 1) = ' LOG10 MAX NORM '
HEADER(2, 1) = ' LOG10 --------------- '
HEADER(3, 1) = ' TIME TIME STEADY CHANGE '
HEADER(4, 1) = ' STEP STEP STATE TO THE '
HEADER(5, 1) = 'NUMBER SIZE RESIDL SOLUTN '
C 123456789_123456789_12345
C 1234 1234 123456
HEADER(1, 2) = 'NEWTON ALGORITHM '
HEADER(2, 2) = '---------------- '
HEADER(3, 2) = 'STEPS '
HEADER(4, 2) = ' JACOBIANS '
HEADER(5, 2) = ' CONDITION REMARK'
IF(0.LT.TEXT .AND. 1.EQ.LEVEL1) THEN
WRITE (TEXT, 10001) ID, ((HEADER(J, K), K = 1, 2), J = 1, 5)
ELSEIF(0.LT.TEXT .AND. 1.LT.LEVEL1) THEN
WRITE (TEXT, 20001) ID
ENDIF
C/// EVALUATE THE STEADY STATE RESIDUAL AT THE INITIAL SOLUTION.
ASSIGN 1100 TO ROUTE
GO TO 7011
1100 CONTINUE
CALL TPNORM (N, YNORM, BUFFER)
C/// LEVEL 1 AND LEVEL 2 PRINTING.
IF(0.LT.TEXT .AND. 1.EQ.LEVEL1) THEN
WRITE (COLUMN(1), '(I6)') STEP
COLUMN(2) = ' '
CALL LOGSTR ('(F6.2)', 6, COLUMN(3), YNORM)
WRITE (TEXT, 10002) (COLUMN(J), J = 1, 3)
ELSEIF(0.LT.TEXT .AND. 1.LT.LEVEL1) THEN
CALL LOGSTR ('(F10.2)', 10, STRING, YNORM)
WRITE (TEXT, 20002) STRING
ENDIF
C///////////////////////////////////////////////////////////////////////
C
C TOP OF THE LOOP OVER THE TIME STEPS.
C
C///////////////////////////////////////////////////////////////////////
C/// SAVE THE LATEST SOLUTION SHOULD NEWTON FAIL; STORE IT FOR USE BY
C/// THE FUNCTION.
CALL COPY (N, X, XSAVE)
ASSIGN 1200 TO ROUTE
GO TO 7021
1200 CONTINUE
C/// TOP OF THE LOOP.
BEGIN = STEP
IF(STEP.EQ.0) THEN
AGE = 0
QMOVE = QNONE
QSTAT = QNONE
TRIALS = 0
ENDIF
1300 CONTINUE
C///////////////////////////////////////////////////////////////////////
C
C USE NEWTON TO SOLVE THE TIME DIFFERENCED EQUATIONS.
C
C///////////////////////////////////////////////////////////////////////
C/// LEVEL 2 PRINTING.
IF(0.LT.TEXT .AND. 1.LT.LEVEL1) THEN
CALL LOGSTR ('(F10.2)', 10, STRING, SIZE1)
WRITE (TEXT, 20003) STRING
ENDIF
C/// CALL NEWTON.
COUNT = 0
EXIST = (0.LT.AGE) .AND. (BEGIN.LT.STEP)
REENT1 = .FALSE.
2100 CONTINUE
C SUBROUTINE NEWTON
C + (ERROR, TEXT, ABOVE, ACCEPT, BELOW, BUFFER, CONDIT, DECIDE,
C + EXIST, FUNCT, JACOB, LEVEL1, LEVEL2, N, REENT, REPORT, RETIRE,
C + S, SHOW, SOLVE, STEPS, STMAX, STRIAL, SUCCES, X, XTRIAL, Y,
C + YNORM, YTRIAL)
CALL NEWTON
+ (ERROR, TEXT, ABOVE, ACCEPT, BELOW, BUFFER, CONDIT, DECIDE,
+ EXIST, FUNCT, JACOB, LEVEL1 - 1, LEVEL2 - 1, N, REENT1, XREPOR,
+ TDAGE, S, SHOW, SOLVE, NUMBER, STMAX, STRIAL, SUCCE1, X,
+ XTRIAL, Y, DUMMY, YTRIAL)
IF(ERROR) GO TO 9006
C/// PASS MOST NEWTON REQUESTS TO THE CALLER. TREAT JACOBIAN REQUESTS
C/// SEPARATELY TO RETAIN THE CONDITION ESTIMATE AND THE COUNT OF
C/// JACOBIANS.
IF(REENT1) THEN
TIME = .TRUE.
IF(DECIDE .AND. 0.EQ.NUMBER) THEN
C FORCE NEWTON TO TAKE AT LEAST ONE STEP
ACCEPT = .FALSE.
GO TO 2100
ELSEIF(JACOB) THEN
GO TO 7031
ELSE
GO TO 7041
ENDIF
ENDIF
C/// DETERMINE THE CHANGE TO X AND EVALUATE THE STEADY STATE RESIDUAL.
IF(.NOT. SUCCE1) GO TO 2400
DO 2200 J = 1, N
BUFFER(J) = X(J) - XSAVE(J)
2200 CONTINUE
CALL TPNORM (N, CHANGE, BUFFER)
ASSIGN 2300 TO ROUTE
GO TO 7011
2300 CONTINUE
CALL TPNORM (N, YNORM, BUFFER)
2400 CONTINUE
C/// DETERMINE THE STATUS OF THE TIME STEP.
IF(SUCCE1) THEN
IF(CHANGE.EQ.0.0) THEN
QSTAT = QNOCH
ELSE
QSTAT = QGOOD
AGE = AGE + 1
STEP = STEP + 1
ENDIF
ELSE
QSTAT = QBAD
ENDIF
C/// IF NEWTON FAILED, THEN RESTORE THE SOLUTION.
IF(.NOT. SUCCE1) THEN
CALL COPY (N, XSAVE, X)
C/// OTHERWISE SAVE THE LATEST SOLUTION SHOULD NEWTON FAIL; STORE IT
C/// FOR USE BY THE FUNCTION.
ELSE
CALL COPY (N, X, XSAVE)
ASSIGN 2500 TO ROUTE
GO TO 7021
ENDIF
2500 CONTINUE
C///////////////////////////////////////////////////////////////////////
C
C REPORT THE OUTCOME OF THE TIME STEP.
C
C///////////////////////////////////////////////////////////////////////
C/// LEVEL 1 PRINTING.
IF(0.LT.TEXT .AND. 1.EQ.LEVEL1) THEN
DO 3100 J = 1, 7
COLUMN(J) = ' '
3100 CONTINUE
C COLUMN 1: NUMBER OF THE TIME STEP
IF(QSTAT.EQ.QGOOD) WRITE (COLUMN(1), '(I6)') STEP
C COLUMN 2: SIZE OF THE TIME STEP
WRITE (COLUMN(2), '(F6.2)') LOG10 (SIZE1)
C COLUMN 3: NORM OF THE STEADY STATE RESIDUAL
IF(QSTAT.EQ.QGOOD .OR. QSTAT.EQ.QNOCH)
+ CALL LOGSTR ('(F6.2)', 6, COLUMN(3), YNORM)
C COLUMN 4: NORM OF THE CHANGE TO THE SOLUTION
IF(QSTAT.EQ.QGOOD) CALL LOGSTR ('(F6.2)', 6, COLUMN(4), CHANGE)
C COLUMN 5: NUMBER OF NEWTON STEPS
WRITE (COLUMN(5), '(I4)') NUMBER
C COLUMN 6: NUMBER OF JACOBIANS
IF(0.LT.COUNT) WRITE (COLUMN(6), '(I4)') COUNT
C COLUMN 7: MAXIMUM CONDITION NUMBER
IF(0.LT.COUNT) CALL LOGSTR ('(F6.2)', 6, COLUMN(7), CSAVE)
C REMARK:
IF(XREPOR.EQ.QBNDS) THEN
REMARK = 'AT BOUNDS'
ELSEIF(XREPOR.EQ.QDVRG) THEN
REMARK = 'DIVERGING'
ELSEIF(XREPOR.EQ.QSMAX) THEN
REMARK = 'STEP LIMIT'
ELSEIF(QSTAT.EQ.QNOCH) THEN
REMARK = 'NO CHANGE'
ELSE
REMARK = ' '
ENDIF
CALL EXTENT (LENGTH, REMARK)
WRITE (TEXT, 10002) COLUMN, REMARK (1 : LENGTH)
C/// LEVEL 2 PRINTING.
ELSEIF(0.LT.TEXT .AND. 1.LT.LEVEL1) THEN
IF(QSTAT.EQ.QGOOD .OR. QSTAT.EQ.QNOCH) THEN
CALL LOGSTR ('(F10.2)', 10, STRING, YNORM)
WRITE (TEXT, 20004) ID, STEP, STRING
ELSE
WRITE (TEXT, 20005) ID
ENDIF
ENDIF
C///////////////////////////////////////////////////////////////////////
C
C DECISION TREE.
C
C///////////////////////////////////////////////////////////////////////
C/// DECIDE WHETHER TO CONTINUE.
IF(QSTAT.EQ.QGOOD) THEN
AGAIN = TOLER.LT.YNORM
ELSE
AGAIN = TRIALS.LT.TRMAX
ENDIF
IF(AGAIN) THEN
C/// CONTINUE WITH THE TIME STEP.
IF(QSTAT.EQ.QGOOD .AND. (AGE.LT.RETIRE .OR. TINC.EQ.1.0))
+ THEN
Q0 = QNONE
QMOVE = QNONE
TRIALS = 0
C/// CHANGE THE TIME STEP.
ELSE
TRIALS = TRIALS + 1
C/// INCREASE THE TIME STEP BY THE SPECIFIED AMOUNT.
IF((QSTAT.EQ.QGOOD .AND. AGE.EQ.RETIRE) .OR.
+ QMOVE.EQ.QNONE .AND. QSTAT.EQ.QNOCH .AND. 1.0.LT.TINC)
+ THEN
IF(0.LT.TEXT .AND. 1.LT.LEVEL1) WRITE (TEXT, 20006)
AGE = 0
QMOVE = QINC
SIZE0 = SIZE1
SIZE1 = MIN (TMAX, SIZE1 * TINC)
C/// REDUCE THE INCREASE BY THE SQUARE ROOT.
ELSEIF(QMOVE.EQ.QINC .AND. QSTAT.EQ.QBAD) THEN
IF(0.LT.TEXT .AND. 1.LT.LEVEL1) WRITE (TEXT, 20007)
AGE = 0
SIZE1 = SQRT (SIZE0 * SIZE1)
C/// REDUCE THE TIME STEP BY THE SPECIFIED AMOUNT.
ELSEIF((QMOVE.EQ.QDEC .OR. QMOVE.EQ.QNONE) .AND.
+ QSTAT.EQ.QBAD .AND. 1.0.LT.TDEC) THEN
IF(0.LT.TEXT .AND. 1.LT.LEVEL1) WRITE (TEXT, 20007)
AGE = 0
QMOVE = QDEC
SIZE0 = SIZE1
SIZE1 = MAX (TMIN, SIZE1 / TDEC)
C/// REDUCE THE REDUCTION BY THE SQUARE ROOT.
ELSEIF(QMOVE.EQ.QDEC .AND. QSTAT.EQ.QNOCH) THEN
IF(0.LT.TEXT .AND. 1.LT.LEVEL1) WRITE (TEXT, 20006)
AGE = 0
SIZE1 = SQRT (SIZE0 * SIZE1)
C/// BOTTOM OF THE TIME STEP SELECTION TREE.
ELSE
ERROR = .TRUE.
GO TO 9005
ENDIF
ENDIF
IF(STEP.LT.BEGIN + DESIRE) GO TO 1300
ENDIF
C///////////////////////////////////////////////////////////////////////
C
C EPILOGUE.
C
C///////////////////////////////////////////////////////////////////////
C/// LEVEL 1 PRINTING.
IF(0.LT.TEXT .AND. 1.EQ.LEVEL1) THEN
DO 3200 J = 1, 7
COLUMN(J) = ' '
3200 CONTINUE
C COLUMN 1: NUMBER OF THE TIME STEP
WRITE (COLUMN(1), '(I6)') BEGIN + DESIRE
C COLUMN 3: NORM OF THE STEADY STATE RESIDUAL
CALL LOGSTR ('(F6.2)', 6, COLUMN(3), TOLER)
REMARK = 'LIMITS'
CALL EXTENT (LENGTH, REMARK)
WRITE (TEXT, 10002) COLUMN, REMARK (1 : LENGTH)
C/// LEVEL 2 PRINTING.
ELSEIF(0.LT.TEXT .AND. 1.LT.LEVEL1) THEN
IF(QSTAT.EQ.QGOOD) THEN
IF(YNORM.LE.TOLER) THEN
WRITE (TEXT, 20008)
ELSE
WRITE (TEXT, 20009)
ENDIF
ELSEIF(QMOVE.EQ.QNONE) THEN
IF(QSTAT.EQ.QBAD .AND. TDEC.LE.1.0) THEN
WRITE (TEXT, 20010)
ELSEIF(QSTAT.EQ.QNOCH .AND. TINC.EQ.1.0) THEN
WRITE (TEXT, 20011)
ENDIF
ELSE
WRITE (TEXT, 20012)
ENDIF
ENDIF
C/// PRINT THE SOLUTION.
IF(0.LT.TEXT .AND. 1.EQ.LEVEL2 .AND. LEVEL2.LE.LEVEL1) THEN
WRITE (TEXT, 20013)
GO TO 7051
ENDIF
3300 CONTINUE
C/// SET THE COMPLETION STATUS FLAGS.
IF(STEP.EQ.BEGIN + DESIRE) THEN
SUCCES = .TRUE.
IF(YNORM.LE.TOLER) THEN
REPORT = QSTDY
ELSE
REPORT = QDONE
ENDIF
ELSEIF(BEGIN.LT.STEP) THEN
SUCCES = .TRUE.
REPORT = XREPOR
ELSE
SUCCES = .FALSE.
REPORT = XREPOR
ENDIF
GO TO 99999
C///////////////////////////////////////////////////////////////////////
C
C REVERSE COMMUNICATION BLOCKS.
C
C///////////////////////////////////////////////////////////////////////
C/// EVALUATE THE STEADY STATE RESIDUAL.
7011 CONTINUE
CALL COPY (N, X, BUFFER)
REENT = .TRUE.
FUNCT = .TRUE.
TIME = .FALSE.
GO TO 99999
C/// STORE THE LATEST SOLUTION FOR USE BY THE FUNCTION.
7021 CONTINUE
CALL COPY (N, X, BUFFER)
REENT = .TRUE.
STORE = .TRUE.
GO TO 99999
C/// PASS REQUESTS FOR JACOBIANS TO THE CALLER.
7031 CONTINUE
REENT = .TRUE.
ASSIGN 7032 TO ROUTE
GO TO 99999
7032 CONTINUE
IF(COUNT.EQ.0) THEN
CSAVE = CONDIT
ELSE
CSAVE = MAX (CONDIT, CSAVE)
ENDIF
COUNT = COUNT + 1
GO TO 2100
C/// PASS NEWTON REQUESTS TO THE CALLER.
7041 CONTINUE
REENT = .TRUE.
ASSIGN 2100 TO ROUTE
GO TO 99999
C/// PRINT THE SOLUTION.
7051 CONTINUE
CALL COPY (N, X, BUFFER)
REENT = .TRUE.
SHOW = .TRUE.
ASSIGN 3300 TO ROUTE
GO TO 99999
C///////////////////////////////////////////////////////////////////////
C
C INFORMATIVE MESSAGES.
C
C///////////////////////////////////////////////////////////////////////
10001 FORMAT
+ (/1X, A9, 'PERFORM A TIME INTEGRATION.'
+ /5(/10X, A36, A25)/)
10002 FORMAT
+ (10X, A6, 3X, A6, 3X, A6, 3X, A6, 3X, A4, 1X, A4, 1X, A6, 3X, A)
20001 FORMAT
+ (/1X, A9, 'BEGIN A TIME INTEGRATION.')
20002 FORMAT
+ (/10X, A10, ' LOG10 MAX NORM STEADY STATE RESIDUAL')
20006 FORMAT
+ (/10X, 'INCREASING THE TIME STEP SIZE.')
20007 FORMAT
+ (/10X, 'DECREASING THE TIME STEP SIZE.')
20003 FORMAT
+ (/10X, 'CALLING NEWTON TO ATTEMPT A TIME STEP.'
+ //10X, A10, ' LOG10 TIME STEP SIZE')
20004 FORMAT
+ (/1X, A9, 'NEWTON COMPLETED THE TIME STEP.'
+ //10X, I10, ' TIME STEP NUMBER'
+ /10X, A10, ' LOG10 MAX NORM STEADY STATE RESIDUAL')
20005 FORMAT
+ (/1X, A9, 'NEWTON DID NOT COMPLETE THE TIME STEP.')
20008 FORMAT
+ (/10X, 'THE STEADY STATE RESIDUAL IS SO SMALL THAT A STEADY'
+ /10X, 'STATE SOLUTION MAY HAVE BEEN OBTAINED.')
20009 FORMAT
+ (/10X, 'THE DESIRED NUMBER OF TIME STEPS HAS BEEN PERFORMED.')
20010 FORMAT
+ (/10X, 'THE TIME STEP IS NOT ALLOWED TO DECREASE.')
20011 FORMAT
+ (/10X, 'THE TIME STEP IS NOT ALLOWED TO INCREASE.')
20012 FORMAT
+ (/10X, 'THE TIME STEP WAS NOT COMPLETED DESPITE REPEATED'
+ /10X, 'CHANGES TO THE TIME STEP SIZE.')
20013 FORMAT
+ (/10X, 'THE SOLUTION:')
C///////////////////////////////////////////////////////////////////////
C
C ERROR MESSAGES.
C
C///////////////////////////////////////////////////////////////////////
9001 IF(TEXT.GT.0) WRITE (TEXT, 99001) ID, DESIRE
GO TO 99999
9002 IF(TEXT.GT.0) WRITE (TEXT, 99002) ID, TDEC, TINC
GO TO 99999
9003 IF(TEXT.GT.0) WRITE (TEXT, 99003) ID, RETIRE
GO TO 99999
9004 IF(TEXT.GT.0) WRITE (TEXT, 99004) ID, STEP
GO TO 99999
9005 IF(TEXT.GT.0) WRITE (TEXT, 99005) ID, AGE, TRIALS, QMOVE, QSTAT
GO TO 99999
9006 IF(TEXT.GT.0) WRITE (TEXT, 99006) ID
GO TO 99999
99001 FORMAT
+ (/1X, A9, 'ERROR. THE DESIRED NUMBER OF TIME STEPS MUST BE'
+ /10X, 'POSITIVE.'
+ //1X, I10, ' DESIRED NUMBER')
99002 FORMAT
+ (/1X, A9, 'ERROR. THE FACTORS BY WHICH THE TIME STEP MAY'
+ /10X, 'CHANGE MUST BE GREATER THAN ONE.'
+ //1X, E10.2, ' INCREASE FACTOR'
+ /1X, E10.2, ' DECREASE FACTOR')
99003 FORMAT
+ (/1X, A9, 'ERROR. THE RETIREMENT AGE OF TIME STEP SIZES MUST'
+ /10X, 'BE POSITIVE.'
+ //1X, I10, ' AGE')
99004 FORMAT
+ (/1X, A9, 'ERROR. THE STEP NUMBER MUST BE ZERO OR POSITIVE.'
+ //1X, I10, ' NUMBER')
99005 FORMAT
+ (/1X, A9, 'ERROR. THE BEHAVIOR OF THE NEWTON ALGORITHM IS NOT'
+ /10X, 'CONSISTENT WITH THE SIZE OF THE TIME STEP.'
+ //1X, I10, ' AGE OF THE PRESENT STEP'
+ /1X, I10, ' TRIALS TO FIND A GOOD STEP'
+ //1X, I10, ' QMOVE'
+ /1X, I10, ' QSTAT')
99006 FORMAT
+ (/1X, A9, 'ERROR. NEWTON ERRS.')
C/// EXIT.
99999 CONTINUE
END
c
SUBROUTINE TIMWOR (VALUE, WORD)
IMPLICIT COMPLEX (A - Z)
C CHARACTER DIGIT*1, WORD*9
C THE CRAY/CTSS FORTLIB I/O ROUTINES REQUIRE THAT CHARACTER STRINGS
C USED AS FORMATS OR INTERNAL FILES HAVE WORD BOUNDARIES. THUS, THE
C STRINGS' LENGTHS MUST BE DIVISIBLE BY 8.
CHARACTER DIGIT*1, WORD*16
C*****precision > double
DOUBLE PRECISION SECOND, SIXTY, VALUE, ZERO
C*****END precision > double
INTEGER HOUR, J, MINUTE
LOGICAL FOUND
C*****precision > single
C REAL SECOND, SIXTY, VALUE, ZERO
C*****END precision > single
C PRECISION-INDEPENDENT CONSTANT
ZERO = 0.0
SIXTY = 60.0
WORD = ' '
IF(VALUE.GT.ZERO) THEN
SECOND = MOD (VALUE, SIXTY)
MINUTE = MOD (INT (VALUE) / 60, 60)
HOUR = INT (VALUE) / 3600
IF(HOUR.LE.9) THEN
WRITE (WORD, 10001, ERR = 0200) HOUR, MINUTE, SECOND
ELSE
C THE LINE BELOW WAS REPLACED WITH A FUNCTIONALLY EQUAIVALENT LINE
C BECAUSE CRAY/CTSS/FORTLIB LACKS THE INTRINSIC FUNCTION IDNINT.
C WRITE (WORD, 10002, ERR = 0200) HOUR, MINUTE, INT(SECOND)
WRITE (WORD, 10002, ERR = 0200)
+ HOUR, MINUTE, NINT (REAL (SECOND))
ENDIF
FOUND = .FALSE.
DO 0100 J = 1, 9
DIGIT = WORD(J : J)
FOUND = FOUND .OR. (DIGIT.NE.' ' .AND. DIGIT.NE.'0'
+ .AND. DIGIT.NE.':')
IF(FOUND) THEN
IF(DIGIT.EQ.' ') WORD(J : J) = '0'
ELSE
WORD(J : J) = ' '
ENDIF
0100 CONTINUE
ENDIF
0200 CONTINUE
10001 FORMAT (I1, ':', I2, ':', F4.1)
10002 FORMAT (I3, ':', I2, ':', I2)
END
c
SUBROUTINE TWOPNT
+ (ERROR, TEXT, LEVEL,
+ ISIZE, IWORK, RSIZE, RWORK,
+ ABOVE, ACTIVE, ADAPT, BELOW, BINARY, BUFFER, COMPS, CONDIT,
+ FUNCT, JACOB, MARK, MESH, PASS, PASSES, PLIMIT, PMAX, POINTS,
+ REENT, REPORT, SAVE, SCALRS, SHOW, SOLVE, SSABS, SSAGE, SSREL,
+ STEPS0, STEPS1, STEPS2, STORE, SUCCES, TDABS, TDAGE, TDEC,
+ TDREL, TIME, TINC, TMAX, TMIN, TOLER0, TOLER1, TOLER2, TSTEP0,
+ TSTEP1, UPDATE, X)
IMPLICIT COMPLEX (A - P, R - Z), INTEGER (Q)
CHARACTER
+ BANNER*60, COLUMN*16, HEADER*60, ID*9, LINE*60, REMARK*60,
+ WORD*16
EXTERNAL
+ COPY, CPUTIM, TPELPS, EXTENT, GRAB, LOGSTR, NEWTON,
+ REFINE, TIMSTP, TIMWOR
INTEGER
+ BINARY, COMPS, DESIRE, EVENT, GMAX, GRID, ISIZE, IWORK, J,
+ JACOBS, K, LABEL, LAST, LENGTH, LEVEL, PASS, PASSES, PLIMIT,
+ PMAX, POINTS, REPORT, RETURN, ROUTE, RSIZE, SCALRS, SIZE,
+ SSAGE, STEP, STEPS, STEPS0, STEPS1, STEPS2, STMAX, TDAGE, TEXT,
+ XREPOR
LOGICAL
+ ACCEPT, ACTIVE, ADAPT, DECIDE, ERROR, EXIST, FIRST, FOUND,
+ FUNCT, JACOB, MARK, REENT, SATISF, SAVE, SHOW, SOLVE, STORE,
+ TIME, UPDATE, XREENT, SUCCES
C*****precision > double
DOUBLE PRECISION
C*****END precision > double
C*****precision > single
C REAL
C*****END precision > single
+ ABOVE, BELOW, BUFFER, CENT, CONDIT, DETAIL, MAXCON, MESH,
+ RWORK, SSABS, SSREL, TDABS, TDEC, TDREL, TEMP, TIMER, TINC,
+ TMAX, TMIN, TOLER, TOLER0, TOLER1, TOLER2, TOTAL, TSTEP0,
+ TSTEP1, X, YNORM, ZERO
PARAMETER (ID = 'TWOPNT: ')
PARAMETER (GMAX = 100)
PARAMETER (TOLER = 1.0E-10)
C REPORT CODES
PARAMETER (QNULL = 0)
C REPORT CODES FROM NEWTON AND TIMSTP
PARAMETER (QDONE = 1, QDVRG = 2, QBNDS = 3, QSTDY = 4, QSMAX = 5)
C VALUES OF REPORT UPON EXIT WITHOUT SUCCESS AND WITHOUT ERROR.
PARAMETER (QNSLVE = 1, QNSPAC = 2)
C LOCATION OF DATA IN ARRAYS DETAIL, EVENT, TIMER, AND TOTAL. THE
C LOCATIONS ARE CHOSEN TO SIMPLIFY WRITE STATEMENTS. DETAIL USES
C ONLY 1 : 8, EVENT USES ONLY 5 : 8, TIMER USES 1 : 9, AND TOTAL
C USES ONLY 2 : 9. IN ADDITION, 2, 3, 4, 10, AND 11 ARE USED AS
C MNEMONIC VALUES FOR QTASK.
PARAMETER
+ (QGRID = 1,
+ QTIMST = 2,
+ QNEWTO = 3,
+ QREFIN = 4,
+ QFUNCT = 5,
+ QJACOB = 6,
+ QSOLVE = 7,
+ QOTHER = 8,
+ QTOTAL = 9,
+ QENTRY = 10,
+ QEXIT = 11)
DIMENSION
+ ABOVE(SCALRS + COMPS * PMAX), ACTIVE(COMPS), BANNER(4),
+ BELOW(SCALRS + COMPS * PMAX), BUFFER(SCALRS + COMPS * PMAX),
+ COLUMN(6), DETAIL(GMAX, QTOTAL), EVENT(GMAX, QTOTAL),
+ HEADER(3, 3), IWORK(ISIZE), LEVEL(2), MARK(PMAX), MESH(PMAX),
+ RWORK(RSIZE), SIZE(GMAX), TIMER(QTOTAL), TOTAL(QTOTAL),
+ X(SCALRS + COMPS * PMAX)
C/// SAVE LOCAL VALUES DURING RETURNS FOR REVERSE COMMUNCIATION.
SAVE
C///////////////////////////////////////////////////////////////////////
C
C RETURN FROM REVERSE COMMUNICATION.
C
C///////////////////////////////////////////////////////////////////////
C/// TURN OFF ALL REVERSE COMMUNICATION FLAGS.
FUNCT = .FALSE.
JACOB = .FALSE.
SAVE = .FALSE.
SHOW = .FALSE.
SOLVE = .FALSE.
STORE = .FALSE.
TIME = .FALSE.
UPDATE = .FALSE.
C/// TURN OFF THE ERROR FLAG.
ERROR = .FALSE.
C/// CONTINUE WHERE THE PROGRAM PAUSED.
IF(REENT) THEN
REENT = .FALSE.
GO TO ROUTE
ENDIF
C/// TURN OFF THE COMPLETION STATUS FLAGS.
REPORT = QNULL
SUCCES = .FALSE.
C///////////////////////////////////////////////////////////////////////
C
C ENTRY BLOCK. INITIALIZE A NEW PROBLEM.
C
C///////////////////////////////////////////////////////////////////////
C/// CHECK THE ARGUMENTS.
ERROR = LEVEL(1).LT.LEVEL(2)
IF(ERROR) GO TO 9001
ERROR = .NOT. (0.LT.PASSES)
IF(ERROR) GO TO 9002
ERROR = .NOT. (0.LE.SCALRS .AND. 0.LE.COMPS
+ .AND. 0.LE.POINTS .AND. POINTS.LE.PMAX
+ .AND. ((0.LT.COMPS) .EQV. (0.LT.POINTS)))
IF(ERROR) GO TO 9003
ERROR = .NOT.
+ (1.LE.COMPS .AND. 1.LE.PMAX .AND.
+ 1.LE.SCALRS + COMPS * PMAX)
IF(ERROR) GO TO 9004
C/// ONE-TIME INITIALIZATION OF CONSTANTS AND VARIABLES.
C DECISION BLOCK PARAMETERS.
QTASK = QENTRY
C PRECISION-INDEPENDENT CONSTANTS.
ZERO = 0
CENT = 100
C STATISTICS POINTER AND ARRAYS.
GRID = 1
DO 1100 K = 1, QTOTAL
TOTAL(K) = ZERO
DO 1100 J = 1, GMAX
DETAIL(J, K) = ZERO
EVENT(J, K) = 0
1100 CONTINUE
C MARKERS FOR NEW POINTS.
DO 1200 K = 1, POINTS
MARK(K) = .FALSE.
1200 CONTINUE
C TIME STEP NUMBER AND SIZE.
STEP = 0
TSTEP1 = TSTEP0
C/// PARTITION THE INTEGER WORK SPACE.
C SUBROUTINE GRAB (ERROR, TEXT, LAST, FIRST, NUMBER)
LAST = 0
C VARY(PMAX)
CALL GRAB (ERROR, TEXT, LAST, QVARY, PMAX)
IF(ERROR) GO TO 9005
C VARY1(PMAX)
CALL GRAB (ERROR, TEXT, LAST, QVARY1, PMAX)
IF(ERROR) GO TO 9005
C VARY2(PMAX)
CALL GRAB (ERROR, TEXT, LAST, QVARY2, PMAX)
IF(ERROR) GO TO 9005
ERROR = .NOT. (LAST.LE.ISIZE)
IF(ERROR) GO TO 9006
C/// PARTITION THE REAL WORK SPACE.
C SUBROUTINE GRAB (ERROR, TEXT, LAST, FIRST, NUMBER)
LAST = 0
C RATIO1(PMAX)
CALL GRAB (ERROR, TEXT, LAST, QRAT1, PMAX)
IF(ERROR) GO TO 9005
C RATIO2(PMAX)
CALL GRAB (ERROR, TEXT, LAST, QRAT2, PMAX)
IF(ERROR) GO TO 9005
C S(SCALRS + COMPS * PMAX)
CALL GRAB (ERROR, TEXT, LAST, QS, SCALRS + COMPS * PMAX)
IF(ERROR) GO TO 9005
C STRIAL(SCALRS + COMPS * PMAX)
CALL GRAB (ERROR, TEXT, LAST, QSTRL, SCALRS + COMPS * PMAX)
IF(ERROR) GO TO 9005
C XSAVE(SCALRS + COMPS * PMAX)
CALL GRAB (ERROR, TEXT, LAST, QXSAV, SCALRS + COMPS * PMAX)
IF(ERROR) GO TO 9005
C XTEMP(SCALRS + COMPS * PMAX)
CALL GRAB (ERROR, TEXT, LAST, QXTMP, SCALRS + COMPS * PMAX)
IF(ERROR) GO TO 9005
C XTRIAL(SCALRS + COMPS * PMAX)
CALL GRAB (ERROR, TEXT, LAST, QXTRL, SCALRS + COMPS * PMAX)
IF(ERROR) GO TO 9005
C Y(SCALRS + COMPS * PMAX)
CALL GRAB (ERROR, TEXT, LAST, QY, SCALRS + COMPS * PMAX)
IF(ERROR) GO TO 9005
C YTRIAL(SCALRS + COMPS * PMAX)
CALL GRAB (ERROR, TEXT, LAST, QYTRL, SCALRS + COMPS * PMAX)
IF(ERROR) GO TO 9005
ERROR = .NOT. (LAST.LE.RSIZE)
IF(ERROR) GO TO 9007
C/// INITIALIZE THE TOTAL TIME STATISTIC.
CALL CPUTIM (TIMER(QTOTAL))
C/// INITIALIZE THE STATISTICS FOR THE FIRST MESH.
IF(GRID.LE.GMAX) SIZE(GRID) = POINTS
IF(GRID.LE.GMAX) CALL CPUTIM (TIMER(QGRID))
C/// SAVE THE INITIAL SOLUTION.
ASSIGN 1300 TO RETURN
GO TO 9911
1300 CONTINUE
C/// PRINT THE BANNER AT ALL PRINT LEVELS.
BANNER(1) = 'TWO-POINT BOUNDARY VALUE PROBLEM SOLVER,'
BANNER(2) =
+ 'VERSION 2.22B OF MARCH 1991 BY DR. JOSEPH F. GRCAR.'
C*****precision > double
BANNER(3) = 'DOUBLE PRECISION'
C*****END precision > double
C*****precision > single
C BANNER(3) = 'SINGLE PRECISION'
C*****END precision > single
IF(ADAPT) THEN
BANNER(4) = 'PERMIT MESH REFINEMENT'
ELSE
BANNER(4) = 'DO NOT REFINE THE MESH'
ENDIF
IF(LEVEL(1) .GE. 1) THEN
IF(0.LT.TEXT) WRITE (TEXT, 10001) ID, BANNER
ENDIF
C/// PRINT THE INITIAL GUESS AT PRINT LEVELS 11, 21, AND 22.
IF(LEVEL(2) .GE. 1) THEN
IF(0.LT.TEXT) WRITE (TEXT, 10002) ID
ASSIGN 1400 TO RETURN
GO TO 9921
ENDIF
1400 CONTINUE
C/// PRINT LEVEL 10 OR 11.
C 123456789_123456789_123456789_123456
C 123456789 123456789 123456789
HEADER(1, 1) = ' LOG10 LOG10 MAX '
HEADER(2, 1) = ' MAX NORM CONDITION '
HEADER(3, 1) = ' ACTIVITY RESIDUAL NUMBER '
C 123456789_123456789_123456
C 123456789_1234 123456789
HEADER(1, 2) = 'POINTS '
HEADER(2, 2) = ' STEPS '
HEADER(3, 2) = ' JACOBIANS REMARK '
IF(LEVEL(1).NE.1) GO TO 1600
IF(0.LT.TEXT) WRITE (TEXT, 10003)
+ ID, ((HEADER(J, K), K = 1, 2), J = 1, 3)
ASSIGN 1500 TO LABEL
GO TO 7100
1500 CONTINUE
CALL EXTENT (LENGTH, REMARK)
IF(0.LT.TEXT) WRITE (TEXT, 10017)
+ COLUMN, REMARK (1 : LENGTH)
1600 CONTINUE
C///////////////////////////////////////////////////////////////////////
C
C DECISION BLOCK. THE PREVIOUS TASK DETERMINES THE NEXT.
C
C///////////////////////////////////////////////////////////////////////
2100 CONTINUE
C/// ENTRY WAS THE PREVIOUS TASK.
IF(QTASK.EQ.QENTRY) THEN
PASS = 1
IF(0.LT.ABS (STEPS0)) THEN
QTASK = QTIMST
DESIRE = ABS (STEPS0)
ELSE
QTASK = QNEWTO
ENDIF
C/// NEWTON WAS THE PREVIOUS TASK.
ELSEIF(QTASK.EQ.QNEWTO) THEN
IF(FOUND) THEN
IF(PASS.LT.PASSES) THEN
PASS = PASS + 1
QTASK = QNEWTO
ELSE
PASS = 1
IF(ADAPT) THEN
QTASK = QREFIN
ELSE
QTASK = QEXIT
SUCCES = .TRUE.
ENDIF
ENDIF
ELSE
IF(STEPS1.GT.0) THEN
QTASK = QTIMST
DESIRE = STEPS1
ELSE
QTASK = QEXIT
REPORT = QNSLVE
SUCCES = .FALSE.
ENDIF
ENDIF
C/// REFINE WAS THE PREVIOUS TASK.
ELSEIF(QTASK.EQ.QREFIN) THEN
IF(FOUND) THEN
PASS = 1
STEP = 0
C RESET THE TIME STEP.
TSTEP1 = TSTEP0
QTASK = QNEWTO
ELSE
QTASK = QEXIT
IF(.NOT. SATISF) REPORT = QNSPAC
SUCCES = SATISF
ENDIF
C/// TIMSTP WAS THE PREVIOUS TASK.
ELSEIF(QTASK.EQ.QTIMST) THEN
IF(FOUND) THEN
IF(STEPS0.LT.0) THEN
QTASK = QEXIT
SUCCES = .TRUE.
ELSE
QTASK = QNEWTO
ENDIF
ELSE
QTASK = QEXIT
REPORT = QNSLVE
SUCCES = .FALSE.
ENDIF
ENDIF
C/// BRANCH TO THE NEXT TASK.
IF(QTASK.EQ.QEXIT) GO TO 3100
IF(QTASK.EQ.QNEWTO) GO TO 4100
IF(QTASK.EQ.QREFIN) GO TO 5100
IF(QTASK.EQ.QTIMST) GO TO 6100
ERROR = .TRUE.
GO TO 9008
C///////////////////////////////////////////////////////////////////////
C
C EXIT BLOCK.
C
C///////////////////////////////////////////////////////////////////////
3100 CONTINUE
C/// COMPLETE STATISTICS FOR THE LAST MESH.
CALL TPELPS (TIMER(QGRID))
IF(GRID.LE.GMAX) DETAIL(GRID, QGRID) = TIMER(QGRID)
C/// PRINT LEVEL 11 OR 21.
IF(LEVEL(2).EQ.1) THEN
IF(0.LT.TEXT) WRITE (TEXT, 10005) ID
ASSIGN 3200 TO RETURN
GO TO 9921
ENDIF
3200 CONTINUE
C/// COMPLETE THE TOTAL TIME STATISTIC.
CALL TPELPS (TIMER(QTOTAL))
TOTAL(QTOTAL) = TIMER(QTOTAL)
C/// REPORT THE TOTAL TIME.
IF(LEVEL(1) .GE. 1) THEN
IF(TOTAL(QTOTAL) .GE. 36000.0) THEN
LINE = 'HOURS : MINUTES : SECONDS'
ELSEIF(TOTAL(QTOTAL) .GE. 3600.0) THEN
LINE = 'HOURS : MINUTES : SECONDS . FRACTION'
ELSEIF(TOTAL(QTOTAL) .GE. 60.0) THEN
LINE = 'MINUTES : SECONDS . FRACTION'
ELSE
LINE = 'SECONDS . FRACTION'
ENDIF
IF(TOTAL(QTOTAL).GT.ZERO) THEN
CALL TIMWOR (TOTAL(QTOTAL), WORD)
IF(0.LT.TEXT) WRITE (TEXT, 10006) ID, WORD, LINE
ENDIF
ENDIF
C/// REPORT THE PERCENT OF TOTAL CPU TIME.
C 123456789_123456
C 123456 123456
HEADER(1, 1) = ' '
HEADER(2, 1) = ' GRID GRID '
HEADER(3, 1) = 'POINTS TOTAL '
C 123456789_123456789_12
C 123456 123456 123456
HEADER(1, 2) = 'ACTIVITY '
HEADER(2, 2) = '-------------------- '
HEADER(3, 2) = 'TIMSTP NEWTON REFINE '
C 123456789_123456789_1234567
C 123456 123456 123456 123456
HEADER(1, 3) = 'EVENT '
HEADER(2, 3) = '---------------------------'
HEADER(3, 3) = 'FUNCTN JACOBN SOLVE OTHER'
IF(LEVEL(1) .GE. 1) THEN
IF(TOTAL(QTOTAL).GT.ZERO) THEN
IF(0.LT.TEXT) WRITE (TEXT, 10007)
+ ID, ((HEADER(J, K), K = 1, 3), J = 1, 3),
+ (SIZE(J), (CENT * DETAIL(J, K) / TOTAL(QTOTAL),
+ K = 1, 8), J = 1, GRID)
IF(GRID.GT.1) THEN
IF(0.LT.TEXT) WRITE (TEXT, 10008)
+ (CENT * TOTAL(K) / TOTAL(QTOTAL), K = 2, 8)
ENDIF
IF(GRID.GT.GMAX) THEN
IF(0.LT.TEXT) WRITE (TEXT, 10009)
ENDIF
ENDIF
ENDIF
C/// REPORT THE AVERAGE CPU TIME.
C 123456789
C 123456
HEADER(1, 1) = ' '
HEADER(2, 1) = ' GRID '
HEADER(3, 1) = 'POINTS '
C 123456789_123456789_12345678
C 1234567 1234567 1234567
HEADER(1, 2) = 'AVERAGE SECONDS '
HEADER(2, 2) = '------------------------- '
HEADER(3, 2) = ' FUNCTN JACOBN SOLVE '
C 123456789_123456789_12345
C 1234567 1234567 1234567
HEADER(1, 3) = 'NUMBER OF EVENTS '
HEADER(2, 3) = '-------------------------'
HEADER(3, 3) = ' FUNCTN JACOBN SOLVE'
IF(LEVEL(1) .GE. 1) THEN
IF(TOTAL(QTOTAL).GT.ZERO) THEN
IF(0.LT.TEXT) WRITE (TEXT, 10010)
+ ID, ((HEADER(J, K), K = 1, 3), J = 1, 3),
+ (SIZE(J), (DETAIL(J, K) / EVENT(J, K), K = 5, 7),
+ (EVENT(J, K), K = 5, 7), J = 1, GRID)
IF(GRID.GT.GMAX) THEN
IF(0.LT.TEXT) WRITE (TEXT, 10011)
ENDIF
ENDIF
ENDIF
C/// REPORT THE COMPLETION STATUS.
IF(LEVEL(1) .GE. 1) THEN
IF(SUCCES) THEN
IF(0.LT.TEXT) WRITE (TEXT, 10012) ID
ELSEIF(REPORT.EQ.QNSLVE) THEN
IF(GRID.EQ.1) THEN
IF(0.LT.TEXT) WRITE (TEXT, 10013) ID
ELSE
IF(0.LT.TEXT) WRITE (TEXT, 10014) ID
ENDIF
ELSEIF(REPORT.EQ.QNSPAC) THEN
IF(0.LT.TEXT) WRITE (TEXT, 10015) ID
ELSE
ERROR = .TRUE.
GO TO 9009
ENDIF
ENDIF
C/// BOTTOM OF THE EXIT BLOCK.
GO TO 99999
C///////////////////////////////////////////////////////////////////////
C
C NEWTON BLOCK.
C
C///////////////////////////////////////////////////////////////////////
4100 CONTINUE
C/// INITIALIZE STATISTICS ON ENTRY TO THE NEWTON BLOCK.
CALL CPUTIM (TIMER(QNEWTO))
FIRST = .TRUE.
JACOBS = 0
MAXCON = ZERO
C/// PRINT LEVEL 20, 21, OR 22 ON ENTRY TO THE NEWTON BLOCK.
IF(LEVEL(1) .GE. 2) THEN
IF(0.LT.TEXT) WRITE (TEXT, 10016) ID
ENDIF
C/// PREPARE TO CALL NEWTON.
C SAVE THE SOLUTION SHOULD NEWTON FAIL.
CALL COPY (SCALRS + COMPS * POINTS, X, RWORK(QXSAV))
EXIST = .FALSE.
C/// CALL NEWTON.
STMAX = 0
XREENT = .FALSE.
4200 CONTINUE
C SUBROUTINE NEWTON
C + (ERROR, TEXT, ABOVE, ACCEPT, BELOW, BUFFER, CONDIT, DECIDE,
C + EXIST, FUNCT, JACOB, LEVEL1, LEVEL2, N, REENT, REPORT, RETIRE,
C + S, SHOW, SOLVE, STEPS, STMAX, STRIAL, SUCCES, X, XTRIAL, Y,
C + YNORM, YTRIAL)
CALL NEWTON
+ (ERROR, TEXT, ABOVE, ACCEPT, BELOW, BUFFER, CONDIT, DECIDE,
+ EXIST, FUNCT, JACOB, LEVEL(1) - 1, LEVEL(2) - 1,
+ SCALRS + COMPS * POINTS, XREENT, XREPOR, SSAGE, RWORK(QS),
+ SHOW, SOLVE, STEPS, STMAX, RWORK(QSTRL), FOUND, X,
+ RWORK(QXTRL), RWORK(QY), YNORM, RWORK(QYTRL))
IF(ERROR) GO TO 9010
C/// SERVICE REQUESTS FROM NEWTON: DECIDE WHETHER THE SOLUTION IS
C/// ACCEPTABLE.
IF(XREENT) THEN
IF(DECIDE) THEN
ACCEPT = .TRUE.
DO 4300 J = 1, SCALRS + COMPS * POINTS
ACCEPT = ACCEPT .AND.
+ ABS (RWORK(QS - 1 + J)) .LE.
+ MAX (SSABS, SSREL * ABS (X(J)))
4300 CONTINUE
GO TO 4200
C/// PASS OTHER REQUESTS FROM NEWTON TO THE CALLER.
ELSE
ASSIGN 4200 TO RETURN
GO TO 9931
ENDIF
ENDIF
C/// REACT TO THE COMPLETION OF NEWTON.
IF(FOUND) THEN
C SAVE THE LATEST SOLUTION.
ASSIGN 4400 TO RETURN
GO TO 9911
ELSE
C RESTORE THE SOLUTION.
CALL COPY (SCALRS + COMPS * POINTS, RWORK(QXSAV), X)
ENDIF
4400 CONTINUE
C/// COMPLETE STATISTICS FOR THE NEWTON BLOCK.
CALL TPELPS (TIMER(QNEWTO))
TOTAL(QNEWTO) = TOTAL(QNEWTO) + TIMER(QNEWTO)
IF(GRID.LE.GMAX)
+ DETAIL(GRID, QNEWTO) = DETAIL(GRID, QNEWTO) + TIMER(QNEWTO)
C/// PRINT LEVEL 10 OR 11 ON EXIT FROM THE NEWTON BLOCK.
IF(LEVEL(1).NE.1) GO TO 4600
ASSIGN 4500 TO LABEL
GO TO 7100
4500 CONTINUE
CALL EXTENT (LENGTH, REMARK)
IF(0.LT.TEXT) WRITE (TEXT, 10017)
+ COLUMN, REMARK (1 : LENGTH)
4600 CONTINUE
C/// PRINT LEVEL 20, 21, OR 22 ON EXIT FROM THE NEWTON BLOCK.
IF(LEVEL(1) .GE. 2) THEN
IF(FOUND) THEN
IF(0.LT.TEXT) WRITE (TEXT, 10018) ID
ELSE
IF(0.LT.TEXT) WRITE (TEXT, 10019) ID
ENDIF
ENDIF
C/// BOTTOM OF THE NEWTON BLOCK.
GO TO 2100
C///////////////////////////////////////////////////////////////////////
C
C REFINE BLOCK.
C
C///////////////////////////////////////////////////////////////////////
5100 CONTINUE
C/// INITIALIZE STATISTICS ON ENTRY TO THE REFINE BLOCK.
CALL CPUTIM (TIMER(QREFIN))
C/// PRINT LEVEL 20, 21, OR 22 ON ENTRY TO THE REFINE BLOCK.
IF(LEVEL(1) .GE. 2) THEN
IF(0.LT.TEXT) WRITE (TEXT, 10020) ID
ENDIF
C/// CALL REFINE.
XREENT = .FALSE.
5200 CONTINUE
C SUBROUTINE REFINE
C + (ERROR, TEXT, ABOVE, ACTIVE, BELOW, BUFFER, COMPS, LEVEL1,
C + LEVEL2, MARK, MESH, N, NEWMSH, PLIMIT, PMAX, POINTS, RATIO1,
C + RATIO2, REENT, SHOW, SUCCES, TOLER0, TOLER1, TOLER2, UPDATE,
C + VARY, VARY1, VARY2, X)
CALL REFINE
+ (ERROR, TEXT, ABOVE(SCALRS + 1), ACTIVE, BELOW(SCALRS + 1),
+ BUFFER(SCALRS + 1), COMPS, LEVEL(1) - 1, LEVEL(2) - 1, MARK,
+ MESH, FOUND, PLIMIT, PMAX, POINTS, RWORK(QRAT1), RWORK(QRAT2),
+ XREENT, SHOW, SATISF, TOLER0, TOLER1, TOLER2, UPDATE,
+ IWORK(QVARY), IWORK(QVARY1), IWORK(QVARY2), X(SCALRS + 1))
IF(ERROR) GO TO 9011
C/// SERVICE REQUESTS FROM REFINE: PASS REQUESTS TO THE CALLER.
IF(XREENT) THEN
IF(0.LT.SCALRS) CALL COPY (SCALRS, X, BUFFER)
ASSIGN 5200 TO RETURN
GO TO 9931
ENDIF
C/// REACT TO THE COMPLETION OF REFINE.
IF(.NOT. FOUND) GO TO 5400
C COMPLETE STATISTICS FOR THE OLD MESH.
CALL TPELPS (TIMER(QGRID))
IF(GRID.LE.GMAX) DETAIL(GRID, QGRID) = TIMER(QGRID)
C INITIALIZE STATISTICS FOR THE NEW MESH.
GRID = GRID + 1
IF(GRID.LE.GMAX) THEN
CALL CPUTIM (TIMER(QGRID))
SIZE(GRID) = POINTS
ENDIF
C SAVE THE LATEST SOLUTION.
ASSIGN 5300 TO RETURN
GO TO 9911
5300 CONTINUE
5400 CONTINUE
C/// COMPLETE STATISTICS FOR THE REFINE BLOCK.
CALL TPELPS (TIMER(QREFIN))
TOTAL(QREFIN) = TOTAL(QREFIN) + TIMER(QREFIN)
IF(GRID.LE.GMAX)
+ DETAIL(GRID, QREFIN) = DETAIL(GRID, QREFIN) + TIMER(QREFIN)
C/// PRINT LEVEL 10 OR 11 ON EXIT FROM THE REFINE BLOCK.
IF(LEVEL(1).NE.1) GO TO 5600
ASSIGN 5500 TO LABEL
GO TO 7100
5500 CONTINUE
CALL EXTENT (LENGTH, REMARK)
IF(0.LT.TEXT) WRITE (TEXT, 10004)
+ COLUMN, REMARK (1 : LENGTH)
5600 CONTINUE
C/// PRINT LEVEL 20, 21, OR 22 ON EXIT FROM THE REFINE BLOCK.
IF(LEVEL(1) .GE. 2) THEN
IF(FOUND) THEN
IF(0.LT.TEXT) WRITE (TEXT, 10021) ID
ELSE
IF(0.LT.TEXT) WRITE (TEXT, 10022) ID
ENDIF
ENDIF
C/// BOTTOM OF THE REFINE BLOCK.
GO TO 2100
C///////////////////////////////////////////////////////////////////////
C
C TIMSTP BLOCK.
C
C///////////////////////////////////////////////////////////////////////
6100 CONTINUE
C/// INITIALIZE STATISTICS ON ENTRY TO THE TIMSTP BLOCK.
CALL CPUTIM (TIMER(QTIMST))
FIRST = .TRUE.
JACOBS = 0
MAXCON = ZERO
STEPS = STEP
C/// PRINT LEVEL 20, 21, OR 22 ON ENTRY TO THE TIMSTP BLOCK.
IF(LEVEL(1) .GE. 2) THEN
IF(0.LT.TEXT) WRITE (TEXT, 10023) ID
ENDIF
C/// CALL TIMSTP.
STMAX = 100
XREENT = .FALSE.
6200 CONTINUE
C SUBROUTINE TIMSTP
C + (ERROR, TEXT, ABOVE, ACCEPT, BELOW, BUFFER, CONDIT, DECIDE,
C + DESIRE, FUNCT, JACOB, LEVEL1, LEVEL2, N, REENT, REPORT,
C + RETIRE, S, SHOW, SIZE1, SOLVE, STEP, STMAX, STORE, STRIAL,
C + SUCCES, TDAGE, TDEC, TIME, TINC, TMAX, TMIN, TOLER, X, XSAVE,
C + XTRIAL, Y, YNORM, YTRIAL)
CALL TIMSTP
+ (ERROR, TEXT, ABOVE, ACCEPT, BELOW, BUFFER, CONDIT, DECIDE,
+ DESIRE, FUNCT, JACOB, LEVEL(1) - 1, LEVEL(2) - 1,
+ SCALRS + COMPS * POINTS, XREENT, XREPOR, STEPS2, RWORK(QS),
+ SHOW, TSTEP1, SOLVE, STEP, STMAX, STORE, RWORK(QSTRL), FOUND,
+ TDAGE, TDEC, TIME, TINC, TMAX, TMIN, TOLER, X, RWORK(QXTMP),
+ RWORK(QXTRL), RWORK(QY), YNORM, RWORK(QYTRL))
IF(ERROR) GO TO 9012
C/// SERVICE REQUESTS FROM TIMSTP: DECIDE WHETHER THE SOLUTION IS
C/// ACCEPTABLE.
IF(XREENT) THEN
IF(DECIDE) THEN
ACCEPT = .TRUE.
DO 6300 J = 1, SCALRS + COMPS * POINTS
ACCEPT = ACCEPT .AND.
+ ABS (RWORK(QS - 1 + J)) .LE.
+ MAX (TDABS, TDREL * ABS (X(J)))
6300 CONTINUE
GO TO 6200
C/// PASS OTHER REQUESTS FROM TIMSTP TO THE CALLER.
ELSE
ASSIGN 6200 TO RETURN
GO TO 9931
ENDIF
ENDIF
C/// REACT TO THE COMPLETION OF TIMSTP.
IF(FOUND) THEN
C SAVE THE LATEST SOLUTION.
ASSIGN 6400 TO RETURN
GO TO 9911
ENDIF
6400 CONTINUE
C/// COMPLETE STATISTICS FOR THE TIMSTP BLOCK.
CALL TPELPS (TIMER(QTIMST))
TOTAL(QTIMST) = TOTAL(QTIMST) + TIMER(QTIMST)
IF(GRID.LE.GMAX)
+ DETAIL(GRID, QTIMST) = DETAIL(GRID, QTIMST) + TIMER(QTIMST)
STEPS = STEP - STEPS
C/// PRINT LEVEL 10 OR 11 ON EXIT FROM THE TIMSTP BLOCK.
IF(LEVEL(1).NE.1) GO TO 6600
ASSIGN 6500 TO LABEL
GO TO 7100
6500 CONTINUE
CALL EXTENT (LENGTH, REMARK)
IF(0.LT.TEXT) WRITE (TEXT, 10017)
+ COLUMN, REMARK (1 : LENGTH)
6600 CONTINUE
C/// PRINT LEVEL 20, 21, OR 22 ON EXIT FROM THE TIMSTP BLOCK.
IF(LEVEL(1) .GE. 2) THEN
IF(FOUND) THEN
IF(0.LT.TEXT) WRITE (TEXT, 10024) ID
ELSE
IF(0.LT.TEXT) WRITE (TEXT, 10025) ID
ENDIF
ENDIF
C/// BOTTOM OF THE TIMSTP BLOCK.
GO TO 2100
C///////////////////////////////////////////////////////////////////////
C
C BLOCK TO PREPARE LOG LINES.
C
C///////////////////////////////////////////////////////////////////////
7100 CONTINUE
DO 7200 J = 1, 6
COLUMN(J) = ' '
7200 CONTINUE
REMARK = ' '
C COLUMN 1: NAME OF THE ACTIVITY
IF(QTASK.EQ.QENTRY) COLUMN(1) = ' START'
IF(QTASK.EQ.QNEWTO) COLUMN(1) = ' NEWTON'
IF(QTASK.EQ.QREFIN) COLUMN(1) = ' REFINE'
IF(QTASK.EQ.QTIMST) COLUMN(1) = ' TIMSTP'
C COLUMN 2: NORM OF THE STEADY-STATE FUNCTION
IF(.NOT. FOUND) GO TO 7400
ASSIGN 7300 TO RETURN
GO TO 9941
7300 CONTINUE
CALL TPNORM (SCALRS + COMPS * POINTS, TEMP, BUFFER)
CALL LOGSTR ('(F9.2)', 9, COLUMN(2), TEMP)
7400 CONTINUE
C COLUMN 3: LARGEST CONDITION NUMBER
IF(QTASK.EQ.QNEWTO .OR. QTASK.EQ.QTIMST) THEN
IF(MAXCON.NE.ZERO)
+ CALL LOGSTR ('(F9.2)', 9, COLUMN(3), MAXCON)
ENDIF
C COLUMN 4: NUMBER OF POINTS
IF(QTASK.EQ.QENTRY .OR. QTASK.EQ.QREFIN) THEN
WRITE (COLUMN(4), '(I4)') POINTS
ENDIF
C COLUMN 5: NUMBER OF STEPS
IF(QTASK.EQ.QNEWTO .OR. QTASK.EQ.QTIMST) THEN
WRITE (COLUMN(5), '(I4)') STEPS
ENDIF
C COLUMN 6: NUMBER OF JACOBIANS
IF(QTASK.EQ.QNEWTO .OR. QTASK.EQ.QTIMST) THEN
WRITE (COLUMN(6), '(I4)') JACOBS
ENDIF
C REMARK
IF(QTASK.EQ.QNEWTO .OR. QTASK.EQ.QTIMST) THEN
IF(XREPOR.EQ.QDVRG) THEN
REMARK = 'DIVERGING'
ELSEIF(XREPOR.EQ.QNULL) THEN
REMARK = ' '
ELSEIF(XREPOR.EQ.QBNDS) THEN
REMARK = 'AT BOUNDS'
ELSEIF(XREPOR.EQ.QDONE) THEN
REMARK = ' '
ELSEIF(XREPOR.EQ.QSMAX) THEN
REMARK = 'STEP LIMIT'
ELSEIF(XREPOR.EQ.QSTDY) THEN
REMARK = 'STEADY STATE'
ELSE
REMARK = '?'
ENDIF
ELSEIF(QTASK.EQ.QREFIN) THEN
IF(FOUND) REMARK = 'NEW MESH'
ENDIF
GO TO LABEL
C///////////////////////////////////////////////////////////////////////
C
C REQUEST REVERSE COMMUNICATION.
C
C///////////////////////////////////////////////////////////////////////
C/// SAVE THE SOLUTION.
9911 CONTINUE
CALL CPUTIM (TIMER(QOTHER))
CALL COPY (SCALRS + COMPS * POINTS, X, BUFFER)
SAVE = .TRUE.
REENT = .TRUE.
ASSIGN 9912 TO ROUTE
GO TO 99999
9912 CONTINUE
CALL TPELPS (TIMER(QOTHER))
TOTAL(QOTHER) = TOTAL(QOTHER) + TIMER(QOTHER)
IF(GRID.LE.GMAX) THEN
DETAIL(GRID, QOTHER) = DETAIL(GRID, QOTHER) + TIMER(QOTHER)
EVENT(GRID, QOTHER) = EVENT(GRID, QOTHER) + 1
ENDIF
GO TO RETURN
C/// PRINT THE LATEST SOLUTION.
9921 CONTINUE
CALL CPUTIM (TIMER(QOTHER))
CALL COPY (SCALRS + COMPS * POINTS, X, BUFFER)
SHOW = .TRUE.
REENT = .TRUE.
ASSIGN 9922 TO ROUTE
GO TO 99999
9922 CONTINUE
CALL TPELPS (TIMER(QOTHER))
TOTAL(QOTHER) = TOTAL(QOTHER) + TIMER(QOTHER)
IF(GRID.LE.GMAX) THEN
DETAIL(GRID, QOTHER) = DETAIL(GRID, QOTHER) + TIMER(QOTHER)
EVENT(GRID, QOTHER) = EVENT(GRID, QOTHER) + 1
ENDIF
GO TO RETURN
C/// PASS REQUESTS FROM NEWTON, REFINE, OR TIMSTP TO THE CALLER.
9931 CONTINUE
C IDENTIFY THE REQUEST. THIS MUST BE SAVED TO GATHER STATISTICS
C AT REENTRY. THE REVERSE COMMUNICATION FLAGS WILL NOT BE SAVED
C BECAUSE THEY ARE CLEARED AT EVERY ENTRY.
IF(FUNCT) THEN
QTYPE = QFUNCT
ELSEIF(JACOB) THEN
QTYPE = QJACOB
ELSEIF(SOLVE) THEN
QTYPE = QSOLVE
ELSE
QTYPE = QOTHER
ENDIF
C COUNT THE NUMBER OF JACOBIANS FOR INCLUSION IN THE LOG LINE.
IF(JACOB) JACOBS = JACOBS + 1
CALL CPUTIM (TIMER(QTYPE))
REENT = .TRUE.
ASSIGN 9932 TO ROUTE
GO TO 99999
9932 CONTINUE
C SAVE THE MAXIMUM CONDITION NUMBER.
IF(QTYPE.EQ.QJACOB) MAXCON = MAX (MAXCON, CONDIT)
CALL TPELPS (TIMER(QTYPE))
TOTAL(QTYPE) = TOTAL(QTYPE) + TIMER(QTYPE)
IF(GRID.LE.GMAX) THEN
DETAIL(GRID, QTYPE) = DETAIL(GRID, QTYPE) + TIMER(QTYPE)
EVENT(GRID, QTYPE) = EVENT(GRID, QTYPE) + 1
ENDIF
GO TO RETURN
C/// EVALUATE THE STEADY-STATE FUNCTION.
9941 CONTINUE
CALL CPUTIM (TIMER(QFUNCT))
CALL COPY (SCALRS + COMPS * POINTS, X, BUFFER)
FUNCT = .TRUE.
TIME = .FALSE.
REENT = .TRUE.
ASSIGN 9942 TO ROUTE
GO TO 99999
9942 CONTINUE
CALL TPELPS (TIMER(QFUNCT))
TOTAL(QFUNCT) = TOTAL(QFUNCT) + TIMER(QFUNCT)
IF(GRID.LE.GMAX) THEN
DETAIL(GRID, QFUNCT) = DETAIL(GRID, QFUNCT) + TIMER(QFUNCT)
EVENT(GRID, QFUNCT) = EVENT(GRID, QFUNCT) + 1
ENDIF
GO TO RETURN
C///////////////////////////////////////////////////////////////////////
C
C INFORMATIVE MESSAGES.
C
C///////////////////////////////////////////////////////////////////////
10001 FORMAT
+ (/1X, A9, A
+ /10X, A
+ //20X, A
+ /20X, A
+ /20X, A)
10002 FORMAT
+ (/1X, A9, 'INITIAL GUESS:')
10003 FORMAT
+ (/1X, A9, 'LOG OF THE BOUNDARY VALUE PROBLEM SOLVER.'
+ //10X, A36, A26
+ /10X, A36, A26
+ /10X, A36, A26)
10004 FORMAT
+ (/10X, A9,
+ 3X, A9,
+ 3X, A9,
+ 3X, A4,
+ 1X, A4,
+ 1X, A4,
+ 3X, A9,
+ 3X, A9)
10005 FORMAT
+ (/1X, A9, 'FINAL SOLUTION:')
10006 FORMAT
+ (/1X, A9, 'TOTAL CPU TIME.'
+ //10X, A9, 2X, A)
10007 FORMAT
+ (/1X, A9, 'PERCENT OF TOTAL CPU TIME.'
+ /3(/10X, A16, A22, A27)
+ //(10X, I6, 2X, F6.1, 1X, 3(1X, F6.1), 1X, 4(1X, F6.1)))
10008 FORMAT
+ (/10X, ' TOTAL', 9X, 3(1X, F6.1), 1X, 4(1X, F6.1))
10009 FORMAT
+ (/10X, 'SPACE LIMITATIONS PREVENT RETAINING TIMING DATA FOR'
+ /10X, 'EVERY GRID. NEVERTHELESS, THE TOTALS REFER TO ALL'
+ /10X, 'GRIDS; THEY ARE NOT COLUMN SUMS.')
10010 FORMAT
+ (/1X, A9, 'AVERAGE CPU TIME.'
+ /3(/10X, A9, A28, A25)
+ //(10X, I6, 3X, F7.3, 2X, F7.2, 2X, F7.3, 1X, 3(2X, I7)))
10011 FORMAT
+ (/10X, 'SPACE LIMITATIONS PREVENT RETAINING TIMING DATA FOR'
+ /10X, 'EVERY GRID.')
10012 FORMAT
+ (/1X, A9, 'SUCCESS. BOUNDARY VALUE PROBLEM SOLVED.')
10013 FORMAT
+ (/1X, A9, 'FAILURE. NO SOLUTION WAS FOUND, NOT EVEN ON THE'
+ /10X, 'FIRST GRID.')
10014 FORMAT
+ (/1X, A9, 'FAILURE. A SOLUTION WAS FOUND FOR SOME GRIDS, BUT'
+ /10X, 'NOT THE LAST.')
10015 FORMAT
+ (/1X, A9, 'FAILURE. THE SOLUTION IS NOT RESOLVED TO THE EXTENT'
+ /10X, 'DESIRED BECAUSE THERE IS NO SPACE FOR MORE POINTS.')
10016 FORMAT
+ (/1X, A9, 'CALLING NEWTON TO SOLVE THE STEADY-STATE PROBLEM.')
10017 FORMAT
+ (10X, A9,
+ 3X, A9,
+ 3X, A9,
+ 3X, A4,
+ 1X, A4,
+ 1X, A4,
+ 3X, A)
10018 FORMAT
+ (/1X, A9, 'NEWTON SOLVED THE STEADY-STATE PROBLEM.')
10019 FORMAT
+ (/1X, A9, 'NEWTON DID NOT SOLVE THE STEADY-STATE PROBLEM.')
10020 FORMAT
+ (/1X, A9, 'CALLING REFINE TO ADAPT THE MESH.')
10021 FORMAT
+ (/1X, A9, 'REFINE PRODUCED A NEW MESH.')
10022 FORMAT
+ (/1X, A9, 'REFINE DID NOT PRODUCE A NEW MESH.')
10023 FORMAT
+ (/1X, A9, 'CALLING TIMSTP TO PERFORM A TIME INTEGRATION.')
10024 FORMAT
+ (/1X, A9, 'TIMSTP COMPLETED THE TIME INTEGRATION.')
10025 FORMAT
+ (/1X, A9, 'TIMSTP DID NOT COMPLETE THE TIME INTEGRATION.')
C///////////////////////////////////////////////////////////////////////
C
C ERROR MESSAGES.
C
C///////////////////////////////////////////////////////////////////////
C GO TO 99999
9001 IF(0.LT.TEXT) WRITE (TEXT, 99001) ID
GO TO 99999
9002 IF(0.LT.TEXT) WRITE (TEXT, 99002) ID, PASSES
GO TO 99999
9003 IF(0.LT.TEXT) WRITE (TEXT, 99003) ID,
+ SCALRS, COMPS, POINTS, PMAX
GO TO 99999
9004 IF(0.LT.TEXT) WRITE (TEXT, 99004) ID,
+ COMPS, PMAX, SCALRS + COMPS * PMAX
GO TO 99999
9005 IF(0.LT.TEXT) WRITE (TEXT, 99005) ID
GO TO 99999
9006 IF(0.LT.TEXT) WRITE (TEXT, 99006) ID, LAST, ISIZE
GO TO 99999
9007 IF(0.LT.TEXT) WRITE (TEXT, 99007) ID, LAST, RSIZE
GO TO 99999
9008 IF(0.LT.TEXT) WRITE (TEXT, 99008) ID
GO TO 99999
9009 IF(0.LT.TEXT) WRITE (TEXT, 99009) ID
GO TO 99999
9010 IF(0.LT.TEXT) WRITE (TEXT, 99010) ID
GO TO 99999
9011 IF(0.LT.TEXT) WRITE (TEXT, 99011) ID
GO TO 99999
9012 IF(0.LT.TEXT) WRITE (TEXT, 99012) ID
GO TO 99999
99001 FORMAT
+ (/1X, A9, 'ERROR. THE PRINTING LEVELS ARE NOT ORDERED.'
+ //10X, I10, ' LEVEL(1) FOR MESSAGES'
+ /10X, I10, ' LEVEL(2) FOR THE SOLUTION')
99002 FORMAT
+ (/1X, A9, 'ERROR. THE NUMBER OF PASSES MUST BE POSITIVE.'
+ //10X, I10, ' PASSES')
99003 FORMAT
+ (/1X, A9, 'ERROR. THE NUMBERS OF VARIABLES ARE INCONSISTENT.'
+ //10X, I10, ' SCALARS'
+ /10X, I10, ' COMPONENTS'
+ /10X, I10, ' POINTS'
+ /10X, I10, ' LIMITING NUMBER OF POINTS')
99004 FORMAT
+ (/1X, A9, 'ERROR. THE ARRAY DIMENSIONS ARE MUST BE POSITIVE.'
+ //10X, I10, ' COMPS'
+ /10X, I10, ' PMAX'
+ /10X, I10, ' SCALRS + COMPS * PMAX')
99005 FORMAT
+ (/1X, A9, 'ERROR. GRAB FAILS.')
99006 FORMAT
+ (/1X, A9, 'ERROR. THE INTEGER WORK SPACE IS TOO SMALL.'
+ //10X, I10, ' NEEDED'
+ /10X, I10, ' AVAILABLE')
99007 FORMAT
+ (/1X, A9, 'ERROR. THE REAL WORK SPACE IS TOO SMALL.'
+ //10X, I10, ' NEEDED'
+ /10X, I10, ' AVAILABLE')
99008 FORMAT
+ (/1X, A9, 'ERROR. UNKNOWN TASK.')
99009 FORMAT
+ (/1X, A9, 'ERROR. UNKNOWN REPORT CODE.')
99010 FORMAT
+ (/1X, A9, 'ERROR. NEWTON ERRS.')
99011 FORMAT
+ (/1X, A9, 'ERROR. REFINE ERRS.')
99012 FORMAT
+ (/1X, A9, 'ERROR. TIMSTP ERRS.')
C/// EXIT.
99999 CONTINUE
END
c