DQED.TXT 09 November 2000 Introduction DQED solves the constrained nonlinear least squares problem: Minimize the sum of squares of MEQUA generally nonlinear equations, f(1:MEQUA)(x) = 0, Equation (1) where x is a set of NVARS unknowns. The vector function with these MEQUA components is called f(x) in the discussion that follows. The components of x may have upper and lower bounds given by the user. In fact, all of the possible cases can be specified: no bounds on X; bounds at one end only; upper and lower bounds. Linear constraints on the unknowns, more general than simple bounds, can also be given. These linear constraints can be of the equality or inequality type: a(L,1) x(1)+ ... + a(L,NVARS) x(NVARS) = y(L), L = 1,...,MCON, Equation (2) with bounds specified on the right hand side values y(L), again given by the user. The constraints can actually be slightly nonlinear. In this case the constraints can be described as: g(L)(x) = y(L), L = 1,...,MCON, Equation (2') where bounds are specified on each y(L). The functions g(L)(x) must be defined for all x in the set described by the simple bounds. Experienced users may wish to turn directly to the examples before reading the subprogram documentation. There is no size relation required for the problem dimensions MEQUA, NVARS, and MCON except that MEQUA and NVARS are both positive, and MCON is nonnegative. DQED will do a decent job of solving most nonlinear least squares problems that can be expressed as Equations (1) and (2) above, provided that continuous derivatives of the functions with respect to the parameters can be computed. This can also include problems where the derivatives must be computed using some form of numerical differentiation. Numerical differentiation is not provided with this software for solving nonlinear least squares problems. The authors also plan to develop methods that will do a much better job of coping with constraints more general than the essentially linear ones indicated above in Equations (2)-(2'). There are nonlinear least squares problems with innocent-looking but highly nonlinear constraints where this package will fail to work. The authors also hope to reduce the overhead required by the software. This high overhead is due primarily to the method used to solve the inner-loop quadratic model problem. The authors recommend that users consider using the option number 14, described below, to suppress use of the quadratic model. The user may find that the software works quite well without the quadratic model. This may be important when the function and derivatives evaluations are not expensive but many individual problems are being solved. There are two fundamental ways to use the subprogram DQED. The most straightforward way is to make one call to DQED and obtain values for the unknowns, x. The user provides a subprogram DQEDEV that gives DQED the values of the functions f(x) and g(x), and the derivative or Jacobian matrices for f(x) and g(x) at each desired point x. This usage is called 'forward communication.' An alternate way to use DQED is to provide an option that allows the user to communicate these values by 'reverse communication.' DQED returns to the user calling program and requests values for f(x) and g(x), and the Jacobian matrices for f(x) and g(x) for a given value of x. This framework is often required in applications that have complicated algorithmic requirements for evaluation of the functions. Examples using both 'forward' and 'reverse' communication are provided. The source code for DQED can be obtained from NETLIB, in the OPT subdirectory. Calling Sequence Explained There are arrays used by the subprogram that must have dimensions equivalent to the following declarations. INTEGER IWORK(LIWORK) DOUBLE PRECISION WORK(LWORK) The array dimensions for the arrays IWORK and WORK can change if either option 14 or option 15 are in use. In the usual case, NPMAX=5, unless the user is using option 15, when NPMAX = user specified value, must be >= 2. In the usual case, NALL = MCON + 2*NVARS + NPMAX + 1, unless the user is using option 14, when NALL = MCON + NVARS + 2. Also, if MCON is zero, then NPLUS = 0 else NPLUS = 3*NALL + 2. In terms of these values defined above, LIWORK >= 3*MCON + 9*NVARS + 4*NPMAX + NALL + 11 LWORK >= NALL*NALL + 4*NALL + NVARS*NPMAX + 33*NVARS + MEQUA*NPMAX + max(MEQUA,NVARS)*NPMAX + 13*NPMAX + 9*MCON + 26 + NPLUS SUBROUTINE DQED(DQEDEV, MEQUA, NVARS, MCON, IND, BL, BU, X, FJ, LDFJ, FNORM, IGO, IOPT, ROPT, IWORK, WORK) DQEDEV Input, EXTERNAL DQEDEV. This is the name of a subprogram that the user will usually supply for evaluation of the values of the constraints and model, and the derivatives of these functions. The user must provide this subprogram unless 'reverse communication' is used. A model for writing the subprogram DQEDEV is provided under the heading "Example 1 Using Forward Communication". Users may find it convenient to modify this model subprogram when writing a subprogram for their own application. If 'reverse communication' is used, the user does not need to write a stub or dummy subroutine named DQEDEV. All that is required is to declare exactly this name in an EXTERNAL statement. The code package has a dummy subroutine DQEDEV that will be used in the linking or load step. "Example 2 Using Reverse Communication", listed below, illustrates this detail. MEQUA Input, INTEGER MEQUA. The number of least squares equations. MEQUA must be positive. NVARS Input, INTEGER NVARS. The number of unknowns or variables. NVARS must be positive. MCON Input, INTEGER MCON. The number of general constraints for the solution, not including simple bounds. MCON must be nonnegative. IND Input, INTEGER IND(NVARS+MCON). IND describes the form of the simple bounds that the components of x are to satisfy. Components numbered 1,...,NVARS are used to describe the form of the simple bounds that the unknown are to satisfy. Components numbered NVARS+1,...,NVARS+MCON are used to describe the form of the general MCON linear constraints. The first NVARS components of IND indicate the type of simple bounds that the solution is to satisfy. The corresponding entries of BL and BU are the bounding value. The only values of IND allowed are 1,2,3 or 4. Other values are errors. Specifically, IND(J)= 1, if x(J)>= BL(J) is required; BU(J) is not used. 2, if x(J)<= BU(J) is required; BL(J) is not used. 3, if x(J)>= BL(J) and x(J)<= BU(J) is required. 4, if no bounds on x(J) are required; BL, BU are not used. General linear constraints of the form shown in Equation (2) require that bounds be given for the linear functions y(L). Specifically, IND(NVARS+L)= 1, if y(L)>= BL(NVARS+L) is required; BU is not needed. 2, if y(L)<= BU(NVARS+L) is required; BL is not needed. 3, if y(L)>= BL(NVARS+L) and y(L)<= BU(NVARS+L) 4, if no bounds on y(L) are required; BL, BU are not used. The values of the bounds for the unknowns may be changed by the user during the evaluation of the functions f(x) and g(x) and their Jacobian matrices. BL Input, REAL BL(NVARS+MCON). For each component J, if IND(J) requires it, BL(J) is the lower value for X(J). BU Input, REAL BU(NVARS+MCON). For each component J, if IND(J) requires it, BU(J) is the upper value for X(J). X Input/output, REAL X(NVARS). The array X contains the NVARS values, x, where the functions f(x) and g(x) and their Jacobian matrices will be evaluated by the subprogram DQED. After the computation has successfully completed, the array X will contain a solution, namely the unknowns of the problem, x. Success is determined by an appropriate value for IGO. On input, the array X must contain a starting guess for the unknowns of the problem, x. The initial values do not need to satisfy the constraints or the bounds. If they do not satisfy the bounds, then the point will be simply projected onto the bounds as a first step. The first linear change to the values of x must satisfy the general constraints. It is here that the assumption of their linearity is utilized. FJ Input/output, REAL FJ(LDFJ,NVARS+1). FJ is used to store the linear constraints of Equation (2) or more generally the Jacobian matrix of the functions g(x) with respect to the variables x, and the Jacobian matrix of the function f(x). The Jacobian matrix of the (linear) constraints is placed in rows 1,...,MCON. The Jacobian matrix of f(x) is placed in rows MCON+1, ...,MCON+MEQUA. Normally the array FJ is assigned values by the user when the nonlinear solver requests evaluations of the constraints g(x) and the function f(x) together with the Jacobian matrices G(x) and J(x). The values of the constraint functions g(L)(x) are placed in the array FJ(L,NVARS+1), L=1,...,MCON. The values of the model functions f(I)(x) are placed in the array at entries FJ(MCON+I,NVARS+1), I=1,...,MEQUA. LDFJ Input, INTEGER LDFJ, the leading or row dimension of the array FJ. LDFJ must be greater than or equal to MEQUA+MCON. FNORM Output, REAL FNORM. This is the value of the Euclidean length or square root of sums of squares of components of the function f(x) after the approximate solution, x, has been found. During the computation it is updated and equals the best value of the length of f(x) that has been found. IGO Input/output, INTEGER IGO. This flag directs user action and informs the user about the type of results obtained by the subprogram. The user may find it convenient to treat the cases abs(IGO) <= 1 the same as IGO=1. This has no effect on the solution process. The user can interrupt the computation at any time, obtaining the best values of the vector x up to this point, by setting IGO to any value .gt. 1 and then return control to DQED. For example if a calculation must be done in a certain length of time, the user can, as the end of the time draws near, set IGO=20 and return control to DQED. It is important that this method be used to stop further computing, rather than simply proceeding. The reason for this is that certain flags in DQED must be reset before any further solving on subsequent problems can take place. The value of IGO .gt. 1 used to interrupt the computation is arbitrary and is the value of IGO returned. If values of IGO =2,...,18 are used to flag this interrupt, they do not mean the same thing as indicated below. For this reason the value IGO=20 is recommended for signaling interrupts in DQED. Another situation that may occur is the request for an evaluation of the functions and derivatives at a point x where these can't be evaluated. If this occurs, set IGO=99 and return control to DQED. This will have the effect of defining the derivatives to be all zero and the functions to be 'large.' Thus a reduction in the trust region around the current best estimate will occur. Assigning the value IGO=99 will not cause DQED to stop computing. Output values of IGO: 0, Place the value of f(x) in FJ(MCON+*,NVARS+1). If 'reverse communication' is being used, CALL DQED again. If 'forward communication' is being used, do a RETURN. 1 or -1, Evaluate the Jacobians for the functions g(x) and f(x) as well as evaluating g(x) and f(x). Use the vector x that is now in the array X as the values where this evaluation will be performed. Place the Jacobian matrix for g(x) in the first MCON rows of FJ. Place the Jacobian matrix for f(x) in rows MCON+1,...,MCON+MEQUA in FJ. Place the value of g(x) in FJ(*,NVARS+1). Place the value of f(x) in FJ(MCON+*,NVARS+1). Users who have complicated functions whose derivatives cannot be computed analytically may want to use the numerical differentiation subroutine JAGC. This is available on the SLATEC library. If 'reverse communication' is being used, CALL DQED again. If 'forward communication' is being used, do a RETURN. A value IGO=-1 flags that that the number of terms in the quadratic model is being restricted by the amount of storage given for that purpose. It is suggested, but it is not required, that additional storage be given for the quadratic model parameters. See the following description of the option array, option number 15, for the way to designate more storage for this purpose. 2, The function f(x) has a length less than TOLF. This is the value for IGO to be expected when an actual zero value of f(x) is anticipated. See the description of the option array for the value. 3, The function f(x) has reached a value that may be a local minimum. However, the bounds on the trust region defining the size of the step are being hit at each step. Thus the situation is suspect. (Situations of this type can occur when the solution is at infinity in some of the components of the unknowns, x. See the description of the option array for ways to avoid this value of output value of IGO. 4, The function f(x) has reached a local minimum. This is the value of IGO that is expected when a nonzero value of f(x) is anticipated. See the description of the option array for the conditions that have been satisfied. 5, The model problem solver has noted a value for the linear or quadratic model problem residual vector length that is >= the current value of the function, which is the Euclidean length of f(x). This situation probably means that the evaluation of f(x) has more uncertainty or noise than is possible to account for in the tolerances used to note a local minimum. The value for x is suspect, but a minimum has probably been found. 6, A small change (absolute) was noted for the vector x. A full model problem step was taken. The condition for IGO=4 may also be satisfied, so that a minimum has been found. However, this test is made before the test for IGO=4. 7, A small change (relative to the length of x) was noted for the vector x. A full model problem step was taken. The condition for IGO=4 may also be satisfied, so that a minimum has been found. However, this test is made before the test for IGO=4. 8, More than ITMAX iterations were taken to obtain the solution. The value obtained for x is suspect, although it is the best set of x values that occurred in the entire computation. See the description of the option array for directions on how to increase this value. Note that the nominal value for ITMAX, 75, is sufficient to solve all of the nonlinear test problems described in Reference (2). 9 through 18, Errors in the usage of the subprogram were noted. The exact condition will be noted using an error processor that prints an informative message unless this printing has been suppressed. A minimum value has not been found for x. The relation between IGO and the error number are IGO=NERR + 8. Here NERR is the identifying number. IOPT Input, INTEGER IOPT(17+*). IOPT must have a dimension of at least 17! The dimension of IOPT may be greater than 17 if many options are specified. IOPT is an array of option indices and option values, which makes it easy for the user to specify many different DQED options. In some cases, the ROPT array may also be used. See the example for details. 99, There are no more options to change. Normally this is the first and only option that a user needs to specify, and it can be simply IOPT(1)=99. The total dimension of IOPT must be at least 17, however. This can lead to a hard-to-find program bug if the dimension is too small. 1, Change the amount of printed output. The next value of IOPT is the print level desired, IPRINT. Any value of IPRINT .gt. 0 gives all the available output. 2, Change the value of ITMAX. The next value of IOPT is the value of ITMAX desired. 3, Pass prior determined bounds for the box containing the initial point. This box is the trust region for the first move from the initial point. The next entry in IOPT points to the place in ROPT where the NVARS values for the edges of the box are found. 4, Change the value of TOLF. The next entry of IOPT points to the place in ROPT where the new value of TOLF is found. 5, Change the value of TOLX. The next entry of IOPT points to the place in ROPT where the new value of TOLX is found. 6, Change the value of TOLD. The next entry of IOPT points to the place in ROPT where the new value of TOLD is found. 7, Change the value of TOLSRN. The next entry of IOPT points to the place in ROPT where the new value of TOLSNR is found. 8, Change the value of TOLP. The next entry of IOPT points to the place in ROPT where the new value of TOLP is found. 9, Change the value of TOLUSE. The next entry of IOPT points to the place in ROPT where the new value of TOLUSE is found. 10, Change the value of COND. The next entry of IOPT points to the place in ROPT where the new value of COND is found. 11, Change the value of LEVEL. The next entry of IOPT is the new value of LEVEL. 12, Pass an option array to the subprogram DQEDGN used as the inner loop solver for the model problem. The next entry of IOPT is the starting location for the option array for DQEDGN within the array IOPT. Thus the option array for DQEDGN must be a part of the array IOPT. 13, Move (or leap) the processing pointer LP for the option array by the next value in IOPT. 14, Change a logical flag that suppresses the use of the quadratic model in the inner loop. Use the next value in IOPT for this flag. If this value = 1, then never use the quadratic model. Just use the linear model. Otherwise, use the quadratic model when appropriate. This option decreases the amount of scratch storage as well as the computing overhead required by the code package. A user may want to determine if the application really requires the use of the quadratic model. If it does not, then use this option to save both storage and computing time. 15, Change NTERMS, the maximum number of array columns that can be used for saving quadratic model data. The value of NTERMS is one more than the maximum number of terms used. Each unit increase for NTERMS increases the required dimension of the array WORK by 2*MEQUA+NVARS. Use the value in IOPT(LP+1) for the new value of NTERMS. Decreasing this value to 2 (its minimum) decreases the amount of storage required by the code package. 16, Change a logical flag so that 'reverse communication' is used instead of 'forward communication.' Example EX01, listed below, uses 'forward communication.' Example EX02, also listed below, uses 'reverse communication.' Use the next value in IOPT for this flag. If this value = 1, then use 'reverse communication.' Otherwise, use 'forward communication.' WARNING: This usage may not work unless the operating system saves variables between subroutine calls to DQED. 17, Do not allow the flag IGO to return with the value IGO=3. This means that convergence will not be claimed unless a full model step is taken. Normal output values will then be IGO = 2, 4, 6 or 7. Use the next value in IOPT for this flag. If this value = 1, then force a full model step. Otherwise, do not force a full model step if small steps are noted. ROPT Input, REAL ROPT(*). Depending on the options specified by the user in IOPT, ROPT may contain real numbers required to set certain DQED options. IWORK Input/output/workspace, INTEGER IWORK(*). IWORK is an integer scratch array that DQED uses for storage of intermediate results. It is important not to modify the contents of this storage during the computation. The array locations IWORK(1) and IWORK(2) must contain the actual lengths of the arrays WORK and IWORK before the call to DQED. These array entries are replaced by the actual amount of storage required for each array. If the amount of storage for either array is too small, an informative error message will be printed, and the value IGO=13 or 14 returned. The user may find it useful to let the subprogram DQED return the amounts of storage required for these arrays. For example set IWORK(1)=1, IWORK(2)=1. The subprogram will return with IGO=13, IWORK(1)=required length of WORK, and IWORK(2)=required length of IWORK. Appropriate error messages will normally be printed. WORK Workspace, REAL WORK(*). WORK is a real scratch array that DQED uses for storage of intermediate results. It is important not to modify the contents of this storage during the computation. The array location IWORK(1) must contain the actual length of WORK before the call to DQED. This array entry is replaced by the actual amount of storage required. If the amount of storage is too small, an informative error message will be printed, and the value IGO=13 or 14 returned. The Option Array In order to use the option array technique to change selected data within a subprogram, it is necessary to understand how this array is processed within the software. Let LP designate the processing pointer that moves to positions of the IOPT array. Initially LP=1, and as each option is noted and changed, the value of LP is updated. The values of IOPT(LP) determine what options get changed. The amount that LP changes is known by the software to be equal to the value two except for two options. These exceptional cases are the last option (=99) and the 'leap' option (=13) which advances LP by the value in IOPT(LP+1). A negative value for IOPT(LP) means that this option is not to be changed. This aids the programmer in using options; often the code for using an option can be in the calling program but a negative value of the option number avoids rewriting code. Glossary of Items Modified by Options. Those items with Nominal Values listed can be changed. The quantity 'eps', used below, is the machine precision parameter. Its value is obtained by a call to the Bell Labs. Port subprogram D1MACH(4). It is machine dependent. The conditions (abs(FB-PV)<=TOLSNR*FB and abs(FC-PV) <= TOLP*FB) and (ABS(FC-FL)<=TOLSNR*FB) together with taking a full model step, must be satisfied before the condition IGO=4 is returned. Decreasing any of the values for TOLF, TOLD, TOLX, TOLSNR, or TOLP will likely increase the number of iterations required for convergence. Names Nominal Definition Values ----- ------- ---------- FC Current value of length of f(x). FB Best value of length of f(x). FL Value of length of f(x) at the previous step. PV Predicted value of length of f(x), after the step is taken, using the approximating model. MIN(1.D-5, TOLF sqrt(eps)) Tolerance for stopping when FC <= TOLF. MIN(1.D-5, TOLD sqrt(eps)) Tolerance for stopping when change to x values has length <= TOLD. MIN(1.D-5, TOLX sqrt(eps)) Tolerance for stopping when change to x values has length <= TOLX*length of x values. TOLSNR 1.D-5 Tolerance used in stopping condition IGO=4. Explained below. TOLP 1.D-5 Tolerance used in stopping condition IGO=4. Explained below. COND 30. Largest condition number to allow when solving for the quadratic model coefficients. Increasing this value may result in more terms being used in the quadratic model. TOLUSE sqrt(eps) A tolerance that is used to avoid values of x in the quadratic model's interpolation of previous points. Decreasing this value may result in more terms being used in the quadratic model. ITMAX 75 The number of iterations to take with the algorithm before giving up and noting it with the value IGO=8. IPRINT 0 Control the level of printed output in the solver. A value of IPRINT .gt. 0 will result in output of information about each iteration. The output unit used is obtained using the Bell Labs. Port subprogram, i. e. I1MACH(2). LEVEL 1 Error processor error level. See the SLATEC library documentation for XERROR for an explanation. NTERMS 5 One more than the maximum number of terms used in the quadratic model. Option Usage Example In the Fortran code fragment that follows, an example is given where we change the value of TOLF and decrease the maximum number of iterations allowed from 75 to 30. In this example the dimensions of IOPT and ROPT must satisfy: double precision ropt(1) integer iopt(5) . . . c c Set the option to change the value of TOLF. c iopt(1)=4 c c The next entry points to the place in ROPT where c the new value of TOLF is located. c iopt(2)=1 c c This is the new value of TOLF. c ropt(1)=1.d-9 c c Change the number of iterations. c iopt(3)=2 c c This next entry is the new value for the maximum number of iterations. c iopt(4)=30 c c This next option is a signal that there are no more options. c iopt(5)=99 . . . call dqed . . . Examples The following complete program units, EX01 and EX02, show how one can use the nonlinear solver for fitting exponential functions to given data. These examples are calculations that match two terms of an exponential series to five given data points. There are some subtle points about exponential fitting that are important to note. First, the signs of the exponential arguments are restricted to be nonpositive. The size of the arguments should not be much larger than the start of the time data (reciprocated). This is the reason the lower bounds are set a bit less than the reciprocal of the time value. In many applications that require exponential modeling this is a natural assumption. The nonlinear solver allows these bounds on the arguments explicitly. In addition, the coefficients are constrained to be nonnegative. These bounds are harder to justify. The idea is to avoid the situation where a coefficient is very large and negative, and the corresponding exponential argument is also large and negative. The resulting contribution to the series may be very small, but its presence is spurious. Finally, the single general linear constraint that keeps the arguments separated (by 0.05 in this example) is used for two purposes. First, it naturally orders these values so that the first one is algebraically largest. Second, this constraint moves the parameters from the local minimum corresponding to the initial values used in the examples. This constraint also retains the validity of the model function h(t) = w*exp(x*t) + y*exp(z*t). Namely, if the arguments are allowed to coalesce to the same value, then the model itself must change. The form of the model can become h(t)=(a+b*t)*exp(c*t) or h(t) = d*exp(e*t). Either one could occur, and the choice is problem dependent. c EX1.F 25 January 1996 c program ex1 c c*********************************************************************** c c EX1 illustrates the use of the Hanson-Krogh nonlinear least squares c solver DQED for fitting two exponentials to data. c c The problem is to find values for the four variables x(1),...,x(4), c which specify the model function c c h(t) = x(1)*exp(x(2)*t) + x(3)*exp(x(4)*t) c c which mininize the sum of the squares of the error between h(t) and c recorded values of h at five values of t: c c t=0.05, 0.1, 0.4, 0.5, and 1.0. c c We also have problem constraints that c c 0 <= x(1) c -25.0 <= x(2) <= 0 c 0 <= x(3) c -25.0 <= x(4) <= 0. c c and a minimal separation requirement between x(2) and x(4) that c can be expressed as c c 0.05 <= x(2) - x(4) c c c Output from Example 1 Program c c model is h(t) = x(1)*exp(-t*x(2)) + x(3)*exp(t*x(4)) c x(1),x(2),x(3),x(4) = c 1.999475 -.999801 .500057 -9.953988 c residual after the fit = 4.2408d-04 c output flag from solver = 4 c integer ldfj integer liwork integer lwork integer mcon integer mequa integer nvars c parameter (liwork=84) parameter (lwork=640) parameter (mcon=1) parameter (mequa=5) parameter (nvars=4) c parameter (ldfj=mcon+mequa) c double precision bl(nvars+mcon) double precision bu(nvars+mcon) double precision fj(ldfj,nvars+1) double precision fnorm integer i integer igo integer ind(nvars+mcon) integer iopt(24) integer iwork(liwork) double precision ropt(1) double precision work(lwork) double precision x(nvars) c external dqedev c c Define the bounding constraints on the variables X(1) through X(4), c and then the linear constraints on the separation condition. c c 1: 0.0 <= X(1) c ind(1)=1 bl(1)=0.0 bu(1)=0.0 c c 2: -25.0 <= X(2) <= 0.0 c ind(2)=3 bl(2)=-25.0 bu(2)=0.0 c c 3: 0.0 <= X(3) c ind(3)=1 bl(3)=0.0 bu(3)=0.0 c c 4: -25.0 <= X(4) <= 0.0 c ind(4)=3 bl(4)=-25.0 bu(4)=0.0 c c 5: 0.05 <= X(2) - X(4) c ind(5)=1 bl(5)=0.05 bu(5)=0.0 c c Define the initial values of the variables. c We don't know anything more, so all variables are set zero. c do i = 1,nvars x(i) = 0.0 enddo c c Tell how much storage we gave the solver. c iwork(1) = lwork iwork(2) = liwork c c No additional options are in use. c iopt(1) = 99 c c Call the program. c call dqed(dqedev,mequa,nvars,mcon,ind,bl,bu,x,fj,ldfj,fnorm,igo, & iopt,ropt,iwork,work) write(*,*)' ' write(*,*)'h(t) = x(1)*exp(-t*x(2)) + x(3)*exp(t*x(4))' write(*,*)' ' write(*,*)'Constrained minimizing X:' write(*,*)' ' write(*,*)(x(i),i=1,nvars) write(*,*)' ' write(*,*)'Residual of the fitted data, FNORM = ',fnorm write(*,*)'DQED output flag IGO = ',igo stop end subroutine dqedev(x,fj,ldfj,igo,iopt,ropt) c c*********************************************************************** c c DQEDEV is the subprogram for evaluating the functions and derivatives c for the nonlinear solver DQED. c c The user problem has MCON constraint functions, c MEQUA least squares equations, and involves NVARS c unknown variables. c c When this subprogram is entered, the general (near) c linear constraint partial derivatives, the derivatives c for the least squares equations, and the associated c function values are placed into the array FJ(*,*). c All partials and functions are evaluated at the point c in X(*). Then the subprogram returns to the calling c program unit. Typically one could do the following c steps: c c step 1. Place the partials of the i-th constraint c function with respect to variable j in the c array FJ(i,j), i=1,...,MCON, j=1,...,NVARS. c c step 2. Place the values of the i-th constraint c equation into FJ(i,NVARS+1). c c step 3. Place the partials of the i-th least squares c equation with respect to variable j in the c array FJ(MCON+i,j), i=1,...,MEQUA, c j=1,...,NVARS. c c step 4. Place the value of the i-th least squares c equation into FJ(MCON+i,NVARS+1). c c step 5. Return to the calling program unit. c integer ldfj c integer mcon integer mequa integer nvars c parameter (mcon=1) parameter (mequa=5) parameter (nvars=4) c double precision f(mequa) double precision fj(ldfj,nvars+1) integer i integer igo integer iopt(*) double precision ropt(*) double precision t(mequa) double precision x(nvars) c data t/0.05,0.10,0.40,0.50,1.00/ data f/2.206d+00,1.994d+00,1.350d+00,1.216d+00,.7358d0/ c c Set the value of the constraint, and the residual functions. c fj(1,5) = x(2) - x(4) do i = 1,mequa fj(mcon+i,5) = x(1)*exp(x(2)*t(i)) + x(3)*exp(x(4)*t(i)) - f(i) enddo c c If IGO is nonzero, compute the derivatives. c if(igo.ne.0)then fj(1,1) = 0.0 fj(1,2) = 1.0 fj(1,3) = 0.0 fj(1,4) = -1.0 do i = 1,mequa fj(mcon+i,1) = exp(x(2)*t(i)) fj(mcon+i,2) = x(1)*t(i)*exp(x(2)*t(i)) fj(mcon+i,3) = exp(x(4)*t(i)) fj(mcon+i,4) = x(3)*t(i)*exp(x(4)*t(i)) enddo endif return end c EX2.F 25 January 1996 c program ex2 c c*********************************************************************** c c EX2 illustrates the use of the Hanson-Krogh nonlinear least squares c solver DQED for fitting two exponentials to data. c c The problem is to find values for the four variables x(1),...,x(4), c which specify the model function c c h(t) = x(1)*exp(x(2)*t) + x(3)*exp(x(4)*t) c c which mininize the sum of the squares of the error between h(t) and c recorded values of h at five values of t: c c t=0.05, 0.1, 0.4, 0.5, and 1.0. c c We also have problem constraints that c c 0 <= x(1) c -25.0 <= x(2) <= 0 c 0 <= x(3) c -25.0 <= x(4) <= 0. c c and a minimal separation requirement between x(2) and x(4) that c can be expressed as c c 0.05 <= x(2) - x(4) c c c Output from Example 2 Program c c model is h(t) = x(1)*exp(-t*x(2)) + x(3)*exp(t*x(4)) c x(1),x(2),x(3),x(4) = c 1.999475 -.999801 .500057 -9.953988 c residual after the fit = 4.2408d-04 c number of evaluations of parameter model = 14 c output flag from solver = 4 c integer ldfj integer liwork integer lwork integer mcon integer mequa integer nvars c parameter (liwork=84) parameter (lwork=640) parameter (mcon=1) parameter (mequa=5) parameter (nvars=4) c parameter (ldfj=mcon+mequa) c double precision bl(nvars+mcon) double precision bu(nvars+mcon) double precision f(mequa) double precision fj(ldfj,nvars+1) double precision fnorm integer i integer igo integer ind(nvars+mcon) integer iopt(24) integer iwork(liwork) integer niters double precision ropt(1) double precision t(mequa) double precision work(lwork) double precision x(nvars) c external dqedev c data t/0.05,0.10,0.40,0.50,1.00/ data f/2.206d+00,1.994d+00,1.350d+00,1.216d+00,.7358d0/ c c Define the bounding constraints on the variables X(1) through X(4), c and then the linear constraints on the separation condition. c c 1: 0.0 <= X(1) c ind(1)=1 bl(1)=0.0 bu(1)=0.0 c c 2: -25.0 <= X(2) <= 0.0 c ind(2)=3 bl(2)=-25.0 bu(2)=0.0 c c 3: 0.0 <= X(3) c ind(3)=1 bl(3)=0.0 bu(3)=0.0 c c 4: -25.0 <= X(4) <= 0.0 c ind(4)=3 bl(4)=-25.0 bu(4)=0.0 c c 5: 0.05 <= X(2) - X(4) c ind(5)=1 bl(5)=0.05 bu(5)=0.0 c c Define the initial values of the variables. c We don't know anything at all, so all variables are set zero. c do i=1,nvars x(i) = 0.0 enddo c c Tell how much storage we gave the solver. c iwork(1) = lwork iwork(2) = liwork c c Initialize the call counter. c niters = 0 c c Use reverse commumication to evaluate the derivatives. c iopt(1)=16 iopt(2)=1 c c No more options. c iopt(3) = 99 10 continue call dqed(dqedev,mequa,nvars,mcon,ind,bl,bu,x,fj,ldfj,fnorm, & igo,iopt,ropt,iwork,work) if (igo.gt.1)then write(*,*)' ' write(*,*)'h(t) = x(1)*exp(-t*x(2)) + x(3)*exp(t*x(4))' write(*,*)' ' write(*,*)'Constrained minimizing X:' write(*,*)' ' write(*,*)(x(i),i=1,nvars) write(*,*)' ' write(*,*)'Residual of the fitted data, FNORM = ',fnorm write(*,*)'Number of model evaluations NITERS = ',niters write(*,*)'DQED output flag IGO = ',igo stop endif c c Count function evaluations. c niters = niters + 1 c c Set the value of the constraint, and the residual functions. c fj(1,5) = x(2) - x(4) do i = 1,mequa fj(mcon+i,5) = x(1)*exp(x(2)*t(i)) + x(3)*exp(x(4)*t(i)) - f(i) enddo c c If IGO is nonzero, compute the derivatives. c if(igo.ne.0)then fj(1,1) = 0.0 fj(1,2) = 1.0 fj(1,3) = 0.0 fj(1,4) = -1.0 do i = 1,mequa fj(mcon+i,1) = exp(x(2)*t(i)) fj(mcon+i,2) = x(1)*t(i)*exp(x(2)*t(i)) fj(mcon+i,3) = exp(x(4)*t(i)) fj(mcon+i,4) = x(3)*t(i)*exp(x(4)*t(i)) enddo endif go to 10 end References J J Dongarra, J R Bunch, C B Moler, G W Stewart, LINPACK User's Guide, Society for Industrial and Applied Mathematics, Philadelphia, PA, 1979. R J Hanson, "Least Squares with Bounds and Linear Constraints," SIAM Journal of Scientific and Statistical Computing, Volume 7, number 3, July 1986, pages 826-834. R B Schnabel, P D Frank, "Tensor Methods for Nonlinear Equations," SIAM Journal of Numerical Analysis, Volume 21, number 5, October 1984, pages 815-843. subroutine dpchek(df,dqedev,fj,iopt,ldfj,nvars,ropt,x,y) DPCHEK carries out a check of the user's jacobian evaluation routine, by comparing each column with a finite difference estimate. DF Workspace, DOUBLE PRECISION DF(NVARS). DQEDEV Input, EXTERNAL DQEDEV, the name of the user written jacobian and function evaluation routine. FJ Workspace, DOUBLE PRECISION FJ(LDFJ,NVARS+1), space to store the jacobian and the function, as required by DQEDEV. IOPT Throughput, INTEGER IOPT(*), parameters to be passed to DQEDEV. LDFJ Input, INTEGER LDFJ, the leading dimension of FJ, which must be at least NVARS. NVARS Input, INTEGER NVARS, the number of variables. ROPT Throughput, DOUBLE PRECISION ROPT(*), parameters to be passed to DQEDEV. X Input, DOUBLE PRECISION X(NVARS), the point at which the jacobian should be evaluated. Y Workspace, DOUBLE PRECISION Y(NVARS). subroutine difcen(fj,func,fx,iopt,ldfj,mcon,mequa,nvars,ropt,x) DIFCEN estimates the jacobian of a set of NEQN nonlinear functions with respect to a set of NVARS variables, using central finite differences, requiring 2*(MCON+MEQUA) function evaluations. FJ Output, DOUBLE PRECISION FJ(LDFJ,NVARS), the estimated jacobian. FUNC Input, EXTERNAL FUNC, the name of the user written function evaluation routine. FUNC should have the form: subroutine func(fx,iopt,mcon,mequa,nvars,ropt,x) and should accept X as input, and return in FX the value of the MEQUA+MCON functions. FX Workspace, DOUBLE PRECISION FX(MEQUA+MCON). IOPT Throughput, INTEGER IOPT(*), parameters to be passed to FUNC. LDFJ Input, INTEGER LDFJ, the leading dimension of FJ, which must be at least MEQUA+MCON. MCON Input, INTEGER MCON, the number of constraints. MEQUA Input, INTEGER MEQUA, the number of nonlinear functions. NVARS Input, INTEGER NVARS, the number of variables. ROPT Throughput, DOUBLE PRECISION ROPT(*), parameters to be passed to FUNC. X Input, DOUBLE PRECISION X(NVARS), the point at which the jacobian should be evaluated. subroutine diffor(fj,func,fx,iopt,ldfj,mcon,mequa,nvars,ropt,x) DIFFOR estimates the jacobian of a set of NEQN nonlinear functions with respect to a set of NVARS variables, using forward finite differences, requiring MCON+MEQUA+1 function evaluations. FJ Output, DOUBLE PRECISION FJ(LDFJ,NVARS), the estimated jacobian. FUNC Input, EXTERNAL FUNC, the name of the user written function evaluation routine. FUNC should have the form: subroutine func(fx,iopt,mcon,mequa,nvars,ropt,x) and should accept X as input, and return in FX the value of the MEQUA+MCON functions. FX Workspace, DOUBLE PRECISION FX(MEQUA+MCON). IOPT Throughput, INTEGER IOPT(*), parameters to be passed to FUNC. LDFJ Input, INTEGER LDFJ, the leading dimension of FJ, which must be at least MEQUA+MCON. MCON Input, INTEGER MCON, the number of constraints. MEQUA Input, INTEGER MEQUA, the number of nonlinear functions. NVARS Input, INTEGER NVARS, the number of variables. ROPT Throughput, DOUBLE PRECISION ROPT(*), parameters to be passed to FUNC. X Input, DOUBLE PRECISION X(NVARS), the point at which the jacobian should be evaluated. subroutine dqedev(x,fj,ldfj,igo,iopt,ropt) DQEDEV evaluates certain functions being treated by DQED, as well as their partial derivatives. The user has NVARS variables X(I), and is trying to minimize the square root of the sum of the squares of MEQUA functions F(I)(X), subject to MCON constraints which have the form BL(I) <= G(I)(X) <= BU(I) where either the left or right bounding inequality may be dropped. X Input, REAL X(*), an array of length NVARS, containing the values of the independent variables at which the functions and partial derivatives should be evaluated. FJ Output, REAL FJ(LDFJ,NVARS+1). If IGO is nonzero, then partial derivatives must be placed in the first NVARS columns of FJ, as follows: Rows I = 1 to MCON, and columns J = 1 to NVARS should contain the partials of the I-th constraint function G(I) with respect to the J-th variable. Rows I=MCON+1 to MCON+MEQUA, and columns J = 1 to NVARS should contain the partials of the (I-MCON)-th nonlinear function F(I-MCON) with respect to the J-th variable. Regardless of the value of IGO, column NVARS+1 must be assigned the values of the constraints and functions, as follows: Rows I = 1 to MCON, column J = NVARS+1, should contain the value of G(I). Rows I=MCON+1 to MCON+MEQUA, column J = NVARS+1, should contain the value of F(I-MCON). LDFJ Input, INTEGER LDFJ, the leading dimension of FJ, which must be at least MCON+MEQUA. IGO Input/output, INTEGER IGO. On input, IGO tells the user whether the partial derivatives are needed. 0, the partial derivatives are not needed. nonzero, the partial derivatives are needed. On output, the user may reset the input value of IGO if one of two situations is encountered: 99, the functions, constraints, or partial derivatives could not be evaluated at the given input point X. Request that DQED reject that point, and try a different one. Any other value, abort the run. IOPT Input, INTEGER IOPT(*), the integer option array. ROPT Input, REAL ROPT(*), the real option array.