# to unbundle, sh this file (in an empty directory) echo README 1>&2 sed >README <<'//GO.SYSIN DD README' 's/^-//' -ve08ad.f is Fortran source for VE08AD, a subroutine for minimizing - partially separable objective functions. Comments near - the beginning explain how to use VE08AD. - -ve08ts.f is Fortran source for an example calling program for VE08AD. - -ve08ts.out contains output obtained by running the test program on an - SGI 4D/240S computer (IEEE arithmetic). - -VE08AD and the test program were written and supplied (with permission -from Harwell) by - - Phillipe L. Toint - FUNDP - 61, Rue de Bruxelles - B-5000 Namur - Belgium - - na.toint@na-net.ornl.gov - phtoint@bnandp51.bitnet - phtoint@cc.fundp.ac.be //GO.SYSIN DD README echo ve08ad.f 1>&2 sed >ve08ad.f <<'//GO.SYSIN DD ve08ad.f' 's/^-//' -*From CUNYVM.CUNY.EDU!phtoint%NEWTON Mon Nov 5 10:35 +01 1990 - SUBROUTINE VE08AD(N, NS, INVAR, NVAR, - * ELFNCT, RANGE, XLOWER, XUPPER, - * X, FX, EPSIL, EBOUND, NGR, NIT, FKNOWN, FLOWBD, - * RELPR, DIFGRD, RESTRT, TESTGX, HESDIF, STMAX, STEPL, ISTATE, - * IPDEVC, IPFREQ, IPWHAT, LWK, WK, INFO, IFLAG) -C -C -C -C -C*********************************************************************** -C -C -C ***************************************************************** -C ***************************************************************** -C ** ** -C ** VE08AD ** -C ** ** -C ** A ROUTINE FOR SOLVING PARTIALLY SEPARABLE MINIMIZATION ** -C ** PROBLEMS WITH POSSIBLE UPPER AND LOWER BOUNDS ON THE ** -C ** VARIABLES. ** -C ** ** -C ***************************************************************** -C ***************************************************************** -C -C ************** DOUBLE PRECISION VERSION ************* -C -C PROGRAMMING: PH.L. TOINT, NAMUR (BELGIUM), 1982. -C PROGRAM VERSION NUMBER: 1.2 (06.07.1983) -C MODIFIED FOR THE HARWELL SUBROUTINE LIBRARY IN NOVEMBER 1983 -C WITH THE ASSISTANCE OF IAIN DUFF (EXT 2670) AND JOHN REID. -C THE ORIGINAL NAME OF THIS PACKAGE WAS PSPMIN. -C -C LAST MINOR UPDATE : 05.03.1985 (PH.L. TOINT) -C -C -C ********************************************************************* -C * * -C * PURPOSE * -C * * -C ********************************************************************* -C -C To minimize an objective function consisting of a sum of -C finite element' functions, each of which involves only a few -C variables or whose second derivative matrix has a low rank for -C other reasons. -C Bounds on the variables and known values may be specified. -C The routine is especially effective on problems involving a -C large number of variables. -C -C The objective function has the form -C -C ns -C F(x) = SUM f ( x ) -C k=1 k k -C -C where -C -C x = ( x , x , ..., x ) -C 1 2 n -C -C and where -C -C x are small subsets of x. -C -C The user is required to write a subroutine that -C calculates each function {f SUB k} and (optionally) its -C gradient -C df -C ns k -C g = SUM --- -C k k=1 dx -C k -C -C When code for the gradient is supplied (and this usually -C results in faster and more reliable execution), there are -C facilities for automatically checking it during the first -C iteration. -C -C There are flexible re-start facilities for problems that do -C not differ too substantially from a problem previously solved. -C Using these can sometimes lead to dramatically reduced -C execution times. -C -C If there is a linear transformation of variables such that -C one or more of the element functions depend on fewer -C transformed variables than the number of original variables -C they involve, this may be specified. This too can lead -C to substantially improved computing time. -C -C There are facilities for the user: -C - to specify the termination criterion, -C - to influence the step size in the line-search procedure, -C - to choose whether to calculate the first approximations to -C the Hessian matrices -C 2 -C d f -C ns k -C H = SUM ---- -C k k=1 2 -C dx -C k -C by finite differences. -C -C ********************************************************************** -C * * -C * SOURCE * -C * * -C ********************************************************************** -C -C VE08 is a routine from the Harwell Subroutrine Library (Harwell -C Laboratory, Harwell, Oxfordshire, United Kingdom). -C -C Programming: Ph. L. TOINT, FUNDP, Namur, Belgium -C December 1983. -C -C -C ********************************************************************** -C * * -C * USE * -C * * -C ********************************************************************** -C -C The Argument List -C ----------------- -C -C VE08AD(N,NS,INVAR,NVAR,ELFNCT,RANGE,XLOWER,XUPPER, -C X,FX,EPSIL,EBOUND,NGR,NIT,FKNOWN,FLOWBD,RELPR, -C DIFGRD,RESTRT,TESTGX,HESDIF,STMAX,STEPL,ISTATE, -C IPDEVC,IPFREQ,IPWHAT,LWK,WK,INFO,IFLAG) -C -C -C N is an INTEGER variable, that must be set by the user to -C the number of variables. It is not altered by the routine. -C Restriction: N > 0. -C -C NS is an INTEGER variable, that must be set by the user to ns, -C the number of element functions. It is not altered by the -C routine. -C Restriction: NS > 0. -C -C INVAR is an INTEGER array containing the indices of the variables -C in the first element, followed by those in the second -C element, etc. See below for an example. -C It is not altered by the routine. -C -C NVAR is an INTEGER array of length at least NS+1, whose k-th -C value is the position of the first variable of element -C function k, in the list INVAR. In addition, NVAR(NS+1) -C should be equal to the total length of INVAR plus one. -C See below for an example. -C It is not altered by the routine. -C -C ELFNCT is the name of a user supplied subroutine that -C computes the values of the element functions at a given point. -C See complete description below. This name should be -C declared EXTERNAL in the calling program. -C -C RANGE is the name of a user supplied subroutine that performs -C various operations related to the true range of the element -C function second derivative matrices. See complete description -C below. This name should be declared EXTERNAL in the -C calling program. -C -C XLOWER is the name of a user supplied DOUBLE PRECISION -C function routine that returns the lower bounds on the bounded -C variables. See complete description below. This name should -C be declared `EXTERNAL in the calling program. -C -C XUPPER is the name of a user supplied DOUBLE PRECISION -C function routine that returns the upper bounds on the bounded -C variables. See complete description below. This name should -C be declared `EXTERNAL in the calling program. -C -C X is a DOUBLE PRECISION array of length N which must be -C set by the user to the value of the variables at the starting -C point. On exit, it contains the values of the variables at the -C best point found (usually the solution). -C -C FX is a DOUBLE PRECISION variable. If the user knows the -C value of the objective function at the starting point, -C he must set its value in FX and set FKNOWN to true. -C Otherwise, the user should set FKNOWN to false and need -C not set FX. On exit, FX contains the value of the -C objective function at the point X. -C -C EPSIL is a DOUBLE PRECISION variable, that must be set by -C the user to a measure of the accuracy (in the sense of the -C Euclidean norm of the projected gradient) required to stop the -C minimization procedure. This value is passed to the -C termination routine VE08SD (see below). It is not altered by -C the routine. -C -C EBOUND is a DOUBLE PRECISION variable, that must be set by -C the user to a tolerance above which a bound is considered to -C be nearly active. If it is negative, the default value -C 10**(-5) is used. This tolerance is used in the -C antizigzaging device of (Bertsekas 1982), and has no -C influence on the precision with which the bound -C constraints will be satisfied at the solution. -C It is not altered by the routine. -C -C NGR is an INTEGER array of length 2. The first element -C NGR(1) must be set by the user to the maximum number of calls -C to ELFNCT allowed for the minimization and is not altered. A -C suitable value will depend on the problem, and may be chosen in -C the first instance between NIT(1) and 3*NIT(1). It may be -C necessary to increase it if the problem is highly nonlinear. -C NGR(2) need not be set by the user and will contain, on exit, -C the actual number of calls to the routine ELFNCT that were made -C during the routine's execution. -C -C NIT is an INTEGER array of length 2. The first element -C NIT(1) must be set by the user to the maximum number of -C iterations that is allowed for the minimization and is not -C altered. The second element NIT(2) need not be be set by the -C user, and will contain, on exit, the actual number of -C minimization iterations performed by the routine. -C -C FKNOWN is a LOGICAL variable. If it is set to true, FX must -C contain the function value at the starting point X, the first -C NS positions of the working vector WK must contain the -C individual values of the element functions at X, and the -C NVAR(NS+1)-1 following positions in WK must contain the -C gradients of the element functions. This option must be used -C only when the gradients of all element functions are -C analytically available. It is not altered by the routine. -C -C FLOWBD is a DOUBLE PRECISION variable, that must be set by -C the user to a lower bound on the optimal value of F(x). -C For instance, if it is known that F(x) is always positive, -C it can help convergence to pass this information to the routine. -C If no bound is known, then a large negative value should be -C used. It is not altered by the routine. -C -C RELPR is a DOUBLE PRECISION variable, that should be set by -C the user to the machine rounding error unit, i.e. the -C smallest positive {epsilon} such that {1+ epsilon} is still -C distinguishable from 1. The user may set RELPR to a negative -C number, and the routine will compute a value for it. If it is -C positive on input, then it is not altered by the routine. -C Otherwise, it contains, on exit, the value computed by the -C routine. -C -C DIFGRD is a DOUBLE PRECISION variable, that should be set -C by the user to the relative step size that is to be used in the -C first gradient estimation by differences. If the user does not -C know a suitable value, he may set it to a negative number, and -C the routine will use the square root of RELPR. It is unaltered -C if positive, and set to the square root of RELPR otherwise. -C -C RESTRT is a LOGICAL variable that must be set to true by the -C user if VE08 is to be restarted and to false otherwise. In -C the first case, approximations to the element Hessians must be -C stored in the vector WK from position NS+NVAR(NS+1) onwards -C (see the description of WK). This preservation of -C second-order information may be rather efficient. Examples for -C which it is useful are for multiple criteria optimization, when -C one wishes to reoptimize a weighted sum of objectives with -C new weights and for a discretized problem with varying levels of -C mesh coarseness. RESTRT is not altered by the routine. -C -C TESTGX is a LOGICAL variable that must be set to true by the user -C if he provides analytical gradients in ELFNCT and wishes them -C to be tested for accuracy at the starting point by comparing -C their values to difference approximations, and to false -C otherwise. If RESTRT is true then TESTGX is reset to false -C and it is as if it were false on entry. Otherwise, it is not -C altered by the routine. -C -C HESDIF is a LOGICAL variable that must be set to true by -C the user if the initial Hessian approximations are to be -C computed by differences in the element gradient values, and to -C false otherwise. This option is only available if gradients are -C calculated in ELFNCT. When this option is not used, the -C element Hessians are initialized to a correctly scaled multiple -C of the identity matrix. If RESTRT is true then HESDIF is -C reset to false and it is as if it were false on entry. -C Otherwise HESDIF is not altered by the routine. -C -C STMAX is a DOUBLE PRECISION variable, that should be set -C by the user to the maximum steplength that is allowed during the -C minimization process. If the user does not know a suitable -C value, STMAX may be set to any negative number; the routine -C will then use a very large value. It is unaltered by the -C routine. -C -C STEPL is a DOUBLE PRECISION array of length 2. The first -C element STEPL(1) should be set by the user to an upper bound -C on the length of the first step taken by the method. If -C such a bound is unknown, STEPL(1) may be set to any -C negative number; the routine will then use a large default -C value. STEPL(1) is not altered by the routine. -C The second element STEPL(2) need not be set by the user, -C and on exit will contain the length of the last step taken -C by the algorithm. -C -C ISTATE is an INTEGER array of length at least N+NS. It must -C be set as follows. For I=1,N, ISTATE(I) must be set to -1 if -C the I-th variable is unconstrained, to 0 if it is fixed, and -C to 1 if it is bounded. The values of fixed variables are never -C changed during the minimization process. (They may be useful, -C for example, when dealing with discretized problems having -C boundary conditions; the variables on the boundary can then be -C fixed, but otherwise treated as ordinary variables, which avoids -C the need for special expressions in the elements touching the -C boundary.) For I=N+1, ...,N+NS, ISTATE(I) must be set to 1 -C if the gradient of the (I-N)th element function is supplied -C by the user, or to -1 if it has to be estimated by differences. -C (The difference approximations are computed by a modification of -C the method proposed by (Stewart 1967). Although efficiently -C computed, the use of approximated gradients can sometimes -C deteriorate the overall performance of VE08; analytical -C gradients are always preferable when available.) ISTATE is -C used as workspace, so must be reset on a second entry. -C -C IPDEVC is an INTEGER variable, that must be set by the user -C to the output device unit number for printing of messages. It -C is not altered by the routine. -C -C IPFREQ is an INTEGER variable, that must be set by the user -C to specify the frequency of output by the routine VE08: -C IPFREQ < 0: no output is generated, -C IPFREQ = 0: output only at first and last iteration, -C IPFREQ > 0: output every IPFREQ iteration. -C This argument is not altered by the routine. -C -C IPWHAT is an INTEGER variable, that must be set by the user to -C specify the amount of output generated when output occurs: -C IPWHAT = 0: iteration count, function value, norm of the -C active gradient and number of function calls, -C IPWHAT = 1: + step selection iterations and linesearch -C parameter, -C IPWHAT = 2: + vector of variables, -C IPWHAT = 3: + gradient vector, -C IPWHAT = 4: + last step of the algorithm, -C IPWHAT = 5: + element Hessian approximations. -C This argument is not altered by the routine. -C -C LWK is an INTEGER array of length 2. The first element LWK(1) -C must be set by the user to the length of the working array WK -C and is not altered. The second element LWK(2) need not be set -C by the user and contains, on exit, the length of the array -C needed to store the approximating element matrices. The length -C of the workspace WK that is required by the routine VE08 -C varies from problem to problem. It essentially depends on the -C amount of storage needed for the element Hessian -C approximations. Let -C -C NV = NVAR(NS+1)-1, -C NE = maximum of NVAR(I+1)-NVAR(I), I=1,NS, -C NQ = maximum of "2*N" and NV, -C -C then the minimum required storage is -C -C LWK(1)=NS+NV+4*N+NQ+3*NE+LWK(2) -C -C The quantity NE in the relation above is the maximum number -C of variables appearing in any element function. -C -C WK is a DOUBLE PRECISION working array of length LWK(1). -C On sucessful termination, it contains a list of the element -C function values, gradients and approximating Hessians and is -C organized as follows : WK(I) contains, for -C -C I=1,...,NS, the values of {f (x)}, i = I and X in X. -C i -C I=NS+1,...,NS+NVAR(NS+1)-1, the NS successive element -C gradient vectors, one after the other, -C I=NS+NVAR(NS+1),...,NS+NVAR(NS+1)+LWK(2)-1, the lower -C triangular parts of the element Hessian approximations, -C stored one after the other in row-wise fashion. -C -C INFO is an INTEGER variable. It need not be set by the user -C and contains, on exit under an error condition, information -C about the error (see below for further details). -C -C IFLAG is an INTEGER variable, that need not be set by the -C user and contains, on exit, the exit condition of the routine. -C If this flag is greater or equal to 11, an error has been -C detected by the routine (see below for further details). -C -C -C Subroutine ELFNCT -C ----------------- -C -C We now describe the calling sequence of the routine that the -C user must provide to evaluate the element functions (and -C possibly their gradients). -C -C It is worth mentioning that the points at which element -C functions are computed are almost always feasible, that is they -C satisfy whatever bounds there are. The only cases when they -C may not is when gradients are to be estimated by differences; -C in this case, the function value may be required at some points -C that fail to satisfy one or more bounds by a small amount. -C -C The routine has the following argument list: -C -C SUBROUTINE ELFNCT(K,X,FX,GX,NDIMK,NS,JFLAG,FMAX,FNOISE) -C -C K is an INTEGER variable, that contains the index of the -C element function to be computed. It must be left unchanged by -C the routine. -C -C X is a DOUBLE PRECISION array of length NDIMK, that -C contains the values of the variables of element K, in the -C order in which they appear in the vector INVAR. It must be -C left unchanged by the routine. -C -C FX is a DOUBLE PRECISION variable, that must be set by -C the routine to the value of the element function at the -C point X. -C -C GX is a DOUBLE PRECISION array of length NDIMK. If -C JFLAG has value 2, the routine must set GX to the components -C -C df -C k -C --- -C dx -C k -C -C of the gradient of the element function at the point X. -C If JFLAG has other values, GX need not be set. -C -C NDIMK is an INTEGER variable, that contains the number of -C variables in element K. It is equal to NVAR(K+1) - NVAR(K) -C and must be left unchanged by the routine. -C -C NS is an INTEGER variable, that contains the total number -C of elements in the problem. It must be left unchanged by the -C routine. -C -C JFLAG is an INTEGER variable, that contains a code to -C describe the information expected on return from the routine -C ELFNCT: -C JFLAG = 1: function value is needed, -C JFLAG = 2: gradient and function values are needed. -C JFLAG < 0: if ELFNCT sets a negative value in JFLAG, -C VE08 will terminate abnormally, with exit -C condition 25 and INFO parameter equal to JFLAG. -C JFLAG = 0: if ELFNCT sets JFLAG to 0, no further calls of -C ELFNCT will be made to complete the current -C evaluation of F and the steplength will be -C halved; this can be useful if F value is already -C too large compared to FMAX, and if the user -C does not want to estimate the remaining element -C functions because of cost. -C -C If none of these events is desired, JFLAG should be left -C unchanged. The possibility of tampering with the linesearch -C using a zero return value for JFLAG should be used with caution, -C especially when no bounds are imposed and analytical -C gradients available, as it may reduce the overall efficiency. -C -C FMAX is a DOUBLE PRECISION variable, that contains the -C maximum function value that will be accepted by the -C algorithm as satisfying the current linesearch criteria. It must -C be left unchanged by the routine. -C -C FNOISE is a DOUBLE PRECISION variable that the routine must -C set to an estimate of the noise (roundoff) present in the -C computation of the element function and element gradient. -C -C An example of the use of ELFNCT is shown below. -C -C -C Subroutine RANGE -C ---------------- -C -C We now turn to describing the way in which the user can -C pass to the algorithm information about a linear -C transformation that reduces the number of independent variables -C upon which each element function depends. For example the -C function -C -C 2 2 -C x + ( x - x ) -C 1 2 3 -C -C depends on the two variables y = x and y = x - x . -C 1 1 2 2 3 -C In such cases, the convergence of the algorithm may be substantially -C enhanced. For each element function k, the user must express -C the dependence in terms of a full-rank matrix U , with fewer rows -C k -C than columns, that maps the variables x , i = INVAR(NVAR(k)) up -C i -C to INVAR(NVAR(k+1)-1), that are involved in element function k -C to the smaller set. In the example of this paragraph -C -C ( 1 0 0 ) -C U = ( 0 1 -1 ) -C k -C -C In such cases, the Hessian matrix can be held as -C -C T -C H = U C U -C k k k k -C -C where C is a square symmetric matrix of order the number of -C k -C rows in U . -C k -C -C The user should therefore provide a routine called RANGE, -C that is called by VE08, and has the following argument list: -C -C SUBROUTINE RANGE(K,MODE,W1,W2,NDIMK,NSUBK,NS) -C -C K is an INTEGER variable, that contains the index of the -C element function. It must be left unchanged by the routine. -C -C MODE is an INTEGER variable, that contains a code for -C the work to be accomplished by the routine : -C MODE=1: the user should provide the vector w2 -C such that -C -C U w1 = w2. -C k -C MODE=2: the user should provide the vector w2 -C with the smallest norm such that -C -C U w2 = w1. -C k -C -C Equivalently, w2 is the result of the -C application of the Moore-Penrose generalized -C inverse of U to w1. -C k -C MODE=3: the user should provide the vector w2 -C such that -C -C T -C U w1 = w2. -C k -C -C MODE=4: the user should provide the vector w2 -C such that -C -C T -C U w2 = w1. -C k -C -C MODE must be left unchanged by the routine. -C -C W1 is a DOUBLE PRECISION array containing the input vector w1. -C Itmust be left unchanged by the routine. -C -C W2 is a DOUBLE PRECISION array that must be set by the -C routine to the result vector w2, related to U and w1 -C k -C as required by the argument MODE. -C When NSUBK=NDIMK, W2 must be set equal to W1. -C -C NDIMK is an INTEGER variable, containing the number of -C variables in the element. It must be left unchanged by the -C routine. -C -C NSUBK is an INTEGER variable that must be set by RANGE to -C the number of rows in U. -C -C NS is an INTEGER variable, containing the number of elements -C of the problem. It must be left unchanged by the routine. -C -C If the user does not know such information it is always -C possible to use the following empty routine RANGE: -C -C SUBROUTINE RANGE(K,MODE,W1,W2,NDIMK,NSUBK,NS) -C DOUBLE PRECISION W1(1),W2(1) -C NSUBK=NDIMK -C DO 1 I=1,NDIMK -C W2(I)=W1(I) -C 1 CONTINUE -C RETURN -C END -C -C The use of this empty routine is nevertheless not -C recommended, especially when the gradients are not available -C analytically, or when the option HESDIF=.FALSE. is used in -C VE08. It may, in these cases, result in much slower -C convergence, or possibly in no convergence at all. -C -C -C Functions XLOWER and XUPPER -C --------------------------- -C -C The fact that the I-th variable of the problem is bounded -C is signalled to the routine VE08 by ISTATE(I) being equal -C to 1. VE08 will then need to known the actual values of the -C bounds on this variable. This information is provided by two -C user-supplied double precision functions XLOWER(I) and XUPPER(I). -C -C XLOWER(I) returns the value of the lower bound on the -C I-th variable, while XUPPER(I) returns the value of the -C upper bound on this variable. Their specification is as -C follows: -C -C DOUBLE PRECISION FUNCTION XLOWER(I) -C DOUBLE PRECISION FUNCTION XUPPER(I) -C -C Remarks: -C -C 1. The functions XLOWER and XUPPER are only called for -C arguments I such that ISTATE(I)=1, i.e. they are only -C called for bounded variables. -C -C 2. When the I-th variable is bounded, both XLOWER and -C XUPPER are called for this variable : each bounded -C variable must be bounded below and above. If the problem -C only incorporates one of these bounds, the other should -C be supplied using a very small or a very large -C constant. -C -C 3. It is more efficient to declare a variable "fixed" -C (ISTATE(I)=0) than to constrain it by equal lower and -C upper bounds. -C -C -C Error Messages -C -------------- -C -C When the exit condition of VE08, i.e. IFLAG, is greater -C than 10, this means that it has detected something going wrong -C in the computation. The routine is terminated with IFLAG set -C to an appropriate error index and INFO to a (sometimes) -C meaningful value. A complete list of these errors, together with -C the associated value of IFLAG and INFO is given below. -C -C IFLAG=11: The value of variable N is not positive. INFO = N. -C IFLAG=12: The machine precision constant RELPR is out of -C range. INFO is meaningless. -C IFLAG=13: The last call to ELFNCT returned a negative -C value for FNOISE. INFO is meaningless. -C IFLAG=14: The number of element functions NS is not positive. -C INFO=NS -C IFLAG=15: The starting position for the internal variables -C of the {k}th element NVAR(k) is not positive. -C INFO = {k}. -C IFLAG=16: The number of variables internal to the k-th -C element is not positive. INFO = k. -C IFLAG=17: NSUBK {GT} NDIMK on return from RANGE for -C element k. INFO = k. -C IFLAG=18: The maximum number of iterations has been -C reached. INFO = maximum number of iterations. -C IFLAG=19: The initial status of the k-th element -C function is incorrect: the only allowed values -C for ISTATE(N+k) are 1 (derivatives available) -C or -1 (derivatives not available). INFO = k. -C IFLAG=20: The initial status of the I-th variable is -C incorrect because the only allowed values for -C ISTATE(I) are 1 (bounded), 0 (fixed), or -1 (free). -C INFO=I. -C IFLAG=21: The accumulated dimension of the elements, -C NVAR(NS+1)-1, is not positive. INFO = NVAR(NS+1)-1. -C IFLAG=22: The index of the I-th variable that appears in -C the vector INVAR is not positive, or greater -C than N. INFO = I. -C IFLAG=23: The total available work space provided in -C the vector WK, of length LWK(1), is too small. -C INFO = minimum necessary length of the vector WK. -C IFLAG=24: The bounds on the I-th variable are inconsistent. -C INFO = I. -C IFLAG=25: The user has stopped the minimization -C procedure by setting the element function flag -C to a negative value. INFO = value of the -C element function flag. -C IFLAG=26: {F()} is unbounded because the I-th variable -C appears linearly and is unbounded. INFO = I. -C IFLAG=27: The directional derivative at the beginning of a -C linesearch is non-negative. This can be caused -C by an incorrect analytical gradient. INFO is meaningless. -C IFLAG=28: The linesearch step became too small after I -C trials. This can be caused by an incorrect -C analytical gradient or a too noisy function or an -C exceptionally non-linear function. (This error is -C controlled by the variable MAXLS, representing the -C maximum number of linesearch trials, which is set at -C the beginning of the code of VE08>). This can also -C be caused near the solution, when the accuracy -C requirement for termination is too high. In this -C case, further progress of the algorithm seems -C unlikely on this machine. INFO = I. -C IFLAG=29: The algorithm was stopped because the -C difference between two successive function values -C along the current search direction is -C insignificant compared to the noise on these -C function values. This type of exit usually -C occurs quite close to the solution, and is mainly -C due to a too high accuracy requirement for -C termination. In this case, further progress of -C the algorithm on this machine is unlikely. INFO -C is meaningless. -C IFLAG=30: An error was detected when scaling the I-th -C initial element Hessian matrix at the first step. This -C is usually caused by an incorrect gradient, or by -C incorrect problem structure specifications. INFO = I. -C IFLAG=31: The last call to the routine VE08SD (see below) -C returned a value for IFLAG outside the range 0 .LE. -C IFLAG .LE. 10. INFO = returned value of IFLAG. -C IFLAG=32: The check on the correctness of the analytical -C gradient at the starting point, has discovered a -C possible error in this computation. INFO is -C meaningless. -C IFLAG=33: The lower bound on the objective function is not -C consistent with the function value at the starting -C point. INFO is meaningless. -C IFLAG=34: The parameter IPWHAT, which controls the amount -C of output required by the user, is incorrect -C (outside the range 0 .LE. IPWHAT .LE. 5). INFO = IPWHAT. -C IFLAG=35: The routine stopped because the norm of the -C computed search direction is smaller than the -C given machine precision times the size of the -C current approximate solution. This error -C usually occurs when too much accuracy is required -C by the user for termination. INFO = iteration -C number when VE08 was stopped. -C IFLAG=36: The routine is stopped because 5 successive very -C large steps ({GT_0.999 times STMAX}) were -C taken. This is interpreted as a sign of -C divergence of the algorithm. This usually happens -C when the problem's minimum is at infinity. -C INFO = iteration number when VE08 was stopped. -C -C The values of IFLAG between 1 and 10 are left free for -C description of normal termination (see the description of the -C routine VE08SD below). -C -C Termination Criterion -C --------------------- -C -C One of the features of VE08AD is to allow the interested -C user to provide a problem-suited stopping criterion. There -C may be cases where the classical tests on the Euclidean length -C of the gradient are rather inadequate. One may wish to use -C another norm, or another stopping rule altogether. This can be -C done by replacing the routine VE08SD by another routine -C VE08SD incorporating the user's stopping criterion. -C -C The routine VE08SD has the following argument list: -C -C SUBROUTINE VE08SD(EPSIL,NIT,NGR,DNGR,ITMAX,NFMAX,GXN,NEGCUR, -C DIFEST,X,FX,GX,DX,DF,N,INFO,IFLAG) -C -C EPSIL is a DOUBLE PRECISION variable, containing a measure -C of the accuracy that is required by the user to terminate the -C minimization method successfully. It is the same value as -C that passed to VE08 by the user (see above). It must not be -C altered by VE08SD. -C -C NIT is an INTEGER variable, that contains the number of -C iterations already performed by VE08. It must not be altered -C by VE08SD. -C -C NGR is an INTEGER variable, that contains the number of -C calls already made to the routine ELFNCT. It must not be -C altered by VE08SD. -C -C DNGR is a DOUBLE PRECISION variable, containing the value -C of NGR divided by NS, i.e. a number representing the total -C number of calls to the complete objective function. It must not -C be altered by VE08SD. -C -C ITMAX is an INTEGER variable, containing the maximum number -C of iterations of VE08 allowed. It must not be altered by -C VE08SD. -C -C NFMAX is an INTEGER variable, containing the maximum number -C of calls to the complete objective function allowed. It -C must not be altered by VE08SD. -C -C GXN is a DOUBLE PRECISION variable, containing the -C Euclidean norm of the orthonormal projection of the -C gradient onto the feasible region. It must not be altered by -C VE08SD. -C -C NEGCUR is a LOGICAL variable, whose value is true if and -C only if some direction of negative curvature was detected -C on the current overall quadratic approximation to the -C objective function. It must not be altered by VE08SD. -C -C DIFEST is a LOGICAL variable, whose value is true if and -C only if some of the element gradients are computed by -C differences. It must not be altered by VE08SD. -C -C X is a DOUBLE PRECISION array of length N, containing the -C current best point found by the routine VE08. It must not be -C altered by VE08SD. -C -C FX is a DOUBLE PRECISION variable, containing the value of the -C objective function at X. It must not be altered by VE08SD. -C -C GX is a DOUBLE PRECISION array of length N, containing -C the overall gradient of the overall objective function at X. -C It must not be altered by VE08SD. -C -C DX is a DOUBLE PRECISION variable, containing the -C Euclidean length of the last step taken by the routine VE08. -C It must not be altered by VE08SD. -C -C DF is a DOUBLE PRECISION variable, containing the last improvement -C in objective function values realized by VE08. It must -C not be altered by VE08SD. -C -C N is an INTEGER variable, that contains the number of -C variables of the problem. It must not be altered by VE08SD. -C -C INFO is an INTEGER variable, that is meaningless on input. -C The output value is only relevant if VE08 is terminated by a -C nonzero output value of IFLAG. In this case, the argument -C INFO of VE08 will return this value to the user. -C -C IFLAG is an INTEGER variable, that is meaningless on input. -C If if is decided to continue the minimization process, the -C value 0 should be returned in IFLAG. If it is decided that -C the minimization is complete, a value in the range 1 .LE. IFLAG -C .LE. 10 should be returned. A return value outside the range -C 0 .LE. IFLAG .LE. 10 will cause abnormal termination of VE08 -C with error message 31. -C -C The default routine supplied with VE08 tests -C 1. if the Euclidean length of the gradient is smaller than -C EPSIL when DIFEST is false, or is the maximum of 10**(-4) -C and EPSIL when DIFEST is true, and if the -C quadratic model looks convex (normal successful exit) -C (IFLAG set to 1), -C 2. or if the maximum number of function calls has -C been exceeded (IFLAG set to 2), -C 3. or if both the length of the last step and the last -C improvement in the function values are insignificant -C (IFLAG set to 3). -C -C If none of these occur, then IFLAG is set to 0. -C -C -C Linesearch Step -C --------------- -C -C Another interesting feature of VE08 is to allow the user -C with some knowledge of his problem to verify if the proposed -C value of the step along the current search direction is -C plausible, and to modify it, if it is not the case. This -C verification can be based on some relevant quantities that are -C passed to the user by VE08. -C -C This verification takes the form of a routine called -C VE08TD, whose argument list is as follows: -C -C SUBROUTINE VE08TD(A,X,FX,GXS,SNORM,S,N) -C -C A is a DOUBLE PRECISION variable, containing the steplength -C proposed by VE08. On output (i.e. when control is returned to -C VE08), it may be set to a new and more realistic stepsize -C for minimizing the objective function starting from X in -C the direction S, i.e. A should be such that F( X + A *S ) -C is close to its minimum with respect to A. -C -C X is a DOUBLE PRECISION array of length N containing the base -C point for the current linesearch. It must not be altered by -C VE08TD. -C -C FX is a DOUBLE PRECISION variable, containing the value -C of the objective function at X. It must not be altered by -C VE08TD. -C -C GXS is a DOUBLE PRECISION variable, containing the -C directional derivative at the base point in X along the -C direction S. It must not be altered by VE08TD. -C -C SNORM is a DOUBLE PRECISION variable, containing the -C Euclidean length of the search direction S. It must not be -C altered by VE08TD. -C -C S is a DOUBLE PRECISION array of length N containing -C the current search direction. It must not be altered by -C VE08TD. -C -C N is an INTEGER variable, containing the number of -C variables of the problem. It must not be altered by VE08TD. -C -C The default routine VE08TD supplied with VE08 is just the empty routine: -C -C SUBROUTINE VE08TD(A,X,FX,GSX,SNORM,S,N) -C DOUBLE PRECISION A,X(1),FX,GSX,SNORM,S(1) -C RETURN -C END -C -C It should be replaced by a more complex routine VE08TD only -C if the user feels his knowledge of the problem is sufficient to -C predict rather accurately the minimum of the objective -C function along any search direction. -C -C -C ********************************************************************* -C * * -C * GENERAL * -C * * -C ********************************************************************* -C -C COMMON: None. -C WORKSPACE: Provided by the user (see argument WK). -C USES: VE08BD, VE08CD, VE08DD, ... VE08TD; ELFNCT, RANGE, -C XUPPER, XLOWER. -C I/O: No input; output on device number IPDEVC. -C RESTRICTIONS: N > 0, NS > 0. -C PORTABILITY: Fortran 77. -C -C METHOD -C ------ -C -C The method that is implemented by routine VE08 is a -C partitioned quasi-Newton algorithm, as decribed in (Griewank -C & Toint 1982a), and whose convergence analysis for the convex -C case can be found in (Griewank & Toint 1982b). The theory -C for the non-convex decomposition is still an open question. -C -C The main idea is that a collection of small matrices -C approximating the Hessian matrices of each f is used and -C k -C updated at every iteration using the BFGS or rank 1 updates, -C depending on the curvature properties of f . -C k -C The step is determined by a trust-region approach that -C uses a truncated conjugate-gradient algorithm, followed by a -C very weak linesearch. The trust-region size is then -C augmented if the reduction obtained in the objective function -C is reasonable when compared with the predicted reduction, and -C reduced otherwise. -C -C The strategy for treating bound constraints is based on the -C usual projection device. For a more detailed description, see -C (Bertsekas 1982). -C -C ACKNOWLEDGEMENTS -C ---------------- -C -C The author of the routine wishes to acknowledge the useful -C discussions with A. Griewank, and the most kind help of J.K. Reid -C and I.S. Duff in setting up this specification sheet and -C modifying the routines for the Harwell Subroutine Library format. -C -C REFERENCES -C ---------- -C -C Bertsekas,D.P. 1982. -C Projected Newton Methods for Optimization Problems with -C Simple Constraints. -C SIAM Journal of Control and Optimization 20(2):221-246. -C -C Griewank, A. and Ph.L. Toint. 1982a. -C Partitioned Variable Metric Updates for Large Structured -C Optimization Problems. -C Numerische Mathematik (39):119-137. -C -C Griewank, A. and Ph.L. Toint. 1982b. -C Local Convergence Analysis for Partitioned Quasi-Newton Updates. -C Numerische Mathematik (39):429-448. -C -C Griewank, A. and Ph.L. Toint. 1982c. -C On the Unconstrained Optimization of Partially Separable Functions. -C In M.J.D. Powell (editor), Nonlinear Optimization 1981. -C Academic Press, New-York. -C -C Stewart, G.W. 1967. -C A Modification of Davidon's Minimization Method to Accept -C Difference Approximations of Derivatives. -C Journal of the ACM 14. -C -C ********************************************************************* -C * * -C * EXAMPLE * -C * * -C ********************************************************************* -C -C We now consider the small example problem, -C -C -C --------------------- --------------------- -C / 2 2 / 2 2 -C minimize V 1 + x + ( x - x ) + V 1 + x + ( x - x ) -C 1 2 3 2 3 4 -C -C subject to the bound x .LE. -1. -C 1 -C -C The solution to this problem is at the point (-1,0,0,0) -C and the optimal function value is {1+ SQRT {2}}. -C -C This function has four variables and two elements, namely -C -C ----------------------- -C / 2 2 -C V 1 + x + ( x - x ) -C 1 2 3 -C -C and -C -C ----------------------- -C / 2 2 -C V 1 + x + ( x - x ) -C 2 3 4 -C -C We can then set N=4 and NS=2. The first element -C involves variables 1, 2 and 3 while the second element -C involves variables 2, 3 and 4; we now build the -C vector INVAR corresponding to the problem : -C -C I 1 2 3 4 5 6 -C INVAR(I) 1 2 3 2 3 4 -C <-- elt 1--><--elt 2--> -C -C In this vector, we now locate the position of the first -C variable of each element, and build the vector NVAR as -C follows: -C -C I 1 2 3 -C NVAR(I) 1 4 7 -C -C Observe that NVAR(NS+1) is set to the total length of -C INVAR plus 1, as required. -C -C The first variable is bounded, so we set ISTATE(1) to -C 1, while the other variables are unbounded, which imposes -C ISTATE(I)=-1 for I=2,4. Moreover, the analytical -C gradients of both element functions are available, so that -C the last 2 components of ISTATE are set to 1. We therefore -C obtain -C -C I 1 2 3 4 5 6 -C ISTATE(I) 1 -1 -1 -1 1 1 -C -C The minimum length required for the workspace WK, i.e. -C LWK is given by -C -C LWK = 2+6+4*4+8+3*3+6 = 47 -C -C This count is obviously much more favourable when the -C dimension increases and when the element Hessians can be -C stored in compact form. -C -C We now give the routine ELFNCT that would describe our -C simple example problem: -C -C SUBROUTINE ELFNCT(K,X,FX,GX,NDIMK,NS,JFLAG,FMAX,FNOISE) -C DOUBLE PRECISION X(NDIMK),FX,GX(NDIMK),FMAX,FNOISE,TEMP -C TEMP=X(2)-X(3) -C FX=DSQRT(1.0D0+X(1)**2+TEMP**2) -C FNOISE=1.0D-25*FX -C IF(JFLAG.LT.2)RETURN -C GX(1)=X(1)/FX -C GX(2)=TEMP/FX -C GX(3)=-GX(2) -C GMAX=DMAX1(DABS(GX(1)),DABS(GX(2))) -C FNOISE=DMAX1(FNOISE,1.0D-25*GMAX) -C RETURN -C END -C -C We now wish to write down the routine RANGE associated with -C our example. First, we observe that both element functions -C may be reduced to dependence on two variables with the matrix -C -C ( 1 0 0 ) -C U = ( ) -C k ( 0 1 -1 ) -C -C As a consequence, it is easy to verify that the -C following routine RANGE satisfies all the above requirements -C -C SUBROUTINE RANGE(K,MODE,W1,W2,NDIMK,NSUBK,NS) -C DOUBLE PRECISION W1(1),W2(1 -C NSUBK=2 -C GO TO (1,2,3,4),MODE -C C -C C 1) U*w1=w2 -C C -C 1 W2(1)=W1(1) -C W2(2)=W1(2)-W1(3) -C RETURN -C C -C C 2) U*w2=w1 -C C -C 2 W2(1)=W1(1) -C W2(2)=0.5*W1(2) -C W2(3)=-W2(2) -C RETURN -C C -C C 3) U'*w1=w2 -C C -C 3 W2(1)=W1(1) -C W2(2)=W1(2) -C W2(3)=-W1(2) -C RETURN -C C -C C 4) U'*w2=w1 -C C -C 4 W2(1)=W1(1) -C W2(2)=W1(2) -C RETURN -C END -C -C As one can see, writing routine RANGE is not difficult. -C -C The last thing we have to set up for our example is the bound -C specification. Only the first variable is bounded, so we do -C not have to distinguish between them for finding the bounds. -C Moreover, it is only bounded above, so we have to build an -C artificial lower bound. This gives the following functions : -C -C DOUBLE PRECISION FUNCTION XLOWER(I) -C XLOWER=-1.0D35 -C RETURN -C END -C -C and -C -C DOUBLE PRECISION FUNCTION XUPPER(I) -C XUPPER=-1.0D0 -C RETURN -C END -C -C This completes the supply of information to VE08. We now -C consider a restart of VE08, after some previous computation. -C Assume that after solving the problem we wish to solve another -C which is the same except that the first element function is -C scaled by 0.9. This will require the stored Hessian -C approximations to be rescaled. The components of WK are -C -C -C WK(1) f -C 1 -C WK(2) f -C 2 -C WK(3) to WK(5) g -C 1 -C WK(6) to WK(8) g -C 2 -C WK(9) to WK(11) C (The Hessians are approximated by -C 1 T -C WK(12) to WK(14) C U C U , see above ) -C 2 k k k -C -C The sequence of instructions preceeding the restart call to -C VE08 is now as follows : -C -C C RESCALE THE APPROXIMATE HESSIANS -C C -C DO 10 I=9,11 -C WK(I)=WK(I)*0.9 -C 10 CONTINUE -C C -C C RESET THE PARAMETERS OF VE08AD -C C -C RESTRT=.TRUE. -C C -C C RESTART CALL TO VE08AD -C C -C ... -C -C -C**************************************************************************** -C**************************************************************************** -C**************************************************************************** -C -C DATA TYPE DECLARATIONS -C ---------------------- -C - DOUBLE PRECISION A, ABMAX, ABMIN, ALMAX, ALMIN, B, BBOUND, BETA, - * BND, BOLD, C, CC, CURV, DBIG, DD, DECR, DECR0, DECR1, DECR2, - * DELDIV, DELTAM, DF, DIFGRD, DII, DNGR, DSMALL, DTOL, DX, DXN, - * DXN1, DXN1P, DXN2, EBOUND, EPSACT, EPSIL, EPSIZE, - * EPSQRT, EPSTEP, FA, FACT, FATMAX, FATMIN, FB, FBOLD, FIMAX, - * FLOWBD, FMAX, FN, FNIMAX, FNOISE, FX, FXI, FXOLD, GAM, GAMHAT, - * GATMAX, GATMIN, GMAXI, GNOIS, GNOISE, - * GSA, GSB, GSBOLD, GSINI, GSOUTI, GS0, GTOL, GXMAX, GXNB, - * PRERDI, PRERED, RAPGS, RAPRED, RELPR, RESN2, RESOL2, RS2OLD, - * STEP, STMAX, TEMP, TEMP2, TEST, TOLCNT, UBOUND, XD, XLOWER, - * XNEW, XNOISE, XSIZE, XU, XUPPER, XX, WKIFXI -C - LOGICAL NEGCUR, FKNOWN, BOUNDS, BEND, TESTGX, LTEMP, - * DIFEST, RESTRT, GRADAV, HESDIF, BACTIV -C -C DIMENSION STATEMENTS -C -------------------- -C - DOUBLE PRECISION STEPL(2), X(1), WK(1) - INTEGER INVAR(1), NVAR(1), ISTATE(1), LWK(2), NGR(2), NIT(2) -C - EXTERNAL ELFNCT,RANGE,XUPPER,XLOWER -C -C SOME IMPLEMENTATION REMARKS -C --------------------------- -C -C THE ROUTINE ASSUMES THAT THE RESULT OF AN UNDERFLOWED -C OPERATION IS SET TO ZERO AND THE CALCULATION CONTINUED. -C -C THE CONSTANT DBIG IS SET TO A VERY LARGE NUMBER STILL -C REPRESENTABLE IN MACHINE WITHOUT OVERFLOW. IT IS USED IN THE -C PROCEDURE AS THE VALUE OF +INFINITY. -C -C THE FOLLOWING VALUE IS ADEQUATE FOR DEC2060 -C DIGITAL VAX -C (IN DOUBLE PRECISION ) -C - DATA DBIG /1.0D+36/ -C -C -C*********************************************************************** -C -C INITIALIZATION. -C -C*********************************************************************** -C -C -C SET ITERATION COUNT TO ZERO -C --------------------------- -C - ITMAX = NIT(1) - NIT(2) = -1 -C -C SET FUNCTION CALLS COUNTS TO ZERO -C --------------------------------- -C - NGRMAX = NGR(1) - NGR(2) = 0 -C -C TEST POSITIVITY OF THE DIMENSION (N) -C ------------------------------------ -C - IF (N.GT.0) GO TO 10 - INFO = N - IFLAG = 11 - GO TO 1090 -C -C SET SOME MACHINE DEPENDENT CONSTANTS FOR STEWART'S ALGORITHM FOR -C ---------------------------------------------------------------- -C ESTIMATING THE GRADIENT BY DIFFERENCES -C -------------------------------------- -C -C 1) MACHINE ROUNDING ERROR UNIT -C - 10 CONTINUE - IF (RELPR.GT.0.0D0) GO TO 30 - RELPR = 0.125D0 - 20 CONTINUE - RELPR = 0.5D0*RELPR - TEST = 1.0D0 + RELPR - IF (TEST.GT.1.0D0) GO TO 20 - RELPR = RELPR + RELPR - 30 CONTINUE - IF (RELPR.LT.1.0D0) GO TO 40 - INFO = IDINT(RELPR) - IFLAG = 12 - GO TO 1090 - 40 CONTINUE - EPSQRT = DSQRT(RELPR) - DSMALL = 1000.0D0/(RELPR*DBIG) -C -C 2) NOISE ON THE NORM OF THE CURRENT POINT -C - XNOISE = 2.0*FLOAT(N)*RELPR -C -C 3) TOLERANCE ON THE ROUNDING ERROR ABOVE WHICH ONE SWITCHES TO -C CENTRAL DIFFERENCES -C - TOLCNT = 2.0D-2 -C -C MAXIMUM NUMBER OF ITERATIONS IN STEP STRATEGY. -C ---------------------------------------------- -C - MAXSTE = N + 2 -C -C UNIFORM BOUND ON THE STEP LENGTH AND FIRST STEP SIZE -C ---------------------------------------------------- -C - DELTAM = STMAX - IF (DELTAM.LE.0.0) DELTAM = DSQRT(DBIG) - STEPL(2) = STEPL(1) - IF (STEPL(1).LE.0.0) STEPL(2) = 0.1*DELTAM -C -C INITIALIZE THE TOLERANCE FOR A DIAGONAL ELEMENT TO BE CONSIDERED -C ---------------------------------------------------------------- -C POSITIVE ENOUGH FOR USE IN THE PRECONDITIONING -C ---------------------------------------------- -C - DTOL = 1.0D-10 -C -C SQUARE OF THE ACCURACY ON THE WEIGHTED RESIDUAL THAT IS REQUIRED FOR -C -------------------------------------------------------------------- -C COMPLETION OF A FULL QUASI-NEWTON STEP (LINEAR SYSTEM) -C ------------------------------------------------------ -C - EPSTEP = 0.01D0 -C -C INITIALIZE THE STEPSIZE FOR THE DIFFERENCE ESTIMATION OF THE FIRST -C ------------------------------------------------------------------ -C GRADIENT -C -------- -C - IF (DIFGRD.LE.0.0) DIFGRD = EPSQRT -C -C INITIALIZE THE ACCURACY ABOVE WHICH A POSSIBLE ERROR IN THE -C ----------------------------------------------------------- -C COMPUTATION OF THE ANALYTICAL GRADIENT IS DETECTED -C -------------------------------------------------- -C - GTOL = 100.0*DIFGRD -C -C INITIALIZE THE "DANGEROUS INDICES" TOLERANCE FOR THE BOUNDS -C ----------------------------------------------------------- -C - IF (EBOUND.LE.0.0) EBOUND = 1.0D-5 -C -C SET THE STEP STRATEGY COUNT TO ZERO -C ----------------------------------- -C - NSTEPS = 0 -C -C SET THE QUANTITIES THAT WILL BE USED IN TESTING FOR DIVERGENCE -C -------------------------------------------------------------- -C -C 1) LOWER BOUND THAT DEFINES A VERY LARGE STEP -C - DELDIV = 0.999*DELTAM -C -C 2) NUMBER OF VERY LARGE STEPS ALLOWED BEFORE DIVERGENCE IS -C IDENTIFIED -C - IDIVX = 5 -C -C 3) SET THE CURRENT NUMBER OF VERY LARGE STEPS TO ZERO -C - IDIV = 0 -C -C MAXIMUM NUMBER OF FUNCTION CALLS IN ONE LINESEARCH. -C --------------------------------------------------- -C - MAXLS = 15 -C -C MAXIMUM NUMBER OF FREE TRIALS FOR THE LINE SEARCH IN THE PRESENCE OF -C -------------------------------------------------------------------- -C BOUNDS -C ------ -C - NRLS = 5 -C -C SET THE ANGULAR DECREASE THAT SHOULD BE ACHEIVED IN THE LINE SEARCH -C ------------------------------------------------------------------- -C - DECR = 0.005 -C -C SET THE NEGATIVE CURVATURE INDICATOR TO FALSE -C --------------------------------------------- -C - NEGCUR = .FALSE. -C -C AVOID TESTING THE ANALYTICAL GRADIENT AND ESTIMATING THE INITIAL -C ---------------------------------------------------------------- -C HESSIANS BY DIFFERENCES IF THIS CALL IS A RESTART -C ------------------------------------------------- -C - TESTGX = TESTGX .AND. .NOT.RESTRT - HESDIF = HESDIF .AND. .NOT.RESTRT -C -C PRESET THE INDICES FOR THE USE OF THE USER'S ROUTINE RANGE -C ---------------------------------------------------------- -C - IUB = 1 -C -C SET THE ERROR FLAGS TO THE "EVERYTHING OK" POSITION -C --------------------------------------------------- -C - IFLAG = 0 - INFO = 0 -C -C -C*********************************************************************** -C -C CHECK THE DIMENSIONS -C -C*********************************************************************** -C -C -C COMPUTE THE TOTAL LENGTH OF STORAGE NEEDED FOR THE ELEMENT HESSIANS -C ------------------------------------------------------------------- -C (MB) AND THE MAXIMAL SUBDIMENSION (NDIMAX). VERIFY ALSO -C ------------------------------------------------------- -C POSITIVITY OF NUMBER OF ELEMENTS (NS), OF THE INPUT PARAMETERS -C -------------------------------------------------------------- -C NVAR AND THE COHERENCE OF THE SUBDIMENSIONS. -C -------------------------------------------- -C - IF (NS.GT.0) GO TO 50 - INFO = NS - IFLAG = 14 - GO TO 1090 - 50 CONTINUE - NDIMAX = 0 - M = NVAR(NS+1) - 1 - LB = M + NS + 1 - MB = 0 - DO 90 I=1,NS - ISUBF = I - IF (NVAR(I).GT.0) GO TO 60 - INFO = I - IFLAG = 15 - GO TO 1090 - 60 CONTINUE - NDIMI = NVAR(I+1) - NVAR(I) - IF (NDIMI.GT.0) GO TO 70 - INFO = I - IFLAG = 16 - GO TO 1090 - 70 CONTINUE - CALL RANGE(ISUBF, IUB, X, WK(LB), NDIMI, NSUBI, NS) - IF (NSUBI.LE.NDIMI .AND. NSUBI.GE.0) GO TO 80 - INFO = I - IFLAG = 17 - GO TO 1090 - 80 CONTINUE - MB = MB + ((NSUBI+1)*NSUBI)/2 - NDIMAX = MAX0(NDIMAX,NDIMI) - 90 CONTINUE - INFO = 0 -C -C CHECK THE CONTENT OF THE VECTOR OF SPECIFICATION ISTATE -C ------------------------------------------------------- -C - DIFEST = .FALSE. - ITEMP = N + NS - DO 120 I=1,ITEMP - IF (RESTRT) GO TO 110 - IF (I.LE.N .OR. ISTATE(I).NE.0) GO TO 100 - INFO = I - IFLAG = 19 - GO TO 1090 - 100 CONTINUE - IF (IABS(ISTATE(I)).LE.1) GO TO 110 - INFO = I - IFLAG = 20 - GO TO 1090 -C -C IS DIFFERENCE ESTIMATION FOR THE GRADIENT NEEDED? -C ------------------------------------------------- -C - 110 CONTINUE - IF (I.GT.N .AND. ISTATE(I).LT.0) DIFEST = .TRUE. - 120 CONTINUE -C -C CHECK THE CONTENT OF THE INPUT VECTOR INVAR, AS WELL AS THE TOTAL -C ----------------------------------------------------------------- -C DIMENSION -C --------- -C - IF (M.GT.0) GO TO 130 - INFO = M - IFLAG = 21 - GO TO 1090 - 130 CONTINUE - DO 140 I=1,M - IF (INVAR(I).GT.0 .AND. INVAR(I).LE.N) GO TO 140 - INFO = I - IFLAG = 22 - GO TO 1090 - 140 CONTINUE -C -C DEFINE THE POINTERS TO THE VARIOUS VECTOR SLICES IN WK. -C ------------------------------------------------------- -C -C 1) VALUES OF THE ELEMENT FUNCTIONS -C - NFXI = 0 - LFXI = 1 -C -C 2) SUCCESSIVE ELEMENT GRADIENTS -C - NGXI2 = NFXI + NS -C -C 3) SUCCESSIVE ELEMENT HESSIAN APPROXIMATIONS -C - NB = NGXI2 + M -C -C 4) STEP IN THE VARIABLES -C - NDX = NB + MB - LDX = NDX + 1 -C -C 5) ASSEMBLED DIAGONAL OF THE COMPLETE HESSIAN -C - NDG = NDX + N - LDG = NDG + 1 -C -C 6) ASSEMBLED GRADIENT VECTOR -C - NGX = NDG + N - LGX = NGX + 1 -C -C 7) THREE WORK SPACE VECTORS -C - NW1 = NGX + N - LW1 = NW1 + 1 -C - NW2 = NW1 + N - LW2 = NW2 + 1 -C - NW3 = NW2 + N - LW3 = NW3 + 1 -C -C 8) STORAGE FOR THE OLD ELEMENT GRADIENTS -C (OVERLAPPING WITH THE 2D AND 3RD WORK VECTORS) -C - NGXI1 = NW2 -C -C 9) THREE WORKING VECTORS IN THE MAXIMUM DIMENSION -C OF THE ELEMENT HESSIANS -C - NP1 = MAX0(NW3+N,NGXI1+M) - LP1 = NP1 + 1 -C - NP2 = NP1 + NDIMAX - LP2 = NP2 + 1 -C - NP3 = NP2 + NDIMAX - LP3 = NP3 + 1 -C -C CHECK THE TOTAL LENGTH OF THE WORK SPACE WK. -C -------------------------------------------- -C - LENGTH = NP3 + NDIMAX - IF (LWK(1).GE.LENGTH) GO TO 150 - INFO = LENGTH - IFLAG = 23 - GO TO 1090 - 150 CONTINUE - LWK(2) = MB -C -C CHECK THE AMOUNT OF PRINTOUT PARAMETER IPWHAT -C --------------------------------------------- -C - IF (IPWHAT.GE.0 .AND. IPWHAT.LE.5) GO TO 160 - IFLAG = 34 - INFO = IPWHAT - GO TO 1090 -C -C -C*********************************************************************** -C -C ASSURE FEASIBILITY OF THE STARTING POINT -C -C*********************************************************************** -C -C - 160 CONTINUE - BOUNDS = .FALSE. - GXNB = 0.0 -C -C PROJECT THE STARTING POINT ON THE FEASIBLE DOMAIN: -C -------------------------------------------------- -C -C 1) LOOP ON THE VARIABLES -C - DO 200 I=1,N - NDXPI = NDX + I - WK(NDXPI) = X(I) - IF (ISTATE(I).LE.0) GO TO 200 -C -C 2) THE I-TH VARIABLE IS BOUNDED -C - IVAR = I - UBOUND = XUPPER(IVAR) - BBOUND = XLOWER(IVAR) -C -C 3) CHECK IF THE BOUNDS ARE COMPATIBLE -C - IF (BBOUND.LE.UBOUND) GO TO 170 - INFO = IVAR - IFLAG = 24 - GO TO 1090 -C -C 4) CHECK IF THE BOUNDS ARE EQUAL. -C IN THIS CASE, FIX THE VARIABLE -C - 170 CONTINUE - IF (BBOUND.LT.UBOUND) GO TO 180 - X(I) = UBOUND - FKNOWN = .FALSE. - ISTATE(IVAR) = 0 - GO TO 200 -C -C 5) CHECK IF THE UPPER BOUND IS SATISFIED -C - 180 CONTINUE - BOUNDS = .TRUE. - IF (X(I).LE.UBOUND) GO TO 190 - X(I) = UBOUND - FKNOWN = .FALSE. - ISTATE(I) = 1 - GO TO 200 -C -C 6) CHECK IF THE LOWER BOUND IS SATISFIED -C - 190 CONTINUE - IF (X(I).GE.BBOUND) GO TO 200 - X(I) = BBOUND - FKNOWN = .FALSE. - ISTATE(I) = 1 -C -C 7) END OF THE LOOP ON THE VARIABLES -C - 200 CONTINUE -C -C -C*********************************************************************** -C -C FIND THE ELEMENT FUNCTIONS WHOSE VARIABLES HAVE BEEN MODIFIED -C -C*********************************************************************** -C -C -C LOOP ON THE ELEMENTS -C -------------------- -C - FKNOWN = FKNOWN .AND. (.NOT.TESTGX) - IF (FKNOWN .AND. .NOT.HESDIF) GO TO 240 - DO 210 I=1,NS - IELF = I - IN = N + I - JL = NVAR(I) - NDIMI = NVAR(I+1) - JL - IF (.NOT.FKNOWN) ISTATE(IN) = ISTATE(IN) + 10 - IF (FKNOWN) CALL VE08CD(IELF, NDIMI, WK(LDX), X, INVAR(JL), - * ISTATE(IN)) - 210 CONTINUE -C -C -C*********************************************************************** -C -C RECOMPUTE THE INITIAL FUNCTION AND GRADIENT VALUES -C IF NOT AVAILABLE ON ENTRY, OR IF INITIAL POINT -C WAS INFEASIBLE. -C -C*********************************************************************** -C - LTEMP = .TRUE. - IF (FKNOWN) FXOLD = FX -C -C LOOP ON THE ELEMENT FUNCTIONS -C ----------------------------- -C - LL = LB - FX = 0.0D0 - DO 220 I=1,NS - JL = NVAR(I) - NDIMI = NVAR(I+1) - JL - IELF = I - IN = N + I - LIG2 = NGXI2 + JL -C -C RECOMPUTE THE ELEMENT FUNCTION, AND POSSIBLY TEST ITS GRADIENT -C -------------------------------------------------------------- - ISTATN = ISTATE(IN) - IFXI = NFXI + I - FXI = WK(IFXI) - CALL VE08DD(IELF, NDIMI, NS, WK(LP1), WK(LP2), WK(LP3), - * WK(LIG2), WK(LDG), X, WK(LL), INVAR(JL), ISTATE(1), - * ISTATN, IPDEVC, IPFREQ, FXI, NGRI, NSUBI, LTEMP, TESTGX, - * HESDIF, RELPR, EPSQRT, DIFGRD, FNOISE, DBIG, GTOL, INFO, - * IFLAG, IRTN, ELFNCT, RANGE) - ISTATE(IN) = ISTATN - IF (IRTN.EQ.1) GO TO 1090 - FX = FX + FXI - WK(IFXI) = FXI - NGR(2) = NGR(2) + NGRI - LL = LL + ((NSUBI+1)*NSUBI)/2 - 220 CONTINUE - -C -C CHECK THE PLAUSIBILITY OF THE LOWER BOUND ON THE FUNCTION -C --------------------------------------------------------- -C - IF (FX.GE.FLOWBD) GO TO 230 - IFLAG = 33 - INFO = 0 - GO TO 1090 - 230 CONTINUE -C -C CHECK IF THE ANALYTICAL GRADIENT PASSED THE TEST SUCCESSFULLY -C ------------------------------------------------------------- -C - IF (LTEMP) GO TO 240 - IFLAG = 32 - INFO = IELF - GO TO 1090 -C -C -C*********************************************************************** -C -C TEST FOR STOPPING AFTER FIRST FUNCTION CALL. -C -C*********************************************************************** -C -C -C ASSEMBLE THE OVERAL INITIAL GRADIENT -C ------------------------------------ -C - 240 CONTINUE - DO 250 I=LGX,NW1 - WK(I) = 0.0 - 250 CONTINUE - DO 260 I=1,NS - IELF = I - JL = NVAR(I) - LIG2 = NGXI2 + JL - NDIMI = NVAR(I+1) - JL - CALL VE08FD(IELF, NDIMI, INVAR(JL), WK(LGX), WK(LIG2), ISTATE) - 260 CONTINUE -C -C DETERMINE THE STATE OF THE VARIABLES W.R.T. THEIR BOUNDS -C -------------------------------------------------------- -C - IF (BOUNDS) CALL VE08RD(X, WK(LGX), ISTATE, N, EBOUND, - * XUPPER, XLOWER) -C -C COMPUTE THE NORM OF THE ACTIVE GRADIENT -C --------------------------------------- -C - GXNB = 0.0 - DO 270 I=1,N - IF (ISTATE(I).GE.0 .AND. ISTATE(I).LE.1) GO TO 270 - NGXPI = NGX + I - GXNB = GXNB + WK(NGXPI)**2 - 270 CONTINUE - GXNB = DSQRT(GXNB) -C -C TEST FOR STOPPING -C ----------------- -C - DNGR = FLOAT(NGR(2))/FLOAT(NS) - DX = DBIG - DF = 0.0 - IF (FKNOWN) DF = FXOLD - FX - CALL VE08SD(EPSIL, NIT, NGR, DNGR, ITMAX, NGRMAX, GXNB, NEGCUR, - * DIFEST, X, FX, WK(LGX), DX, DF, N, INFO, IFLAG) - IF (IFLAG.GE.0 .AND. IFLAG.LE.10) GO TO 280 - INFO = IFLAG - IFLAG = 31 - GO TO 1090 -C -C -C*********************************************************************** -C -C WRITE INITIAL DATA -C -C*********************************************************************** -C -C -C - 280 CONTINUE - IF (IPFREQ.GE.0) CALL VE08QD(NIT, FX, GXNB, NGR, DNGR, 1.0D0, - * NSTEPS, BOUNDS, X, WK(LGX), WK(LDX), WK(LB), IPWHAT, IPDEVC, N, - * NS, M, NVAR, WK(LP1), WK(LP2), RANGE) - IF (IFLAG.NE.0) GO TO 1090 -C -C -C*********************************************************************** -C -C COMPUTE THE INITIAL APPROXIMATIONS TO THE ELEMENT HESSIANS -C -C*********************************************************************** -C -C -C AVOID REINITIALIZATION OF THE ELEMENT HESSIANS IF A RESTART OCCURED -C ------------------------------------------------------------------- -C - IF (RESTRT) GO TO 300 -C -C LOOP ON THE ELEMENTS -C -------------------- -C - LL = LB - DO 290 I=1,NS - IELF = I - IN = N + I - NDIMI = NVAR(I+1) - NVAR(I) - CALL VE08ED(IELF, NDIMI, NS, WK(LP1), WK(LP2), WK(LL), - * ISTATE(IN), NSUBI, HESDIF, RANGE) - LL = LL + ((NSUBI+1)*NSUBI)/2 - 290 CONTINUE -C -C -C*********************************************************************** -C -C ASSEMBLE THE DIAGONAL OF THE OVERAL APPROXIMATION B IN WK(LDG) -C -C*********************************************************************** -C -C -C INITIALIZE THE DIAGONAL TO ZERO -C ------------------------------- -C - 300 CONTINUE - DO 310 I=LDG,NGX - WK(I) = 0.0 - 310 CONTINUE - IF (MB.EQ.0) GO TO 330 -C -C ACCUMULATE THE DIAGONAL BY CONSIDERING ELEMENT AFTER ELEMENT -C ------------------------------------------------------------ -C - L = LB - DO 320 I=1,NS - JL = NVAR(I) - NDIMI = NVAR(I+1) - JL - IELF = I - CALL VE08GD(IELF, NDIMI, NS, WK(LP1), WK(LP2), WK(LP3), - * WK(LDG), WK(L), INVAR(JL), NSUBI, RESTRT, HESDIF, RANGE) - L = L + ((NSUBI+1)*NSUBI)/2 - 320 CONTINUE -C -C -C*********************************************************************** -C -C ELIMINATE THE VARIABLES THAT APPEAR ONLY LINEARLY IN THE -C OBJECTIVE FUNCTION BY SETTING THEM TO THEIR BOUNDS -C -C*********************************************************************** -C -C - 330 CONTINUE - NIT(2) = 0 - XSIZE = DSMALL - DX = 0.0 -C -C LOOP ON THE VARIABLES -C --------------------- -C - DO 350 I=1,N - IDX = NDX + I - TEMP = X(I) - WK(IDX) = TEMP -C -C COMPUTE THE SIZE OF THE CURRENT POINT -C ------------------------------------- -C - XSIZE = DMAX1(XSIZE,DABS(TEMP)) -C -C IS THE I-TH VARIABLE A LINEAR VARIABLE? IF THIS IS NOT THE CASE, -C ---------------------------------------------------------------- -C CONTINUE THE LOOP ON THE VARIABLES -C ---------------------------------- -C - NDGPI = NDG + I - IF (WK(NDGPI).NE.0.0) GO TO 350 -C -C YES: CHECK IF IT IS BOUNDED. IF IT IS NOT, ERROR RETURN -C ------------------------------------------------------- -C - IF (ISTATE(I).GE.0) GO TO 340 - INFO = I - IFLAG = 26 - GO TO 1090 -C -C IF IT IS FIXED, CONTINUE THE LOOP ON THE VARIABLES -C -------------------------------------------------- -C - 340 CONTINUE - IF (ISTATE(I).EQ.0) GO TO 350 -C -C ELSE, SET IT TO ITS LOWER OR UPPER BOUND, DEPENDING ON THE SIGN OF -C ------------------------------------------------------------------ -C THE I-TH COMPONENT OF THE OVERAL GRADIENT VECTOR. THEN, CONSIDER -C ---------------------------------------------------------------- -C IT AS A ORDINARY FIXED VARIABLE -C ------------------------------- -C - IGX = NGX + I - IVAR = I - FKNOWN = .FALSE. - IF (WK(IGX).LT.0.0) X(I) = XUPPER(IVAR) - IF (WK(IGX).GT.0.0) X(I) = XLOWER(IVAR) - DX = DX + (X(I)-WK(IDX))**2 - ISTATE(I) = 0 - 350 CONTINUE -C -C COMPUTE THE LOWER BOUND ON THE NORM OF THE NEXT STEP -C ---------------------------------------------------- -C - EPSIZE = RELPR*XSIZE -C -C DETERMINE IF THE APPROXIMATE SOLUTION HAS BEEN CHANGED, DUE TO THE -C ------------------------------------------------------------------ -C MODIFICATION OF PURELY LINEAR VARIABLES -C --------------------------------------- -C -C 1) LOOP ON THE ELEMENTS -C - IF ((.NOT.BOUNDS) .OR. DX.LE.RELPR) GO TO 420 - DO 360 I=1,NS - JL = NVAR(I) - NDIMI = NVAR(I+1) - JL - IN = N + I - IELF = I - CALL VE08CD(IELF, NDIMI, WK(LDX), X, INVAR(JL), ISTATE(IN)) -C -C 5) END OF THE LOOP ON THE ELEMENTS -C - 360 CONTINUE -C -C RECOMPUTE THE OBJECTIVE FUNCTION -C -------------------------------- -C - FXOLD = FX - FMAX = FX -C -C 1) LOOP ON THE ELEMENT FUNCTIONS -C - FX = 0.0D0 - DO 370 I=1,NS - IFXI = NFXI + I - IN = N + I - JL = NVAR(I) - NDIMI = NVAR(I+1) - JL - IELF = I - LIG2 = NGXI2 + JL - ISTATN = ISTATE(IN) - WKIFXI = WK(IFXI) - CALL VE08HD(IELF, NDIMI, NS, WK(LP1), WK(LP2), WK(LP3), X, - * WK(LIG2), INVAR(JL), ISTATE, ISTATN, WKIFXI, NGRI, - * DIFGRD, FNOISE, DBIG, FMAX, INFO, IFLAG, IRTN, ELFNCT, RANGE) - ISTATE(IN) = ISTATN - WK(IFXI) = WKIFXI - IF (IRTN.EQ.1) GO TO 1090 - FX = FX + WK(IFXI) - NGR(2) = NGR(2) + NGRI - 370 CONTINUE -C -C -C*********************************************************************** -C -C TEST AGAIN FOR STOPPING, AFTER MODIFICATION OF THE LINEAR VARIABLES -C -C*********************************************************************** -C -C -C ASSEMBLE THE OVERAL GRADIENT -C ---------------------------- -C - DO 380 I=LGX,NW1 - WK(I) = 0.0 - 380 CONTINUE - DO 390 I=1,NS - IELF = I - JL = NVAR(I) - NDIMI = NVAR(I+1) - JL - LIG2 = NGXI2 + JL - CALL VE08FD(IELF, NDIMI, INVAR(JL), WK(LGX), WK(LIG2), ISTATE) - 390 CONTINUE -C -C DETERMINE THE STATE OF THE VARIABLES W.R.T. THEIR BOUNDS -C -------------------------------------------------------- -C - IF (BOUNDS) CALL VE08RD(X, WK(LGX), ISTATE, N, EBOUND, - * XUPPER, XLOWER) -C -C COMPUTE THE NORM OF THE ACTIVE GRADIENT -C --------------------------------------- -C - GXNB = 0.0 - DO 400 I=1,N - IF (ISTATE(I).GE.0 .AND. ISTATE(I).LE.1) GO TO 400 - NGXPI = NGX + I - GXNB = GXNB + WK(NGXPI)**2 - 400 CONTINUE - GXNB = DSQRT(GXNB) -C -C TEST FOR STOPPING -C ----------------- -C - DNGR = FLOAT(NGR(2))/FLOAT(NS) - DX = DSQRT(DX) - DF = FXOLD - FX - CALL VE08SD(EPSIL, NIT, NGR, DNGR, ITMAX, NGRMAX, GXNB, NEGCUR, - * DIFEST, X, FX, WK(LGX), DX, DF, N, INFO, IFLAG) - IF (IFLAG.GE.0 .AND. IFLAG.LE.10) GO TO 410 - INFO = IFLAG - IFLAG = 31 - GO TO 1090 -C -C -C*********************************************************************** -C -C WRITE THE ACTUAL DATA -C -C*********************************************************************** -C -C -C - 410 CONTINUE - IF (IPFREQ.GE.0) CALL VE08QD(NIT, FX, GXNB, NGR, DNGR, 1.0D0, - * NSTEPS, BOUNDS, X, WK(LGX), WK(LDX), WK(LB), IPWHAT, IPDEVC, N, - * NS, M, NVAR, WK(LP1), WK(LP2), RANGE) - IF (IFLAG.NE.0) GO TO 1090 -C -C -C*********************************************************************** -C -C MAIN LOOP OF THE ITERATIVE MINIMIZATION PROCEDURE -C -C*********************************************************************** -C -C - 420 CONTINUE - RAPRED = 0.0D0 - DXN = 1.0D0 - RESN2 = DXN - DO 1080 NITI=1,ITMAX - RS2OLD = DSQRT(RESN2)/DXN - NIT(2) = NITI -C -C -C*********************************************************************** -C -C STEP SELECTION : INEXACT TRUST REGION STRATEGY WITH PRECONDITIONING -C ( USING TRUNCATED CONJUGATE GRADIENTS ) -C -C*********************************************************************** -C -C -C INITIALIZE THE INITIAL RESIDUAL AS THE OPPOSITE OF THE GRADIENT -C --------------------------------------------------------------- -C AND THE STEP TO ZERO. SOLVE ALSO FOR THE COMPONENTS -C --------------------------------------------------- -C CORRESPONDING TO FIXED VARIABLES OR VARIABLES THAT ARE IN -C --------------------------------------------------------- -C THE DANGEROUS SET. -C ------------------ -C - NEGCUR = .FALSE. - RESN2 = 0.0 - GXMAX = 0.0 - DO 440 J=1,N - IGX = NGX + J - IDX = NDX + J - IW1 = NW1 + J - IW3 = NW3 + J - IDG = NDG + J - WK(IDX) = 0.0 - GXMAX = DMAX1(GXMAX,DABS(WK(IGX))) - IF (ISTATE(J).LT.0 .OR. ISTATE(J).GE.3) GO TO 430 - IF (ISTATE(J).NE.0) WK(IDX) = -WK(IGX)/DABS(WK(IDG)) - WK(IGX) = 0.0 - WK(IW1) = 0.0 - WK(IW3) = 0.0 - GO TO 440 - 430 CONTINUE - WK(IGX) = -WK(IGX) -C -C COMPUTE THE INVERSE OF THE PRECONDITIONING DIAGONAL. -C ---------------------------------------------------- -C - WK(IW3) = 1.0 - DII = DABS(WK(IDG)) - IF (DII.GE.DTOL) WK(IW3) = 1.0/DII -C -C PRECONDITION THE INITIAL RESIDUAL TO OBTAIN THE FIRST SEARCH -C ------------------------------------------------------------ -C DIRECTION -C --------- -C - WK(IW1) = WK(IW3)*WK(IGX) -C -C COMPUTE THE SQUARE OF THE NORM OF THE PRECONDITIONED RESIDUAL -C ------------------------------------------------------------- -C - RESN2 = RESN2 + WK(IW1)*WK(IGX) - 440 CONTINUE -C -C DEFINE THE ACCURACY THAT WILL BE REQUIRED ON THE PRECONDITIONED -C --------------------------------------------------------------- -C RESIDUAL TO STOP THE CONJUGATE GRADIENT ROUTINE -C ----------------------------------------------- -C - EPSACT = DMIN1(EPSTEP,EPSTEP*DSQRT(GXNB))*RESN2 - TEMP = DABS(RAPRED-1.0D0) - TEMP = 0.1*DMAX1(TEMP,RS2OLD) - IF (TEMP.LE.EPSQRT) EPSACT = DMIN1(EPSACT,TEMP) -C -C MAIN LOOP OF THE CONJUGATE GRADIENT ROUTINE THAT DEFINES -C -------------------------------------------------------- -C THE STEP STRATEGY -C ----------------- -C - DO 480 KSTEP=1,MAXSTE - KSTEPS = KSTEP -C -C COMPUTE THE PRODUCT OF THE CURRENT APPROXIMATION B TIMES THE SEARCH -C ------------------------------------------------------------------- -C DIRECTION WK(LW1) AND STORE THE RESULT IN WK(LW2) -C ------------------------------------------------- -C -C 1) SET THE VECTOR WK(LW2) TO ZERO -C - DO 450 K=LW2,NW3 - WK(K) = 0.0 - 450 CONTINUE -C -C 2) LOOP ON THE ELEMENTS OF B -C - L = LB - DO 460 I=1,NS - JL = NVAR(I) - NDIMI = NVAR(I+1) - JL - IELF = I - CALL VE08ID(IELF, NDIMI, NS, WK(LP1), WK(LP2), WK(LW1), - * WK(LW2), WK(L), INVAR(JL), ISTATE, NSUBI, RANGE) - L = L + ((NSUBI+1)*NSUBI)/2 - 460 CONTINUE -C -C COMPUTE THE CURVATURE TERM -C -------------------------- -C - CALL VE08MD(CURV,WK(LW1),WK(LW2),N) -C -C TEST THE SIGN OF THE CURVATURE. IF IT IS NEGATIVE, SET NEGCUR TO -C ---------------------------------------------------------------- -C .TRUE. AND BRANCH TO CHECK THE SIZE OF THE STEP W.R.T. TRUST -C ------------------------------------------------------------ -C REGION BOUDARY -C -------------- -C - IF (CURV.GT.0.0) GO TO 470 - NEGCUR = .TRUE. - GO TO 490 -C -C POSITIVE CURVATURE: COMPUTE THE OPTIMAL STEPSIZE -C ------------------------------------------------ -C - 470 CONTINUE - STEP = RESN2/CURV -C -C COMPUTE THE NEW POINT -C --------------------- -C - RESOL2 = RESN2 - CALL VE08ND(WK(LDX), WK(LW1), STEP, WK(LGX), WK(LW2), - * WK(LW3), RESN2, N) -C -C TEST FOR STOPPING THE CONJUGATE GRADIENT BECAUSE THE REQUIRED -C ------------------------------------------------------------- -C ACCURACY ON THE PRECONDITIONED RESIDUAL HAS BEEN REACHED -C -------------------------------------------------------- -C - IF (RESN2.LE.EPSACT) GO TO 490 -C -C COMPUTE THE NEXT SEARCH DIRECTION OF THE CONJUGATE GRADIENT -C ----------------------------------------------------------- -C - BETA = RESN2/RESOL2 - CALL VE08OD(WK(LW1), WK(LW2), WK(LW1), BETA, N) -C -C END OF THE MAIN LOOP OF THE CONJUGATE GRADIENT -C ---------------------------------------------- -C - 480 CONTINUE -C -C CHECK FEASIBILITY OF THE STEP W.R.T. THE TRUST REGION BOUNDARY -C -------------------------------------------------------------- -C -C 1) CHECK THE LENGTH OF THE PRESENT DIRECTION -C - 490 CONTINUE - CALL VE08MD(XX,WK(LDX),WK(LDX),N) - DXN = DSQRT(XX) - -C 2) THE NORM OF THE DIRECTION IS TOO SMALL : EXIT FROM VE08AD -C WITH AN ERROR MESSAGE -C - IF (DXN.GE.EPSIZE) GO TO 500 - IF (CURV.LT.0.0 .AND. STEPL(2).GE.EPSIZE) GO TO 500 - INFO = NIT(2) - IFLAG = 35 - GO TO 1090 -C -C 3) THE DIRECTION IS TOO SHORT, BUT THE MODEL IS CONVEX : -C NO RESCALING IS NECESSARY -C - 500 CONTINUE - IF (DXN.LE.STEPL(2) .AND. (.NOT.NEGCUR)) GO TO 530 -C -C 4) IF THE DIRECTION IS TOO LONG, SHORTEN IT TO THE BOUNDARY -C - IF (DXN.LE.STEPL(2)) GO TO 520 - FACT = STEPL(2)/DXN - DO 510 I=1,N - IDX = NDX + I - WK(IDX) = WK(IDX)*FACT -510 CONTINUE - DXN = STEPL(2) - GO TO 530 -C -C 5) COMPUTE THE STEP SIZE TO THE BOUNDARY -C - 520 CONTINUE - IF (DXN.EQ.STEPL(2)) GO TO 530 - CALL VE08MD(DD,WK(LW1),WK(LW1),N) - IF (DD.LT.RELPR) GO TO 530 - CALL VE08MD(XD,WK(LDX),WK(LW1),N) - XD = XD/DD - DD = ((STEPL(2)-DXN)*(STEPL(2)+DXN))/DD - XU = XD + DSQRT(XD*XD+DD) - STEP = DD/XU -C -C 6) SET THE NEW DIRECTION -C - CALL VE08OD(WK(LDX), WK(LDX), WK(LW1), STEP, N) - DXN = STEPL(2) -C -C UPDATE THE COUNT OF ITERATIONS IN THE STEP STRATEGY -C --------------------------------------------------- -C - 530 CONTINUE - NSTEPS = NSTEPS + KSTEPS -C -C -C*********************************************************************** -C -C LINESEARCH -C -C*********************************************************************** -C -C -C MODIFY THE SEARCH DIRECTION WHEN MINIMIZING IN A SUBSPACE -C --------------------------------------------------------- -C - BACTIV = .FALSE. - IF (.NOT.BOUNDS) GO TO 560 -C -C 1) LOOP ON THE VARIABLES -C - DO 550 I=1,N - IDX = NDX + I - NGXPI = NGX + I - WK(NGXPI) = WK(IDX) -C -C 2) IF THE VARIABLE IS NOT BOUNDED, CONSIDER NEXT ONE -C - IF (ISTATE(I).NE.1) GO TO 550 - BACTIV = .TRUE. -C -C 3) SET THE CORRESPONDING SEARCH DIRECTION COMPONENT TO ZERO -C - IVAR = I - UBOUND = XUPPER(IVAR) - IF (X(I).NE.UBOUND) GO TO 540 - WK(IDX) = DMIN1(WK(IDX),0.0D0) - GO TO 550 - 540 CONTINUE - WK(IDX) = DMAX1(WK(IDX),0.0D0) -C -C 4) END OF THE LOOP ON THE VARIABLES -C - 550 CONTINUE -C -C COMPUTE THE INITIAL DIRECTIONAL DERIVATIVE -C ------------------------------------------ -C - 560 CONTINUE - GSB = 0.0 - GSOUTI = 0.0 -C -C 1) LOOP ON THE ELEMENTS -C - DO 580 I=1,NS - JL = NVAR(I) - JU = NVAR(I+1) - 1 -C -C 2) ACCUMULATE THE I-TH DIRECTIONAL DERIVATIVE -C - DO 570 J=JL,JU - IVAR = INVAR(J) - ISVAR = ISTATE(IVAR) - IF (ISVAR.EQ.0) GO TO 570 - NGXI2J = NGXI2 + J - NDXIVA = NDX + IVAR - GSB = GSB + WK(NGXI2J)*WK(NDXIVA) - IF (.NOT.BOUNDS) GO TO 570 - NGXIVA = NGX + IVAR - IF (ISVAR.EQ.1 .OR. ISVAR.EQ.2) GO TO 570 - GSOUTI = GSOUTI + WK(NGXI2J)*WK(NGXIVA) - 570 CONTINUE -C -C 3) END OF THE LOOP ON THE ELEMENTS -C - 580 CONTINUE -C -C TEST IF INITIAL SLOPE IS NEGATIVE. IF IT IS NOT THE CASE, STOP -C -------------------------------------------------------------- -C WITH AN ERROR MESSAGE -C --------------------- -C - IF (GSB.LT.0.0) GO TO 590 - INFO = 0 - IFLAG = 27 - GO TO 1090 -C -C COMPUTE THE LENGTH OF THE SEARCH DIRECTION, DUE TO THE -C ------------------------------------------------------ -C POSSIBLE MODIFICATION INTRODUCED BY THE PROJECTION -C -------------------------------------------------- -C ON THE ACTIVE BOUNDS -C -------------------- -C - 590 CONTINUE - DXN1 = DXN - DO 600 I=1,N - NGXPI = NGX + I - WK(NGXPI) = X(I) - 600 CONTINUE - IF (.NOT.BACTIV) GO TO 620 - DXN2 = 0.0 - DO 610 I=1,N - NDXPI = NDX + I - DXN2 = DXN2 + WK(NDXPI)**2 - 610 CONTINUE - DXN1 = DSQRT(DXN2) -C -C CHECK IF THE DIRECTION IS NOT TOO SMALL. IT IS IS -C ------------------------------------------------- -C BELOW EPSIZE IN NORM, EXIT FROM VE08AD WITH AN ERROR -C ---------------------------------------------------- -C MESSAGE -C ------- -C - 620 CONTINUE - IF (DXN1.GE.EPSIZE) GO TO 630 - INFO = NIT(2) - IFLAG = 35 - GO TO 1090 -C -C PRESET THE SLOPE OF THE LINEAR UPPER BOUND ON THE -C ------------------------------------------------- -C FUNCTION VALUES -C --------------- - 630 CONTINUE - DECR0 = DECR*GSOUTI - DECR1 = DECR*GSB - DECR2 = DMIN1(DECR0,DECR1) -C -C INITIALIZE SOME LINESEARCH VARIABLES -C ------------------------------------ -C -C 1) LINESEARCH PARAMETER AT ORIGIN -C - B = 0.0 -C -C 2) FUNCTION VALUE AND DIRECTIONAL DERIVATIVE AT ORIGIN -C - FB = FX - GS0 = GSB - BOLD = 0.0 - FBOLD = FX - GSBOLD = GS0 -C -C 3) MAXIMUM AND MINIMUM LINESEARCH PARAMETER -C - ALMAX = DELTAM/DXN1 - TEMP = FLOWBD-FX - IF (DABS(DECR2).GT.1.D0) THEN - ALMAX = DMIN1(TEMP/DECR2,ALMAX) - ELSE IF (DABS(TEMP).LE.DABS(DECR2)*DBIG) THEN - ALMAX = DMIN1(TEMP/DECR2,ALMAX) - END IF - ALMIN = 0.0 -C -C 4) PRESET THE INDICATOR THAT SAYS IF A BOUND HAS BEEN ENCOUNTERED -C DURING THE LINESEARCH -C - BEND = .FALSE. -C -C SET THE FIRST ESTIMATE OF THE FIRST VALUE OF THE LINESEARCH -C ----------------------------------------------------------- -C PARAMETER -C --------- -C - A = 1.0D0 -C -C BRANCHING POINT FOR THE LINESEARCH MAIN ITERATION -C ------------------------------------------------- -C - DO 900 IKLS=1,MAXLS - KLS = IKLS -C -C LET THE USER VERIFY THE PLAUSIBILITY OF THE STEP LENGTH -C ------------------------------------------------------- -C - IF (.NOT.BEND) CALL VE08TD(A, X, FX, GS0, DXN1, WK(LDX), N) -C -C TEST IF A IS SUFFICIENTLY INSIDE THE PERMITTED INTERVAL -C ------------------------------------------------------- -C - TEMP = DMAX1(1.01*ALMIN,RELPR) - IF (A.LT.TEMP .OR. A.GT.0.99*ALMAX) A = 0.5*(ALMIN+ALMAX) - IF (A.GT.0.4*DELTAM .AND. ALMIN.GT.0.) A = DMIN1(A,10.*ALMIN) -C -C TEST FOR LINESEARCH FAILURE -C --------------------------- -C - TEMP = DMIN1(ALMAX-ALMIN,DABS(A-B),A) - IF (TEMP*DXN1.GT.EPSIZE) GO TO 640 - INFO = KLS - IFLAG = 28 - GO TO 1090 -C -C COMPUTE THE MAXIMUM ADMISSIBLE FUNCTION VALUE -C --------------------------------------------- -C - 640 CONTINUE - FMAX = FX + DECR1*A - IF ((BACTIV .AND. IKLS.GE.NRLS) .OR. BEND) FMAX = FX -C -C BUILD THE NEW POINT WHERE THE FUNCTION SHOULD BE EVALUATED -C ---------------------------------------------------------- -C (IN WK(LGX)) -C ------------ -C -C 1) BUILD THE NEW POINT -C - DO 670 I=1,N - IGX = NGX + I - NW1PI = NW1 + I - WK(NW1PI) = WK(IGX) - XNEW = X(I) - IF (ISTATE(I).EQ.0) GO TO 660 - NDXPI = NDX + I - XNEW = X(I) + A*WK(NDXPI) - IF (ISTATE(I).LT.0) GO TO 660 -C -C 2) PROJECT ON THE FEASIBLE REGION -C - IVAR = I - UBOUND = XUPPER(IVAR) - IF (XNEW.LE.UBOUND) GO TO 650 - XNEW = UBOUND - BEND = .TRUE. - GO TO 660 - 650 CONTINUE - BBOUND = XLOWER(IVAR) - IF (XNEW.GE.BBOUND) GO TO 660 - XNEW = BBOUND - BEND = .TRUE. -C -C 3) BUILD THE NEW POINT IN WK(LGX) -C - 660 CONTINUE - WK(IGX) = XNEW - 670 CONTINUE - BACTIV = BACTIV .OR. BEND -C -C FIND THE ELEMENT FUNCTIONS WHOSE VARIABLES HAVE BEEN MODIFIED -C ------------------------------------------------------------- -C - DO 680 I=1,NS - IELF = I - IN = N + I - JL = NVAR(I) - NDIMI = NVAR(I+1) - JL - CALL VE08CD(IELF, NDIMI, WK(LW1), WK(LGX), INVAR(JL), - * ISTATE(IN)) - 680 CONTINUE -C -C EVALUATE THE OBJECTIVE FUNCTION AT THE NEW POINT WK(LGX) -C -------------------------------------------------------- -C -C -C 1) INITIALIZE THE OVERAL FUNCTION VALUE AND DIRECTIONAL DERIVATIVE -C - FN = 0.0 - FA = 0.0 - GSA = 0.0 - GSINI = 0.0 - FIMAX = -DBIG - FNIMAX = RELPR -C -C 2) LOOP ON THE ELEMENT FUNCTIONS -C - DO 690 I=1,NS - JL = NVAR(I) - NDIMI = NVAR(I+1) - JL - LIG1 = NGXI1 + JL - LIG2 = NGXI2 + JL - IELF = I - IFXI = NFXI + I - IN = N + I - ISTATN = ISTATE(IN) - WKIFXI = WK(IFXI) - CALL VE08JD(IELF, NDIMI, NS, WK(LP1), WK(LGX), X, WK(LIG1), - * WK(LIG2), INVAR(JL), ISTATE, ISTATN, NGRI, DIFEST, - * GSINI, A, FNOISE, FMAX, WKIFXI, GSA, KLS, MAXLS, INFO, - * IFLAG, IRTN, ELFNCT) - ISTATE(IN) = ISTATN - WK(IFXI) = WKIFXI - IF (IRTN.EQ.1) GO TO 1090 - IF (IRTN.EQ.2) GO TO 790 - FIMAX = DMAX1(DABS(WK(IFXI)),FIMAX) - FNIMAX = DMAX1(FNIMAX,FNOISE) - FA = FA + WK(IFXI) - FN = FN + FNOISE - NGR(2) = NGR(2) + NGRI - 690 CONTINUE -C -C TEST SATISFACTION OF THE LINE SEARCH CRITERIA -C --------------------------------------------- -C -C 1) IF BOUNDS ARE PRESENT, DECIDE UPON THE MAXIMUM ADMISSIBLE -C FUNCTION VALUE -C - LTEMP = BEND .OR. (BACTIV .AND. IKLS.GE.NRLS) - IF (.NOT.BACTIV) GO TO 700 - FMAX = FX + DECR0*A + DECR*GSINI - IF (IKLS.GE.NRLS) GO TO 700 - TEMP = FX + DECR1*A - FMAX = DMIN1(FMAX,TEMP) -C -C 2) CHECK IF THE GRADIENT IS AVAILABLE -C - 700 CONTINUE - GRADAV = .NOT.DIFEST - IF (GRADAV) GO TO 720 -C -C 3) CRITERION WITHOUT GRADIENT -C - IF (FA.LE.FMAX) GO TO 730 - ITEMP1 = N + 1 - ITEMP2 = N + NS - DO 710 IN=ITEMP1,ITEMP2 - IF (ISTATE(IN).GE.8) ISTATE(IN) = ISTATE(IN) - 10 - 710 CONTINUE - GO TO 750 -C -C 4) CRITERION WITH GRADIENT -C - 720 CONTINUE - RAPGS = GSA/GS0 - IF (FA.LE.FMAX .AND. (LTEMP .OR. RAPGS.LE.0.9)) GO TO 910 - GO TO 750 -C -C ESTIMATE THE FINAL GRADIENT BY DIFFERENCES -C ------------------------------------------ -C - 730 CONTINUE -C -C 1) LOOP ON THE ELEMENT FUNCTIONS -C - GSA = 0.0 - GSINI = 0.0 - L = LB - DO 740 I=1,NS - JL = NVAR(I) - NDIMI = NVAR(I+1) - JL - LIG1 = NGXI1 + JL - LIG2 = NGXI2 + JL - IELF = I - IN = N + I - IFXI = NFXI + I - ISTATN = ISTATE(IN) - CALL VE08KD(IELF, NDIMI, NS, WK(LP1), WK(LP2), WK(LP3), - * WK(LIG1), WK(LIG2), WK(LGX), X, WK(L), ISTATN, - * INVAR(JL), ISTATE, NSUBI, NGRI, WK(IFXI), FNOISE, XNOISE, - * EPSQRT, RELPR, TOLCNT, A, DXN1, GSA, GSINI, GXMAX, INFO, - * IFLAG, IRTN, ELFNCT, RANGE) - ISTATE(IN) = ISTATN - IF (IRTN.EQ.1) GO TO 1090 - NGR(2) = NGR(2) + NGRI - IF (NSUBI.EQ.0) GO TO 740 - L = L + ((NSUBI+1)*NSUBI)/2 - 740 CONTINUE -C -C 2) FINAL CURVATURE TEST -C - RAPGS = GSA/GS0 - IF (RAPGS.LE.0.9 .OR. LTEMP) GO TO 910 -C -C 4) TEST NOT SATISFIED: BRANCH FOR A NEW LINESEARCH STEP -C - GSB = GSBOLD - B = BOLD - GSBOLD = GSA - BOLD = A - FBOLD = FA - GRADAV = .TRUE. -C -C TEST TO CHECK THAT THE DIFFERENCES IN FUNCTION VALUES ARE ABOVE THE -C ------------------------------------------------------------------- -C NOISE -C ----- -C - 750 CONTINUE - IF (IKLS.EQ.1 .OR. (BEND .AND. IKLS.LE.NRLS)) GO TO 760 - TEMP = DABS(FB-FA) - IF (TEMP.GT.FN .AND. TEMP.GT.FIMAX*FLOAT(NS)*RELPR) GO TO - * 760 - INFO = IKLS - IFLAG = 29 - GO TO 1090 -C -C TEST FOR MAXIMUM NUMBER OF LINESEARCH ITERATIONS -C ------------------------------------------------ -C - 760 CONTINUE - IF (IKLS.GE.MAXLS) GO TO 900 -C -C BERTSEKAS STRATEGY FOR SHORTENING THE STEP IN THE PRESENCE OF -C ------------------------------------------------------------- -C UNDETECTED ACTIVE BOUNDS -C ------------------------ -C - IF (.NOT.BOUNDS .OR. IKLS.NE.NRLS) GO TO 790 -C -C 1) COMPUTE THE MAXIMUM FEASIBLE STEPSIZE -C - GAMHAT = DELTAM - DO 780 K=1,N - IF (ISTATE(K).NE.3) GO TO 780 - IVAR = K - NDXPK = NDX + K - TEMP = WK(NDXPK) - IF (DABS(TEMP).LT.1.D-15) GO TO 780 - IF (TEMP.GT.0.0) BND = XUPPER(IVAR) - IF (TEMP.LT.0.0) BND = XLOWER(IVAR) - TEMP2 = BND - X(K) - IF (DABS(TEMP).GE.1.0) GO TO 770 - IF (TEMP2.GE.TEMP*GAMHAT) GO TO 780 - GAMHAT = TEMP2/TEMP - 770 CONTINUE - GAM = TEMP2/TEMP - GAMHAT = DMIN1(GAM,GAMHAT) - 780 CONTINUE -C -C 2) TEST IT W.R.T. THE BOUNDS ALREADY AVAILABLE -C - IF (GAMHAT.GE.ALMAX .AND. (BEND .OR. DIFEST)) GO TO 810 - IF (GAMHAT.GE.ALMAX) GO TO 790 - C = GAMHAT - ALMAX = 1.02*GAMHAT - ALMIN = 0.0 - BEND = .TRUE. - GO TO 890 -C -C ORDER THE TWO DATA POINTS -C ------------------------- -C - 790 CONTINUE - IF (A.LT.B) GO TO 800 - ABMIN = B - ABMAX = A - FATMIN = FB - FATMAX = FA - IF (.NOT.GRADAV) GO TO 810 - GATMAX = GSA - GATMIN = GSB - GO TO 810 - 800 CONTINUE - ABMIN = A - ABMAX = B - FATMIN = FA - FATMAX = FB - IF (.NOT.GRADAV) GO TO 810 - GATMAX = GSB - GATMIN = GSA -C -C UPDATE THE UPPER AND LOWER BOUNDS ON THE STEP -C --------------------------------------------- -C -C 1) IF THE CALCULATED FUNCTION VALUE EXCEEDS ITS MAXIMUM ALLOWED -C VALUE FMAX, REDUCE THE INTERVAL TO (0,A) -C - 810 CONTINUE - IF (FA.GT.FMAX) ALMAX = DMIN1(A,ALMAX) -C -C 2) IF THE GRADIENT AT A IS UNAVAILABLE, WE KNOW NOTHING MORE. -C HENCE, BRANCH FOR A NEW BISECTION STEP BETWEEN ALMIN -C AND ALMAX -C - IF (.NOT.GRADAV .OR. BEND) GO TO 880 -C -C 3) ELSE, THE GRADIENT IS AVAILABLE AT ABMIN AND ABMAX. -C CHECK THE SLOPE AT ABMIN. -C - LTEMP = .TRUE. - IF (GATMIN.LT.0.0) GO TO 820 -C -C 4) THE SLOPE AT THE FIRST POINT IS POSITIVE : THE MINIMUM -C MUST LIE IN (ALMIN,ABMIN) -C - ALMAX = DMIN1(ALMAX,ABMIN) - GO TO 870 -C -C 5) THE SLOPE AT THE FIRST POINT IS NEGATIVE. -C TEST THE SIZE OF THE FUNCTION AT THE FIRST POINT -C - 820 CONTINUE - IF (FATMIN.LE.(FX+DECR1*ABMIN)) GO TO 830 -C -C 6) THE FUNCTION AT THE FIRST POINT IS ABOVE ITS UPPER LIMIT. -C THE STEP MUST LIE BEFORE ABMIN. -C - ALMAX = DMIN1(ALMAX,ABMIN) - GO TO 870 -C -C 7) THE FUNCTION AT THE FIRST POINT IS BELOW ITS UPPER LIMIT. -C ABMIN MUST BE A LOWER BOUND ON THE STEP -C - 830 CONTINUE - ALMIN = DMAX1(ALMIN,ABMIN) -C -C 8) TEST THE SIGN OF THE SLOPE AT THE SECOND POINT -C - IF (GATMAX.LT.0.0) GO TO 840 -C -C 7) THE SLOPE AT THE SECOND POINT IS POSITIVE (AND NEGATIVE -C AT THE FIRST) : THE MINIMUM MUST LIE BEFORE ABMAX -C - ALMAX = DMIN1(ALMAX,ABMAX) - GO TO 870 -C -C 8) THE SLOPE AT THE SECOND POINT IS NEGATIVE. -C TEST THE SIZE OF THE FUNCTION AT THE SECOND POINT. -C - 840 CONTINUE - IF (FATMAX.LT.FATMIN) GO TO 850 -C -C 9) THE FUNCTION VALUE AT THE SECOND POINT IS ABOVE ITS UPPER LIMIT. -C THE MINIMUM MUST LIE BEFORE ABMAX, BUT THE INTERPOLATION WONT -C WORK : BRANCH FOR A BISECTION BETWEEN ALMIN AND ABMAX -C - ALMAX = DMIN1(ALMAX,ABMAX) - GO TO 880 -C -C 10) THE FUNCTION VALUE AT THE SECOND POINT IS BELOW ITS UPPER LIMIT. -C TEST IT AGAIN AGAINST THE MORE SEVERE BOUND -C - 850 CONTINUE - IF (FATMAX.LE.(FX+DECR1*ABMAX)) GO TO 860 -C -C 11) THE FUNCTION AT THE SECOND POINT IS ABOVE THE LINEAR LIMIT, -C BUT BELOW THE FIRST POINT. THE STEP MUST LIE IN -C (ABMIN,ABMAX) -C - ALMAX = DMIN1(ALMAX,ABMAX) - GO TO 870 -C -C 12) THE FUNCTION AT THE SECOND POINT IS SUFFICIENTLY LOW : -C THE MINIMUM MUST BE FURTHER ON THE SEARCH DIRECTION -C - 860 CONTINUE - LTEMP = .FALSE. - ALMIN = ABMAX -C -C USE CUBIC EXTRAPOLATION -C ----------------------- -C - 870 CONTINUE - C = GATMIN + GATMAX - 3.0*(FATMAX-FATMIN)/(ABMAX-ABMIN) - CC = DMAX1(C*C-GATMIN*GATMAX,0.0D0) - CC = DSQRT(CC) - TEMP = GATMAX - GATMIN + CC + CC - IF (TEMP.EQ.0.0) GO TO 880 - C = ABMIN + (ABMAX-ABMIN)*(C-GATMIN+CC)/TEMP -C -C IF THE MINIMUM LIES FURTHER ON THE SEARCH DIRECTION, MAKE SURE THE -C ------------------------------------------------------------------ -C STEP IS LARGE ENOUGH -C -------------------- -C - IF (LTEMP) GO TO 890 - IF (C.LT.1.01*ALMIN) C = ALMIN + ALMIN - GO TO 890 -C -C HALVE THE STEP -C -------------- -C - 880 CONTINUE - C = 0.5D0*(ALMAX+ALMIN) -C -C SWITCH VALUES -C ------------- -C - 890 CONTINUE - B = A - A = C - FB = FA - GSB = GSA -C -C END OF THE MAIN LOOP OF THE LINE SEARCH -C --------------------------------------- - 900 CONTINUE - INFO = MAXLS - IFLAG = 28 - GO TO 1090 -C -C EXIT OF THE LINESEARCH: COMPUTE THE NEW POINT, THE TRUE DISPLACEMENT -C -------------------------------------------------------------------- -C AND CORRECT THE STEP NORMS AND THE DIRECTIONAL DERIVATIVE -C --------------------------------------------------------- -C ALSO COMPUTE THE LOWER BOUND ON THE NEXT STEP -C --------------------------------------------- -C - 910 CONTINUE - DXN2 = 0.0 - XSIZE = DSMALL - DO 920 I=1,N - IGX = NGX + I - IDX = NDX + I - WK(IDX) = WK(IGX) - X(I) - DXN2 = DXN2 + WK(IDX)**2 - TEMP = WK(IGX) - X(I) = TEMP - XSIZE = DMAX1(XSIZE,DABS(TEMP)) - 920 CONTINUE - FXOLD = FX - FX = FA - DXN1P = DSQRT(DXN2) - FACT = DXN1P/DXN1 - DXN1 = DXN1P - DXN = FACT*DXN - EPSIZE = RELPR*XSIZE - GSA = A*GSA - IF (BEND) GO TO 930 - GS0 = A*GS0 - GO TO 960 -C -C RECOMPUTE THE INITIAL DIRECTIONAL DERIVATIVE IF BENDING HAS OCCURED -C ------------------------------------------------------------------- - 930 CONTINUE - GS0 = 0.0 - DO 950 I=1,NS - JL = NVAR(I) - JU = NVAR(I+1) - 1 - DO 940 J=JL,JU - IVAR = INVAR(J) - IF (ISTATE(IVAR).EQ.0) GO TO 940 - NDXVAR=NDX+IVAR - NGXI2J=NGXI2+J - GS0 = GS0 + WK(NGXI2J)*WK(NDXVAR) - 940 CONTINUE - 950 CONTINUE -C -C -C*********************************************************************** -C -C TEST FOR STOPPING THE MAIN MINIMIZATION ITERATION -C -C*********************************************************************** -C -C -C TEST FOR DIVERGENCE -C ------------------- -C - 960 CONTINUE - IF (DXN.LT.DELDIV) GO TO 970 - IDIV = IDIV + 1 - IF (IDIV.LT.IDIVX) GO TO 980 - INFO = NIT(2) - IFLAG = 36 - GO TO 1090 - 970 CONTINUE - IDIV = 0 -C ASSEMBLE THE OVERAL GRADIENT -C ---------------------------- -C - 980 CONTINUE - DO 990 I=LGX,NW1 - WK(I) = 0.0 - 990 CONTINUE - DO 1000 I=1,NS - IELF = I - JL = NVAR(I) - NDIMI = NVAR(I+1) - JL - LIG1 = NGXI1 + JL - CALL VE08FD(IELF, NDIMI, INVAR(JL), WK(LGX), WK(LIG1), ISTATE) - 1000 CONTINUE -C -C DETERMINE THE STATE OF THE VARIABLES W.R.T. THEIR BOUNDS -C -------------------------------------------------------- -C - IF (BOUNDS) CALL VE08RD(X, WK(LGX), ISTATE, N, EBOUND, - * XUPPER, XLOWER) -C -C COMPUTE THE NORM OF THE ACTIVE GRADIENT -C --------------------------------------- -C - GXNB = 0.0 - DO 1010 I=1,N - IF (ISTATE(I).GE.0 .AND. ISTATE(I).LE.1) GO TO 1010 - NGXPI=NGX+I - GXNB = GXNB + WK(NGXPI)**2 - 1010 CONTINUE - GXNB = DSQRT(GXNB) -C -C TEST FOR STOPPING -C ----------------- -C - DNGR = FLOAT(NGR(2))/FLOAT(NS) - DF = FXOLD - FX - CALL VE08SD(EPSIL, NIT, NGR, DNGR, ITMAX, NGRMAX, GXNB, NEGCUR, - * DIFEST, X, FX, WK(LGX), DXN1, DF, N, INFO, IFLAG) - IF (IFLAG.GE.0 .AND. IFLAG.LE.10) GO TO 1020 - INFO = IFLAG - IFLAG = 31 - GO TO 1090 -C -C -C*********************************************************************** -C -C WRITE THE ACTUAL DATA -C -C*********************************************************************** -C -C - 1020 CONTINUE - IF (IPFREQ.LT.0) GO TO 1030 - MIP = 1 - IF (IPFREQ.NE.0) MIP = MOD(NIT(2),IPFREQ) - IF (MIP.EQ.0 .OR. IFLAG.NE.0) CALL VE08QD(NIT, FX, GXNB, NGR, - * DNGR, A, NSTEPS, BOUNDS, X, WK(LGX), WK(LDX), WK(LB), IPWHAT, - * IPDEVC, N, NS, M, NVAR, WK(LP1), WK(LP2), RANGE) - 1030 CONTINUE - IF (IFLAG.NE.0) GO TO 1090 -C -C -C*********************************************************************** -C -C UPDATE THE ELEMENT HESSIANS APPROXIMATIONS -C -C*********************************************************************** -C -C -C INITIALIZE THE PREDICTED FUNCTION REDUCTION AND THE SQUARE OF THE -C ----------------------------------------------------------------- -C NORM OF THE OVERAL RESIDUAL -C --------------------------- -C - PRERED = -GSA - RESN2 = 0.0 - GNOIS = FNIMAX - IF (DIFEST) GNOIS = DSQRT(GNOIS) -C -C LOOP ON THE ELEMENTS -C -------------------- -C - LL = LB - DO 1050 I=1,NS - JL = NVAR(I) - NDIMI = NVAR(I+1) - JL - LIG1 = NGXI1 + JL - NIG1 = LIG1 - 1 - NIG2 = NGXI2 + JL - 1 -C -C BUILD THE DIFFERENCES IN GRADIENT AND SHIFT THE NEW VALUE -C --------------------------------------------------------- -C - GMAXI = 0.0 - DO 1040 J=1,NDIMI - IGXI1 = NIG1 + J - IGXI2 = NIG2 + J - TEMP = WK(IGXI1) - TEMP2 = WK(IGXI2) - WK(IGXI2) = TEMP - WK(IGXI1) = TEMP - TEMP2 - GMAXI = DMAX1(GMAXI,DABS(TEMP)+DABS(TEMP2)) - 1040 CONTINUE - GNOISE = GMAXI*GNOIS -C -C UPDATE -C ------ -C - IELF = I - IN = N + I - LTEMP = NIT(2).GT.1 .OR. RESTRT .OR. (HESDIF .AND. - * ISTATE(IN).GT.0) - CALL VE08LD(IELF, NDIMI, NS, WK(LP1), WK(LP2), WK(LP3), - * WK(LIG1), WK(LDX), WK(LDG), WK(LL), RESN2, ISTATE(IN), - * INVAR(JL), KU, PRERDI, LTEMP, RELPR, GNOISE, INFO, IFLAG, - * IRTN, RANGE) - IF (IRTN.EQ.1) GO TO 1090 - PRERED = PRERED + PRERDI - LL = LL + KU - 1050 CONTINUE -C -C -C*********************************************************************** -C -C TRUST REGION BOUNDS MANAGEMENT -C -C*********************************************************************** -C -C -C COMPUTE THE RATIO OF THE PREDICTED REDUCTION AND THE TRUE REDUCTION -C ------------------------------------------------------------------- -C - IF (DF.LE.RELPR) GO TO 1060 - RAPRED = PRERED/DF -C -C IF THIS RATIO IS TOO SMALL, HALVE THE MAXIMUM STEPLENGTH -C -------------------------------------------------------- -C - IF (RAPRED.LE.10.0) GO TO 1070 - 1060 CONTINUE - STEPL(2) = 0.5*DXN - GO TO 1080 -C -C ELSE, DOUBLE IT -C --------------- -C - 1070 CONTINUE - STEPL(2) = DXN + DXN - STEPL(2) = DMIN1(STEPL(2),DELTAM) -C -C -C*********************************************************************** -C -C END OF THE MAIN ITERATION OF THE WHOLE MINIMIZATION PROCESS -C -C*********************************************************************** -C -C - 1080 CONTINUE - IFLAG = 18 - INFO = NIT(2) -C -C -C*********************************************************************** -C -C EXITS OF THE MAIN ROUTINE AND RETURN TO USER'S CONTROL -C -C*********************************************************************** -C -C -C DETERMINE WHAT TYPE OF EXIT IS REQUIRED -C --------------------------------------- -C - 1090 CONTINUE - IF (IPFREQ.LT.0) RETURN - DNGR = FLOAT(NGR(2))/FLOAT(NS) - IPWERR = MIN0(2,IPWHAT) - IF (IFLAG.GT.10) CALL VE08QD(NIT, FX, GXNB, NGR, DNGR, A, NSTEPS, - * BOUNDS, X, WK(LGX), WK(LDX), WK(LB), IPWERR, IPDEVC, N, NS, M, - * NVAR, WK(LP1), WK(LP2), RANGE) - WRITE (IPDEVC,99999) - IF (IFLAG.GT.10) WRITE (IPDEVC,99998) - GO TO (1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, - * 1100, 1110, 1120, 1130, 1140, 1150, 1160, 1170, 1180, 1190, - * 1200, 1210, 1220, 1230, 1240, 1250, 1260, 1270, 1280, 1290, - * 1300, 1310, 1320, 1330, 1340, 1350, 1360), IFLAG -C -C EXIT BECAUSE USER'S TERMINATION CRITERION (VE08SD) HAS BEEN -C ----------------------------------------------------------- -C SATISFIED -C --------- -C - 1100 CONTINUE - WRITE (IPDEVC,99997) IFLAG - RETURN -C -C ERROR EXITS -C ----------- -C -C 1) IFLAG=11 -C - 1110 CONTINUE - WRITE (IPDEVC,99996) INFO - RETURN -C -C 2) IFLAG=12 -C - 1120 CONTINUE - WRITE (IPDEVC,99995) RELPR - RETURN -C -C 3) IFLAG=13 -C - 1130 CONTINUE - WRITE (IPDEVC,99994) FNOISE - RETURN -C -C 4) IFLAG=14 -C - 1140 CONTINUE - WRITE (IPDEVC,99993) INFO - RETURN -C -C 5) IFLAG=15 -C - 1150 CONTINUE - WRITE (IPDEVC,99992) INFO, INFO, NVAR(INFO) - RETURN -C -C 6) IFLAG=16 -C - 1160 CONTINUE - WRITE (IPDEVC,99991) INFO, NDIMI - RETURN -C -C 7) IFLAG=17 -C - 1170 CONTINUE - WRITE (IPDEVC,99990) NSUBI, INFO - RETURN -C -C 8) IFLAG=18 -C - 1180 CONTINUE - WRITE (IPDEVC,99989) INFO - RETURN -C -C 9) IFLAG=19 -C - 1190 CONTINUE - IELF = INFO - N - WRITE (IPDEVC,99988) IELF, INFO, ISTATE(INFO) - RETURN -C -C 10) IFLAG=20 -C - 1200 CONTINUE - WRITE (IPDEVC,99987) INFO, INFO, ISTATE(INFO) - RETURN -C -C 11) IFLAG=21 -C - 1210 CONTINUE - WRITE (IPDEVC,99986) INFO - RETURN -C -C 12) IFLAG=22 -C - 1220 CONTINUE - WRITE (IPDEVC,99985) INFO, INVAR(INFO) - RETURN -C -C 13) IFLAG=23 -C - 1230 CONTINUE - WRITE (IPDEVC,99984) LWK(1), INFO - RETURN -C -C 14) IFLAG=24 -C - 1240 CONTINUE - WRITE (IPDEVC,99983) INFO, UBOUND, BBOUND - RETURN -C -C 15) IFLAG=25 -C - 1250 CONTINUE - WRITE (IPDEVC,99982) INFO - RETURN -C -C 16) IFLAG=26 -C - 1260 CONTINUE - WRITE (IPDEVC,99981) INFO - RETURN -C -C 17) IFLAG=27 -C - 1270 CONTINUE - WRITE (IPDEVC,99980) - RETURN -C -C 18) IFLAG=28 -C - 1280 CONTINUE - WRITE (IPDEVC,99979) INFO - RETURN -C -C 19) IFLAG=29 -C - 1290 CONTINUE - WRITE (IPDEVC,99978) - RETURN -C -C 20) IFLAG=30 -C - 1300 CONTINUE - WRITE (IPDEVC,99977) INFO - RETURN -C -C 21) IFLAG=31 -C - 1310 CONTINUE - WRITE (IPDEVC,99976) INFO - RETURN -C -C 22) IFLAG=32 -C - 1320 CONTINUE - WRITE (IPDEVC,99975) - RETURN -C -C 23) IFLAG=33 -C - 1330 CONTINUE - WRITE (IPDEVC,99974) FX, FLOWBD - RETURN -C -C 24) IFLAG=34 -C - 1340 CONTINUE - WRITE (IPDEVC,99973) IPWHAT - RETURN -C -C 25) IFLAG=35 -C - 1350 CONTINUE - WRITE (IPDEVC,99972) RELPR - RETURN -C -C 26) IFLAG=36 -C - 1360 CONTINUE - WRITE (IPDEVC,99971) IDIVX, DELDIV - RETURN -99999 FORMAT (/5X, 50H************* EXIT OF ROUTINE VE08AD ************* - * //) -99998 FORMAT (7X, 14HERROR RETURN .) -99997 FORMAT (5X, 46HTERMINATION CRITERION OF USER SATISFIED (FLAG=, - * I2, 1H)) -99996 FORMAT (5X, 24HTHE TOTAL DIMENSION N = , I5, 2H .) -99995 FORMAT (5X, 29HTHE MACHINE PRECISION RELPR (, D14.7, 8H) IS OUT, - * 10H OF RANGE.) -99994 FORMAT (5X, 44HTHE NOISE IN THE OBJECTIVE FUNCTION FNOISE (, - * D14.7, 14H) IS NEGATIVE.) -99993 FORMAT (5X, 36HTHE NUMBER OF ELEMENT FUNCTION NS = , I5, 6H IS NO, - * 11HN POSITIVE.) -99992 FORMAT (5X, 38HTHE STARTING POSITION OF THE VARIABLES, 8H INTERNA, - * 9HL TO THE /8X, I5, 28H-TH ELEMENT FUNCTION (NVAR(, I5, 4H) = , - * I5, 2H) /8X, 16HIS NON POSITIVE.) -99991 FORMAT (5X, 40HTHE NUMBER OF VARIABLES INTERNAL TO THE , I5, - * 21H-TH ELEMENT FUNCTION /8X, 9H(NDIMI = , I5, 15H) IS NON POSITI, - * 3HVE.) -99990 FORMAT (5X, 26HTHE SUBDIMENSION (NSUBI = , I5, 16H) ASSOCIATED TO - * /8X, 4HTHE , I5, 35H-TH ELEMENT FUNCTION BY THE ROUTINE, - * 12H OF THE USER/8X, 20HRANGE IS IMPOSSIBLE.) -99989 FORMAT (5X, 24HALGORITHM STOPPED AFTER , I5, 12H ITERATIONS.) -99988 FORMAT (5X, 26HTHE INITIAL STATUS OF THE , I3, 15H-TH ELEMENT FUN, - * 5HCTION/8X, 8H(ISTATE(, I5, 4H) = , I2, 21H) IS NOT IN THE ALLOW, - * 9HED RANGE.) -99987 FORMAT (5X, 26HTHE INITIAL STATUS OF THE , I5, 12H-TH VARIABLE/8X, - * 8H(ISTATE(, I5, 4H) = , I5, 30H) IS NOT IN THE ALLOWED RANGE.) -99986 FORMAT (5X, 49HTHE TOTAL OF THE ELEMENT DIMENSIONS M=NVAR(NS+1)=, - * I5, 1H.) -99985 FORMAT (5X, 17HTHE INDEX OF THE , I4, 25H-TH VARIABLE DESCRIBED IN - * /8X, 7HINVAR (, I5, 2H)=, I5, 16H) IS IMPOSSIBLE.) -99984 FORMAT (5X, 37HTHE TOTAL AVAILABLE WORKSPACE LWK(1)=, I6, 4H IS , - * 10HTOO SMALL./8X, 22HIT SHOULD BE AT LEAST , I6, 1H.) -99983 FORMAT (5X, 18HTHE BOUNDS ON THE , I5, 12H-TH VARIABLE, 7H ARE IN, - * 11HCONSISTENT./8X, 14H(UPPER BOUND =, D14.7, 16H AND LOWER BOUND, - * 2H =, D14.7, 2H.)) -99982 FORMAT (5X, 50HTHE USER HAS STOPPED THE MINIMIZATION PROCEDURE BY/ - * 8X, 32HSETTING THE SUBFUNCTION FLAG TO , I3, 1H.) -99981 FORMAT (5X, 4HTHE , I5, 29H-TH VARIABLE IS PURELY LINEAR, - * 15H AND UNBOUNDED.) -99980 FORMAT (5X, 50HTHE DIRECTIONAL DERIVATIVE AT THE BEGINNING OF THE/ - * 8X, 27HLINESEARCH IS NON-NEGATIVE.) -99979 FORMAT (5X, 25HLINESEARCH FAILURE AFTER , I2, 8H TRIALS.) -99978 FORMAT (5X, 46HTHE ALGORITHM WAS STOPPED BECAUSE THE FUNCTION, - * 12H DIFFERENCES/8X, 30HALONG THE CURRENT DIRECTION OF, 7H SEARCH, - * 18H ARE INSIGNIFICANT/8X, 15HCOMPARED TO THE, 15H NOISE ON THE F, - * 15HUNCTION VALUES.) -99977 FORMAT (5X, 24HTHE PRODUCT SBS FOR THE , I5, 17H-TH ELEMENT HESSI, - * 2HAN/8X, 29HIS TOO SMALL FOR BFGS UPDATE.) -99976 FORMAT (5X, 51HINCORRECT FLAG IFLAG RETURNED FROM THE LAST CALL TO - * /8X, 50HTHE ROUTINE VE08SD. IT SHOULD BE BETWEEN 0 AND 10,/8X, - * 17HAND WAS EQUAL TO , I3, 2H .) -99975 FORMAT (5X, 51HTHE MINIMIZATION PROCEDURE HAS BEEN STOPPED BECAUSE - * /8X, 53HA POSSIBLE ERROR HAS BEEN DETECTED IN THE COMPUTATION/8X, - * 25HOF THE ELEMENT GRADIENTS.) -99974 FORMAT (5X, 51HINCONSISTENT LOWER BOUND ON THE OBJECTIVE FUNCTION. - * /8X, 44HTHE FUNCTION VALUE AT THE STARTING POINT IS , D14.7/8X, - * 25HWHILE THE LOWER BOUND IS , D14.7, 1H.) -99973 FORMAT (5X, 37HTHE PARAMETER IPWHAT IS OUT OF RANGE./8X, 6HITS VA, - * 7HLUE IS , I3, 30H AND SHOULD BE 0,1,2,3,4 OR 5.) -99972 FORMAT (5X, 38HTHE ROUTINE VE08AD WAS STOPPED BECAUSE/8X, - * 44HTHE NORM OF THE COMPUTED STEP WAS NEGLIGIBLE/8X,10H( LT RELPR, - * 12H*NORM(X) = , D9.2, 2H)./8X, - * 39HTHIS ERROR USUALLY OCCURS WHEN TOO MUCH/ - * 8X, 37HACCURACY IS REQUIRED FOR TERMINATION.) -99971 FORMAT (5X, 38HTHE ROUTINE VE08AD WAS STOPPED BECAUSE/8X, I1, - * 28H SUCCESSIVE VERY LARGE STEPS/8X, 5H( GT , D9.2, 8H ) WERE , - * 19HTAKEN.. THE ROUTINE/8X, 21HIS PROBABLY DIVERGING) - END - - SUBROUTINE VE08BD(FNOISE, IFFLAG, IFLAG, INFO, IRTN) -C -C -C ********************************************************* -C * * -C * CONTROL OF THE ERROR DURING FUNCTION EVALUATION * -C * * -C ********************************************************* -C -C -C - DOUBLE PRECISION FNOISE -C -C -C IS FNOISE REALISTIC? -C -------------------- -C - IRTN = 0 - IF (FNOISE.GE.0.0) GO TO 10 - INFO = IDINT(FNOISE) - IFLAG = 13 - IRTN = 1 - RETURN -C -C DOES THE USER WANT TO STOP THE ROUTINE? -C --------------------------------------- -C - 10 CONTINUE - IF (IFFLAG.GE.0) RETURN - INFO = IFFLAG - IFLAG = 25 - IRTN = 1 - RETURN - END - - SUBROUTINE VE08CD(I, NDIMI, X, XNEW, INVAR, ISELF) -C -C -C ************************************************** -C * * -C * DETECT IF ELEMENT I IS TO BE RECOMPUTED * -C * * -C ************************************************** -C -C - DOUBLE PRECISION X(1), XNEW(1) - INTEGER INVAR(1) -C -C LOOP ON THE VARIABLES OF THE I-TH ELEMENT FUNCTION -C -------------------------------------------------- -C - DO 10 J=1,NDIMI - IVAR = INVAR(J) - IF (X(IVAR).EQ.XNEW(IVAR)) GO TO 10 -C -C THE IVAR-TH VARIABLE HAS BEEN MODIFIED: ASSURE THE I-TH ELEMENT -C --------------------------------------------------------------- -C FUNCTION WILL BE RECOMPUTED -C --------------------------- -C - IF (ISELF.LE.8) ISELF = ISELF + 10 - RETURN -C -C END OF THE LOOP ON THE VARIABLES OF THE I-TH ELEMENT -C ---------------------------------------------------- -C - 10 CONTINUE - RETURN - END - - SUBROUTINE VE08DD(IELF, NDIMI, NS, P1, P2, P3, GX, DG, X, B, - * INVAR, ISTATE, ISELF, IPDEVC, IPFREQ, FXI, NGRI, NSUBI, GXOK, - * TESTGX, HESDIF, RELPR, EPSQRT, DIFGRD, FNOISE, DBIG, GTOL, - * INFO, IFLAG, IRTN, ELFNCT, RANGE) -C -C -C ********************************************************* -C * * -C * EVALUATION OF THE IELF-THE ELEMENT FUNCTION, * -C * ITS GRADIENT AND POSSIBLY ITS HESSIAN AT * -C * THE INITIAL POINT X. * -C * * -C ********************************************************* -C -C - DOUBLE PRECISION DBIG, DIF, DIFGRD, DIF2, RELPR, EPSQRT, ERR, - * FNOISE, FXI, FXIP, GTOL, TEMP - LOGICAL GXOK, TESTGX, HESDIF - EXTERNAL ELFNCT,RANGE -C -C - DOUBLE PRECISION P1(1), P2(1), P3(1), GX(1), X(1), DG(1), B(1) - INTEGER INVAR(1), ISTATE(1) -C -C -C CHECK IF THE I-TH ELEMENT FUNCTION HAS TO BE RECOMPUTED. IF THIS IS -C -------------------------------------------------------------------- -C THE CASE, BRANCH TO POSSIBLE HESSIAN ESTIMATION -C ----------------------------------------------- -C - IRTN = 0 - IUTMB = 4 - IUMB = 2 - NGRI = 0 - IF (ISELF.GE.8) GO TO 10 - GO TO 90 -C -C THE I-TH ELEMENT FUNCTION HAS TO BE RECOMPUTED -C ---------------------------------------------- -C - 10 CONTINUE - ISELF = ISELF - 10 - IFFLAG = 2 - IF (ISELF.LT.0) IFFLAG = 1 -C -C PREPARE THE DATA POINT FOR THE I-TH ELEMENT FUNCTION IN P1 -C ---------------------------------------------------------- -C - DO 20 I=1,NDIMI - INVARI = INVAR(I) - P1(I) = X(INVARI) - 20 CONTINUE -C -C EVALUATE THE I-TH ELEMENT FUNCTION -C ---------------------------------- -C - CALL ELFNCT(IELF, P1, FXI, GX, NDIMI, NS, IFFLAG, DBIG, FNOISE) - NGRI = NGRI + 1 - CALL VE08BD(FNOISE, IFFLAG, IFLAG, INFO, JRTN) - IF (JRTN.EQ.1) GO TO 140 -C -C ESTIMATION OF THE GRADIENT BY DIFFERENCE -C ---------------------------------------- -C -C 1) IS ESTIMATION NEEDED? -C - IF (ISELF.GE.0 .AND. (.NOT.TESTGX)) GO TO 90 -C -C 2) YES: LOOP ON THE SUBDIMENSIONS TO PREPARE THE PERTURBED DATA -C POINT -C - JJ = 1 - DO 80 J=1,NDIMI - INVARJ = INVAR(J) - IF (ISTATE(INVARJ).EQ.0) GO TO 70 - TEMP = P1(J) - IF (DABS(TEMP).GT.RELPR) GO TO 30 - P1(J) = TEMP + EPSQRT - GO TO 40 - 30 CONTINUE - P1(J) = TEMP*(1.0+DIFGRD) - 40 CONTINUE - DIF = P1(J) - TEMP -C -C 3) RE-EVALUATE THE I-TH ELEMENT FUNCTION AT PERTURBED POINT -C - CALL ELFNCT(IELF, P1, FXIP, P2, NDIMI, NS, IFFLAG, DBIG, FNOISE) - NGRI = NGRI + 1 - CALL VE08BD(FNOISE, IFFLAG, IFLAG, INFO, JRTN) - IF (JRTN.EQ.1) GO TO 140 -C -C 4) ESTIMATE THE J-TH COMPONENT OF THE GRADIENT AND TEST IT W.R.T. IT -C ANALYTICAL VALUE IF NEEDED -C - TEMP = (FXIP-FXI)/DIF - IF (ISELF.GE.0) GO TO 50 - GX(JJ) = TEMP - GO TO 60 - 50 CONTINUE - ERR = DABS(GX(JJ)-TEMP)/(1.0+DABS(TEMP)) - IF (ERR.LE.GTOL) GO TO 60 - IF (IPFREQ.GE.0) WRITE (IPDEVC,99999) J, IELF, GX(JJ), TEMP - GXOK = .FALSE. - 60 CONTINUE -C -C 5) RESET DATA POINT -C - P1(J) = P1(J) - DIF -C -C 6) END OF THE LOOP ON THE SUBDIMENSIONS -C - 70 CONTINUE - JJ = JJ + 1 - 80 CONTINUE -C -C ESTIMATE THE INITIAL I-TH ELEMENT HESSIAN BY DIFFERENCES -C -------------------------------------------------------- -C -C 1) GET THE SUBDIMENSION AND PROJECT THE ELEMENT GRADIENT ON THE -C RANGE -C - 90 CONTINUE - IF (.NOT.HESDIF) RETURN - CALL RANGE(IELF, IUTMB, GX, P3, NDIMI, NSUBI, NS) - IF (NSUBI.EQ.0 .OR. ISELF.LT.0) RETURN - L = 0 - DIF2 = DIFGRD + DIFGRD -C -C 2) LOOP ON THE SUBDIMENSION -C - NSUBIT = NSUBI - DO 130 J=1,NSUBI -C -C 3) BUILD THE PERTURBED DATA POINT -C - DO 100 K=1,NSUBI - P1(K) = 0.0 - 100 CONTINUE - P1(J) = 1.0 - CALL RANGE(IELF, IUMB, P1, P2, NDIMI, NSUBIT, NS) - TEMP = DIF2 - DO 110 K=1,NDIMI - INVARK = INVAR(K) - DG(K) = TEMP*P2(K) + X(INVARK) - 110 CONTINUE -C -C 4) COMPUTE THE NEW PERTURBED GRADIENT, PROJECTED ON THE RANGE -C - IFFLAG = 2 - CALL ELFNCT(IELF, DG, TEMP, P1, NDIMI, NS, IFFLAG, DBIG, FNOISE) - NGRI = NGRI + 1 - CALL VE08BD(FNOISE, IFFLAG, IFLAG, INFO, JRTN) - IF (JRTN.EQ.1) GO TO 140 -C -C 5) COMPUTE THE J-TH ROW OF THE I-TH ELEMENT HESSIAN -C - CALL RANGE(IELF, IUTMB, P1, P2, NDIMI, NSUBIT, NS) - DO 120 K=1,J - L = L + 1 - B(L) = (P2(K)-P3(K))/DIF2 - 120 CONTINUE -C -C 6) END OF THE LOOP ON THE SUBDIMENSION -C - 130 CONTINUE - RETURN -C -C ERROR EXIT -C ---------- -C - 140 IRTN = 1 - RETURN -99999 FORMAT (5X, 36HPOSSIBLE ERROR IN THE COMPUTATION OF, 9H THE ANAL, - * 16HYTICAL GRADIENT /5X, 4HTHE , I3, 21H-TH COMPONENT OF THE , - * I5, 24H-TH ELEMENT GRADIENT HAS/5X, 22HA COMPUTED VALUE OF , - * D14.7, 4H AND/5X, 22HAN ESTIMATED VALUE OF , D14.7, 1H./) - END - - SUBROUTINE VE08ED(IELF, NDIMI, NS, P1, P2, B, ISELF, NSUBI, - * HESDIF, RANGE) -C -C -C ************************************************ -C * * -C * INITIALIZES THE IELF-TH ELEMENT HESSIAN * -C * TO THE IDENTITY MATRIX * -C * * -C ************************************************ -C -C - LOGICAL HESDIF - EXTERNAL RANGE -C -C - DOUBLE PRECISION P1(1), P2(1), B(1) -C -C -C GET THE TRUE SUBDIMENSION NSUBI -C ------------------------------- -C - L = 1 - IUB = 1 - CALL RANGE(IELF, IUB, P1, P2, NDIMI, NSUBI, NS) -C -C INITIALIZE THE I-TH ELEMENT HESSIAN -C ----------------------------------- -C - IF (NSUBI.EQ.0 .OR. (HESDIF .AND. ISELF.GT.0)) RETURN - B(1) = 1.0 - IF (NSUBI.LE.1) RETURN - DO 20 J=2,NSUBI - JJ = J - 1 - DO 10 K=1,JJ - L = L + 1 - B(L) = 0.0 - 10 CONTINUE - L = L + 1 - B(L) = 1.0 - 20 CONTINUE - RETURN - END - - SUBROUTINE VE08FD(IELF, NDIMI, INVAR, GX, GI, ISTATE) -C -C -C ************************************************** -C * * -C * THIS ROUTINE ACCUMULATES THE OVERAL GRADIENT * -C * FROM THE IELF-TH ELEMENT GRADIENT * -C * * -C ************************************************** -C -C -C -C - DOUBLE PRECISION GX(1), GI(1) - INTEGER INVAR(1), ISTATE(1) -C -C - DO 10 I=1,NDIMI - IVAR = INVAR(I) - IF (ISTATE(IVAR).NE.0) GX(IVAR) = GX(IVAR) + GI(I) - 10 CONTINUE - RETURN - END - - SUBROUTINE VE08GD(IELF, NDIMI, NS, P1, P2, P3, DG, B, INVAR, - * NSUBI, RESTRT, HESDIF, RANGE) -C -C -C ************************************************************ -C * * -C * ASSEMBLY OF THE DIAGONAL OF THE FIRST ELEMENT HESSIAN * -C * APPROXIMATIONS * -C * * -C ************************************************************ -C -C - LOGICAL RESTRT, HESDIF - EXTERNAL RANGE -C -C - DOUBLE PRECISION P1(1), P2(1), P3(1), DG(1), B(1) - INTEGER INVAR(1) -C -C -C -C LOOP ON THE SUBDIMENSION -C ------------------------ -C - IUB = 1 - IUTB = 3 - DO 40 J=1,NDIMI -C -C PREPARE THE J-TH SUB-BASIS VECTOR -C --------------------------------- -C - DO 10 K=1,NDIMI - P1(K) = 0.0 - 10 CONTINUE - P1(J) = 1.0 -C -C PROJECT ON THE RANGE -C -------------------- -C - CALL RANGE(IELF, IUB, P1, P2, NDIMI, NSUBI, NS) - IF (NSUBI.EQ.0) RETURN -C -C WHEN RESTARTING, COMPUTE THE PRODUCT WITH THE STARTING APPROXIMATION -C -------------------------------------------------------------------- -C TO THE ELEMENT HESSIANS: -C ------------------------ -C - IF (.NOT.(RESTRT .OR. HESDIF)) GO TO 30 - CALL VE08PD(B, P2, P3, NSUBI) - DO 20 K=1,NSUBI - P2(K) = P3(K) - 20 CONTINUE -C -C RESTORE TO FULL SUBDIMENSION -C ---------------------------- -C - 30 CONTINUE - CALL RANGE(IELF, IUTB, P2, P1, NDIMI, NSUBI, NS) -C -C ACCUMULATE IN THE CORRESPONDING DIAGONAL ELEMENT -C ------------------------------------------------ -C - IVAR = INVAR(J) - DG(IVAR) = DG(IVAR) + P1(J) - 40 CONTINUE - RETURN - END - - SUBROUTINE VE08HD(IELF, NDIMI, NS, P1, P2, P3, X, GXI, INVAR, - * ISTATE, ISELF, FXI, NGRI, DIFGRD, FNOISE, DBIG, FMAX, INFO, - * IFLAG, IRTN, ELFNCT, RANGE) -C -C -C ************************************************************* -C * * -C * RECOMPUTATION OF THE IELF-TH ELEMENT FUNCTION AND * -C * GRADIENT AFTER MODIFICATION OF THE LINEAR VARIABLES * -C * * -C ************************************************************* -C -C - DOUBLE PRECISION DBIG, DERJ, DIFGRD, FMAX, FNOISE, FXFWD, FXI - EXTERNAL ELFNCT, RANGE -C -C - DOUBLE PRECISION P1(1), P2(1), P3(1), X(1), GXI(1) - INTEGER INVAR(1), ISTATE(1) -C -C - IRTN = 0 - IUB = 1 - IUTB = 3 - IUMB = 2 - NGRI = 0 -C -C CHECK IF THE I-TH ELEMENT FUNCTION HAS TO BE RECOMPUTED. -C -------------------------------------------------------- -C - IF (ISELF.LT.8) RETURN -C -C THE I-TH ELEMENT FUNCTION HAS TO BE RECOMPUTED -C ---------------------------------------------- -C - ISELF = ISELF - 10 -C -C PREPARE THE DATA POINT FOR THE I-TH ELEMENT FUNCTION IN P1 -C ---------------------------------------------------------- -C - DO 10 I=1,NDIMI - INVARI = INVAR(I) - P1(I) = X(INVARI) - 10 CONTINUE -C -C EVALUATE THE I-TH ELEMENT FUNCTION -C ---------------------------------- -C - IFFLAG = 2 - IF (ISELF.LT.0) IFFLAG = 1 - CALL ELFNCT(IELF, P1, FXI, GXI, NDIMI, NS, IFFLAG, FMAX, FNOISE) - NGRI = NGRI + 1 - CALL VE08BD(FNOISE, IFFLAG, IFLAG, INFO, JRTN) - IF (JRTN.EQ.1) GO TO 60 -C -C ESTIMATE THE ELEMENT GRADIENT BY DIFFERENCES -C -------------------------------------------- -C -C 1) CHECK IF ESTIMATION IS NEEDED -C - IF (ISELF.GT.0) RETURN -C -C 2) GET THE TRUE SUBDIMENSION -C - CALL RANGE(IELF, IUB, GXI, P3, NDIMI, NSUBI, NS) - IF (NSUBI.EQ.0) RETURN -C -C 3) LOOP ON THE SUBDIMENSION -C - NSUBIT = NSUBI - DO 40 J=1,NSUBI - INVARJ = INVAR(J) - IF (NSUBI.EQ.NDIMI .AND. ISTATE(INVARJ).EQ.0) GO TO 40 -C -C 4) SET DATA POINT -C - DO 20 K=1,NSUBI - P2(K) = 0.0 - 20 CONTINUE - P2(J) = 1.0 - CALL RANGE(IELF, IUTB, P2, P1, NDIMI, NSUBIT, NS) - DO 30 K=1,NDIMI - INVARK = INVAR(K) - P2(K) = P1(K)*DIFGRD + X(INVARK) - 30 CONTINUE -C -C 5) COMPUTE FUNCTION VALUE AT PERTURBED DATA POINT -C - CALL ELFNCT(IELF, P2, FXFWD, P1, NDIMI, NS, IFFLAG, DBIG, - * FNOISE) - NGRI = NGRI + 1 - CALL VE08BD(FNOISE, IFFLAG, IFLAG, INFO, JRTN) - IF (JRTN.EQ.1) GO TO 60 -C -C 6) ESTIMATE THE J-TH COMPONENT OF THE ELEMENT GRADIENT IN -C SUBDIMENSION -C - DERJ = (FXFWD-FXI)/DIFGRD - P3(J) = DERJ - P3(J) - 40 CONTINUE -C -C 7) RESTORE THE ELEMENT GRADIENT TO FULL SUBDIMENSION -C - CALL RANGE(IELF, IUMB, P3, P2, NDIMI, NSUBI, NS) - DO 50 J=1,NDIMI - GXI(J) = GXI(J) + P2(J) - INVARJ = INVAR(J) - IF (ISTATE(INVARJ).EQ.0) GXI(J) = 0.0 - 50 CONTINUE - RETURN - 60 CONTINUE - IRTN = 1 - RETURN - END - - SUBROUTINE VE08ID(IELF, NDIMI, NS, P1, P2, W1, W2, B, INVAR, - * ISTATE, NSUBI, RANGE) -C -C -C ********************************************************** -C * * -C * COMPUTATION OF THE PRODUCT OF THE IELF-TH ELEMENT * -C * HESSIAN WITH W1 INTO W2. * -C * * -C ********************************************************** -C -C -C -C - DOUBLE PRECISION W1(1), W2(1), B(1), P1(1), P2(1) - INTEGER INVAR(1), ISTATE(1) - EXTERNAL RANGE -C -C - IUB = 1 - IUTB = 3 -C -C BUILD THE CORRESPONDING DIRECTION VECTOR IN SUBDIMENSION -C -------------------------------------------------------- -C - DO 10 I=1,NDIMI - INVARI = INVAR(I) - P1(I) = W1(INVARI) - 10 CONTINUE -C -C PROJECT ON THE RANGE OF THE I-TH ELEMENT HESSIAN -C ------------------------------------------------ -C - CALL RANGE(IELF, IUB, P1, P2, NDIMI, NSUBI, NS) - IF (NSUBI.EQ.0) RETURN -C -C CHECK IF PRODUCT IS NEEDED -C -------------------------- -C - DO 20 J=1,NSUBI - IF (P2(J).NE.0.0) GO TO 30 - 20 CONTINUE - RETURN -C -C COMPUTE THE PRODUCT ON THE RANGE -C -------------------------------- -C - 30 CONTINUE - CALL VE08PD(B, P2, P1, NSUBI) -C -C RESTORE PRODUCT TO FULL SUBDIMENSION -C ------------------------------------ -C - CALL RANGE(IELF, IUTB, P1, P2, NDIMI, NSUBI, NS) -C -C ACCUMULATE THE EXPANDED PRODUCT IN THE VECTOR W2 -C ------------------------------------------------ -C - DO 40 J=1,NDIMI - IVAR = INVAR(J) - IF (ISTATE(IVAR).NE.0) W2(IVAR) = W2(IVAR) + P2(J) - 40 CONTINUE - RETURN - END - - SUBROUTINE VE08JD(IELF, NDIMI, NS, P1, X, XOLD, GI1, GI2, INVAR, - * ISTATE, ISELF, NGRI, DIFEST, GSINI, A, FNOISE, FMAX, FXI, GSA, - * IKLS, MAXLS, INFO, IFLAG, IRTN, ELFNCT) -C -C -C ***************************************************** -C * * -C * COMPUTATION OF THE VALUE OF THE IELF-TH ELEMENT * -C * FUNCTION AT THE POINT X, AND OF ITS GRADIENT, * -C * IF ANALYTICAL GRADIENTS ARE AVAILABLE. * -C * * -C ***************************************************** -C -C - DOUBLE PRECISION A, FMAX, FNOISE, FXI, GSA, GSINI, TEMP - LOGICAL DIFEST - EXTERNAL ELFNCT -C -C - DOUBLE PRECISION P1(1), GI1(1), GI2(1), X(1), XOLD(1) - INTEGER INVAR(1), ISTATE(1) -C -C - IRTN = 0 - NGRI = 0 - IF (ISELF.GE.8) GO TO 20 - DO 10 I=1,NDIMI - GI1(I) = GI2(I) - 10 CONTINUE - GO TO 40 -C -C PREPARE THE DATA POINT -C ---------------------- -C - 20 CONTINUE - ITEMP = ISELF - 10 - IF (.NOT.DIFEST) ISELF = ITEMP - IFFLAG = 2 - IF (ITEMP.LT.0) IFFLAG = 1 - DO 30 I=1,NDIMI - INVARI = INVAR(I) - P1(I) = X(INVARI) - 30 CONTINUE -C -C COMPUTE THE NEW VALUE OF THE I-TH ELEMENT FUNCTION -C -------------------------------------------------- -C - CALL ELFNCT(IELF, P1, FXI, GI1, NDIMI, NS, IFFLAG, FMAX, FNOISE) - NGRI = NGRI + 1 - CALL VE08BD(FNOISE, IFFLAG, IFLAG, INFO, JRTN) - IF (JRTN.EQ.1) GO TO 60 -C -C TEST IF OVERSHOOT IS DETECTED BY THE USER. IN THIS CASE, STOP -C ------------------------------------------------------------- -C EVALUATION AND BRANCH TO HALVE THE STEP -C --------------------------------------- -C - IF (IFFLAG.EQ.0 .AND. IKLS.LT.MAXLS) IRTN = 2 - IF (IFFLAG.EQ.0 .AND. IKLS.LT.MAXLS) RETURN -C -C ACCUMULATE THE NEW DIRECTIONAL DERIVATIVE -C ----------------------------------------- -C - 40 CONTINUE - IF (DIFEST) RETURN - DO 50 J=1,NDIMI - IVAR = INVAR(J) - ISVAR = ISTATE(IVAR) - IF (ISVAR.EQ.0) GO TO 50 - TEMP = X(IVAR)-XOLD(IVAR) - GSA = GSA + TEMP*GI1(J)/A - IF (ISVAR.LT.0 .OR. ISVAR.GT.2) GO TO 50 - GSINI = GSINI + TEMP*GI2(J) - 50 CONTINUE - RETURN - 60 CONTINUE - IRTN = 1 - RETURN - END - - SUBROUTINE VE08KD(IELF, NDIMI, NS, P1, P2, P3, G1, G2, X, XOLD, - * B, ISELF, INVAR, ISTATE, NSUBI, NGRI, FXI, FNOISE, XNOISE, - * EPSQRT, RELPR, TOLCNT, A, DXN1, GSA, GSINI, GXMAX, INFO, - * IFLAG, IRTN, ELFNCT, RANGE) -C -C -C ************************************************************* -C * * -C * COMPUTATION OF ONE ELEMENT GRADIENT BY DIFFERENCES, * -C * USING STEWART'S ALGORITHM. * -C * * -C ************************************************************* -C -C -C - DOUBLE PRECISION A, ALPHA, DALPHA, DBIG, DDENOM, DELTA, DERJ, - * DGAMMA, DISCR, DNUM, DPHI, DXN1, RELPR, EPSQRT, ETA, - * FNOISE, FXBWD, FXFWD, FXI, GAMMA, GAM2, GSA, GSINI, GXMAX, - * TEMP, TEST, TOLCNT, XI, XNOISE - EXTERNAL ELFNCT, RANGE -C -C - DOUBLE PRECISION P1(1), P2(1), P3(1), G1(1), G2(1), X(1), - * XOLD(1), B(1) - INTEGER ISTATE(1), INVAR(1) -C -C - IRTN = 0 - IUB = 1 - IUTB = 3 - IUMB = 2 - NGRI = 0 -C -C PROJECT THE ELEMENT GRADIENT ON THE RANGE -C ----------------------------------------- -C - CALL RANGE(IELF, IUB, G2, G1, NDIMI, NSUBI, NS) -C -C IF THE GRADIENT IS KNOWN ALREADY, COPY ITS VALUE -C ------------------------------------------------ -C - IF (ISELF.GE.8) GO TO 20 - DO 10 I=1,NDIMI - G1(I) = G2(I) - 10 CONTINUE - GO TO 210 - 20 CONTINUE - ISELF = ISELF - 10 - IF (ISELF.GE.0) RETURN - IF (NSUBI.EQ.0) GO TO 190 -C -C PREPARE THE DATA POINT -C ---------------------- -C - IFFLAG = 1 - DO 30 I=1,NDIMI - INVARI = INVAR(I) - P1(I) = X(INVARI) - 30 CONTINUE -C -C SET THE PARAMETERS FOR STEWART'S ALGORITHM -C ------------------------------------------ -C - DPHI = DMAX1(DABS(FXI),FNOISE,RELPR) - XI = 0.0 - DO 40 I=1,NDIMI - XI = XI + P1(I)**2 - 40 CONTINUE - XI = DSQRT(XI) -C -C LOOP ON THE TRUE SUBDIMENSION -C ----------------------------- -C - LL = 0 - NSUBIT = NSUBI - DO 150 J=1,NSUBI - LL = LL + J - INVARJ = INVAR(J) - IF (NSUBI.EQ.NDIMI .AND. ISTATE(INVARJ).EQ.0) GO TO 150 - ALPHA = B(LL) - DALPHA = DABS(ALPHA) - IF (DALPHA.GE.EPSQRT) GO TO 50 - DELTA = EPSQRT - TEST = 1.0 - GO TO 90 - 50 CONTINUE - GAMMA = G1(J) - DGAMMA = DABS(GAMMA) -C -C COMPUTE ETA, THE ACCURACY FACTOR -C -------------------------------- -C - ETA = DMAX1(FNOISE,(DGAMMA*XI*XNOISE)/DPHI) -C -C TEST IF CENTRAL DIFFERENCE IS NEEDED BECAUSE OF CANCELLATION ERRORS -C ------------------------------------------------------------------- -C IN THE COMPUTATION OF THE OVERAL GRADIENT -C ---------------------------------------- -C - IF (GXMAX.LE.100.0*EPSQRT*DGAMMA) GO TO 80 -C -C COMPUTE THE VALUE OF ACCURACY CONDITION -C --------------------------------------- -C - TEST = DALPHA*DPHI*ETA - GAM2 = GAMMA**2 - IF (GAM2.LE.TEST .AND. DGAMMA.GE.(100.*EPSQRT)) GO TO 60 -C -C FIRST ALTERNATIVE -C ----------------- -C - DELTA = 2.0*DSQRT(ETA*DPHI/DALPHA) - DNUM = DALPHA*DELTA - GO TO 70 -C -C SECOND ALTERNATIVE -C ------------------ -C - 60 CONTINUE - DELTA = 2.0*(ETA*DPHI*DGAMMA/ALPHA**2)**(1.0/3.0) - DNUM = GAMMA + GAMMA -C -C NEWTON'S REFINEMENT -C ------------------- -C - 70 CONTINUE - DDENOM = 3.0*DALPHA*DELTA + 4.0*DGAMMA - IF (DABS(DDENOM).GT.RELPR) DELTA = DELTA*(1.0-DNUM/DDENOM) - DELTA = DSIGN(DELTA,GAMMA) - DELTA = DSIGN(DELTA,ALPHA) -C -C TEST IF CENTRAL DIFFERENCE ARE NEEDED. IF YES, COMPUTE THE NEW DELTA -C -------------------------------------------------------------------- -C - TEST = 1.0 - IF (DABS(ALPHA*DELTA).LT.TOLCNT*DGAMMA .AND. DELTA.LE.A*DXN1) - * GO TO 90 - 80 CONTINUE - DISCR = GAM2 + 4.0*DALPHA*DPHI*ETA/TOLCNT - DELTA = -DGAMMA/DALPHA + DSQRT(DISCR)/DALPHA - TEST = -1.0 - 90 CONTINUE - TEMP = DMAX1(DABS(DELTA),0.001*EPSQRT) - DELTA = DSIGN(TEMP,DELTA) -C -C FORWARD DIFFERENCE: COMPUTE THE PERTURBED POINT -C ----------------------------------------------- -C - DO 100 K=1,NSUBI - P2(K) = 0.0 - 100 CONTINUE - P2(J) = 1.0 - CALL RANGE(IELF, IUTB, P2, P1, NDIMI, NSUBIT, NS) - DO 110 K=1,NDIMI - P2(K) = P1(K)*DELTA - INVARK = INVAR(K) - P1(K) = X(INVARK) + P2(K) - 110 CONTINUE -C -C FORWARD DIFFERENCE: COMPUTE THE PERTURBED VALUE OF THE ELEMENT -C -------------------------------------------------------------- -C FUNCTION -C -------- -C - CALL ELFNCT(IELF, P1, FXFWD, P3, NDIMI, NS, IFFLAG, DBIG, - * FNOISE) - NGRI = NGRI + 1 - CALL VE08BD(FNOISE, IFFLAG, IFLAG, INFO, JRTN) - IF (JRTN.EQ.1) GO TO 230 -C -C BACKWARD DIFFERENCE: COMPUTE THE SECOND PERTURBED DATA POINT -C ------------------------------------------------------------ -C - IF (TEST.GT.0.0) GO TO 130 - DO 120 K=1,NDIMI - INVARK = INVAR(K) - P1(K) = X(INVARK) - P2(K) - 120 CONTINUE -C -C BACKWARD DIFFERENCE: COMPUTE THE ELEMENT FUNCTION VALUE AT SECOND -C ----------------------------------------------------------------- -C PERTURBED DATA POINT -C -------------------- -C - CALL ELFNCT(IELF, P1, FXBWD, P3, NDIMI, NS, IFFLAG, DBIG, - * FNOISE) - NGRI = NGRI + 1 - CALL VE08BD(FNOISE, IFFLAG, IFLAG, INFO, JRTN) - IF (JRTN.EQ.1) GO TO 230 -C -C CENTRAL DIFFERENCE: ESTIMATE THE J-TH COMPONENT OF GRADIENT IN -C -------------------------------------------------------------- -C SUBDIMENSION -C ------------ -C - DERJ = (FXFWD-FXBWD)/(DELTA+DELTA) - GO TO 140 -C -C FORWARD DIFFERENCE: ESTIMATE THE J-TH COMPONENT OF GRADIENT IN -C -------------------------------------------------------------- -C SUBDIMENSION -C ------------ -C - 130 CONTINUE - DERJ = (FXFWD-FXI)/DELTA -C -C MODIFY THE J-TH COMPONENT OF THE PROJECTED GRADIENT -C --------------------------------------------------- -C - 140 CONTINUE - G1(J) = DERJ - G1(J) -C -C END OF THE LOOP ON THE SUBDIMENSIONS -C ------------------------------------ -C - 150 CONTINUE -C -C RESTORE THE ELEMENT GRADIENT TO THE FULL SUBDIMENSION -C ----------------------------------------------------- -C - IF (NDIMI.EQ.NSUBI) GO TO 170 - DO 160 I=1,NSUBI - P1(I) = G1(I) - 160 CONTINUE - CALL RANGE(IELF, IUMB, P1, G1, NDIMI, NSUBI, NS) - 170 CONTINUE - DO 180 J=1,NDIMI - G1(J) = G2(J) + G1(J) - INVARJ = INVAR(J) - IF (ISTATE(INVARJ).EQ.0) G1(J) = 0.0 - 180 CONTINUE - GO TO 210 - 190 CONTINUE - DO 200 I=1,NDIMI - G1(I) = G2(I) - 200 CONTINUE -C -C ACCUMULATE THE DIRECTIONAL DERIVATIVE -C ------------------------------------- -C - 210 CONTINUE - DO 220 J=1,NDIMI - IVAR = INVAR(J) - ISVAR = ISTATE(IVAR) - IF (ISVAR.EQ.0) GO TO 220 - TEMP = X(IVAR)-XOLD(IVAR) - GSA = GSA + TEMP*G1(J)/A - IF (ISVAR.LT.0 .OR. ISVAR.GT.2) GO TO 220 - GSINI = GSINI + TEMP*G2(J) - 220 CONTINUE - RETURN - 230 CONTINUE - IRTN = 1 - RETURN - END - - SUBROUTINE VE08LD(IELF, NDIMI, NS, P1, P2, P3, Y, S, D, B, RESN2, - * ISELF, INVAR, KU, PRERED, SCAL1, RELPR, GNOISE, INFO, IFLAG, - * IRTN, RANGE) -C -C -C -C **************************************************** -C * * -C * UPDATING OF THE IELF-TH ELEMENT HESSIAN * -C * * -C **************************************************** -C -C -C - DOUBLE PRECISION DSBSI, DYSI, RELPR, GNOISE, PRERED, Q, RESN2, - * RNI1, RNI2, RSI, RSII, SBSI, SNI2, TEMP, YSI - LOGICAL SCAL1 - EXTERNAL RANGE -C -C - DOUBLE PRECISION P1(1), P2(1), P3(1), Y(1), S(1), D(1), B(1) - INTEGER INVAR(1) -C -C -C - IRTN = 0 - PRERED = 0.0 - IUTB = 3 - IUTMB = 4 - IUB = 1 -C -C BUILD THE PRODUCT OF THE DIFFERENCE IN GRADIENTS TIMES THE STEP(YSI) -C -------------------------------------------------------------------- -C -C 1) BUILD THE UNPROJECTED DIFFERENCE IN P1 -C - DO 10 I=1,NDIMI - P1(I) = Y(I) - 10 CONTINUE -C -C 2) PROJECT IT ON THE RANGE -C - CALL RANGE(IELF, IUTMB, P1, P2, NDIMI, NSUBI, NS) - KU = (NSUBI*(NSUBI+1))/2 - IF (NSUBI.EQ.0) RETURN - IF (NSUBI.EQ.NDIMI) GO TO 30 -C -C 3) STORE THE RESULT IN Y -C - DO 20 I=1,NSUBI - Y(I) = P2(I) - 20 CONTINUE -C -C 4) BUILD THE UNPROJECTED STEP IN P1 -C - 30 CONTINUE - DO 40 I=1,NDIMI - INVARI = INVAR(I) - P1(I) = S(INVARI) - 40 CONTINUE -C -C 5) PROJECT THE STEP ON THE RANGE -C - CALL RANGE(IELF, IUB, P1, P2, NDIMI, NSUBI, NS) -C -C 6) COMPUTE YSI -C - YSI = 0.0 - DO 50 J=1,NSUBI - YSI = YSI + P2(J)*Y(J) - 50 CONTINUE -C -C COMPUTE THE PRODUCT OF THE PROJECTED STEP TIMES THE ELEMENT HESSIAN -C ------------------------------------------------------------------- -C -C -C 1) COMPUTE THE PRODUCT ON THE RANGE -C - CALL VE08PD(B, P2, P1, NSUBI) -C -C 2) COMPUTE THE PRODUCT OF THE RESULT WITH THE PROJECTED STEP (SBSI) -C AND THE SQUARE OF THE NORM OF THE PROJECTED STEP -C - SBSI = 0.0 - SNI2 = 0.0 - RNI2 = 0.0 - DO 60 J=1,NSUBI - SNI2 = SNI2 + P2(J)**2 - SBSI = SBSI + P1(J)*P2(J) -C -C 3) UPDATE THE SQUARE OF THE RESIDUAL NORM -C - RNI2 = RNI2 + (Y(J)-P1(J))**2 - 60 CONTINUE - RESN2 = RESN2 + RNI2 -C -C 4) UPDATE THE VALUE OF THE PREDICTED REDUCTION IN THE -C OBJECTIVE FUNCTION -C - PRERED = YSI - 0.5*SBSI -C -C CHECK IF THE NOISE ON THE DIFFERENCE IN GRADIENT IS BELOW THE -C ------------------------------------------------------------- -C NORM OF THE RESIDUAL. IF IT IS NOT THE CASE, SKIP UPDATING -C ---------------------------------------------------------- -C - RNI1 = DSQRT(RNI2) - IF (RNI1.LT.GNOISE) RETURN -C -C SELECT BFGS OR RK1 UPDATE, DEPENDING ON CONVEXITY OF THE ELEMENT -C ---------------------------------------------------------------- -C FUNCTION -C -------- -C - IF (YSI.GT.0.0 .AND. SBSI.GT.0.0 .AND. IABS(ISELF).LE.1) GO TO 120 - ISELF = ISIGN(2,ISELF) -C -C 1) THE ELEMENT FUNCTION IS NOT CONVEX: APPLY RK1 UPDATE. BUILD THE -C INNER PRODUCT OF THE RESIDUAL TIMES THE PROJECTED STEP, AND THE -C NORM OF THESE VECTORS -C - RSI = YSI - SBSI - DO 70 J=1,NSUBI - P1(J) = Y(J) - P1(J) - 70 CONTINUE -C -C 2) TEST IF THE DENOMINATOR IS BIG ENOUGH. IF THIS IS NOT THE CASE, -C AVOID UPDATING -C - TEMP = DABS(RSI) - IF (TEMP.LE.RELPR .OR. TEMP.LE.(0.1*DSQRT(SNI2)*RNI1)) RETURN -C -C 3) APPLY RK1 UPDATE -C - L = 0 - RSII = 1.0D0/RSI - DO 80 J=1,NSUBI - P2(J) = P1(J)*RSII - 80 CONTINUE - DO 100 J=1,NSUBI - DO 90 K=1,J - L = L + 1 - B(L) = B(L) + P1(J)*P2(K) - 90 CONTINUE - 100 CONTINUE -C -C 4) UPDATE THE DIAGONAL ACCORDINGLY -C - CALL RANGE(IELF, IUTB, P1, P2, NDIMI, NSUBI, NS) - DO 110 J=1,NDIMI - IVAR = INVAR(J) - D(IVAR) = D(IVAR) + P2(J)*P2(J)/RSI - 110 CONTINUE - RETURN -C -C 5) THE ELEMENT FUNCTION IS CONVEX: TEST FOR SCALING AT FIRST -C ITERATION -C - 120 CONTINUE - IF (SCAL1) GO TO 160 -C -C 6) SCALING OF THE ELEMENT HESSIAN -C - IF (SBSI.LE.RELPR .AND. YSI.LE.RELPR) RETURN - IF (SBSI.GT.RELPR) GO TO 130 - INFO = IELF - IFLAG = 30 - GO TO 260 - 130 CONTINUE - Q = YSI/SBSI - DO 140 J=1,KU - B(J) = Q*B(J) - 140 CONTINUE - DO 150 J=1,NSUBI - P1(J) = Q*P1(J) - 150 CONTINUE - SBSI = YSI - Q = Q - 1.0 -C -C 7) APPLY BFGS UPDATE -C - 160 CONTINUE - L = 0 - DYSI = 1./DSQRT(YSI) - DSBSI = 1./DSQRT(SBSI) - DO 170 J=1,NSUBI - P3(J) = Y(J)*DYSI - P2(J) = P1(J)*DSBSI - 170 CONTINUE - DO 190 J=1,NSUBI - DO 180 K=1,J - L = L + 1 - B(L) = B(L) + P3(J)*P3(K) - P2(J)*P2(K) - 180 CONTINUE - 190 CONTINUE -C -C 8) UPDATE THE DIAGONAL ACCORDINGLY -C - CALL RANGE(IELF, IUTB, P1, P2, NDIMI, NSUBI, NS) - DO 200 J=1,NDIMI - P1(J) = (P2(J)**2)/SBSI - 200 CONTINUE - CALL RANGE(IELF, IUTB, Y, P2, NDIMI, NSUBI, NS) - DO 210 J=1,NDIMI - Y(J) = (P2(J)**2)/YSI - P1(J) - 210 CONTINUE - NDIMIT = NDIMI - DO 250 J=1,NDIMI - IVAR = INVAR(J) - IF (SCAL1) GO TO 240 - DO 220 K=1,NDIMI - P1(K) = 0.0 - 220 CONTINUE - P1(J) = 1.0 - IF (NSUBI.EQ.NDIMI) GO TO 230 - CALL RANGE(IELF, IUB, P1, P2, NDIMIT, NSUBI, NS) - CALL RANGE(IELF, IUTB, P2, P1, NDIMIT, NSUBI, NS) - 230 CONTINUE - D(IVAR) = D(IVAR) + Q*P1(J) - 240 CONTINUE - D(IVAR) = D(IVAR) + Y(J) - 250 CONTINUE - RETURN - 260 CONTINUE - IRTN = 1 - RETURN - END - - SUBROUTINE VE08MD(RES,X,Y,N) -C -C -C ********************************************** -C * * -C * INNER PRODUCT OF X AND Y IN RES * -C * * -C ********************************************** -C - - INTEGER N - DOUBLE PRECISION RES, X(N), Y(N) - - - RES = 0.0 - DO 10 I=1,N - RES = RES + X(I)*Y(I) -10 CONTINUE - RETURN - END - - SUBROUTINE VE08ND(DX, W1, STEP, GX, W2, W3, RES, N) -C -C -C ********************************************** -C * * -C * FIRST INNER LOOP OF CONJUGATE GRADIENT * -C * * -C ********************************************** -C -C - DOUBLE PRECISION RES, STEP -C -C - DOUBLE PRECISION DX(1), W1(1), W2(1), W3(1), GX(1) -C -C - RES = 0.0 - DO 10 I=1,N - DX(I) = DX(I) + STEP*W1(I) - GX(I) = GX(I) - STEP*W2(I) - W2(I) = GX(I)*W3(I) - RES = RES + W2(I)*GX(I) - 10 CONTINUE - RETURN - END - - SUBROUTINE VE08OD(R, X, Y, BETA, N) -C -C -C ************************************************ -C * * -C * SECOND INNER LOOP OF CONJUGATE GRADIENTS * -C * * -C ************************************************ -C -C - DOUBLE PRECISION BETA -C -C - DOUBLE PRECISION R(1), X(1), Y(1) -C -C - DO 10 I=1,N - R(I) = X(I) + BETA*Y(I) - 10 CONTINUE - RETURN - END - - SUBROUTINE VE08PD(H, X, Y, NSUBI) -C -C -C ********************************************************* -C * * -C * THIS ROUTINE COMPUTES THE PRODUCT H*X=Y, WHERE * -C * H IS AN ELEMENT HESSIAN MATRIX OF DIMENSION * -C * NSUBI*NSUBI AND WHERE X AND Y AND TWO VECTORS * -C * OF DIMENSION NSUBI. * -C * * -C ********************************************************* -C -C -C - DOUBLE PRECISION H(1), X(1), Y(1) -C -C -C -C RESET THE RESULT VECTOR Y TO ZERO -C --------------------------------- -C - DO 10 I=2,NSUBI - Y(I) = 0.0 - 10 CONTINUE -C -C COMPUTE THE PRODUCT -C ------------------- -C - L = 1 - Y(1) = H(1)*X(1) - IF (NSUBI.LT.2) RETURN - DO 30 I=2,NSUBI - ITEMP = I - 1 - DO 20 J=1,ITEMP - L = L + 1 - Y(I) = Y(I) + H(L)*X(J) - Y(J) = Y(J) + H(L)*X(I) - 20 CONTINUE - L = L + 1 - Y(I) = Y(I) + H(L)*X(I) - 30 CONTINUE - RETURN - END - - SUBROUTINE VE08QD(KIT, FX, GXN, NGR, DNGR, ALPHA, ITSTEP, BOUNDS, - * X, GX, S, B, IP1, NLP, N, NS, M, NVAR, P1, P2, RANGE) -C -C -C **************************************************** -C * * -C * PRINTING SUBROUTINE FOR PARTIALLY SEPARABLE * -C * MINIMIZATION * -C * * -C **************************************************** -C -C - DOUBLE PRECISION ALPHA, DNGR, FX, GXN - INTEGER KIT(2), NGR(2) - LOGICAL BOUNDS - EXTERNAL RANGE -C -C -C THE PRINTOUT GENERATED BY THIS ROUTINE IS GOVERNED BY THE VALUE -C OF THE ARGUMENTS IP1 AS FOLLOWS : -C -C IP1 DETERMINE THE ITEMS TO BE PRINTED : -C -C IP1=0 THE ITERATION COUNT (KIT(2)), THE FUNCTION VALUE (FX), -C THE NORM OF THE GRADIENT (GXN) AND THE NUMBER OF -C OBJECTIVE FUNCTION SUBROUTINE CALLS (NGR(2),DNGR) -C IP1=1 AS WITH IP1=0 + STEP SELECTION ITERATIONS (ITSTEP) -C AND THE LINESEARCH PARAMETER (ALPHA) -C IP1=2 AS WITH IP1=1 + THE VECTOR OF VARIABLES (X(K),K=1,N) -C IP1=3 AS WITH IP1=2 + THE GRADIENT VECTOR (GX(K),K=1,N) -C IP1=4 AS WITH IP1=3 + THE RECENT CHANGE IN THE VARIABLES -C (S(K),K=1,N) -C IP1=5 AS WITH IP1=4 + THE VALUES OF THE NONZERO ELEMENTS IN -C THE APPROXIMATE HESSIAN MATRIX (B(K),K=1,M) -C -C -C THE ARRAYS APPEARING IN THE ARGUMENT LIST ARE DIMENSIONED BY - DOUBLE PRECISION X(1), GX(1), S(1), B(1), P1(1), P2(1) - INTEGER NVAR(1) -C -C PRINT THE REQUIRED ITEMS -C - WRITE (NLP,99999) KIT(2), DNGR, NGR(2) - IF (.NOT.BOUNDS) WRITE (NLP,99998) FX, GXN - IF (BOUNDS) WRITE (NLP,99997) FX, GXN - IF (IP1.LE.0) RETURN - WRITE (NLP,99996) ITSTEP, ALPHA - IF (IP1.LE.1) RETURN - WRITE (NLP,99995) - WRITE (NLP,99994) (X(K),K=1,N) - IF (IP1.LE.2) RETURN - WRITE (NLP,99993) - WRITE (NLP,99994) (GX(K),K=1,N) - IF (IP1.LE.3) RETURN - IF (KIT(2).GE.0) WRITE (NLP,99992) - IF (KIT(2).LT.0) WRITE (NLP,99991) - WRITE (NLP,99994) (S(K),K=1,N) - IF (IP1.LE.4) RETURN - L = 0 - NST = NS - DO 20 J=1,NS - WRITE (NLP,99989) J - JL = NVAR(J) - NDIMI = NVAR(J+1) - JL - IELF = J - CALL RANGE(IELF, 1, P1, P2, NDIMI, NSUBI, NST) - IF (NSUBI.EQ.0) GO TO 20 - DO 10 K=1,NSUBI - LP1 = L + 1 - LPK = L + K - WRITE (NLP,99988) (B(KK),KK=LP1,LPK) - L = LPK - 10 CONTINUE - 20 CONTINUE - WRITE (NLP,99990) - RETURN -C -C FORMATS -C -99999 FORMAT (///4X, 15H**** ITERATION , I4, 11H ( AFTER , F9.2, 2H (, - * I7, 17H) FUNCTION CALLS)) -99998 FORMAT (//5X, 17HFUNCTION VALUE = , D23.15//5X, 14HNORM OF THE GR, - * 9HADIENT = , D16.7) -99997 FORMAT (//5X, 17HFUNCTION VALUE = , D23.15//5X, 14HNORM OF THE PR, - * 19HOJECTED GRADIENT = , D16.7) -99996 FORMAT (/5X, 25HSTEP STRATEGY ITERATIONS , I7, 5X, 10HLINESEARCH, - * 11H PARAMETER , D10.3) -99995 FORMAT (/5X, 22HAPPROXIMATE SOLUTION X/) -99994 FORMAT (5X, 5D15.7) -99993 FORMAT (/5X, 34HGRADIENT AT APPROXIMATE SOLUTION X/) -99992 FORMAT (/5X, 26HLAST STEP OF THE ALGORITHM/) -99991 FORMAT (/5X, 23HORIGINAL STARTING POINT/) -99990 FORMAT (/5X, 50(1H*)//) -99989 FORMAT (/5X, 10HPROJECTED , I4, 19H-TH ELEMENT HESSIAN/) -99988 FORMAT (5X, 10D12.2) - END - - SUBROUTINE VE08RD(X, GX, ISTATE, N, EBOUND, XUPPER, XLOWER) -C -C -C ************************************************************ -C * * -C * THIS ROUTINE REDEFINES THE STATE OF THE BOUNDED * -C * VARIABLES WITH RESPECT TO THEIR BOUND, ACCORDING * -C * TO A STRATEGY PROPOSED BY BERTSEKAS. * -C * * -C ************************************************************ -C -C - DOUBLE PRECISION DBI, EBOUND, EPSIT, UBI, WI, XUPPER, XLOWER - EXTERNAL XUPPER, XLOWER -C -C -C - DOUBLE PRECISION X(1), GX(1) - INTEGER ISTATE(1) -C -C -C*********************************************************************** -C -C DETERMINE THE TOLERANCE ASSOCIATED WITH THE PRESENT ITERATE -C -C*********************************************************************** -C -C -C LOOP ON THE VARIABLES -C --------------------- -C - EPSIT = 0.0 - DO 30 I=1,N - IF (ISTATE(I).LE.0) GO TO 30 - IVAR = I -C -C BUILD THE GRADIENT STEP FROM THE PRESENT ITERATE -C ------------------------------------------------ -C - WI = X(I) - GX(I) -C -C CHECK GRADIENT STEP W.R.T. THE UPPER BOUND -C ------------------------------------------ -C - UBI = XUPPER(IVAR) - IF (WI.LE.UBI) GO TO 10 - WI = UBI - GO TO 20 -C -C CHECK GRADIENT STEP W.R.T. THE LOWER BOUND -C ------------------------------------------ -C - 10 CONTINUE - DBI = XLOWER(IVAR) - IF (WI.LT.DBI) WI = DBI -C -C ACCUMULATE THE SQUARE OF THE TOLERANCE -C -------------------------------------- -C - 20 CONTINUE - EPSIT = EPSIT + (X(I)-WI)**2 -C -C END OF THE LOOP ON THE VARIABLES: DEFINE THE FINAL TOLERANCE -C ------------------------------------------------------------ -C - 30 CONTINUE - EPSIT = DSQRT(EPSIT) - EPSIT = DMIN1(EPSIT,EBOUND) -C -C -C*********************************************************************** -C -C REDEFINE THE STATE OF EACH BOUNDED VARIABLE -C -C*********************************************************************** -C -C -C LOOP ON THE BOUNDED VARIABLES -C ----------------------------- -C - DO 60 I=1,N - IF (ISTATE(I).LE.0) GO TO 60 - IVAR = I - ISTATE(I) = 3 -C -C CHECK THE I-TH VARIABLE W.R.T. ITS UPPER BOUND -C ---------------------------------------------- -C - UBI = XUPPER(IVAR) - WI = UBI - X(I) - IF (WI.GT.EPSIT) GO TO 40 - IF (GX(I).LT.0.0) GO TO 50 -C -C CHECK THE I-TH VARIABLE W.R.T. ITS LOWER BOUND -C ---------------------------------------------- -C - 40 CONTINUE - DBI = XLOWER(IVAR) - WI = X(I) - DBI - IF (WI.GT.EPSIT) GO TO 60 - IF (GX(I).LE.0.0) GO TO 60 -C -C THE I-TH VARIABLE SHOULD BE IN THE DANGEROUS SET. SET THE VALUE -C --------------------------------------------------------------- -C OF ISTATE ACCORDINGLY -C --------------------- -C - 50 CONTINUE - IF (WI.LE.0.0) ISTATE(I) = 1 - IF (WI.GT.0.0) ISTATE(I) = 2 -C -C END OF THE LOOP ON THE BOUNDED VARIABLES -C ---------------------------------------- -C - 60 CONTINUE - RETURN - END - - SUBROUTINE VE08SD(EPSIL, NIT, NGR, DNGR, ITMAX, NFMAX, GXN, - * NEGCUR, DIFEST, X, FX, GX, DX, DF, N, INFO, IFLAG) -C -C -C ********************************************************* -C * * -C * THIS ROUTINE IMPLEMENTS THE CLASSICAL STOPPING * -C * CRITERION FOR MINIMIZATION ROUTINES. * -C * * -C * NORM USED: ORDINARY 2-NORM. * -C * * -C ********************************************************* -C -C - DOUBLE PRECISION DF, DNGR, DNGRMX, DX, EPSIL, EPS2, EPS3 - DOUBLE PRECISION FX, GXN, EPS - LOGICAL NEGCUR, DIFEST -C - DOUBLE PRECISION X(1), GX(1) - INTEGER NGR(2), NIT(2) -C -C -C -C INITIALIZE THE MAIN ACCURACY CONSTANTS -C -------------------------------------- -C - EPS = EPSIL - IF (DIFEST) EPS = DMAX1(EPSIL,1.0D-4) - EPS2 = EPS*EPS - EPS3 = 0.001*EPS2 -C -C SET RETURN FLAG TO DISABLED -C --------------------------- -C - IFLAG = 0 -C -C TEST FOR SUFFICIENT ACCURACY -C ---------------------------- -C - IF (GXN.GT.EPS .OR. NEGCUR) GO TO 10 - IFLAG = 1 - RETURN -C -C CHECK IF THE MAXIMUM NUMBER OF FUNCTION CALLS HAS BEEN EXCEEDED -C --------------------------------------------------------------- -C - 10 CONTINUE - DNGRMX = FLOAT(NFMAX) - IF (DNGR.LT.DNGRMX) GO TO 20 - IFLAG = 2 - RETURN -C -C CHECK IF PROGRESS OF THE ALGORITHM IS SUFFICIENT -C ------------------------------------------------ -C - 20 CONTINUE - IF (DF.GT.EPS3 .OR. DX.GT.EPS3) RETURN - IFLAG = 3 - RETURN - END - - SUBROUTINE VE08TD(A, X, FX, GXS, SNORM, S, N) -C -C -C *************************************************************** -C * * -C * THIS SUBROUTINE ALLOWS THE USER TO VERIFY THE PLAUSIBILTY * -C * OF THE VALUE TO BE ATTRIBUTED TO THE LINESEARCH * -C * PARAMETER. IF THE VALUE PROPOSED IN A IS NOT SUITABLE, * -C * THE USER SHOULD RETURN (IN A) A MORE APPROPRIATE ONE. * -C * * -C * THE LINESEARCH IS STARTING FROM THE POINT X IN DIRECTION * -C * S (OF NORM SNORM). THE INITIAL FUNCTION AND DIRECTIONAL * -C * DERIVATIVE VALUES ARE GIVEN IN FX AND GXS. THE NUMBER * -C * OF VARIABLES IS N. * -C * * -C *************************************************************** -C - DOUBLE PRECISION X(1), S(1), A, FX, GXS, SNORM -C -C -C - RETURN - END //GO.SYSIN DD ve08ad.f echo ve08ts.f 1>&2 sed >ve08ts.f <<'//GO.SYSIN DD ve08ts.f' 's/^-//' -*From CUNYVM.CUNY.EDU!phtoint%NEWTON Mon Nov 5 10:34 +01 1990 -C -C ********************************************* -C * * -C * TEST OF MINIMIZATION ROUTINE VE08AD * -C * * -C ********************************************* -C -C - DOUBLE PRECISION DIFGRD, EBOUND, EPSIL, FLOWBD, FX, RELPR, STMAX - LOGICAL FKNOWN, TESTGX, RESTRT, HESDIF -C - DOUBLE PRECISION STEPL(2), X(100), WK(1000) - INTEGER ISTATE(300), NVAR(101), INVAR(700), LWK(2), NGR(2), NIT(2) -C - COMMON /TEST/ ITEST -C - EXTERNAL ELFTS, RGETS, XLWTS, XUPTS -C -C LOOP OVER THE TEST PROBLEMS -C - NPROB = 6 - DO 30 ITEMP=1,NPROB -C -C READ IN SOME PARAMETERS -C ----------------------- -C -C 1) INDEX OF TEST PROBLEM CONSIDERED -C - ITEST = ITEMP -C -C 2) DIMENSION -C - N = 50 - IF (ITEMP.EQ.4) N = 3 - IF (ITEMP.EQ.6) N = 10 -C -C 3) AVAILABILITY OF ANALYTICAL DERIVATIVES -C - IDER = MOD(ITEMP,2) -C -C 4) OUTPUT PARAMETERS (IPFREQ, IPWHAT) -C - IPFREQ = 0 - IF (ITEMP.EQ.1) IPFREQ = 1 - IPWHAT = 0 - IF (ITEMP.EQ.4) IPWHAT = 5 -C -C 5) ESTIMATION OF FIRST HESSIAN MATRIX AT THE STARTING POINT ? -C - HESDIF = .FALSE. - IF (ITEMP.EQ.3 .OR. ITEMP.EQ.4) HESDIF = .TRUE. -C -C 6) ACCURACY -C - EPSIL = 0.0000001 - IF (ITEMP.EQ.5) EPSIL = 0.0 -C -C SET THE LENGTH OF THE WORK SPACE -C -------------------------------- -C - LWK(1) = 1000 -C -C DEFINE OUTPUT DEVICE -C -------------------- -C - IPDEVC = 6 -C -C SET THE MAXIMUM NUMBER OF FUNCTION CALLS -C ---------------------------------------- -C - NGR(1) = 1000 - NIT(1) = 500 -C -C DEFINE TOLERANCE FOR ANTI-ZIGZAG DEVICE (BERTSEKAS) -C --------------------------------------------------- -C - EBOUND = 1.0D-4 -C -C INITIALIZE THE LOGICAL VARIABLE THAT ASKS FOR TESTING OF THE -C ------------------------------------------------------------ -C ANALYTICAL GRADIENT -C ------------------- -C - TESTGX = .TRUE. -C -C SET THE RESTART PARAMETER TO FALSE -C ---------------------------------- -C - RESTRT = .FALSE. -C -C DEFINE STEPSIZE FOR DIFFERENCE ESTIMATION OF THE FIRST GRADIENT -C --------------------------------------------------------------- -C - DIFGRD = 1.0D-7 -C -C DEFINE THE MACHINE PRECISION -C ---------------------------- -C - RELPR = -1.0 - IF (ITEMP.EQ.6) RELPR = 1.0D-7 -C -C DEFINE THE OVERAL MAXIMUM STEPLENGTH AND THE MAXIMUM STEP -C --------------------------------------------------------- -C AT THE FIRST ITERATION -C ---------------------- -C - STMAX = 1.0D+10 - STEPL(1) = 1.D+5 -C -C INITIALIZE THE PROCESS FOR THE CONSIDERED TEST PROBLEM -C ------------------------------------------------------ -C - CALL INIELF(INVAR, NVAR, ISTATE, X, FLOWBD, FKNOWN, N, NS) -C -C IF THE DERIVATIVES ARE NOT AVAILABLE, MODIFY ISTATE -C --------------------------------------------------- -C - IF (IDER.NE.0) GO TO 20 - ITEMP1 = N + 1 - ITEMP2 = N + NS - DO 10 I=ITEMP1,ITEMP2 - ISTATE(I) = -1 - 10 CONTINUE - 20 CONTINUE -C -C MINIMIZE -C -------- -C - CALL VE08AD(N, NS, INVAR, NVAR, ELFTS, RGETS, XLWTS, XUPTS, - * X, FX, EPSIL, EBOUND, NGR, NIT, FKNOWN, FLOWBD, RELPR, DIFGRD, - * RESTRT, TESTGX, HESDIF, STMAX, STEPL, ISTATE, IPDEVC, IPFREQ, - * IPWHAT, LWK, WK, INFO, IFLAG) -C -C END OF TEST -C ----------- -C - 30 CONTINUE - STOP - END - SUBROUTINE INIELF(INVAR, NVAR, ISTATE, X, FLOWBD, FKNOWN, N, NS) -C -C -C ****************************************************** -C * * -C * INITIALIZATION OF THE IELF-TH ELEMENT FUNCTION * -C * * -C ****************************************************** -C -C - DOUBLE PRECISION DN, FLOWBD, XTEMP -C - DOUBLE PRECISION X(1) - INTEGER INVAR(1), NVAR(1), ISTATE(1) - LOGICAL FKNOWN -C - COMMON /TEST/ ITEST -C -C -C DECIDE WHAT FUNCTION IS USED -C ---------------------------- -C - GO TO (50, 50, 50, 10, 20, 80), ITEST -C -C FIRST TRIAL TRIVIAL PROBLEM (LOWER BOUND) -C ----------------------------------------- -C - 10 CONTINUE - FKNOWN = .FALSE. - N = 3 - NS = 2 - X(1) = 10.0 - X(2) = 4.0 - X(3) = 10.0 - NVAR(1) = 1 - NVAR(2) = 2 - NVAR(3) = 4 - INVAR(1) = 1 - INVAR(2) = 2 - INVAR(3) = 3 - ISTATE(1) = 1 - ISTATE(2) = -1 - ISTATE(3) = -1 - ISTATE(4) = 1 - ISTATE(5) = 1 - FLOWBD = 0.0 - RETURN -C -C EXTENDED ROSENBROCK -C ------------------- -C - 20 CONTINUE - FKNOWN = .FALSE. - NS = N - 1 - FLOWBD = 0.0 - DO 30 I=1,NS - IP = I + I - 1 - NVAR(I) = IP - INVAR(IP) = I - INVAR(IP+1) = I + 1 - NPI = N + I - ISTATE(NPI) = 1 - 30 CONTINUE - NVAR(NS+1) = NS + NS + 1 - ITEMP = -1 - DN = FLOAT(N) + 1.0 - DO 40 I=1,N - ISTATE(I) = ITEMP - X(I) = -1.0 - 40 CONTINUE - RETURN -C -C BROYDEN TRIDIAGONAL NONLINEAR SYSTEM (BROY3) -C -------------------------------------------- -C - 50 CONTINUE - FKNOWN = .FALSE. - FLOWBD = 0.0 - NS = N - 2 - IF (NS.LE.0) STOP - DO 60 I=1,NS - IP = 3*I - 2 - NVAR(I) = IP - INVAR(IP) = I - INVAR(IP+1) = I + 1 - INVAR(IP+2) = I + 2 - NPI = N + I - ISTATE(NPI) = 1 - 60 CONTINUE - NVAR(NS+1) = 3*NS + 1 - XTEMP = -1.0 - ITEMP = 1 - ITEMP1 = N - 1 - DO 70 I=2,ITEMP1 - X(I) = XTEMP - ISTATE(I) = ITEMP - 70 CONTINUE - X(1) = 0.0 - ISTATE(1) = 0 - X(N) = 0.0 - ISTATE(N) = 0 - RETURN -C -C BROYDEN BANDED FUNCTION (BROYB) -C ----------------------- -C - 80 CONTINUE - FKNOWN = .FALSE. - FLOWBD = 0.0 - NS = N - IP = 1 - DO 100 I=1,NS - NVAR(I) = IP - ILOW = MAX0(1,I-5) - IUP = MIN0(N,I+1) - DO 90 J=ILOW,IUP - INVAR(IP) = J - IP = IP + 1 - 90 CONTINUE - NPI = N + I - ISTATE(NPI) = 1 - 100 CONTINUE - NVAR(NS+1) = IP - DO 110 I=1,N - ISTATE(I) = -1 - X(I) = -1.0 - 110 CONTINUE - RETURN - END - SUBROUTINE RGETS(K, MODE, W1, W2, NDIMK, NSUBK, NS) -C -C -C * -C * TRANSFERS ON THE RANGE OF THE HESSIAN OF THE K-TH ELEMENT -C * -C -C MODE=1 <=> U*W1=W2 -C MODE=2 <=> U*W2=W1 -C MODE=3 <=> U'*W1=W2 -C MODE=4 <=> U'*W2=W1 -C -C - DOUBLE PRECISION W1(1), W2(1) -C -C - COMMON /TEST/ ITEST -C -C -C DECIDE WHAT TEST PROBLEM IS CONSIDERED -C -------------------------------------- -C - GO TO (50, 50, 50, 30, 10, 10), ITEST -C -C NO NULLSPACE AT ALL -C ------------------- -C - 10 CONTINUE - DO 20 I=1,NDIMK - W2(I) = W1(I) - 20 CONTINUE - NSUBK = NDIMK - RETURN -C -C CORRECT NULLSPACE FOR TEST PROBLEM 1 -C ------------------------------------ -C - 30 CONTINUE - IF (K.GT.1) GO TO 40 - NSUBK = 0 - RETURN - 40 CONTINUE - NSUBK = 2 - W2(1) = W1(1) - W2(2) = W1(2) - RETURN -C -C BROYDEN TRIDIAGONAL NONLINEAR SYSTEM (BROY3) -C -------------------------------------------- -C - 50 CONTINUE - GO TO (60, 100, 80, 120), MODE -C -C 1) UW1=W2C -C - 60 CONTINUE - IF (K.EQ.1 .OR. K.EQ.NS) GO TO 70 - NSUBK = 2 - W2(1) = W1(1) + W1(3) + W1(3) - W2(2) = W1(2) - RETURN - 70 CONTINUE - NSUBK = 3 - W2(1) = W1(1) - W2(2) = W1(2) - W2(3) = W1(3) - RETURN -C -C 2) U'W1=W2 -C - 80 CONTINUE - IF (K.EQ.1 .OR. K.EQ.NS) GO TO 90 - W2(1) = W1(1) - W2(2) = W1(2) - W2(3) = W2(1) + W2(1) - RETURN - 90 CONTINUE - W2(1) = W1(1) - W2(2) = W1(2) - W2(3) = W1(3) - RETURN -C -C 3) UW2=W1 -C - 100 CONTINUE - IF (K.EQ.1 .OR. K.EQ.NS) GO TO 110 - W2(1) = 0.2*W1(1) - W2(2) = W1(2) - W2(3) = W2(1) + W2(1) - RETURN - 110 CONTINUE - W2(1) = W1(1) - W2(2) = W1(2) - W2(3) = W1(3) - RETURN -C -C 4) U'W2=W1 -C - 120 CONTINUE - IF (K.EQ.1 .OR. K.EQ.NS) GO TO 130 - NSUBK = 2 - W2(1) = W1(1) - W2(2) = W1(2) - RETURN - 130 CONTINUE - NSUBK = 3 - W2(1) = W1(1) - W2(2) = W1(2) - W2(3) = W1(3) - RETURN - END - DOUBLE PRECISION FUNCTION XLWTS(IVAR) -C -C -C *************************************** -C * * -C * LOWER BOUNDS ON THE VARIABLES * -C * * -C *************************************** -C - COMMON /TEST/ ITEST -C -C -C DECIDE WHICH TEST IS CONDIDERED -C ------------------------------- -C - GO TO (30, 30, 30, 20, 10, 10), ITEST -C -C NO LOWER BOUND -C -------------- -C - 10 CONTINUE - XLWTS = -1.0D20 - RETURN -C -C POSITIVE VARIABLE -C ----------------- -C - 20 CONTINUE - XLWTS = 0.0 - RETURN -C -C BOUNDED BROYDEN TRIDIAGONAL PROBLEM -C ----------------------------------- -C - 30 CONTINUE - XLWTS = 0.65 - RETURN - END - DOUBLE PRECISION FUNCTION XUPTS(IVAR) -C -C -C *************************************** -C * * -C * UPPER BOUNDS ON THE VARIABLES * -C * * -C *************************************** -C - COMMON /TEST/ ITEST -C -C -C DECIDE WHICH TEST IS CONDIDERED -C ------------------------------- -C - GO TO (20, 20, 20, 10, 10, 10), ITEST -C -C NO UPPER BOUND -C -------------- -C - 10 CONTINUE - XUPTS = 1.0D+20 - RETURN -C -C BOUNDED BROYDEN TRIDIAGONAL PROBLEM -C ----------------------------------- -C - 20 CONTINUE - XUPTS = 0.71 - RETURN - END - SUBROUTINE ELFTS(K, X, FX, GX, NDIMK, NS, JFLAG, FMAX, - * FNOISE) -C -C -C ************************************************** -C * * -C * COMPUTATION OF THE ELEMENT FUNCTION VALUE * -C * * -C ************************************************** -C -C - DOUBLE PRECISION FMAX, FNOISE, FX, TEMP, XX -C - DOUBLE PRECISION X(1), GX(1) -C - COMMON /TEST/ ITEST -C -C -C DECIDE WHICH TEST IS CONDIDERED -C ------------------------------- -C - GO TO (40, 40, 40, 10, 30, 50), ITEST -C -C FIRST TRIVIAL TEST PROBLEM -C -------------------------- -C - 10 CONTINUE - IF (K.GT.1) GO TO 20 - FX = X(1) - FNOISE = 0.0 - IF (JFLAG.LT.2) RETURN - GX(1) = 1.0 - RETURN - 20 CONTINUE - FX = 0.5*(X(1)-X(2))**2 + X(1)**2 - FNOISE = 0.0 - IF (JFLAG.LT.2) RETURN - GX(1) = X(1) - X(2) + X(1) + X(1) - GX(2) = X(2) - X(1) - RETURN -C -C EXTENDED ROSENBROCK FUNCTION -C ---------------------------- -C - 30 CONTINUE - XX = X(2) - X(1)**2 - FX = 100.0*XX*XX + (X(1)-1.0)**2 - FNOISE = 0.0 - IF (JFLAG.LT.2) RETURN - GX(1) = -400.0*XX*X(1) + 2.0*(X(1)-1.0) - GX(2) = 200.0*XX - RETURN -C -C BROYDEN TRIDIAGONAL NONLINEAR SYSTEM -C ------------------------------------ -C - 40 CONTINUE - TEMP = (3.0-X(2)-X(2))*X(2) - X(1) - X(3) - X(3) + 1.0 - FX = TEMP*TEMP - TEMP = TEMP + TEMP - FNOISE = 0.0 - IF (JFLAG.LT.2) RETURN - GX(1) = -TEMP - GX(2) = TEMP*(3.0-4.0*X(2)) - GX(3) = -TEMP - TEMP - RETURN -C -C BROYDEN BANDED FUNCTION (BROYB) -C ------------------------------- -C - 50 CONTINUE - ILOW = MAX0(1,K-5) - IUP = MIN0(NS,K+1) - FX = 1.0 - DO 70 J=ILOW,IUP - JJ = J - ILOW + 1 - IF (J.EQ.K) GO TO 60 - FX = FX - X(JJ)*(1.0+X(JJ)) - GO TO 70 - 60 CONTINUE - FX = FX + X(JJ)*(2.0+5.0*X(JJ)*X(JJ)) - 70 CONTINUE - TEMP = FX + FX - FX = FX*FX - FNOISE = 0.0 - IF (JFLAG.LT.2) RETURN - DO 80 J=ILOW,IUP - JJ = J - ILOW + 1 - GX(JJ) = 0.0 - 80 CONTINUE - DO 90 J=ILOW,IUP - JJ = J - ILOW + 1 - IF (J.EQ.K) GX(JJ) = TEMP*(2.0+15.0*X(JJ)*X(JJ)) - IF (J.NE.K) GX(JJ) = -TEMP*(1.0+2.0*X(JJ)) - 90 CONTINUE - RETURN - END //GO.SYSIN DD ve08ts.f echo ve08ts.out 1>&2 sed >ve08ts.out <<'//GO.SYSIN DD ve08ts.out' 's/^-//' - - - - **** ITERATION -1 ( AFTER 3.96 ( 190) FUNCTION CALLS) - - - FUNCTION VALUE = 0.387020000000000D+01 - - NORM OF THE PROJECTED GRADIENT = 0.7196225D+01 - - - - **** ITERATION 1 ( AFTER 4.96 ( 238) FUNCTION CALLS) - - - FUNCTION VALUE = 0.246460256000000D+01 - - NORM OF THE PROJECTED GRADIENT = 0.8726189D+00 - - - - **** ITERATION 2 ( AFTER 5.88 ( 282) FUNCTION CALLS) - - - FUNCTION VALUE = 0.243145821034848D+01 - - NORM OF THE PROJECTED GRADIENT = 0.1551483D+00 - - - - **** ITERATION 3 ( AFTER 6.79 ( 326) FUNCTION CALLS) - - - FUNCTION VALUE = 0.243050506851653D+01 - - NORM OF THE PROJECTED GRADIENT = 0.2122874D-01 - - - - **** ITERATION 4 ( AFTER 7.71 ( 370) FUNCTION CALLS) - - - FUNCTION VALUE = 0.243048146951896D+01 - - NORM OF THE PROJECTED GRADIENT = 0.5156672D-02 - - - - **** ITERATION 5 ( AFTER 8.63 ( 414) FUNCTION CALLS) - - - FUNCTION VALUE = 0.243047999376789D+01 - - NORM OF THE PROJECTED GRADIENT = 0.4939792D-03 - - - - **** ITERATION 6 ( AFTER 9.54 ( 458) FUNCTION CALLS) - - - FUNCTION VALUE = 0.243047998035647D+01 - - NORM OF THE PROJECTED GRADIENT = 0.2001071D-03 - - - - **** ITERATION 7 ( AFTER 10.46 ( 502) FUNCTION CALLS) - - - FUNCTION VALUE = 0.243047997839509D+01 - - NORM OF THE PROJECTED GRADIENT = 0.2673273D-04 - - - - **** ITERATION 8 ( AFTER 11.38 ( 546) FUNCTION CALLS) - - - FUNCTION VALUE = 0.243047997834873D+01 - - NORM OF THE PROJECTED GRADIENT = 0.7055696D-05 - - - - **** ITERATION 9 ( AFTER 12.29 ( 590) FUNCTION CALLS) - - - FUNCTION VALUE = 0.243047997834545D+01 - - NORM OF THE PROJECTED GRADIENT = 0.1430553D-05 - - - - **** ITERATION 10 ( AFTER 13.21 ( 634) FUNCTION CALLS) - - - FUNCTION VALUE = 0.243047997834530D+01 - - NORM OF THE PROJECTED GRADIENT = 0.2340654D-06 - - - - **** ITERATION 11 ( AFTER 14.13 ( 678) FUNCTION CALLS) - - - FUNCTION VALUE = 0.243047997834529D+01 - - NORM OF THE PROJECTED GRADIENT = 0.1346918D-07 - - ************* EXIT OF ROUTINE VE08AD ************* - - - TERMINATION CRITERION OF USER SATISFIED (FLAG= 1) - - - - **** ITERATION -1 ( AFTER 3.96 ( 190) FUNCTION CALLS) - - - FUNCTION VALUE = 0.387020000000000D+01 - - NORM OF THE PROJECTED GRADIENT = 0.7196223D+01 - - - - **** ITERATION 7 ( AFTER 23.46 ( 1126) FUNCTION CALLS) - - - FUNCTION VALUE = 0.243047997839556D+01 - - NORM OF THE PROJECTED GRADIENT = 0.2682123D-04 - - ************* EXIT OF ROUTINE VE08AD ************* - - - TERMINATION CRITERION OF USER SATISFIED (FLAG= 1) - - - - **** ITERATION -1 ( AFTER 6.00 ( 288) FUNCTION CALLS) - - - FUNCTION VALUE = 0.387020000000000D+01 - - NORM OF THE PROJECTED GRADIENT = 0.7196225D+01 - - - - **** ITERATION 13 ( AFTER 42.08 ( 2020) FUNCTION CALLS) - - - FUNCTION VALUE = 0.243047997834529D+01 - - NORM OF THE PROJECTED GRADIENT = 0.1299014D-08 - - ************* EXIT OF ROUTINE VE08AD ************* - - - TERMINATION CRITERION OF USER SATISFIED (FLAG= 1) - - - - **** ITERATION -1 ( AFTER 2.50 ( 5) FUNCTION CALLS) - - - FUNCTION VALUE = 0.440000000000000D+02 - - NORM OF THE PROJECTED GRADIENT = 0.6403125D+01 - - STEP STRATEGY ITERATIONS 0 LINESEARCH PARAMETER 0.100D+01 - - APPROXIMATE SOLUTION X - - 0.1000000D+02 0.4000000D+01 0.1000000D+02 - - GRADIENT AT APPROXIMATE SOLUTION X - - 0.1000000D+01 0.2000001D+01 0.6000001D+01 - - ORIGINAL STARTING POINT - - 0.1000000D+02 0.4000000D+01 0.1000000D+02 - - PROJECTED 1-TH ELEMENT HESSIAN - - - PROJECTED 2-TH ELEMENT HESSIAN - - 0.10D+02 - 0.40D+01 0.19D-15 - - ************************************************** - - - - - - **** ITERATION 0 ( AFTER 3.00 ( 6) FUNCTION CALLS) - - - FUNCTION VALUE = 0.340000000000000D+02 - - NORM OF THE PROJECTED GRADIENT = 0.6324556D+01 - - STEP STRATEGY ITERATIONS 0 LINESEARCH PARAMETER 0.100D+01 - - APPROXIMATE SOLUTION X - - 0.0000000D+00 0.4000000D+01 0.1000000D+02 - - GRADIENT AT APPROXIMATE SOLUTION X - - 0.0000000D+00 0.2000001D+01 0.6000001D+01 - - LAST STEP OF THE ALGORITHM - - 0.1000000D+02 0.4000000D+01 0.1000000D+02 - - PROJECTED 1-TH ELEMENT HESSIAN - - - PROJECTED 2-TH ELEMENT HESSIAN - - 0.10D+01 - 0.00D+00 0.10D+01 - - ************************************************** - - - - - - **** ITERATION 6 ( AFTER 12.50 ( 25) FUNCTION CALLS) - - - FUNCTION VALUE = 0.273162791442526D-11 - - NORM OF THE PROJECTED GRADIENT = 0.3498870D-05 - - STEP STRATEGY ITERATIONS 10 LINESEARCH PARAMETER 0.100D+01 - - APPROXIMATE SOLUTION X - - 0.0000000D+00 0.1412197D-06 -0.2187593D-05 - - GRADIENT AT APPROXIMATE SOLUTION X - - 0.0000000D+00 0.2611275D-05 -0.2328806D-05 - - LAST STEP OF THE ALGORITHM - - 0.0000000D+00 0.1314405D-02 0.9426197D-03 - - PROJECTED 1-TH ELEMENT HESSIAN - - - PROJECTED 2-TH ELEMENT HESSIAN - - 0.30D+01 - -0.10D+01 0.10D+01 - - ************************************************** - - - - ************* EXIT OF ROUTINE VE08AD ************* - - - TERMINATION CRITERION OF USER SATISFIED (FLAG= 1) - - - - **** ITERATION -1 ( AFTER 3.00 ( 147) FUNCTION CALLS) - - - FUNCTION VALUE = 0.197960000000000D+05 - - NORM OF THE GRADIENT = 0.8389755D+04 - - - - **** ITERATION 28 ( AFTER 38.20 ( 1872) FUNCTION CALLS) - - - FUNCTION VALUE = 0.742022288973514D-29 - - NORM OF THE GRADIENT = 0.1134824D-12 - - ************* EXIT OF ROUTINE VE08AD ************* - - - ERROR RETURN . - THE ROUTINE VE08AD WAS STOPPED BECAUSE - THE NORM OF THE COMPUTED STEP WAS NEGLIGIBLE - ( LT RELPR*NORM(X) = , 0.22D-15). - THIS ERROR USUALLY OCCURS WHEN TOO MUCH - ACCURACY IS REQUIRED FOR TERMINATION. - - - - **** ITERATION -1 ( AFTER 6.40 ( 64) FUNCTION CALLS) - - - FUNCTION VALUE = 0.360000000000000D+03 - - NORM OF THE GRADIENT = 0.8147639D+03 - - - - **** ITERATION 29 ( AFTER 380.30 ( 3803) FUNCTION CALLS) - - - FUNCTION VALUE = 0.820822406926453D-11 - - NORM OF THE GRADIENT = 0.3789656D-04 - - ************* EXIT OF ROUTINE VE08AD ************* - - - TERMINATION CRITERION OF USER SATISFIED (FLAG= 1) //GO.SYSIN DD ve08ts.out