PROGRAM HEART C REVISED 890628-0930, 950228-1300 C REVISED YYMMDD-HHMM C ILLUSTRATE THE USE OF THE DQED, HANSON-KROGH NONLINEAR LEAST C SQUARES SOLVER: SOLVING A CERTAIN HEART DIPOLE MOMENT EQUATION. C REFERENCE: A NEW NONLINEAR EQUATIONS TEST PROBLEM, RICE UNIV. C DEPT. OF MATH. REPORT 83-16, 6/83, 6/85, BY J. DENNIS, C JR., DAVID GAY, AND PHUONG AHN VU. C C COMMENTS ON THE MODULES DISTRIBUTED WITH THIS PACKAGE. C HEART - A TEST PROGRAM THAT READS AN INPUT FILE, CALLS C THE NONLINEAR TENSOR-QUADRATIC MODEL SOLVER, AND C WRITES THE OUTPUT FILE. OUTPUT RESULTS WILL VARY C IN ENVIRONMENTS DIFFERENT FROM THE PC/AT WITH THE C LAHEY F77L FORTRAN 77 COMPILER, v. 3.0. C DQED, DQEDMN, DQEDGN, DQEDIP, DQEDEV - C THE NONLINEAR SOLVING SOFTWARE. C DBOCLS, DBOLS, DBOLSM - C LINEARLY CONSTRAINED LINEAR LEAST SQUARES SOLVER. C DMOUT, DVOUT, IVOUT - C MATRIX-VECTOR PRINTING SUBPROGRAMS. C XERROR, XERRWV, CHRCNT - C A SMALL AND PARTIAL IMPLEMENTATION OF THE SLATEC C LIBRARY ERROR PROCESSING SUBPROGRAMS. DISCARD THESE C IF YOUR LOCAL LIBRARY HAS THIS PACKAGE. C D1MACH, I1MACH - C THE BELL LABS./PORT LIBRARY ENVIRONMENT INQUIRY C SUBPROGRAMS. DISCARD THESE IF YOUR LOCAL LIBRARY C HAS THIS PACKAGE. THE CURRENT VERSION IS SET FOR C THE LAHEY F77L FORTRAN 77 COMPILER ON THE PC/AT. C YOU MUST MODIFY THESE SUBROUTINES FOR YOUR OWN SYSTEM C IF THEY ARE NOT AVAILABLE IN YOUR LOCAL LIBRARY. C DGECO, DGEFA, DGESL - LINPACK LINEAR EQUATION SOLVING SOFTWARE. C (NOT INCLUDED IN THIS PACKAGE.) C DCOPY, DASUM, DDOT, DROTG, DAXPY, IDAMAX, DROT, DSWAP, DNRM2, C DSCAL - C BLAS, BASIC LINEAR ALGEBRA SUBPROGRAMS. C (NOT INCLUDED IN THIS PACKAGE.) C C EDIT on 950228-1300: ***LINE BELOW IS THE START OF THE INPUT FILE. Make this *** a separate file and remove the * in column 1 for all lines. *EXP. 0121C *-.807 -.021 -2.379 -3.64 -10.541 -1.961 -51.551 21.053 *-.074 -.733 0.013 -.034 -3.632 3.632 -.289 .289 *EXP. 0121B *-.809 -.021 -2.04 -.614 -6.903 -2.934 -26.328 18.639 *-.056 -.753 .026 -.047 -2.991 2.991 -.568 .568 *EXP. 0121A *-.816 -.017 -1.826 -.754 -4.839 -3.259 -14.023 15.467 *-.041 -.775 0.03 -.047 -2.565 2.565 -.754 .754 *EXP. 791226 *-.69 -.044 -1.57 -1.31 -2.65 2.0 -12.6 9.48 *-.3 -.39 .3 -.344 -1.2 2.69 1.59 -1.5 *EXP. 791129 *.485 -.0019 -.0581 .015 .105 .0406 .167 -.399 *.299 .186 -.0273 .0254 -.474 .474 -.0892 .0892 ***LINE ABOVE IS END OF THE INPUT FILE. C ***LINE BELOW IS THE START OF THE OUTPUT FILE. * PROBLEM TITLE IS = EXP. 0121C * HEART DIPOLE MODEL, * X(1),...,X(8) = * -0.692099 -0.114901 0.710964 -0.731964 * -0.815740 3.426596 1.184719 -2.332848 * RESIDUAL AFTER THE FIT = 1.5561E-10 * OUTPUT FLAG FROM SOLVER = 2 * NUMBER OF FUNCTION, DERIVATIVE EVALUATIONS = 19 * * * PROBLEM TITLE IS = EXP. 0121B * HEART DIPOLE MODEL, * X(1),...,X(8) = * 0.009035 -0.818035 -0.000445 -0.020555 * 2.773429 2.529477 -14.800972 0.522047 * RESIDUAL AFTER THE FIT = 7.7707E-08 * OUTPUT FLAG FROM SOLVER = 6 * NUMBER OF FUNCTION, DERIVATIVE EVALUATIONS = 19 * * * PROBLEM TITLE IS = EXP. 0121A * HEART DIPOLE MODEL, * X(1),...,X(8) = * 0.003100 -0.819100 -0.000224 -0.016776 * 2.681514 2.250216 -20.241705 0.797098 * RESIDUAL AFTER THE FIT = 9.6877E-13 * OUTPUT FLAG FROM SOLVER = 2 * NUMBER OF FUNCTION, DERIVATIVE EVALUATIONS = 20 * * * PROBLEM TITLE IS = EXP. 791226 * HEART DIPOLE MODEL, * X(1),...,X(8) = * -0.311627 -0.378373 0.328244 -0.372244 * -1.282227 2.494300 1.554866 -1.384638 * RESIDUAL AFTER THE FIT = 3.3973E-08 * OUTPUT FLAG FROM SOLVER = 6 * NUMBER OF FUNCTION, DERIVATIVE EVALUATIONS = 4 * * * PROBLEM TITLE IS = EXP. 791129 * HEART DIPOLE MODEL, * X(1),...,X(8) = * 0.491321 -0.006321 0.000098 -0.001998 * -0.100315 0.122657 -0.020718 -4.023518 * RESIDUAL AFTER THE FIT = 1.0066E-13 * OUTPUT FLAG FROM SOLVER = 2 * NUMBER OF FUNCTION, DERIVATIVE EVALUATIONS = 18 ***LINE ABOVE IS THE END OF OUTPUT FILE. C C DIMENSION FOR THE NONLINEAR SOLVER. DOUBLE PRECISION FJ(08,9),BL(10),BU(10),X(8),ROPT(001),WA(2000) DOUBLE PRECISION SIGMA(8),Y(08),DF(08),T,D1MACH,RNORM INTEGER IND(10),IOPT(28),IWA(150) INTEGER LDFJ,LWA,LIWA,NVARS,NIN,I,J,NFEV,MCON,MEQUA,NOUT INTEGER IGO,K,MODE CHARACTER*12 TITLE LOGICAL WANTPR COMMON /SIGMA/SIGMA COMMON /MODE/MODE,NFEV EXTERNAL DQEDHD DATA LDFJ,LWA,LIWA/08,2000,150/ C C TO CHECK PARTIAL DERIVATIVES NUMERICALLY SET WANTPR=.TRUE. WANTPR = .FALSE. NVARS = 8 10 CONTINUE C C READ THE PROBLEM TITLE. NIN = 5 READ (NIN,'(A)',END=60) TITLE C C READ IN THE SUMS OF VALUES AND INITIAL PARAMETER VALUES. READ (NIN,*,END=60) SIGMA,X C C TEST THE PARTIAL DERIVATIVE COMPUTATION. IF (WANTPR) THEN T = SQRT(D1MACH(4)) DO 30 J = 1,NVARS CALL DCOPY(NVARS,X,1,Y,1) Y(J) = X(J) + T*X(J) CALL DQEDHD(Y,FJ,LDFJ,IGO,IOPT,ROPT) CALL DCOPY(NVARS,FJ(1,NVARS+1),1,DF,1) Y(J) = X(J) - T*X(J) CALL DQEDHD(Y,FJ,LDFJ,IGO,IOPT,ROPT) DO 20 I = 1,NVARS DF(I) = (DF(I)-FJ(I,NVARS+1))/ (2.*T*X(J)) 20 CONTINUE CALL DQEDHD(X,FJ,LDFJ,IGO,IOPT,ROPT) CALL DVOUT(NVARS,DF,'(///,'' NUMERICAL DERIVATIVE'')',-4) CALL IVOUT(1,J,'('' = VARIABLE NUMBER'')',-4) CALL DVOUT(NVARS,FJ(1,J),'('' ANALYTIC PARTIAL'')',-4) 30 CONTINUE END IF C C SOLVE THE PROBLEM WITH THE FIRST TWO EQUATIONS AS CONSTRAINTS. DO 50 K = 1,1 C C THE LINEAR EQUATIONS WILL NOT BE CONSTRAINTS IF MODE=0. C CONVERGENCE IS NORMALLY FASTER IF THE LINEAR EQUATIONS ARE USED C AS CONSTRAINTS, MODE=1. MODE = K NFEV = 0 C TELL HOW MUCH STORAGE THE SOLVER HAS. IWA(1) = LWA IWA(2) = LIWA C C SET UP PRINT OPTION. (NOT USED HERE). IOPT(1) = -1 IOPT(2) = 1 C C SET UP FOR LINEAR MODEL WITHOUT ANY QUADRATIC TERMS. (NOT USED). IOPT(3) = -14 IOPT(4) = 1 C C DO NOT ALLOW CONVERGENCE TO BE CLAIMED ON SMALL STEPS. IOPT(5) = 17 IOPT(6) = 1 C C ALLOW UP TO NVARS = 8 QUADRATIC MODEL TERMS. IOPT(7)=15 IOPT(8)=NVARS C C CHANGE CONDITION NUMBER FOR QUADRATIC MODEL DEGREE. IOPT(9)=10 IOPT(10)=1 ROPT(01)=1.D4 C C NO MORE OPTIONS. IOPT(11) = 99 IF (MODE.EQ.1) THEN C SOLVE WITH TWO CONSTRAINT EQUATIONS. MCON = 2 IND(NVARS+1) = 3 IND(NVARS+2) = 3 BL(NVARS+1) = SIGMA(1) BU(NVARS+1) = SIGMA(1) BL(NVARS+2) = SIGMA(2) BU(NVARS+2) = SIGMA(2) ELSE C SOLVE WITH ALL REGRESSION EQUATIONS. MCON = 0 END IF MEQUA = NVARS - MCON C C ALL VARIABLES ARE OTHERWISE FREE. DO 40 I = 1,NVARS IND(I) = 4 40 CONTINUE CALL DQED(DQEDHD,MEQUA,NVARS,MCON,IND,BL,BU,X,FJ,LDFJ,RNORM, . IGO,IOPT,ROPT,IWA,WA) NOUT = 6 WRITE (NOUT,9031) TITLE WRITE (NOUT,9001) (X(J),J=1,NVARS) WRITE (NOUT,9011) RNORM WRITE (NOUT,9021) IGO,NFEV 50 CONTINUE GO TO 10 60 STOP 9001 FORMAT (' HEART DIPOLE MODEL,',/,' X(1),...,X(8) = ',/,4F12.6,/, . 4F12.6) 9011 FORMAT (' RESIDUAL AFTER THE FIT = ',1PE12.4) 9021 FORMAT (' OUTPUT FLAG FROM SOLVER =',17X,I6,/, . ' NUMBER OF FUNCTION, DERIVATIVE EVALUATIONS =',08X,I6,//) 9031 FORMAT (' PROBLEM TITLE IS = ',A) END SUBROUTINE DQEDHD(X,FJ,LDFJ,IGO,IOPT,ROPT) C REVISED 860113-0930 C REVISED YYMMDD-HHMM C THIS IS THE SUBPROGRAM FOR EVALUATING THE FUNCTIONS C AND DERIVATIVES 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 IF(IGO.NE.0) THEN C STEP 1. PLACE THE PARTIALS OF THE I-TH CONSTRAINT C FUNCTION WITH REPECT TO VARIABLE J IN THE C ARRAY FJ(I,J), I=1,...,MCON, J=1,...,NVARS. C END IF C STEP 2. PLACE THE VALUES OF THE I-TH CONSTRAINT C EQUATION INTO FJ(I,NVARS+1). C IF(IGO.NE.0) THEN C STEP 3. PLACE THE PARTIALS OF THE I-TH LEAST SQUARES C EQUATION WITH RESPECT TO VARIABLE J IN THE C ARRAY FJ(I,J), I=1,...,MEQUA, C J=1,...,NVARS. C END IF C STEP 4. PLACE THE VALUE OF THE I-TH LEAST SQUARES C EQUATION INTO FJ(I,NVARS+1). C STEP 5. RETURN TO THE CALLING PROGRAM UNIT. INTEGER LDFJ,NVARS,I,J,NFEV,MODE DOUBLE PRECISION FJ(LDFJ,*),X(*),ROPT(*),ONE,TWO,THREE,ZERO DOUBLE PRECISION SIGMA(08) COMMON /SIGMA/SIGMA COMMON /MODE/MODE,NFEV INTEGER IGO,IOPT(*) DATA ONE,TWO,THREE/1.D0,2.D0,3.D0/,ZERO/0.D0/ C C HEART DIPOLE EQUATIONS. FIRST TWO EQUATIONS MAY BE CONSTRAINTS. C C DEFINE THE CONSTRAINTS AND THEIR DERIVATIVES. NVARS = 8 DO 20 J = 1,NVARS DO 10 I = 1,2 FJ(I,J) = ZERO 10 CONTINUE 20 CONTINUE FJ(1,1) = ONE FJ(1,2) = ONE FJ(1,NVARS+1) = X(1) + X(2) FJ(2,3) = ONE FJ(2,4) = ONE FJ(2,NVARS+1) = X(3) + X(4) C C DEFINE THE REMAINING NONLINEAR EQUATIONS AND THEIR DERIVATIVES. C (START EQUATION ONE). FJ(3,1) = X(5) FJ(3,2) = X(6) FJ(3,3) = -X(7) FJ(3,4) = -X(8) FJ(3,5) = X(1) FJ(3,6) = X(2) FJ(3,7) = -X(3) FJ(3,8) = -X(4) FJ(3,NVARS+1) = X(1)*X(5) + X(2)*X(6) - X(3)*X(7) - X(4)*X(8) - . SIGMA(3) C (END EQUATION ONE). C (START EQUATION TWO). FJ(4,1) = X(7) FJ(4,2) = X(8) FJ(4,3) = X(5) FJ(4,4) = X(6) FJ(4,5) = X(3) FJ(4,6) = X(4) FJ(4,7) = X(1) FJ(4,8) = X(2) FJ(4,NVARS+1) = X(1)*X(7) + X(2)*X(8) + X(3)*X(5) + X(4)*X(6) - . SIGMA(4) C (END EQUATION TWO). C (START EQUATION THREE). FJ(5,1) = X(5)**2 - X(7)**2 FJ(5,2) = X(6)**2 - X(8)**2 FJ(5,3) = -TWO*X(5)*X(7) FJ(5,4) = -TWO*X(6)*X(8) FJ(5,5) = TWO* (X(1)*X(5)-X(3)*X(7)) FJ(5,6) = TWO* (X(2)*X(6)-X(4)*X(8)) FJ(5,7) = -TWO* (X(1)*X(7)+X(3)*X(5)) FJ(5,8) = -TWO* (X(2)*X(8)+X(4)*X(6)) FJ(5,NVARS+1) = X(1)*FJ(5,1) + X(2)*FJ(5,2) + X(3)*FJ(5,3) + . X(4)*FJ(5,4) - SIGMA(5) C (END EQUATION THREE). C (START EQUATION FOUR). C FJ(6,1) = TWO*X(5)*X(7) C FJ(6,2) = TWO*X(6)*X(8) C FJ(6,3) = X(5)**2 - X(7)**2 C FJ(6,4) = X(6)**2 - X(8)**2 FJ(6,1) = -FJ(5,3) FJ(6,2) = -FJ(5,4) FJ(6,3) = FJ(5,1) FJ(6,4) = FJ(5,2) C FJ(6,5) = TWO* (X(3)*X(5)+X(1)*X(7)) C FJ(6,6) = TWO* (X(4)*X(6)+X(2)*X(8)) C FJ(6,7) = -TWO* (X(3)*X(7)-X(1)*X(5)) C FJ(6,8) = -TWO* (X(4)*X(8)-X(2)*X(6)) FJ(6,5) = -FJ(5,7) FJ(6,6) = -FJ(5,8) FJ(6,7) = FJ(5,5) FJ(6,8) = FJ(5,6) FJ(6,NVARS+1) = X(1)*FJ(6,1) + X(2)*FJ(6,2) + X(3)*FJ(6,3) + . X(4)*FJ(6,4) - SIGMA(6) C (END EQUATION FOUR). C (START EQUATION FIVE). FJ(7,1) = X(5)* (X(5)**2-THREE*X(7)**2) FJ(7,2) = X(6)* (X(6)**2-THREE*X(8)**2) FJ(7,3) = X(7)* (X(7)**2-THREE*X(5)**2) FJ(7,4) = X(8)* (X(8)**2-THREE*X(6)**2) FJ(7,5) = THREE* (X(1)*FJ(5,1)+X(3)*FJ(5,3)) FJ(7,6) = THREE* (X(2)*FJ(5,2)+X(4)*FJ(5,4)) FJ(7,7) = -THREE* (X(1)*FJ(6,1)+X(3)*FJ(6,3)) FJ(7,8) = -THREE* (X(2)*FJ(6,2)+X(4)*FJ(6,4)) FJ(7,NVARS+1) = X(1)*FJ(7,1) + X(2)*FJ(7,2) + X(3)*FJ(7,3) + . X(4)*FJ(7,4) - SIGMA(7) C (END EQUATION FIVE). C (START EQUATION SIX). C FJ(8,1) = -X(7)* (X(7)**2-THREE*X(5)**2) C FJ(8,2) = -X(8)* (X(8)**2-THREE*X(6)**2) C FJ(8,3) = X(5)* (X(5)**2-THREE*X(7)**2) C FJ(8,4) = X(6)* (X(6)**2-THREE*X(8)**2) FJ(8,1) = -FJ(7,3) FJ(8,2) = -FJ(7,4) FJ(8,3) = FJ(7,1) FJ(8,4) = FJ(7,2) C FJ(8,5) = THREE* (X(1)*FJ(6,1)+X(3)*FJ(6,3)) C FJ(8,6) = THREE* (X(2)*FJ(6,2)+X(4)*FJ(6,4)) C FJ(8,7) = THREE* (X(1)*FJ(5,1)+X(3)*FJ(5,3)) C FJ(8,8) = THREE* (X(2)*FJ(5,2)+X(4)*FJ(5,4)) FJ(8,5) = THREE* (X(1)*FJ(6,1)+X(3)*FJ(6,3)) FJ(8,6) = THREE* (X(2)*FJ(6,2)+X(4)*FJ(6,4)) FJ(8,7) = FJ(7,5) FJ(8,8) = FJ(7,6) FJ(8,NVARS+1) = X(1)*FJ(8,1) + X(2)*FJ(8,2) + X(3)*FJ(8,3) + . X(4)*FJ(8,4) - SIGMA(8) C (END EQUATION SIX). IF (MODE.NE.1) THEN FJ(1,NVARS+1) = FJ(1,NVARS+1) - SIGMA(1) FJ(2,NVARS+1) = FJ(2,NVARS+1) - SIGMA(2) END IF NFEV = NFEV + 1 RETURN END SUBROUTINE DQED(DQEDEV,MEQUA,NVARS,MCON,IND,BL,BU,X,FJAC,LDFJAC, . FNORM,IGO,IOPT,ROPT,IWA,WA) C***BEGIN PROLOGUE DQED C***DATE WRITTEN 851210 (YYMMDD) C***REVISION DATE 870204 (YYMMDD) C***CATEGORY NO. K1b,K1b1a2,K1b2a C***KEYWORDS NONLINEAR LEAST SQUARES, SIMPLE BOUNDS, C LINEAR CONSTRAINTS C***AUTHOR HANSON, R. J., SNLA C KROGH, F. T., JPL-CIT C***PURPOSE SOLVE NONLINEAR LEAST SQUARES AND NONLINEAR C EQUATIONS. USER PROVIDES SIMPLE BOUNDS, LINEAR C CONSTRAINTS AND EVALUATION CODE FOR THE FUNCTIONS. C***ROUTINES CALLED DQEDMN, DQEDEV, CHRCNT, XERRWV C***LONG DESCRIPTION C SUBROUTINE DQED (DQEDEV, MEQUA, NVARS, MCON, IND, BL, BU, X, C * FJ, LDFJ, RNORM, IGO, IOPT, ROPT, C * IWORK, WORK) C C C Table of Sections C ----------------- C 1. Introduction C ------------ C 2. Calling Sequence Explained C ------- -------- --------- C 3. Remarks on the Usage Examples C ------- -- --- ----- -------- C 4. Error Messages for DQED() C -------------------------- C 5. References C ---------- C C 1. Introduction C ------------ C This software package is available in both single and double C precision. The double precision version (type REAL*8) is C described below. For the REAL version of the C documentation substitute 'REAL' for 'DOUBLE PRECISION' in the C type statements. Change the names of the subprograms: 'DQED()' C to 'SQED()', 'DQEDEV()' to 'SQEDEV()', and 'D1MACH' to 'R1MACH.' C C The Fortran subprogram, DQED(), solves the constrained nonlinear C least squares problem: C C Minimize the sum of squares of MEQUA (generally nonlinear) C equations, C C f (x) = 0, I=1,...,MEQUA Eq. (1) C I C where x is a (vector) set of NVARS unknowns. (The vector C function with these MEQUA components is called f(x) in the C discussion that follows.) The components of x may have upper C and lower bounds given by the user. (In fact all of the C possible cases, no bounds, bounds at one end only, or upper and C lower bounds can be specified.) Linear constraints on the C unknowns, more general than simple bounds, can also be given. C These constraints can be of the equality or inequality type: C C a x + ... + a x = y , L = 1,...,MCON, C L1 1 L,NVARS NVARS L C Eq. (2) C C with bounds specified on the y , again given by the user. The C L C constraints can actually be slightly nonlinear. In this case C the constraints can be described as: C C g (x) = y , L = 1,...,MCON, Eq. (2') C L L C where bounds are specified on each y . The functions g (x) must C L L C be defined for all x in the set described by the simple bounds. C Experienced users may wish to turn directly to Examples 1 and 2, C listed below, before reading the subprogram documentation. C There is no size relation required for the problem dimensions C MEQUA, NVARS, and MCON except that MEQUA and NVARS are both C positive, and MCON is nonnegative. C C This code package will do a decent job of solving most nonlinear C least squares problems that can be expressed as Eqs. (1) and (2) C above, provided that continuous derivatives of the functions C with respect to the parameters can be computed. This can also C include problems where the derivatives must be computed using C some form of numerical differentiation. Numerical C differentiation is not provided with this software for solving C nonlinear least squares problems. Refer to the subprogram C JACG for numerical differentiation. (Note: D. Salane has this C submitted to TOMS. It is not included here.) C C The authors also plan to develop methods that will do a much C better job of coping with constraints more general than the C essentially linear ones indicated above in Eqs. (2)-(2'). There C are nonlinear least squares problems with innocent looking but C highly nonlinear constraints where this package will fail to C work. The authors also hope to reduce the overhead required by C the software. This high overhead is due primarily to the method C used to solve the inner-loop quadratic model problem. The C authors recommend that users consider using the option number C 14, described below, to suppress use of the quadratic model. The C user may find that the software works quite well without the C quadratic model. This may be important when the function and C derivatives evaluations are not expensive but many individual C problems are being solved. C C There are two fundamental ways to use the subprogram DQED(). C The most staightforward way is to make one Fortran CALL to the C subprogram and obtain values for the unknowns, x. The user C provides a subprogram DQEDEV(), described below, that gives the C subprogram DQED() values of the functions f(x) and g(x), and the C derivative or Jacobian matrices for f(x) and g(x) at each C desired point x. This usage is called 'forward communication.' C An alternate way to use the subprogram is to provide an option C that allows the user to communicate these values by 'reverse C communication.' The subprogram returns to the calling program C unit and requests values for f(x) and g(x), and the Jacobian C matrices for f(x) and g(x) for a given value of x. (This C framework is often required in applications that have C complicated algorithmic requirements for evaluation of the C functions.) An example using both 'forward' and 'reverse' C communication is provided below (see Remarks on the Usage C Examples) for least squares fitting of two exponential functions C to five data points. C C 2. Calling Sequence Explained C ------- -------- --------- C There are arrays used by the subprogram that must have C dimensions equivalent to the following declarations. C C INTEGER MEQUA, NVARS, MCON, LDFJ, IGO C INTEGER IND(NVARS+MCON), IOPT(LIOPT), IWORK(LIWORK) C C DOUBLE PRECISION BL(NVARS+MCON), BU(NVARS+MCON), X(NVARS), RNORM, C *ROPT(LROPT), FJ(LDFJ,NVARS+1), WORK(LWORK) C C EXTERNAL DQEDEV C The array dimensions must satisfy the bounds: C LIOPT .ge. Number required for options in use. C LROPT .ge. Number required for options in use. C LDFJ .ge. MEQUA+MCON, C The array dimensions for the arrays IWORK(*) and WORK(*) can C change if either option 14 or option 15 are in use. For use in C the formulas, define: C C MC=MCON C ME=MEQUA C NV=NVARS C MX=MAX(MEQUA,NVARS) C If the user is not using option 15, then C NT=5. C If the user is using option 15, then C NT=new number, must be .ge. 2. C If the user is not using option 14, then C NA=MC+2*NV+NT. C If the user is using option 14, then C NA=MC+NV+1. C C C In terms of these values defined above, C LIWORK .ge. 3*MC+9*NV+4*NT+NA+10 C LWORK .ge. NA*(NA+4)+NV*(NT+33)+(ME+MX+14)*NT+9*MC+26 C C The subprogram DQEDEV must be declared in a Fortran EXTERNAL C statement: C C EXTERNAL DQEDEV C C Initialize the values of the parameters: C MEQUA, NVARS, MCON, IND(*), BL(*), BU(*), X(*), LDFJ, C IOPT(*), IWORK(1), IWORK(2), C C CALL DQED (DQEDEV, MEQUA, NVARS, MCON, IND, BL, BU, X, C * FJ, LDFJ, RNORM, IGO, IOPT, ROPT, C * IWORK, WORK) C C Subprogram parameters: C C DQEDEV (Input) C ----- C This is the name of a subprogram that the user will usually C supply for evaluation of the values of the constraints and C model, and the derivatives of these functions. The user must C provide this subprogram unless 'reverse communication' is used. C A model for writing the subprogram DQEDEV() is provided under C the heading Example 1 Using Forward Communication, listed below. C Users may find it convenient to modify this model subprogram C when writing a subprogram for their own application. If C 'reverse communication' is used, the user does not need to write C a stub or dummy subroutine named DQEDEV(). All that is required C is to declare exactly this name in an EXTERNAL statement. The C code package has a dummy subroutine DQEDEV() that will be used C in the linking or load step. Example 2 Using Reverse C Communication, listed below, illustrates this detail. C C MEQUA, NVARS, MCON (Input) C ------------------ C Respectively they are: The number of least squares equations, C the number of unknowns or variables, and the number of general C constraints for the solution, not including simple bounds. The C values of MEQUA and NVARS must be positive; the value of MCON C must be nonnegative. Other values for these parameters are C errors. C C IND(*),BL(*),BU(*) (Input) C ------------------ C These arrays describe the form of the simple bounds that the C components of x are to satisfy. Components numbered 1,...,NVARS C are used to describe the form of the simple bounds that the C unknown are to satisfy. Components numbered C NVARS+1,...,NVARS+MCON are used to describe the form of the C general MCON linear constraints. The first NVARS components of C IND(*) indicate the type of simple bounds that the solution is C to satisfy. The corresponding entries of BL(*) and BU(*) are C the bounding value. The only values of IND(*) allowed are 1,2,3 C or 4. Other values are errors. Specifically: C C IND(J)=1, if x .ge. BL(J) is required; BU(J) is not used. C J C =2, if x .le. BU(J) is required; BL(J) is not used. C J C =3, if x .ge. BL(J) and x .le. BU(J) is required. C J J C =4, if no bounds on x are required; C J C BL(*),BU(*) are not used. C General linear constraints of the form shown in Eq. (2) require C that bounds be given for the linear functions y . Specifically: C L C C IND(NVARS+L)=1, if y .ge. BL(NVARS+L) is required; BU(*) is not C L C needed. C C =2, if y .le. BU(NVARS+L) is required; BL(*) is not C L C needed. C =3, if y .ge. BL(NVARS+L) and y .le. BU(NVARS+L) C L L C C =4, if no bounds on y are required; C L C BL(*),BU(*) are not used. C C The values of the bounds for the unknowns may be changed by the C user during the evaluation of the functions f(x) and g(x) and C their Jacobian matrices. C C X(*),FJ(*,*),LDFJ (Input and Output, except LDFJ which is Input) C ----------------- C The array X(*) contains the NVARS values, x, where the C functions f(x) and g(x) and their Jacobian matrices will be C evaluated by the subprogram DQED(). After the computation has C successfully completed, the array X(*) will contain a solution, C namely the unknowns of the problem, x. (Success is determined C by an appropriate value for IGO. This parameter is described C below.) Initially the array X(*) must contain a starting guess C for the unknowns of the problem, x. The initial values do not C need to satisfy the constraints or the bounds. If they do not C satisfy the bounds, then the point will be simply projected onto C the bounds as a first step. The first linear change to the C values of x must satisfy the general constraints. (It is here C that the assumption of their linearity is utilized.) C C The Fortran two-dimensional array FJ(*,*) is used to store the C linear constraints of Eq. (2) (or more generally the Jacobian C matrix of the functions g(x) with respect to the variables x), C and the Jacobian matrix of the function f(x). The Jacobian C matrix of the (linear) constraints is placed in rows 1,...,MCON. C The Jacobian matrix of f(x) is placed in rows MCON+1, ..., C MCON+MEQUA. The parameter LDFJ is the leading or row dimension C of the array FJ(*,*). Normally the array FJ(*,*) is assigned C values by the user when the nonlinear solver requests C evaluations of the constraints g(x) and the function f(x) C together with the Jacobian matrices G(x) and J(x). The values C of the constraint functions g (x) are placed in the array C L C FJ(L,NVARS+1), L=1,...,MCON. The values of the model functions C f (x) are placed in the array at entries FJ(MCON+I,NVARS+1), C I C I=1,...,MEQUA. Note that the second dimension of FJ(*,*) must C be at least NVARS+1 in value. C C RNORM (Output) C ----- C This is the value of the Euclidean length or square root of sums C of squares of components of the function f(x) after the C approximate solution, x, has been found. During the computation C it is updated and equals the best value of the length of f(x) C that has been found. C C IGO (Output; it can be an Input if interrupting the code) C --- C This flag directs user action and informs the user about the C type of results obtained by the subprogram. The user may find C it convenient to treat the cases abs(IGO) .le. 1 the same as C IGO=1. This has no effect on the solution process. C C The user can interrupt the computation at any time, obtaining C the best values of the vector x up to this point, by setting C IGO to any value .gt. 1 and then return control to DQED(). For C example if a calculation must be done in a certain length of C time, the user can, as the end of the time draws near, set C IGO=20 and return control to DQED(). It is important that this C method be used to stop further computing, rather than simply C proceeding. The reason for this is that certain flags in DQED() C must be reset before any further solving on subsequent problems C can take place. The value of IGO .gt. 1 used to interrupt the C computation is arbitrary and is the value of IGO returned. If C values of IGO =2,...,18 are used to flag this interrupt, they do C not mean the same thing as indicated below. For this reason the C value IGO=20 is recommended for signaling interrupts in DQED(). C Another situation that may occur is the request for an C evaluation of the functions and derivatives at a point x where C these can't be evaluated. If this occurs, set IGO=99 and return C control to DQED(). This will have the effect of defining the C derivatives to be all zero and the functions to be 'large.' C Thus a reduction in the trust region around the current best C estimate will occur. Assigning the value IGO=99 will not cause C DQED() to stop computing. C C =0 Place the value of f(x) in FJ(MCON+*,NVARS+1). If C 'reverse communication' is being used, CALL DQED() again. If C 'forward communication' is being used, do a RETURN. C C =1 or (-1) Evaluate the Jacobians for the functions g(x) C and f(x) as well as evaluating g(x) and f(x). Use the vector x C that is now in the array X(*) as the values where this C evaluation will be performed. Place the Jacobian matrix C for g(x) in the first MCON rows of FJ(*,*). Place the Jacobian C matrix for f(x) in rows MCON+1,...,MCON+MEQUA in FJ(*,*). Place C the value of g(x) in FJ(*,NVARS+1). Place the value of f(x) in C FJ(MCON+*,NVARS+1). C C (Users who have complicated functions whose derivatives cannot be C computed analytically may want to use the numerical differentiation C subroutine JAGC. This is available on the SLATEC library.) C C If 'reverse communication' is being used, CALL DQED() again. C If 'forward communication' is being used, do a RETURN. C C A value IGO=(-1) flags that that the number of terms in the C quadratic model is being restricted by the amount of storage C given for that purpose. It is suggested, but it is not C required, that additional storage be given for the quadratic C model parameters. See the following description of The Option C Array, option number 15, for the way to designate more storage C for this purpose. C C =2 The function f(x) has a length less than TOLF. This is C the value for IGO to be expected when an actual zero value of C f(x) is anticipated. See the description of The Option Array C for the value. C C =3 The function f(x) has reached a value that may be a C local minimum. However, the bounds on the trust region C defining the size of the step are being hit at each step. Thus C the situation is suspect. (Situations of this type can occur C when the solution is at infinity in some of the components of C the unknowns, x. See the description of The Option Array for C ways to avoid this value of output value of IGO. C C =4 The function f(x) has reached a local minimum. This is C the value of IGO that is expected when a nonzero value of f(x) C is anticipated. See the description of The Option Array for the C conditions that have been satisfied. C C =5 The model problem solver has noted a value for the C linear or quadratic model problem residual vector length that is C .ge. the current value of the function, i.e. the Euclidean C length of f(x). This situation probably means that the C evaluation of f(x) has more uncertainty or noise than is C possible to account for in the tolerances used to note a local C minimum. The value for x is suspect, but a minimum has probably C been found. C C =6 A small change (absolute) was noted for the vector x. A C full model problem step was taken. The condition for IGO=4 may C also be satisfied, so that a minimum has been found. However, C this test is made before the test for IGO=4. C C =7 A small change (relative to the length of x) was noted C for the vector x. A full model problem step was taken. The C condition for IGO=4 may also be satisfied, so that a minimum has C been found. However, this test is made before the test for C IGO=4. C C =8 More than ITMAX iterations were taken to obtain the C solution. The value obtained for x is suspect, although it is C the best set of x values that occurred in the entire C computation. See the description of The Option Array for C directions on how to increase this value. (Note that the C nominal value for ITMAX, 75, is sufficient to solve all of the C nonlinear test problems described in Ref. (2).) C C =9-18 Errors in the usage of the subprogram were noted. C The exact condition will be noted using an error processor that C prints an informative message unless this printing has been C suppressed. A minimum value has not been found for x. The C relation between IGO and the error number are IGO=NERR + 8. C Here NERR is the identifying number. See below, Error Messages C for DQED(). C The Option Array C --- ------ ----- C Glossary of Items Modified by Options. Those items with Nominal C Values listed can be changed. C C Names Nominal Definition C Values C ----- ------- ---------- C FC Current value of length of f(x). C FB Best value of length of f(x). C FL Value of length of f(x) at the C previous step. C PV Predicted value of length of f(x), C after the step is taken, using the C approximating model. C The quantity 'eps', used below, is the machine precision C parameter. Its value is obtained by a call to the Bell Labs. C Port subprogram D1MACH(4). It is machine dependent. C MIN(1.D-5, C TOLF sqrt(eps)) Tolerance for stopping when C FC .le. TOLF. C MIN(1.D-5, C TOLD sqrt(eps)) Tolerance for stopping when C change to x values has length C .le. TOLD. C MIN(1.D-5, C TOLX sqrt(eps)) Tolerance for stopping when C change to x values has length C .le. TOLX*length of x values. C TOLSNR 1.D-5 Tolerance used in stopping C condition IGO=4. Explained below. C TOLP 1.D-5 Tolerance used in stopping C condition IGO=4. Explained below. C C (The conditions (abs(FB-PV).le.TOLSNR*FB and abs(FC-PV) .le. C TOLP*FB) and (ABS(FC-FL).le.TOLSNR*FB) together with taking C a full model step, must be satisfied before the condition IGO=4 C is returned. Decreasing any of the values for TOLF, TOLD, TOLX, C TOLSNR, or TOLP will likely increase the number of iterations C required for convergence.) C C COND 30. Largest condition number to allow C when solving for the quadratic C model coefficients. Increasing C this value may result in more C terms being used in the quadratic C model. C TOLUSE sqrt(eps) A tolerance that is used to avoid C values of x in the quadratic C model's interpolation of previous C points. Decreasing this value may C result in more terms being used in C the quadratic model. C ITMAX 75 The number of iterations to take C with the algorithm before giving C up and noting it with the value C IGO=8. C IPRINT 0 Control the level of printed C output in the solver. A value C of IPRINT .gt. 0 will result in C output of information about each C iteration. The output unit used C is obtained using the Bell Labs. C Port subprogram, i. e. I1MACH(2). C LEVEL 1 Error processor error level. See C the SLATEC library documentation C for XERROR() for an explanation. C NTERMS 5 One more than the maximum number C of terms used in the quadratic C model. C C IOPT(*) (Input) C ------- C In order to use the option array technique to change selected C data within a subprogram, it is necessary to understand how this C array is processed within the software. Let LP designate the C processing pointer that moves to positions of the IOPT(*) array. C Initially LP=1, and as each option is noted and changed, the C value of LP is updated. The values of IOPT(LP) determine what C options get changed. The amount that LP changes is known by the C software to be equal to the value two except for two options. C These exceptional cases are the last option (=99) and the 'leap' C option (=13) which advances LP by the value in IOPT(LP+1). A C negative value for IOPT(LP) means that this option is not to be C changed. This aids the programmer in using options; often the C code for using an option can be in the calling program but a C negative value of the option number avoids rewriting code. C C Option Usage Example C ------ ----- ------- C In the Fortran code fragment that follows, an example is given C where we change the value of TOLF and decrease the maximum C number of iterations allowed from 75 to 30. C In this example the dimensions of IOPT(*) and ROPT(*) must C satisfy: C C DOUBLE PRECISION ROPT(01) C INTEGER IOPT(005) C . C . C . C C SET THE OPTION TO CHANGE THE VALUE OF TOLF. C C IOPT(01)=4 C C C THE NEXT ENTRY POINTS TO THE PLACE IN ROPT(*) WHERE C C THE NEW VALUE OF TOLF IS LOCATED. C C IOPT(02)=1 C C THIS IS THE NEW VALUE OF TOLF. THE SPECIFIC VALUE C C 1.D-9 IS USED HERE ONLY FOR ILLUSTRATION. C C ROPT(01)=1.D-9 C C C CHANGE THE NUMBER OF ITERATIONS. C C IOPT(03)=2 C C C THIS NEXT ENTRY IS THE NEW VALUE FOR THE MAXIMUM NUMBER OF C C ITERATIONS. C C IOPT(04)=30 C C C THIS NEXT OPTION IS A SIGNAL THAT THERE ARE NO MORE C C OPTIONS. C C IOPT(05)=99 C . C . C . C CALL DQED() C . C . C . C Option Values Explanation C ------ ------ ----------- C =99 There are no more options to change. C Normally this is the first and only C option that a user needs to specify, C and it can be simply IOPT(01)=99. The C total dimension of IOPT(*) must be at C least 17, however. This can lead to a C hard-to-find program bug if the dimension C is too small. C C = 1 Change the amount of printed output. C The next value of IOPT(*) is the print C level desired, IPRINT. Any value of C IPRINT .gt. 0 gives all the available C output. C C = 2 Change the value of ITMAX. The next value C of IOPT(*) is the value of ITMAX desired. C C = 3 Pass prior determined bounds for the box C containing the initial point. This box is the C trust region for the first move from the initial C point. The next entry in IOPT(*) points to C the place in ROPT(*) where the NVARS values for C the edges of the box are found. C C = 4 Change the value of TOLF. The next entry of C IOPT(*) points to the place in ROPT(*) where the C new value of TOLF is found. C C = 5 Change the value of TOLX. The next entry of C IOPT(*) points to the place in ROPT(*) where the C new value of TOLX is found. C C = 6 Change the value of TOLD. The next entry of C IOPT(*) points to the place in ROPT(*) where the C new value of TOLD is found. C C = 7 Change the value of TOLSRN. The next entry of C IOPT(*) points to the place in ROPT(*) where the C new value of TOLSNR is found. C C = 8 Change the value of TOLP. The next entry of C IOPT(*) points to the place in ROPT(*) where the C new value of TOLP is found. C C = 9 Change the value of TOLUSE. The next entry of C IOPT(*) points to the place in ROPT(*) where the C new value of TOLUSE is found. C C =10 Change the value of COND. The next entry of C IOPT(*) points to the place in ROPT(*) where the C new value of COND is found. C C =11 Change the value of LEVEL. The next entry of C IOPT(*) is the new value of LEVEL. C C =12 Pass an option array to the subprogram DQEDGN() C used as the inner loop solver for the C model problem. The next entry of IOPT(*) is the C starting location for the option array for C DQEDGN() within the array IOPT(*). Thus the C option array for DQEDGN() must be a part of C the array IOPT(*). C C =13 Move (or leap) the processing pointer LP for the C option array by the next value in IOPT(*). C C =14 Change a logical flag that suppresses the C use of the quadratic model in the inner C loop. Use the next value in IOPT(*) for C this flag. If this value = 1, then never C use the quadratic model. (Just use the C linear model). Otherwise, use the quadratic C model when appropriate. This option decreases C the amount of scratch storage as well as the C computing overhead required by the code package. C A user may want to determine if the application C really requires the use of the quadratic model. C If it does not, then use this option to save C both storage and computing time. C C =15 Change, NTERMS, the maximum number of array C columns that can be used for saving quadratic C model data. (The value of NTERMS is one more C than the maximum number of terms used.) Each C unit increase for NTERMS increases the required C dimension of the array WORK(*) by 2*MEQUA+NVARS. C Use the value in IOPT(LP+1) for the new value C of NTERMS. Decreasing this value to 2 (its C minimum) decreases the amount of storage C required by the code package. C C =16 Change a logical flag so that 'reverse C communication' is used instead of 'forward C communication.' Example EX01, listed below, C uses 'forward communication.' Example EX02, C also listed below, uses 'reverse communication.' C Use the next value in IOPT(*) for C this flag. If this value = 1, then C use 'reverse communication.' Otherwise, C use 'forward communication.' WARNING: This C usage may not work unless the operating system C saves variables between subroutine calls to DQED. C C =17 Do not allow the flag IGO to return with the C value IGO=3. This means that convergence will C not be claimed unless a full model step is taken. C Normal output values will then be IGO = 2,4,6 or 7. C Use the next value in IOPT(*) for this flag. If C this value = 1, then force a full model step. C Otherwise, do not force a full model step if small C steps are noted. C C IWORK(*), WORK(*) (Input and Output) C ---------------- C These are scratch arrays that the software uses for storage of C intermediate results. It is important not to modify the C contents of this storage during the computation. C C The array locations IWORK(1) and IWORK(2) must contain the C actual lengths of the arrays WORK(*) and IWORK(*) before the C call to the subprogram. These array entries are replaced by the C actual amount of storage required for each array. If the amount C of storage for either array is too small, an informative error C message will be printed, and the value IGO=13 or 14 returned. C C The user may find it useful to let the subprogram DQED() return C the amounts of storage required for these arrays. For example C set IWORK(1)=1, IWORK(2)=1. The subprogram will return with C IGO=13, IWORK(1)=required length of WORK(*), and C IWORK(2)=required length of IWORK(*). (Appropriate error C messages will normally be printed.) C C 3. Remarks on the Usage Examples C ------- -- --- ----- -------- C The following complete program units, EX01 and EX02, show how C one can use the nonlinear solver for fitting exponential C functions to given data. These examples are calculations that C match two terms of an exponential series to five given data C points. There are some subtle points about exponential fitting C that are important to note. First, the signs of the C exponential arguments are restricted to be nonpositive. C The size of the arguments should not be much larger than the start C of the time data (reciprocated). This is the reason the lower C bounds are set a bit less than the reciprocal of the time value. C In many applications that require exponential modeling this is a C natural assumption. The nonlinear solver allows these bounds C on the arguments explicitly. In addition, the coefficients are C constrained to be nonnegative. These bounds are harder to C justify. The idea is to avoid the situation where a coefficient C is very large and negative, and the corresponding exponential C argument is also large and negative. The resulting contribution C to the series may be very small, but its presence is spurious. C Finally, the single general linear constraint that keeps the C arguments separated (by 0.05 in this example) is used for two C purposes. First, it naturally orders these values so that the C first one is algebraically largest. Second, this constraint C moves the parameters from the local minimum corresponding to the C initial values used in the examples. This constraint also C retains the validity of the model function h(t) = w*exp(x*t) + C y*exp(z*t). Namely, if the arguments are allowed to coalesce to C the same value, then the model itself must change. The form of C the model can become h(t)=(a+b*t)*exp(c*t) or h(t) = d*exp(e*t). C Either one could occur, and the choice is problem dependent. C Example 1 Using Forward Communication C --------- ----- ------- ------------- C PROGRAM EX01 C CC Illustrate the use of the Hanson-Krogh nonlinear least CC squares solver for fitting two exponentials to data. CC CC The problem is to find the four variables x(1),...,x(4) CC that are in the model function CC CC h(t) = x(1)*exp(x(2)*t) + x(3)*exp(x(4)*t) CC There are values of h(t) given at five values of t, CC t=0.05, 0.1, 0.4, 0.5, and 1.0. CC We also have problem constraints that x(2), x(4) .le. 0, x(1), CC x(3) .ge. 0, and a minimal separation of 0.05 between x(2) and CC x(4). Nothing more about the values of the parameters is known CC except that x(2),x(4) are approximately .ge. 1/min t. CC Thus we have no further knowledge of their values. CC For that reason all of the initial values are set to zero. CC CC Dimension for the nonlinear solver. C DOUBLE PRECISION FJ(6,5),BL(5),BU(5),X(4),ROPT(001),WA(640) CC EDIT on 950228-1300: C DOUBLE PRECISION RNORM C INTEGER IND(5),IOPT(24),IWA(084) C C EXTERNAL DQEDEX C C DATA LDFJ,LWA,LIWA/6,640,084/ C C MCON = 1 C MEQUA = 5 C NVARS = 4 CC Define the constraints for variables. C BL(1) = 0. C BL(2) = -25. C BU(2) = 0. C BL(3) = 0. C BL(4) = -25. C BU(4) = 0. CC Define the constraining value (separation) for the arguments. C BL(5) = 0.05 CC Define all of the constraint indicators. C IND(1) = 1 C IND(2) = 3 C IND(3) = 1 C IND(4) = 3 C IND(5) = 1 CC Define the initial values of the variables. CC We don't know anything more, so all variables are set zero. C DO 10 J = 1,NVARS C X(J) = 0.D0 C 10 CONTINUE CC Tell how much storage we gave the solver. C IWA(1) = LWA C IWA(2) = LIWA CC No additional options are in use. C IOPT(01) = 99 C CALL DQED(DQEDEX,MEQUA,NVARS,MCON,IND,BL,BU,X,FJ,LDFJ,RNORM,IGO, C . IOPT,ROPT,IWA,WA) C NOUT = 6 C WRITE (NOUT,9001) (X(J),J=1,NVARS) C WRITE (NOUT,9011) RNORM C WRITE (NOUT,9021) IGO C C STOP C C 9001 FORMAT (' 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) = ',/,4F12.6) C 9011 FORMAT (' RESIDUAL AFTER THE FIT = ',1PD12.4) C 9021 FORMAT (' OUTPUT FLAG FROM SOLVER =',17X,I6) C END C SUBROUTINE DQEDEX(X,FJ,LDFJ,IGO,IOPT,ROPT) CC This is the subprogram for evaluating the functions CC and derivatives for the nonlinear solver, DQED. CC CC The user problem has MCON constraint functions, CC MEQUA least squares equations, and involves NVARS CC unknown variables. CC CC When this subprogram is entered, the general (near) CC linear constraint partial derivatives, the derivatives CC for the least squares equations, and the associated CC function values are placed into the array FJ(*,*). CC All partials and functions are evaluated at the point CC in X(*). Then the subprogram returns to the calling CC program unit. Typically one could do the following CC steps: CC CC step 1. Place the partials of the i-th constraint CC function with respect to variable j in the CC array FJ(i,j), i=1,...,MCON, j=1,...,NVARS. CC step 2. Place the values of the i-th constraint CC equation into FJ(i,NVARS+1). CC step 3. Place the partials of the i-th least squares CC equation with respect to variable j in the CC array FJ(MCON+i,j), i=1,...,MEQUA, CC j=1,...,NVARS. CC step 4. Place the value of the i-th least squares CC equation into FJ(MCON+i,NVARS+1). CC step 5. Return to the calling program unit. C DOUBLE PRECISION FJ(LDFJ,*),X(*),ROPT(*) C DOUBLE PRECISION T(5),F(5) C INTEGER IOPT(*) C C DATA T/0.05,0.10,0.40,0.50,1.00/ C DATA F/2.206D+00,1.994D+00,1.350D+00,1.216D+00,.7358D0/ C C DATA MCON,MEQUA,NVARS/1,5,4/ C CC Define the derivatives of the constraint with respect to the x(j). C FJ(1,1) = 0.D0 C FJ(1,2) = 1.D0 C FJ(1,3) = 0.D0 C FJ(1,4) = -1.D0 CC Define the value of this constraint. C FJ(1,5) = X(2) - X(4) CC Define the derivatives and residuals for the data model. C DO 10 I = 1,MEQUA C E1 = EXP(X(2)*T(I)) C E2 = EXP(X(4)*T(I)) C FJ(MCON+I,1) = E1 C FJ(MCON+I,2) = X(1)*T(I)*E1 C FJ(MCON+I,3) = E2 C FJ(MCON+I,4) = X(3)*T(I)*E2 C FJ(MCON+I,5) = X(1)*E1 + X(3)*E2 - F(I) C 10 CONTINUE C RETURN C END C Output from Example 1 Program C ------ ---- --------- ------- 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 C C Example 2 Using Reverse Communication C --------- ----- ------- ------------- C PROGRAM EX02 C CC Illustrate the use of the Hanson-Krogh nonlinear least CC squares solver for fitting two exponentials to data. CC CC The problem is to find the four variables x(1),...,x(4) CC that are in the model function CC CC h(t) = x(1)*exp(x(2)*t) + x(3)*exp(x(4)*t) CC There are values of h(t) given at five values of t, CC t=0.05, 0.1, 0.4, 0.5, and 1.0. CC We also have problem constraints that x(2), x(4) .le. 0, x(1), CC x(3) .ge. 0, and a minimal separation of 0.05 between x(2) and CC x(4). Nothing more about the values of the parameters is known CC except that x(2),x(4) are approximately .ge. 1/min t. CC Thus we have no further knowledge of their values. CC For that reason all of the initial values are set to zero. CC CC Dimension for the nonlinear solver. C DOUBLE PRECISION FJ(6,5),BL(5),BU(5),X(4),ROPT(001),WA(640) CC EDIT on 950228-1300: C DOUBLE PRECISION RNORM C INTEGER IND(5),IOPT(24),IWA(084) C DOUBLE PRECISION T(5),F(5) C C EXTERNAL DQEDEV C C DATA LDFJ,LWA,LIWA/6,640,084/ C C DATA T/0.05,0.10,0.40,0.50,1.00/ C DATA F/2.206D+00,1.994D+00,1.350D+00,1.216D+00,.7358D0/ C C MCON = 1 C MEQUA = 5 C NVARS = 4 CC Define the constraints for variables. C BL(1) = 0. C BL(2) = -25. C BU(2) = 0. C BL(3) = 0. C BL(4) = -25. C BU(4) = 0. CC Define the constraining value (separation) for the arguments. C BL(5) = 0.05 CC Define all of the constraint indicators. C IND(1) = 1 C IND(2) = 3 C IND(3) = 1 C IND(4) = 3 C IND(5) = 1 CC Define the initial values of the variables. CC We don't know anything at all, so all variables are set zero. C DO 10 J = 1,NVARS C X(J) = 0.D0 C 10 CONTINUE CC Tell how much storage we gave the solver. C IWA(1) = LWA C IWA(2) = LIWA C NITERS = 0 CC TELL HOW MUCH STORAGE WE GAVE THE SOLVER. C IWA(1) = LWA C IWA(2) = LIWA CC USE REVERSE COMMUMICATION TO EVALUATE THE DERIVATIVES. C IOPT(01)=16 C IOPT(02)=1 CC NO MORE OPTIONS. C IOPT(03) = 99 C 20 CONTINUE C CALL DQED(DQEDEV,MEQUA,NVARS,MCON,IND,BL,BU,X,FJ,LDFJ,RNORM, C .IGO,IOPT, ROPT,IWA,WA) C IF (IGO.GT.1) GO TO 40 CC COUNT FUNCTION EVALUATIONS. C NITERS = NITERS + 1 CC DEFINE THE DERIVATIVES OF THE CONSTRAINT WITH RESPECT TO THE X(J). C FJ(1,1) = 0.D0 C FJ(1,2) = 1.D0 C FJ(1,3) = 0.D0 C FJ(1,4) = -1.D0 CC DEFINE THE VALUE OF THIS CONSTRAINT. C FJ(1,5) = X(2) - X(4) CC DEFINE THE DERIVATIVES AND RESIDUALS FOR THE DATA MODEL. C DO 30 I = 1,MEQUA C E1 = EXP(X(2)*T(I)) C E2 = EXP(X(4)*T(I)) C FJ(MCON+I,1) = E1 C FJ(MCON+I,2) = X(1)*T(I)*E1 C FJ(MCON+I,3) = E2 C FJ(MCON+I,4) = X(3)*T(I)*E2 C FJ(MCON+I,5) = X(1)*E1 + X(3)*E2 - F(I) C 30 CONTINUE C GO TO 20 C C 40 CONTINUE C NOUT = 6 C WRITE (NOUT,9001) (X(J),J=1,NVARS) C WRITE (NOUT,9011) RNORM C WRITE (NOUT,9021) NITERS, IGO C C 9001 FORMAT (' 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) = ',/,4F12.6) C 9011 FORMAT (' RESIDUAL AFTER THE FIT = ',1PD12.4) C 9021 FORMAT (' NUMBER OF EVALUATIONS OF PARAMETER MODEL =',I6,/, C . ' OUTPUT FLAG FROM SOLVER =',17X,I6) C STOP C END C Output from Example 2 Program C ------ ---- --------- ------- 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 C 4. Error Messages for DQED() C -------------------------- C 'DQED. VALUE OF MEQUA=NO. OF EQUAS. MUST .GT.0. NOW = (I1).' C NERR = 01 C IGO=9 C C 'DQED. VALUE OF NVARS=NO. OF EQUAS. MUST .GT.0. NOW = (I1).' C NERR = 02 C IGO=10 C C 'DQED. VALUE OF MCON=NO. OF EQUAS. MUST .GE.0. NOW = (I1).' C NERR = 03 C IGO=11 C C 'DQED. INVALID OPTION PROCESSED. I1=IOPT(*) ENTRY. I2=IOPT(I1).' C NERR = 04 C IGO=12 C C 'DQED. WA(*) STORAGE SHORT. I1=AMOUNT NEEDED. I2=AMOUNT GIVEN.' C NERR = 05 C IGO=13 C C 'DQED. IWA(*) STORAGE SHORT. I1=AMOUNT NEEDED. I2=AMOUNT GIVEN.' C NERR = 06 C IGO=14 C C 'DQEDMN. INVALID OPTION PROCESSED. I1=IOPT(*) ENTRY. I2=IOPT(I1). C C NERR=07 C IGO=15 C C 'DQEDIP. INVALID OPTION PROCESSED. I1=IOPT(*) ENTRY. I2=IOPT(I1).' C C NERR=08 C IGO=16 C C 'DQED. THE EVALUATOR PROGRAM DQEDEV MUST BE WRITTEN BY THE USER.' C NERR=09 C IGO=17 C C 'DQED. BOUND INDICATORS MUST BE 1-4. NOW I1=J, I2=IND(I1).' C NERR=10 C IGO=18 C C 5. References C ---------- C***REFERENCES C Dongarra, J. J., Bunch, J. R., Moler, C. B, Stewart, G. W., C LINPACK User's Guide, Soc. Indust. and Appl. Math, Phil., C PA, (1979). C Hanson, R. J., "Least Squares with Bounds and Linear C Constraints," SIAM J. Sci. Stat. Comput., vol. 7, no. 3, July, C (1986), p. 826-834. C Schnabel, R. B., Frank, P. D, "Tensor Methods for Nonlinear C Equations," SIAM J. Num. Anal., vol. 21, no. 5, Oct., (1984), C p. 815-843. C***END PROLOGUE DQED C REVISED 870204-1100 C REVISED YYMMDD-HHMM DOUBLE PRECISION BL(*),BU(*),X(*),FJAC(LDFJAC,*) DOUBLE PRECISION ROPT(*),WA(*) DOUBLE PRECISION FNORM REAL RDUM INTEGER IND(*),IOPT(*),IWA(*) LOGICAL NOQUAD CHARACTER XMESS*128 EXTERNAL DQEDEV C--PROCEDURES-- C -NAME------TYPE--------ARGS------CLASS----- C C DQEDEV 6 EXTERNAL C DQEDMN 35 SUBROUTINE C CHRCNT 2 SUBROUTINE C XERRWV 10 SUBROUTINE C DQED: C Name Memory Status Type Argument Uses and comments. C Status C ---- ------------- ---- -------- ------------------ C BL DUMMY-ARG REAL ADJ-ARY Lower bounds C BU DUMMY-ARG REAL ADJ-ARY Upper bounds C FJAC DUMMY-ARG REAL ADJ-ARY Jacobian array C FNORM DUMMY-ARG REAL Norm at solution C I /S$A$V$E/ SAV INTEGER Dummy loop variable C IDUM /S$A$V$E/ SAV INTEGER Dummy for error pack C IFLAG /S$A$V$E/ SAV INTEGER Gate for reverse comm C IGO DUMMY-ARG INTEGER Directs user action C IGOOK /S$A$V$E/ SAV INTEGER Internal gate for errors C IIWAAV /S$A$V$E/ SAV INTEGER Length claimed for IWA C IND DUMMY-ARG INTEGER ADJ-ARY Bound indicators C IOPT DUMMY-ARG INTEGER ADJ-ARY Option array C IWA DUMMY-ARG INTEGER ADJ-ARY Work array C IWAAV /S$A$V$E/ SAV INTEGER Length claime for WA C J /S$A$V$E/ SAV INTEGER Dummy loop variable C MILAST /S$A$V$E/ SAV INTEGER Last integer in IWA C MIND /S$A$V$E/ SAV INTEGER Point to start IND C MINDB /S$A$V$E/ SAV INTEGER Point to start INDB C MPJ /S$A$V$E/ SAV INTEGER Point to start PJ C MQC /S$A$V$E/ SAV INTEGER Point to start QC C MUT /S$A$V$E/ SAV INTEGER Point to start UT C MWA /S$A$V$E/ SAV INTEGER Point to start WA C MWJ /S$A$V$E/ SAV INTEGER Point to start WJ C MWLAST /S$A$V$E/ SAV INTEGER Last value in WA C MXB /S$A$V$E/ SAV INTEGER Point to start XB C MXP /S$A$V$E/ SAV INTEGER Point to start XP C MZP /S$A$V$E/ SAV INTEGER Point to start ZP C NALL /S$A$V$E/ SAV INTEGER Sum of dimensions C KP /S$A$V$E/ SAV INTEGER Dummy option loop pointer C LDFJAC DUMMY-ARG INTEGER Row dimension of FJAC C LEVEL /S$A$V$E/ SAV INTEGER Error processor status C LP /S$A$V$E/ SAV INTEGER Dummy option loop pointer C LPDIFF /S$A$V$E/ SAV INTEGER Dummy option loop diff C MB /S$A$V$E/ SAV INTEGER Point to start B C MBB /S$A$V$E/ SAV INTEGER Point to start BB C MBLB /S$A$V$E/ SAV INTEGER Point to start BLB C MBUB /S$A$V$E/ SAV INTEGER Point to start BUB C MCON DUMMY-ARG INTEGER Number of constraints C MDX /S$A$V$E/ SAV INTEGER Point to start DX C MDXL /S$A$V$E/ SAV INTEGER Point to start DXL C MEQUA DUMMY-ARG INTEGER Numer of least squares eqs C MGR /S$A$V$E/ SAV INTEGER Point to start GR C NERR /S$A$V$E/ SAV INTEGER Error processor number C NMESS /S$A$V$E/ SAV INTEGER Length of error message C NOQUAD /S$A$V$E/ SAV LOGICAL Flag, suppress quad model C NPMAX /S$A$V$E/ SAV INTEGER Max number of quad terms C NVARS DUMMY-ARG INTEGER Number of unknowns C RDUM /S$A$V$E/ SAV REAL Dummy variable, error proc C ROPT DUMMY-ARG REAL ADJ-ARY Option data C WA DUMMY-ARG REAL ADJ-ARY Working array C X DUMMY-ARG REAL ADJ-ARY Values of the variables C XMESS /S$A$V$E/ SAV CHAR*128 Hold error message C DATA IFLAG/0/ C***FIRST EXECUTABLE STATEMENT DQED ASSIGN 40 TO IGOOK IF (IFLAG.EQ.0) THEN NOQUAD = .FALSE. LEVEL = 1 IF (MEQUA.LE.0) THEN XMESS = . 'DQED. VALUE OF MEQUA=NO. OF EQUAS. MUST .GT.0. NOW = (I1).' CALL CHRCNT(XMESS,NMESS) NERR = 01 IGO = 9 CALL XERRWV(XMESS,NMESS,NERR,LEVEL,1,MEQUA,IDUM,0,RDUM, . RDUM) ASSIGN 50 TO IGOOK * END IF * IF (NVARS.LE.0) THEN XMESS = . 'DQED. VALUE OF NVARS=NO. OF EQUAS. MUST .GT.0. NOW = (I1).' CALL CHRCNT(XMESS,NMESS) NERR = 02 IGO = 10 CALL XERRWV(XMESS,NMESS,NERR,LEVEL,1,NVARS,IDUM,0,RDUM, . RDUM) ASSIGN 50 TO IGOOK * END IF * IF (MCON.LT.0) THEN XMESS = . 'DQED. VALUE OF MCON=NO. OF EQUAS. MUST .GE.0. NOW = (I1).' CALL CHRCNT(XMESS,NMESS) NERR = 03 IGO = 11 CALL XERRWV(XMESS,NMESS,NERR,LEVEL,1,MCON,IDUM,0,RDUM, . RDUM) ASSIGN 50 TO IGOOK * END IF * DO 10 J = 1,NVARS + MCON I = IND(J) GO TO (10,10,10,10),I * XMESS = . 'DQED. BOUND INDICATORS MUST BE 1-4. NOW I1=J, I2= IND(I1).' CALL CHRCNT(XMESS,NMESS) NERR = 10 IGO = 18 CALL XERRWV(XMESS,NMESS,NERR,LEVEL,2,J,IND(J),0,RDUM,RDUM) ASSIGN 50 TO IGOOK * 10 CONTINUE NPMAX = 05 C LOOK THROUGH THE OPTION ARRAY FOR A CHANGE TO NPMAX, C THE AMOUNT OF ARRAY STORAGE USED FOR THE QUADRATIC PARAMETERS. LP = 1 LPDIFF = 0 20 CONTINUE LP = LP + LPDIFF LPDIFF = 2 KP = IOPT(LP) JP = ABS(KP) IF (JP.EQ.99) THEN IF (KP.GT.0) GO TO 30 END IF C THIS IS THE ONLY OPTION WHERE THE PROCESSING POINTER C MUST BE CHANGED FROM THE VALUE 2. IF (JP.EQ.13) LPDIFF = IOPT(LP+1) C FOUND A CHANGE TO THE ARRAY SIZE FOR THE QUADRATIC MODEL. IF (JP.EQ.15) THEN IF (KP.GT.0) NPMAX = IOPT(LP+1) END IF C SEE IF THE QUADRATIC MODEL IS SUPPRESSED. C THIS REQUIRES LESS STORAGE IN THE USER PROGRAM. IF (JP.EQ.14) THEN IF (KP.GT.0) NOQUAD = IOPT(LP+1) .EQ. 1 END IF * IF (JP.LT.1 .OR. JP.GT.17) THEN C SAW AN OPTION (OR GARBAGE) THAT IS NOT ON THE LIST. XMESS = . 'DQED. INVALID OPTION PROOCESSED. I1=IOPT(*) ENTRY. I2=IOPT(I1).' NERR = 04 IGO = 12 CALL CHRCNT(XMESS,NMESS) CALL XERRWV(XMESS,NMESS,NERR,LEVEL,2,LP,IOPT(LP),0,RDUM, . RDUM) ASSIGN 50 TO IGOOK * END IF * GO TO 20 * 30 CONTINUE END IF * IF (NOQUAD) THEN NALL = MCON + NVARS + 2 * ELSE NALL = MCON + 2*NVARS + NPMAX + 1 END IF C C COMPUTE POINTERS INTO WORK SPACE FOR VARIOUS ARRAYS C REQUIRED IN MAIN PROGRAMS. MDX = 1 MXB = MDX + NALL + 2 IF(MCON.GT.0)MXB=MXB+NALL+2 MB = MXB + NVARS MBB = MB + NVARS MBLB = MBB + NVARS IF(MCON.GT.0)MBLB=MBLB+NALL MBUB = MBLB + NALL IF(MCON.GT.0)MBUB=MBUB+NALL MZP = MBUB + NALL MXP = MZP + MEQUA*NPMAX MQC = MXP + NVARS*NPMAX MWJ = MQC + MAX(MEQUA,NVARS)*NPMAX MPJ = MWJ + (NALL)* (NALL+1) MGR = MPJ + NVARS + 1 MDXL = MGR + NVARS MWA = MDXL + NVARS + NVARS MWLAST = MWA + 9* (MCON+1) + 13* (2*NVARS+NPMAX+1) C MINDB = 3 MIND = MINDB + NALL + NVARS MILAST = MIND + 3* (MCON+1) + 4* (2*NVARS+NPMAX+1) C CHECK LENGTHS OF ARRAYS ONCE PER PROBLEM. IF (IFLAG.EQ.0) THEN IWAAV = IWA(1) IIWAAV = IWA(2) C RETURN THE ACTUAL AMOUNTS OF STORAGE REQD. FOR WA(*), IWA(*). IWA(1) = MWLAST IWA(2) = MILAST IF (IWAAV.LT.MWLAST) THEN XMESS = . 'DQED. WA(*) STORAGE SHORT. I1=AMOUNT NEEDED. I2=AMOUNT GIVEN.' NERR = 05 IGO = 13 CALL CHRCNT(XMESS,NMESS) CALL XERRWV(XMESS,NMESS,NERR,LEVEL,2,MWLAST,IWAAV,0,RDUM, . RDUM) ASSIGN 50 TO IGOOK * END IF * IF (IIWAAV.LT.MILAST) THEN XMESS = . 'DQED. IWA(*) STORAGE SHORT. I1=AMOUNT NEEDED. I2=AMOUNT GIVEN.' NERR = 06 IGO = 14 CALL CHRCNT(XMESS,NMESS) CALL XERRWV(XMESS,NMESS,NERR,LEVEL,2,MILAST,IIWAAV,0,RDUM, . RDUM) ASSIGN 50 TO IGOOK * END IF * IFLAG = 1 END IF * GO TO IGOOK * 40 CONTINUE * CALL DQEDMN(DQEDEV,MEQUA,NVARS,MCON,IND,BL,BU,X,FJAC,LDFJAC,FNORM, . IGO,IOPT,ROPT,IWA(MIND),WA(MWA),WA(MDX),WA(MXB), . WA(MB),WA(MBB),WA(MBLB),WA(MBUB),IWA(MINDB),NPMAX, . WA(MZP),WA(MXP),WA(MQC),MAX(MEQUA,NVARS),WA(MPJ), . WA(MWJ),NALL,WA(MGR),WA(MDXL)) 50 CONTINUE IF (IGO.GT.1) IFLAG = 0 RETURN C TOTAL WORKING STORAGE (FOR DQEDGN) IN WA(*)= C 9*NALL + 4*NVARS = 9*MCON + 13*NVARS. C TOTAL WORKING STORAGE (FOR DQEDGN) IN IWA(*)= C 3*NALL + NV = 3*MCON +4*NVARS. C IN THE ABOVE FORMULA, NVARS (FOR DQEDGN) IS .LE. 3*NVARS C IN TERMS OF THE USER'S VARIABLES. END SUBROUTINE DQEDMN(DQEDEV,MEQUA,NVARS,MCON,IND,BL,BU,X,FJAC,LDFJAC, . FB,IGO,IOPT,ROPT,IWA,WA,DX,XB,B,BB,BLB,BUB,INDB, . NPMAX,ZP,XP,QC,MDQC,PJ,WJ,LDWJ,GR,DXL) C***BEGIN PROLOGUE DQEDMN C***REFER TO DQED C***ROUTINES CALLED DGECO,DGESL,IDAMAX,DNRM2,I1MACH, C DQEDEV,DQEDGN,D1MACH,DAXPY,DSCAL,DCOPY, C DVOUT,DDOT,CHRCNT,XERRWV C***END PROLOGUE DQEDMN C REVISED 861216-1100 C REVISED YYMMDD-HHMM DOUBLE PRECISION BL(*),BU(*),X(*),FJAC(LDFJAC,*) DOUBLE PRECISION DX(*),XB(*),B(*),BB(*) DOUBLE PRECISION BLB(*),BUB(*),ROPT(*),WA(*),PJ(*) DOUBLE PRECISION WJ(LDWJ,*),GR(*),DXL(*) C ARRAYS TO HOLD VALUES FOR 2-ND ORDER TERMS . DOUBLE PRECISION ZP(MEQUA,NPMAX),XP(NVARS,NPMAX),QC(MDQC,NPMAX) C SCALARS: DOUBLE PRECISION AJN,ALB,ALFAC,ALPHA,AUB,BBOOST,BOLD,CHG DOUBLE PRECISION CHGFAC,COLNRM,COND,COSL,COSM,COSQ,DFN,DXNRM DOUBLE PRECISION C1516,FB,FC,FL,PV,PVL,RC,RG,RB,PB,PD DOUBLE PRECISION RCOND,SEMIBG,T,TOLD,TOLF,TOLP,TOLSNR,TOLUSE DOUBLE PRECISION TOLX,TT,TWO,T2,ZERO,ONE,ZN,SS,SC DOUBLE PRECISION D1MACH,DDOT,DNRM2 INTEGER IND(*),INDB(*),IOPT(*),IWA(*) LOGICAL RETREA,TERM,FULNWT,USEQ,NEWBST,PASSB,NOQUAD,NEWOPT LOGICAL LINMOD,REVERS,USEQL,MUSTCN,JACTRI CHARACTER XMESS*128 DATA IFLAG/0/ C--PROCEDURES-- C -NAME------TYPE--------ARGS------CLASS----- C C DGECO 6 SUBROUTINE C DGESL 6 SUBROUTINE C IDAMAX INTEGER 3 FUNCTION C DNRM2 REAL 3 FUNCTION C DQEDEV 6 DUMMY-SUBR C DQEDGN 15 SUBROUTINE C I1MACH INTEGER 1 FUNCTION C D1MACH REAL 1 FUNCTION C DAXPY 6 SUBROUTINE C DSCAL 4 SUBROUTINE C DCOPY 5 SUBROUTINE C SVOUT 4 SUBROUTINE C DDOT REAL 5 FUNCTION C CHRCNT 2 SUBROUTINE C XERRWV 10 SUBROUTINE C DQEDMN: C GLOSSARY OF VARIABLES. NOTATION: C DUMMY-ARG A DUMMY ARGUMENT, THAT IS AN ARGUMENT TO THIS PROG. UNIT. C /S$A$V$E/ SAV DENOTES THAT THIS VARIABLE IS LOCAL TO THE ROUTINE C AND IS SAVED BETWEEN CALLS TO IT. C INTEGER, REAL, DOUBLE PRECISION, LOGICAL, CHARACTER C THE TYPES OF THE VARIABLES. C ADJ-ARR AN ADJUSTABLE ARRAY, THAT IS AN ARGUMENT TO THIS PROG. UNIT. C NAME MEMORY STATUS TYPE ARGUMENT USES AND COMMENTS. C STATUS C ---- ------------- ---- -------- ------------------ C AJN /S$A$V$E/ SAV REAL NORM GRAD VECTOR C ALB /S$A$V$E/ SAV REAL TEMP BOUND VALUE C ALFAC /S$A$V$E/ SAV REAL TRUST REGION FACTOR C ALPHA /S$A$V$E/ SAV REAL TRUST REGION FACTOR C AUB /S$A$V$E/ SAV REAL TEMP BOUND VALUE C B DUMMY-ARG REAL ADJ-ARY CURRENT TRUST BOUNDS C BB DUMMY-ARG REAL ADJ-ARY TRUST BOUNDS AT BEST C BBOOST /S$A$V$E/ SAV REAL FACTOR TO BOOST BOUNDS C BL DUMMY-ARG REAL ADJ-ARY USER LOWER BOUNDS C BLB DUMMY-ARG REAL ADJ-ARY MODEL LOWER BOUNDS C BOLD /S$A$V$E/ SAV REAL TEMP LAST BOUND VALUE C BU DUMMY-ARG REAL ADJ-ARY USER UPPER BOUNDS C BUB DUMMY-ARG REAL ADJ-ARY MODEL UPPER BOUNDS C CHG /S$A$V$E/ SAV REAL TRUST REGION FACTOR C CHGFAC /S$A$V$E/ SAV REAL TRUST REGION FACTOR C COLNRM /S$A$V$E/ SAV REAL TEMP JACOBIAN COL NORM C COND /S$A$V$E/ SAV REAL MAX COND NUMBER, QUAD TERM C COSL /S$A$V$E/ SAV REAL COSINE, GRAD AND LIN STEP C COSM /S$A$V$E/ SAV REAL COSINE, LIN AND QUAD STEP C COSQ /S$A$V$E/ SAV REAL COSINE, GRAD AND QUAD STEP C C1516 /S$A$V$E/ SAV REAL THE CONSTANT 15/16 C DFN /S$A$V$E/ SAV REAL NORM J*DX/NORM DEL F C DX DUMMY-ARG REAL ADJ-ARY THE CHANGE, X=X-DX C DXL DUMMY-ARG REAL ADJ-ARY A CANDIDATE DX, LIN MODEL C DXNRM /S$A$V$E/ SAV REAL NORM OF DX C FB DUMMY-ARG REAL NORM AT THE BEST X C FC /S$A$V$E/ SAV REAL NORM AT THE CURRENT X C FJAC DUMMY-ARG REAL ADJ-ARY JACOBIAN ARRAY C FL /S$A$V$E/ SAV REAL NORM AT THE LAST X C FULNWT /S$A$V$E/ SAV LOGICAL FLAG, TOOK FULL STEP C GR DUMMY-ARG REAL ADJ-ARY GRADIENT VECTOR C I /S$A$V$E/ SAV INTEGER DUMMY LOOP VARIABLE C ICASE /S$A$V$E/ SAV INTEGER GATE VARIABLE C IFLAG /S$A$V$E/ SAV INTEGER INTERNAL FIRST TIME FLAG C IGO DUMMY-ARG INTEGER DIRECT USER, PROGRAM ACTION C NERR /S$A$V$E/ SAV INTEGER ERROR PROCESSOR NUMBER C NEWBST /S$A$V$E/ SAV LOGICAL FLAG, GOT A NEW BEST C NEWOPT /S$A$V$E/ SAV LOGICAL FLAG, SET AN OPTION C NIT /S$A$V$E/ SAV INTEGER ITERS, QUAD MODEL SOLVING C NMESS /S$A$V$E/ SAV INTEGER LENGTH OF ERROR MESSAGE C NOQUAD /S$A$V$E/ SAV LOGICAL FLAG, SUPPRESS QUAD MODEL C NOUT /S$A$V$E/ SAV INTEGER UNIT NUMBER, SAMPLE OUTPUT C NP /S$A$V$E/ SAV INTEGER POTENTIAL QUAD TERMS + 1 C NPMAX DUMMY-ARG INTEGER MAX NUMBER QUAD TERMS C NT /S$A$V$E/ SAV INTEGER MIN(NVARS+1,MEQUA-1) C NTTERM /S$A$V$E/ SAV INTEGER NUMBER QUAD TERMS USED C NV /S$A$V$E/ SAV INTEGER NUMBER OF VARIABLES, MODEL C NVARS DUMMY-ARG INTEGER NUMBER OF USER VARIABLES C ONE /S$A$V$E/ SAV REAL THE NUMBER 1. C PASSB /S$A$V$E/ SAV LOGICAL FLAG, USER GAVE TRUST BNDS C PB /S$A$V$E/ SAV REAL PREDICTED RESIDUAL, AT BEST C PD /S$A$V$E/ SAV REAL 1.5(FB+PB(PB/FB)) - 4PB C PJ DUMMY-ARG REAL ADJ-ARY J**T*F AND J*DX AT X C PV /S$A$V$E/ SAV REAL PREDICTED RESIDUAL, CURRENT C PVL /S$A$V$E/ SAV REAL PREDICTED RESIDUAL, LIN MOD C QC DUMMY-ARG REAL ADJ-ARY QUAD MODEL COEFFICIENTS C RB /S$A$V$E/ SAV REAL ALPHA AFTER NOT A BEST C RC /S$A$V$E/ SAV REAL RETREAT COEFFICIENT C RCOND /S$A$V$E/ SAV REAL RECIPROCAL, QUAD TERM COEFF C RDUM /S$A$V$E/ SAV REAL DUMMY VARIABLE, ERROR PROC C RETREA /S$A$V$E/ SAV LOGICAL FLAG, WILL RETREAT C REVERS /S$A$V$E/ SAV LOGICAL FLAG, REVERSE COMMUNICATION C IGOELM /S$A$V$E/ SAV INTEGER BACK FROM LIN MOD EVALUATE C IGOEQM /S$A$V$E/ SAV INTEGER BACK FROM QUAD MOD EVALUATE C IGOTFC /S$A$V$E/ SAV INTEGER BACK FROM TEST FOR CONVERGE C IGOW /S$A$V$E/ SAV INTEGER DIRECT, SOLVE QUAD MODEL C IND DUMMY-ARG INTEGER ADJ-ARY USER BOUND INDICATORS C INDB DUMMY-ARG INTEGER ADJ-ARY INTERNAL BOUND INDICATORS C IOPT DUMMY-ARG INTEGER ADJ-ARY USER OPTION ARRAY C IPLS /S$A$V$E/ SAV INTEGER C IPRINT /S$A$V$E/ SAV INTEGER WANT OUTPUT IF .GT. 0 C ITERS /S$A$V$E/ SAV INTEGER NUMBER OF ITERATIONS C ITMAX /S$A$V$E/ SAV INTEGER MAX NUMBER OF ITERATIONS C IWA DUMMY-ARG INTEGER ADJ-ARY WORKING ARRAY C J /S$A$V$E/ SAV INTEGER DUMMY LOOP VARIABLE C JACTRI /S$A$V$E/ SAV LOGICAL FLAG, JACOBIAN IS TRIANGLE C JK /S$A$V$E/ SAV INTEGER TEMP LOOP VARIABLE C JP /S$A$V$E/ SAV INTEGER OPTION ARRAY POINTER C K /S$A$V$E/ SAV INTEGER DUMMY LOOP VARIABLE C KL /S$A$V$E/ SAV INTEGER DIRECT STATE OF DQED C KP /S$A$V$E/ SAV INTEGER OPTION ARRAY POINTER C L /S$A$V$E/ SAV INTEGER DUMMY LOOP VARIABLE C LDFJAC DUMMY-ARG INTEGER ROW DIMENSION OF FJAC C LDWJ DUMMY-ARG INTEGER ROW DIMENSION OF WJ C LEVEL /S$A$V$E/ SAV INTEGER ERROR PROC RESPONSE LEVEL C LINMOD /S$A$V$E/ SAV LOGICAL FLAG, SOLVING LINEAR PROBLEM C LK /S$A$V$E/ SAV INTEGER MIN(MEQUA,NVARS+1) C LP /S$A$V$E/ SAV INTEGER OPTION ARRAY POINTER C LPDIFF /S$A$V$E/ SAV INTEGER OPTION ARRAY INCREMENT C MCON DUMMY-ARG INTEGER NUMBER, LINEAR CONSTRAINTS C MCONST /S$A$V$E/ SAV INTEGER NUMBER, MODEL CONSTRAINTS C MDQC DUMMY-ARG INTEGER ROW DIM OF QC(,) C ME /S$A$V$E/ SAV INTEGER NUMBER, MODEL EQUATIONS C MEQUA DUMMY-ARG INTEGER NUMBER, USER EQUATIONS C MK /S$A$V$E/ SAV INTEGER 0 OR MIN(MEQUA,NVARS+NP+1) C MUSTCN /S$A$V$E/ SAV LOGICAL FLAG, MUST TAKE FULL STEP C NALL /S$A$V$E/ SAV INTEGER MCON+NVARS C RG /S$A$V$E/ SAV REAL RETREAT GAUGE C ROPT DUMMY-ARG REAL ADJ-ARY USER PASSED OPTION DATA C SEMIBG /S$A$V$E/ SAV REAL CONSTANT 1.D+10 C T /S$A$V$E/ SAV REAL TEMP VARIABLE C TERM /S$A$V$E/ SAV LOGICAL FLAG, STOP THIS PROBLEM C TOLD /S$A$V$E/ SAV REAL REL TOLERANCE ON CHANGE C TOLF /S$A$V$E/ SAV REAL TOLERANCE ON FUNCTION C TOLP /S$A$V$E/ SAV REAL TOLERANCE FOR MIN FLAG C TOLSNR /S$A$V$E/ SAV REAL TOLERANCE FOR MIN FLAG C TOLUSE /S$A$V$E/ SAV REAL TOLERANCE FOR QUAD MODEL C TOLX /S$A$V$E/ SAV REAL ABS TOLERANCE ON CHANGE C TT /S$A$V$E/ SAV REAL TEMP VARIABLE C TWO /S$A$V$E/ SAV REAL CONSTANT 2. C T2 /S$A$V$E/ SAV REAL REL STEP WITHIN TRUST REG C USEQ /S$A$V$E/ SAV LOGICAL USE QUAD MODEL THIS STEP C USEQL /S$A$V$E/ SAV LOGICAL USED QUAD MODEL LAST STEP C UT DUMMY-ARG REAL ADJ-ARY WORKING ARRAY C WA DUMMY-ARG REAL ADJ-ARY WORKING ARRAY C WJ DUMMY-ARG REAL ADJ-ARY MODEL JACOBIAN ARRAY C X DUMMY-ARG REAL ADJ-ARY SOLUTION ARRAY C XB DUMMY-ARG REAL ADJ-ARY BEST VALUES OF X C XMESS /S$A$V$E/ SAV CHAR*128 TEMP FOR ERROR MESSAGE C XP DUMMY-ARG REAL ADJ-ARY WORKING ARRAY C ZERO /S$A$V$E/ SAV REAL CONSTANT 0. C ZN /S$A$V$E/ SAV REAL NORM J*DX/NORM DF C ZP DUMMY-ARG REAL ADJ-ARY WORKING ARRAY C IF (IFLAG.NE.0) GO TO 50 * LK = MIN(MEQUA,NVARS+1) NT = MIN(NVARS+1,MEQUA-1) ZERO = 0.D0 ONE = 1.D0 TWO = 2.D0 C DO(PROCESS OPTION ARRAY) GO TO 1100 * 10 CONTINUE C DO(INITIALIZE OTHER VALUES) GO TO 1030 * 20 CONTINUE C SET SO X(*)-DX(*) UPDATE WILL BE CORRECT FIRST TIME. DX(1) = ZERO CALL DCOPY(NVARS,DX,0,DX,1) K = 0 C D1MACH(2)="INFINITY" ON THIS MACHINE. FB = D1MACH(2) DXNRM = FB FL = ZERO C C MODEL PROBLEM RESIDUAL. PV = ZERO PVL = ZERO RETREA = .FALSE. FULNWT = .FALSE. TERM = .FALSE. C DO FOREVER 30 CONTINUE ITERS = ITERS + 1 IF (RETREA) THEN C MUST RETREAT TO BEST X VALUES. CALL DCOPY(NVARS,XB,1,X,1) K = 0 KL = -1 FL = FB * ELSE KL = K DO 40 J = 1,NVARS X(J) = X(J) - DX(J) 40 CONTINUE END IF * IF (TERM) THEN IFLAG = 0 C EXIT FOREVER GO TO 840 * END IF * IFLAG = 1 IGO = 1 IF (NP.EQ.NPMAX-1 .AND. NP.LT.NVARS) IGO = -1 IF (REVERS) THEN C THERE ARE TWO POSSIBLE WAYS TO GET FUNCTION AND DERIVATIVE C VALUES FROM THE USER. THE OPTIONAL WAY IS REVERSE COMMUNICATION. C DO(RETURN TO USER PROGRAM UNIT) GO TO 1020 * ELSE C THE NOMINAL WAY IS THROUGH FORWARD COMMUNICATION. CALL DQEDEV(X,FJAC,LDFJAC,IGO,IOPT,ROPT) END IF * 50 CONTINUE C IF IGO HAS BEEN CHANGED BY THE USER TO A VALUE .GT. 1, THEN C THIS IS AN ABORT SIGNAL. STOP UNLESS IT = 99. IF (IGO.EQ.99) THEN C IF IGO = 99 THE EVALUATION CAN'T BE PERFORMED. C WE FORCE A RETREAT AND RESTART IN THIS CASE. DO 70 I = MCON + 1, MCON + MEQUA FJAC(I,NVARS+1) = FC DO 60 J = 1,NVARS FJAC(I,J) = ZERO 60 CONTINUE 70 CONTINUE C A RETREAT IS FORCED TO OCCUR WITH THIS ASSIGNMENT. RETREA = .TRUE. * END IF * FC = DNRM2(MEQUA,FJAC(MCON+1,NVARS+1),1) * IF (IGO.GT.1 .AND. IGO.NE.99) THEN IFLAG = 0 CALL DCOPY(NVARS,XB,1,X,1) C DO(RETURN TO USER PROGRAM UNIT) GO TO 1020 * END IF C SAVE PAST FUNCTION AND VARIABLE VALUES. C DO NOT UPDATE THE PAST POINTS UNLESS THERE IS A C SIGNIFICANT CHANGE IN THE X(*) VALUES. IF (NP.GE.0) THEN IF (DXNRM.GT.TOLUSE*DNRM2(NVARS,X,1)) THEN LP = NVARS IF ( .NOT. NOQUAD) NP = MIN(NP,NPMAX-1,LP) + 1 DO 150 J = NP - 1,1,-1 C SAVE THE PAST VALUES OF THE VARIABLES. C DIFFERENCES ARE LATER COMPUTED USING THESE VALUES. CALL DCOPY(NVARS,XP(1,J),1,XP(1,J+1),1) C SAVE THE PAST FUNCTION VALUES IN ZP( , ). C DIFFERENCES ARE LATER COMPUTED USING THESE VALUES. CALL DCOPY(MEQUA,ZP(1,J),1,ZP(1,J+1),1) 150 CONTINUE END IF * END IF C PUT IN THE PRESENT VALUES OF THE VARIABLES. CALL DCOPY(NVARS,X,1,XP(1,1),1) C PUT IN THE PRESENT VALUES OF THE FUNCTIONS. CALL DCOPY(MEQUA,FJAC(MCON+1,NVARS+1),1,ZP(1,1),1) C THIS STATEMENT HAS THE EFFECT OF A FIRST TIME FLAG. NP = MAX(NP,0) C COMPUTE THE COSINES OF THE PAST MOVES WITH THE MOST CURRENT MOVE. DO 170 L = 2,NP DO 160 J = 1,NVARS QC(J,L) = XP(J,L) - XP(J,1) 160 CONTINUE 170 CONTINUE L = 3 180 CONTINUE C DO WHILE(L.LE.NP) IF ( .NOT. (L.LE.NP)) GO TO 200 C CALCULATE THE DIRECTION COSINES OF THE PAST MOVES. T = DDOT(NVARS,QC(1,2),1,QC(1,L),1) TT = DNRM2(NVARS,QC(1,2),1)*DNRM2(NVARS,QC(1,L),1) IF (TT.GT.ZERO) THEN T = T/TT * ELSE T = ONE END IF * IF (IPRINT.GT.0) THEN WRITE (NOUT, . '('' PAST MOVE NUMBER, COSINE OF MOVE'',I3,2X,F6.2)') L - 2, . T END IF * IF(ABS(T) .GT. .98) THEN C DISCARD PAST INFORMATION ASSOCIATED WITH THIS MOVE IF CLOSE TO C A PAST MOVE. DO 190 J = L,NP - 1 CALL DCOPY(MEQUA,ZP(1,J+1),1,ZP(1,J),1) CALL DCOPY(NVARS,XP(1,J+1),1,XP(1,J),1) CALL DCOPY(NVARS,QC(1,J+1),1,QC(1,J),1) 190 CONTINUE NP = NP - 1 C CYCLE WHILE GO TO 180 * END IF * L = L + 1 C END WHILE GO TO 180 * 200 CONTINUE C COMPUTE FUNCTION DIFFERENCES IN QC( , ). DO 220 J = 1,NP - 1 DO 210 I = 1,MEQUA QC(I,J+1) = ZP(I,J+1) - ZP(I,1) 210 CONTINUE 220 CONTINUE C NOW HAVE F(PAST)-F(CURRENT) IN QC( , ), COLS. 2,...,NP USED. C COMPUTE NORM OF DIFFERENCE OF FUNCTION VALUES. IF (NP.GT.1) THEN DFN = DNRM2(MEQUA,QC(1,2),1) * ELSE DFN = ZERO END IF * DO 240 I = 1,NP - 1 DO 230 J = 1,NVARS C NEXT ADD PRODUCT OF JACOBIAN AND PAST VARIABLE DIFFERENCES. CALL DAXPY(MEQUA,- (XP(J,I+1)-XP(J,1)),FJAC(MCON+1,J),1, . QC(1,I+1),1) 230 CONTINUE 240 CONTINUE C DO FOREVER 250 CONTINUE C COMPUTE THE SYMMETRIC MATRIX WHOSE ENTRIES ARE THE C SQUARES OF THE DOT PRODUCTS OF THE PAST DIRECTIONS. C THIS MATRIX IS REQUIRED TO OBTAIN THE QUADRATIC TERMS C ASSOCIATED WITH INTERPOLATING TO PAST FUNCTION VALUES. DO 280 L = 2,NP DO 270 J = L,NP T = ZERO DO 260 I = 1,NVARS T = T + (XP(I,J)-XP(I,1))* (XP(I,L)-XP(I,1)) 260 CONTINUE WJ(L-1,J-1) = T WJ(J-1,L-1) = T 270 CONTINUE 280 CONTINUE C COMPUTE NORM OF REMAINDER INCLUDING LINEAR TERMS, C USING THE LAST MOVE. USEQ = NP .GT. 1 .AND. .NOT. RETREA ZN = ONE IF (NP.GT.1) THEN ZN = DNRM2(MEQUA,QC(1,2),1) C COMPUTE RATIO OF Z TERMS TO CURRENT F VALUE.. IF (USEQ) USEQ = (ZN .GT.1.D-4*DFN .AND. ZN .LT. DFN*.75D0) . .OR. USEQL IF (DFN.GT.ZERO) THEN ZN = ZN/DFN * ELSE ZN = ONE END IF * IF (IPRINT.GT.0) THEN CALL DVOUT(1,ZN,'('' RATIO OF Z TERM TO PAST DF NORM'')', . 4) END IF C SCALE THE MATRIX (MATRIX := D*MATRIX*D, WHERE D**2 = RECIPROCAL C OF THE DIAGONAL TERMS OF THE MATRIX. DO 290 I = 1,NP - 1 DXL(I) = WJ(I,I) IF (DXL(I).EQ.ZERO) THEN NP = I GO TO 250 * ELSE DXL(I) = ONE/DXL(I) END IF * 290 CONTINUE DO 310 I = 1,NP - 1 DO 300 J = 1,NP - 1 WJ(I,J) = (WJ(I,J)*DXL(I))* (WJ(I,J)*DXL(J)) 300 CONTINUE 310 CONTINUE C USE THE LINPACK ROUTINES DGECO(), DGESL() TO OBTAIN C THE COEFFICIENTS OF THE QUADRATIC TERMS, ONE ROW AT C A TIME. CALL DGECO(WJ,LDWJ,NP-1,IWA,RCOND,WA) IF (IPRINT.GT.0) WRITE (NOUT, . '('' RCOND FROM DGECO() = '',2X, 1PD15.4)') RCOND IF (COND*RCOND.LT.ONE) THEN NP = NP - 1 C CYCLE FOREVER GO TO 250 * END IF * DO 340 I = 1,MEQUA C COPY A ROW OF THE INTERPOLATED DATA TO A WORKING ARRAY. C USE THIS ARRAY TO OBTAIN A ROW OF THE QUADRATIC TERMS. CALL DCOPY(NP-1,QC(I,2),MDQC,WA,1) C SCALE THE RIGHT HAND SIDE DATA. DO 320 J = 1,NP - 1 WA(J) = WA(J)*DXL(J) 320 CONTINUE CALL DGESL(WJ,LDWJ,NP-1,IWA,WA,0) C RESCALE THE SOLUTION DATA. DO 330 J = 1,NP - 1 WA(J) = WA(J)*DXL(J) 330 CONTINUE C THE FOLLOWING SIGN CHANGE COMES FROM A CHANGE C OF SIGN IN THE INNER LOOP MODEL PROBLEM. CALL DSCAL(NP-1,-TWO,WA,1) CALL DCOPY(NP-1,WA,1,QC(I,2),MDQC) 340 CONTINUE END IF * 350 CONTINUE C EXIT FOREVER C END FOREVER C NOW HAVE THE QUADRATIC TERMS COMPUTED. C NEXT WILL TRIANGULARIZE THE JACOBIAN TO SAVE SPACE C WHEN USING THE QUADRATIC MODEL. IF (JACTRI) THEN C CONSTRUCT AND THEN APPLY PLANE ROTATIONS C TO ACHIEVE UPPER TRIANGULAR FORM. THIS LOOP C AFFECTS THE JACOBIAN AND RIGHT HAND SIDE. DO 370 J = 1,NT DO 360 I = J + 1,MEQUA CALL DROTG(FJAC(MCON+J,J),FJAC(MCON+I,J),SC,SS) CALL DROT(NVARS-J+1,FJAC(MCON+J,J+1),LDFJAC, . FJAC(MCON+I,J+1),LDFJAC,SC,SS) C NOW APPLY THE TRANSFORMATION TO THE QUADRATIC TERMS. CALL DROT(NP-1,QC(J,2),MDQC, . QC(I,2),MDQC,SC,SS) FJAC(MCON+I,J) = ZERO 360 CONTINUE 370 CONTINUE C NOW WE FINISH TRIANGULARIZING THE QUADRATIC TERMS. C NOTE THAT THIS DOES NOT AFFECT THE RIGHT HAND SIDE. DO 390 L = 1,NP - 1 DO 395 I=NVARS+L+2,MEQUA CALL DROTG(QC(NVARS+L+1,L+1),QC(I,L+1),SC,SS) CALL DROT(NP-L-1,QC(NVARS+L+1,MIN(L+2,NPMAX)),MDQC, . QC(I,MIN(L+2,NPMAX)),MDQC,SC,SS) QC(I,L+1) = ZERO 395 CONTINUE 390 CONTINUE END IF C COMPUTE CURRENT NORM OF J**T*F(X). DO 400 J = 1,NVARS IF (JACTRI) THEN JK = J * ELSE JK = MEQUA END IF * PJ(J) = DDOT(JK,FJAC(MCON+1,J),1,FJAC(MCON+1,NVARS+1),1) 400 CONTINUE AJN = DNRM2(NVARS,PJ,1) C SAVE J**T*F FOR DIRECTION TESTING WITH LINEAR AND QUADRATIC MOVES. IF (AJN.GT.ZERO) CALL DSCAL(NVARS,ONE/AJN,PJ,1) CALL DCOPY(NVARS,PJ,1,GR,1) NEWBST = FC .LT. FB IF (NEWBST) K = 0 IF (K.EQ.0) THEN PB = ZERO PD = D1MACH(2) C WANT TO POSITION AT BEST X VALUES. IF ( .NOT. RETREA) THEN FB = FC END IF * GO TO (410,430,450),2 - KL * GO TO 470 C CASE 1 410 CONTINUE C IMMEDIATELY GOT A NEW BEST X. RG = ZERO IF (T2.LE.0.25D0) THEN BBOOST = ONE CHG = MAX(4.D0*T2,.1D0) END IF DO 420 J = 1,NVARS BB(J) = CHG*BB(J) 420 CONTINUE T = .25D0/ (ALFAC-1.D0) ALPHA = (ZN+ALFAC*T)/ (ZN+T) ALFAC = 1.5*ALPHA BBOOST = MIN(1.5D0*ALPHA*BBOOST,SEMIBG) GO TO 490 C CASE 2 430 CONTINUE USEQL = .FALSE. RG = ZERO C AT THE INITIAL X. ALFAC = 256.D0 DO 440 J = 1,NVARS IF ( .NOT. PASSB) THEN BB(J) = -X(J) END IF * IF (BB(J).EQ.ZERO) THEN IF (JACTRI) THEN JK = J * ELSE JK = MEQUA END IF * COLNRM = DNRM2(JK,FJAC(MCON+1,J),1) IF (COLNRM.NE.ZERO) THEN BB(J) = DDOT(JK,FJAC(MCON+1,J),1, . FJAC(MCON+1,NVARS+1),1) BB(J) = -MAX(ABS(BB(J))/COLNRM/COLNRM,FC/COLNRM) * ELSE BB(J) = -ONE END IF * END IF * XB(J) = X(J) B(J) = BB(J) 440 CONTINUE ALPHA = ONE BBOOST = 0.5D0 C EXIT IF GO TO 520 C CASE 3 450 CONTINUE C RETREAT TO BEST X. IF (ALFAC.NE.256.D0) THEN ALPHA = MIN(4.D0/ALFAC,.25D0) ALFAC = 1.25D0 * ELSE ALPHA = .25D0*ALPHA END IF * BBOOST = .25D0 USEQL = .FALSE. DO 460 J = 1,NVARS C C THE NEXT LINES HAVE THE EFFECT OF REDUCING THE BOUNDS C AT THE CURRENT BEST TO ABOUT 1/16 THEIR CURRENT VALUES C PROVIDED THE CURRENT BOUNDS ARE RELATIVELY SMALL. IF C THE CURRENT BOUNDS ARE RELATIVELY LARGE, THE BOUNDS AT C THE BEST ARE LEFT ABOUT THE SAME. T = ABS(B(J)) TT = ABS(BB(J)) T = (T+.25D0*TT)/ (T+4.D0*TT) B(J) = T*BB(J) BB(J) = B(J) DX(J) = ZERO 460 CONTINUE C EXIT IF GO TO 520 C CASE OTHER 470 CONTINUE C NOT IMMEDIATELY A BEST X. RB = ZERO DO 480 J = 1,NVARS RB = .125D0*MAX(RB,ABS((XB(J)-X(J))/BB(J))) 480 CONTINUE ALPHA = RB ALFAC = TWO BBOOST = (1.D0+RG)/ (.25D0+2.D0*RG) C END CASE 490 CONTINUE DO 500 J = 1,NVARS DX(J) = XB(J) - X(J) IF (DX(J).EQ.ZERO) THEN B(J) = ALPHA*BB(J) * ELSE XB(J) = X(J) B(J) = SIGN(ALPHA*BB(J),DX(J)) + BBOOST*DX(J) END IF * BB(J) = SIGN(MIN(SQRT(D1MACH(2)),ABS(B(J))),B(J)) 500 CONTINUE * ELSE C MUST MAKE SURE THAT PD GETS SET TO A REASONABLE VALUE. C COMPUTE A GAUGE FOR RETREATING IF DID NOT C GET A NEW BEST. ALPHA = (.5D0*ZN+1.D0)/ (ZN+1.D0) CHG = ALPHA*CHG IF (K.EQ.1) THEN CHG = MIN(CHG,T2) CHG = MAX(CHG,.1D0) PB = PV PD = 1.5D0* (FB+PB* (PB/FB)) - 4.D0*PB END IF * DO 510 J = 1,NVARS B(J) = CHG*B(J) IF (K.EQ.1) BB(J) = B(J) 510 CONTINUE END IF * 520 CONTINUE C DO(TEST FOR CONVERGENCE) ASSIGN 530 TO IGOTFC GO TO 980 * 530 CONTINUE IF (TERM) THEN IFLAG = 0 C EXIT FOREVER GO TO 840 * END IF * K = K + 1 C SOLVE MODEL BOUNDED PROBLEM. DO 590 J = 1,NVARS IF (B(J).LT.ZERO) THEN ALB = B(J) IF (DX(J).EQ.ZERO) THEN C THIS CASE IS REQD. TO AVOID USING BUB(*) AT THE INITIAL PT. AUB = -C1516*ALB * ELSE AUB = MIN(-C1516*ALB,-DX(J)+BUB(J)) END IF * ELSE AUB = B(J) IF (DX(J).EQ.ZERO) THEN ALB = -C1516*AUB * ELSE ALB = MAX(-C1516*AUB,-DX(J)+BLB(J)) END IF * END IF * THIS NEXT CODE, ENDING WITH ***, POINTS THE BOX TOWARDS THE BEST * VALUE OF X WHEN NOT AT A NEW BEST. IF (K.GE.2) THEN IF (XB(J).GT.X(J)) THEN BUB(J) = BUB(J)*.25D0 IF (X(J)-BLB(J).GT.XB(J)) THEN BLB(J) = BLB(J)*.75D0 * ELSE BLB(J) = MIN(BLB(J),.75D0* (X(J)-XB(J))) END IF * ELSE BLB(J) = BLB(J)*.25D0 IF (X(J)-BUB(J).LT.XB(J)) THEN BUB(J) = BUB(J)*.75D0 * ELSE BUB(J) = MAX(BUB(J),.75D0* (X(J)-XB(J))) END IF * END IF * END IF *** C RESTRICT THE STEP FURTHER IF USER GIVES BOUNDS. ICASE = IND(J) GO TO (540,550,560,570),ICASE C CASE 1 540 AUB = MIN(AUB,X(J)-BL(J)) GO TO 580 C CASE 2 550 ALB = MAX(ALB,X(J)-BU(J)) GO TO 580 C CASE 3 560 AUB = MIN(AUB,X(J)-BL(J)) ALB = MAX(ALB,X(J)-BU(J)) GO TO 580 C CASE 4 570 CONTINUE C END CASE 580 BLB(J) = ALB C THIS NEXT LINE IS TO GUARANTEE THAT THE LOWER BOUND C IS .LE. THE UPPER BOUND. AUB = MAX(AUB,ALB) BUB(J) = AUB INDB(J) = 3 590 CONTINUE C COMPUTE JACOBIAN*DX AND COMPARE NORM WITH CURRENT FUNCTION. IF (NP.GT.1) THEN PJ(1) = ZERO CALL DCOPY(LK,PJ,0,PJ,1) DO 600 J = 1,NVARS IF (JACTRI) THEN JK = J * ELSE JK = MEQUA END IF * CALL DAXPY(JK,DX(J),FJAC(MCON+1,J),1,PJ,1) 600 CONTINUE T = DNRM2(LK,PJ,1) C THIS TEST SAYS TO USE THE QUADRATIC MODEL IF C THE LAST STEP IS APPROXIMATELY IN THE NULL SPACE C OF THE JACOBIAN. USEQ = USEQ .OR. T .LT. DFN*0.75D0 IF (DFN.GT.ZERO) THEN DFN = T/DFN * ELSE DFN = ZERO END IF * END IF * IF (IPRINT.GT.0) THEN CALL DVOUT(1,DFN,'('' RATIO OF J*DX NORM TO PAST DF NORM'')', . 4) END IF C CHECK IF QUAD. MODEL IS BEING SUPPRESSED. USEQ = USEQ .AND. .NOT. NOQUAD C START THE PROCESS USING THE LINEAR MODEL. LINMOD = .TRUE. MCONST = MCON CA: DO FOREVER 610 CONTINUE C COMPUTE THE REQUIRED DIMENSIONS OF THE MODEL PROBLEM. IF (LINMOD) THEN MK = 0 ME = MIN(MEQUA,NVARS+1) C SET THE INITIAL VALUES FOR THE LINEAR MODEL PROBLEM. DX(1) = ZERO CALL DCOPY(NVARS,DX,0,DX,1) * ELSE IF (USEQ) THEN MK = MIN(MEQUA,NVARS+NP+1) ME = NVARS + MK CALL DCOPY(MK,FJAC(MCON+1,NVARS+1),1,DX(NVARS+1),1) * ELSE C EXIT FOREVER(A) GO TO 730 * END IF * NV = NVARS + MK C NOTE THAT THE RESIDUALS ARE FREE VARIABLES. DO 620 I = NVARS + 1,NV INDB(I) = 4 620 CONTINUE NIT = 0 C THE JACOBIAN, RIGHT SIDE, QUAD. TERMS ARE AN C UPPER TRAPAZOIDAL DATA ARRAY. THIS WILL MAKE SOLVING C THE MODEL PROBLEM MORE EFFICIENT. C DO FOREVER 630 CONTINUE C CALL A SIMPLIFIED VERSION OF THE ALGORITHM TO SOLVE C THE MODEL PROBLEM. THE FUNCTION AND JACOBIAN C ARE COMPUTED FOR THIS SUBPROBLEM DURING THE REVERSE C COMMUNICATION REQUESTS. CALL DQEDGN(ME,NV,MCONST,INDB,BLB,BUB,DX,WJ,LDWJ,PV,IGOW, . IOPT(IPLS),ROPT,IWA,WA) C CHECK FOR AN ERROR THAT WAS SEEN IN THE LOW-LEVEL C NONLINEAR SOLVER. IF (IGOW.GT.7) THEN IGO = IGOW IFLAG = 0 C DO(RETURN TO USER PROGRAM UNIT) GO TO 1020 * END IF C CLEAR OUT THE WJ(*,*) ARRAY THAT HOLDS C THE JACOBIAN FOR THE INNER LOOP PROBLEM. WJ(1,1) = ZERO CALL DCOPY(LDWJ* (NV+1),WJ,0,WJ,1) IF (USEQ .AND. .NOT. LINMOD) WJ(MCONST+1,NVARS+1) = ONE C PUT IN A UNIT MATRIX FOR THE PARTIALS C WITH RESPECT TO THE RESIDUALS. CALL DCOPY(MK,WJ(MCONST+1,NVARS+1),0,WJ(MCONST+1,NVARS+1),LDWJ+1) DO 640 I = 1,MCON C THE FORM OF THE UPDATE BEING COMPUTED IS X(*)-DX(*). C THE VALUE OF DX IS ITSELF COMPUTED AS THE SOLUTION C TO A NONLINEAR PROBLEM WITH UPDATES OF THE FORM C DX(*)-D(DX)(*). CALL DCOPY(NVARS,FJAC(I,1),LDFJAC,WJ(I,1),LDWJ) WJ(I,NV+1) = FJAC(I,NVARS+1) - DDOT(NVARS,DX,1,WJ(I,1),LDWJ) 640 CONTINUE C SEE IF USER HAS GIVEN GENERAL CONSTRAINTS. C USE THESE CONSTRAINTS TO PLACE EQUIVALENT CONSTRAINTS C ON THE CHANGES BEING COMPUTED. DO 650 J = 1,MCON BLB(NV+J) = BL(J+NVARS) BUB(NV+J) = BU(J+NVARS) INDB(NV+J) = IND(J+NVARS) 650 CONTINUE C DO(EVALUATE LINEAR MODEL) ASSIGN 660 TO IGOELM GO TO 940 * 660 CONTINUE C DO(EVALUATE QUADRATIC MODEL) ASSIGN 670 TO IGOEQM GO TO 850 * 670 CONTINUE IF (IGOW.GT.1) THEN PV = DNRM2(ME,WJ(MCONST+1,NV+1),1) IF (LINMOD) THEN PVL = PV CALL DCOPY(NVARS,DX,1,DXL,1) LINMOD = .FALSE. C CYCLE FOREVER(A) GO TO 610 C IF THE PREDICTED NORM IS GREATER THAN THE CURRENT C RESIDUAL NORM, DROP THE QUADRATIC MODEL AND USE THE C LINEAR MODEL. ELSE IF ((PV.GE.FC) .AND. USEQ) THEN IF (IPRINT.GT.0) WRITE (NOUT, . '('' ABANDON QUAD. MODEL.'')') USEQ = .FALSE. END IF * C EXIT FOREVER(A) GO TO 730 * END IF C FOR EITHER CASE TRANSFER THE JACOBIAN FOR THE MODEL C PROBLEM. THE TRANSPOSE OF THIS MATRIX IS THE PARTIALS C WITH RESPECT TO THE RESIDUALS. DO 680 J = 1,NVARS CALL DCOPY(MK,WJ(MCONST+1,J),1,WJ(MCONST+MK+J,NVARS+1),LDWJ) 680 CONTINUE C NOW UPDATE THE RESIDUALS FOR BOTH SETS OF MODEL EQUATIONS. C IN ROWS 1,...,MK THIS INVOLVES ADDING DX(NVARS+I) TO ROW I. C FOR ROWS MK+1,...,ME THIS REQUIRES ADDING MULTIPLES OF THE C COLS. OF THE TRANSPOSED JACOBIAN. DO 690 I = 1,MK T = DX(NVARS+I) WJ(MCONST+I,NV+1) = WJ(MCONST+I,NV+1) + T 690 CONTINUE C SYMMETRIZE THE SECOND DERIVATIVE MATRIX. THIS C IS NOT REQUIRED WHEN THE MODEL IS LINEAR, BUT IT C DOES NOT HURT THEN. IF (USEQ .AND. .NOT. LINMOD) THEN DO 710 J = 1,NVARS DO 700 I = J,NVARS WJ(MCONST+MK+I,J) = WJ(MCONST+MK+J,I) 700 CONTINUE 710 CONTINUE END IF C COMPUTE RESIDUALS ON THE EQUATIONS K*R = 0. DO 720 J = NVARS + 1,NV CALL DAXPY(NVARS,DX(J),WJ(MCONST+MK+1,J),1, . WJ(MCONST+MK+1,NV+1),1) 720 CONTINUE NIT = NIT + 1 GO TO 630 C END FOREVER C END FOREVER (A) 730 CONTINUE C COMPUTE THE ANGLES BETWEEN THE LINEAR AND QUADRATIC STEP. C TAKE THE ONE, IF THERE IS A CHOICE, CLOSEST TO THE GRADIENT. C IF THE QUADRATIC MOVE IS QUITE CLOSE TO THE GRADIENT, TAKE C THAT MOVE IN PREFERENCE TO THE LINEAR MOVE. COSL = DDOT(NVARS,GR,1,DXL,1) T = DNRM2(NVARS,DXL,1) IF (T.GT.ZERO) COSL = COSL/T COSQ = -ONE COSM = -ONE TT = ZERO IF (USEQ) THEN COSQ = DDOT(NVARS,GR,1,DX,1) TT = DNRM2(NVARS,DX,1) IF (TT.GT.ZERO) COSQ = COSQ/TT C COMPUTE THE COSINE OF THE ANGLE BETWEEN THE QUAD. AND C LINEAR MOVES. IF (T.GT.ZERO .AND. TT.GT.ZERO) COSM = DDOT(NVARS,DX,1,DXL,1)/ . T/TT END IF * IF (IPRINT.GT.0) WRITE (NOUT, .'('' COS OF QUAD. MOVE AND GRAD., COS OF LIN. MOVE AND GRAD., COS . OF EACH MOVE'',/, .'' FLAG FOR TRYING QUAD. MOVE.'',/, .1P3D12.4,L6)') COSQ,COSL,COSM,USEQ IF (IPRINT.GT.0) WRITE (NOUT, .'('' LENGTH OF QUAD., THEN LINEAR MOVES'',/, .1P2D12.4)') TT,T C CHOOSE MOVE PARTIALLY BASED ON ANGLE MOVES MAKE WITH EACH OTHER. USEQ = USEQ .AND. (COSM.GT.ZERO .OR. COSL.LT.COSM) .AND. COSQ .GT. . ZERO USEQL = USEQ * IF ( .NOT. USEQ) THEN PV = PVL CALL DCOPY(NVARS,DXL,1,DX,1) NTTERM = 0 * ELSE NTTERM = NP - 1 END IF C TEST FOR NOISE IN MODEL PROBLEM SOLN. TERM = (PV.GE.FC) .AND. .NOT. RETREA .AND. .NOT. USEQ TERM = TERM .AND. MCON .EQ. 0 TERM = .FALSE. IF (TERM) THEN IF (IPRINT.GT.0) THEN WRITE (NOUT,9021) PV,FC END IF C VALUE MEANS MODEL RES. .GE. NONLINEAR FUNCTION VALUE. IGO = 5 IFLAG = 0 C EXIT FOREVER GO TO 840 * END IF * RC = ZERO IF (PV.GT.PB) RC = 4.D0* (PV-PB)/ (PV+PD) C C IF USING A QUADRATIC MODEL AND RETREATING SEEMS TO BE C NECESSARY, SEE IF RETREATING WOULD BE NEEDED WITH A C LINEAR MODEL. ONLY THEN RETREAT. IF (RC.LE.ONE .OR. .NOT. USEQ) GO TO 750 C DO(EVALUATE LINEAR MODEL) NV = NVARS ASSIGN 740 TO IGOELM GO TO 940 * 740 CONTINUE PVL = DNRM2(MIN(MEQUA,NVARS+1),WJ(MCONST+1,NV+1),1) IF (PVL.GT.PB) RC = 4.D0* (PVL-PB)/ (PVL+PD) 750 CONTINUE RG = MAX(RG,RC) IF (IPRINT.GT.0) THEN WRITE (NOUT,9011) ITERS,FC,PV,RC,AJN,K,KL,FB,ALPHA,BBOOST, . NIT,USEQ,NTTERM WRITE (NOUT,9001) ' X=', (X(J),J=1,NVARS) WRITE (NOUT,9001) ' DX=', (DX(J),J=1,NVARS) WRITE (NOUT,9001) ' B=', (B(J),J=1,NVARS) WRITE (NOUT,9001) ' LB=', (BLB(J),J=1,NALL) WRITE (NOUT,9001) ' UB=', (BUB(J),J=1,NALL) WRITE (NOUT,'('' ///END OF ITERATION.///'')') END IF * RETREA = RC .GT. 1 IF ( .NOT. RETREA) THEN CHG = ONE T2 = ZERO DO 810 J = 1,NVARS BOLD = B(J) T = DX(J)/BOLD ALB = TWO AUB = TWO C IF USER GIVES BOUNDS, AND THESE BOUNDS ARE HIT, C DO NOT DETRACT FROM DECLARING A FULL NEWTON STEP. ICASE = IND(J) GO TO (760,770,780,790),ICASE C CASE 1 760 ALB = (X(J)-BL(J))/BOLD AUB = -SEMIBG GO TO 800 C CASE 2 770 AUB = (X(J)-BU(J))/BOLD ALB = -SEMIBG GO TO 800 C CASE 3 780 ALB = (X(J)-BL(J))/BOLD AUB = (X(J)-BU(J))/BOLD GO TO 800 C CASE 4 790 CONTINUE ALB = -SEMIBG AUB = -SEMIBG C END CASE 800 CONTINUE IF (T.EQ.ONE) THEN T2 = ONE B(J) = BOLD + BOLD CHG = CHG*CHGFAC * ELSE IF (ABS(T).LT..25D0 .AND. DX(J).NE.ZERO) THEN B(J) = SIGN(.25D0*BOLD,DX(J)) + 3.D0*DX(J) * ELSE B(J) = SIGN(BOLD,DX(J)) END IF * END IF C THIS TEST AVOIDS THE USER BOUNDS IN DECLARING A NEWTON STEP. IF (ABS(ALB-T).GE..01D0*ABS(T) .AND. ABS(AUB-T).GE..01D0* . ABS(T)) THEN IF (T.GT.ZERO) THEN T2 = MAX(T2,T) * ELSE T2 = MAX(T2,-T/C1516) END IF * END IF * 810 CONTINUE FULNWT = T2 .LT. .99D0 DXNRM = ABS(DX(IDAMAX(NVARS,DX,1))) C DO BLOCK C TEST FOR SMALL ABSOLUTE CHANGE IN X VALUES. TERM = DXNRM .LT. TOLD .AND. FULNWT IF (TERM) THEN IGO = 6 C VALUE MEANS CHANGE (IN PARAMETERS) WAS SMALL AND A C FULL STEP (NOT HITTING TRUST CONSTRAINTS) TAKEN. C EXIT BLOCK GO TO 820 * END IF * TERM = DXNRM .LT. DNRM2(NVARS,X,1)*TOLX .AND. FULNWT TERM = TERM .AND. (ITERS.GT.1) IF (TERM) THEN IGO = 7 C VALUE MEANS RELATIVE CHANGE IN PARAMETERS WAS SMALL AND A C FULL STEP (NOT HITTING CONSTRAINTS WITH AT LEAST 2 ITERATIONS) C WAS TAKEN. C EXIT BLOCK GO TO 820 * END IF C EXIT IF GO TO 830 C END BLOCK 820 CONTINUE C CYCLE FOREVER GO TO 30 * END IF * 830 CONTINUE FL = FC C END FOREVER GO TO 30 * 840 CONTINUE C DO(RETURN TO USER PROGRAM UNIT) GO TO 1020 C PROCEDURE(EVALUATE QUADRATIC MODEL) 850 CONTINUE C IF THE MODEL IS GENUINELY QUADRATIC, ADD IN THE EXTRA C TERMS AND COMPUTE THE SECOND DERIVATIVE INFORMATION. IF (USEQ .AND. .NOT. LINMOD) THEN C COMPUTE THE DOT PRODUCT OF CURRENT PROPOSED STEP AND C PAST DIRECTIONS REPRESENTED IN THE MODEL. DO 870 L = 1,NP - 1 T = ZERO DO 860 J = 1,NVARS T = T + DX(J)* (XP(J,L+1)-XP(J,1)) 860 CONTINUE PJ(L) = T 870 CONTINUE C STORAGE LAYOUT, WITH K = J**T, OF WJ(*,*). C [J : I : F+R ] C [H : K : K*R ] C ADD IN THE QUADRATIC TERMS FOR THE FUNCTION. DO 880 L = 1,NP - 1 JK = MIN(NVARS+L+1,MEQUA) CALL DAXPY(JK,.5D0*PJ(L)**2,QC(1,L+1),1,WJ(MCONST+1,NV+1), . 1) 880 CONTINUE C ADD THE LINEAR TERMS TO THE INNER LOOP JACOBIAN. DO 900 L = 1,NP - 1 JK = MIN(NVARS+L+1,MEQUA) DO 890 J = 1,NVARS CALL DAXPY(JK,PJ(L)* (XP(J,L+1)-XP(J,1)),QC(1,L+1),1, . WJ(MCONST+1,J),1) 890 CONTINUE 900 CONTINUE C COMPUTE THE UPPER TRIANGULAR PART OF THE SECOND C DERIVATIVE TERMS. DO 930 I = 1,NVARS DO 920 J = I,NVARS DO 910 L = 1,NP - 1 JK = MIN(NVARS+L+1,MEQUA) WJ(MCONST+MK+I,J) = WJ(MCONST+MK+I,J) + . (XP(J,L+1)-XP(J,1))* . (XP(I,L+1)-XP(I,1))* . DDOT(JK,DX(NVARS+1),1,QC(1,L+1), . 1) 910 CONTINUE 920 CONTINUE 930 CONTINUE END IF C END PROCEDURE GO TO IGOEQM C PROCEDURE(EVALUATE LINEAR MODEL) 940 CONTINUE C TRANSFER THE JACOBIAN THAT WOULD RESULT FROM C USING JUST A LINEAR MODEL. DO 950 J = 1,NVARS IF (JACTRI) THEN JK = J * ELSE JK = MEQUA END IF * CALL DCOPY(JK,FJAC(MCON+1,J),1,WJ(MCONST+1,J),1) 950 CONTINUE C TRANSFER THE PRESENT VALUES OF THE FUNCTION. CALL DCOPY(MIN(MEQUA,NVARS+1),FJAC(MCON+1,NVARS+1),1, . WJ(MCONST+1,NV+1),1) C CHANGE SIGN FOR THE MODEL PROBLEM. DO 960 I = 1,MIN(MEQUA,NVARS+1) WJ(MCONST+I,NV+1) = -WJ(MCONST+I,NV+1) 960 CONTINUE C COMPUTE THE LINEAR TERM OF THE MODEL. DO 970 J = 1,NVARS IF (JACTRI) THEN JK = J * ELSE JK = MEQUA END IF * CALL DAXPY(JK,DX(J),WJ(MCONST+1,J),1,WJ(MCONST+1,NV+1),1) 970 CONTINUE C END PROCEDURE GO TO IGOELM C PROCEDURE(TEST FOR CONVERGENCE) 980 CONTINUE TERM = ITERS .GE. ITMAX IF (TERM) THEN IGO = 8 C VALUE MEANS THAT MAX. NUMBER OF ALLOWED ITERATIONS TAKEN. C EXIT PROCEDURE GO TO 1000 * END IF C TEST FOR SMALL FUNCTION NORM. TERM = FC .LE. TOLF C IF HAVE CONSTRAINTS MUST ALLOW AT LEAST ONE MOVE. TERM = TERM .AND. (MCON.EQ.0 .OR. ITERS.GT.1) IF (TERM) THEN IGO = 2 C VALUE MEANS FUNCTION NORM WAS SMALL. C EXIT PROCEDURE GO TO 1000 * END IF C DO(TEST FOR NO CHANGE) GO TO 1010 * 990 CONTINUE TERM = TERM .AND. .NOT. RETREA IF (TERM) THEN IGO = 3 C VALUE MEANS THE FUNCTION IS PROBABLY REACHING A LOCAL MINIMUM C BUT MOVES ARE STILL HITTING TRUST REGION CONSTRAINTS. IF (FULNWT) IGO = 4 C VALUE MEANS THAT FUNCTION IS REACHING A LOCAL MINIMUM C AND MOVES ARE NOT HITTING THE TRUST REGION CONSTRAINTS. IF (IGO.EQ.3) TERM = TERM .AND. .NOT. MUSTCN C EXIT PROCEDURE GO TO 1000 * END IF C END PROCEDURE 1000 CONTINUE GO TO IGOTFC C PROCEDURE(TEST FOR NO CHANGE) 1010 CONTINUE TERM = (ABS(FB-PV).LE.TOLSNR*FB) .AND. (ABS(FC-PV).LE.FB*TOLP) TERM = TERM .AND. (ABS(FC-FL).LE.FB*TOLSNR) TERM = TERM .AND. (ABS(PVL-PV).LE.FB*TOLSNR) C END PROCEDURE GO TO 990 C PROCEDURE(RETURN TO USER PROGRAM UNIT) 1020 CONTINUE RETURN C PROCEDURE(INITIALIZE OTHER VALUES) 1030 CONTINUE C THE NUMBER OF PAST DIFFERENCES USED IN THE QUADRATIC MODEL. NP = 0 C IF NO MORE EQUATIONS THAN VARIABLES, NO NEED TO C PRETRIANGULARIZE THE JACOBIAN MATRIX. JACTRI = MEQUA .GT. NVARS C MAKE SURE THAT VARIABLES SATISFY CONSTRAINTS. C GENERALLY THIS MAY TAKE A CALL TO DBOCLS(). C AS LONG AS THE FUNCTIONS ARE DEFINED AT POINTS C THAT DO NOT SATISFY THE CONSTRAINTS, THE FIRST C ALGORITHM STEP WILL BRING IT ONTO THE CONSTRAINTS. C DO 1090 J = 1,NVARS GO TO (1040,1050,1060,1070),IND(J) C CASE 1 1040 X(J) = MAX(X(J),BL(J)) GO TO 1080 C CASE 2 1050 X(J) = MIN(X(J),BU(J)) GO TO 1080 C CASE 3 1060 X(J) = MAX(X(J),BL(J)) X(J) = MIN(X(J),BU(J)) GO TO 1080 C CASE 4 1070 CONTINUE C END CASE 1080 CONTINUE 1090 CONTINUE ITERS = 0 NALL = MCON + NVARS CHGFAC = TWO** (-ONE/REAL(NVARS)) C1516 = 15.D0/16.D0 SEMIBG = 1.D10 C END PROCEDURE GO TO 20 C PROCEDURE(PROCESS OPTION ARRAY) 1100 CONTINUE IPRINT = 0 C I1MACH(2)=STANDARD OUTPUT UNIT. NOUT = I1MACH(2) C D1MACH(4)=RELPR=MACHINE REL. PREC. T = D1MACH(4) TOLF = SQRT(T) TOLUSE = TOLF TOLX = TOLF TOLD = TOLF TOLSNR = 1.D-5 TOLP = 1.D-5 COND = 30. ITMAX = 75 LEVEL = 1 IPLS = 0 PASSB = .FALSE. NOQUAD = .FALSE. REVERS = .FALSE. MUSTCN = .FALSE. LP = 1 LPDIFF = 0 C DO FOREVER 1110 CONTINUE LP = LP + LPDIFF LPDIFF = 2 KP = IOPT(LP) NEWOPT = KP .GT. 0 JP = ABS(KP) C SEE IF THIS IS THE LAST OPTION.. IF (JP.EQ.99) THEN IF (NEWOPT) THEN C THE POINTER TO THE START OF OPTIONS FOR THE INNER LOOP C SOLVER MUST SATISFY THE REQUIREMENTS FOR THAT OPTION ARRAY. IF (IPLS.EQ.0) IPLS = LP C EXIT FOREVER GO TO 1120 * ELSE LPDIFF = 1 C CYCLE FOREVER GO TO 1110 * END IF * END IF C CHANGE PRINT OPTION. IF (JP.EQ.1) THEN IF (NEWOPT) IPRINT = IOPT(LP+1) C CYCLE FOREVER GO TO 1110 * END IF C SEE IF MAX. NUMBER OF ITERATIONS CHANGING. IF (JP.EQ.2) THEN IF (NEWOPT) ITMAX = IOPT(LP+1) C CYCLE FOREVER GO TO 1110 * END IF C SEE IF BOUNDS FOR THE TRUST REGION ARE BEING PASSED. IF (JP.EQ.3) THEN IF (NEWOPT) THEN CALL DCOPY(NVARS,ROPT(IOPT(LP+1)),1,BB,1) PASSB = .TRUE. END IF C CYCLE FOREVER GO TO 1110 * END IF C CHANGE TOLERANCE ON THE LENGTH OF THE RESIDUALS. IF (JP.EQ.4) THEN IF (NEWOPT) TOLF = ROPT(IOPT(LP+1)) C CYCLE FOREVER GO TO 1110 * END IF C CHANGE TOLERANCE ON THE NORM OF THE RELATIVE C CHANGE TO THE PARAMETERS. IF (JP.EQ.5) THEN IF (NEWOPT) TOLX = ROPT(IOPT(LP+1)) C CYCLE FOREVER GO TO 1110 * END IF C CHANGE TOLERANCE ON ABSOLUTE CHANGE TO THE PARAMETERS. IF (JP.EQ.6) THEN IF (NEWOPT) TOLD = ROPT(IOPT(LP+1)) C CYCLE FOREVER GO TO 1110 * END IF * IF (JP.EQ.7) THEN C CHANGE TOLERANCE FOR RELATIVE AGREEMENT BETWEEN C BEST FUNCTION NORM, LAST FUNCTION NORM AND THE C CURRENT FUNCTION NORM. IF (NEWOPT) TOLSNR = ROPT(IOPT(LP+1)) C CYCLE FOREVER GO TO 1110 * END IF * IF (JP.EQ.8) THEN C CHANGE TOLERANCE FOR AGREEMENT BETWEEN PREDICTED C VALUE OF RESIDUAL NORM AND THE PREVIOUS VALUE OF C THE RESIDUAL NORM. IF (NEWOPT) TOLP = ROPT(IOPT(LP+1)) C CYCLE FOREVER GO TO 1110 * END IF * IF (JP.EQ.9) THEN C CHANGE TOLERANCE SUCH THAT RELATIVE CHANGES IN THE C VALUES OF THE PARAMETERS IMPLY THAT THE PREVIOUS C VALUE OF THE FUNCTION WILL NOT BE USED IN THE C QUADRATIC MODEL. IF (NEWOPT) TOLUSE = ROPT(IOPT(LP+1)) C CYCLE FOREVER GO TO 1110 * END IF * IF (JP.EQ.10) THEN C CHANGE THE LARGEST CONDITION NUMBER TO ALLOW WHEN C SOLVING FOR THE QUADRATIC COEFFICIENTS OF THE MODEL. IF (NEWOPT) COND = ROPT(IOPT(LP+1)) C CYCLE FOREVER GO TO 1110 * END IF C CHANGE THE PRINT LEVEL IN THE ERROR PROCESSOR. IF (JP.EQ.11) THEN IF (NEWOPT) LEVEL = IOPT(LP+1) C CYCLE FOREVER GO TO 1110 * END IF C PASS AN OPTION ARRAY TO THE CONSTRAINED LINEAR SOLVER. C THIS OPTION IS A POINTER TO THE START OF THE OPTION C ARRAY FOR THE SUBPROGRAM. IF (JP.EQ.12) THEN IF (NEWOPT) IPLS = IOPT(LP+1) C CYCLE FOREVER GO TO 1110 * END IF C MOVE THE PROCESSING POINTER BY THE VALUE IN THE C NEXT ENTRY OF THE OPTION ARRAY. THIS DEVICE IS C INCLUDED SO THAT PASSING OPTIONS TO LOWER LEVEL C SUBROUTINES IS EASY TO DO. IF (JP.EQ.13) THEN IF (NEWOPT) LPDIFF = IOPT(LP+1) C CYCLE FOREVER GO TO 1110 * END IF C OPTION TO SUPPRESS USING THE QUADRATIC MODEL, EVER. IF (JP.EQ.14) THEN IF (NEWOPT) NOQUAD = IOPT(LP+1) .EQ. 1 C CYCLE FOREVER GO TO 1110 * END IF C MORE STORAGE WAS GIVEN FOR THE QUADRATIC MODEL ARRAYS. C THIS OPTION WAS PROCESSED BY THE INTERFACE UNIT. C IF(JP.EQ.15)CYCLE FOREVER IF (JP.EQ.15) GO TO 1110 C USE FORWARD COMMUNICATION TO GET THE DERIVATIVES C AND FUNCTION VALUES. IF (JP.EQ.16) THEN IF (NEWOPT) REVERS = IOPT(LP+1) .EQ. 1 C CYCLE FOREVER GO TO 1110 * END IF C FORCE A FULL NEWTON STEP WHEN NEAR THE MINIMUM. C DO NOT ALLOW CONVERGENCE CLAIMS WHEN HITTING BOUNDS. IF (JP.EQ.17) THEN IF (NEWOPT) MUSTCN = IOPT(LP+1) .EQ. 1 C CYCLE FOREVER GO TO 1110 * END IF C SAW AN OPTION (OR GARBAGE) THAT IS NOT ON THE LIST. XMESS = .'DQEDMN. INVALID OPTION PROCESSED. I1=IOPT(*) ENTRY. I2=IOPT(I1).' NERR = 07 IGO = 15 CALL CHRCNT(XMESS,NMESS) CALL XERRWV(XMESS,NMESS,NERR,LEVEL,2,LP,IOPT(LP),0,RDUM,RDUM) IFLAG = 0 C DO(RETURN TO USER PROGRAM UNIT) GO TO 1020 C END FOREVER C END PROCEDURE 1120 CONTINUE GO TO 10 * 9001 FORMAT (A4,1P10D12.4/ (4X,10D12.4)) 9011 FORMAT ('0ITER.=',I3,' FC=',1PD10.4,' PV=',1PD10.4,03X,' RC=', . 1PD10.4,' J**T*F=',1PD10.4,/,' K=',I4,' KL=',I4,/10X,' FB=', . 1PD10.4,' AL=',1PD10.4,' BB=',1PD12.4/' INNER ITERATIONS =',I5, . ' USE QUAD. MODEL?=',L5,' NUM. OF TERMS =',I5) 9021 FORMAT (' MODEL RESIDUAL.GE.CURRENT F. QUITTING.',1P2D12.5) END SUBROUTINE DQEDGN(MEQUA,NVARS,MCON,IND,BL,BU,X,FJAC,LDFJAC,FNORM, . IGO,IOPT,ROPT,IWA,WA) C***BEGIN PROLOGUE DQEDGN C***REFER TO DQED C***ROUTINES CALLED DQEDIP C***END PROLOGUE DQEDGN C DQEDGN: C GLOSSARY OF VARIABLES. NOTATION: C DUMMY-ARG A dummy argument, that is an argument to this prog. unit. C /S$A$V$E/ SAV Denotes that this variable is local to the routine C and is saved between calls to it. C INTEGER, REAL, DOUBLE PRECISION, LOGICAL, CHARACTER C The types of the variables. C ADJ-ARR An adjustable array, that is an argument to this prog. unit. C Name Memory Status Type Argument Uses and comments. C Status C ---- ------------- ---- -------- ------------------ C BL DUMMY-ARG REAL ADJ-ARY Model lower bounds C BU DUMMY-ARG REAL ADJ-ARY Model upper bounds C FJAC DUMMY-ARG REAL ADJ-ARY Model Jacobian array C FNORM DUMMY-ARG REAL Model residual norm C IGO DUMMY-ARG INTEGER direct model action C IND DUMMY-ARG INTEGER ADJ-ARY Model bound indicators C IOPT DUMMY-ARG INTEGER ADJ-ARY Option array C IWA DUMMY-ARG INTEGER ADJ-ARY Working array C LDFJAC DUMMY-ARG INTEGER Row dim of FJAC(*,*) C MB /S$A$V$E/ SAV INTEGER Pointer to B(*) C MBB /S$A$V$E/ SAV INTEGER Pointer to BB(*) C MBLB /S$A$V$E/ SAV INTEGER Pointer to BLB(*) C MBUB /S$A$V$E/ SAV INTEGER Pointer to BLB(*) C MCON DUMMY-ARG INTEGER Number, model constraints C MDX /S$A$V$E/ SAV INTEGER Pointer to DX(*) C MEQUA DUMMY-ARG INTEGER Number, model equations C MINDB /S$A$V$E/ SAV INTEGER Pointer to INDB(*) C MIWA /S$A$V$E/ SAV INTEGER Pointer to IWA(*) C MWA /S$A$V$E/ SAV INTEGER Pointer to WA(*) C MXB /S$A$V$E/ SAV INTEGER Pointer to XB(*) C NALL /S$A$V$E/ SAV INTEGER NVARS+MEQUA C NVARS DUMMY-ARG INTEGER Number, user variables C ROPT DUMMY-ARG REAL ADJ-ARY Option array data C WA DUMMY-ARG REAL ADJ-ARY Working array C C REVISED 870204-1100 C REVISED YYMMDD-HHMM DOUBLE PRECISION BL(*),BU(*),X(*),FJAC(LDFJAC,*),FNORM DOUBLE PRECISION ROPT(*),WA(*) INTEGER IND(*),IOPT(*),IWA(*) C ALLOCATE BLOCKS OF WORKING STORAGE TO LOGICAL ARRAYS. NALL = MCON + NVARS MDX = 1 MXB = MDX + 2*NALL + 2 MB = MXB + NVARS MBB = MB + NVARS MBLB = MBB + NVARS MBUB = MBLB + NALL MWA = MBUB + NALL C MINDB = 1 MIWA = MINDB + NALL + NVARS CALL DQEDIP(MEQUA,NVARS,MCON,IND,BL,BU,X,FJAC,LDFJAC,FNORM,IGO, . IOPT,ROPT,IWA(MIWA),WA(MWA),WA(MDX),WA(MXB),WA(MB), . WA(MBB),WA(MBLB),WA(MBUB),IWA(MINDB)) C THESE DEFINE THE AMOUNT OF STORAGE FOR THE DOUBLE PRECISION AND C INTEGER WORK ARRAYS, WA(*) AND IWA(*). MWA = MWA + 6*NVARS + 5*MCON MIWA = MIWA + 2*NALL C TOTAL WORKING STORAGE IN WA(*)= C 9*NALL + 4*NVARS = 9*MCON + 13*NVARS. C TOTAL WORKING STORAGE IN IWA(*)= C 3*NALL + NV = 3*MCON +4*NVARS. RETURN END SUBROUTINE DQEDIP(MEQUA,NVARS,MCON,IND,BL,BU,X,FJAC,LDFJAC,FB,IGO, . IOPT,ROPT,IWA,WA,DX,XB,B,BB,BLB,BUB,INDB) C***BEGIN PROLOGUE DQEDIP C***REFER TO DQED C***ROUTINES CALLED DBOCLS,IDAMAX,DCOPY,DNRM2,D1MACH, C CHRCNT,XERRWV C***END PROLOGUE DQEDIP C DBOCLS 15 SUBROUTINE C IDAMAX INTEGER 3 FUNCTION C DCOPY 5 SUBROUTINE C DNRM2 REAL 3 FUNCTION C D1MACH REAL 1 FUNCTION C CHRCNT 2 SUBROUTINE C XERRWV 10 SUBROUTINE C DQEDIP: C GLOSSARY OF VARIABLES. NOTATION: C DUMMY-ARG A dummy argument, that is an argument to this prog. unit. C /S$A$V$E/ SAV Denotes that this variable is local to the routine C and is saved between calls to it. C INTEGER, REAL, DOUBLE PRECISION, LOGICAL, CHARACTER C The types of the variables. C ADJ-ARR An adjustable array, that is an argument to this prog. unit. C Name Memory Status Type Argument Uses and comments. C Status C ---- ------------- ---- -------- ------------------ C ALB /S$A$V$E/ SAV REAL C ALFAC /S$A$V$E/ SAV REAL C ALPHA /S$A$V$E/ SAV REAL C AUB /S$A$V$E/ SAV REAL C B DUMMY-ARG REAL ADJ-ARY C BB DUMMY-ARG REAL ADJ-ARY C BBOOST /S$A$V$E/ SAV REAL C BL DUMMY-ARG REAL ADJ-ARY C BLB DUMMY-ARG REAL ADJ-ARY C BOLD /S$A$V$E/ SAV REAL C BU DUMMY-ARG REAL ADJ-ARY C BUB DUMMY-ARG REAL ADJ-ARY C CHG /S$A$V$E/ SAV REAL C CHGFAC /S$A$V$E/ SAV REAL C COLNRM /S$A$V$E/ SAV REAL C C1516 /S$A$V$E/ SAV REAL C DX DUMMY-ARG REAL ADJ-ARY C IPLS /S$A$V$E/ SAV INTEGER C IPRINT /S$A$V$E/ SAV INTEGER C ITERS /S$A$V$E/ SAV INTEGER C ITMAX /S$A$V$E/ SAV INTEGER C IWA DUMMY-ARG INTEGER ADJ-ARY C J /S$A$V$E/ SAV INTEGER C JP /S$A$V$E/ SAV INTEGER C K /S$A$V$E/ SAV INTEGER C KL /S$A$V$E/ SAV INTEGER C KP /S$A$V$E/ SAV INTEGER C LDFJAC DUMMY-ARG INTEGER C LEVEL /S$A$V$E/ SAV INTEGER C EDIT on 950228-1300. REMOVE references to: C LINPRB /S$A$V$E/ SAV INTEGER C LP /S$A$V$E/ SAV INTEGER C LPDIFF /S$A$V$E/ SAV INTEGER C MCON DUMMY-ARG INTEGER C MEQUA DUMMY-ARG INTEGER C MODE /S$A$V$E/ SAV INTEGER C NALL /S$A$V$E/ SAV INTEGER C NERR /S$A$V$E/ SAV INTEGER C NEWBST /S$A$V$E/ SAV LOGICAL C NEWOPT /S$A$V$E/ SAV LOGICAL C NMESS /S$A$V$E/ SAV INTEGER C NOUT /S$A$V$E/ SAV INTEGER C NVARS DUMMY-ARG INTEGER C ONE /S$A$V$E/ SAV REAL C PASSB /S$A$V$E/ SAV LOGICAL C DXNRM /S$A$V$E/ SAV REAL C FB DUMMY-ARG REAL C FC /S$A$V$E/ SAV REAL C FJAC DUMMY-ARG REAL ADJ-ARY C FL /S$A$V$E/ SAV REAL C FULNWT /S$A$V$E/ SAV LOGICAL C GVAL /S$A$V$E/ SAV REAL C ICASE /S$A$V$E/ SAV INTEGER C IFLAG /S$A$V$E/ SAV INTEGER C IGO DUMMY-ARG INTEGER C IGOIOV /S$A$V$E/ SAV INTEGER C IGOPOA /S$A$V$E/ SAV INTEGER C IGOTFC /S$A$V$E/ SAV INTEGER C IGOTNC /S$A$V$E/ SAV INTEGER C IND DUMMY-ARG INTEGER ADJ-ARY C INDB DUMMY-ARG INTEGER ADJ-ARY C IOPT DUMMY-ARG INTEGER ADJ-ARY C PB /S$A$V$E/ SAV REAL C PD /S$A$V$E/ SAV REAL C PV /S$A$V$E/ SAV REAL C RB /S$A$V$E/ SAV REAL C RDUM /S$A$V$E/ SAV REAL C RETREA /S$A$V$E/ SAV LOGICAL C RG /S$A$V$E/ SAV REAL C RNORMC /S$A$V$E/ SAV REAL C ROPT DUMMY-ARG REAL ADJ-ARY C SEMIBG /S$A$V$E/ SAV REAL C T /S$A$V$E/ SAV REAL C TERM /S$A$V$E/ SAV LOGICAL C TOLD /S$A$V$E/ SAV REAL C TOLF /S$A$V$E/ SAV REAL C TOLP /S$A$V$E/ SAV REAL C TOLSNR /S$A$V$E/ SAV REAL C TOLX /S$A$V$E/ SAV REAL C TWO /S$A$V$E/ SAV REAL C T2 /S$A$V$E/ SAV REAL C WA DUMMY-ARG REAL ADJ-ARY C X DUMMY-ARG REAL ADJ-ARY C XB DUMMY-ARG REAL ADJ-ARY C XMESS /S$A$V$E/ SAV CHAR*128 C ZERO /S$A$V$E/ SAV REAL C C REVISED 870204-1100 C REVISED YYMMDD-HHMM DOUBLE PRECISION BL(*),BU(*),X(*),FJAC(LDFJAC,*) DOUBLE PRECISION DX(*),XB(*),B(*),BB(*) DOUBLE PRECISION BLB(*),BUB(*),ROPT(*),WA(*) DOUBLE PRECISION ALB,ALFAC,ALPHA,AUB,BBOOST DOUBLE PRECISION BOLD,CHG,CHGFAC,COLNRM DOUBLE PRECISION DXNRM,C1516,FB,FC,FL DOUBLE PRECISION RG,RB,PB,PD DOUBLE PRECISION GVAL,ONE,PV DOUBLE PRECISION RNORMC,SEMIBG,T,TOLD,TOLF,TOLP DOUBLE PRECISION TOLSNR,TOLX,TWO,T2,ZERO DOUBLE PRECISION D1MACH,DNRM2 REAL RDUM INTEGER IND(*),INDB(*),IOPT(*),IWA(*) CHARACTER XMESS*128 C OPTIONS.. C 1 SET THE PRINTED OUTPUT OFF/ON. REQUIRES TWO ENTRIES C IN IOPT(*). IF IOPT(*+1)=0, NO PRINTING; =1 PRINT. C (THIS PRINTOUT SHOWS VARIOUS QUANTITIES ASSOCIATED C WITH THE NONLINEAR PROBLEM BEING SOLVED. GOOD FOR C DEBUGGING A PROBLEM IN CASE OF TROUBLE. C 2 SET THE MAXIMUM NUMBER OF INTERATIONS. REQUIRES TWO ENTRIES C IN IOPT(*). USE IOPT(*+1)= MAXIMUM NUMBER OF ITERATIONS. C 3 PASS INITIAL BOUNDS FOR THE TRUST REGION TO THE NONLINEAR C SOLVER. REQUIRES TWO ENTRIES IN IOPT(*). USE IOPT(*+1) AS A C POINTER INTO ROPT(*) FOR START OF THE NVARS BOUNDS. C EDIT on 950228-1300: C LOGICAL RETREA,TERM,FULNWT,PASSB,NEWBST,NEWOPT,LINPRB LOGICAL RETREA,TERM,FULNWT,PASSB,NEWBST,NEWOPT DATA IFLAG/0/ C--PROCEDURES-- C -NAME------TYPE--------ARGS------CLASS----- C C DBOCLS 15 SUBROUTINE C IDAMAX INTEGER 3 FUNCTION C DCOPY 5 SUBROUTINE C DNRM2 REAL 3 FUNCTION C D1MACH REAL 1 FUNCTION C CHRCNT 2 SUBROUTINE C XERRWV 10 SUBROUTINE GO TO (50),IFLAG * ZERO = 0.D0 ONE = 1.D0 TWO = 2.D0 C DO(PROCESS OPTION ARRAY) ASSIGN 10 TO IGOPOA GO TO 470 * 10 CONTINUE C DO(INITIALIZE OTHER VALUES) ASSIGN 20 TO IGOIOV GO TO 450 * 20 CONTINUE C SET SO X(*)-DX(*) UPDATE WILL BE CORRECT FIRST TIME. DX(1) = ZERO CALL DCOPY(NVARS,DX,0,DX,1) K = 0 C D1MACH(2)="INFINITY" ON THIS MACHINE. FB = D1MACH(2) DXNRM = FB FL = ZERO C C LINEAR PROBLEM RESIDUAL. PV = ZERO RETREA = .FALSE. FULNWT = .FALSE. TERM = .FALSE. 30 CONTINUE IF ( .NOT. RETREA) ITERS = ITERS + 1 IF (RETREA) THEN C MUST RETREAT TO BEST X VALUES. K = 0 KL = -1 FL = FB CALL DCOPY(NVARS,XB,1,X,1) * ELSE KL = K DO 40 J = 1,NVARS X(J) = X(J) - DX(J) 40 CONTINUE IF (TERM) THEN IFLAG = 0 GO TO 390 * END IF * END IF * IGO = 1 IFLAG = 1 C DO(RETURN TO USER PROGRAM UNIT) GO TO 440 * 50 CONTINUE FC = DNRM2(MEQUA,FJAC(MCON+1,NVARS+1),1) C DO(TEST FOR CONVERGENCE) ASSIGN 60 TO IGOTFC GO TO 400 * 60 CONTINUE IF (TERM) THEN IFLAG = 0 GO TO 390 * END IF * NEWBST = FC .LT. FB .OR. (MCON.GT.0 .AND.ITERS.EQ.2) IF (NEWBST) K = 0 IF (K.EQ.0) THEN RG = ZERO PB = ZERO PD = D1MACH(2) C WANT TO POSITION AT BEST X VALUES. FB = FC C DO CASE(2-KL,3) GO TO (70,90,110),2 - KL * GO TO 120 * 70 CONTINUE C CASE 1 C IMMEDIATELY GOT A NEW BEST X. IF (T2.LE.0.25D0) THEN BBOOST = ONE CHG = MAX(4.D0*T2,.1D0) END IF * DO 80 J = 1,NVARS BB(J) = CHG*BB(J) 80 CONTINUE C THIS CODE FOR ALPHA HAS THE FOLLOWING EFFECT. C IF FC .EQ. PV, THEN ALPHA=ALFAC. C IF FC**2 .EQ. PV*FL THEN ALPHA=2.-1./ALFAC C IF FC**2 IS MUCH LARGER THAN PV*FL, THEN ALPHA=1. T = FC - PV IF (T.EQ.ZERO) THEN ALPHA = ALFAC * ELSE ALPHA = (PV* (FL-PV))/ (FC+PV)/ (ALFAC-ONE) ALPHA = (ABS(T)+ALFAC*ALPHA)/ (ABS(T)+ALPHA) END IF * ALFAC = 1.5D0*ALPHA BBOOST = MIN(1.5D0*ALPHA*BBOOST,SEMIBG) GO TO 140 * 90 CONTINUE C CASE 2 C AT THE INITIAL X. ALFAC = 256.D0 DO 100 J = 1,NVARS IF ( .NOT. PASSB) BB(J) = -X(J) IF (BB(J).EQ.ZERO) THEN COLNRM = DNRM2(MEQUA,FJAC(MCON+1,J),1) IF (COLNRM.NE.ZERO) BB(J) = -FC/COLNRM END IF * IF (BB(J).EQ.ZERO) BB(J) = -ONE XB(J) = X(J) B(J) = BB(J) 100 CONTINUE ALPHA = ONE BBOOST = 0.5D0 C EXIT IF GO TO 170 * 110 CONTINUE C CASE 3 C RETREAT TO BEST X. IF (ALFAC.NE.256.D0) THEN ALPHA = MIN(ONE/ALFAC,.25D0) ALFAC = 1.25D0 * ELSE ALPHA = .25D0*ALPHA END IF * BBOOST = .25D0 GO TO 140 * 120 CONTINUE C CASE OTHER C NOT IMMEDIATELY A BEST X. RB = ZERO DO 130 J = 1,NVARS RB = MAX(RB,ABS((XB(J)-X(J))/BB(J))) 130 CONTINUE ALPHA = RB ALFAC = TWO BBOOST = (8.D0/7.D0+RG)/ (2.D0/7.D0+RG) C END CASE 140 CONTINUE DO 150 J = 1,NVARS DX(J) = XB(J) - X(J) IF (DX(J).EQ.ZERO) THEN B(J) = ALPHA*BB(J) * ELSE XB(J) = X(J) B(J) = SIGN(ALPHA*BB(J),DX(J)) + BBOOST*DX(J) END IF * BB(J) = SIGN(MIN(SQRT(D1MACH(2)),ABS(B(J))),B(J)) 150 CONTINUE * ELSE C COMPUTE A GAUGE FOR RETREATING IF DID NOT C GET A NEW BEST. IF (K.EQ.1) THEN PB = PV * PD = .5D0* (FB-PB)* ((FB+PB)/FB) PD = 1.5D0* (FB+PB* (PB/FB)) - 4.D0*PB END IF * ALPHA = (.5D0*FC+FL)/ (FC+FL) CHG = MIN(ALPHA*CHG,T2) CHG=MAX(CHG,.1D0) DO 160 J = 1,NVARS B(J) = CHG*B(J) IF (K.EQ.1) BB(J) = B(J) 160 CONTINUE END IF * 170 CONTINUE C DO(TEST FOR CONVERGENCE) ASSIGN 180 TO IGOTFC GO TO 400 * 180 CONTINUE IF (TERM) THEN IFLAG = 0 GO TO 390 * END IF * K = K + 1 C SOLVE LINEAR BOUNDED PROBLEM. DO 240 J = 1,NVARS IF (B(J).LT.ZERO) THEN ALB = B(J) IF (DX(J).EQ.ZERO) THEN C THIS CASE IS REQD. TO AVOID USING BUB(*) AT THE INITIAL PT. AUB = -C1516*ALB * ELSE AUB = MIN(-C1516*ALB,-DX(J)+BUB(J)) END IF * ELSE AUB = B(J) IF (DX(J).EQ.ZERO) THEN ALB = -C1516*AUB * ELSE ALB = MAX(-C1516*AUB,-DX(J)+BLB(J)) END IF * END IF C Edit on 950228-1300: C IF(LINPRB) THEN C AUB=D1MACH(2) C ALB=-AUB C END IF C RESTRICT THE STEP FURTHER IF USER GIVES BOUNDS. ICASE = IND(J) C DO CASE(ICASE,4) GO TO (190,200,210,220),ICASE * 190 CONTINUE C CASE 1 AUB = MIN(AUB,X(J)-BL(J)) GO TO 230 * 200 CONTINUE C CASE 2 ALB = MAX(ALB,X(J)-BU(J)) GO TO 230 * 210 CONTINUE C CASE 3 AUB = MIN(AUB,X(J)-BL(J)) ALB = MAX(ALB,X(J)-BU(J)) GO TO 230 * 220 CONTINUE C CASE 4 C END CASE 230 CONTINUE BLB(J) = ALB C THIS NEXT LINE IS TO GUARANTEE THAT THE LOWER BOUND C IS .LE. THE UPPER BOUND. AUB = MAX(AUB,ALB) BUB(J) = AUB INDB(J) = 3 240 CONTINUE C SEE IF USER HAS GIVEN GENERAL CONSTRAINTS. DO 300 J = NVARS + 1,NALL ICASE = IND(J) GVAL = FJAC(J-NVARS,NVARS+1) C DO CASE(ICASE,4) GO TO (250,260,270,280),ICASE * 250 CONTINUE C CASE 1 BLB(J) = - (GVAL-BL(J)) INDB(J) = 1 GO TO 290 * 260 CONTINUE C CASE 2 BUB(J) = - (GVAL-BU(J)) INDB(J) = 2 GO TO 290 * 270 CONTINUE C CASE 3 BLB(J) = - (GVAL-BL(J)) BUB(J) = - (GVAL-BU(J)) INDB(J) = 3 GO TO 290 * 280 CONTINUE C CASE 4 INDB(J) = 4 C END CASE 290 CONTINUE 300 CONTINUE C SOLVE THE LEAST SQUARES PROBLEM WITH BOUNDS AND LINEAR C CONSTRAINTS. THESE BOUNDS CAN COME FROM THE USER OR C THE ALGORITHM. CALL DBOCLS(FJAC,LDFJAC,MCON,MEQUA,NVARS,BLB,BUB,INDB,IOPT(IPLS), . DX,RNORMC,PV,MODE,WA,IWA) IF (IPRINT.GT.0) THEN WRITE (NOUT,9011) ITERS,FC,PV,K,KL,FB,ALPHA,BBOOST WRITE (NOUT,9001) ' +X=', (X(J),J=1,NVARS) WRITE (NOUT,9001) ' +XB=', (XB(J),J=1,NVARS) WRITE (NOUT,9001) ' +DX=', (DX(J),J=1,NALL) WRITE (NOUT,9001) ' + B=', (B(J),J=1,NVARS) WRITE (NOUT,9001) ' +LB=', (BLB(J),J=1,NALL) WRITE (NOUT,9001) ' +UB=', (BUB(J),J=1,NALL) WRITE (NOUT,'('' +///END OF ITERATION.///'')') END IF C TEST FOR NOISE IN LINEAR PROBLEM SOLN. TERM=MCON.EQ.0.AND.(PV.GE.FC) TERM=.FALSE. IF (TERM) THEN IF (IPRINT.GT.0) THEN WRITE (NOUT,9021) PV,FC END IF * CALL DCOPY(NVARS,XB,1,X,1) IGO = 4 IFLAG = 0 GO TO 390 * END IF * RG = MAX(RG, (PV-PB)/PD) IF ( .NOT. RETREA) THEN CHG = ONE T2 = ZERO DO 360 J = 1,NVARS BOLD = B(J) T = DX(J)/BOLD C IF USER GIVES BOUNDS, AND THESE BOUNDS ARE HIT, C DO NOT DETRACT FROM DECLARING A FULL NEWTON STEP. ICASE = IND(J) GO TO (310,320,330,340),ICASE C CASE 1 310 ALB = (X(J)-BL(J))/BOLD AUB = -SEMIBG GO TO 350 C CASE 2 320 AUB = (X(J)-BU(J))/BOLD ALB = -SEMIBG GO TO 350 C CASE 3 330 ALB = (X(J)-BL(J))/BOLD AUB = (X(J)-BU(J))/BOLD GO TO 350 C CASE 4 340 CONTINUE ALB = -SEMIBG AUB = -SEMIBG C END CASE 350 CONTINUE IF (T.EQ.ONE) THEN T2 = ONE B(J) = BOLD + BOLD CHG = CHG*CHGFAC * ELSE IF (ABS(T).LT..25D0 .AND. DX(J).NE.ZERO) THEN B(J) = SIGN(.25D0*BOLD,DX(J)) + 3.D0*DX(J) * ELSE B(J) = SIGN(BOLD,DX(J)) END IF * END IF C THIS TEST AVOIDS THE USER BOUNDS IN DECLARING A NEWTON STEP. IF (ABS(ALB-T).GE..01D0*ABS(T) .AND. ABS(AUB-T).GE..01D0* . ABS(T)) THEN IF (T.GT.ZERO) THEN T2 = MAX(T2,T) * ELSE T2 = MAX(T2,-T/C1516) END IF * END IF * 360 CONTINUE FULNWT = T2 .LT. .99D0 FL = FC DXNRM = ABS(DX(IDAMAX(NVARS,DX,1))) C DO BLOCK C TEST FOR SMALL ABSOLUTE CHANGE IN X VALUES. TERM = DXNRM .LT. TOLD .AND. FULNWT IF (TERM) THEN IGO = 5 C EXIT BLOCK GO TO 370 * END IF * TERM = DXNRM .LT. DNRM2(NVARS,X,1)*TOLX .AND. FULNWT TERM = TERM .AND. (ITERS.GT.1) IF (TERM) THEN IGO = 6 C EXIT BLOCK GO TO 370 * END IF C EXIT IF GO TO 380 C END BLOCK 370 CONTINUE GO TO 30 * END IF * 380 CONTINUE GO TO 30 * 390 CONTINUE C DO(RETURN TO USER PROGRAM UNIT) GO TO 440 C PROCEDURE(TEST FOR CONVERGENCE) 400 CONTINUE C TEST FOR SMALL FUNCTION NORM. TERM = FC .LE. TOLF .OR. TERM C IF HAVE CONSTRAINTS MUST ALLOW AT LEAST ONE MOVE. TERM = TERM .AND. (MCON.EQ.0 .OR. ITERS.GT.1) IF (TERM) THEN IGO = 2 C EXIT PROCEDURE GO TO 420 * END IF C DO(TEST FOR NO CHANGE) ASSIGN 410 TO IGOTNC GO TO 430 * 410 CONTINUE TERM = TERM .AND. .NOT. RETREA IF (TERM) THEN IGO = 3 C EXIT PROCEDURE GO TO 420 * END IF * TERM = ITERS .GE. ITMAX IF (TERM) THEN IGO = 7 END IF C END PROCEDURE 420 CONTINUE GO TO IGOTFC C PROCEDURE(TEST FOR NO CHANGE) 430 CONTINUE T = SQRT(MAX(ZERO, (FL-PV)* (FL+PV))) TERM = (ABS(FB-FC).LE.TOLSNR*FB) .AND. (T.LE.PV*TOLP) TERM = TERM .AND. (ABS(FC-FL).LE.FB*TOLSNR) TERM = TERM .AND. FULNWT C END PROCEDURE GO TO IGOTNC C PROCEDURE(RETURN TO USER PROGRAM UNIT) 440 CONTINUE RETURN C END PROCEDURE C PROCEDURE(INITIALIZE OTHER VALUES) 450 CONTINUE C Edit on 950228-1300: C LINPRB = IGO.EQ.0 ITERS = 0 NALL = MCON + NVARS CHGFAC = TWO** (-ONE/REAL(NVARS)) C1516 = 15.D0/16.D0 SEMIBG = 1.D10 C MAKE SURE THAT VARIABLES SATISFY THE BOUNDS AND CONSTRAINTS. DO 460 J = 1,NALL BLB(J) = BL(J) BUB(J) = BU(J) INDB(J) = IND(J) 460 CONTINUE C END PROCEDURE GO TO IGOIOV C PROCEDURE(PROCESS OPTION ARRAY) 470 CONTINUE IPRINT = 0 C I1MACH(2)=STANDARD OUTPUT UNIT. NOUT = I1MACH(2) C D1MACH(4)=RELPR=MACHINE REL. PREC. T = D1MACH(4) TOLF = T TOLX = TOLF TOLD = TOLF TOLSNR = 1.D-3 TOLP = 1.D-3 ITMAX = 18 PASSB = .FALSE. LEVEL = 1 IPLS = 0 LPDIFF = 0 LP = 1 480 CONTINUE LP = LP + LPDIFF LPDIFF = 2 KP = IOPT(LP) NEWOPT = KP .GT. 0 JP = ABS(KP) C SEE IF THIS IS THE LAST OPTION.. IF (JP.EQ.99) THEN IF (NEWOPT) THEN C THE POINTER TO THE START OF OPTIONS FOR THE LINEAR C SOLVER MUST SATISFY THE REQUIREMENTS FOR THAT OPTION ARRAY. IF (IPLS.EQ.0) IPLS = LP GO TO 490 * ELSE LPDIFF = 1 GO TO 480 * END IF * END IF C CHANGE PRINT OPTION. IF (JP.EQ.1) THEN IF (NEWOPT) IPRINT = IOPT(LP+1) GO TO 480 * END IF C SEE IF MAX. NUMBER OF ITERATIONS CHANGING. IF (JP.EQ.2) THEN IF (NEWOPT) ITMAX = IOPT(LP+1) GO TO 480 * END IF C SEE IF BOUNDS FOR THE TRUST REGION ARE BEING PASSED. IF (JP.EQ.3) THEN IF (NEWOPT) THEN CALL DCOPY(NVARS,ROPT(IOPT(LP+1)),1,BB,1) PASSB = .TRUE. END IF * GO TO 480 * END IF C CHANGE TOLERANCE ON THE LENGTH OF THE RESIDUALS. IF (JP.EQ.4) THEN IF (NEWOPT) TOLF = ROPT(IOPT(LP+1)) GO TO 480 * END IF C CHANGE TOLERANCE ON THE NORM OF THE RELATIVE C CHANGE TO THE PARAMETERS. IF (JP.EQ.5) THEN IF (NEWOPT) TOLX = ROPT(IOPT(LP+1)) GO TO 480 * END IF C CHANGE TOLERANCE ON ABSOLUTE CHANGE TO THE PARAMETERS. IF (JP.EQ.6) THEN IF (NEWOPT) TOLD = ROPT(IOPT(LP+1)) GO TO 480 * END IF * IF (JP.EQ.7) THEN C CHANGE TOLERANCE FOR RELATIVE AGREEMENT BETWEEN C BEST FUNCTION NORM, LAST FUNCTION NORM AND THE C CURRENT FUNCTION NORM. IF (NEWOPT) TOLSNR = ROPT(IOPT(LP+1)) GO TO 480 * END IF * IF (JP.EQ.8) THEN C CHANGE TOLERANCE FOR AGREEMENT BETWEEN PREDICTED C VALUE OF RESIDUAL NORM AND THE PREVIOUS VALUE OF C THE RESIDUAL NORM. IF (NEWOPT) TOLP = ROPT(IOPT(LP+1)) GO TO 480 * END IF C CHANGE THE PRINT LEVEL IN THE ERROR PROCESSOR. IF (JP.EQ.9) THEN IF (NEWOPT) LEVEL = IOPT(LP+1) GO TO 480 * END IF C PASS AN OPTION ARRAY TO THE CONSTRAINED LINEAR SOLVER. C THIS OPTION IS A POINTER TO THE START OF THE OPTION C ARRAY FOR THE SUBPROGRAM. IF (JP.EQ.10) THEN IF (NEWOPT) IPLS = IOPT(LP+1) GO TO 480 * END IF C MOVE THE PROCESSING POINTER BY THE VALUE IN THE C NEXT ENTRY OF THE OPTION ARRAY. THIS DEVICE IS C INCLUDED SO THAT PASSING OPTIONS TO LOWER LEVEL C SUBROUTINES IS EASY TO DO. IF (JP.EQ.11) THEN IF (NEWOPT) LPDIFF = IOPT(LP+1) GO TO 480 * END IF C SAW AN OPTION (OR GARBAGE) THAT IS NOT ON THE LIST. XMESS = .'DQEDIP. INVALID OPTION PROCESSED. I1=IOPT(*) ENTRY. I2=IOPT(I1).' NERR = 08 IGO = 16 CALL CHRCNT(XMESS,NMESS) CALL XERRWV(XMESS,NMESS,NERR,LEVEL,2,LP,IOPT(LP),0,RDUM,RDUM) IFLAG = 0 C DO(RETURN TO USER PROGRAM UNIT) GO TO 440 C END PROCEDURE 490 CONTINUE GO TO IGOPOA * 9001 FORMAT (A4,1P10D12.4/ (4X,10D12.4)) 9011 FORMAT ('0+ITER.=',I3,' FC=',1PD10.4,' PV=',1PD10.4,' K=',I4, . ' KL=',I4,' FB=',1PD10.4/12X,'AL=',1PD14.4,' BB=',1PD14.4) 9021 FORMAT (' LINEAR RESIDUAL.GE.CURRENT F. QUITTING.',1P2D12.5) END SUBROUTINE DQEDEV(X,FJ,LDFJ,IGO,IOPT,ROPT) C***BEGIN PROLOGUE DQEDEV C***REFER TO DQED C***ROUTINES CALLED XERROR,CHRCNT C***END PROLOGUE DQEDEV C REVISED 870204-1100 C REVISED YYMMDD-HHMM C THIS IS THE SUBPROGRAM FOR EVALUATING THE FUNCTIONS C AND DERIVATIVES 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 IF(IGO.NE.0) THEN C STEP 1. PLACE THE PARTIALS OF THE I-TH CONSTRAINT C FUNCTION WITH REPECT TO VARIABLE J IN THE C ARRAY FJ(I,J), I=1,...,MCON, J=1,...,NVARS. C END IF C STEP 2. PLACE THE VALUES OF THE I-TH CONSTRAINT C EQUATION INTO FJ(I,NVARS+1). C IF(IGO.NE.0) THEN 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 END IF C STEP 4. PLACE THE VALUE OF THE I-TH LEAST SQUARES C EQUATION INTO FJ(MCON+I,NVARS+1). C STEP 5. RETURN TO THE CALLING PROGRAM UNIT. C DQEDEV: C GLOSSARY OF VARIABLES. NOTATION: C DUMMY-ARG A dummy argument, that is an argument to this prog. unit. C /S$A$V$E/ SAV Denotes that this variable is local to the routine C and is saved between calls to it. C INTEGER, REAL, DOUBLE PRECISION, LOGICAL, CHARACTER C The types of the variables. C ADJ-ARR An adjustable array, that is an argument to this prog. unit. C Name Memory Status Type Argument Uses and comments. C Status C ---- ------------- ---- -------- ------------------ C FJ DUMMY-ARG REAL ADJ-ARY C IGO DUMMY-ARG INTEGER C IOPT DUMMY-ARG INTEGER ADJ-ARY C LDFJ DUMMY-ARG INTEGER C NERR INTEGER C NMESS INTEGER C ROPT DUMMY-ARG REAL ADJ-ARY C X DUMMY-ARG REAL ADJ-ARY C XMESS CHAR*128 DOUBLE PRECISION FJ(LDFJ,*),X(*),ROPT(*) INTEGER IGO,IOPT(*) CHARACTER XMESS*128 XMESS = . 'DQED. THE EVALUATOR PROGRAM DQEDEV MUST BE WRITTEN BY THE USER.' NERR = 09 IGO = 17 C C THE INTENT HERE IS THAT THE EVALUATOR WILL NOT RETURN C FROM THE ERROR PROCESSOR CALL. CALL CHRCNT(XMESS,NMESS) CALL XERROR(XMESS,NMESS,NERR,01) RETURN END SUBROUTINE DBOCLS(W,MDW,MCON,MROWS,NCOLS,BL,BU,IND,IOPT,X,RNORMC, * RNORM,MODE,RW,IW) C***BEGIN PROLOGUE DBOCLS C***DATE WRITTEN 821220 (YYMMDD) C***REVISION DATE 880722 (YYMMDD) C***CATEGORY NO. K1A2A,G2E,G2H1,G2H2 C***KEYWORDS BOUNDS,CONSTRAINTS,INEQUALITY,LEAST SQUARES,LINEAR C***AUTHOR HANSON, R. J., SNLA C***PURPOSE Solve the bounded and constrained least squares C problem consisting of solving the equation C E*X = F (in the least squares sense) C subject to the linear constraints C C*X = Y. C***DESCRIPTION C C This subprogram solves the bounded and constrained least squares C problem. The problem statement is: C C Solve E*X = F (least squares sense), subject to constraints C C*X=Y. C C In this formulation both X and Y are unknowns, and both may C have bounds on any of their components. This formulation C of the problem allows the user to have equality and inequality C constraints as well as simple bounds on the solution components. C C This constrained linear least squares subprogram solves E*X=F C subject to C*X=Y, where E is MROWS by NCOLS, C is MCON by NCOLS. C C The user must have dimension statements of the form C C DIMENSION W(MDW,NCOLS+MCON+1), BL(NCOLS+MCON), BU(NCOLS+MCON), C * X(2*(NCOLS+MCON)+2+NX), RW(6*NCOLS+5*MCON) C INTEGER IND(NCOLS+MCON), IOPT(17+NI), IW(2*(NCOLS+MCON)) C C (here NX=number of extra locations required for the options; NX=0 C if no options are in use. Also NI=number of extra locations C for options 1-9. C C INPUT C ----- C C ------------------------- C W(MDW,*),MCON,MROWS,NCOLS C ------------------------- C The array W contains the (possibly null) matrix [C:*] followed by C [E:F]. This must be placed in W as follows: C [C : *] C W = [ ] C [E : F] C The (*) after C indicates that this data can be undefined. The C matrix [E:F] has MROWS rows and NCOLS+1 columns. The matrix C is C placed in the first MCON rows of W(*,*) while [E:F] C follows in rows MCON+1 through MCON+MROWS of W(*,*). The vector F C is placed in rows MCON+1 through MCON+MROWS, column NCOLS+1. The C values of MDW and NCOLS must be positive; the value of MCON must C be nonnegative. An exception to this occurs when using option 1 C for accumulation of blocks of equations. In that case MROWS is an C OUTPUT variable only, and the matrix data for [E:F] is placed in C W(*,*), one block of rows at a time. See IOPT(*) contents, option C number 1, for further details. The row dimension, MDW, of the C array W(*,*) must satisfy the inequality: C C If using option 1, C MDW .ge. MCON + max(max. number of C rows accumulated, NCOLS) C C If using option 8, MDW .ge. MCON + MROWS. C Else, MDW .ge. MCON + max(MROWS, NCOLS). C C Other values are errors, but this is checked only when using C option=2. The value of MROWS is an output parameter when C using option number 1 for accumulating large blocks of least C squares equations before solving the problem. C See IOPT(*) contents for details about option 1. C C ------------------ C BL(*),BU(*),IND(*) C ------------------ C These arrays contain the information about the bounds that the C solution values are to satisfy. The value of IND(J) tells the C type of bound and BL(J) and BU(J) give the explicit values for C the respective upper and lower bounds on the unknowns X and Y. C The first NVARS entries of IND(*), BL(*) and BU(*) specify C bounds on X; the next MCON entries specify bounds on Y. C C 1. For IND(J)=1, require X(J) .ge. BL(J); C IF J.gt.NCOLS, Y(J-NCOLS) .ge. BL(J). C (the value of BU(J) is not used.) C 2. For IND(J)=2, require X(J) .le. BU(J); C IF J.gt.NCOLS, Y(J-NCOLS) .le. BU(J). C (the value of BL(J) is not used.) C 3. For IND(J)=3, require X(J) .ge. BL(J) and C X(J) .le. BU(J); C IF J.gt.NCOLS, Y(J-NCOLS) .ge. BL(J) and C Y(J-NCOLS) .le. BU(J). C (to impose equality constraints have BL(J)=BU(J)= C constraining value.) C 4. For IND(J)=4, no bounds on X(J) or Y(J-NCOLS) are required. C (the values of BL(J) and BU(J) are not used.) C C Values other than 1,2,3 or 4 for IND(J) are errors. In the case C IND(J)=3 (upper and lower bounds) the condition BL(J) .gt. BU(J) C is an error. The values BL(J), BU(J), J .gt. NCOLS, will be C changed. Significant changes mean that the constraints are C infeasible. (Users must make this decision themselves.) C The new values for BL(J), BU(J), J .gt. NCOLS, define a C region such that the perturbed problem is feasible. If users C know that their problem is feasible, this step can be skipped C by using option number 8 described below. C C ------- C IOPT(*) C ------- C This is the array where the user can specify nonstandard options C for DBOCLS( ). Most of the time this feature can be ignored by C setting the input value IOPT(1)=99. Occasionally users may have C needs that require use of the following subprogram options. For C details about how to use the options see below: IOPT(*) CONTENTS. C C Option Number Brief Statement of Purpose C ------ ------ ----- --------- -- ------- C 1 Return to user for accumulation of blocks C of least squares equations. The values C of IOPT(*) are changed with this option. C The changes are updates to pointers for C placing the rows of equations into position C for processing. C 2 Check lengths of all arrays used in the C subprogram. C 3 Column scaling of the data matrix, [C]. C [E] C 4 User provides column scaling for matrix [C]. C [E] C 5 Provide option array to the low-level C subprogram DBOLS( ). C {Provide option array to the low-level C subprogram DBOLSM( ) by imbedding an C option array within the option array to C DBOLS(). Option 6 is now disabled.} C 7 Move the IOPT(*) processing pointer. C 8 Do not preprocess the constraints to C resolve infeasibilities. C 9 Do not pretriangularize the least squares matrix. C 99 No more options to change. C C ---- C X(*) C ---- C This array is used to pass data associated with options 4,5 and C 6. Ignore this parameter (on input) if no options are used. C Otherwise see below: IOPT(*) CONTENTS. C C C OUTPUT C ------ C C ----------------- C X(*),RNORMC,RNORM C ----------------- C The array X(*) contains a solution (if MODE .ge.0 or .eq.-22) for C the constrained least squares problem. The value RNORMC is the C minimum residual vector length for the constraints C*X - Y = 0. C The value RNORM is the minimum residual vector length for the C least squares equations. Normally RNORMC=0, but in the case of C inconsistent constraints this value will be nonzero. C The values of X are returned in the first NVARS entries of X(*). C The values of Y are returned in the last MCON entries of X(*). C C ---- C MODE C ---- C The sign of MODE determines whether the subprogram has completed C normally, or encountered an error condition or abnormal status. A C value of MODE .ge. 0 signifies that the subprogram has completed C normally. The value of mode (.ge. 0) is the number of variables C in an active status: not at a bound nor at the value zero, for C the case of free variables. A negative value of MODE will be one C of the cases (-57)-(-41), (-37)-(-22), (-19)-(-2). Values .lt. -1 C correspond to an abnormal completion of the subprogram. These C error messages are in groups for the subprograms DBOCLS(), C DBOLSM(), and DBOLS(). An approximate solution will be returned C to the user only when max. iterations is reached, MODE=-22. C C ----------- C RW(*),IW(*) C ----------- C These are working arrays. (normally the user can ignore the C contents of these arrays.) C C IOPT(*) CONTENTS C ------- -------- C The option array allows a user to modify some internal variables C in the subprogram without recompiling the source code. A central C goal of the initial software design was to do a good job for most C people. Thus the use of options will be restricted to a select C group of users. The processing of the option array proceeds as C follows: a pointer, here called LP, is initially set to the value C 1. At the pointer position the option number is extracted and C used for locating other information that allows for options to be C changed. The portion of the array IOPT(*) that is used for each C option is fixed; the user and the subprogram both know how many C locations are needed for each option. The value of LP is updated C for each option based on the amount of storage in IOPT(*) that is C required. A great deal of error checking is done by the C subprogram on the contents of the option array. Nevertheless it C is still possible to give the subprogram optional input that is C meaningless. For example option 4 uses the locations C X(NCOLS+IOFF),...,X(NCOLS+IOFF+NCOLS-1) for passing scaling data. C The user must manage the allocation of these locations. C C 1 C - C This option allows the user to solve problems with a large number C of rows compared to the number of variables. The idea is that the C subprogram returns to the user (perhaps many times) and receives C new least squares equations from the calling program unit. C Eventually the user signals "that's all" and a solution is then C computed. The value of MROWS is an output variable when this C option is used. Its value is always in the range 0 .le. MROWS C .le. NCOLS+1. It is the number of rows after the C triangularization of the entire set of equations. If LP is the C processing pointer for IOPT(*), the usage for the sequential C processing of blocks of equations is C C C IOPT(LP)=1 C Move block of equations to W(*,*) starting at C the first row of W(*,*). C IOPT(LP+3)=# of rows in the block; user defined C C The user now calls DBOCLS( ) in a loop. The value of IOPT(LP+1) C directs the user's action. The value of IOPT(LP+2) points to C where the subsequent rows are to be placed in W(*,*). Both of C these values are first defined in the subprogram. The user C changes the value of IOPT(LP+1) (to 2) as a signal that all of C the rows have been processed. C C C .