C ALGORITHM 681, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 16, NO. 2, PP. 152-157. Installation Distributed with the package are four source code files, two sample input files and the sample output file, and this documentation file. While the input and output files are explained below, the source files consist of 1) INTARI.FOR -- the interval arithmetic subpackage; 2) STKMAN.FOR -- the stack management subpackage; 3) INTBIS.FOR -- the core routines for the package; and 4) LPKD1M.FOR -- the machine constant routine D1MACH and the LINPACK routines used in the package. Installation involves (i) making sure input and output units are correct, and possibly supplying appropriate attachments; (ii) possibly adjusting the size of workspace parameters; (iii) setting a machine-dependent constant in SIMINI and selecting the machine constants routine D1MACH; and (iv) making sure the package has access to the appropriate LINPACK routines. Input occurs from two files, which we term the configuration file and the equation file. Parameter NUNITC in program GENBIS defines the configuration file unit; it is set to unit 20 by default. Parameter NUNITI in GENBIS defines the equation file unit; it is set to 5 by default. All input occurs in GENBIS and INPUT. Output occurs to a single output file, whose unit is defined with parameter NUNITO in GENBIS; the default unit number is 6. All output occurs in GENBIS, OUTPUT, and ERROUT. Workspace parameters all occur in GENBIS. They include the maximum allowable number of variables and equations, the maximum number of possible root-containing boxes, and the maximum height of a stack. The package does error checking for these values, so they need not be altered on the initial installation. The default maximum dimension is 10; the maximum stack height is 200 and the maximum number in the list is 50. With these values, the executable code for an IBM-PC takes about 310K, depending on the compiler. There are five common blocks. Parameters for two of these may need to be changed in several routines if the package is to be used for more than ten equations and unknowns; this should not be too difficult. See details in GENBIS and INPUT. The installer must modify two machine-dependent package routines. The most important of these routines is SIMINI, where data statement for MAXERR may need to be altered. See the in- line documentation. The second routine is ERRHND, in which the system defaults for handling underflows, etc. may need to be changed. The routine D1MACH also must be properly installed and linked with the package. See the documentation within that routine. THE IMPORTANCE OF CORRECT INSTALLATION OF D1MACH AND SIMINI CANNOT BE OVEREMPHASIZED. This is because, if the machine constants are not properly set, then the code could indicate that there are no roots in a region which actually contains roots. One of the authors should be consulted if there are any doubts. INTBIS includes the test input file TOMS86SH.DT1, the test configuration file CONFIG.BIS, and the corresponding output file TOMS86SH.OU1. Results should vary little from machine to machine. This test takes about 4 seconds of CPU time on an IBM-3090. When compiling and linking INTBIS, one may wish to separate the two subroutines FTESTH and HNSNG. Also, the package calls the LINPACK routines DGEFA and DGESL and the LINPACK or BLA routines DAXPY, DCOPY, DDOT, DNRM2, DSCAL, and IDAMAX. Finally, note that the entry point to the main routine is labeled GENBIS; this may be changed since other programs by that name may exist. However, the documentation throughout the package refers to the main routine as GENBIS. A machine-executable code for the IBM-PC/AT compatible family of machines is available on an MS-DOS 360KB diskette from one of the authors. PROGRAM GENBIS C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C Written by: C C R. Baker Kearfott C C and C C Manuel Novoa III C C Department of Mathematics C U.S.L. Box 4-1010 C Lafayette, LA 70504 C C September 29, 1987 C C Part of the generalized bisection package. C C*********************************************************************** C C Function -- C C This is the main driver program, which allocates and indexes C the workspace vectors. It calls the routine ROOTS, which is C structured on the general algorithm in the paper R. B. Kearfott, C 'Abstract generalized bisection and a cost bound', Math. Comput. 49, C 179 (July, 1987), pp. 187-202. C C*********************************************************************** C C Parameter statements -- C PARAMETER (MN = 10) PARAMETER (MMAXDP = 200) PARAMETER (MMAXLS = 50) PARAMETER (MLINFO = 15) C PARAMETER (NUNITI = 5) PARAMETER (NUNITO = 6) PARAMETER (NUNITC = 20) C PARAMETER (MDWORK = 8*MN+9*MN**2) PARAMETER (MIWORK = MN) PARAMETER (MLWORK = 1) C C*********************************************************************** C C Parameter descriptions -- C C MN is the maximum number of variables and equations in any C system processed. C C MMAXDP is the maximum allowable depth of the stack of boxes yet to be C considered by the algorithm. It is also the maximum allowable C depth of the binary search tree. C C MMAXLS MMAXLS - 2 is the maximum number of root-containing boxes which C can be stored. Depending on DLSFLG (see below), this could C include boxes which have been removed due to possible C redundancies. C C MLINFO is the length of the vector INFO, which is set in ROOTS C and which is used in OUTPUT. C C NUNITI is the unit number of the file (device)from which the system, C the initial box, and other initial information will be read. C C NUNITO is the unit number of the file (device) to which all output C will be directed. C C NUNITC is the unit number of the file (device) from which certain C configuration parameters will be read. C C MDWORK is the length of the double precision work vector. C C MIWORK is the length of the integer work vector. C C MLWORK is the length of the logical work vector. C C*********************************************************************** C C Major variable declarations -- C INTEGER NCASE DOUBLE PRECISION CURBOX(2,MN), CURPT(MN) DOUBLE PRECISION STACK(2,MN,MMAXDP) INTEGER STKLVL(MMAXDP) LOGICAL DPIVOT(MMAXDP) INTEGER IPIVOT(MMAXDP) INTEGER TSTITR(MMAXDP) C DOUBLE PRECISION BOXES(2,MN,MMAXLS), POINTS(MN,MMAXLS) INTEGER PTR(2,MMAXLS) LOGICAL UNUSED(MMAXLS) INTEGER BXINFO(5,MMAXLS) LOGICAL PINFO1(MMAXDP,MMAXLS) INTEGER PINFO2(2,MMAXDP,MMAXLS) C INTEGER NINLST, INDBOX(MMAXLS) C DOUBLE PRECISION DWORK(MDWORK) INTEGER IWORK(MIWORK) LOGICAL LWORK(MLWORK) C INTEGER INFO(MLINFO) INTEGER N INTEGER ADJFLG, MAXFT DOUBLE PRECISION EPS, EPSF INTEGER DLSFLG, PRTCON C C*********************************************************************** C C Major variable descriptions (See ROOTS and other routines for more C detailed information.) -- C C NCASE is the number of different problems to be run. (See the C routine INPUT.) C C CURBOX is an array for an interval N-vector (box). C (It is used to pass the initial region to ROOTS.) C C CURPT is space for a point N-vector. C C STACK is an array for the stack of boxes the algorithm must C still process. C C STKLVL, DPIVOT, IPIVOT, and TSTITR are associated with STACK. C The algorithm puts auxiliary information in these arrays. C STKLVL is indexed on the stack depth, while DPIVOT, C IPIVOT, and TSTITR are indexed on the depth of the binary C tree. C C BOXES is an array to contain the lists of root-containing boxes and C possibly other lists of boxes. C C POINTS, PTR, UNUSED, BXINFO, PINFO1, and PINFO2 are associated C with BOXES. The algorithm stores auxiliary information C in these arrays. C C NINLST Upon return from OUTPUT, NINLST is the number of C root-containing boxes in the list. C C INDBOX Upon return from OUTPUT, INDBOX(I) contains the index such C that BOXES(1,J,INDBOX(I)) and BOXES(2,J,INDBOX(I)) are the left C and right endpoints, respectively, of the J-th coordinate C interval of the I-th root-containing box, for I between C 1 and NINLST. C C DWORK is the double precision work array. C C IWORK is the integer work array. C C LWORK is the logical work array. C C INFO contains information about the lists of boxes, etc. C See the documentation in ROOTS, INPUT, OUTPUT, or ERROUT C for more information. C C N is the number of variables and equations in the system C presently under consideration. C C ADJFLG signals whether or not to implement the expansion/deletion C steps (steps 4 of the 1987 Math. Comput. paper (ibid.)): C C If ADJFLG = 1 then expansion/deletion steps are done. C C If ADJFLG is not equal to 1, then expansion/deletion steps C are omitted. C C MAXFT is the maximum allowable number of calls to the root-inclusion C test (FTEST). C C EPS is the minimum allowable width (in the norm defined by DIAMCP) C of an unresolved box. C C EPSF is the range tolerance; if it is certain that each component C of the vector function has absolute value less than EPSF C within a box, then that box is signalled as containing a C root. C C DLSFLG signals whether or not boxes deleted from the list during C the expansion/deletion process are stored in another C list. C C If DLSFLG = 1 then such boxes are stored. C C If DLSFLG is not equal to 1, then such boxes are NOT stored. C C (This flag is only significant if ADJFLG = 1. C C PRTCON controls the amount of information the routine OUTPUT C prints. (See OUTPUT for information.) C C IPTFMT controls the format for printing of floating point C numbers in the OUTPUT routine. See the OUTPUT routine C for details. C (OUTPUT) C C*********************************************************************** C C Common block declarations -- none C C However, users should be aware of the common blocks EQUAT, COEFFS, C CONFG1, and CONFG2 in INPUT. EQUAT and COEFFS also occur in C POLFUN, POLJAC, POLFSC, AND SCLFSC, while CONFG1 and CONFG2 also occur C in HNSNG. The block CONFG2 also occurs in CHKLST and DELLST, while C the block CONFG1 also occurs in PVSLCT. C Additionally the routines SIMINI and RNDOUT share machine constants C in common block MACH1. C C*********************************************************************** C C Fortran-supplied functions and subroutines -- none C C*********************************************************************** C C Package-supplied functions and subroutines -- C C ERRHND, ERROUT, INPUT, OUTPUT, ROOTS, SIMINI C EXTERNAL FTESTH, POLFSC, POLFUN, POLJAC, POLJSC C C*********************************************************************** C C User-supplied functions and subroutines -- none C C*********************************************************************** C C I/O functions -- C C This routine gets the number of problems to be solved from UNITI C in I5 format. C C This routine outputs a message concerning which interval arithmetic C routines it employs to unit UNITO. C C*********************************************************************** C C Internal constant declarations -- C INTEGER MAXDP, MAXLST, UNITI, UNITO, UNITC INTEGER LINFO, LNDWRK, LNIWRK, LNLWRK C C*********************************************************************** C C Internal constant descriptions -- C C MAXDP has the same value as parameter MMAXDP. C C MAXLST has the same value as parameter MMAXLS. C C UNITI has the same value as parameter NUNITI. C C UNITO has the same value as parameter NUNITO. C C UNITC has the same value as parameter NUNITC. C C LINFO has the same value as parameter MLINFO. C C LNDWRK has the same value as parameter MDWORK. C C LNIWRK has the same value as parameter MIWORK. C C LNLWRK has the same value as parameter MLWORK. C C*********************************************************************** C C Beginning of executable statements -- C MAXDP = MMAXDP MAXLST = MMAXLS C LNDWRK = MDWORK LNIWRK = MIWORK LNLWRK = MLWORK C LINFO = MLINFO C UNITI = NUNITI UNITO = NUNITO UNITC = NUNITC C C Call the system-dependent error handling setup routine. C CALL ERRHND C C*********************************************************************** C C If the user has supplied the interval arithmetic routines, then C the following lines of code should not be commented. They should C otherwise be commented. C C WRITE(UNITO,1030) C WRITE(UNITO,1030) 'The user has supplied the interval arithmetic' C WRITE(UNITO,1030) 'routines.' C WRITE(UNITO,1030) C C When using the interval arithmetic routines supplied with the C package, the following lines of code should not be commented. C They should otherwise be commented. C CALL SIMINI WRITE(UNITO,1030) WRITE(UNITO,1030) 'The interval arithmetic routines supplied' WRITE(UNITO,1030) 'with this package will be used.' WRITE(UNITO,1030) C C*********************************************************************** C C The following dummy assignment is required for systems which C require adjustable array dimensions to be assigned before C entry into subroutines -- C N = 1 C C Input the number of problems to be solved. C READ(UNITI,*) NCASE C C Input the system, the initial region, etc. C DO 10 I = 1,NCASE C CALL INPUT(N,MN,UNITI,UNITO,UNITC,CURBOX,EPS,EPSF,ADJFLG,MAXFT, 1 DLSFLG,PRTCON,IPTFMT,INFO(1),INFO(2)) IF (INFO(1).NE.0) GOTO 800 C C Call the main generalized bisection algorithm. C CALL ROOTS(N,MAXLST,MAXDP,MAXFT,EPS,EPSF, 1 ADJFLG,DLSFLG,FTESTH,POLFUN,POLJAC, 2 CURBOX,CURPT,LINFO,INFO, 3 STACK,STKLVL,DPIVOT,IPIVOT,TSTITR, 4 BOXES,POINTS,PTR,UNUSED,BXINFO,PINFO1,PINFO2, 5 LNDWRK,DWORK,LNIWRK,IWORK,LNLWRK,LWORK) C C Take appropriate action if an error occurred during the root-finding C process. C 800 CALL ERROUT(N,MAXDP,MAXLST,MAXFT,UNITO,CURBOX,CURPT,LINFO, 1 INFO,STACK,STKLVL,DPIVOT,IPIVOT,TSTITR) C C Output the results (roots and diagnostics). C CALL OUTPUT(POLFUN,POLFSC,POLJSC,N,MAXDP,MAXLST,MAXFT, 1 EPS,EPSF,PRTCON,IPTFMT,UNITO,ADJFLG,DLSFLG,LINFO, 2 INFO,BOXES,POINTS,PTR,UNUSED,BXINFO,PINFO1,PINFO2, 3 NINLST,INDBOX,LNDWRK,DWORK,LNIWRK,IWORK) C 10 CONTINUE C STOP C 1030 FORMAT(1X,A,I5) C END C*********************************************************************** C*********************************************************************** SUBROUTINE ERRHND C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C Written by: C C R. Baker Kearfott C C and C C Manuel Novoa III C C Department of Mathematics C U.S.L. Box 4-1010 C Lafayette, LA 70504 C C September 29, 1987 C C Part of the generalized bisection package C C*********************************************************************** C C Called by -- GENBIS C C*********************************************************************** C C Function -- C C This machine-dependent routine is meant to access operating system C parameters which determine handling of arithmetic errors. C In particular -- C C (1) Underflows should be allowed. C C*********************************************************************** C C FORTRAN-SUPPLIED FUNCTIONS AND SUBROUTINES -- C C FOR THE IBM VM/CMS ENVIRONMENT -- C C EXTERNAL ERRSET C C*********************************************************************** C C Beginning of executable statements -- C C Set error handling for exponent underflow to no messages printed and C an unlimited number of occurrences. C C statement below is for IBM VM/CMS -- C C CALL ERRSET(208,300,-1,1,0,0) C RETURN END C*********************************************************************** C*********************************************************************** SUBROUTINE INPUT(N,MN,UNITI,UNITO,UNITC,BOX,EPS,EPSF,ADJFLG,MAXFT, 1 DLSFLG,PRTCON,IPTFMT,ERRFLG,ERRVAL) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C PARAMETER (MN2 = 10) PARAMETER (MT = 30) C C Written by: C C R. Baker Kearfott C C and C C Manuel Novoa III C C Department of Mathematics C U.S.L. Box 4-1010 C Lafayette, LA 70504 C C Formats are taken partially from Alexander Morgan. (See below). C C September 29, 1987 C C Part of the generalized bisection package C C*********************************************************************** C C Called by -- GENBIS C C*********************************************************************** C C Function -- C C This routine reads in the dimension and the coefficients of C the polynomial system from the main input file. It also reads C in tolerances, parameters which define the algorithm, and flags C which control printing from a separate configuration file. C C Depending on a configuration parameter, this routine may echo C the input. The routine always outputs a title associated with C the problem. C C This routine performs some error checking, and sets error flags C as appropriate. C C*********************************************************************** C C Argument declarations -- C INTEGER N, MN INTEGER UNITI, UNITO, UNITC DOUBLE PRECISION BOX(2,N), EPS, EPSF INTEGER ADJFLG, MAXFT, DLSFLG, PRTCON, IPTFMT, ERRFLG, ERRVAL C C*********************************************************************** C C Argument descriptions -- (INPUT = set on entry and not alterable) C (OUTPUT = to be set by the routine) C (I/O = set on entry but alterable) C C N is the actual number of equations and variables. C (OUTPUT) C C MN is the maximum number of equations as defined in the C parameter statement in GENBIS. C (INPUT) C C UNITI is the unit number of the file (device) from which the C coefficients and exponents of the equations and the C endpoint intervals of the initial box will be read. C (INPUT) C C UNITO is the unit number of the file (device) to which all output C is directed. C (INPUT) C C UNITC is the unit number of the file (device) from which C tolerances and certain other configuration parameters will c be read. C (INPUT) C C BOX will contain the initial region to be examined. C (an N-dimensional interval vector) C (OUTPUT) C C EPS is the domain tolerance, and is related to the smallest C box to be produced by bisection. C (OUTPUT) C C EPSF is the range tolerance. C (OUTPUT) C C ADJFLG is set to 1 if expansion steps are to be implemented, and C is set to some other value otherwise. (See Algorithm 2.5 C in the paper R. B. Kearfott, 'Abstract generalized C bisection and a cost bound', Math. Comput. 49, 179 C (July, 1987), pp. 187-202. C (OUTPUT) C C MAXFT is the maximum number of calls to the routine which C tests boxes for roots. C (OUTPUT) C C DLSFLG is only used if expansion steps are implemented. If it C is set to 1, the algorithm produces a list of possible C root-containing boxes which have been removed from the C primary list because they are adjacent to or intersect C other such boxes. C (OUTPUT) C C PRTCON is a print control flag for the OUTPUT routine. Check C See the OUTPUT routine documentation for details. C (OUTPUT) C C IPTFMT controls the format for printing of floating point C numbers in the OUTPUT routine. See the OUTPUT routine C for details. C (OUTPUT) C C ERRFLG signals errors for the error-printing routine ERROUT. C If ERRFLG = 0 on return, then no monitored error has C occurred. C (OUTPUT) C C ERRVAL will contain additional information about the error, C if one has occurred. C (OUTPUT) C C*********************************************************************** C C Internal variable declarations -- C CHARACTER*1 TITLE(60) C INTEGER IECHO INTEGER J, K, L, NT C C*********************************************************************** C C Internal variable descriptions -- C C TITLE identifies the input file which has been accessed. C C IECHO controls echoing of input data within this routine. (See C the description of the configuration file below.) C C J, K, and L are loop indices. C C NT temporarily stores the number of terms in equations. C C*********************************************************************** C C Common block declarations -- C C In the common blocks below, MN2 denotes the number of equations, C and MT denotes the maximum number of terms in any single C equation. C INTEGER NUMT, KDEG COMMON/EQUAT/NUMT(MN2),KDEG(MN2,MN2,MT) C DOUBLE PRECISION A COMMON/COEFS/A(MN2,MT) C C /EQUAT/ passes exponents of variables from INPUT to C POLFUN, POLJAC, POLFSC, AND POLJSC. C C /COEFS/ passes coefficient values from INPUT to POLFUN, C POLJAC, POLFSC, AND POLJSC. C INTEGER ITRFLG, PIVFLG, JACFLG COMMON /CONFG1/ ITRFLG, PIVFLG, JACFLG C DOUBLE PRECISION TOL1, TOL2, TOL3 COMMON /CONFG2/ TOL1, TOL2, TOL3 C C /CONFG1/ passes control flags to the interval Newton C method routine HNSNG and to the routine PVSLCT. C These flags are described below with the C configuration file. C C /CONFG2/ passes tolerances to the interval Newton method C routine HNSNG and to the routines CHKLST and C DELLST. These tolerances are described below C with the configuration file. C C*********************************************************************** C C Fortran-supplied functions and subroutines -- none C C*********************************************************************** C C Package-supplied functions and subroutines -- none C C*********************************************************************** C C User-supplied functions and subroutines -- none C C*********************************************************************** C C I/O functions -- C C INPUT: Is read from the configuration file UNITC and the C main input file UNITI. C C INPUT from UNITC: C C Input of the following variables occurs in the following order. C C DLSFLG -- See argument descriptions above. C C PRTCON -- See argument descriptions above. C C ITRFLG -- This is presently unused, but may control iteration C of the interval Newton method, if that routine is C modified in the future. C C PIVFLG -- This controls the way the index of the coordinate C interval to be bisected is chosen. C C If PIVFLG = 0, then the coordinate interval of C maximum width is chosen. C C If PIVFLG = 1, then the coordinate interval of C maximum width according to a certain scaling scheme C is chosen. C C If PIVFLG = 2, then the coordinate interval C corresponding to maximum 'smear' in the interval C Jacobian matrix is chosen. (See the explanation in the C paper accompanying this program.) C C JACFLG -- This controls the way the preconditioner matrix Y for C the interval Newton method is computed. If A is an C approximation to the value of the Jacobian matrix at C the midpoint of the box, then we mathematically C set Y to be the inverse of A. If F' is the interval C Jacobian matrix, then Y F' is the preconditioned C interval Jacobian matrix. If we numerically compute C Y as ( A-transpose)**(-1) )-transpose, then the C residual Y A - I may be smaller than if we compute Y C directly as A**(-1). C C If JACFLG = 0 then we compute Y in the ordinary way. C C If JACFLG = 1 then we compute Y from the transpose. C C ADJFLG -- See argument descriptions above. C C TOL1 -- This controls iteration of the interval Newton method. C If the ratio of volumes of the new box to the old C box is less than TOL1 but greater than TOL2, then C routine HNSNG does another iteration of the interval C Newton method after re-evaluating the Jacobian matrix. C C TOL2 -- If the ration of volumes of the new box to the old C box in the interval Newton method is less than TOL2, C then routine HNSNG does another iteration of the C interval Newton method with the same interval C Jacobian matrix as before. C C TOL3 -- TOL3 is presently unused, but is available to C HNSNG, CHKLST, and DELLST. C C EPS -- is the domain tolerance. C C EPSF -- is the range tolerance. C C MAXFT -- is the maximum allowable number of calls to the C root-inclusion test routine HNSNG. C C IECHO -- If IECHO = 0, then this routine does not echo any C of the input. C C If IECHO = 1, then this routine echos the configuration C parameters to the output unit UNITO. C C If IECHO = 2, then this routine echos both the C configuration parameters and the polynomial coefficients C and exponents to the output unit UNITO. C C IPTFMT -- See argument descriptions above. C C A sample UNITI input file follows. (Delete the first column to use C as a template.) C C 1 DLSFLG (0 no del. list, 1 del. list) C 2 PRTCON (1 ordinary, 2 ext., 3 1+path) C 1 ITRFLG (presently unused) C 2 PIVFLG (0 wid., 1 scl. wid., 2 max. sm.) C 1 JACFLG (0 ordinary, 1 transpose) C 1 ADJFLG (0 no exp./del., 1 exp./del.) C .6d0 TOL1 (good vol. ratio. in iteration) C -.4d0 TOL2 (extremely good vol. rat. in it.) C 0d0 TOL3 (presently unused) C 1.0D-05 EPSMIN (domain tolerance times 16) C 1.0D-10 EPSF (range tolerance) C 10000 MAXFT (max. no. calls to incl.test) C 2 IECHO (0 no echo, 1 config., 2 all) C 1 IPTFMT (0 lg, 1 med, 2 sht, 3 132) C C INPUT from UNITI: C C TITLE -- a 60 character alphanumeric problem description. C C N -- the number of equations and variables. C C The following three sets of input variables are grouped by equation C and term, as in the sample below: C C NUMTRM(I) -- number of terms in the I-th equation, C I=1 to N. C C DEG(I,K,J) -- degree of the K-th variable in the J-th term C of the I-th equation, I=1 to N, K=1 to N, J=1 to C NUMTRM(I). C C A(I,J) -- coefficient of the J-th term of the I-th equation, C I=1 to N, J=1 to NUMTRM(I). C C A sample UNITI input file follows. (Delete the first column to use C as a template.) Note that the first record contains the number of C datasets, and is read in the calling program. The remaining records C of this data file contain a sample data set. If NSETS > 1, then C additional data sets would be placed directly after the first one C (without spaces), and would follow the same format. C C 1 NSETS CTITLE: CUBIC - PARABOLA (TOMS 1986 -- PROBLEM 1)......... C 2 N C 3 NUMTRM(1) C 3 DEG(1,1,1) C 0 DEG(1,2,1) C 4.0D 00 A(1,1) C 1 DEG(1,1,2) C 0 DEG(1,2,2) C -3.0D 00 A(1,2) C 0 DEG(1,1,3) C 1 DEG(1,2,3) C -1.0D 00 A(1,3) C 2 NUMTRM(2) C 2 DEG(2,1,1) C 0 DEG(2,2,1) C 1.0D 00 A(2,1) C 0 DEG(2,1,2) C 1 DEG(2,2,2) C -1.0D 00 A(2,2) C -2.0D0 D0(1,1) C 2.0D0 D0(2,1) C -2.0D0 D0(1,2) C 2.0D0 D0(2,2) C C The last 2*N lines of this input file represent the endpoints C of the intervals which define the initial box. The other lines C of the input file are identical in format and content to the C lines of the input file for the program INPTAT on p. 366 of C Alexander Morgan, Solving Polynomial Systems using Continuation C for Engineering and Scientific Problems, Prentice-Hall, Englewood C Cliffs, NJ, 1987. Morgan refers to this as 'tableau' input of C the polynomial system. C C OUTPUT: If PRTCON > 0, then TITLE is output to UNITO. The input C and configuration parameters are also possibly echoed C to UNITO, depending on the value of IECHO from the C configuration file. C C*********************************************************************** C C INTERNAL CONSTANT DECLARATIONS -- NONE C C*********************************************************************** C C Beginning of executable statements -- C C Initialize error flag to 'no error'. C ERRFLG = 0 C C Check for compatibility between MN and common block dimensions. C C Input configuration parameters. C READ(UNITC,1040) DLSFLG READ(UNITC,1040) PRTCON READ(UNITC,1040) ITRFLG READ(UNITC,1040) PIVFLG READ(UNITC,1040) JACFLG READ(UNITC,1040) ADJFLG READ(UNITC,1030) TOL1 READ(UNITC,1030) TOL2 READ(UNITC,1030) TOL3 READ(UNITC,1030) EPS READ(UNITC,1030) EPSF READ(UNITC,1040) MAXFT READ(UNITC,1040) IECHO READ(UNITC,1040) IPTFMT C REWIND UNITC C C INPUT STATEMENTS FOR THE SYSTEM FOLLOW. C READ(UNITI,1010) TITLE READ(UNITI,1020) N C IF (N .GT.MN2) THEN ERRFLG = 14 ERRVAL = N RETURN END IF C IF (N.GT.MN) THEN ERRFLG = 3 ERRVAL = N RETURN END IF C DO 30 J= 1, N READ(UNITI,1020) NUMT(J) NT = NUMT(J) IF (NT .GT. MT) THEN ERRFLG = 15 ERRVAL = J RETURN END IF DO 20 K = 1, NT DO 10 L = 1, N READ(UNITI,1020) KDEG(J,L,K) 10 CONTINUE READ(UNITI,1030) A(J,K) 20 CONTINUE 30 CONTINUE C C*********************************************************************** C C Input the initial box. C C DO 40 J = 1, N READ(UNITI,1030) BOX(1,J) READ(UNITI,1030) BOX(2,J) 40 CONTINUE C C Output section. C IF (PRTCON.GT.0 .OR. IECHO.GE.1) THEN WRITE(UNITO,2030) 1 '************************************************************' WRITE(UNITO,2050) WRITE(UNITO,2050) WRITE(UNITO,2050) WRITE(UNITO,2030) 'Generalized bisection package' WRITE(UNITO,2040) TITLE END IF C IF (IECHO.GE.2) THEN WRITE(UNITO,2030) 'Number of equations = ', N WRITE(UNITO,2050) C WRITE(UNITO,2020) 'Initial region:',((BOX(I,J),I=1,2),J=1,N) WRITE(UNITO,2050) DO 70 J = 1, N WRITE(UNITO,2050) WRITE(UNITO,2050) WRITE(UNITO,2060) J,NUMT(J) NT = NUMT(J) DO 60 K = 1, NT DO 50 L = 1, N WRITE(UNITO,2070) J, L, K, KDEG(J,L,K) 50 CONTINUE WRITE(UNITO,2080) J, K, A(J,K) 60 CONTINUE 70 CONTINUE WRITE(UNITO,2050) WRITE(UNITO,2050) END IF C IF (IECHO.GE.1) THEN WRITE(UNITO,2050) WRITE(UNITO,2010) 'CONFIGURATION PARAMETERS:' WRITE(UNITO,2050) IF (ADJFLG.EQ.1) THEN WRITE(UNITO,2010) 'Expansion steps are implemented.' ELSE WRITE(UNITO,2010) 'Expansion steps are not implemented.' END IF WRITE(UNITO,2050) WRITE(UNITO,2030) 'ITRFLG: ',ITRFLG WRITE(UNITO,2030) 'PIVFLG: ',PIVFLG WRITE(UNITO,2030) 'JACFLG: ',JACFLG WRITE(UNITO,2010) 'TOL1 : ',TOL1 WRITE(UNITO,2010) 'TOL2 : ',TOL2 WRITE(UNITO,2010) 'TOL3 : ',TOL3 WRITE(UNITO,2010) 'EPS : ',EPS WRITE(UNITO,2010) 'EPSF : ',EPSF WRITE(UNITO,2030) 'MAXFT : ',MAXFT WRITE(UNITO,2030) 'IECHO : ',IECHO WRITE(UNITO,2030) 'IPTFMT: ',IPTFMT WRITE(UNITO,2030) 'ADJFLG: ',ADJFLG WRITE(UNITO,2050) WRITE(UNITO,2050) END IF C C*********************************************************************** C C Input formats -- C 1010 FORMAT(60A1) 1020 FORMAT(I5) 1030 FORMAT(D22.15) 1040 FORMAT(I8) C C Output formats -- C 2010 FORMAT(1X,A,1X,D30.22) 2020 FORMAT(1X,A,50(/2(2X,D30.22))) 2030 FORMAT(1X,A,I8) 2040 FORMAT(1X,60A1) 2050 FORMAT(1X) 2060 FORMAT(' NUMT(',I2,')=',I5) 2070 FORMAT(' KDEG(',I2,',',I2,',',I2,')=',I5) 2080 FORMAT(' A(',I2,',',I2,')=',D22.15) C RETURN END C*********************************************************************** C*********************************************************************** SUBROUTINE OUTPUT(FUNC,FUNCSC,JACSC,N,MAXDP,MAXLST,MAXFT, 1 EPS,EPSF,PRTCON,IPTFMT,UNITO,ADJFLG,DLSFLG,LINFO,INFO,BOXES, 2 POINTS,PTR,UNUSED,BXINFO,PINFO1,PINFO2,NINLST,INDBOX, 3 LNDWRK,DWORK,LNIWRK,IWORK) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C Written by: C C R. Baker Kearfott C C and C C Manuel Novoa III C C Department of Mathematics C U.S.L. Box 4-1010 C Lafayette, LA 70504 C C September 29, 1987 C C Part of the generalized bisection package C C*********************************************************************** C C Called by -- GENBIS C C*********************************************************************** C C Function -- C C This is the primary output routine for the generalized bisection C package. Most output except for error messages occurs from this C routine. C C*********************************************************************** C C Argument declarations -- C EXTERNAL FUNC, FUNCSC, JACSC INTEGER N, MAXDP, MAXLST, MAXFT INTEGER PRTCON, IPTFMT, UNITO INTEGER ADJFLG, DLSFLG, LINFO, INFO(LINFO) DOUBLE PRECISION BOXES(2,N,MAXLST), POINTS(N,MAXLST) INTEGER PTR(2,MAXLST) LOGICAL UNUSED(MAXLST) INTEGER BXINFO(5,MAXLST) LOGICAL PINFO1(MAXDP,MAXLST) INTEGER PINFO2(2,MAXDP,MAXLST) INTEGER NINLST, INDBOX(MAXLST) INTEGER LNDWRK DOUBLE PRECISION DWORK(LNDWRK) INTEGER LNIWRK INTEGER IWORK(LNIWRK) C C*********************************************************************** C C Argument descriptions -- (INPUT = set on entry and not alterable) C (OUTPUT = to be set by the routine) C (I/O = set on entry but alterable) C C FUNC is the name of the external routine which performs interval C function evaluations. GENBIS uses POLFUN. C C FUNCSC is the name of an external routine which does point function C evaluations. GENBIS uses POLFSC. C C JACSC is the name of an external routine which dies point Jacobian C matrix evaluations. GENBIS uses POLJSC. C C N is the number of equations and variables. C (INPUT) C C MAXDP is the maximum allowable depth of the stack of boxes yet to be C considered. It is also the maximum allowable depth of the C binary search tree. C (INPUT) C C MAXLST MAXLST - 2 is the maximum number of root-containing boxes which C can be stored. C (INPUT) C C MAXFT is the maximum allowable number of calls to the root-inclusion C test (FTESTH). C (INPUT) C C EPS is the minimum allowable width (in the norm defined by DIAMCP) C of an unresolved box. It is also the domain tolerance in C the scalar Newton method. C (INPUT) C C EPSF is the range tolerance; if it is certain that each component C of the vector function has absolute value less than EPSF C within a box, then that box is signalled as containing a C root. It is also the range tolerance in the scalar Newton C method. C (INPUT) C C PRTCON controls the amount of information printed in this routine. C C If PRTCON = 0, then this routine does not print; it merely C computes the number of boxes in the linked C list and returns their indices. C C If PRTCON = 1, then this routine prints the C root-containing boxes, approximate roots, C residuals, and a list of boxes which were C removed from the list due to possible C redundancies. C C If PRTCON = 2, then this routine also prints the boxes' C pointers in the linked list, the reasons C the boxes were added to the list, etc. C C If PRTCON = 3, then this routine also prints a table of C which boxes intersected boxes in the list, C as well as information about which coordinates C were bisected and the number of iterations of C the interval Newton method. C C (INPUT) C C IPTFMT controls the format for printing of floating point C numbers in this routine. C C If IPTFMT = 0, then floating point results are printed C with 22 digits in 80 column format. C C If IPTFMT = 1, then floating point results are printed C with 9 digits in 80 column format. C C If IPTFMT = 2, then floating point results are printed C with 2 digits in 80 column format. C C If IPTFMT = 3, then floating point results are printed C with 3 digits in 132 column format. This C setting is most appropriate for larger C systems. C C (INPUT) C C UNITO is the unit number of the file (device) to which all output C is directed. C (INPUT) C C ADJFLG signals whether or not to implement the expansion/deletion C steps (steps 4 of the 1987 Math. Comput. paper (ibid.)): C C If ADJFLG = 1 then expansion/deletion steps are done. C C If ADJFLG is not equal to 1, then expansion/deletion steps C are omitted. C C (INPUT) C C DLSFLG signals whether or not boxes deleted from the list during C the expansion/deletion process are stored in another C list. C (INPUT) C C LINFO is the length of the information vector INFO. C (INPUT) C C INFO is a vector containing information about the lists, as well as C other statistics: C C INFO(1) = error type. (ERRFLG : See subroutine ERROUT.) C C INFO(2) = additional error information (ERRVAL). C C INFO(3) = current depth of the stack. C C INFO(4) = the number of boxes which have been examined. C C INFO(5) = the number of boxes in list 1 (which contains C root-containing boxes). C C INFO(6) = the number of boxes in list 2, if it exists. C (See the explanation of DLSFLG above). C Otherwise, it is the number of boxes deleted in C expansion/deletion steps. C C INFO(7) = the header note for list 1. C C INFO(8) = the header node for list 2. C C INFO(9) = the maximum depth reached in the stack. C C INFO(10) = the maximum level reached in the binary tree. C C INFO(11) = the number of calls to the root inclusion test. C C INFO(12) = the number of interval function calls made to C possibly determine that the zero vector is not C in the range. C C INFO(13) = the number of function calls made for the C interval Newton method. C C INFO(14) = the number of interval and scalar Jacobian calls. C C INFO(15) = the number of expanded boxes. C C (INPUT) C C BOXES contains the lists of root-containing boxes, boxes which C have been deleted due to redundancies, and possibly other C lists of boxes. C (INPUT) C C POINTS contains points corresponding to the boxes in BOXES. These C points may simply be midpoints of the boxes, or approximations C to roots. C (I/O) C C PTR contains information about the 'links' between nodes in the C linked lists in BOXES and POINTS. C (I/O) C C UNUSED UNUSED(I) contains the state of the storage areas BOXES(*,*,i) C and POINTS(*,*,I); UNUSED(I) = 'TRUE' means that these C areas are not in use. C (I/O) C C BXINFO, PINFO1, and PINFO2 are storage areas whose entries C correspond to entries in BOXES, POINTS, and UNUSED. C Additional information about the root-containing boxes is C placed in these arrays. C (I/O) C C NINLST Upon return from OUTPUT, NINLST is the number of C root-containing boxes in the list. C (OUTPUT) C C INDBOX Upon return from OUTPUT, INDBOX(I) contains the index such C that BOXES(1,J,INDBOX(I)) and BOXES(2,J,INDBOX(I)) are the left C and right endpoints, respectively, of the J-th coordinate C interval of the I-th root-containing box, for I between C 1 and NINLST. C (OUTPUT) C C LNDWRK is the dimension of the double precision work vector. C (INPUT) C C DWORK is the double precision work vector. C (INPUT) C C LNIWRK is the length of the integer work vector. C (INPUT) C C IWORK is the integer work vector. C C*********************************************************************** C C Internal variable declarations -- C INTEGER ERRFLG, ERRVAL, DEPTH INTEGER NLEAVE, NLIST1, NLIST2, LIST(2), DMAX, LMAX INTEGER NFTCAL, NFCAL1, NFCAL2, NJCALL, NADJ DOUBLE PRECISION RESID INTEGER I, J, K, COUNT INTEGER CURPTR, BOXPTR LOGICAL PTFLG INTEGER IERR INTEGER IFINT, ISCLF, IDEND, IJAC, ISCX, IV INTEGER IIPVT, IIEND C INTEGER I1000, I1010, I1020, I1030, I1040, I1050, I1060 C C*********************************************************************** C C Internal variable descriptions -- C C ERRFLG Error type: 0 if no error, another value otherwise C (See subroutine ERROUT for more detailed information.) C C ERRVAL Additional error information C (See subroutine ERROUT for more detailed information.) C C DEPTH The current depth of the stack of boxes yet to be examined. C C NLEAVE The number of boxes for which the algorithm has reached a C conclusion about whether they contain roots. C C NLIST1 The number of boxes in list 1 C (List 1 is the list of root-containg boxes.) C C NLIST2 The number of boxes in list 2, if it exists; otherwise the C number of boxes deleted in expansion/deletion steps C (See the description of argument DLSFLG.) C C LIST The pointers to the header nodes of the linked lists C C DMAX The maximum depth reached in the stack of boxes yet to be C considered C C LMAX The maximum level reached in the binary search tree C corresponding to the bisection process C C NFTCAL The number of calls made to FTEST C C NFCAL1 the number of interval function calls made to possibly C determine that the zero vector is not in the range C C NFCAL2 the number of function calls made for the interval Newton C method C C NJCALL the number of interval and scalar Jacobian calls C C NADJ the number of expanded boxes C C RESID a temporary variable C C I, J, and K are loop indices. C C COUNT is an index for the boxes in a list as they are printed. C C CURPTR and BOXPTR are temporary pointers. C C PTFLG is set to .TRUE. if the intersection of two boxes has zero C volume, and is set to .FALSE. otherwise. C C IERR indicates if the classical Newton method has converged, C when precisely finding a root in a root-containing box. C C IFINT, ISCLF, IDEND, IJAC, ISCX, and IV are pointers into the C double precision work array. C C IIPVT and IIEND are pointers into the integer work array. C C I1000, I1010, I1020, I1030, I1040, I1050, and I1060 are variable C format statement numbers. I1010 corresponds to floating C point interval output, while I1020 corresponds to C non-interval floating point output. C C*********************************************************************** C C Common block declarations -- none C C*********************************************************************** C C Fortran-supplied functions and subroutines -- none C C*********************************************************************** C C Package-supplied functions and subroutines -- C C FUNC (POLFUN), FUNCSC (POLFSC), C JACSC (POLJSC), NEWTON C C DCOPY, DNRM2 (from LINPACK) C C*********************************************************************** C C User-supplied functions and subroutines -- none C C*********************************************************************** C C I/O functions -- C C Output is to unit UNITO. See the explanations of PRTCON and IPTFMT C above. C C*********************************************************************** C C Internal constant declarations -- C INTEGER EOLIST, PREV, NEXT LOGICAL LEFT, RIGHT C C*********************************************************************** C C Internal constant descriptions -- C C EOLIST THE END-OF-LIST VALUE (0) FOR THE LINKED LIST STRUCTURE C C PREV In a linked list, PTR(PREV,*) is the pointer to the box C preceeding the box with pointer '*'. C C NEXT In a linked list, PTR(NEXT,*) is the pointer to the box C following the box with pointer '*'. C C LEFT is the logical value associated in this program with the C left box formed from bisection. The left box is that box C formed by replacing the right endpoint of the bisected C interval by its midpoint.) C C RIGHT is the logical value associated in this program with the C right box formed from bisection. The right box is that box C formed by replacing the left endpoint of the bisected C interval by its midpoint.) C C*********************************************************************** C DATA EOLIST/0/, PREV/1/, NEXT/2/ DATA LEFT/.TRUE./, RIGHT/.FALSE./ C C Beginning of executable statements -- C C Initialize pointers, etc. C IFINT = 1 ISCLF = IFINT + 2*N IJAC = ISCLF + N IDEND = IJAC + N*N -1 C ISCX = IFINT IV = ISCX + N C IIPVT = 1 IIEND = IIPVT + N-1 C ASSIGN 1000 TO I1000 ASSIGN 1030 TO I1030 ASSIGN 1040 TO I1040 ASSIGN 1050 TO I1050 ASSIGN 1060 TO I1060 IF (IPTFMT .EQ. 0) THEN ASSIGN 1020 TO I1010 ASSIGN 1020 TO I1020 ELSE IF (IPTFMT .EQ. 1) THEN ASSIGN 1011 TO I1010 ASSIGN 1021 TO I1020 ELSE IF (IPTFMT .EQ. 2) THEN ASSIGN 1012 TO I1010 ASSIGN 1022 TO I1020 ELSE IF (IPTFMT .EQ. 3) THEN ASSIGN 1013 TO I1010 ASSIGN 1023 TO I1020 END IF C ERRFLG = INFO(1) ERRVAL = INFO(2) DEPTH = INFO(3) NLEAVE = INFO(4) NLIST1 = INFO(5) NLIST2 = INFO(6) LIST(1) = INFO(7) LIST(2) = INFO(8) DMAX = INFO(9) LMAX = INFO(10) NFTCAL = INFO(11) NFCAL1 = INFO(12) NFCAL2 = INFO(13) NJCALL = INFO(14) NADJ = INFO(15) C C*********************************************************************** C C Output statistics. C IF (PRTCON.GT.0) THEN WRITE(UNITO,I1030) 'After ROOTS:' WRITE(UNITO,I1030) ' Number of roots found:',NLIST1 WRITE(UNITO,I1030) ' Number of leaves:',NLEAVE IF (ADJFLG.EQ.1) THEN WRITE(UNITO,I1030) ' Number of adjusted boxes:',NADJ WRITE(UNITO,I1030) 1 ' Number of boxes deleted due to redundancy:',NLIST2 ELSE WRITE(UNITO,I1030) 1 'Algorithm was configured to not adjust or delete boxes.' WRITE(UNITO,I1030) 'Number of undetermined boxes:', NLIST2 END IF WRITE(UNITO,I1030) ' Maximum binary tree level reached:',LMAX WRITE(UNITO,I1030) ' Maximum stack depth reached:',DMAX WRITE(UNITO,I1030) ' Number of ftest calls:',NFTCAL WRITE(UNITO,I1030) 1 ' Number of interval function calls to determine the range:', 2 NFCAL1 WRITE(UNITO,I1030) ' Number of Jacobian calls:', NJCALL WRITE(UNITO,I1030) 1 ' Number of function calls to do the interval Newton method:' 2 , NFCAL2 WRITE(UNITO,I1000) WRITE(UNITO,I1000) END IF C C Output boxes. C DO 200 I = 1, 2 C IF ((I.EQ.1).AND.(NLIST1.EQ.0)) THEN IF (ADJFLG.EQ.1) THEN IF (PRTCON.GT.0) WRITE(UNITO,I1030) 1 'No roots were found.' GOTO 200 ELSE IF (PRTCON.GT.0) THEN WRITE(UNITO,I1030) 1 'The test signalled no boxes as root-containing.' END IF GOTO 200 END IF ELSE IF (I.EQ.1) THEN IF (PRTCON.GT.0) WRITE(UNITO,I1030) 1 'LIST OF ROOT-CONTAINING BOXES FOLLOWS:' ELSE IF ((I.EQ.2).AND.(NLIST2.EQ.0).AND.(ADJFLG.EQ.1)) THEN IF (PRTCON.GT.0) THEN WRITE(UNITO,I1000) WRITE(UNITO,I1000) WRITE(UNITO,I1030) 'No boxes were deleted from the list.' END IF GOTO 200 ELSE IF ((I.EQ.2).AND.(DLSFLG.EQ.1).AND.(ADJFLG.EQ.1)) THEN IF (PRTCON.GT.0) THEN WRITE(UNITO,I1000) WRITE(UNITO,I1000) WRITE(UNITO,I1030) 'LIST OF DELETED BOXES FOLLOWS:' END IF ELSE IF ((I.EQ.2).AND.(ADJFLG.NE.1)) THEN IF (PRTCON.GT.0) THEN WRITE(UNITO,I1000) WRITE(UNITO,I1000) WRITE(UNITO,I1030) 1 'The test did not reject the following boxes, but their' WRITE(UNITO,I1030) 1 'diameters were smaller than EPS. NOTE: Diameters are' WRITE(UNITO,I1030) 1 'computed with the scaled norm defined in DIAMCP.' END IF ELSE GOTO 200 END IF C COUNT = 0 CURPTR = LIST(I) 100 CURPTR = PTR(NEXT,CURPTR) IF (CURPTR.EQ.EOLIST) THEN IF (I.EQ.1) NINLST = COUNT GOTO 200 END IF COUNT = COUNT + 1 INDBOX(COUNT) = CURPTR IF (PRTCON.GT.0) THEN WRITE(UNITO,I1000) WRITE(UNITO,I1000) WRITE(UNITO,I1030) 'BOX NUMBER ',COUNT END IF IF (PRTCON.GT.1) THEN WRITE(UNITO,I1030) ' Pointer in the linked list',CURPTR END IF C IF (IDEND.LE.LNDWRK .AND. IIEND.LE.LNIWRK) THEN CALL DCOPY (N,POINTS(1,CURPTR),1,DWORK(ISCX),1) CALL NEWTON (N,DWORK(ISCX),DWORK(IJAC),DWORK(IV),EPSF, 1 EPS,20,NITR,FUNCSC,JACSC,IWORK(IIPVT),IERR) IF (IERR.EQ.0) THEN IF (PRTCON.GT.0) WRITE (UNITO,I1030) 1 'Point Newton method succeeded with NITR =', NITR CALL DCOPY(N,DWORK(ISCX),1,POINTS(1,CURPTR),1) ELSE IF (PRTCON.GT.0) WRITE (UNITO,I1030) 1 'POINT NEWTON METHOD FAILED WITH IERR =', IERR END IF END IF C IF (PRTCON.GT.0) THEN C WRITE(UNITO,I1000) WRITE(UNITO,I1010) 1 'Containing intervals for the coordinates:', 2 ((BOXES(J,K,CURPTR),J=1,2),K=1,N) C WRITE(UNITO,I1020) 'Approximate root:', 1 (POINTS(K,CURPTR),K=1,N) C IF (IJAC-1.LE.LNDWRK) THEN CALL FUNC(N,BOXES(1,1,CURPTR),POINTS(1,CURPTR), 1 DWORK(IFINT), DWORK(ISCLF)) RESID = DNRM2(N,DWORK(ISCLF),1) WRITE(UNITO,I1020) 1 'Euclidean norm of residual at approximate root:',RESID WRITE(UNITO,I1010) 'Interval residual vector:', 1 ((DWORK(IFINT+J+2*K-3),J=1,2),K=1,N) WRITE(UNITO,I1020) 'Residual vector at approx. root:', 1 (DWORK(ISCLF+K-1), K=1,N) END IF END IF C IF (PRTCON.GT.1) THEN WRITE(UNITO,I1030) 'Resolved box number:',BXINFO(5,CURPTR) WRITE(UNITO,I1030) 'Box added at level:',BXINFO(2,CURPTR) IF (BXINFO(1,CURPTR).EQ.10) THEN WRITE(UNITO,I1030) ' Reason: | image | < EPSF' ELSE IF (BXINFO(1,CURPTR).EQ.11) THEN WRITE(UNITO,I1030) ' Reason: test signalled unique root' ELSE IF (BXINFO(1,CURPTR).EQ.18) THEN WRITE(UNITO,I1030) ' Reason: underwent expansion step' ELSE WRITE(UNITO,I1030) ' Reason: ',BXINFO(1,CURPTR) END IF IF ((I.EQ.2).AND.(DLSFLG.EQ.1).AND.(ADJFLG.EQ.1)) THEN WRITE(UNITO,I1030) 'Reason box deleted:',BXINFO(3,CURPTR) WRITE(UNITO,I1030) 'Pointer of intersecting box:', 1 BXINFO(4,CURPTR) END IF END IF IF (PRTCON.GT.2) THEN WRITE(UNITO,I1050) LEFT, RIGHT WRITE(UNITO,I1060) (PINFO2(1,J,CURPTR),PINFO1(J,CURPTR), 1 PINFO2(2,J,CURPTR),J=1,BXINFO(2,CURPTR)) END IF GOTO 100 C 200 CONTINUE C C*********************************************************************** C C Print intersection table for boxes in list 2. C IF (PRTCON.LE.3) RETURN WRITE(UNITO,I1000) WRITE(UNITO,I1000) WRITE(UNITO,I1030) 'Intersection table of boxes in second list' WRITE(UNITO,I1000) BOXPTR = LIST(2) 300 BOXPTR = PTR(NEXT,BOXPTR) IF (BOXPTR.EQ.EOLIST) RETURN WRITE(UNITO,I1030) 'Boxptr = ',BOXPTR CURPTR = LIST(2) 310 CURPTR = PTR(NEXT,CURPTR) IF (CURPTR.EQ.BOXPTR) GOTO 310 IF (CURPTR.EQ.EOLIST) GOTO 300 PTFLG = .TRUE. DO 320 I = 1, N IF (BOXES(2,I,CURPTR).LT.BOXES(1,I,BOXPTR)) GOTO 310 IF (BOXES(1,I,CURPTR).GT.BOXES(2,I,BOXPTR)) GOTO 310 IF (BOXES(2,I,CURPTR).GT.BOXES(1,I,BOXPTR)) PTFLG = .FALSE. IF (BOXES(1,I,CURPTR).LT.BOXES(2,I,BOXPTR)) PTFLG = .FALSE. 320 CONTINUE IF (PTFLG) THEN WRITE(UNITO,I1040) 1 ' intersects BOXPTR ',CURPTR,' with zero volume.' ELSE WRITE(UNITO,I1040) ' Intersects BOXPTR ',CURPTR END IF GOTO 310 C C*********************************************************************** C C Format statements -- C 1000 FORMAT(1X) C 1011 FORMAT(1X,A,25(/2(2X,2(2X,D17.9)))) 1012 FORMAT(1X,A,50(/3(2X,2(2X,D10.2)))) 1013 FORMAT(1X,A,100(/5(2X,2(2X,D30.22)))) C 1020 FORMAT(1X,A,50(/2(2X,D30.22))) 1021 FORMAT(1X,A,25(/4(2X,D17.9))) 1022 FORMAT(1X,A,50(/6(2X,D10.2))) 1023 FORMAT(1X,A,100(/10(2X,D30.22))) C 1030 FORMAT(1X,A,I8) C 1040 FORMAT(1X,A,I8,A) 1050 FORMAT(2X,'Left is ',L1,' and right is ',L1) 1060 FORMAT(20(/' / ',6(I3,L1,I3,' / '))) C END C*********************************************************************** C*********************************************************************** SUBROUTINE ROOTS(N,MAXLST,MAXDP,MAXFT,EPS,EPSF,ADJFLG,DLSFLG, 1 FTEST,FUNC,JAC,CURBOX,CURPT,LINFO,INFO,STACK,STKLVL, 2 DPIVOT,IPIVOT,TSTITR,BOXES,POINTS,PTR,UNUSED,BXINFO, 3 PINFO1,PINFO2,LNDWRK,DWORK,LNIWRK,IWORK,LNLWRK,LWORK) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C Written by: C C R. Baker Kearfott C C and C C Manuel Novoa III C C Department of Mathematics C U.S.L. Box 4-1010 C Lafayette, LA 70504 C C September 29, 1987 C C Part of the generalized bisection package. C C*********************************************************************** C C Called by -- GENBIS C C*********************************************************************** C C Function -- C C This routine controls the overall generalized bisection C scheme. In the comments, steps are labelled as in the paper C paper, R. B. Kearfott, 'Abstract generalized bisection and a cost C bound', Math. Comput. 49, 179 (July, 1987), pp. 187-202. C C*********************************************************************** C C Argument declarations -- C INTEGER N, MAXLST, MAXDP, MAXFT DOUBLE PRECISION EPS, EPSF INTEGER ADJFLG, DLSFLG EXTERNAL FTEST, FUNC, JAC DOUBLE PRECISION CURBOX(2,N), CURPT(N) INTEGER LINFO, INFO(LINFO) DOUBLE PRECISION STACK(2,N,MAXDP) INTEGER STKLVL(MAXDP) LOGICAL DPIVOT(MAXDP) INTEGER IPIVOT(MAXDP), TSTITR(MAXDP) DOUBLE PRECISION BOXES(2,N,MAXLST), POINTS(N,MAXLST) INTEGER PTR(2,MAXLST) LOGICAL UNUSED(MAXLST) INTEGER BXINFO(5,MAXLST) LOGICAL PINFO1(MAXDP,MAXLST) INTEGER PINFO2(2,MAXDP,MAXLST) INTEGER LNDWRK DOUBLE PRECISION DWORK(LNDWRK) INTEGER LNIWRK, IWORK(LNIWRK), LNLWRK LOGICAL LWORK(LNLWRK) C C*********************************************************************** C C Argument descriptions -- (INPUT = set on entry and not alterable) C (OUTPUT = to be set by the routine) C (I/O = set on entry but alterable) C C N is the number of equations and variables. C (INPUT) C C MAXLST MAXLST - 2 is the maximum number of root-containing boxes which C can be stored. C (INPUT) C C MAXDP is the maximum allowable depth of the stack of boxes yet to be C considered. It is also the maximum allowable depth of the C binary search tree. C (INPUT) C C MAXFT is the maximum allowable number of calls to the root-inclusion C test (FTESTH). C (INPUT) C C EPS is the minimum allowable width (in the norm defined by DIAMCP) C of an unresolved box. C (INPUT) C C EPSF is the range tolerance; if it is certain that each component C of the vector function has absolute value less than EPSF C within a box, then that box is signalled as containing a C root. C (INPUT) C C ADJFLG signals whether or not to implement the expansion/deletion C steps (steps 4 of the 1987 Math. Comput. paper (ibid.)): C C If ADJFLG = 1 then expansion/deletion steps are done. C C If ADJFLG is not equal to 1, then expansion/deletion steps C are omitted. C C (INPUT) C C DLSFLG signals whether or not boxes deleted from the list during C the expansion/deletion process are stored in another C list. C C If DLSFLG = 1 then such boxes are stored. C C If DLSFLG is not equal to 1, then such boxes are NOT stored. C C (This flag is only significant if ADJFLG = 1. C C (INPUT) C C CURBOX is the initial region, on entry. It holds the current box C to be considered during the computation. C (I/O) C C CURPT is storage space for an N-vector (of numbers). C (I/O) C C LINFO is the length of the information vector 'INFO'. C (INPUT) C C INFO is a vector containing information about the lists, as well as C other statistics: C C INFO(1) = error type. (ERRFLG : See subroutine ERROUT.) C C INFO(2) = additional error information (ERRVAL). C C INFO(3) = current depth of the stack. C C INFO(4) = the number of boxes which have been examined. C C INFO(5) = the number of boxes in list 1 (which contains C root-containing boxes). C C INFO(6) = the number of boxes in list 2, if it exists. C (See the explanation of DLSFLG above). C Otherwise, it is the number of boxes deleted in C expansion/deletion steps. C C INFO(7) = the header note for list 1. C C INFO(8) = the header node for list 2. C C INFO(9) = the maximum depth reached in the stack. C C INFO(10) = the maximum level reached in the binary tree. C C INFO(11) = the number of calls to the root inclusion test. C C INFO(12) = the number of interval function calls made to C possibly determine that the zero vector is not C in the range. C C INFO(13) = the number of function calls made for the C interval Newton method. C C INFO(14) = the number of interval and scalar Jacobian calls. C C INFO(15) = the number of expanded boxes. C C (OUTPUT) C C STACK is temporary storage space for the stack of boxes yet to be C examined. C (I/O) C C STKLVL is an integer stack whose entries correspond to the entries in C STACK. In it is stored the level (in the binary tree C from bisection) at which the corresponding box in STACK C occurred. C (I/O) C C DPIVOT, IPIVOT, and TSTITR are for additional stack information C associated with boxes stored in STACK. Their entries C correspond to entries in STACK. C (I/O) C C BOXES contains the lists of root-containing boxes. (There are C possibly two such linked lists; see DLSFLG above.) C C (I/O) C C POINTS contains points corresponding to the boxes in BOXES. These C points may simply be midpoints of the boxes, or approximations C to roots. C (I/O) C C PTR contains information about the 'links' between nodes in the C linked lists in BOXES and POINTS. C (I/O) C C UNUSED UNUSED(I) contains the state of the storage areas BOXES(*,*,i) C and POINTS(*,*,I); UNUSED(I) = 'TRUE' means that these C areas are not in use. C (I/O) C C BXINFO, PINFO1, and PINFO2 are storage areas whose entries C correspond to entries in BOXES, POINTS, and UNUSED. C Additional information about the root-containing boxes is C placed in these arrays. C (I/O) C C LNDWRK is the length of the double precision work vector DWORK. C (INPUT) C C DWORK is double precision work space (used in roots, FTEST, etc.) C (I/O) C C LNIWRK is the length of the integer work vector IWORK. C (INPUT) C C IWORK is integer work space (used in ROOTS, FTEST, etc.) C (I/O) C C LNLWRK is the length of the logical work vector LWORK. C (INPUT) C C LWORK is logical work space (used in ROOTS, FTEST, etc.) C (I/O) C C*********************************************************************** C C Internal variable declarations -- C INTEGER ERRFLG, ERRVAL, DEPTH, LEVEL INTEGER NLEAVE, NLIST1, NLIST2, LIST(2), DMAX, LMAX INTEGER NFTCAL, NFCAL1, NFCAL2, NJCALL, NADJ DOUBLE PRECISION DIAM LOGICAL UNKNWN, SIGRT, DP INTEGER IP, RETCON INTEGER TPTR1, TPTR2 INTEGER I, NT C C*********************************************************************** C C Internal variable descriptions -- C C ERRFLG Error type: 0 if no error, another value otherwise C (See subroutine ERROUT for more detailed information.) C C ERRVAL Additional error information C (See subroutine ERROUT for more detailed information.) C C DEPTH The current depth of the stack of boxes yet to be examined. C C LEVEL The current level in the binary search tree C C NLEAVE The number of boxes for which the algorithm has reached a C conclusion about whether they contain roots. C C NLIST1 The number of boxes in list 1 C (List 1 is the list of root-containg boxes.) C C NLIST2 The number of boxes in list 2, if it exists; otherwise the C number of boxes deleted in expansion/deletion steps C (See the description of argument DLSFLG.) C C LIST The pointers to the header nodes of the linked lists C C DMAX The maximum depth reached in the stack of boxes yet to be C considered C C LMAX The maximum level reached in the binary search tree C corresponding to the bisection process C C NFTCAL The number of calls made to FTEST C C NFCAL1 the number of interval function calls made to possibly C determine that the zero vector is not in the range C C NFCAL2 the number of function calls made for the interval Newton C method C C NJCALL the number of interval and scalar Jacobian calls C C NADJ The number of expanded boxes C C DIAM The diameter of the current box C C UNKNWN Set to .TRUE. in FTEST to signal 'unknown' for the C root-inclusion test; set to .FALSE. in FTEST if the C root-inclusion test is true or false. C C SIGRT Set to .TRUE. in FTEST if the root-inclusion test signals a C unique root; set to .FALSE. in FTEST if the root-inclusion C test signals no root C (This flag is meaningful only if UNKNWN = .FALSE.) C C DP is set to .TRUE. if the left box formed from bisection is to C be considered next, and is set to .FALSE. if the right box is C to be considered next. (The left box is that box formed by C replacing the right endpoint of the bisected interval by its C midpoint.) (This choice is made in PVSLCT.) C C IP Set to the index of the coordinate PVSLCT has chosen to be C bisected C C RETCON is the return condition from FTEST. Its values are as C follows. C C 0: The preconditioner matrix could not be computed. C C 1: The root-inclusion test was inconclusive. C C 2: The diameter of the box is less than the domain C tolerance EPS/16. C C 10: Each component of the function has absolute value less C than the range tolerance EPSF over the entire box. C C 11: The box contains a unique root. C C 20: Zero is not in the interval value of the function. C C 21: The image of the box under the interval Newton method C has null intersection with the box. C C TPTR1 and TPTR2 are temporary pointers. C C I A loop index C C NT Temporary storage C C*********************************************************************** C C Common block declarations -- none C C*********************************************************************** C C Fortran-supplied functions and subroutines -- C C MAX C C*********************************************************************** C C Package-supplied functions and subroutines -- C C BISECT, CHKLST, DELLST, DIAMCP, EXPAND, XINFO C C From the linked list and stack management subpackage -- C C ADDBOX, ALLOC, POP, and PUSH C C From LINPACK (BLAS) -- C C DCOPY C C The root inclusion test -- C C (FTEST) (The package-supplied routine is FTESTH; the user C may substitute his or her own routine.) C C*********************************************************************** C C User-supplied functions and subroutines -- none C C*********************************************************************** C C I/O functions -- none C C*********************************************************************** C C Internal constant declarations -- C INTEGER NTTWO, EOLIST DOUBLE PRECISION R, EPSBOX C C*********************************************************************** C C Internal constant descriptions -- C C NTTWO 2*N C C EOLIST The end-of-list value (0) for the linked list structure C C R A parameter used in the expansion/deletion process C C EPSBOX The minimum allowable width (in the norm defined DIAMCP) of C an unresolved box C C*********************************************************************** C DATA EOLIST/0/ C C Beginning of executable statements -- C C Step 0 (housekeeping) C NTTWO = N * 2 R = 0.5D0 IF (ADJFLG.EQ.1) THEN EPSBOX = EPS * (R/2D0)**2 ELSE EPSBOX = EPS END IF C ERRFLG = 0 ERRVAL = 0 LEVEL = 1 DEPTH = 0 C NLEAVE = 0 NLIST1 = 0 NLIST2 = 0 DMAX = 0 LMAX = 0 NFTCAL = 0 NFCAL1 = 0 NFCAL2 = 0 NJCALL = 0 NADJ = 0 C C Check lengths of workspace vectors. C IF (LNDWRK.LT.N) THEN ERRFLG = 2 ERRVAL = N GOTO 900 END IF C C Set all nodes in the linked list pool to 'free'. C DO 10 I = 1, MAXLST UNUSED(I) = .TRUE. 10 CONTINUE C C Create linked lists. (Allocate header nodes.) C DO 40 I = 1, 2 IF ((I.EQ.2).AND.(ADJFLG.EQ.1).AND.(DLSFLG.NE.1)) GOTO 40 CALL ALLOC(TPTR1,MAXLST,UNUSED) IF (TPTR1.EQ.EOLIST) THEN ERRFLG = 20 ERRVAL = MAXLST GOTO 900 END IF LIST(I) = TPTR1 PTR(1,TPTR1) = EOLIST PTR(2,TPTR1) = EOLIST 40 CONTINUE C C*********************************************************************** C C Step 1 (Initialization phase) C DO 110 I = 1, N CURPT(I) = (CURBOX(1,I) + CURBOX(2,I)) / 2D0 110 CONTINUE CALL DIAMCP(N,CURBOX,DIAM) C GOTO 330 C C*********************************************************************** C C Step 2 (Subdivision phase) C 200 CONTINUE C CALL BISECT(N,CURBOX,CURPT,IP,DP,DWORK(1),ERRFLG,ERRVAL) IF (ERRFLG.NE.0) GOTO 900 C IF (LEVEL.LE.MAXDP) THEN IPIVOT(LEVEL) = IP DPIVOT(LEVEL) = DP END IF C LEVEL = LEVEL + 1 CALL PUSH(N,DWORK(1),MAXDP,DEPTH,LEVEL,STACK,STKLVL,ERRFLG,ERRVAL) IF (ERRFLG.NE.0) GOTO 900 C C*********************************************************************** C C Step 3 (Test phase and storage of roots) C 300 CONTINUE C C Step 3(a) C CALL DIAMCP(N,CURBOX,DIAM) C C Step 3(b) C IF ((ADJFLG.EQ.1).AND.(DIAM.LT.(R/2D0)*EPS)) THEN CALL CHKLST(N,MAXLST,LIST,CURBOX,BOXES,PTR,R,EPS, 1 BXINFO,TPTR1) IF (TPTR1.NE.EOLIST) THEN NLIST2 = NLIST2 + 1 IF (DLSFLG.EQ.1) THEN CALL ALLOC(TPTR2,MAXLST,UNUSED) IF (TPTR2.EQ.EOLIST) THEN ERRFLG = 20 ERRVAL = MAXLST GOTO 900 END IF CALL DCOPY(NTTWO,CURBOX,1,BOXES(1,1,TPTR2),1) CALL DCOPY(N,CURPT,1,POINTS(1,TPTR2),1) CALL ADDBOX(LIST(2),TPTR2,PTR) CALL XINFO(TPTR2,19,LEVEL,1,TPTR1,NLEAVE,MAXDP,MAXLST, 1 BXINFO,PINFO1,PINFO2,DPIVOT,IPIVOT,TSTITR) END IF GOTO 500 END IF END IF C C Step 3(c) C 330 CONTINUE CALL FTEST(N,CURBOX,CURPT,DIAM,ERRFLG,ERRVAL,UNKNWN,SIGRT, 1 FUNC,JAC,EPSBOX,EPSF,RETCON,NFCAL1,NFCAL2,NJCALL, 2 MAXDP,DEPTH,LEVEL,STACK,STKLVL, 3 IP,DP,DPIVOT,IPIVOT,TSTITR, 4 LNDWRK,DWORK,LNIWRK,IWORK,LNLWRK,LWORK) NFTCAL = NFTCAL + 1 IF (ERRFLG.NE.0) GOTO 900 IF (NFTCAL.GT.MAXFT) THEN ERRFLG = 1 ERRVAL = MAXFT GOTO 900 END IF C C Step 3(d) C IF (UNKNWN.AND.(DIAM.GE.EPSBOX)) GOTO 200 C C Step 3(e) C IF ((.NOT.UNKNWN).AND.(.NOT.SIGRT)) GOTO 500 C C Step 3(f) C IF ((.NOT.UNKNWN).AND.SIGRT) THEN NLIST1 = NLIST1 + 1 CALL ALLOC(TPTR2,MAXLST,UNUSED) IF (TPTR2.EQ.EOLIST) THEN ERRFLG = 20 ERRVAL = MAXLST GOTO 900 END IF CALL DCOPY(NTTWO,CURBOX,1,BOXES(1,1,TPTR2),1) CALL DCOPY(N,CURPT,1,POINTS(1,TPTR2),1) CALL ADDBOX(LIST(1),TPTR2,PTR) CALL XINFO(TPTR2,RETCON,LEVEL,2,TPTR1,NLEAVE,MAXDP,MAXLST, 1 BXINFO,PINFO1,PINFO2,DPIVOT,IPIVOT,TSTITR) * IF ((ADJFLG.EQ.1).AND.(DIAM.LT.(R/2D0)*EPS)) THEN NT = NLIST1 CALL DELLST(N,MAXLST,NLIST1,TPTR2,LIST,BOXES,PTR,R,EPS, 1 DLSFLG,BXINFO,UNUSED) NLIST2 = NLIST2 + (NT - NLIST1) END IF GOTO 500 END IF C C*********************************************************************** C C Steps 4(a) AND 4(c) C CALL ALLOC(TPTR2,MAXLST,UNUSED) IF (TPTR2.EQ.EOLIST) THEN ERRFLG = 20 ERRVAL = MAXLST GOTO 900 END IF IF (ADJFLG.EQ.1) THEN NADJ = NADJ + 1 NLIST1 = NLIST1 + 1 CALL EXPAND(N,R,CURBOX,BOXES(1,1,TPTR2)) CALL DCOPY(N,CURPT,1,POINTS(1,TPTR2),1) CALL ADDBOX(LIST(1),TPTR2,PTR) ELSE NLIST2 = NLIST2 + 1 CALL DCOPY(NTTWO,CURBOX,1,BOXES(1,1,TPTR2),1) CALL DCOPY(N,CURPT,1,POINTS(1,TPTR2),1) CALL ADDBOX(LIST(2),TPTR2,PTR) END IF CALL XINFO(TPTR2,18,LEVEL,0,0,NLEAVE,MAXDP,MAXLST,BXINFO, 1 PINFO1,PINFO2,DPIVOT,IPIVOT,TSTITR) C C Step 4(b) C IF ((ADJFLG.EQ.1).AND.(NLIST1.GT.1)) THEN NT = NLIST1 CALL DELLST(N,MAXLST,NLIST1,TPTR2,LIST,BOXES,PTR,R,EPS, 1 DLSFLG,BXINFO,UNUSED) NLIST2 = NLIST2 + (NT - NLIST1) END IF C C*********************************************************************** C C Step 5 (Backtrack) C 500 CONTINUE LMAX = MAX(LMAX,LEVEL) DMAX = MAX(DMAX,DEPTH) C C Steps 5(a) - 5(d) C NLEAVE = NLEAVE + 1 IF (DEPTH.NE.0) THEN CALL POP(N,CURBOX,MAXDP,DEPTH,LEVEL,STACK,STKLVL,ERRFLG,ERRVAL) IF (ERRFLG.NE.0) GOTO 900 DO 520 I = 1, N CURPT(I) = (CURBOX(1,I) + CURBOX(2,I)) / 2D0 520 CONTINUE IF ((LEVEL-1.GE.1).AND.(LEVEL-1.LE.MAXDP)) THEN DPIVOT(LEVEL-1) = .NOT.DPIVOT(LEVEL-1) END IF GOTO 300 END IF C C*********************************************************************** C C Set information vector for return 900 CONTINUE C INFO(1) = ERRFLG INFO(2) = ERRVAL INFO(3) = DEPTH INFO(4) = NLEAVE INFO(5) = NLIST1 INFO(6) = NLIST2 INFO(7) = LIST(1) INFO(8) = LIST(2) INFO(9) = DMAX INFO(10) = LMAX INFO(11) = NFTCAL INFO(12) = NFCAL1 INFO(13) = NFCAL2 INFO(14) = NJCALL INFO(15) = NADJ C RETURN END C*********************************************************************** C*********************************************************************** SUBROUTINE BISECT(N,BOX1,CURPT,IP,DP,BOX2,ERRFLG,ERRVAL) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C Written by: C C R. Baker Kearfott C C and C C Manuel Novoa III C C Department of Mathematics C U.S.L. Box 4-1010 C Lafayette, LA 70504 C C September 29, 1987 C C Part of the generalized bisection package C C*********************************************************************** C C Called by -- ROOTS C C*********************************************************************** C C Function -- C C This routine bisects the box in BOX1 by bisecting the IP-th C coordinate interval. It then stores the half dictated by DP C back in BOX1 and the other half in BOX2. It also stores the midpoint C of the new BOX1 in CURPT. C C*********************************************************************** C C Argument declarations -- C INTEGER N DOUBLE PRECISION BOX1(2,N), CURPT(N) INTEGER IP LOGICAL DP DOUBLE PRECISION BOX2(2,N) INTEGER ERRFLG, ERRVAL C C*********************************************************************** C C Argument descriptions -- (INPUT = set on entry and not alterable) C (OUTPUT = to be set by the routine) C (I/O = set on entry but alterable) C C N is the number of variables and equations. C (INPUT) C C BOX1 On entry, BOX1 contains the box to be bisected. C On return, BOX1 contains the half of the bisected box C which DP defines to be the next current box. C (I/O) C C CURPT On return, CURPT is set to the midpoint of the box in BOX1. C (OUTPUT) C C IP is the index of the coordinate interval to be bisected. C (INPUT) C C DP indicates the half of the bisected box (left or right) to be C the next current box. (The left box is that box formed by C replacing the right endpoint of the bisected interval by its C midpoint.) (This choice is made in PVSLCT.) C (INPUT) C C BOX2 On return, BOX2 contains the half of the bisected box to be C placed on the stack for future examination. C (OUTPUT) C C ERRFLG On return, ERRFLG is set to 0 if no error has occured, and is C otherwise set to some other value. (See subroutine ERROUT C for more detailed information.) C (OUTPUT) C C ERRVAL contains additional error information. (See subroutine ERROUT C for more detailed information.) C (OUTPUT) C C*********************************************************************** C C Internal variable declarations -- C DOUBLE PRECISION TMIDPT INTEGER I C C*********************************************************************** C C Internal variable descriptions -- C C TMIDPT is a temporary midpoint value. C C I is a loop index. C C*********************************************************************** C C Common block declarations -- none C C*********************************************************************** C C Fortran-supplied functions and subroutines -- none C C*********************************************************************** C C Package-supplied functions and subroutines -- none C C DCOPY (from LINPACK) C C*********************************************************************** C C User-supplied functions and subroutines -- none C C*********************************************************************** C C I/O functions -- none C C*********************************************************************** C C Internal constant declarations -- C INTEGER NTTWO LOGICAL LEFT, RIGHT C C*********************************************************************** C C Internal constant descriptions -- C C NTTWO 2*N C C LEFT the logical value associated in this program with the C left box formed from bisection. (See DP above.) C C RIGHT the logical value associated in this program with the C right box formed from bisection. (See DP above.) C C*********************************************************************** C DATA LEFT/.TRUE./, RIGHT/.FALSE./ C C Beginning of executable statements -- C NTTWO = 2 * N C C Check to see if coordinate index selected for bisection is within C bounds. C IF ((IP.LE.0).OR.(IP.GT.N)) THEN ERRFLG = 32 ERRVAL = IP RETURN END IF C C Compute the midpoint of the interval to be bisected. C TMIDPT = (BOX1(1,IP) + BOX1(2,IP)) / 2D0 C C Make sure the interval to be bisected does not have endpoints C which are adjacent machine numbers. C IF ((TMIDPT.LE.BOX1(1,IP)).OR.(TMIDPT.GE.BOX1(2,IP))) THEN ERRFLG = 33 ERRVAL = IP RETURN END IF C C Perform the bisection. C CALL DCOPY(NTTWO,BOX1,1,BOX2,1) IF (DP.EQV.LEFT) THEN C WANT LEFT BOX TO BE NEW CURRENT BOX BOX1(2,IP) = TMIDPT BOX2(1,IP) = TMIDPT ELSE C WANT RIGHT BOX TO BE NEW CURRENT BOX BOX1(1,IP) = TMIDPT BOX2(2,IP) = TMIDPT END IF C C Compute the midpoint of the new current box. C DO 10 I = 1, N CURPT(I) = (BOX1(1,I) + BOX1(2,I)) / 2D0 10 CONTINUE C RETURN END C*********************************************************************** C*********************************************************************** SUBROUTINE CHKLST(N,MAXLST,LIST,CURBOX,BOXES,PTR,R,EPS,BXINFO, 1 CURPTR) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C Written by: C C R. Baker Kearfott C C and C C Manuel Novoa III C C Department of Mathematics C U.S.L. Box 4-1010 C Lafayette, LA 70504 C C September 29, 1987 C C Part of the generalized bisection package C C*********************************************************************** C C Called by -- ROOTS C C*********************************************************************** C C Function -- C C The general purpose of this subroutine is to determine when to C remove a redundantly listed box from the list of root containing C boxes. Such redundancy could occur when a root happens C to occur near a boundary of two or more boxes, when the C Jacobian matrix is nearly singular at the root, or when the C interval approximation to the Jacobian matrix is poor. The C removal process will work well in the first case, but functions C heuristically in the other two cases. C C Specifically this subroutine returns a pointer to the first C box which intersects the box in CURBOX and which has a sufficiently C small diameter, if such a box exists. If there is no such C intersecting box, this routine returns the end-of-list pointer. C C*********************************************************************** C C Argument declarations -- C INTEGER N, MAXLST, LIST(2) DOUBLE PRECISION CURBOX(2,N), BOXES(2,N,MAXLST) INTEGER PTR(2,MAXLST) DOUBLE PRECISION R, EPS INTEGER BXINFO(5,MAXLST), CURPTR C C*********************************************************************** C C Argument descriptions -- (INPUT = set on entry and not alterable) C (OUTPUT = to be set by the routine) C (I/O = set on entry but alterable) C C N is the number of equations and variables. C (INPUT) C C MAXLST MAXLST - 2 is the maximum number of root-containing boxes which C can be stored. C (INPUT) C C LIST contains the pointers to the header nodes of the linked lists. C (INPUT) C C CURBOX contains the current box, to be compared with boxes in the C list. C (INPUT) C C BOXES is an array which holds the coordinates of the boxes in C the lists. C (INPUT) C C PTR contains information about the locations of boxes in the C array BOXES. C (INPUT) C C R is a parameter used in the expansion/deletion process C and here to determine when diameters are sufficiently C small. C (INPUT) C C EPS is the minimum allowable width (in the norm defined by C DIAMCP) of an unresolved box. C (INPUT) C C BXINFO is a storage areas whose entries correspond to entries C in BOXES. Additional information about the root-containing C boxes is placed here. C (INPUT) C C CURPTR On exit, CURPTR is set to 0 (the end-of-list pointer) if no C box in the list intersected the box in CURBOX; otherwise, C CURPTR is set to the pointer of the first box in the list C which intersected CURBOX. C (OUTPUT) C C*********************************************************************** C C Internal variable declarations -- C DOUBLE PRECISION DIAM1, DIAM2 INTEGER I C C*********************************************************************** C C Internal variable descriptions -- C C DIAM1 is the diameter of the current box. C C DIAM2 is a temporary variable for the diameters of boxes in the C list. C C I is a loop index. C C*********************************************************************** C C Common block declarations -- C DOUBLE PRECISION TOL1, TOL2, TOL3 COMMON /CONFG2/ TOL1, TOL2, TOL3 C C /CONFG2/ passes tolerances from INPUT to this routine and to C the routines HNSNGT and DELLST. These tolerances C are described in INPUT with the configuration file. C C*********************************************************************** C C Fortran-supplied functions and subroutines -- none C C*********************************************************************** C C Package-supplied functions and subroutines -- C C DIAMCP C C*********************************************************************** C C User-supplied functions and subroutines -- none C C*********************************************************************** C C I/O functions -- none C C*********************************************************************** C C Internal constant declarations -- C DOUBLE PRECISION TOL, REPSD2 INTEGER EOLIST, PREV, NEXT C C*********************************************************************** C C Internal constant descriptions -- C C TOL the intersection tolerance C C REPSD2 R*EPS/2 C C EOLIST the end-of-list value (0) for the linked list structure C C PREV PTR(PREV,*) is the pointer to the box preceeding the C box with pointer '*' in the linked list. C C NEXT PTR(NEXT,*) is the pointer to the box following the C box with pointer '*' in the linked list. C C*********************************************************************** C DATA EOLIST/0/, PREV/1/, NEXT/2/ C C Beginning of executable statements -- C REPSD2 = R * EPS / 2D0 C TOL = 10D0 * EPS TOL = TOL3 CALL DIAMCP(N,CURBOX,DIAM1) CURPTR = LIST(1) 10 CURPTR = PTR(NEXT,CURPTR) IF (CURPTR.EQ.EOLIST) GOTO 20 C C Remove comment in next line if only checking expanded boxes. C C IF (BXINFO(1,CURPTR).NE.18) GOTO 10 CALL DIAMCP(N,BOXES(1,1,CURPTR),DIAM2) IF (DIAM1+DIAM2.GT.REPSD2) GOTO 10 DO 15 I = 1, N IF (CURBOX(2,I).LT.BOXES(1,I,CURPTR)-TOL) GOTO 10 IF (CURBOX(1,I).GT.BOXES(2,I,CURPTR)+TOL) GOTO 10 15 CONTINUE 20 RETURN END C*********************************************************************** C*********************************************************************** SUBROUTINE DELLST(N,MAXLST,NINLST,BOXPTR,LIST,BOXES,PTR,R,EPS, 1 DLSFLG,BXINFO,UNUSED) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C Written by: C C R. Baker Kearfott C C and C C Manuel Novoa III C C Department of Mathematics C U.S.L. Box 4-1010 C Lafayette, LA 70504 C C September 29, 1987 C C Part of the generalized bisection package C C*********************************************************************** C C Called by -- ROOTS C C*********************************************************************** C C Function -- C C This subroutine removes those boxes from the list of root-containing C boxes which intersect a specified box and which have sufficiently C small diameters. C C The general purpose is to reduce redundancy in the list of C boxes. Such redundancy could occur when a root happens C to occur near a boundary of two or more boxes, when the C Jacobian matrix is nearly singular at the root, or when the C interval approximation to the Jacobian matrix is poor. The C removal process will work well in the first case, but functions C heuristically in the other two cases. C C*********************************************************************** C C Argument declarations -- C INTEGER N, MAXLST, NINLST, BOXPTR, LIST(2) DOUBLE PRECISION BOXES(2,N,MAXLST) INTEGER PTR(2,MAXLST) DOUBLE PRECISION R, EPS INTEGER DLSFLG, BXINFO(5,MAXLST) LOGICAL UNUSED(MAXLST) C C*********************************************************************** C C Argument descriptions -- (INPUT = set on entry and not alterable) C (OUTPUT = to be set by the routine) C (I/O = set on entry but alterable) C C N is the number of equations and variables. C (INPUT) C C MAXLST MAXLST - 2 is the maximum number of root-containing boxes which C can be stored. C (INPUT) C C NINLST is the number of root-containing boxes presently in the list. C (INPUT) C C BOXPTR is the pointer to the box whose intersection with other C boxes in the list is to be checked. C (INPUT) C C LIST contains the pointers to the header nodes of the linked lists. C (INPUT) C C BOXES is an array which holds the coordinates of the boxes in C the lists. C (INPUT) C C PTR contains information about the locations of boxes in the C array BOXES. C (INPUT) C C R is a parameter used in the expansion/deletion process C and here to determine when diameters are sufficiently C small. C (INPUT) C C EPS is the minimum allowable width (in the norm defined by C DIAMCP) of an unresolved box. C (INPUT) C C DLSFLG is equal to 1 if boxes removed from the list are to be C stored in another list. If DLSFLG is not equal to 1, then C such boxes are simply discarded. C (INPUT) C C BXINFO is a storage areas whose entries correspond to entries C in BOXES. Additional information about the root-containing C boxes is placed here. C (INPUT) C C UNUSED UNUSED(I) contains the state of the storage areas BOXES(*,*,i) C and BXINFO(*,I); UNUSED(I) = 'TRUE' means that these C areas are not in use. C (I/O) C C*********************************************************************** C C Internal variable declarations -- C DOUBLE PRECISION DIAM1, DIAM2 INTEGER CURPTR, I, TPTR C C*********************************************************************** C C Internal variable descriptions -- C C DIAM1 is the diameter of the box whose intersection with other C boxes is to be checked. C C DIAM2 is a temporary variable for diameters of boxes in the list. C C CURPTR is a temporary pointer. C C I is a loop index. C C TPTR is a temporary pointer. C C*********************************************************************** C C Common block declarations -- C DOUBLE PRECISION TOL1, TOL2, TOL3 COMMON /CONFG2/ TOL1, TOL2, TOL3 C C /CONFG2/ passes tolerances from INPUT to this routine and to C the routines CHKLST and HNSNG. These tolerances C are described in INPUT with the configuration file. C C*********************************************************************** C C Fortran-supplied functions and subroutines -- none C C*********************************************************************** C C Package-supplied functions and subroutines -- C C ADDBOX, DELBOX, FREE, DIAMCP C C*********************************************************************** C C user-supplied functions and subroutines -- none C C*********************************************************************** C C I/O functions -- none C C*********************************************************************** C C Internal constant declarations -- C DOUBLE PRECISION TOL, REPSD2 INTEGER EOLIST, PREV, NEXT C C*********************************************************************** C C Internal constant descriptions -- C C TOL the intersection tolerance C C REPSD2 R*EPS/2 C C EOLIST the end-of-list value (0) for the linked list structure C C PREV PTR(PREV,*) is the pointer to the box preceeding the C box with pointer '*' in the linked list. C C NEXT PTR(NEXT,*) is the pointer to the box following the C box with pointer '*' in the linked list. C C*********************************************************************** C DATA EOLIST/0/, PREV/1/, NEXT/2/ C C Beginning of executable statements -- C REPSD2 = R * EPS / 2D0 C TOL = 10D0 * EPS TOL = TOL3 CALL DIAMCP(N,BOXES(1,1,BOXPTR),DIAM1) IF (DIAM1.GE.REPSD2) THEN CURPTR = EOLIST GOTO 20 END IF CURPTR = LIST(1) 10 CURPTR = PTR(NEXT,CURPTR) 12 IF (CURPTR.EQ.BOXPTR) GOTO 10 IF (CURPTR.EQ.EOLIST) GOTO 20 CALL DIAMCP(N,BOXES(1,1,CURPTR),DIAM2) IF (DIAM1+DIAM2.GT.REPSD2) GOTO 10 DO 15 I = 1, N IF (BOXES(2,I,BOXPTR).LT.BOXES(1,I,CURPTR)-TOL) GOTO 10 IF (BOXES(1,I,BOXPTR).GT.BOXES(2,I,CURPTR)+TOL) GOTO 10 15 CONTINUE NINLST = NINLST - 1 TPTR = PTR(NEXT,CURPTR) CALL DELBOX(CURPTR,PTR) IF (DLSFLG.EQ.1) THEN CALL ADDBOX(LIST(2),CURPTR,PTR) BXINFO(3,CURPTR) = 3 BXINFO(4,CURPTR) = BOXPTR ELSE CALL FREE(CURPTR,MAXLST,UNUSED) END IF CURPTR = TPTR GOTO 12 20 RETURN C END C*********************************************************************** C*********************************************************************** SUBROUTINE DIAMCP(N,X,DIAM) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C Written by: C C R. Baker Kearfott C C and C C Manuel Novoa III C C Department of Mathematics C U.S.L. Box 4-1010 C Lafayette, LA 70504 C C September 29, 1987 C C Part of the generalized bisection package C C*********************************************************************** C C CALLED BY -- CHKLST, DELLST, FTESTH, and ROOTS. C C*********************************************************************** C C Function -- C C This routine computes the diameter of the N-box X (in a scaled C infinity norm) and returns the value in DIAM. C C*********************************************************************** C C Argument declarations -- C INTEGER N DOUBLE PRECISION X(2,N) DOUBLE PRECISION DIAM C C*********************************************************************** C C Argument descriptions -- (INPUT = set on entry and not alterable) C (OUTPUT = to be set by the routine) C (I/O = set on entry but alterable) C C N is the number of variables and equations. C (INPUT) C C X X(1,I) is the left endpoint of the I-th coordinate interval C of the box, and X(2,I) is the right endpoint of the I-th C coordinate interval of the box, for I between 1 and N. C (INPUT) C C DIAM is set to the diameter of X in a scaled infinity norm. C (OUTPUT) C C*********************************************************************** C C Internal variable declarations -- C INTEGER I DOUBLE PRECISION WIDTH, BIG C C*********************************************************************** C C Common block declarations -- none C C*********************************************************************** C C Fortran-supplied functions and subroutines -- C C MAX, ABS C C*********************************************************************** C C Package-supplied functions and subroutines -- none C C*********************************************************************** C C User-supplied functions and subroutines -- none C C*********************************************************************** C C I/O functions -- none C C*********************************************************************** C C Internal constant declarations -- none C C*********************************************************************** C C Beginning of executable statements -- C DIAM = 0D0 DO 10 I = 1, N WIDTH = X(2,I) - X(1,I) BIG = MAX(1D0,ABS(X(1,I)),ABS(X(2,I))) DIAM = MAX (DIAM,WIDTH/BIG) 10 CONTINUE C RETURN END C*********************************************************************** C*********************************************************************** SUBROUTINE DVSBIN(X,D,N,A,ICASE,R1,R2) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C C Written by: C C R. Baker Kearfott C C and C C Manuel Novoa III C C Department of Mathematics C U.S.L. Box 4-1010 C Lafayette, LA 70504 C C September 29, 1987 C C Part of the generalized bisection package C C*********************************************************************** C C Called by -- HNSNG C C*********************************************************************** C C Function -- C C This routine performs the following portions of an interval C Gauss-Seidel step: C C (1) the extended-value interval division; C C (2) subtraction of the quotient (1) from a point value; C C and C C (3) intersection of the difference (2) with the appropriate C component of the original box. C C The number of resulting intervals is returned in ICASE. These C intervals are stored as follows: C C If ICASE = 0 then no intervals are stored; C C if ICASE = 1 then the interval is stored in R1; C C and C C if ICASE = 2 then one interval is stored in R1 and the other C interval is stored in R2. C C*********************************************************************** C C Argument declarations -- C DOUBLE PRECISION X, D(2), N(2), A(2) INTEGER ICASE DOUBLE PRECISION R1(2), R2(2) C C*********************************************************************** C C Argument descriptions -- (INPUT = set on entry and not alterable) C (OUTPUT = to be set by the routine) C (I/O = set on entry but alterable) C C X is the point from which the interval is to be subtracted. C (INPUT) C C D is the denominator of the interval quotient. C (INPUT) C C N is the numerator of the interval quotient. C (INPUT) C C A is the appropriate component interval of the box to which C the interval Gauss-Seidel is to be applied. C (INPUT) C C ICASE indicates the number of intervals returned. C (OUTPUT) C C R1 will contain the first resultant interval. C (OUTPUT) C C R2 will contain the second resultant interval. C (OUTPUT) C C*********************************************************************** C C Internal variable declarations -- C INTEGER XCASE DOUBLE PRECISION B(2) C C*********************************************************************** C C Internal variable descriptions -- C C XCASE indicates the type of the current temporary C extended-value result of any of the extended-value C interval arithmetic functions (XDIV, XINT, XSCLSB). C See the documentation of these routines for more C information. C C B is a temporary interval. C C*********************************************************************** C C Common block declarations -- none C C*********************************************************************** C C Fortran-supplied functions and subroutines -- none C C*********************************************************************** C C Package-supplied functions and subroutines -- C C XDIV, XINT, XSCLSB C C*********************************************************************** C C User-supplied functions and subroutines -- none C C*********************************************************************** C C I/O functions -- none C C*********************************************************************** C C Internal constant declarations -- none C C*********************************************************************** C C Beginning of executable statements -- C C Save the interval A in case it is overwritten before needed. C B(1) = A(1) B(2) = A(2) C C Do the extended-value interval division. C CALL XDIV(XCASE,N,D,R1,R2) C C If the result is a single interval, do the subtraction and C the intersection, then return. C IF (XCASE.NE.5) THEN CALL XSCLSB(XCASE,X,R1,R1) CALL XINT(XCASE,B,R1,R1) ICASE = XCASE RETURN END IF C C Since the result of the division is two intervals, do the C subtraction and intersection with the first one and place C the result in R1. C XCASE = 3 CALL XSCLSB(XCASE,X,R1,R1) CALL XINT(XCASE,B,R1,R1) ICASE = XCASE XCASE = 2 C C Now do the subtraction and intersection with the second interval. C If the corresponding result with the first interval was empty, C place the result using the second interval in R1. Otherwise, C place the result in R2 (unless that result is empty). C IF (ICASE.EQ.0) THEN CALL XSCLSB(XCASE,X,R2,R1) CALL XINT(XCASE,B,R1,R1) ICASE = XCASE RETURN ELSE CALL XSCLSB(XCASE,X,R2,R2) CALL XINT(XCASE,B,R2,R2) ICASE = ICASE + XCASE RETURN END IF C END C*********************************************************************** C*********************************************************************** SUBROUTINE ERROUT(N,MAXDP,MAXLST,MAXFT,UNITO,CURBOX,CURPT,LINFO, 1 INFO,STACK,STKLVL,DPIVOT,IPIVOT,TSTITR) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C Written by: C C R. Baker Kearfott C C and C C Manuel Novoa III C C Department of Mathematics C U.S.L. Box 4-1010 C Lafayette, LA 70504 C C September 29, 1987 C C Part of the generalized bisection package C C*********************************************************************** C C Called by -- GENBIS C C*********************************************************************** C C Function -- C C This routine prints all error messages. It will also, depending C on the error, output the boxes which had not yet been resolved at C the time the error occurred. This is done so that the user does not C need to run the program again for the entire region if it has C terminated for reasons such as stack overflow, no free space in the C linked list, maximu, number of FTEST calls exceeded, etc. C C*********************************************************************** C C Argument declarations -- C INTEGER N, MAXDP, MAXLST, MAXFT INTEGER UNITO DOUBLE PRECISION CURBOX(2,N), CURPT(N) INTEGER LINFO, INFO(LINFO) DOUBLE PRECISION STACK(2,N,MAXDP) INTEGER STKLVL(MAXDP) LOGICAL DPIVOT(MAXDP) INTEGER IPIVOT(MAXDP), TSTITR(MAXDP) C C*********************************************************************** C C Argument descriptions -- (INPUT = set on entry and not alterable) C (OUTPUT = to be set by the routine) C (I/O = set on entry but alterable) C C N is the number of equations and variables. C (INPUT) C C MAXDP is the maximum allowable depth of the stack of boxes yet to be C considered. It is also the maximum allowable depth of the C binary search tree. C (INPUT) C C MAXLST MAXLST - 2 is the maximum number of root-containing boxes which C can be stored. C (INPUT) C C MAXFT is the maximum allowable number of calls to the root-inclusion C test (FTESTH). C (INPUT) C C UNITO is the unit number of the file (device) to which all output C will be directed. C (INPUT) C C CURBOX Space for an N-dimensional interval vector (box) used C to store the box being examined at any given time. C (INPUT) C C CURPT Space for a point in N dimensions, which usually contains C either the midpoint of CURBOX or some approximation to a C root in CURBOX. C (INPUT) C C LINFO is the number of spaces in the information vector INFO. C (INPUT) C C INFO is a vector containing information about the lists, as well as C other statistics: C C INFO(1) = error type. (ERRFLG : See subroutine ERROUT.) C C INFO(2) = additional error information (ERRVAL). C C INFO(3) = current depth of the stack. C C INFO(4) = the number of boxes which have been examined. C C INFO(5) = the number of boxes in list 1 (which contains C root-containing boxes). C C INFO(6) = the number of boxes in list 2, if it exists. C (See the explanation of DLSFLG above). C Otherwise, it is the number of boxes deleted in C expansion/deletion steps. C C INFO(7) = the header note for list 1. C C INFO(8) = the header node for list 2. C C INFO(9) = the maximum depth reached in the stack. C C INFO(10) = the maximum level reached in the binary tree. C C INFO(11) = the number of calls to the root inclusion test. C C INFO(12) = the number of interval function calls made to C possibly determine that the zero vector is not C in the range. C C INFO(13) = the number of function calls made for the C interval Newton method. C C INFO(14) = the number of interval and scalar Jacobian calls. C C INFO(15) = the number of expanded boxes. C C (INPUT) C C STACK is temporary storage space for the stack of boxes yet to be C examined. C (INPUT) C C STKLVL is an integer stack whose entries correspond to the entries in C STACK. In it is stored the level (in the binary tree C from bisection) at which the corresponding box in STACK C occurred. C (INPUT) C C DPIVOT, IPIVOT, and TSTITR are for additional stack information C associated with boxes stored in STACK. Their entries C correspond to entries in STACK. C (INPUT) C C*********************************************************************** C C Internal variable declarations -- C INTEGER ERRFLG, ERRVAL, DEPTH INTEGER NLEAVE, NLIST1, NLIST2, LIST(2), DMAX, LMAX INTEGER NFTCAL, NFCAL1, NFCAL2, NJCALL, NADJ INTEGER J, K, LEVEL C C*********************************************************************** C C Internal variable descriptions -- C C ERRFLG Error type: 0 if no error, another value otherwise C (See error messages below for details.) C C ERRVAL Additional error information C (See error messages below for details.) C C DEPTH The current depth of the stack of boxes which still need to C be considered. C C NLEAVE The number of boxes for which the algorithm has reached a C conclusion about whether they contain roots. C C NLIST1 The number of boxes in list 1 C (List 1 is the list of root-containg boxes.) C C NLIST2 The number of boxes in list 2, if it exists; otherwise the C number of boxes deleted in expansion/deletion steps C (See the description of argument DLSFLG.) C C LIST The pointers to the header nodes of the linked lists C C DMAX The maximum depth reached in the stack of boxes yet to be C considered C C LMAX The maximum level reached in the binary search tree C corresponding to the bisection process C C NFTCAL The number of calls made to FTEST C C NFCAL1 the number of interval function calls made to possibly C determine that the zero vector is not in the range C C NFCAL2 the number of function calls made for the interval Newton C method C C NJCALL The number of interval and scalar Jacobian calls C C NADJ The number of expanded boxes C C J and K are loop indices. C C LEVEL a temporary value used to store the level of a box in the C binary search tree C C*********************************************************************** C C Common block declarations -- none C C*********************************************************************** C C Fortran-supplied functions and subroutines -- none C C*********************************************************************** C C Package-supplied functions and subroutines -- C C POP C C*********************************************************************** C C User-supplied functions and subroutines -- none C C*********************************************************************** C C I/O functions -- C C Output: This routine outputs a message corresponding to the value C if ERRFLG to UNITO. If an error occurs, it may also output C additional information about the program's state. Such C information may include the box being examined at the C time the error occurred, as well as the boxes remaining C on the stack. C C*********************************************************************** C C Internal constant declarations -- none C C*********************************************************************** C C Beginning of executable statements -- C ERRFLG = INFO(1) ERRVAL = INFO(2) DEPTH = INFO(3) NLEAVE = INFO(4) NLIST1 = INFO(5) NLIST2 = INFO(6) LIST(1) = INFO(7) LIST(2) = INFO(8) DMAX = INFO(9) LMAX = INFO(10) NFTCAL = INFO(11) NFCAL1 = INFO(12) NFCAL2 = INFO(13) NJCALL = INFO(14) NADJ = INFO(15) C C*********************************************************************** C C Error messages C IF (ERRFLG.EQ.0) THEN WRITE(UNITO,1030) 'SUCCESSFUL COMPLETION OF BINARY SEARCH.' GOTO 900 C ELSE IF (ERRFLG.EQ.1) THEN WRITE(UNITO,1030) 'THE MAXIMUM NUMBER OF CALLS TO FTEST',ERRVAL WRITE(UNITO,1030) 1 'HAS BEEN EXCEEDED. CHANGE THE VALUE OF MAXFT IN YOUR' WRITE(UNITO,1030) 1 'INPUT FILE AND RERUN, RUN THE PROGRAM FOR EACH OF THE' WRITE(UNITO,1030) 1 'BOXES NOT YET EXAMINED, INCREASE THE DOMAIN TOLERANCE EPS' WRITE(UNITO,1030) 1 'OR THE RANGE TOLERANCE EPSF, OR TRY TO DETERMINE PROPERTIES' WRITE(UNITO,1030) 1 'OF THE SYSTEM WHICH MAKE IT DIFFICULT TO SOLVE.' C ELSE IF (ERRFLG.EQ.2) THEN WRITE(UNITO,1030) 'ROOTS REQUIRES A MINIMUM OF ',ERRVAL WRITE(UNITO,1030) 'FOR THE DIMENSION OF DWORK. CHANGE' WRITE(UNITO,1030) 'THE PARAMETER MDWORK IN GENBIS TO THIS' WRITE(UNITO,1030) 'VALUE AND RECOMPILE.' GOTO 900 C ELSE IF (ERRFLG.EQ.3) THEN WRITE(UNITO,1030) 'THE DIMENSION OF THE SYSTEM ',ERRVAL WRITE(UNITO,1030) 'IS GREATER THAN THE MAXIMUM ALLOWED' WRITE(UNITO,1030) 'FOR IN THE PROGRAM. CHANGE PARAMETER MN' WRITE(UNITO,1030) 'IN GENBIS AND RECOMPILE.' GOTO 900 C ELSE IF (ERRFLG.EQ.4) THEN WRITE(UNITO,1030) 'THE TEST ROUTINE REQUIRES ',ERRVAL WRITE(UNITO,1030) 'AS THE DIMENSION OF THE DOUBLE PRECISION' WRITE(UNITO,1030) 'WORK ARRAY. CHANGE PARAMETER MDWORK IN' WRITE(UNITO,1030) 'GENBIS TO THIS VALUE AND RECOMPILE.' GOTO 900 C ELSE IF (ERRFLG.EQ.5) THEN WRITE(UNITO,1030) 'THE TEST ROUTINE REQUIRES ',ERRVAL WRITE(UNITO,1030) 'AS THE DIMENSION OF THE INTEGER' WRITE(UNITO,1030) 'WORK ARRAY. CHANGE PARAMETER MIWORK IN' WRITE(UNITO,1030) 'GENBIS TO THIS VALUE AND RECOMPILE.' GOTO 900 C ELSE IF (ERRFLG.EQ.6) THEN WRITE(UNITO,1030) 'THE TEST ROUTINE REQUIRES ',ERRVAL WRITE(UNITO,1030) 'AS THE DIMENSION OF THE LOGICAL' WRITE(UNITO,1030) 'WORK ARRAY. CHANGE PARAMETER MLWORK IN' WRITE(UNITO,1030) 'GENBIS TO THIS VALUE AND RECOMPILE.' GOTO 900 C ELSE IF (ERRFLG.EQ.10) THEN WRITE(UNITO,1030) 'THE MAXIMUM DEPTH OF THE STACK,',ERRVAL WRITE(UNITO,1030) 1 'HAS BEEN REACHED. CHANGE PARAMETER MMAXDP IN GENBIS' WRITE(UNITO,1030) 1 'AND RECOMPILE, RUN THE PROGRAM FOR EACH OF THE BOXES' WRITE(UNITO,1030) 1 'NOT YET EXAMINED, OR TRY TO DETERMINE PROPERTIES' WRITE(UNITO,1030) 1 'SYSTEM WHICH MAKE IT DIFFICULT TO SOLVE.' C ELSE IF (ERRFLG.EQ.11) THEN WRITE(UNITO,1030) 'AN ATTEMPT HAS BEEN MADE TO POP SOMETHING' WRITE(UNITO,1030) 'FROM THE STACK WHEN IT WAS EMPTY.' GOTO 900 C ELSE IF (ERRFLG.EQ.12) THEN WRITE(UNITO,1030) 'AN ILLEGAL VALUE ',ERRVAL WRITE(UNITO,1030) 'FOR THE STACK DEPTH HAS OCCURRED IN PUSH.' WRITE(UNITO,1030) 'THIS IS A FATAL ERROR.' GOTO 900 C ELSE IF (ERRFLG.EQ.13) THEN WRITE(UNITO,1030) 'AN ILLEGAL VALUE ',ERRVAL WRITE(UNITO,1030) 'FOR THE STACK DEPTH HAS OCCURRED IN POP.' WRITE(UNITO,1030) 'THIS IS A FATAL ERROR.' GOTO 900 C ELSE IF (ERRFLG.EQ.14) THEN WRITE(UNITO,1030) 'THE DIMENSION N =', N WRITE(UNITO,1030) 1 'IS GREATER THAN PARAMETER MN2 IN INPUT, POLFUN, POLJAC,' WRITE(UNITO,1030) 1 'POLFPT, AND SCLFPT.' WRITE(UNITO,1030) 1 'THIS CAUSES A DIMENSION PROBLEM IN COMMON BLOCKS EQUAT' WRITE(UNITO,1030) 1 'AND COEFS IN THESE ROUTINES. MAKE MN2 EQUAL IN ALL OF' WRITE(UNITO,1030) 1 'THESE ROUTINES AND GREATER THAN OR EQUAL TO N.' GOTO 900 C ELSE IF (ERRFLG.EQ.15) THEN WRITE(UNITO,1030) 1 'THE NUMBER OF TERMS IN EQUATION ', ERRVAL WRITE(UNITO,1030) 1 'IS GREATER THAN PARAMETER MT IN INPUT, POLFUN, POLJAC,' WRITE(UNITO,1030) 1 'POLFPT, AND SCLFPT.' WRITE(UNITO,1030) 1 'THIS CAUSES A DIMENSION PROBLEM IN COMMON BLOCKS EQUAT' WRITE(UNITO,1030) 1 'AND COEFS IN THESE ROUTINES. MAKE MT EQUAL IN ALL OF' WRITE(UNITO,1030) 1 'THESE ROUTINES AND GREATER THAN OR EQUAL TO THE MAXIMUM' WRITE(UNITO,1030) 1 'NUMBER OF TERMS IN ANY EQUATION.' GOTO 900 C ELSE IF (ERRFLG.EQ.20) THEN WRITE(UNITO,1030) 'TOTAL REQUESTS FOR LINKED LIST NODES HAS' WRITE(UNITO,1030) 'EXCEEDED THE AVAILABLE NUMBER ',ERRVAL WRITE(UNITO,1030) 'TRY ONE OF THE FOLLOWING:' WRITE(UNITO,1030) ' (1) INCREASE VALUE OF MMAXLS IN THE' WRITE(UNITO,1030) ' DRIVER AND RECOMPILE.' WRITE(UNITO,1030) ' (2) IF EXPANSION STEPS ARE IMPLEMENTED' WRITE(UNITO,1030) ' SET DLSFLG TO 0 IN THE CONFIGURATION' WRITE(UNITO,1030) ' FILE.' WRITE(UNITO,1030) ' (3) RUN THE PROGRAM FOR EACH OF THE' WRITE(UNITO,1030) ' BOXES NOT EXAMINED.' WRITE(UNITO,1030) ' (4) DECREASE THE DOMAIN TOLERANCE EPS' WRITE(UNITO,1030) ' OR THE RANGE TOLERANCE EPSF.' WRITE(UNITO,1030) ' (5) DETERMINE PROPERTIES OF' WRITE(UNITO,1030) ' THE SYSTEM WHICH MAKE IT' WRITE(UNITO,1030) ' DIFFICULT TO SOLVE.' C ELSE IF (ERRFLG.EQ.30) THEN WRITE(UNITO,1030) 'AN ILLEGAL VALUE FOR PIVFLG ',ERRVAL WRITE(UNITO,1030) 'WAS GIVEN IN THE CONFIGURATION FILE.' WRITE(UNITO,1030) 'CHANGE IT AND RERUN THE PROGRAM (OR CHANGE' WRITE(UNITO,1030) 'THE ROUTINE PVSLCT, RECOMPILE, AND RERUN).' C ELSE IF (ERRFLG.EQ.31) THEN WRITE(UNITO,1030) 1 'PVSLCT SELECTED THE ILLEGAL COORDINATE ',ERRVAL WRITE(UNITO,1030) 'TO BE BISECTED.' C ELSE IF (ERRFLG.EQ.32) THEN WRITE(UNITO,1030) 'AN ILLEGAL COORDINATE ',ERRVAL WRITE(UNITO,1030) 'WAS PASSED TO THE ROUTINE BISECT.' C ELSE IF (ERRFLG.EQ.33) THEN WRITE(UNITO,1030) 'THE COORDINATE SELECTED TO BISECT ',ERRVAL WRITE(UNITO,1030) 'IS AN INTERVAL WITH ADJACENT MACHINE' WRITE(UNITO,1030) 'NUMBERS FOR ENDPOINTS.' C ELSE IF (ERRFLG.EQ.40) THEN WRITE(UNITO,1030) 'AN ILLEGAL VALUE WAS RETURNED FROM FTEST' WRITE(UNITO,1030) 'FOR THE VARIABLE RETCON :',ERRVAL C ELSE WRITE(UNITO,1030) 'ERROR: RETURN FROM ROOTS WITH' WRITE(UNITO,1030) ' ERRFLG = ',ERRFLG WRITE(UNITO,1030) ' ERRVAL = ',ERRVAL END IF C C*********************************************************************** C C OUTPUT CURRENT BOX AND DUMP STACK C WRITE(UNITO,1000) WRITE(UNITO,1020) 'CURRENT BOX:',((CURBOX(J,K),J=1,2),K=1,N) WRITE(UNITO,1000) WRITE(UNITO,1030) 'BOXES FROM THE STACK:' WRITE(UNITO,1000) 100 IF (DEPTH.EQ.0) GOTO 890 CALL POP(N,CURBOX,MAXDP,DEPTH,LEVEL,STACK,STKLVL,ERRFLG,ERRVAL) WRITE(UNITO,1020) ' ',((CURBOX(J,K),J=1,2),K=1,N) WRITE(UNITO,1000) GOTO 100 C 890 WRITE(UNITO,1030) 'END OF LIST OF STACK BOXES' WRITE(UNITO,1000) C 900 CONTINUE RETURN C 1000 FORMAT(1X) 1010 FORMAT(1X,A,1X,D30.22) 1015 FORMAT(1X,A,/2X,D30.22) 1020 FORMAT(1X,A,50(/2(2X,D30.22))) 1025 FORMAT(1X,A,50(/2X,D30.22)) 1030 FORMAT(1X,A,I8) 1040 FORMAT(I8) 1050 FORMAT(D22.15) 1080 FORMAT('+',2X,I3,A,I3) 1082 FORMAT(2X,'LEFT IS ',L1,' AND RIGHT IS ',L1) 1085 FORMAT(20(/' | ',6(I3,L1,I3,' | '))) 1090 FORMAT(D30.22) C END C*********************************************************************** C*********************************************************************** SUBROUTINE EXPAND(N,R,DL,DLB) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C Written by: C C R. Baker Kearfott C C and C C Manuel Novoa III C C Department of Mathematics C U.S.L. Box 4-1010 C Lafayette, LA 70504 C C September 29, 1987 C C Part of the generalized bisection package C C*********************************************************************** C C Called by -- ROOTS C C*********************************************************************** C C Function -- C C This routine performs step 4 of Algorithm 2.5 in R. B. Kearfott, C 'Abstract generalized bisection and a cost bound', Math. Comput. 49, C 179 (July, 1987), pp. 187-202. C C The purpose of this algorithm is to reduce redundancies and make C the algorithm more efficient when roots happen to occur near C boundaries of boxes. This algorithm is moderately successful at C this. Future versions may be even more so, if they take account C of the convergence rate of the interval Newton method to expand C boxes before they become small. C C*********************************************************************** C C Argument declarations -- C INTEGER N DOUBLE PRECISION R DOUBLE PRECISION DL(2,N) DOUBLE PRECISION DLB(2,N) C C*********************************************************************** C C Argument descriptions -- (INPUT = set on entry and not alterable) C (OUTPUT = to be set by the routine) C (I/O = set on entry but alterable) C C N is the number of equations and variables. C (INPUT) C C R is the number R defined in Definition 2.2 (vi) of Kearfott C (ibid.). C (INPUT) C C DL DL(1,I) is the lower limit of the I-th coordinate interval C and DL(2,I) is the upper limit of the I-th coordinate C interval, for I between 1 and N. C (INPUT) C C DLB is an array with the same structure as DL, but which will C represent the expanded box in Algorithm 2.5 of Kearfott, C (ibid.). C (OUTPUT) C C*********************************************************************** C C Internal variable declarations -- C DOUBLE PRECISION CI INTEGER I C C*********************************************************************** C C Common block declarations -- none C C*********************************************************************** C C Fortran-supplied functions and subroutines -- none C C*********************************************************************** C C Package-supplied functions and subroutines -- none C C*********************************************************************** C C User-supplied functions and subroutines -- none C C*********************************************************************** C C I/O functions -- none C C*********************************************************************** C C Internal constant declarations -- none C C*********************************************************************** C C Beginning of executable statements -- C DO 10 I = 1, N CI = (DL(1,I)+DL(2,I)) / 2D0 DLB(1,I) = CI + (2D0/R) * (DL(1,I) - CI) DLB(2,I) = CI + (2D0/R) * (DL(2,I) - CI) 10 CONTINUE C RETURN END C*********************************************************************** C*********************************************************************** SUBROUTINE FTESTH(N,X,XPT,DIAM,ERRFLG,ERRVAL,UNKNWN,SIGRT, 1 FUNC,JAC,EPSBOX,EPSF,RETCON,NFCAL1,NFCAL2,NJCALL, 2 MAXDP,DEPTH,LEVEL,STACK,STKLVL, 3 IP,DP,DPIVOT,IPIVOT,TSTITR, 4 LNDWRK,DWORK,LNIWRK,IWORK,LNLWRK,LWORK) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C Written by: C C R. Baker Kearfott C C and C C Manuel Novoa III C C Department of Mathematics C U.S.L. Box 4-1010 C Lafayette, LA 70504 C C September 29, 1987 C C Part of the generalized bisection package. C C*********************************************************************** C C Called by -- ROOTS C C*********************************************************************** C C Function -- C C This subroutine partitions the work storage for use in the C arrays in the root inclusion test routine HNSNG (to allow C changes in HNSNG without necessarily modifying the overall C algorithm ROOTS). It also calls the routine PVSLCT; that routine C selects the coordinate direction to be bisected next, and also C indicates which of the two resulting boxes is most likely to C contain roots. C C*********************************************************************** C C Argument declarations -- C INTEGER N DOUBLE PRECISION X(2,N), XPT(N), DIAM INTEGER ERRFLG, ERRVAL LOGICAL UNKNWN, SIGRT EXTERNAL FUNC, JAC DOUBLE PRECISION EPSBOX, EPSF INTEGER RETCON, NFCAL1, NFCAL2, NJCALL, MAXDP, DEPTH INTEGER LEVEL DOUBLE PRECISION STACK(2,N,MAXDP) INTEGER STKLVL(MAXDP) INTEGER IP LOGICAL DP LOGICAL DPIVOT(MAXDP) INTEGER IPIVOT(MAXDP), TSTITR(MAXDP) INTEGER LNDWRK DOUBLE PRECISION DWORK(LNDWRK) INTEGER LNIWRK, IWORK(LNIWRK), LNLWRK LOGICAL LWORK(LNLWRK) C C*********************************************************************** C C Argument descriptions -- (INPUT = set on entry and not alterable) C (OUTPUT = to be set by the routine) C (I/O = set on entry but alterable) C C N is the number of variables and equations. C (INPUT) C C X On entry, this array contains the box to be tested for C roots. On exit, this array either contains the box which C was to be tested or else a sub-box in which all roots must C lie. C (I/O) C C XPT This array is used to store a point approximation to C a root in X (Usually, this is the midpoint of X). C (I/O) C C DIAM is the diameter of the box stored in X. C (I/O) C C ERRFLG is the error type. FTESTH sets ERRFLG to 0 if there is no C error, but sets ERRFLG to another value if certain errors have C occurred. (See subroutine ERROUT for more detailed C information.) C (OUTPUT) C C ERRVAL contains the index of an error which has occurred. C (See subroutine ERROUT for more detailed information.) C (OUTPUT) C C UNKNWN Is set to .TRUE. to signal 'unknown' for the C root-inclusion test, and is set to .FALSE. if the C root-inclusion test is true or false. C (OUTPUT) C C SIGRT is set to .TRUE. if the root-inclusion test signals a C unique root, and is set to .FALSE. if the root-inclusion C test signals no root. (This flag is meaningful only if C UNKNWN = .FALSE.) C (OUTPUT) C C FUNC is the name of the routine to compute the interval values C of the function. (The system routine is POLFUN.) C (INPUT) C C JAC is the name of the routine to compute the interval values C of the Jacobian matrix. (The system routine is POLJAC.) C (INPUT) C C EPSBOX is the minimum allowable width (in the norm defined DIAMCP) of C an unresolved box. C (INPUT) C C EPSF is the range tolerance; if it is certain that each component C of the vector function has absolute value less than EPSF C within a box, then that box is signalled as containing a C root. C (INPUT) C C RETCON is the return condition from FTEST. Its values are as C follows. C C 0: The preconditioner matrix could not be computed. C C 1: The root-inclusion test was inconclusive. C C 2: The diameter of the box is less than the domain C tolerance EPS/16. C C 10: Each component of the function has absolute value less C than the range tolerance EPSF over the entire box. C C 11: The box contains a unique root. C C 20: Zero is not in the interval value of the function. C C 21: The image of the box under the interval Newton method C has null intersection with the box. C C (OUTPUT) C C NFCAL1 is the number of interval function calls made to possibly C determine that the zero vector is not in the range. C (I/O) C C NFCAL2 is the number of function calls made for the interval Newton C method. C (I/O) C C NJCALL is the number of interval and scalar Jacobian calls C (I/O) C C MAXDP is the maximum depth of the stack of boxes yet to be C considered. C (INPUT) C C DEPTH is the current depth of the stack of boxes yet to be examined. C (I/O) C C LEVEL is the current level in the binary search tree. C (I/O) C C STACK is temporary storage space for the stack of boxes yet to be C examined. C (I/O) C C STKLVL is an integer stack whose entries correspond to the entries in C STACK. In it is stored the level (in the binary tree C from bisection) at which the corresponding box in STACK C occurred. C (I/O) C C IP is set to the index of the coordinate PVSLCT has chosen to be C bisected. C (OUTPUT) C C DP is set to .TRUE. if the left box formed from bisection is to C be considered next, and is set to .FALSE. if the right box is C to be considered next. (The left box is that box formed by C replacing the right endpoint of the bisected interval by its C midpoint.) (This choice is made in PVSLCT.) C (OUTPUT) C C DPIVOT, IPIVOT, and TSTITR are for additional stack information C associated with boxes stored in STACK. Their entries C are indexed on the level in the binary tree instead of the C depth of the stack. C (I/O) C C LNDWRK is the length of the double precision work array. C (INPUT) C C DWORK is the double precision work array. C (I/O) C C LNIWRK is the length of the integer work array. C (INPUT) C C IWORK is the integer work array. C (I/O) C C LNLWRK is the length of the logical work array. C (INPUT) C C LWORK is the logical work array. C (I/O) C C*********************************************************************** C C Internal variable declarations -- C INTEGER NUMITR C C*********************************************************************** C C Internal variable descriptions -- C C NUMITR is the number of iterations of the interval Newton method. C C*********************************************************************** C C Common block declarations -- none C C*********************************************************************** C C Fortran-supplied functions and subroutines -- none C C*********************************************************************** C C Package-supplied functions and subroutines -- HNSNG, PVSLCT C C*********************************************************************** C C User-supplied functions and subroutines -- none C C*********************************************************************** C C I/O functions -- none C C*********************************************************************** C C Internal constant declarations -- C INTEGER ISCLF, ISCLJC, IFINT, IJACIN, IYMAT INTEGER IWIV0, IWIV1, IWIV2, ITMPM1, ITMPM2 INTEGER IDEND, IIPVT, IIEND, ILWORK, ILEND INTEGER LUNKNW, UUNKNW, LTRUE, UTRUE, LFALSE, UFALSE C C*********************************************************************** C C INTERNAL CONSTANT DESCRIPTIONS -- C C ISCLF a pointer to the starting address of the scalar function C value in the double precision work array C C ISCLJC a pointer to the starting address of the scalar Jacobian C matrix value in the double precision work array C C IFINT a pointer to the starting address of the interval function C value in the double precision work array C C IJACIN a pointer to the starting address of the interval Jacobian C matrix value in the double precision work array. C C IYMAT a pointer to the starting address of the prescaling matrix C in the double precision work array C C IWIV0, IWIV1, and IWIV2 are pointers to starting addresses of work C arrays for interval N-vectors in the double precision work C array. C C ITMPM1 and ITMPM2 are pointers to starting addresses of work C arrays for intervan N by N matrices in the double precision C work arrays. C C IDEND a pointer to the last used location in the double C precision work array C C IIPVT a pointer to the starting address of an integer work C vector in the integer work array C C IIEND a pointer to the last used location in the integer work C array C C ILWORK a pointer to an address in the logical work array C (currently not used) C C ILEND a pointer to the la