cgxbfgr.for
 
c
c
C

File-name TWOPNT.htm 251096

C TWOPNT C C

TWO POINT BOUNDARY VALUE PROBLEM SOLVER

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 C#### JCL 15.02.12 change ITIME(8) to ITIME(20) for consistency with C#### declaration in MCLOCK. Can cause failure on 64-bit C#### INTEGER ITIME(8) 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#### JKW 09.12.93 Changed ELAPSE to TPELPS 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 C#### JKW 09.12.93 Changed name of NORM to TPNORM for PHOENICS 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 C#### JKW 09.12.93 Changed name of NORM to TPNORM for PHOENICS 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 C#### JKW 09.12.93 Changed name of NORM to TPNORM for PHOENICS 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 C#### JKW 09.12.93 Changed name of NORM to TPNORM for PHOENICS 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 C#### JKW 09.12.93 Changed name of NORM to TPNORM for PHOENICS 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#### JKW 09.12.93 Changed name of NORM to TPNORM 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 C#### JKW 09.12.93 Changed name of NORM to TPNORM for PHOENICS 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 C#### JKW 09.12.93 Changed name of NORM to TPNORM for PHOENICS CALL TPNORM (N, CHANGE, BUFFER) ASSIGN 2300 TO ROUTE GO TO 7011 2300 CONTINUE C#### JKW 09.12.93 Changed name of NORM to TPNORM for PHOENICS 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 C#### JKW 09.12.93 Changed ELAPSE to TPELPS + 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. C#### JKW 09.12.93 Changed ELAPSE to TPELPS 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. C#### JKW 09.12.93 Changed ELAPSE to TPELPS 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. C#### JKW 09.12.93 Changed ELAPSE to TPELPS 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. C#### JKW 09.12.93 Changed ELAPSE to TPELPS 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. C#### JKW 09.12.93 Changed ELAPSE to TPELPS 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. C#### JKW 09.12.93 Changed ELAPSE to TPELPS 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 C#### JKW 09.12.93 Changed name of NORM to TPNORM for PHOENICS 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 C#### JKW 09.12.93 Changed ELAPSE to TPELPS 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 C#### JKW 09.12.93 Changed ELAPSE to TPELPS 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) C#### JKW 09.12.93 Changed ELAPSE to TPELPS 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 C#### JKW 09.12.93 Changed ELAPSE to TPELPS 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