C C C********** COPYRIGHT 1996 EDWIN H. KAUFMAN JR., DAVID J. LEEMING, C********** GERALD D. TAYLOR C********** THE AUTHORS GRATEFULLY ACKNOWLEDGE THE ASSISTANCE OF C********** CENTRAL MICHIGAN UNIVERSITY, THE UNIVERSITY OF VICTORIA C********** (CANADA), AND COLORADO STATE UNIVERSITY. C********** PERMISSION TO USE, COPY, MODIFY, AND DISTRIBUTE THIS C********** SOFTWARE FOR ANY PURPOSE WITHOUT FEE IS HEREBY GRANTED, C********** PROVIDED THAT THIS ENTIRE NOTICE IS INCLUDED IN ANY C********** SOFTWARE WHICH IS OR INCLUDES A COPY OR MODIFICATION OF C********** THIS SOFTWARE AND IN ALL COPIES OF THE SUPPORTING C********** DOCUMENTATION FOR SUCH SOFTWARE. C********** THIS SOFTWARE IS PROVIDED "AS IS", WITHOUT EXPRESS OR C********** IMPLIED WARRANTY. IN PARTICULAR, NEITHER THE AUTHORS NOR C********** THEIR UNIVERSITIES MAKE ANY REPRESENTATION OR WARRANTY OF C********** ANY KIND CONCERNING THE MERCHANTIBILITY OF THIS SOFTWARE C********** OR ITS FITNESS FOR ANY PARTICULAR PURPOSE. C C C WELCOME TO CONMAX, VERSION 1.7 C C THIS PACKAGE WAS LAST REVISED ON DECEMBER 4, 1996. C C SUCCESS STORIES AND BUG REPORTS MAY BE SENT TO EDWIN H. KAUFMAN, JR. C BY E-MAIL AT 32ZEJ7N@CMUVM.CSV.CMICH.EDU C ALTHOUGH WE DO NOT PROMISE TO ATTEMPT TO FIX BUGS, WE MAY LOOK AT C THEM AS TIME PERMITS. C C THIS PACKAGE CONSISTS OF THREE PARTS: C C (I) THE USERS GUIDE YOU ARE READING, C C (II) A SAMPLE DRIVER PROGRAM AND SUBROUTINE FNSET, AND C C (III) THE CONMAX SUBPROGRAMS. C C THIS PACKAGE IS SELF-CONTAINED. ONCE THE USER HAS SET (ONCE) THE C THREE MACHINE DEPENDENT CONSTANTS AND THE PRECISION AS DESCRIBED IN (A) C BELOW, THE PACKAGE CAN BE COMPILED AND EXECUTED TO RUN THE EXAMPLE IN C THE SAMPLE DRIVER PROGRAM AND SUBROUTINE FNSET. DETAILED INSTRUCTIONS C FOR APPLYING THE PROGRAM TO OTHER PROBLEMS ARE GIVEN IN THIS USERS GUIDE. C C CONMAX CONSISTS OF TWO PROGRAMS FOR SOLVING THE PROBLEM C C MINIMIZE W C C SUBJECT TO C C ABS(FUN(I) - CONFUN(I,1)) .LE. W IF ICNTYP(I)=2, C C CONFUN(I,1) .LE. W IF ICNTYP(I)=1, C C CONFUN(I,1) .LE. 0.0 IF ICNTYP(I)=-1 OR -2, C C WHERE IF ICNTYP(I)=-1 THE CONSTRAINT IS LINEAR (I.E. THE LEFT SIDE C CONSISTS OF A LINEAR COMBINATION OF THE PARAMETERS IN THE VECTOR ARRAY C PARAM PLUS A CONSTANT) AND IF ICNTYP(I)=-2 THE CONSTRAINT MAY BE C NONLINEAR. C C THUS WE ARE CHOOSING THE PARAMETERS TO MINIMIZE THE LEFT SIDES OF C THE TYPE 2 AND TYPE 1 (I.E. PRIMARY) CONSTRAINTS SUBJECT TO THE C TYPE -1 AND TYPE -2 (I.E. STANDARD) CONSTRAINTS. C C IF THERE ARE NO PRIMARY CONSTRAINTS THE PROGRAM WILL ATTEMPT TO C FIND A FEASIBLE POINT (THAT IS, A VECTOR PARAM FOR WHICH THE TYPE -1 C AND TYPE -2 CONSTRAINTS ARE SATISFIED WITHIN A TOLERANCE TOLCON C DESCRIBED BELOW) WHICH IS RELATIVELY CLOSE TO THE USER SUPPLIED INITIAL C APPROXIMATION, THEN RETURN. C C TO USE THE PACKAGE THE USER MUST DO THREE THINGS: C C (A) SET THE THREE MACHINE DEPENDENT CONSTANTS AND THE PRECISION. C C IF DOUBLE PRECISION IS DESIRED AND THE NETLIB-SUPPLIED SUBPROGRAMS C I1MACH, D1MACH, AND R1MACH ARE BEING USED, THEN THIS HAS BEEN TAKEN C CARE OF ALREADY. C C IF DOUBLE PRECISION IS DESIRED BUT THE NETLIB-SUPPLIED SUBPROGRAMS C I1MACH, D1MACH, AND R1MACH ARE NOT BEING USED, THEN REPLACE THE CS C IN COLUMN 1 OF THE SUBPROGRAMS I1MACH AND D1MACH INCLUDED IN THIS C PACKAGE BY BLANKS (EXCEPT REPLACE CC BY C) AND ADJUST THE ASSIGNMENT C STATEMENTS IN THESE SUBPROGRAMS TO SET C C I1MACH(1) = (THE INPUT UNIT NUMBER FOR THE MACHINE BEING USED), C C I1MACH(2) = (THE OUTPUT UNIT NUMBER FOR THE MACHINE BEING USED), C C D1MACH(3) = B**(-ITT), WHERE B IS THE BASE FOR FLOATING POINT NUMBERS C AND ITT IS THE NUMBER OF BASE B DIGITS IN THE MANTISSA. C C IF ONE WISHES TO CHANGE THE PRECISION TO SINGLE PRECISION (FOR C EXAMPLE), THEN ONLY THREE THINGS NEED TO BE DONE: C C (1) REPLACE ALL THE STATEMENTS ONE=1.0D0 IN THE PACKAGE BY ONE=1.0, C C (2) PLACE A C IN COLUMN 1 OF ALL THE STATEMENTS IN THE PACKAGE C IMPLICIT REAL*8 (A-H,O-Z), C C (3) REPLACE ALL OCCURENCES OF D1MACH IN THE PACKAGE BY R1MACH C IF THE NETLIB SUPPLIED SUBPROGRAMS I1MACH, D1MACH, AND R1MACH ARE C BEING USED, AND OTHERWISE MERELY CHANGE THE DEFINITION OF C D1MACH(3) IN SUBPROGRAM D1MACH TO REFLECT SINGLE PRESISION. C C (WE NOTE HERE THAT THE ONLY ACTIVE WRITE STATEMENTS IN THIS PACKAGE C ARE IN THE SAMPLE DRIVER PROGRAM, BUT SOME OF THOSE WHICH HAVE BEEN C COMMENTED OUT (ALONG WITH THE STATEMENTS NWRIT=I1MACH(2)) ELSEWHERE C IN THE PROGRAM COULD PROVE USEFUL, ESPECIALLY THE STATEMENT C WRITE(NWRIT,1100)... IN SUBROUTINE CONMAX.) C C (B) CREATE A DRIVER (I.E. MAIN) PROGRAM WHICH DIMENSIONS THE ARRAYS C FUN, PTTBL, IWORK, WORK, PARAM, AND ERROR, SETS THE INPUT VARIABLES C IOPTN, NPARM, NUMGR, ITLIM, IFUN, IPTB, INDM, LIWRK, LWRK, PARAM C (AND POSSIBLY ALSO FUN AND PTTBL), CONTAINS THE STATEMENT C C CALL CONMAX(IOPTN,NPARM,NUMGR,ITLIM,FUN,IFUN,PTTBL,IPTB, C *INDM,IWORK,LIWRK,WORK,LWRK,ITER,PARAM,ERROR) C C AND PRINTS THE OUTPUT. IN ADDITION, WITH CERTAIN SETTINGS OF IOPTN C (SEE BELOW) THE USER WILL ALSO SUPPLY INPUT VALUES IN IWORK(1) AND C WORK(1), AND/OR IWORK(2), AND/OR WORK(2); SINCE IWORK AND WORK ARE C CHANGED BY THE PROGRAM, IN THIS CASE THE USER WILL NEED TO RESET THESE C VALUES EACH TIME CONMAX IS CALLED. C C (C) CREATE A SUBROUTINE FNSET (DESCRIBED LATER) FOR INPUT OF FUNCTION C DATA AND THE VECTOR ARRAY ICNTYP, WHICH SPECIFIES THE TYPE OF THE C CONSTRAINTS. C C THE VARIABLES IN THE CALLING SEQUENCE FOR CONMAX ARE THE FOLLOWING. C C IOPTN (INPUT) THIS IS THE OPTION SWITCH, WHICH SHOULD BE SET TO C 0 UNLESS ONE OR MORE OF THE EXTRA OPTIONS DESCRIBED BELOW C IS USED. C C NPARM (INPUT) THIS IS THE NUMBER OF PARAMETERS IN THE PROBLEM. C (THEY ARE STORED IN PARAM--SEE BELOW.) C C NUMGR (INPUT) THIS IS THE TOTAL NUMBER OF CONSTRAINTS. C C ITLIM (INPUT) THIS IS THE LIMIT ON THE NUMBER OF ITERATIONS, I.E. C THE LIMIT ON THE NUMBER OF TIMES THE PROGRAM REDUCES W. IF C ITLIM IS SET TO 0 THE PROGRAM WILL COMPUTE THE ERRORS FOR C THE INITIAL APPROXIMATION AND STOP WITHOUT CHECKING C FEASIBILITY. C C FUN (OPTIONAL INPUT) (VECTOR ARRAY OF DIMENSION IFUN) THIS IS C A VECTOR ARRAY IN WHICH DATA OR FUNCTION VALUES IN TYPE 2 C CONSTRAINTS (SEE ABOVE) CAN BE STORED. FUN(I) NEED NOT BE C ASSIGNED A VALUE IF IT IS NOT GOING TO BE USED. C C IFUN (INPUT) THIS IS THE DIMENSION OF FUN IN THE DRIVER PROGRAM. C IT MUST BE .GE. THE LARGEST INDEX I FOR WHICH FUN(I) IS USED C UNLESS NO FUN(I) IS USED, IN WHICH CASE IT MUST BE .GE. 1. C C PTTBL (OPTIONAL INPUT) (MATRIX ARRAY OF DIMENSION (IPTB,INDM)) C ROW I OF THIS ARRAY NORMALLY CONTAINS A POINT USED IN THE ITH C CONSTRAINT. THE ENTRIES IN ROW I NEED NOT BE ASSIGNED VALUES IF C SUCH A POINT IS NOT USED IN THE ITH CONSTRAINT. C (EXAMPLE: IF THE LEFT SIDE OF CONSTRAINT I IS A POLYNOMIAL IN C ONE INDEPENDENT VARIABLE, THEN THE VALUE OF THE INDEPENDENT C VARIABLE SHOULD BE IN PTTBL(I,1), AND THE COEFFICIENTS SHOULD BE C IN PARAM.) C PTTBL CAN ALSO BE USED TO PASS OTHER INFORMATION FROM THE DRIVER C PROGRAM TO SUBROUTINE FNSET. C C IPTB (INPUT) THIS IS THE FIRST DIMENSION OF PTTBL IN THE DRIVER C PROGRAM. IT MUST BE .GE. THE LARGEST SUBSCRIPT I FOR WHICH A C VALUE PTTBL(I,J) IS USED, AND MUST BE .GE. 1 IF NO SUCH VALUES C ARE USED. C C INDB (INPUT) THIS IS THE SECOND DIMENSION OF PTTBL IN THE DRIVER C PROGRAM. IT MUST BE .GE. THE LARGEST SUBSCRIPT J FOR WHICH A C VALUE PTTBL(I,J) IS USED, AND MUST BE .GE. 1 IF NO SUCH VALUES C ARE USED. C C IWORK (INPUT) (VECTOR ARRAY OF DIMENSION LIWRK) THIS IS AN INTEGER C WORK ARRAY. THE USER NEED NOT PLACE ANY VALUES IN IT, EXCEPT C POSSIBLY CERTAIN OPTIONAL INFORMATION AS DESCRIBED BELOW. C C LIWRK (INPUT) THIS IS THE DIMENSION OF IWORK IN THE DRIVER PROGRAM. C IT MUST BE AT LEAST 7*NUMGR + 7*NPARM + 3. IF NOT, CONMAX WILL C RETURN WITH THIS MINIMUM VALUE MULTIPLIED BY -1 AS A WARNING. C C WORK (INPUT) (VECTOR ARRAY OF DIMENSION LWRK) THIS IS A FLOATING C POINT WORK ARRAY. THE USER NEED NOT PLACE ANY VALUES IN IT, C EXCEPT POSSIBLY CERTAIN OPTIONAL INFORMATION AS DESCRIBED BELOW. C C LWRK (INPUT) THIS IS THE DIMENSION OF WORK IN THE DRIVER PROGRAM. C IT MUST BE AT LEAST 2*NPARM**2 + 4*NUMGR*NPARM + 11*NUMGR + C 27*NPARM + 13. IF NOT, CONMAX WILL RETURN WITH THIS MINIMUM C VALUE MULTIPLIED BY -1 AS A WARNING. C C ITER (OUTPUT) THIS IS THE NUMBER OF ITERATIONS PERFORMED BY CONMAX, C INCLUDING THOSE USED IN ATTEMPTING TO GAIN FEASIBILITY, C UNTIL EITHER IT CAN NO LONGER IMPROVE THE SITUATION OR THE C ITERATION LIMIT IS REACHED. IF ITER=ITLIM IT IS POSSIBLE C THAT THE PROGRAM COULD FURTHER REDUCE W IF RESTARTED C (POSSIBLY WITH THE NEW PARAMETERS). ITER=-1 IS A SIGNAL C THAT TYPE -1 FEASIBILITY COULD NOT BE ACHIEVED, IN THIS CASE C ERROR WILL CONTAIN THE VALUES COMPUTED USING THE USER C SUPPLIED INITIAL APPROXIMATION. ITER=-2 IS A SIGNAL THAT C TYPE -1 FEASIBILITY WAS ACHIEVED BUT TYPE -2 FEASIBILITY C COULD NOT BE ACHIEVED, IN THIS CASE VALUES IN ERROR C CORRESPONDING TO TYPE 1 OR TYPE 2 CONSTRAINTS WILL BE ZERO. C C PARAM (INPUT AND OUTPUT) (VECTOR ARRAY OF DIMENSION AT LEAST NPARM C IN THE DRIVER PROGRAM) THE USER SHOULD PLACE AN INITIAL GUESS C FOR THE PARAMETERS IN PARAM, AND ON OUTPUT PARAM WILL CONTAIN C THE BEST PARAMETERS CONMAX HAS BEEN ABLE TO FIND. IF THE C INITIAL PARAM IS NOT FEASIBLE THE PROGRAM WILL FIRST TRY TO C FIND A FEASIBLE PARAM. C C ERROR (OUTPUT) (VECTOR ARRAY OF DIMENSION AT LEAST NUMGR + 3 IN THE C DRIVER PROGRAM) FOR I=1,...,NUMGR, CONMAX WILL PLACE IN C ERROR(I) THE ERROR IN CONSTRAINT I (DEFINED TO BE THE VALUE C OF THE LEFT SIDE OF CONSTRAINT I, EXCEPT WITHOUT THE ABSOLUTE C VALUE IN TYPE 2 CONSTRAINTS). FURTHER, C C ERROR(NUMGR+1) WILL BE THE (PRINCIPAL) ERROR NORM, THAT IS, THE C MAXIMUM VALUE OF THE LEFT SIDES OF THE TYPE 2 (INCLUDING THE C ABSOLUTE VALUE) AND TYPE 1 CONSTRAINTS. C C ERROR(NUMGR+2) WILL BE THE MAXIMUM VALUE OF THE LEFT SIDES OF THE C TYPE -1 CONSTRAINTS, OR 0.0 IF THERE ARE NO TYPE -1 C CONSTRAINTS. EXCEPT FOR ROUNDOFF ERROR AND SMALL TOLERANCES C IN SOME SUBROUTINES THIS VALUE WILL NORMALLY BE .LE. 0.0, AND C IT WILL NOT BE ALLOWED TO BE .GT. TOLCON IN THE MAIN PART OF C THE PROGRAM. C C ERROR(NUMGR+3) WILL BE THE MAXIMUM VALUE OF THE LEFT SIDES OF THE C TYPE -2 CONSTRAINTS, OR 0.0 IF THERE ARE NO TYPE -2 C CONSTRAINTS. THIS VALUE SHOULD BE .LE. TOLCON, SINCE THE C PROGRAM WILL NOT EVEN ATTEMPT TO COMPUTE VALUES FOR THE C TYPE 2 AND TYPE 1 CONSTRAINTS OTHERWISE (EXCEPT FOR VALUES C CORRESPONDING TO THE INITIAL PARAMETERS PLACED IN PARAM BY C THE USER). THE USER CAN USE THIS FEATURE TO INSERT TYPE -2 C OR -1 CONSTRAINTS TO KEEP THE PARAMETERS AWAY FROM VALUES C WHERE A TYPE 2 OR TYPE 1 CONSTRAINT IS UNDEFINED. C C C THE SUBPROGRAM FNSET CREATED BY THE USER MUST HAVE AS ITS FIRST THREE C STATEMENTS C C SUBROUTINE FNSET(NPARM,NUMGR,PTTBL,IPTB,INDM,PARAM,IPT, C *INDFN,ICNTYP,CONFUN) C C IMPLICIT REAL*8 (A-H,O-Z) C C DIMENSION PTTBL(IPTB,INDM),PARAM(NPARM),ICNTYP(NUMGR), C *CONFUN(NUMGR,NPARM+1) C C THE FIRST EIGHT VARIABLES IN THE CALLING SEQUENCE FOR FNSET ARE FOR C INPUT TO FNSET, WITH THE FIRST FIVE VARIABLES BEING EXACTLY AS THE C USER SET THEM IN THE DRIVER PROGRAM. IF THE TEN THOUSANDS DIGIT OF C IOPTN WAS SET TO 0, FNSET SHOULD BE WRITTEN TO PLACE THE APPROPRIATE C VALUES IN ICNTYP AND CONFUN USING THE PARAMETERS IN PARAM, AS FOLLOWS. C C ICNTYP(IPT) = THE TYPE OF THE IPT(TH) CONSTRAINT (I.E. 2, 1, -1, C OR -2), OR THE USER CAN SET ICNTYP(IPT)=0 AS A SIGNAL TO IGNORE C CONSTRAINT IPT. C C CONFUN(IPT,1) = THE APPROPRIATE VALUE AS DISCUSSED ABOVE. (THIS CAN C BE LEFT UNDEFINED IF ICNTYP(IPT)=0.) C C IF INDFN=1 (WHICH IS THE ONLY POSSIBILITY OTHER THAN INDFN=0) THEN IN C ADDITION TO THE ABOVE, FOR J=1,...,NPARM, FNSET SHOULD COMPUTE C C CONFUN(IPT,J+1) = THE VALUE OF THE PARTIAL DERIVATIVE WITH RESPECT C TO PARAM(J) OF THE FUNCTION WHOSE VALUE WAS COMPUTED IN C CONFUN(IPT,1) (UNLESS ICNTYP(IPT)=0, IN WHICH CASE THESE VALUES C NEED NOT BE COMPUTED). C C C THE USER HAS SEVERAL EXTRA OPTIONS WHICH ARE ACTIVATED BY SETTING C IOPTN TO A VALUE OTHER THAN 0; MORE THAN ONE AT A TIME CAN BE USED. C IN PARTICULAR, C C IF THE ONES DIGIT OF IOPTN IS 0, THEN THE USER SHOULD GIVE FORMULAS C IN FNSET FOR COMPUTING THE PARTIAL DERIVATIVES WHEN INDFN=1 AS DESCRIBED C ABOVE. C C IF THE USER SETS THE ONES DIGIT OF IOPTN TO 1, THEN INDFN WILL ALWAYS C BE 0 WHEN FNSET IS CALLED, AND THE PROGRAM WILL AUTOMATICALLY C APPROXIMATE THE PARTIAL DERIVATIVES WHEN REQUIRED USING THE FOLLOWING C CENTERED DIFFERENCE FORMULA: C PARTIAL DERIVATIVE WITH RESPECT TO THE JTH VARIABLE (WHERE 1 .LE. J C .LE. NPARM) OF THE FUNCTION WHOSE VALUE IS COMPUTED IN CONFUN(IPT,1) C (WHERE 1 .LE. IPT .LE. NUMGR) IS APPROXIMATELY THE VALUE OF THIS C FUNCTION WHEN THE JTH COMPONENT OF PARAM IS INCREASED BY DELT, MINUS C THE VALUE OF THIS FUNCTION WHEN THE JTH COMPONENT OF PARAM IS C DECREASED BY DELT, ALL DIVIDED BY 2.0*DELT, WHERE DELT = SQRT(B**(-ITT)), C WHERE B IS THE BASE FOR FOR FLOATING POINT NUMBERS AND ITT IS THE NUMBER C OF BASE B DIGITS IN THE MANTISSA. C C IF THE TENS DIGIT OF IOPTN IS 0, THE PROGRAM WILL NOT GIVE UP C UNTIL EITHER AN ITERATION FAILS TO PRODUCE A REDUCTION ABS(ENCHG) OF C MORE THAN 100.0*B**(-ITT) IN THE PRINCIPAL ERROR NORM, OR ITLIM C ITERATIONS HAVE BEEN USED. C C IF THE TENS DIGIT OF IOPTN IS 1, THE PROGRAM WILL ALSO GIVE UP IF C ABS(ENCHG) .LT. ENCSM FOR LIMSM CONSECUTIVE STEPS IN THE MAIN PART OF C THE PROGRAM (I.E. AFTER TYPE -1 AND TYPE -2 FEASIBILITY HAVE BOTH BEEN C ACHIEVED). THIS OPTION MAY BE USEFUL IN SAVING SOME TIME BY C FORESTALLING A LONG STRING OF ITERATIONS AT THE END OF A RUN WITH ONLY C A TINY IMPROVEMENT IN EACH ONE. ENCSM AND LIMSM ARE TRANSMITTED TO C CONMAX IN WORK(1) AND IWORK(1) RESPECTIVELY. WORK(1) AND IWORK(1) DO C NOT NEED TO BE ASSIGNED VALUES IF THE TENS DIGIT OF IOPTN IS 0. C C IF THE HUNDREDS DIGIT OF IOPTN IS 0 OR 2, THEN THE INTERNAL PARAMETER C NSTEP WILL BE GIVEN THE DEFAULT VALUE 1. NSTEP IS THE NUMBER OF C RUNGE-KUTTA STEPS USED IN EACH RK ITERATION. C C IF THE HUNDREDS DIGIT OF IOPTN IS 1 OR 3, THEN THE VALUE OF NSTEP USED C WILL BE THE POSITIVE INTEGER VALUE PLACED IN IWORK(2) BY THE USER IN THE C DRIVER PROGRAM. SETTING NSTEP LARGER THAN 1 MAY ALLOW THE PROGRAM TO C SUCCEED ON DIFFICULT PROBLEMS WHERE THE CONVERGENCE WOULD BE EXTREMELY C SLOW WITH NSTEP=1, BUT IT WILL GENERALLY CAUSE MORE COMPUTER TIME TO BE C USED IN EACH RK ITERATION. IWORK(2) DOES NOT NEED TO BE ASSIGNED A C VALUE IF THE HUNDREDS DIGIT OF IOPTN IS 0 OR 2. (NSTEP IS SOMETIMES C CALLED THE OVERDRIVE PARAMETER.) C C IF THE HUNDREDS DIGIT OF IOPTN IS 0 OR 1, THEN THE INTERNAL PARAMETER C TOLCON WILL BE GIVEN THE DEFAULT VALUE SQRT(B**(-ITT)), WHERE B IS THE C BASE FOR FLOATING POINT NUMBERS AND ITT IS THE NUMBER OF BASE B DIGITS C IN THE MANTISSA. C C IF THE HUNDREDS DIGIT OF IOPTN IS 2 OR 3, THEN THE VALUE OF TOLCON USED C WILL BE THE VALUE PLACED IN WORK(2) BY THE USER IN THE DRIVER PROGRAM. C THIS VALUE SHOULD ALWAYS BE NONNEGATIVE. IF THERE ARE NO TYPE -2 OR -1 C CONSTRAINTS THEN THE SETTING OF TOLCON WILL HAVE NO EFFECT, BUT IF C THERE ARE TYPE -2 OR -1 CONSTRAINTS THEN IN GENERAL SMALLER VALUES OF C TOLCON MAY GIVE MORE ACCURACY IN THE FINAL ANSWER, BUT MAY SLOW DOWN C OR PREVENT CONVERGENCE, WHILE LARGER VALUES OF TOLCON MAY HAVE THE C REVERSE EFFECT. IN PARTICULAR, IF THE TYPE -2 AND -1 CONSTRAINTS C CANNOT BE SATISFIED SUMULTANEOUSLY WITH STRICT INEQUALITY (THIS CASE C COULD OCCUR, FOR EXAMPLE, IF AN EQUALITY CONSTRAINT G = 0.0 WAS C ENTERED AS THE TWO INEQUALITY CONSTRAINTS G .LE. 0.0 AND C -G .LE. 0.0), THEN SETTING TOLCON=0.0 WILL ALMOST CERTAINLY CAUSE THE C PROGRAM TO FAIL, SINCE AT EACH ITERATION THE PROGRAM WILL NOT ACCEPT C PROSPECTIVE NEW VALUES OF THE PARAMETERS UNLESS IT CAN CORRECT THEM C BACK INTO THE RELAXED FEASIBLE REGION WHERE CONFUN(IPT,1) .LE. TOLCON C FOR ALL THE TYPE -2 AND -1 CONSTRAINTS. C C IF THE THOUSANDS DIGIT OF IOPTN IS 0, THE PROGRAM WILL USE THE RK METHOD C (WHICH INVOLVES APPLYING FOURTH ORDER RUNGE-KUTTA TO A SYSTEM OF C DIFFERENTIAL EQUATIONS), THEN IF THIS FAILS IT WILL TRY TO REDUCE C W WITH AN SLP STEP (I.E. SUCCESSIVE LINEAR PROGRAMMING WITH A TRUST C REGION), THEN GO BACK TO RK IF THE SLP STEP REDUCES W. C C IF THE THOUSANDS DIGIT OF IOPTN IS 1, THE PROGRAM WILL USE SLP STEPS ONLY. C IN GENERAL, IN SOME PROBLEMS SLP WORKS VERY WELL, AND IN THOSE C PROBLEMS IT WILL USUALLY BE FASTER THAN RK BECAUSE IT REQUIRES FEWER C GRADIENT EVALUATIONS THAN RK, BUT IN OTHER PROBLEMS THE CONVERGENCE C OF SLP MAY BE EXCRUCIATINGLY SLOW, AND THE MORE POWERFUL RK METHOD C MAY BE MUCH FASTER. C C IF THE THOUSANDS DIGIT OF IOPTN IS 2 THE PROGRAM WILL USE THE RK METHOD C ONLY, QUITTING WHEN RK CAN NO LONGER PRODUCE AN IMPROVEMENT. THIS C MAY GIVE A LITTLE LESS ACCURACY THAN SETTING THE THOUSANDS DIGIT TO 0, C BUT MAY SAVE SIGNIFICANT COMPUTER TIME IN SOME CASES. C C IF THE TEN THOUSANDS DIGIT OF IOPTN IS 0, THEN FNSET SHOULD BE WRITTEN AS C DESCRIBED ABOVE. C C IF THE USER SETS THE TEN THOUSANDS DIGIT OF IOPTN TO 1, THEN FNSET SHOULD BE C WRITTEN AS DESCRIBED ABOVE EXCEPT THAT THE COMPUTATIONS SHOULD BE DONE C FOR ALL IPT=1,..,NUMGR INSTEAD OF FOR A SINGLE GIVEN VALUE OF IPT. C THIS OPTION MAY SAVE COMPUTER TIME IN PROBLEMS WHERE A LARGE PART OF C THE COMPUTATION IS THE SAME FOR DIFFERENT VALUES OF IPT, SINCE IT C AVOIDS UNNECESSARY REPITITION OF THIS COMMON COMPUTATION BY HAVING C THE LOOP OVER IPT IN FNSET INSTEAD OF OUTSIDE FNSET. C IF THE TEN THOUSANDS DIGIT OF IOPTN IS 1, EVEN MORE TIME CAN OFTEN BE C SAVED IF FNSET IS WRITTEN SO THAT ALL CONSTRAINTS ARE COMPUTED IF IPT=0, C BUT ONLY THE STANDARD (I.E. TYPE -1 OR -2) CONSTRAINTS ARE COMPUTED IF C IPT=-1. NOTE THAT IF THE TEN THOUSANDS DIGIT OF IOPTN IS 0, THEN IPT C WILL BE POSITIVE WHENEVER FNSET IS CALLED, INDICATING THAT ONLY C CONSTRAINT IPT SHOULD BE COMPUTED. C THE DRAWBACK OF USING THIS OPTION IS THAT IN GENERAL SOME CONSTRAINT C VALUES AND DERIVATIVES WILL BE COMPUTED UNNECESSARILY, WHICH COULD COST C TIME IN SOME PROBLEMS. C C IF THE COMPILER USED FOR THIS PROGRAM WILL NOT ACCEPT VARIABLE C DIMENSIONING AS USED HERE, ONE CAN (1) CHANGE ALL THE DIMENSIONS TO C THE APPROPRIATE INTEGER CONSTANTS, WHERE IN THE DIMENSION STATEMENTS C NPARM AND NUMGR CAN BE CONSISTENTLY REPLACED BY FIXED INTEGERS OF C EQUAL OR GREATER VALUE, (2) FOR EACH SUBPROGRAM EXCEPT FNSET, DELETE C ALL ARRAY NAMES FROM THE ARGUMENT LISTS IN THE SUBROUTINE OR FUNCTION C STATEMENT EXCEPT FUN, PTTBL, PARAM, ERROR, IWORK, AND WORK, AND DELETE C THE CORRESPONDING ELEMENTS IN ALL CALL STATEMENTS AND FUNCTION C SUBPROGRAM REFERENCES, AND (3) USE AN EQUIVALENCE STATEMENT AFTER C THE DIMENSION STATEMENT IN EACH SUBPROGRAM WITH ARRAYS DELETED FROM C THE ARGUMENT LIST TO ASSOCIATE THE FIRST ELEMENT OF EACH DELETED C ARRAY WITH THE APPROPRIATE ELEMENT IN IWORK OR WORK USING THE C INFORMATION IN FUNCTION SUBPROGRAM ILOC. C C FOR FURTHER INFORMATION ABOUT THE ALGORITHM AND PROGRAM, SEE THE PAPER C C AN ODE-BASED APPROACH TO NONLINEARLY CONSTRAINED MINIMAX PROBLEMS, BY C E. H. KAUFMAN, JR., D. J. LEEMING, AND G. D. TAYLOR, NUMERICAL C ALGORITHMS (1995), PP. 25-37. C C NOTE THAT SINCE THAT PAPER WAS WRITTEN, TOLCON WAS DELETED FROM THE C ARGUMENT LIST OF CONMAX. C C C A SAMPLE DRIVER PROGRAM AND SUBROUTINE FNSET ARE INCLUDED IN THIS C PACKAGE, AND WE NOW DISCUSS THE EXAMPLE USED THEREIN; FOR MORE C INFORMATION, SEE THE COMMENTS IN THESE TWO ROUTINES. C C THE EXAMPLE IS TO CHOOSE (DOUBLE PRECISION) PARAMETERS A, B, C, AND D C TO MINIMIZE C C MAX{ MAX (|F(X,Y) - (AX+BY+C)/(DX+1)|, |(AX+BY+C)/(DX+1)|) : C (X,Y) IN Z} C C SUBJECT TO THE CONSTRAINTS C C DX+1 .GE. EPS FOR (X,Y) IN Z, C C AND THE FIRST PARTIAL DERIVATIVE OF (AX+BY+C)/(DX+1) WITH RESPECT TO C X IS .LE. 0.0 AT (X,Y) = (0.0,0.0). C C HERE WE ARE TAKING Z = {(0.0,0.0), (0.0,1.0), (-2.0/3.0,1.0/3.0), C (1.0,1.0), (1.0,2.0)}, EPS = .0001, F(0.0,0.0) = .5, F(0.0,1.0) = 1.0, C F(-2.0/3.0,1.0/3.0) = -1.0, F(1.0,1.0) =1.5, F(1.0,2.0) =-1.0, C AND TAKING THE INITIAL GUESSES FOR THE PARAMETERS TO BE C A = B = C = D = 0.0 C C TO USE CONMAX, WE WRITE THIS PROBLEM AS THE OPTIMIZATION PROBLEM C C MINIMIZE W, SUBJECT TO C C |F(X,Y) - (AX+BY+C)/(DX+1)| .LE. W, (X,Y) IN Z, (TYPE 2) C C (AX+BY+C)/(DX+1) .LE. W, (X,Y) IN Z, (TYPE 1) C C -(AX+BY+C)/(DX+1) .LE. W, (X,Y) IN Z, (TYPE 1) C C -DX - 1 + EPS .LE. 0, (X,Y) IN Z, (TYPE -1) C C A - CD .LE. 0. (TYPE -2) C C ONE CAN PROVE THAT THE UNIQUE BEST VALUES ARE A = -B = C = D = 1.0, C WITH W (= ENORM = ERROR(NUMGR+1)) = 1.0. ONE CAN ALSO PROVE THAT C THERE ARE LOCAL SOLUTIONS WHICH ARE NOT GLOBAL SOLUTIONS, THAT IS, C SOLUTIONS FOR WHICH W CANNOT BE REDUCED BY SMALL CHANGES IN A, B, C, C AND D, BUT FOR WHICH A, B, C, D CAN BE FOUND SATISFYING THE CONSTRAINTS C AND GIVING SMALLER W. SOME SUCH SOLUTIONS ARE GIVEN BY A = B = D = 0.0, C C = 0.25, WITH W = 1.25, AND OTHER CHOICES WHERE THE RATIONAL FUNCTION C REDUCES TO THE CONSTANT 0.25 AND THE COEFFICIENTS SATISFY THE CONSTRAINTS. C WHEN THE PROGRAM IS RUN, ANY OF THESE SOLUTIONS MAY BE FOUND (UP TO A C SMALL DISCREPANCY DUE IN PART TO ROUNDOFF), DEPENDING ON THE ACCURACY OF C THE COMPUTER BEING USED. C C THE OUTPUT FOR THE SAMPLE DRIVER AND FNSET FOLLOWS. THIS WAS RUN ON C A MAC SE (WITH ONE MEGABYTE OF RAM) WITH D1MACH(3) SET TO 16.0D0**(-14) C USING THE ABSOFT MACFORTRAN COMPILER (VERSION 2.4); RUNS WITH A C DIFFERENT MACHINE AND/OR A DIFFERENT D1MACH(3) AND/OR A DIFFERENT C COMPILER COULD PRODUCE DIFFERENT RESULTS, ESPECIALLY CONSIDERING THE C POSSIBILITY OF LOCAL SOLUTIONS WHICH ARE NOT GLOBALY BEST. C C IOPTN IS 0 NPARM IS 4 NUMGR IS 21 C C ITLIM IS 100 IFUN IS 5 IPTB IS 6 INDM IS 2 C C THE FUNCTION VALUES ARE C 0.50000E+00 0.10000E+01 -0.10000E+01 0.15000E+01 C -0.10000E+01 C C THE POINTS ARE C 0.00000E+00 0.00000E+00 C 0.00000E+00 0.10000E+01 C -0.66667E+00 0.33333E+00 C 0.10000E+01 0.10000E+01 C 0.10000E+01 0.20000E+01 C C EPS IS 0.10000E-03 C C THE INITIAL PARAMETERS ARE C C 0.00000000000000E+00 0.00000000000000E+00 0.00000000000000E+00 C C 0.00000000000000E+00 C C *****AFTER CONMAX ITER IS 13 LIWRK IS 178 LWRK IS 720 C C THE FINAL PARAMETERS ARE C C -0.90123453614326E-02 -0.12378422071992E-12 0.25000000000006E+00 C C -0.36049381446237E-01 C C THE ERROR NORMS ARE C C 0.1250000000000E+01 -0.9638506185538E+00 0.1288396472843E-12 C C THE ERRORS ARE C C 0.24999999999994E+00 0.75000000000006E+00 -0.12499999999999E+01 C 0.12499999999999E+01 -0.12499999999999E+01 0.25000000000006E+00 C -0.25000000000006E+00 0.24999999999994E+00 -0.24999999999994E+00 C 0.24999999999994E+00 -0.24999999999994E+00 0.25000000000006E+00 C -0.25000000000006E+00 0.24999999999994E+00 -0.24999999999994E+00 C -0.99990000000000E+00 -0.99990000000000E+00 -0.10239329209642E+01 C -0.96385061855376E+00 -0.96385061855376E+00 0.12883964728427E-12 C C C WE NOW DESCRIBE FOUR SUBUNITS OF CONMAX WHICH THE USER CAN CALL C DIRECTLY IF IT IS DESIRED TO ACOMPLISH CERTAIN TASKS MORE SIMPLY C AND QUICKLY THAN SETTING UP THE PROBLEMS FOR CONMAX. IN EACH CASE A C DRIVER PROGRAM (AND A SUBROUTINE FNSET IF NECESSARY) IS SUPPLIED, C WHICH CAN BE USED BY REPLACING THE C IN COLUMN 1 OF EACH LINE BY A C BLANK (EXCEPT CC IS REPLACED BY C) AND APPROPRIATELY MODIFYING THOSE C STATEMENTS WHICH ARE INDICATED TO BE USER SETTABLE. THE OPEN STATEMENT C MAY ALSO NEED TO BE CHANGED, DEPENDING ON THE MACHINE. OUTPUT PRODUCED C BY THE SAMPLE DRIVER ON THE MAC SE IS ALSO PROVIDED. ALTHOUGH CERTAIN C ADVANTAGES (SUCH AS REDUCED STORAGE) COULD BE GAINED BY REWRITING SOME C OF THE SUBPROGRAMS IN THE CONMAX PACKAGE, FOR SIMPLICITY THE DRIVER C PROGRAMS ARE SET UP SO THAT THE SUBPROGRAMS IN THIS PACKAGE CAN BE USED C AS IS, WITH NO CHANGES OTHER THAN IN THE DRIVER PROGRAMS AND (IN CASES C (A) AND (B) BELOW) THE SUBROUTINES FNSET. THUS IF DESIRED ONE COULD C COMPILE AND LOAD THE ENTIRE CONMAX PACKAGE (EXCEPT FOR DRIVERS AND C SUBROUTINES FNSET) AND USE IT AS A LIBRARY. FOR SIMPLICITY WE ALSO DO C NOT DESCRIBE HERE CERTAIN ADDITIONAL OPTIONS, SUCH AS A HOT START C OPTION IN (C) AND (D), BUT MORE INFORMATION CAN BE OBTAINED FROM THE C COMMENTS IN THE SUBPROGRAMS INVOLVED IF DESIRED. C C C (A) MULLERS METHOD DERIVATIVE FREE REAL ROOT FINDING C C (SUBPROGRAMS INVOLVED: MULLER, FNSET (USER SUPPLIED), ILOC, D1MACH, C ERCMP1) C C GIVEN A FUNCTION F OF ONE VARIABLE (WHERE F(X) IS COMPUTED IN SUBROUTINE C FNSET AS CONFUN(1,1), WITH X = PARAM(1)), A NONNEGATIVE TOLERANCE TOLCON, C TWO POINTS (P1,F1) AND (PROCOR,EMIN) WITH F1 = F(P1), EMIN = F(PROCOR), C P1 .LT. PROCOR, F1 .GT. TOLCON, AND EMIN .LT. -TOLCON, THE PROGRAM WILL C ATTEMPT TO LOCATE A NEW PROCOR WITH NEW EMIN = F(PROCOR) AND ABS(EMIN) C .LE. TOLCON. IF IT FAILS TO DO THIS, THE PROGRAM WILL RETURN WITH C PROCOR = THE LEFTMOST ABSCISSA FOUND WITH EMIN = F(PROCOR) .LT. -TOLCON. C NOTE: IF INSTEAD OF F1 .GT. TOLCON AND EMIN .LT. -TOLCON WE START WITH C F1 .LT. -TOLCON AND EMIN .GT. TOLCON, THIS PROGRAM CAN STILL BE USED C BY REPLACING F BY -F BEFORE RUNNING THE PROGRAM. C THE SOLUTION PROCEDURE IS A MODIFICATION OF THE FOLLOWING: BISECT THE C INTERVAL [P1,PROCOR] TO GET A THIRD POINT, PASS A QUADRATIC POLYNOMIAL C THROUGH THE THREE POINTS AND USE ITS UNIQUE ZERO IN [P1,PROCOR] TO C REPLACE P1 OR PROCOR, MAINTAINING THE CONDITIONS F(LEFT ENDPOINT) .GT. C TOLCON AND F(RIGHT ENDPOINT) .LT. -TOLCON, AND CONTINUE UNTIL A SOLUTION C IS FOUND OR F HAS BEEN COMPUTED 5 TIMES OR THE INTERVAL LENGTH FALLS C BELOW 100.0*B**(-ITT). THIS PROCEDURE MAY BE ESPECIALLY USEFUL IN CASES C WHERE F IS EXPENSIVE TO COMPUTE SINCE IT MAINTAINS A SHRINKING INTERVAL C ABOUT THE SOLUTION, HAS A HIGHER ORDER OF CONVERGENCE THAN THE REGULA C FALSI METHOD, AND REQUIRES NO DERIVATIVES. THE FOLLOWING SAMPLE DRIVER C PROGRAM AND SUBROUTINE FNSET ARE SET UP TO FIND PROCOR IN [-4.0D0,2.0D0] C WITH ABS(F(PROCOR)) .LE. 0.001D0, WHERE F(X) = 2.0D0**(-X) - 0.5D0 C (THE EXACT SOLUTION IS PROCOR = 1.0D0, EMIN = 0.0D0). C C SAMPLE DRIVER AND FNSET FOR (A) MULLERS METHOD C C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION DVEC(1),FUN(1),PTTBL(1,1),ZWORK(1),ERR1(4), C *PARWRK(1),IWORK(17),WORK(6) C OPEN(6,FILE='MULOUT') C DVEC(1)=1.0D0 C ZWORK(1)=0.0D0 C IWORK(16)=-2 CC*****BEGIN USER SETTABLE STATEMENTS 1 OF 2 C TOLCON=1.0D-3 C P1=-2.0D0 C F1=3.5D0 C PROCOR=2.0D0 C EMIN=-0.25D0 CC*****END USER SETTABLE STATEMENTS 1 OF 2 C WRITE(6,100)TOLCON,P1,F1,PROCOR,EMIN C 100 FORMAT(/10H TOLCON IS,E22.13//16H INITIALLY P1 IS,E22.13, C *7H F1 IS,E22.13//10H PROCOR IS,E22.13,9H EMIN IS,E22.13) C CALL MULLER(0,1,1,DVEC,FUN,1,PTTBL,1,1,ZWORK,TOLCON,0, C *IWORK,17,WORK,6,PARWRK,ERR1,P1,F1,PROCOR,EMIN) C WRITE(6,200)PROCOR,EMIN C 200 FORMAT(/23H AFTER MULLER PROCOR IS,E22.13//8H EMIN IS, C *E22.13) C STOP C END C SUBROUTINE FNSET(NPARM,NUMGR,PTTBL,IPTB,INDM,PARAM,IPT, C *INDFN,ICNTYP,CONFUN) C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION PTTBL(IPTB,INDM),PARAM(NPARM),ICNTYP(NUMGR), C *CONFUN(NUMGR,NPARM+1) CC*****BEGIN USER SETTABLE STATEMENTS 2 OF 2 C CONFUN(1,1)=2.0D0**(-PARAM(1))-0.5D0 CC*****END USER SETTABLE STATEMENTS 2 OF 2 C RETURN C END C C OUTPUT FOR (A) MULLERS METHOD C C TOLCON IS 0.1000000000000E-02 C C INITIALLY P1 IS -0.2000000000000E+01 F1 IS 0.3500000000000E+01 C C PROCOR IS 0.2000000000000E+01 EMIN IS -0.2500000000000E+00 C C AFTER MULLER PROCOR IS 0.9993746852789E+00 C C EMIN IS 0.2167645412207E-03 C C C (B) ONE-DIMENSIONAL DERIVATIVE-FREE QUADRATIC SEARCH FOR A POSITIVE C LOCAL MINIMUM C C (SUBPROGRAMS INVOLVED: SEARSL, FNSET (USER SUPPLIED), ILOC, D1MACH, C ERCMP1, RCHMOD, CORRCT, SEARCR, MULLER, WOLFE, CONENR, HOUSE, C DOTPRD, REFWL; ONLY THE FIRST FIVE OF THESE ARE ACTUALLY USED) C C GIVEN A FUNCTION F OF ONE VARIABLE (WHERE F(X) IS COMPUTED IN SUBROUTINE C FNSET AS CONFUN(1,1), WITH X = PARAM(1)) AND A POSITIVE NUMBER PROJCT, C THE PROGRAM WILL ATTEMPT TO LOCATE A NEW PROJCT WITH EMIN = F(PROJCT), C WITH THE NEW PROJCT APPROXIMATELY GIVING A LOCAL MINIMUM OF F IN C [(OLD PROJCT)/1024, 1024*(OLD PROJCT)]. IF IT FAILS TO DO THIS, THE C PROGRAM WILL RETURN WITH PROJCT = THE ABSCISSA FOUND WITH SMALLEST C EMIN = F(PROJCT). ON OUTPUT, NSRCH WILL BE THE NUMBER OF EVALUATIONS C OF F THAT WERE DONE. THE SOLUTION PROCEDURE IS A MODIFICATION OF THE C FOLLOWING: COMPUTE F(PROJCT/2), F(PROJCT), AND F(2*PROJCT), IF THE C CONDITIONS F(MIDDLE POINT) .LE. F(LEFT POINT) AND F(MIDDLE POINT) .LE. C F(RIGHT POINT) ARE NOT BOTH SATISFIED THEN TRY TO GET THIS BY COMPUTING C F AT SMALLER (OR LARGER) POINTS AT MOST 3 MORE TIMES; ONCE THE C CONDITIONS ARE SATISFIED, ASSUMING THE POINTS ARE NOT TOO CLOSE TO BEING C COLLINEAR, PASS A QUADRATIC POLYNOMIAL THROUGH THE THREE POINTS AND USE C ITS UNIQUE MINIMUM IN THE INTERVAL TO REPLACE ONE OF THE ENDPOINTS WHILE C MAINTAINING THE TWO CONDITIONS, CONTINUING UNTIL F HAS BEEN COMPUTED 4 C MORE TIMES OR THE INTERVAL LENGTH FALLS BELOW 100.0*B**(-ITT) OR THE C POINTS BECOME NEARLY COLLINEAR. THE FOLLOWING SAMPLE DRIVER AND C SUBROUTINE FNSET ARE SET UP TO APPROXIMATE A SOLUTION OF THE LINE SEARCH C PROBLEM OF MINIMIZING G((6.0D0,2.0D0) + PROJCT*(-2.0D0,-1.0D0)), WHERE C G(U,V) = 3.0D0*ABS(U) + 2.0D0*ABS(V). WE START WITH PROJCT = 1.0D0, C THEN RUN THE PROGRAM AGAIN STARTING WITH THE RESULT OF THE FIRST RUN. C (THE EXACT SOLUTION IS PROJCT = 3.0D0, EMIN = 2.0D0; CONVERGENCE IS C RATHER SLOW, MAINLY BECAUSE F IS NOT DIFFERENTIABLE AT THE MINIMUM.) C C SAMPLE DRIVER AND FNSET FOR (B) ONE-DIMENSIONAL SEARCH C C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION X(2),FUN(1),PTTBL(1,1),PARAM(1),ERROR(4),IACT(1), C *IWORK(17),WORK(42),ERR1(4),PARPRJ(1),PARSER(1) C OPEN(6,FILE='SEROUT') C SPCMN=D1MACH(3) C TOL1=100.0D0*SPCMN C TOLCON=SQRT(SPCMN) C IACT(1)=1 C IWORK(7)=1 C PRJLIM=1.0D0/SPCMN C PARAM(1)=0.0D0 C X(1)=1.0D0 CC*****BEGIN USER SETTABLE STATEMENTS 1 OF 2 C PROJCT=1.0D0 CC*****END USER SETTABLE STATEMENTS 1 OF 2 C WRITE(6,100)PROJCT C 100 FORMAT(/20H INITIALLY PROJCT IS,E22.13) C CALL SEARSL(0,1,1,PRJLIM,TOL1,X,FUN,1,PTTBL,1,1,PARAM,ERROR, C *2.0D0,1,IACT,0,1.0D0,TOLCON,2.0D0,0,0,IWORK,17,WORK,42, C *ERR1,PARPRJ,PROJCT,EMIN,EMIN1,PARSER,NSRCH) C WRITE(6,200)PROJCT,EMIN,NSRCH C 200 FORMAT(/23H AFTER SEARSL PROJCT IS,E22.13//8H EMIN IS, C *E22.13,10H NSRCH IS,I4) C STOP C END C SUBROUTINE FNSET(NPARM,NUMGR,PTTBL,IPTB,INDM,PARAM,IPT, C *INDFN,ICNTYP,CONFUN) C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION PTTBL(IPTB,INDM),PARAM(NPARM),ICNTYP(NUMGR), C *CONFUN(NUMGR,NPARM+1) CC*****BEGIN USER SETTABLE STATEMENTS 2 OF 2 C U=6.0D0+PARAM(1)*(-2.0D0) C V=2.0D0+PARAM(1)*(-1.0D0) C CONFUN(1,1)=3.0D0*ABS(U)+2.0D0*ABS(V) CC*****END USER SETTABLE STATEMENTS 2 OF 2 C RETURN C END C C OUTPUT 1 FOR (B) ONE-DIMENSIONAL SEARCH C C INITIALLY PROJCT IS 0.1000000000000E+01 C C AFTER SEARSL PROJCT IS 0.2956250000000E+01 C C EMIN IS 0.2175000000000E+01 NSRCH IS 8 C C OUTPUT 2 FOR (B) ONE-DIMENSIONAL SEARCH C C INITIALLY PROJCT IS 0.2956250000000E+01 C C AFTER SEARSL PROJCT IS 0.3010732751581E+01 C C EMIN IS 0.2085862012645E+01 NSRCH IS 7 C C C (C) FREE-VARIABLE INEQUALITY-CONSTRAINED LINEAR PROGRAMMING C C (SUBPROGRAMS INVOLVED: SLNPRO, D1MACH, SJELIM) C C GIVEN POSITIVE INTEGERS M AND N WITH M .GE. N, AND A MATRIX V, THIS C PROGRAM WILL ATTEMPT TO SOLVE THE LINEAR PROGRAMMING PROBLEM C MAXIMIZE -V(M+1,1)*X(1)-...-V(M+1,N)*X(N) C SUBJECT TO V(I,1)*X(1)+...+V(I,N)*X(N) .LE. V(I,N+1) FOR I=1,...,M. C ON OUTPUT, INDIC WILL BE AN ERROR FLAG WHOSE VALUE WILL BE 0 FOR A C NORMAL SOLUTION, NEGATIVE FOR A POSSIBLY INACCURATE SOLUTION, AND C POSITIVE FOR A PROBABLE FAILURE. THE METHOD IS AN ENHANCED VERSION OF C THAT IN C AVDEYEVA, L. I. AND ZUKHOVITSKIY, S. I., LINEAR AND CONVEX PROGRAMMING, C SAUNDERS, PHILADELPHIA, 1966. C THE FOLLOWING SAMPLE DRIVER IS SET UP TO MAXIMIZE -Y, SUBJECT TO C X + Y .GE. 2.0D0 (I.E. -X - Y .LE. -2.0D0) AND X - Y .LE. 4.0DO. C (THE EXACT SOLUTION IS X = 3.0D0, Y = -1.0D0, WITH MAXIMUM OBJECTIVE C FUNCTION VALUE 1.0D0.) C C SAMPLE DRIVER FOR (C) LINEAR PROGRAMMING C C IMPLICIT REAL*8 (A-H,O-Z) CC THE MINIMUM DIMENSIONS ARE V(NUMGR+2*NPARM+1,NPARM+2), X(NPARM+1), CC IYCCT(NPARM+1), Y(NUMGR+2*NPARM), IXRCT(NUMGR+2*NPARM), CC IYRCT(NUMGR+2*NPARM), WHERE NPARM .GE. N AND NUMGR .GE. M. THE FIRST CC DIMENSION OF V MUST BE EXACTLY NUMGR+2*NPARM+1. CC*****BEGIN USER SETTABLE STATEMENTS 1 OF 1 C DIMENSION V(7,4),X(3),IYCCT(3),Y(6),IXRCT(6),IYRCT(6) C N=2 C M=2 C V(1,1)=-1.0D0 C V(1,2)=-1.0D0 C V(1,3)=-2.0D0 C V(2,1)=1.0D0 C V(2,2)=-1.0D0 C V(2,3)=4.0D0 C V(3,1)=0.0D0 C V(3,2)=1.0D0 C NPARM=N C NUMGR=M CC*****END USER SETTABLE STATEMENTS 1 OF 1 C MP1=M+1 C NP1=N+1 C V(MP1,NP1)=0.0D0 C IYRCT(1)=-1 C OPEN(6,FILE='LPOUT') C WRITE(6,100)M,N C 100 FORMAT(/10H THERE ARE,I5,17H CONSTRAINTS AND,I5, C *11H VARIABLES//32H THE CONSTRAINT COEFFICIENTS AND, C *16H RIGHT SIDES ARE) C DO 300 I=1,M C WRITE(6,200)(V(I,J),J=1,NP1) C 200 FORMAT(/(3E22.13)) C 300 CONTINUE C WRITE(6,400)(V(MP1,J),J=1,N) C 400 FORMAT(/40H THE NEGATIVES OF THE OBJECTIVE FUNCTION, C *17H COEFFICIENTS ARE//(3E22.13)) C CALL SLNPRO(V,M,N,IYRCT,Y,IXRCT,IYCCT,NPARM,NUMGR,X,INDIC) C WRITE(6,500)INDIC,(X(I),I=1,N) C 500 FORMAT(/31H AFTER SLNPRO THE ERROR FLAG IS,I4, C *19H THE VARIABLES ARE//(3E22.13)) C WRITE(6,600)V(MP1,NP1) C 600 FORMAT(/32H THE OBJECTIVE FUNCTION VALUE IS,E22.13) C STOP C END C C OUTPUT FOR (C) LINEAR PROGRAMMING C C THERE ARE 2 CONSTRAINTS AND 2 VARIABLES C C THE CONSTRAINT COEFFICIENTS AND RIGHT SIDES ARE C C -0.1000000000000E+01 -0.1000000000000E+01 -0.2000000000000E+01 C C 0.1000000000000E+01 -0.1000000000000E+01 0.4000000000000E+01 C C THE NEGATIVES OF THE OBJECTIVE FUNCTION COEFFICIENTS ARE C C 0.0000000000000E+00 0.1000000000000E+01 C C AFTER SLNPRO THE ERROR FLAG IS 0 THE VARIABLES ARE C C 0.3000000000000E+01 -0.1000000000000E+01 C C THE OBJECTIVE FUNCTION VALUE IS 0.1000000000000E+01 C C C (D) LEAST-DISTANCE QUADRATIC PROGRAMMING C C (SUBPROGRAMS INVOLVED: WOLFE, D1MACH, ILOC, CONENR, HOUSE, DOTPRD, C REFWL) C C GIVEN POSITIVE INTEGERS M AND N, AND A MATRIX PMAT, THIS PROGRAM WILL C ATTEMPT TO LOCATE AN N-DIMENSIONAL POINT WPT IN THE POLYHEDRON C DETERMINED BY THE INEQUALITIES C PMAT(1,J)*WPT(1)+...+PMAT(N,J)*WPT(N)+PMAT(N+1,J) .LE. 0.0 FOR J=1,...,M C WHOSE DISTANCE WDIST FROM THE ORIGIN IS MINIMIZED. ON OUTPUT, JFLAG C WILL BE AN ERROR FLAG WHOSE VALUE WILL BE 0 FOR A NORMAL SOLUTION AND C POSITIVE FOR A LIKELY FAILURE. THE METHOD IS AN ENHANCED VERSION OF C THAT IN C WOLFE, PHILIP, FINDING THE NEAREST POINT IN A POLYTOPE, MATHEMATICAL C PROGRAMMING 11 (1976), 128-149. C THE FOLLOWING SAMPLE DRIVER IS SET UP TO FIND THE NEAREST POINT TO C THE ORIGIN IN THE POLYHEDRON DEFINED BY X + Y + Z .GE. 2.0D0 (I.E. C -X - Y - Z + 2.0D0 .LE. 0.0D0) AND Z .GE. 1.0D0 (I.E. -Z + 1.0D0 .LE. C 0.0D0). (THE EXACT SOLUTION IS (X,Y,Z) = (0.5D0,0.5D0,1.0D0) WITH C DISTANCE SQRT(1.5D0) FROM THE ORIGIN.) C C SAMPLE DRIVER FOR (D) LEAST-DISTANCE QUADRATIC PROGRAMMING C C IMPLICIT REAL*8 (A-H,O-Z) CC THE MINIMUM DIMENSIONS ARE PMAT(NPARM+1,NUMGR), PMAT1(NPARM+1,NUMGR), CC WPT(NPARM), ICOR(NPARM+1), R(NPARM+1), PTNR(NPARM+1), COEF(NUMGR), CC WCOEF(NUMGR), IWORK(4*NUMGR+5*NPARM+3), WORK(2*NPARM**2+4*NUMGR*NPARM CC +9*NUMGR+22*NPARM+10), WHERE NPARM .GE. N AND NUMGR .GE. M. THE FIRST CC DIMENSIONS OF PMAT AND PMAT1 MUST BE EXACTLY NPARM+1. CC*****BEGIN USER SETTABLE STATEMENTS 1 OF 1 C DIMENSION PMAT(4,2),PMAT1(4,2),WPT(3),ICOR(4),R(4),PTNR(4), C *COEF(2),WCOEF(2),IWORK(26),WORK(136) C N=3 C M=2 C PMAT(1,1)=-1.0D0 C PMAT(2,1)=-1.0D0 C PMAT(3,1)=-1.0D0 C PMAT(4,1)=2.0D0 C PMAT(1,2)=0.0D0 C PMAT(2,2)=0.0D0 C PMAT(3,2)=-1.0D0 C PMAT(4,2)=1.0D0 C NPARM=N C NUMGR=M CC*****END USER SETTABLE STATEMENTS 1 OF 1 C LIWRK=4*NUMGR+5*NPARM+3 C LWRK=2*NPARM**2+4*NUMGR*NPARM+9*NUMGR+22*NPARM+10 C MP1=M+1 C NP1=N+1 C OPEN(6,FILE='QPOUT') C WRITE(6,100)M,N C 100 FORMAT(/10H THERE ARE,I5,16H BOUNDARIES AND,I5, C *12H DIMENSIONS//30H THE BOUNDARY COEFFICIENTS ARE) C DO 300 J=1,M C WRITE(6,200)(PMAT(I,J),I=1,NP1) C 200 FORMAT(/(3E22.13)) C 300 CONTINUE C CALL WOLFE(N,M,PMAT,0,S,NCOR,ICOR,IWORK,LIWRK,WORK,LWRK, C *R,COEF,PTNR,PMAT1,NPARM,NUMGR,WCOEF,WPT,WDIST,NMAJ,NMIN, C *JFLAG) C WRITE(6,400)JFLAG,WDIST,(WPT(I),I=1,N) C 400 FORMAT(/30H AFTER WOLFE THE ERROR FLAG IS,I4// C *32H THE DISTANCE FROM THE ORIGIN IS,E22.13// C *26H THE POINT HAS COORDINATES//(3E22.13)) C STOP C END C C OUTPUT FOR (D) LEAST-DISTANCE QUADRATIC PROGRAMMING C C THERE ARE 2 BOUNDARIES AND 3 DIMENSIONS C C THE BOUNDARY COEFFICIENTS ARE C C -0.1000000000000E+01 -0.1000000000000E+01 -0.1000000000000E+01 C 0.2000000000000E+01 C C 0.0000000000000E+00 0.0000000000000E+00 -0.1000000000000E+01 C 0.1000000000000E+01 C C AFTER WOLFE THE ERROR FLAG IS 0 C C THE DISTANCE FROM THE ORIGIN IS 0.1224744871392E+01 C C THE POINT HAS COORDINATES C C 0.5000000000000E+00 0.5000000000000E+00 0.1000000000000E+01 C C C**********END OF USERS MANUAL FOR CONMAX. C C C THIS IS A TEST DRIVER PROGRAM FOR CONMAX. FOR A DESCRIPTION OF C CONMAX, PLEASE SEE THE CONMAX USERS GUIDE, WHICH APPEARS AT THE C BEGINNING OF THIS PACKAGE. FOR MORE INFORMATION ABOUT THE C EXAMPLE WHICH IS SET UP IN THIS DRIVER PROGRAM AND IN SUBROUTINE C FNSET, PLEASE SEE THE COMMENTS IN THESE TWO ROUTINES AS WELL AS C THE COMMENTS IN THE AFOREMENTIONED USERS GUIDE. C C THIS TEST DRIVER PROGRAM AND FNSET ARE SET UP TO CHOOSE REAL (DOUBLE C PRECISION) PARAMETERS A, B, C, AND D TO MINIMIZE C C MAX{ MAX (|F(X,Y) - (AX+BY+C)/(DX+1)|, |(AX+BY+C)/(DX+1)|) : C (X,Y) IN Z} C C SUBJECT TO THE CONSTRAINTS C C DX+1 .GE. EPS FOR (X,Y) IN Z, C C AND THE FIRST PARTIAL DERIVATIVE OF (AX+BY+C)/(DX+1) WITH RESPECT TO C X IS .LE. 0.0 AT (X,Y) = (0.0,0.0). C C HERE WE ARE TAKING Z = {(0.0,0.0), (0.0,1.0), (-2.0/3.0,1.0/3.0), C (1.0,1.0), (1.0,2.0)}, EPS = .0001, F(0.0,0.0) = .5, F(0.0,1.0) = 1.0, C F(-2.0/3.0,1.0/3.0) = -1.0, F(1.0,1.0) =1.5, F(1.0,2.0) =-ONE, C AND TAKING THE INITIAL GUESSES FOR THE PARAMETERS TO BE C A = B = C = D = 0.0 C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION FUN(5),PTTBL(6,2),IWORK(178),WORK(720),PARAM(4), *ERROR(24) C C*****MAC INSERT C OPEN(6,FILE='TSTOT') C*****END MAC INSERT C C SET MACHINE AND PRECISION DEPENDENT CONSTANTS. ONE=1.0D0 ZERO=ONE-ONE TWO=ONE+ONE THREE=ONE+TWO NREAD=I1MACH(1) NWRIT=I1MACH(2) C C SET PARAMETERS FOR CONMAX. C C SET IOPTN=0 SINCE NO EXTRA OPTIONS ARE TO BE USED. IOPTN=0 C C SET NPARM=4 SINCE THERE ARE 4 PARAMETERS (VARIABLES) A, B, C, AND D. NPARM=4 C C SET NUMGR=21 SINCE THERE ARE 21 CONSTRAINTS (5 OF TYPE 2, 10 OF TYPE 1, C 5 OF TYPE -1, AND ONE OF TYPE -2). NUMGR=21 C C SET ITLIM=100, OR WHATEVER LIMIT IS DESIRED ON THE NUMBER OF ITERATIONS. ITLIM=100 C C SET IFUN=5 (OR GREATER) SINCE THERE ARE 5 TYPE 2 CONSTRAINTS, AND USE C THIS NUMBER AS THE DIMENSION OF FUN ABOVE. IFUN=5 C C SET IPTB=6 (OR GREATER) AND SET INDM=2 (OR GREATER), SINCE THE GREATEST C FIRST SUBSCRIPT WE WILL USE IN PTTBL IS 6, AND THE GREATEST SECOND C SUBSCRIPT WE WILL USE IN PTTBL IS 2. WE USE THESE NUMBERS TO DIMENSION C PTTBL ABOVE. NOTE THAT IT IS ESSENTIAL THAT THE FIRST DIMENSION OF C PTTBL BE EXACTLY IPTB. IPTB=6 INDM=2 C C SET LIWRK=178 (OR GREATER), AND USE THIS NUMBER TO DIMENSION LIWRK C ABOVE, BECAUSE OF THE COMPUTATION 7*21 + 7*4 + 3 = 178 (SEE CONMAX C USERS GUIDE FOR AN EXPLANATION OF THIS AND OF LWRK BELOW). LIWRK=178 C C SET LWRK=720 (OR GREATER), AND USE THIS NUMBER TO DIMENSION LWRK ABOVE, C BECAUSE OF THE COMPUTATION 2*4**2 + 4*21*4 + 11*21 + 27*4 + 13 = 720. LWRK=720 C C THE DIMENSION OF PARAM ABOVE MUST BE NPARM (I.E. 4) OR GREATER, AND THE C DIMENSION OF ERROR ABOVE MUST BE NUMGR+3 (I.E. 24) OR GREATER. C C SET THE VALUES OF THE FUNCTION F. FUN(1)=ONE/TWO FUN(2)=ONE FUN(3)=-ONE FUN(4)=THREE/TWO FUN(5)=-ONE C C SET THE COORDINATES OF THE FIVE POINTS. PTTBL(1,1)=ZERO PTTBL(1,2)=ZERO PTTBL(2,1)=ZERO PTTBL(2,2)=ONE PTTBL(3,1)=-TWO/THREE PTTBL(3,2)=ONE/THREE PTTBL(4,1)=ONE PTTBL(4,2)=ONE PTTBL(5,1)=ONE PTTBL(5,2)=TWO C C PUT EPS IN PTTBL(6,1) FOR TRANSMITTAL TO FNSET. THIS IS WHY WE NEEDED C IPTB TO BE AT LEAST 6. EPS=(5*TWO)**(-4) PTTBL(6,1)=EPS C C SET THE INITIAL GUESSES FOR THE PARAMETERS. PARAM(1)=ZERO PARAM(2)=ZERO PARAM(3)=ZERO PARAM(4)=ZERO C C WRITE THE INITIAL DATA. WRITE(NWRIT,100)IOPTN,NPARM,NUMGR,ITLIM,IFUN,IPTB,INDM 100 FORMAT(/9H IOPTN IS,I5,10H NPARM IS,I4,10H NUMGR IS,I4// *9H ITLIM IS,I5,9H IFUN IS,I3,9H IPTB IS,I3,9H INDM IS,I3) WRITE(NWRIT,200)(FUN(I),I=1,5) 200 FORMAT(/24H THE FUNCTION VALUES ARE/(4E15.5)) WRITE(NWRIT,300) 300 FORMAT(/15H THE POINTS ARE) DO 500 I=1,5 WRITE(NWRIT,400)(PTTBL(I,J),J=1,INDM) 400 FORMAT(2E15.5) 500 CONTINUE WRITE(NWRIT,600)EPS 600 FORMAT(/7H EPS IS,E15.5) WRITE(NWRIT,700)(PARAM(J),J=1,NPARM) 700 FORMAT(/27H THE INITIAL PARAMETERS ARE/(/3E23.14)) C C NOW CALL CONMAX. CALL CONMAX(IOPTN,NPARM,NUMGR,ITLIM,FUN,IFUN,PTTBL,IPTB, *INDM,IWORK,LIWRK,WORK,LWRK,ITER,PARAM,ERROR) C C WRITE THE OUTPUT. C NOTE THAT WE HAVE DEFERRED WRITING LIWRK AND LWRK UNTIL AFTER CALLING C CONMAX SINCE CONMAX WILL CHANGE THEM TO THE NEGATIVE OF THE SMALLEST C ALLOWABLE VALUES AND RETURN IF THEY WERE TOO SMALL. WRITE(NWRIT,800)ITER,LIWRK,LWRK 800 FORMAT(/26H *****AFTER CONMAX ITER IS,I4,10H LIWRK IS,I5, *9H LWRK IS,I6) WRITE(NWRIT,900)(PARAM(J),J=1,NPARM) 900 FORMAT(/25H THE FINAL PARAMETERS ARE/(/3E23.14)) WRITE(NWRIT,1000)ERROR(NUMGR+1),ERROR(NUMGR+2),ERROR(NUMGR+3), *(ERROR(I),I=1,NUMGR) 1000 FORMAT(/20H THE ERROR NORMS ARE//3E23.13// *15H THE ERRORS ARE//(3E23.14)) STOP END SUBROUTINE FNSET(NPARM,NUMGR,PTTBL,IPTB,INDM,PARAM,IPT, *INDFN,ICNTYP,CONFUN) C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION PTTBL(IPTB,INDM),PARAM(NPARM),ICNTYP(NUMGR), *CONFUN(NUMGR,NPARM+1) C C THIS IS THE SUBROUTINE FNSET FOR THE EXAMPLE DISCUSSED IN THE CONMAX C USERS GUIDE. C C SET PRECISION DEPENDENT CONSTANTS. ONE=1.0D0 ZERO=ONE-ONE C C WE BREAK FNSET INTO SECTIONS BASED ON THE VALUE OF IPT, THAT IS, ON C WHICH CONSTRAINT IS BEING SET. IF(IPT-5)300,300,100 100 IF(IPT-15)600,600,200 200 IF(IPT-20)1100,1100,1300 C C HERE IPT .LE. 5 AND WE SET A CONSTRAINT OF THE FORM C ABS(F(X,Y) - (AX+BY+C)/(DX+1)) .LE. W. C NOTE THAT SINCE THIS IS A TYPE 2 CONSTRAINT WE DO NOT NEED TO DEAL C WITH THE ABSOLUTE VALUE OR THE F(X,Y) HERE. 300 ICNTYP(IPT)=2 A=PARAM(1) B=PARAM(2) C=PARAM(3) D=PARAM(4) X=PTTBL(IPT,1) Y=PTTBL(IPT,2) P=A*X+B*Y+C Q=D*X+ONE CONFUN(IPT,1)=P/Q IF(INDFN)400,400,500 C 400 RETURN C C HERE IPT .LE. 5 AND INDFN=1, AND WE SET THE PARTIAL DERIVATIVES. 500 CONFUN(IPT,2)=X/Q CONFUN(IPT,3)=Y/Q CONFUN(IPT,4)=ONE/Q CONFUN(IPT,5)=-P*X/(Q*Q) RETURN C C HERE 6 .LE. IPT .LE. 15 AND IF IPT IS EVEN WE SET THE CONSTRAINT C (AX+BY+C)/(DX+1) .LE. W, WHICH IS HALF OF THE CONSTRAINT C ABS((AX+BY+C)/(DX+1)) .LE. W, WHILE IF IPT IS ODD WE SET THE CONSTRAINT C -(AX+BY+C)/(DX+1) .LE. W, WHICH IS THE OTHER HALF OF THE CONSTRAINT C ABS((AX+BY+C)/(DX+1)) .LE. W. 600 ICNTYP(IPT)=1 II=(IPT-4)/2 A=PARAM(1) B=PARAM(2) C=PARAM(3) D=PARAM(4) X=PTTBL(II,1) Y=PTTBL(II,2) P=A*X+B*Y+C Q=D*X+ONE IREM=IPT-4-2*II IF(IREM)700,700,900 C C HERE 6 .LE. IPT .LE. 15 AND IPT IS EVEN. 700 CONFUN(IPT,1)=P/Q IF(INDFN)400,400,800 800 CONFUN(IPT,2)=X/Q CONFUN(IPT,3)=Y/Q CONFUN(IPT,4)=ONE/Q CONFUN(IPT,5)=-P*X/(Q*Q) RETURN C C HERE 6 .LE. IPT .LE. 15 AND IPT IS ODD. 900 CONFUN(IPT,1)=-P/Q IF(INDFN)400,400,1000 1000 CONFUN(IPT,2)=-X/Q CONFUN(IPT,3)=-Y/Q CONFUN(IPT,4)=-ONE/Q CONFUN(IPT,5)=P*X/(Q*Q) RETURN C C HERE 16 .LE. IPT .LE. 20 AND WE SET A CONSTRAINT OF THE FORM C -DX - 1.0 + EPS .LE. 0.0 1100 ICNTYP(IPT)=-1 D=PARAM(4) EPS=PTTBL(6,1) II=IPT-15 X=PTTBL(II,1) CONFUN(IPT,1)=-D*X-ONE+EPS IF(INDFN)400,400,1200 1200 CONFUN(IPT,2)=ZERO CONFUN(IPT,3)=ZERO CONFUN(IPT,4)=ZERO CONFUN(IPT,5)=-X RETURN C C HERE IPT=21 AND WE SET THE CONSTRAINT C (PARTIAL DERIVATIVE OF (AX+BY+C)/(DX+1) WITH RESPECT TO X AT C (X,Y) = (0.0,0.0)) .LE. 0.0, C I.E. A - CD .LE. 0.0 1300 ICNTYP(IPT)=-2 A=PARAM(1) C=PARAM(3) D=PARAM(4) CONFUN(IPT,1)=A-C*D IF(INDFN)400,400,1400 1400 CONFUN(IPT,2)=ONE CONFUN(IPT,3)=ZERO CONFUN(IPT,4)=-D CONFUN(IPT,5)=-C RETURN END C C C*****END OF SAMPLE LISTING, START OF CONMAX SUBPROGRAMS C C SUBROUTINE CONMAX(IOPTN,NPARM,NUMGR,ITLIM,FUN,IFUN,PTTBL, *IPTB,INDM,IWORK,LIWRK,WORK,LWRK,ITER,PARAM,ERROR) C C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION FUN(IFUN),PTTBL(IPTB,INDM),ERROR(NUMGR+3), *PARAM(NPARM),IWORK(LIWRK),WORK(LWRK) C C C********** COPYRIGHT 1996 EDWIN H. KAUFMAN JR., DAVID J. LEEMING, C********** GERALD D. TAYLOR C********** THE AUTHORS GRATEFULLY ACKNOWLEDGE THE ASSISTANCE OF C********** CENTRAL MICHIGAN UNIVERSITY, THE UNIVERSITY OF VICTORIA C********** (CANADA), AND COLORADO STATE UNIVERSITY. C********** PERMISSION TO USE, COPY, MODIFY, AND DISTRIBUTE THIS C********** SOFTWARE FOR ANY PURPOSE WITHOUT FEE IS HEREBY GRANTED, C********** PROVIDED THAT THIS ENTIRE NOTICE IS INCLUDED IN ANY C********** SOFTWARE WHICH IS OR INCLUDES A COPY OR MODIFICATION OF C********** THIS SOFTWARE AND IN ALL COPIES OF THE SUPPORTING C********** DOCUMENTATION FOR SUCH SOFTWARE. C********** THIS SOFTWARE IS PROVIDED "AS IS", WITHOUT EXPRESS OR C********** IMPLIED WARRANTY. IN PARTICULAR, NEITHER THE AUTHORS NOR C********** THEIR UNIVERSITIES MAKE ANY REPRESENTATION OR WARRANTY OF C********** ANY KIND CONCERNING THE MERCHANTIBILITY OF THIS SOFTWARE C********** OR ITS FITNESS FOR ANY PARTICULAR PURPOSE. C C C PLEASE SEE THE USERS GUIDE FOR CONMAX AT THE BEGINNING OF THIS C PACKAGE FOR MORE INFORMATION ABOUT THE USE OF THESE SUBPROGRAMS. C C C CHECK TO SEE IF THE DIMENSIONS LIWRK AND LWRK ARE LARGE ENOUGH. IF C EITHER IS NOT, REPLACE IT BY THE NEGATIVE OF ITS CORRECT MINIMUM VALUE C AND RETRUN. JIWRK=7*NUMGR+7*NPARM+3 JWRK=2*NPARM**2+4*NUMGR*NPARM+11*NUMGR+27*NPARM+13 IF(LIWRK-JIWRK)10,20,20 10 LIWRK=-JIWRK IF(LWRK-JWRK)30,40,40 20 IF(LWRK-JWRK)30,50,50 30 LWRK=-JWRK 40 RETURN C C SET MACHINE AND PRECISION DEPENDENT CONSTANTS. 50 ONE=1.0D0 ZERO=ONE-ONE TWO=ONE+ONE FOUR=TWO+TWO TEN=FOUR+FOUR+TWO C NWRIT=I1MACH(2) SPCMN=D1MACH(3) C C INITIALIZE SOME OTHER PARAMETERS. NPAR1=NPARM+1 ISUCC=0 ITER=0 ITERSL=0 ITLIM1=ITLIM ENCHG=ZERO ILC02=ILOC(2,NPARM,NUMGR) ILC06=ILOC(6,NPARM,NUMGR) ILC08=ILOC(8,NPARM,NUMGR) ILC11=ILOC(11,NPARM,NUMGR) ILC12=ILOC(12,NPARM,NUMGR) ILC13=ILOC(13,NPARM,NUMGR) ILC14=ILOC(14,NPARM,NUMGR) ILC15=ILOC(15,NPARM,NUMGR) ILC17=ILOC(17,NPARM,NUMGR) ILC20=ILOC(20,NPARM,NUMGR) ILC21=ILOC(21,NPARM,NUMGR) ILC22=ILOC(22,NPARM,NUMGR) ILC24=ILOC(24,NPARM,NUMGR) ILC25=ILOC(25,NPARM,NUMGR) ILC26=ILOC(26,NPARM,NUMGR) ILC27=ILOC(27,NPARM,NUMGR) ILC29=ILOC(29,NPARM,NUMGR) ILC30=ILOC(30,NPARM,NUMGR) ILC31=ILOC(31,NPARM,NUMGR) ILC33=ILOC(33,NPARM,NUMGR) ILC35=ILOC(35,NPARM,NUMGR) ILC40=ILOC(40,NPARM,NUMGR) ILC42=ILOC(42,NPARM,NUMGR) ILC44=ILOC(44,NPARM,NUMGR) ILC46=ILOC(46,NPARM,NUMGR) C C IF THE TENS DIGIT OF IOPTN IS 1, SET KNTSM TO 0 AND GET ENCSM C FROM WORK(1) AND LIMSM FROM IWORK(1). IOPTEN=(IOPTN-(IOPTN/100)*100)/10 IF(IOPTEN)53,53,52 52 KNTSM=0 ENCSM=WORK(1) LIMSM=IWORK(1) C C IF THE HUNDREDS DIGIT OF IOPTN IS 1 OR 3, SET NSTEP = IWORK(2), C AND OTHERWISE SET NSTEP TO ITS DEFAULT VALUE OF 1. 53 IOPHUN=(IOPTN-(IOPTN/1000)*1000)/100 IF(IOPHUN-(IOPHUN/2)*2)55,55,54 54 NSTEP=IWORK(2) GO TO 56 55 NSTEP=1 C C IF THE HUNDREDS DIGIT OF IOPTN IS 2 OR 3, SET TOLCON = WORK(2), C AND OTHERWISE SET TOLCON TO ITS DEFAULT VALUE OF SQRT(SPCMN). 56 IF(IOPHUN-2)58,57,57 57 TOLCON=WORK(2) GO TO 60 58 TOLCON=SQRT(SPCMN) C C IN THIS VERSION OF CONMAX WE SET THE LINEAR CONSTRAINT TOLERANCE C EQUAL TO THE NONLINEAR CONSTRAINT TOLERANCE. 60 TOLLIN=TOLCON C C SET IRK=1 IF THE THOUSANDS DIGIT OF IOPTN IS 0 AND OTHERWISE SET IRK=0. IOPTHO=(IOPTN-(IOPTN/10000)*10000)/1000 IF(IOPTHO)100,100,120 100 IRK=1 GO TO 200 120 IRK=0 C C COMPUTE THE TEN THOUSANDS DIGIT OF IOPTN FOR LATER USE. 200 IOPTTH=(IOPTN-(IOPTN/100000)*100000)/10000 C C SET IPHSE=-1 TO INDICATE WE HAVE NOT CHECKED TYPE -1 FEASIBILITY YET. IPHSE=-1 C SET RCHDWN = THE NUMBER OF LENGTHS OF PROJCT IN RKSACT (OR NUMBER OF C LENGTHS OF BNDLGT IN SETU1) WE WILL GO BELOW ERROR(NUMGR+1) TO DECLARE C A PRIMARY CONSTRAINT TO BE ACTIVE. RCHDWN=TWO RCHDNK=RCHDWN C SET RCHIN = THE NUMBER OF LENGTHS OF PROJCT (OR BNDLGT) WE WILL GO C BELOW 0.0 TO DECLARE A TYPE -2 CONSTRAINT TO BE ACTIVE. RCHIN=TWO C SET A NORMAL VALUE FOR NUMLIM FOR USE IN SLPCON. NUMLIM=11 C C END OF PRELIMINARY SECTION. THE STATEMENTS ABOVE THIS POINT WILL NOT C BE EXECUTED AGAIN IN THIS CALL TO CONMAX. C C C CALL ERCMP1 WITH ICNUSE=0 TO COMPUTE THE ERRORS, ERROR NORMS, AND ICNTYP. C WE TAKE IPHSE AS 0 SO ALL CONSTRAINTS WILL BE COMPUTED BY FNSET IN CASE C THE TEN THOUSANDS DIGIT OF IOPTN IS 1. C THIS IS ONE OF ONLY TWO PLACES IN THE PROGRAM WHERE WE CALL ERCMP1 WITH C ICNUSE=0, THE OTHER BEING STATEMENT 1415 BELOW.. 500 CALL ERCMP1(IOPTN,NPARM,NUMGR,FUN,IFUN,PTTBL,IPTB,INDM,PARAM, *0,0,IWORK,LIWRK,WORK(ILC08),IWORK(ILC17),IPMAX,ISMAX,ERROR) C IF ITLIM=0 WE RETURN. IF(ITLIM)510,510,520 510 RETURN C C COMPUTE ITYP2, ITYP1, ITYPM1, AND ITYPM2 AS THE NUMBER OF CONSTRAINTS OF C TYPE 2 (I.E. PRIMARY, ABS(FUN(I)-CONFUN(I,1)) .LE. W) OR 1 (I.E. PRIMARY, C CONFUN(I,1) .LE. W) OR -1 (I.E. STANDARD LINEAR, CONFUN(I,1) .LE. 0.0) C OR -2 (I.E. STANDARD NONLINEAR) RESPECTIVELY. 520 ITYP2=0 ITYP1=0 ITYPM1=0 ITYPM2=0 C NOTE THAT ARRAYS NOT IN THE CALLING SEQUENCE FOR CONMAX ARE ACCESSED C THROUGH THEIR LOCATION IN IWORK OR WORK. CONMAX IS THE ONLY C SUBROUTINE IN WHICH THIS IS NECESSARY. DO 900 I=1,NUMGR II=ILC17-1+I C HERE IWORK(II)=ICNTYP(I). IF(IWORK(II))600,900,550 550 IF(IWORK(II)-1)585,585,570 570 ITYP2=ITYP2+1 GO TO 900 585 ITYP1=ITYP1+1 GO TO 900 600 IF(IWORK(II)+1)800,700,700 700 ITYPM1=ITYPM1+1 GO TO 900 800 ITYPM2=ITYPM2+1 900 CONTINUE C C COMPUTE THE ERROR NORMS. ENORM IS THE PRINCIPAL ERROR NORM. 1000 ENORM=ERROR(NUMGR+1) ENOR2=ERROR(NUMGR+2) ENOR3=ERROR(NUMGR+3) C C WRITE ITER, ISUCC, IRK, ENCHG, AND THE ERROR NORMS. 1050 CONTINUE C1050 WRITE(NWRIT,1100)ITER,ISUCC,IRK,ENCHG,ENORM,ENOR2,ENOR3 C WRITE(9,1100)ITER,ISUCC,IRK,ENCHG,ENORM,ENOR2,ENOR3 C1100 FORMAT(/8H ITER IS,I5,10H ISUCC IS,I4,8H IRK IS,I4, C *10H ENCHG IS,E24.14/9H ENORM IS,E24.14,10H ENOR2 IS,E24.14/ C *9H ENOR3 IS,E24.14) C C C THE NEXT SECTION DETERMINES WHETHER WE WILL TERMINATE DUE TO ITERATION C COUNT, AND IF SO FOR OUTPUT PURPOSES IT MODIFIES ITER (OR TWO OF THE C ERROR NORMS IF THE FAILURE IS DUE TO INABILITY TO GAIN TYPE -2 C FEASIBILITY). C C IF IOPTEN=1 AND WE HAVE DONE AT LEAST ONE ITERATION IN THE MAIN PART C OF CONMAX, WE WILL GIVE UP IF ABS(ENCHG) HAS BEEN LESS THAN ENCSM FOR C LIMSM CONSECUTIVE MAIN ITERATIONS (INCLUDING THIS ONE). IF(IOPTEN-1)1118,1106,1118 1106 IF(IPHSE)1118,1108,1118 1108 IF(ITER)1118,1118,1110 1110 IF(-ENCHG-ENCSM)1114,1112,1112 1112 KNTSM=0 GO TO 1118 1114 KNTSM=KNTSM+1 IF(KNTSM-LIMSM)1118,1200,1200 C 1118 IF(ITER-ITLIM1)1300,1120,1120 C C HERE ITER = ITLIM1, SO WE RETURN. 1120 IF(IPHSE)1140,1200,1200 C C HERE WE HAVE FAILED TO ACHIEVE TYPE -2 FEASIBILITY AND WE SET ITER=-2 C AS A WARNING, PUT ERROR(NUMGR+1) IN ITS PROPER LOCATION, SET C ERROR(NUMGR+1) = 0.0 SINCE THE PRIMARY CONSTRAINTS WERE NOT COMPUTED, C AND RETURN. NOTE THAT WE CANNOT HAVE IPHSE=-1 HERE SINCE THAT WOULD C IMPLY ITER=0, THUS ITLIM=ITLIM1=0, IN WHICH CASE WE WOULD HAVE C TERMINATED EARLIER. 1140 ITER=-2 ERROR(NUMGR+3)=ERROR(NUMGR+1) ERROR(NUMGR+1)=ZERO C WRITE(6,1150) C1150 FORMAT(43H ***WARNING NONLINEAR STANDARD FEASIBILITY, C *16H NOT ACHIEVED***) RETURN C C FOR OUTPUT PURPOSES REPLACE ITER BY ITER + ITLIM - ITLIM1, THE TRUE C NUMBER OF ITERATIONS COUNTING INITIALIZATION. ITLIM - ITLIM1 WILL BE C THE NUMBER OF ITERATIONS NEEDED TO GAIN TYPE -2 FEASIBILITY. WORK C DONE TO GAIN TYPE -1 FEASIBILITY IS NOT COUNTED AS AN ITERATION. 1200 ITER=ITER+ITLIM-ITLIM1 C 1205 RETURN C C HERE ITER .LT. ITLIM1. IF IPHSE = 0 OR -2 HERE WE GO INTO THE C ITERATIVE PHASE OF CONMAX. 1300 IF(IPHSE+1)1450,1302,1450 C C C HERE IPHSE=-1 AND WE CHECK TYPE -1 FEASIBILITY, TRY TO REGAIN IT IF C WE DONT HAVE IT, CHECK TYPE -2 FEASIBILITY, AND SET UP FOR TYPE -2 C FEASIBILITY ITERATIONS IF WE DONT HAVE IT. THE STATEMENTS FROM HERE C DOWN TO THE TRIPLE BLANK LINE WILL BE EXECUTED AT MOST ONCE. C C NOTE THAT ENOR2=0.0 IF THERE ARE NO TYPE -1 CONSTRAINTS. 1302 IF(ENOR2-TOLLIN)1304,1304,1316 C C HERE WE HAD TYPE -1 FEASIBILITY INITIALLY. 1304 IF(ENOR3-TOLCON)1444,1444,1430 C C HERE WE DO NOT HAVE TYPE -1 FEASIBILITY SO WE TRY TO GET IT. C WE WILL NEED TO TELL DERST TO COMPUTE THE VALUES OF THE LEFT SIDES C OF THE TYPE -1 CONSTRAINTS WITH THE VARIABLES EQUAL TO ZERO (I.E. C THE CONSTANT TERMS IN THE CONSTRAINTS), SO WE SET PARWRK TO THE C ZERO VECTOR TO CARRY THE MESSAGE. 1316 DO 1324 J=1,NPARM JJ=ILC27-1+J C HERE WORK(JJ) = PARWRK(J). WORK(JJ)=ZERO 1324 CONTINUE IF(IOPTTH)1328,1328,1326 C HERE IOPTTH=1 AND WE CALL DERST WITH IPT=-1 TO PUT ALL THE STANDARD C CONSTRAINT AND DERIVATIVE VALUES IN CONFUN. C WE SET IPT=-1 TO TELL DERST IT NEED ONLY COMPUTE STANDARD CONSTRAINTS. 1326 IPT=-1 CALL DERST(IOPTN,NPARM,NUMGR,PTTBL,IPTB,INDM,WORK(ILC27),IPT, *WORK(ILC24),WORK(ILC35),IWORK(ILC22),WORK(ILC08)) C 1328 M=0 DO 1350 I=1,NUMGR II=ILC17-1+I C HERE WE CONSIDER ONLY TYPE -1 CONSTRAINTS. THERE MUST BE AT LEAST C ONE OF THESE, SINCE OTHERWISE WE WOULD NOT BE HERE ATTEMPTING TO C GAIN TYPE -1 FEASIBILITY. C HERE IWORK(II)=ICNTYP(I). IF(IWORK(II)+1)1350,1330,1350 1330 M=M+1 IF(IOPTTH)1332,1332,1335 C HERE IOPTTH=0 AND WE HAVE NOT YET CALLED DERST TO PUT CONSTRAINT I C AND ITS DERIVATIVES IN CONFUN, SO WE DO IT NOW. 1332 IPT=I CALL DERST(IOPTN,NPARM,NUMGR,PTTBL,IPTB,INDM,WORK(ILC27),IPT, * WORK(ILC24),WORK(ILC35),IWORK(ILC22),WORK(ILC08)) C COPY THE DERIVATIVES INTO PMAT FOR USE BY WOLFE. 1335 DO 1340 L=1,NPARM L1=ILC29-1+L+(M-1)*NPAR1 L2=ILC08-1+I+L*NUMGR C HERE WORK(L1)=PMAT(L,M) AND WORK(L2)=CONFUN(I,L+1). WORK(L1)=WORK(L2) 1340 CONTINUE C C NOW THE ITH CONSTRAINT (WHICH IS ALSO THE MTH TYPE -1 CONSTRAINT) HAS C THE FORM PMAT(1,M)*Z1+...+PMAT(NPARM,M)*ZNPARM + CONFUN(I,1) .LE. C 0.0. WE MAKE THE CHANGE OF VARIABLES ZZ = Z - PARAM TO TRANSLATE THE C ORIGIN TO PARAM. THE ITH CONSTRAINT WILL THEN HAVE THE FORM C PMAT(1,M)*ZZ1+...+PMAT(NPARM,M)*ZZNPARM + (CONFUN(I,1) + PMAT(1,M)* C PARAM(1)+...+PMAT(NPARM,M)*PARAM(NPARM)) .LE. 0.0. AFTER WOLFE FINDS C THE CLOSEST POINT TO THE ORIGIN IN THE POLYHEDRON DEFINED BY THE NEW C CONSTRAINTS, WE WILL ADD PARAM TO TRANSLATE BACK TO THE POINT WE WANT. L1=ILC29-1+NPAR1+(M-1)*NPAR1 L2=ILC08-1+I C HERE WORK(L1)=PMAT(NPAR1,1) AND WORK(L2)=CONFUN(I,1). WORK(L1)=WORK(L2) DO 1345 L=1,NPARM L2=ILC29-1+L+(M-1)*NPAR1 C HERE WORK(L1)=PMAT(NPAR1,1) AND WORK(L2)=PMAT(L,M). WORK(L1)=WORK(L1)+WORK(L2)*PARAM(L) 1345 CONTINUE 1350 CONTINUE C CALL WOLFE WITH ISTRT=0 TO COMPUTE THE SOLUTION IN THE ZZ COORDINATE C SYSTEM FROM SCRATCH. CALL WOLFE(NPARM,M,WORK(ILC29),0,S,NCOR,IWORK(ILC15),IWORK,LIWRK, *WORK,LWRK,WORK(ILC33),WORK(ILC06),WORK(ILC31),WORK(ILC30),NPARM, *NUMGR,WORK(ILC40),WORK(ILC42),WDIST,NMAJ,NMIN,JFLAG) IF(JFLAG)1365,1365,1355 C C HERE WE HAVE FAILED TO ACHIEVE TYPE -1 FEASIBILITY. WE SET ITER=-1 C AS A WARNING AND RETURN. 1355 ITER=-1 C WRITE(NWRIT,1360) C1360 FORMAT(40H ***WARNING LINEAR STANDARD FEASIBILITY, C *16H NOT ACHIEVED***) RETURN C C HERE JFLAG .LE. 0 AND WE PUT PARAM+WPT IN PARWRK TO CHECK WHETHER C THE TYPE -1 CONSTRAINTS ARE NOW FEASIBLE WITHIN TOLLIN. 1365 DO 1370 J=1,NPARM J1=ILC27-1+J J2=ILC42-1+J C HERE WORK(J1)=PARWRK(J) AND WORK(J2)=WPT(J). WORK(J1)=PARAM(J)+WORK(J2) 1370 CONTINUE C FOR USE IN ERCMP1 WE SET JCNTYP(I)=-1 IF ICNTYP(I)=-1 AND SET C JCNTYP(I)=0 OTHERWISE. DO 1385 I=1,NUMGR II=ILC17-1+I JJ=ILC21-1+I C HERE IWORK(II)=ICNTYP(I) AND IWORK(JJ)=JCNTYP(I). IF(IWORK(II)+1)1380,1375,1380 1375 IWORK(JJ)=-1 GO TO 1385 1380 IWORK(JJ)=0 1385 CONTINUE C CALL ERCMP1 WITH ICNUSE=1. CALL ERCMP1(IOPTN,NPARM,NUMGR,FUN,IFUN,PTTBL,IPTB,INDM, *WORK(ILC27),1,IPHSE,IWORK,LIWRK,WORK(ILC08),IWORK(ILC21), *IPMAX,ISMAX,WORK(ILC11)) I1=ILC11-1+(NUMGR+2) C HERE WORK(I1)=ERR1(NUMGR+2). IF(WORK(I1)-TOLLIN)1390,1390,1355 C C HERE WE HAVE ACHIEVED TYPE -1 FEASIBILITY. WE REPLACE PARAM WITH C PARWRK. 1390 DO 1395 J=1,NPARM JJ=ILC27-1+J C HERE WORK(JJ)=PARWRK(J). PARAM(J)=WORK(JJ) 1395 CONTINUE II=ILC11-1+NUMGR+2 C HERE WORK(II)=ERR1(NUMGR+2). C WRITE(NWRIT,1397)WORK(II),(PARAM(J),J=1,NPARM) C1397 FORMAT(48H TYPE -1 FEASIBILITY ACHIEVED. ERR1(NUMGR+2) IS, C *E15.5,10H PARAM IS/(4E20.12)) C C IF THERE ARE TYPE -2 CONSTRAINTS, SET JCNTYP AS ICNTYP WITH ALL BUT -2 C VALUES ZEROED OUT AND CALL ERCMP1 WITH ICNUSE=1 TO CHECK TYPE -2 C FEASIBILITY. WE CANNOT SIMPLY CHECK THE OLD ENOR3 HERE SINCE PARAM HAS C BEEN CHANGED. IF THERE ARE NO TYPE -2 CONSTRAINTS WE WILL AUTOMATICALLY C HAVE TYPE -2 FEASIBILITY. IF(ITYPM2)1415,1415,1398 1398 DO 1410 I=1,NUMGR II=ILC17-1+I JJ=ILC21-1+I C HERE IWORK(II)=ICNTYP(I) AND IWORK(JJ)=JCNTYP(I). IF(IWORK(II)+1)1400,1405,1405 1400 IWORK(JJ)=-2 GO TO 1410 1405 IWORK(JJ)=0 1410 CONTINUE CALL ERCMP1(IOPTN,NPARM,NUMGR,FUN,IFUN,PTTBL,IPTB,INDM, *PARAM,1,IPHSE,IWORK,LIWRK,WORK(ILC08),IWORK(ILC21),IPMAX, *ISMAX,WORK(ILC11)) II=ILC11-1+NUMGR+3 C HERE WORK(II)=ERR1(NUMGR+3). IF(WORK(II)-TOLCON)1415,1415,1430 C C HERE WE HAVE BOTH TYPE -1 AND TYPE -2 FEASIBILITY, BUT PARAM WAS C CHANGED IN GETTING TYPE -1 FEASIBILITY, SO WE CALL ERCMP1 C WITH ICNUSE=0 (ICNUSE=1 WOULD WORK ALSO SINCE ICNTYP HAS NOT BEEN C CHANGED HERE) TO GET THE NEW ERROR VECTOR. 1415 CALL ERCMP1(IOPTN,NPARM,NUMGR,FUN,IFUN,PTTBL,IPTB,INDM, *PARAM,0,IPHSE,IWORK,LIWRK,WORK(ILC08),IWORK(ILC17),IPMAX, *ISMAX,ERROR) GO TO 1444 C C HERE WE HAVE TYPE -1 FEASIBILITY BUT NOT TYPE -2 FEASIBILITY. WE SET C UP FOR THE TYPE -2 FEASIBILITY ITERATIONS, IN WHICH TYPE 1 AND TYPE C 2 CONSTRAINTS ARE IGNORED AND TYPE -2 CONSTRAINTS ARE TREATED AS C TYPE 1 CONSTRAINTS, EXCEPT WE WILL SWITCH OVER TO NORMAL ITERATIONS C ONCE WE CAN FORCE W .LE. TOLCON. THUS WE SET THE INDICATOR IPHSE TO C -2, RESET ICNTYP(I) TO 1 IF IT WAS -2, LEAVE IT AT -1 IF IT WAS -1, C AND SET IT TO 0 OTHERWISE, RESET ITYP2, ITYP1, AND ITYPM2, AND CALL C ERCMP1 WITH ICNUSE=1 TO PUT THE PROPER VALUES IN ERROR. 1430 IPHSE=-2 DO 1439 I=1,NUMGR II=ILC17-1+I C HERE IWORK(II)=ICNTYP(I). IF(IWORK(II)+1)1433,1439,1436 1433 IWORK(II)=1 GO TO 1439 1436 IWORK(II)=0 1439 CONTINUE C SAVE ITYP2 AND ITYP1. ITYP2K=ITYP2 ITYP1K=ITYP1 ITYP2=0 ITYP1=ITYPM2 ITYPM2=0 CALL ERCMP1(IOPTN,NPARM,NUMGR,FUN,IFUN,PTTBL,IPTB,INDM, *PARAM,1,IPHSE,IWORK,LIWRK,WORK(ILC08),IWORK(ILC17),IPMAX, *ISMAX,ERROR) GO TO 1450 C C HERE WE HAVE BOTH TYPE -1 AND TYPE -2 FEASIBILITY, AND WE C SET IPHSE=0 AND GO INTO THE MAIN PART OF CONMAX (UNLESS THERE WERE C NO TYPE 1 OR TYPE 2 CONSTRAINTS, IN WHICH CASE WE RETURN). 1444 IPHSE=0 IF(ITYP1+ITYP2)1205,1205,1450 C C END OF INITIAL FEASIBILITY CHECKING, TYPE -1 FEASIBILITY WORK, AND C TYPE -2 SETUP. THE BLOCK OF STATEMENTS FROM HERE UP TO THE C PRECEDING DOUBLE BLANK LINE WILL NOT BE EXECUTED AGAIN. C C C 1450 IF(IRK)1475,1475,1500 C C HERE IRK IS 0 OR -1 AND WE DO AN SLP STEP. IF SLPCON CANNOT REDUCE THE C PRINCIPAL ERROR NORM ENORM = ERROR(NUMGR+1) BY MORE THAN 100.0*B**(-ITT) C THEN IT WILL LEAVE PARAM AND ERROR UNCHANGED. 1475 CALL SLPCON(IOPTN,NPARM,NUMGR,FUN,IFUN,PTTBL,IPTB,INDM, *TOLCON,RCHIN,IRK,ITYPM1,ITYPM2,IWORK(ILC17),RCHDWN,NUMLIM,ITERSL, *PRJSLP,WORK(ILC12),IWORK(ILC20),WORK(ILC44),MACT1,IWORK(ILC14), *IWORK(ILC21),IPHSE,ENCHG,IWORK,LIWRK,WORK,LWRK,WORK(ILC26), *ISUCC,PARAM,ERROR) GO TO 1600 C C HERE IRK IS 1 OR 2 AND WE DO AN RK STEP. IF RKCON CANNOT REDUCE THE C PRINCIPAL ERROR NORM ENORM = ERROR(NUMGR+1) BY MORE THAN 100.0*B**(-ITT) C THEN IT WILL LEAVE PARAM AND ERROR UNCHANGED. 1500 CALL RKCON(IOPTN,NPARM,NUMGR,FUN,IFUN,PTTBL,IPTB,INDM, *TOLCON,RCHIN,ITER,IRK,ITYP2,ITYP1,ITYPM1,ITYPM2,IWORK(ILC17), *PROJCT,RCHDWN,NSTEP,IPHSE,ENCHG,ENC1,WORK(ILC29),WORK(ILC12), *IWORK,LIWRK,WORK,LWRK,IWORK(ILC13),WORK(ILC02),WORK(ILC25), *WORK(ILC26),WORK(ILC46),WORK(ILC11),WORK(ILC08),ISUCC,PARAM, *ERROR) C 1600 IF(ISUCC)1700,1700,2100 C HERE THE RK OR SLP STEP REDUCED ERROR(NUMGR+1) BY MORE THAN C 100.0*B**(-ITT), AND WE INCREMENT ITER. 1700 ITER=ITER+1 C C IF EITHER IPHSE=0, OR IPHSE=-2 AND ERROR(NUMGR+1) .GT. TOLCON, WE GO C ON AS USUAL TO SET UP ANOTHER STEP WITH THE SAME IPHSE. IF(IPHSE)1710,1790,1790 1710 IF(ERROR(NUMGR+1)-TOLCON)1720,1720,1790 C C HERE IPHSE=-2 AND ERROR(NUMGR+1) .LE. TOLCON, SO WE HAVE JUST ACHIEVED C TYPE -2 FEASIBILITY. WE WILL SET IPHSE=0, AND IF THERE ARE ANY C PRIMARY CONSTRAINTS WE WILL RESET ITER, ITERSL, AND ITLIM1 (SINCE C ITER=0 AND ITERSL=0 HAVE MEANINGS TO RKCON AND SLPCON RESPECTIVELY), C RESET RCHIN AND RCHDWN, AND GO BACK TO THE FIRST ERCMP1 CALL TO C RESTORE ERROR AND ICNTYP (ITYP1, ITYP2, ITYPM1, AND ITYPM2 WILL ALSO C BE RESTORED). 1720 IPHSE=0 IF(ITYP1K+ITYP2K)1205,1205,1730 1730 ITLIM1=ITLIM-ITER ITER=0 ITERSL=0 RCHIN=RCHDWN RCHDWN=RCHDNK GO TO 500 C 1790 IF(IRK)1800,1900,2000 C C HERE WE HAD AN SLP SUCCESS AND WE ARE GOING TO TRY RK AGAIN, SO WE SET C IRK=2 TO WARN RKCON THAT THE SUCCESS CAME FROM SLP. 1800 IRK=2 C HERE WE HAD AN SLP SUCCESS AND WE INCREMENT ITERSL = THE NUMBER OF SLP C SUCCESSES SINCE THE LAST SUCCESSFUL RK STEP (IF ANY). ITERSL IS NEEDED C IN SUBROUTINE BNDSET (CALLED BY SLPCON). 1900 ITERSL=ITERSL+1 GO TO 1000 C C HERE IRK IS 1 OR 2, SO WE JUST HAD AN RK SUCCESS. WE RESET IRK AND C ITERSL. 2000 IRK=1 ITERSL=0 GO TO 1000 C C HERE RKCON OR SLPCON FAILED TO SIGNIFICANTLY REDUCE THE PRINCIPAL ERROR C NORM. IF WE JUST TRIED SLP WE QUIT, AND IF WE JUST TRIED RK WE ATTEMPT C AN SLP STEP UNLESS IOPTHO = 2, IN WHICH CASE WE QUIT. 2100 IF(IRK)2300,2300,2150 2150 IF(IOPTHO-2)2200,2300,2200 2200 IRK=-1 GO TO 1050 C C IF IPHSE=-2 HERE WE WILL SET ITER=-2 AS A WARNING AND CHANGE C ERROR(NUMGR+1) AND ERROR(NUMGR+3) BEFORE RETURNING. OTHERWISE WE WILL C HAVE IPHSE=0 AND WE WILL ADJUST ITER BEFORE RETURNING. 2300 IF(IPHSE)1140,1200,1200 END C FUNCTION I1MACH(I) CC CC THIS IS THE FIRST OF TWO FUNCTION SUBPROGRAMS IN WHICH THE USER SETS CC MACHINE DEPENDENT CONSTANTS. IT SETS THE INPUT AND OUTPUT UNIT CC NUMBERS. CC CC I1MACH(1) IS THE INPUT UNIT NUMBER, AND I1MACH(2) IS THE OUTPUT CC UNIT NUMBER. C IF(I-1)1,1,2 C 1 I1MACH=5 C RETURN C 2 I1MACH=6 C RETURN C END C FUNCTION D1MACH(I) C CC***BEGIN PROLOGUE D1MACH CC***ROUTINES CALLED (NONE) CC***PURPOSE THIS IS THE SECOND OF TWO FUNCTION SUBPROGRAMS IN CC WHICH THE USER MUST SET MACHINE DEPENDENT CONSTANTS. CC IT SETS THE PRECISION DEPENDENT CONSTANT CC CC D1MACH(3) = B**(-ITT) CC CC WHERE B IS THE BASE FOR FLOATING POINT NUMBERS, AND CC ITT IS THE NUMBER OF BASE B DIGITS IN THE MANTISSA. CC***REMARK TO CONVERT THIS PROGRAM FROM DOUBLE PRECISION TO SINGLE CC PRECISION (FOR EXAMPLE), ON MANY MACHINES ONE NEED ONLY CC RESET D1MACH(3), PUT A C IN COLUMN 1 OF ALL THE STATEMENTS CC CC IMPLICIT REAL*8 (A-H,O-Z) CC CC AND CONVERT THE DRIVER PROGRAM AND SUBROUTINE FNSET TO CC SINGLE PRECISION. CC C IMPLICIT REAL*8 (A-H,O-Z) CC C IF(I-3)100,200,100 CC C 100 RETURN CC C 200 D1MACH=16.0D0**(-14) C RETURN C END FUNCTION ILOC(IARR,NPARM,NUMGR) C C THIS FUNCTION SUBPROGRAM RETURNS THE SUBSCRIPT OF THE FIRST ELEMENT OF C ARRAY IARR RELATIVE TO IWORK (IF THE ARRAY IS INTEGER, I.E. 13 .LE. C IARR .LE. 23) OR RELATIVE TO WORK (IF THE ARRAY IS FLOATING POINT, I.E. C 1 .LE. IARR .LE. 12 OR 24 .LE. IARR .LE. 48). C GO TO (10,20,30,40,50,60,70,80,90,100,110,120,130,140,150,160, *170,180,190,200,210,220,230,240,250,260,270,280,290,300,310,320, *330,340,350,360,370,380,390,400,410,420,430,440,450,460,470, *480),IARR C C 1 AA(NPARM+1,NPARM+1) (OPPOSITE V, Y; STARTS AT V STARTING POINT) 10 ILOC=3*NUMGR*NPARM+6*NUMGR+11*NPARM+8 RETURN C C 2 ACTDIF(NUMGR) 20 ILOC=1 RETURN C C 3 B(NPARM+1) (OPPOSITE V, Y; FOLLOWS AA) 30 ILOC=NPARM**2+3*NUMGR*NPARM+6*NUMGR+13*NPARM+9 RETURN C C 4 BETA(NPARM+1) (OPPOSITE V, Y; FOLLOWS B) 40 ILOC=NPARM**2+3*NUMGR*NPARM+6*NUMGR+14*NPARM+10 RETURN C C 5 BNDKP(NPARM) (FOLLOWS ACTDIF) 50 ILOC=NUMGR+1 RETURN C C 6 COEF(NUMGR) 60 ILOC=NUMGR+NPARM+1 RETURN C C 7 COFBND(NPARM) 70 ILOC=2*NUMGR+NPARM+1 RETURN C C 8 CONFUN(NUMGR,NPARM+1) (OPPOSITE PMAT1) 80 ILOC=2*NUMGR+2*NPARM+1 RETURN C C 9 D(NPARM+1) (OPPOSITE V, Y; FOLLOWS BETA) 90 ILOC=NPARM**2+3*NUMGR*NPARM+6*NUMGR+15*NPARM+11 RETURN C C 10 DVEC(NPARM) (FOLLOWS CONFUN) 100 ILOC=NUMGR*NPARM+3*NUMGR+2*NPARM+1 RETURN C C 11 ERR1(NUMGR+3) 110 ILOC=NUMGR*NPARM+3*NUMGR+3*NPARM+1 RETURN C C 12 FUNTBL(NUMGR,NPARM+1) 120 ILOC=NUMGR*NPARM+4*NUMGR+3*NPARM+4 RETURN C C 13 IACT(NUMGR) 130 ILOC=1 RETURN C C 14 IACT1(NUMGR) 140 ILOC=NUMGR+1 RETURN C C 15 ICOR(NPARM+1) 150 ILOC=2*NUMGR+1 RETURN C C 16 ICOR1(NPARM+1) (DOES NOT APPEAR IN PROGRAM BY NAME) 160 ILOC=2*NUMGR+NPARM+2 RETURN C C 17 ICNTYP(NUMGR) 170 ILOC=2*NUMGR+2*NPARM+3 RETURN C C 18 IXRCT(NUMGR+2*NPARM) 180 ILOC=3*NUMGR+2*NPARM+3 RETURN C C 19 IYCCT(NPARM+1) (OPPOSITE KPIVOT) 190 ILOC=4*NUMGR+4*NPARM+3 RETURN C C 20 IYRCT(NUMGR+2*NPARM) 200 ILOC=4*NUMGR+5*NPARM+4 RETURN C C 21 JCNTYP(NUMGR) 210 ILOC=5*NUMGR+7*NPARM+4 RETURN C C 22 KCNTYP(NUMGR) 220 ILOC=6*NUMGR+7*NPARM+4 RETURN C C 23 KPIVOT(NPARM+1) (OPPOSITE IYCCT) 230 ILOC=4*NUMGR+4*NPARM+3 RETURN C C 24 PARAM1(NPARM) (FOLLOWS FUNTBL) 240 ILOC=2*NUMGR*NPARM+5*NUMGR+3*NPARM+4 RETURN C C 25 PARPRJ(NPARM) 250 ILOC=2*NUMGR*NPARM+5*NUMGR+4*NPARM+4 RETURN C C 26 PARSER(NPARM) 260 ILOC=2*NUMGR*NPARM+5*NUMGR+5*NPARM+4 RETURN C C 27 PARWRK(NPARM) 270 ILOC=2*NUMGR*NPARM+5*NUMGR+6*NPARM+4 RETURN C C 28 PICOR(NPARM+1,NPARM+1) (OPPOSITE V, Y; FOLLOWS D) 280 ILOC=NPARM**2+3*NUMGR*NPARM+6*NUMGR+16*NPARM+12 RETURN C C 29 PMAT(NPARM+1,NUMGR) (FOLLOWS PARWRK) 290 ILOC=2*NUMGR*NPARM+5*NUMGR+7*NPARM+4 RETURN C C 30 PMAT1(NPARM+1,NUMGR) (OPPOSITE CONFUN) 300 ILOC=2*NUMGR+2*NPARM+1 RETURN C C 31 PTNR(NPARM+1) (FOLLOWS PMAT) 310 ILOC=3*NUMGR*NPARM+6*NUMGR+7*NPARM+4 RETURN C C 32 PTNRR(NPARM+1) 320 ILOC=3*NUMGR*NPARM+6*NUMGR+8*NPARM+5 RETURN C C 33 R(NPARM+1) 330 ILOC=3*NUMGR*NPARM+6*NUMGR+9*NPARM+6 RETURN C C 34 SAVE(NPARM+1) 340 ILOC=3*NUMGR*NPARM+6*NUMGR+10*NPARM+7 RETURN C C 35 V(NUMGR+2*NPARM+1,NPARM+2) (WITH Y, OPPOSITE AA, B, BETA, D, C PICOR, ZWORK) 350 ILOC=3*NUMGR*NPARM+6*NUMGR+11*NPARM+8 RETURN C C 36 VDER(NPARM) (FOLLOWS Y) 360 ILOC=2*NPARM**2+4*NUMGR*NPARM+9*NUMGR+18*NPARM+10 RETURN C C 37 VDERN(NPARM) 370 ILOC=2*NPARM**2+4*NUMGR*NPARM+9*NUMGR+19*NPARM+10 RETURN C C 38 VDERS(NPARM) 380 ILOC=2*NPARM**2+4*NUMGR*NPARM+9*NUMGR+20*NPARM+10 RETURN C C 39 VEC(NPARM+1) 390 ILOC=2*NPARM**2+4*NUMGR*NPARM+9*NUMGR+21*NPARM+10 RETURN C C 40 WCOEF(NUMGR) 400 ILOC=2*NPARM**2+4*NUMGR*NPARM+9*NUMGR+22*NPARM+11 RETURN C C 41 WCOEF1(NUMGR) (DOES NOT APPEAR IN THE PROGRAM BY NAME) 410 ILOC=2*NPARM**2+4*NUMGR*NPARM+10*NUMGR+22*NPARM+11 RETURN C C 42 WPT(NPARM) 420 ILOC=2*NPARM**2+4*NUMGR*NPARM+11*NUMGR+22*NPARM+11 RETURN C C 43 WVEC(NPARM) 430 ILOC=2*NPARM**2+4*NUMGR*NPARM+11*NUMGR+23*NPARM+11 RETURN C C 44 X(NPARM+1) 440 ILOC=2*NPARM**2+4*NUMGR*NPARM+11*NUMGR+24*NPARM+11 RETURN C C 45 XKEEP(NPARM+1) 450 ILOC=2*NPARM**2+4*NUMGR*NPARM+11*NUMGR+25*NPARM+12 RETURN C C 46 XRK(NPARM+1) 460 ILOC=2*NPARM**2+4*NUMGR*NPARM+11*NUMGR+26*NPARM+13 RETURN C C 47 Y(NUMGR+2*NPARM) (WITH V, OPPOSITE AA, B, BETA, D, PICOR, C ZWORK; FOLLOWS V) 470 ILOC=2*NPARM**2+4*NUMGR*NPARM+8*NUMGR+16*NPARM+10 RETURN C C 48 ZWORK(NPARM) (OPPOSITE V, Y; FOLLOWS PICOR) 480 ILOC=2*NPARM**2+3*NUMGR*NPARM+6*NUMGR+18*NPARM+13 RETURN END SUBROUTINE DERST(IOPTN,NPARM,NUMGR,PTTBL,IPTB,INDM,PARAM, *IPT,PARAM1,V,KCNTYP,CONFUN) C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION PTTBL(IPTB,INDM),PARAM(NPARM),PARAM1(NPARM), *V(NUMGR+2*NPARM+1,NPARM+2),KCNTYP(NUMGR), *CONFUN(NUMGR,NPARM+1) C C THIS SUBROUTINE USES FNSET TO COMPUTE CONFUN(I,1) AND THE PARTIAL C DERIVATIVES OF THE FUNCTION WHOSE VALUE IS IN CONFUN(I,1) FOR C CERTAIN VALUE(S) OF I. NOTE THAT WE DO NOT WANT THE ICNTYP COMPUTED C BY FNSET TO OVERRIDE THE ICNTYP (OR JCNTYP) CARRIED INTO THIS C SUBROUTINE IN ICNTYP, SO WE USE KCNTYP WHEN WE CALL FNSET. (THE C ICNTYP COMPUTED BY FNSET WAS STORED EARLIER THROUGH A CALL TO ERCMP1 C FROM CONMAX.) C C IF THE ONES DIGIT OF IOPTN IS 0, WE CALL FNSET WITH INDFN=1 TO DO THE C COMPUTATIONS DIRECTLY USING FORMULAS SUPPLIED BY THE USER. IOPONE=IOPTN-(IOPTN/10)*10 IF(IOPONE)100,100,200 100 CALL FNSET (NPARM,NUMGR,PTTBL,IPTB,INDM,PARAM,IPT,1,KCNTYP, *CONFUN) RETURN C C HERE THE ONES DIGIT OF IOPTN IS 1, AND WE APPROXIMATE THE PARTIAL C DERIVATIVES USING CENTERED DIFFERENCE APPROXIMATIONS. C 200 IOPTTH=(IOPTN-(IOPTN/100000)*100000)/10000 C C SET PRECISION DEPENDENT CONSTANTS. SPCMN=D1MACH(3) DELT=SQRT(SPCMN) DELT2=DELT+DELT IF(IOPTTH)300,300,700 C C HERE IOPONE=1 AND IOPTTH=0, AND WE WORK ONLY WITH CONSTRAINT IPT, C WHERE IPT WILL BE AN INTEGER BETWEEN 1 AND NUMGR. C L WILL BE THE INDEX OF THE VARIABLE WITH RESPECT TO WHICH WE ARE C COMPUTING THE PARTIAL DERIVATIVE. 300 DO 500 L=1,NPARM C C SET PARAM1 EQUAL TO PARAM, ECXEPT WITH ITS LTH COMPONENT INCREASED C BY DELT. DO 400 J=1,NPARM PARAM1(J)=PARAM(J) 400 CONTINUE PARAM1(L)=PARAM(L)+DELT C C NOW CALL FNSET WITH INDFN=0 TO PLACE THE FUNCTION IN CONSTRAINT C IPT EVALUATED AT POINT PARAM1 IN CONFUN(IPT,1). CALL FNSET(NPARM,NUMGR,PTTBL,IPTB,INDM,PARAM1,IPT,0, * KCNTYP,CONFUN) UP=CONFUN(IPT,1) C C SET PARAM1 EQUAL TO PARAM, ECXEPT WITH ITS LTH COMOPONENT DECREASED C BY DELT, AND CALL FNSET AGAIN. PARAM1(L)=PARAM(L)-DELT CALL FNSET(NPARM,NUMGR,PTTBL,IPTB,INDM,PARAM1,IPT,0, * KCNTYP,CONFUN) C C NOW WE CAN COMPUTE THE CENTERED-DIFFERENCE APPROXIMATION TO THE PARTIAL C DERIVATIVE OF THE FUNCTION IN CONSTRAINT IPT WITH RESPECT TO THE LTH C VARIABLE AT THE POINT PARAM. THIS BELONGS IN CONFUN(IPT,L+1), AND C WE COULD PUT IT THERE NOW IF THE USER FOLLOWED DIRECTIONS AND DID NOT C CHANGE CONFUN(IPT,L+1) (SINCE INDFN=0) IN LATER FNSET CALLS, BUT TO C BE SAFE WE TEMPORARILY STORE IT IN V(L,1). C NOTE THAT V IS USED ELSEWHERE IN THE PROGRAM, BUT HERE IT IS JUST A C WORK ARRAY, WHILE THE WORK ARRAY PARAM1 IS NOT USED ELSEWHERE IN C THE PROGRAM. V(L,1)=(UP-CONFUN(IPT,1))/DELT2 500 CONTINUE C C NOW COMPUTE THE VALUE OF THE FUNCTION AT PARAM, AND THEN PUT THE C EARLIER-COMPUTED PARTIAL DERIVATIVES INTO CONFUN. CALL FNSET(NPARM,NUMGR,PTTBL,IPTB,INDM,PARAM,IPT,0, *KCNTYP,CONFUN) DO 600 L=1,NPARM CONFUN(IPT,L+1)=V(L,1) 600 CONTINUE RETURN C C HERE IOPONE=1 AND IOPTTH=1, AND EACH TIME FNSET IS CALLED IT WILL C COMPUTE VALUES FOR THE FUNCTIONS IN THE LEFT SIDES OF ALL CONSTRAINTS C (EXCEPT THOSE WHERE FNSET SETS ICNTYP(I)=0) IF IPT=0, AND WILL COMPUTE C VALUES FOR THE FUNCTIONS IN THE LEFT SIDES OF ALL STANDARD (I.E. TYPE C -1 OR -2) CONSTRAINTS IF IPT=-1. C WE FIRST SAVE IPT IN CASE THE USER CHANGES IT IN A FNSET CALL; WE WILL C RESTORE IT AFTER EACH FNSET CALL. 700 IPTKP=IPT NPAR1=NPARM+1 C C WE WILL COMPUTE APPROXIMATIONS TO PARTIAL DERIVATIVES FOR THOSE C CONSTRAINTS WHICH FNSET IS ASKED BY IPT TO COMPUTE. TO DETERMINE WHICH C THESE ARE WE ZERO OUT KCNTYP; AFTER A FNSET CALL, THE DESIRED C CONSTRAINTS WILL BE THE CONSTRAINTS K WITH KCNTYP(K) .NE. 0 IF IPT=0, C OR THE CONSTRAINTS K WITH KCNTYP(K) .LT. 0 IF IPT=-1. DO 800 K=1,NUMGR KCNTYP(K)=0 800 CONTINUE C C NOW FOLLOW BASICALLY THE SAME PROCEDURES AS IN THE IOPTTH=0 CASE DONE C ABOVE. DO 1800 L=1,NPARM DO 900 J=1,NPARM PARAM1(J)=PARAM(J) 900 CONTINUE PARAM1(L)=PARAM(L)+DELT CALL FNSET(NPARM,NUMGR,PTTBL,IPTB,INDM,PARAM1,IPT,0, * KCNTYP,CONFUN) IPT=IPTKP DO 1300 K=1,NUMGR IF(IPT)1100,1000,1000 1000 IF(KCNTYP(K))1200,1300,1200 1100 IF(KCNTYP(K))1200,1300,1300 C C SAVE THE UPPER NUMBERS IN COLUMN NPARM+1 OF V. 1200 V(K,NPAR1)=CONFUN(K,1) 1300 CONTINUE C C REVISE PARAM1 AND CALL FNSET AGAIN. PARAM1(L)=PARAM(L)-DELT CALL FNSET(NPARM,NUMGR,PTTBL,IPTB,INDM,PARAM1,IPT,0, * KCNTYP,CONFUN) IPT=IPTKP DO 1700 K=1,NUMGR IF(IPT)1500,1400,1400 1400 IF(KCNTYP(K))1600,1700,1600 1500 IF(KCNTYP(K))1600,1700,1700 C C STORE THE APPROXIMATE PARTIAL DERIVATIVES WITH RESPECT TO THE LTH C VARIABLE IN THE LTH COLUMN OF V. 1600 V(K,L)=(V(K,NPAR1)-CONFUN(K,1))/DELT2 1700 CONTINUE 1800 CONTINUE C CALL FNSET AGAIN TO COMPUTE THE VALUES OF THE FUNCTIONS AT POINT C PARAM, AND THEN PUT THE EARLIER-COMPUTED PARTIAL DERIVATIVES INTO C CONFUN. CALL FNSET(NPARM,NUMGR,PTTBL,IPTB,INDM,PARAM,IPT,0, *KCNTYP,CONFUN) DO 2300 K=1,NUMGR IF(IPT)2000,1900,1900 1900 IF(KCNTYP(K))2100,2300,2100 2000 IF(KCNTYP(K))2100,2300,2300 2100 DO 2200 L=1,NPARM CONFUN(K,L+1)=V(K,L) 2200 CONTINUE 2300 CONTINUE RETURN END SUBROUTINE SLPCON(IOPTN,NPARM,NUMGR,FUN,IFUN,PTTBL,IPTB, *INDM,TOLCON,RCHIN,IRK,ITYPM1,ITYPM2,ICNTYP,RCHDWN,NUMLIM,ITERSL, *PRJSLP,FUNTBL,IYRCT,X,MACT1,IACT1,JCNTYP,IPHSE,ENCHG,IWORK, *LIWRK,WORK,LWRK,PARSER,ISUCC,PARAM,ERROR) C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION FUN(IFUN),PTTBL(IPTB,INDM),ICNTYP(NUMGR), *FUNTBL(NUMGR,NPARM+1),IYRCT(NUMGR+2*NPARM),X(NPARM+1), *IACT1(NUMGR),PARAM(NPARM),ERROR(NUMGR+3),JCNTYP(NUMGR), *PARSER(NPARM),IWORK(LIWRK),WORK(LWRK) C C SET MACHINE AND PRECISION DEPENDENT CONSTANTS. C NWRIT=I1MACH(2) ONE=1.0D0 ZERO=ONE-ONE TWO=ONE+ONE FOUR=TWO+TWO TEN=FOUR+FOUR+TWO SPCMN=D1MACH(3) BIG=ONE/SPCMN TOL1=TEN*TEN*SPCMN TOL2=TEN*SPCMN ILC05=ILOC(5,NPARM,NUMGR) ILC07=ILOC(7,NPARM,NUMGR) ILC08=ILOC(8,NPARM,NUMGR) ILC11=ILOC(11,NPARM,NUMGR) ILC13=ILOC(13,NPARM,NUMGR) ILC18=ILOC(18,NPARM,NUMGR) ILC19=ILOC(19,NPARM,NUMGR) ILC25=ILOC(25,NPARM,NUMGR) ILC35=ILOC(35,NPARM,NUMGR) ILC45=ILOC(45,NPARM,NUMGR) ILC47=ILOC(47,NPARM,NUMGR) NUMIN=0 ISUCC=0 ENORM=ERROR(NUMGR+1) NPAR1=NPARM+1 NG3=NUMGR+3 C IF ITERSL=0, SET IYRCT(1)=-1 FOR USE IN SETU1 AND TO TELL SLNPRO NOT C TO TRY TO USE INFORMATION FROM A PREVIOUS VERTEX. IF(ITERSL)50,50,300 50 IYRCT(1)=-1 C C CALL BNDSET TO SET (OR RESET) THE COEFFICIENT CHANGE BOUNDS. 300 CALL BNDSET(NPARM,X,ITERSL,NUMIN,PRJSLP,WORK(ILC07),WORK(ILC45), *WORK(ILC05)) C C CALL SETU1 TO SET UP FOR SLNPRO AND, IF NUMIN=0, TO DETERMINE C WHICH CONSTRAINTS ARE ACTIVE AND STORE FUNCTION AND GRADIENT VALUES C FOR THEM IN FUNTBL. 400 CALL SETU1(IOPTN,NUMGR,NPARM,NUMIN,RCHIN,PTTBL,IPTB,INDM, *FUN,IFUN,FUNTBL,WORK(ILC07),PARAM,ICNTYP,RCHDWN,ERROR,MACT1, *IACT1,BNDLGT,IYRCT,IPHSE,IWORK,LIWRK,WORK,LWRK,WORK(ILC08), *IWORK(ILC13),WORK(ILC35),M) C C SET UNIT (FOR USE IN RCHMOD) EQUAL TO THE VALUE OF BNDLGT AFTER C SETU1 IS CALLED WITH NUMIN=0. IF(NUMIN)500,500,1000 500 UNIT=BNDLGT C C CALL SLNPRO TO COMPUTE A SEARCH DIRECTION X. 1000 CALL SLNPRO(WORK(ILC35),M,NPAR1,IYRCT,WORK(ILC47), *IWORK(ILC18),IWORK(ILC19),NPARM,NUMGR,X,INDIC) C C IF INDIC .GT. 0 THEN SLNPRO FAILED TO PRODUCE AN X, AND IF WE HAVE C REACHED THE SLPCON ITERATION LIMIT WE RETURN WITH THE WARNING C ISUCC=1. IF(INDIC)1300,1300,1800 C C HERE SLNPRO SUCCEEDED AND WE SET PRJSLP=1.0 INITIALLY FOR SEARSL. 1300 PRJSLP=ONE C C WE NOW WISH TO DETERMINE PRJLIM = THE SMALLER OF 1.0/SPCMN AND C THE LARGEST VALUE OF PRJSLP FOR WHICH THE LINEAR STANDARD CONSTRAINTS C ARE SATISFIED FOR THE PARAMETER VECTOR PARAM+PRJSLP*X. THIS C WILL GIVE AN UPPER BOUND FOR LINE SEARCHING. NOTE THAT IN C THEORY WE SHOULD HAVE PRJLIM .GE. 1.0 SINCE THE LINEAR STANDARD C CONSTRAINTS SHOULD BE SATISFIED FOR PRJSLP=0.0 AND PRJSLP=1.0, BUT C ROUNDOFF ERROR COULD AFFECT THIS A LITTLE. IF THERE ARE NO C LINEAR STANDARD CONSTRAINTS, WE SET PRJLIM=1.0/SPCMN. 1400 PRJLIM=BIG C*****INSERT TO MAKE SEARCHING LESS VIOLENT. C PRJLIM=TWO C*****END INSERT IF(ITYPM1)1430,1430,1405 1405 DO 1425 I=1,NUMGR IF(ICNTYP(I)+1)1425,1407,1425 C WE WISH TO HAVE SUMMATION (FUNTBL(I,J+1)*(PARAM(J)+PRJSLP*X(J))) C + C(I) .LE. 0.0 FOR I=1,...,NUMGR, ICNTYP(I) = -1, C WHERE THE ITH CONSTRAINT APPLIED TO PARAM SAYS C SUMMATION (FUNTBL(I,J+1)*PARAM(J)) + C(I) .LE. 0.0, SO C(I) IS THE C CONSTANT TERM ON THE LEFT SIDE OF LINEAR CONSTRANT I. C THUS FOR I=1,...,NUMGR, ICNTYP(I) = -1, WE WANT PRJLIM*SS .LE. SSS, C WHERE SS = SUMMATION (FUNTBL(I,J+1)*X(J)) AND SSS = -C(I) - C SUMMATION (FUNTBL(I,J+1)*PARAM(J)) = -FUNTBL(I,1). 1407 SS=ZERO DO 1410 J=1,NPARM SS=SS+FUNTBL(I,J+1)*X(J) 1410 CONTINUE C IF SS .LT. 10.0*SPCMN THIS CONSTRAINT WILL NOT PUT A SIGNIFICANT C RESTRICTION ON PRJSLP. IF(SS-TOL2)1425,1415,1415 C HERE SS .GE. 10.0*SPCMN AND WE COMPARE SSS/SS AGIANST PRJLIM. 1415 QUOTS=-FUNTBL(I,1)/SS IF(PRJLIM-QUOTS)1425,1425,1420 1420 PRJLIM=QUOTS 1425 CONTINUE C DO NOT ALLOW A PRJSLP SMALLER THAN TOL1. 1430 IF(PRJSLP-TOL1)1440,1470,1470 1440 PRJSLP=TOL1 C CALL SEARSL TO DO A LINE SEARCH IN DIRECTION X. 1470 CALL SEARSL(IOPTN,NUMGR,NPARM,PRJLIM,TOL1,X,FUN,IFUN,PTTBL, *IPTB,INDM,PARAM,ERROR,RCHDWN,MACT1,IACT1,IPHSE,UNIT, *TOLCON,RCHIN,ITYPM1,ITYPM2,IWORK,LIWRK,WORK,LWRK,WORK(ILC11), *WORK(ILC25),PRJSLP,EMIN,EMIN1,PARSER,NSRCH) C C COMPUTE THE ERROR NORM CHANGE ENCHG. ENCHG=EMIN-ENORM C C IF WE HAVE AN IMPROVEMENT IN THE ERROR NORM ENORM OF MORE THAN TOL1 C WE UPDATE PARAM AND ERROR AND RETURN WITH ISUCC=0, INDICATING SUCCESS. C OTHERWISE WE CHECK TO SEE IF WE HAVE REACHED THE SLPCON ITERATION C LIMIT, AND IF SO WE RETURN WITH ISUCC=1, INDICATING FAILURE. IF(ENCHG+TOL1)1600,1800,1800 C C HERE WE HAD AN IMPROVEMENT IN THE ERROR NORM ENORM OF MORE THAN TOL1. 1600 DO 1700 J=1,NPARM PARAM(J)=PARSER(J) 1700 CONTINUE CALL ERCMP1(IOPTN,NPARM,NUMGR,FUN,IFUN,PTTBL,IPTB,INDM, *PARAM,1,IPHSE,IWORK,LIWRK,WORK(ILC08),ICNTYP,IPMAX, *ISMAX,ERROR) RETURN C C HERE WE DID NOT OBTAIN AN IMPROVED ERROR NORM SO WE RETURN WITH THE C WARNING ISUCC=1 IF WE HAVE DONE NUMLIN ITERATIONS IN SLPCON. 1800 IF(NUMIN-NUMLIM)2000,1900,1900 1900 ISUCC=1 RETURN C C HERE WE DID NOT OBTAIN AN IMPROVED ERROR NORM BUT WE HAVE NOT YET DONE C NUMLIM ITERATIONS IN SLPCON SO WE INCREMENT NUMIN, SET IYRCT(1)=-1 TO C TELL SLNPRO NOT TO TRY TO USE INFORMATION FROM THE PREVIOUS FAILED C VERTEX, AND GO BACK TO CALL BNDSET AND TRY ANOTHER ITERATION WITH C A DIFFERENT TRUST REGION. 2000 NUMIN=NUMIN+1 IYRCT(1)=-1 GO TO 300 END SUBROUTINE BNDSET(NPARM,X,ITERSL,NUMIN,PRJSLP,COFBND,XKEEP, *BNDKP) C C THIS SUBROUTINE SETS THE BOUNDS ON THE COEFFICIENT CHANGES IN C SLNPRO. C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION X(NPARM+1),COFBND(NPARM),XKEEP(NPARM+1),BNDKP(NPARM) C C SET MACHINE AND PRECISION DEPENDENT CONSTANTS FOR BNDSET. ONE=1.0D0 TWO=ONE+ONE FOUR=TWO+TWO TEN=FOUR+FOUR+TWO C NWRIT=I1MACH(2) SPCMN=D1MACH(3) C C SET INITIAL PARAMETERS. FACT1, FACT3A, FACT3B, CHLM1, AND CHLM2 C SHOULD BE BETWEEN 0.0 AND 1.0, WHILE FACT2 SHOULD BE .GT. 1.0. FACT1=(ONE+TWO)/FOUR FACT2=TWO FACT3A=ONE/TEN FACT3B=ONE/(TEN*TEN) FACT4=TWO/TEN CHLM1=ONE/TEN CHLM2=(FOUR+FOUR)/TEN TSTPRJ=ONE/TWO-ONE/(TEN*TEN*TEN) EPSIL=TEN*TEN*SPCMN EPSIL1=(ONE+ONE/(TEN*TEN*TEN))*EPSIL C BND IS THE INITIAL BOUND ON ALL COEFFICIENT CHANGES. BND=TWO/(TEN*TEN) C END OF SETTING MACHINE AND PRECISION DEPENDENT CONSTANTS FOR BNDSET. C IF(NUMIN-1)100,2000,2100 100 IF(ITERSL-1)200,400,600 C C HERE NUMIN=0 AND ITERSL=0, SO WE ARE IN THE FIRST BNDSET CALL SINCE THE C LAST RK SUCCESS (IF ANY), SO WE SET INITIAL BOUNDS. 200 DO 300 J=1,NPARM COFBND(J)=BND 300 CONTINUE RETURN C C HERE NUMIN=0 AND ITERSL=1, SO THE LAST BNDSET CALL RESULTED IN C THE FIRST SUCCESSFUL PRINCIPAL ERROR NORM IMPROVEMENT, C AND SO WE SAVE COFBND IN BNDKP AND X IN XKEEP. WE WILL NOT C CHANGE COFBND HERE. 400 DO 500 J=1,NPARM XKEEP(J)=X(J) BNDKP(J)=COFBND(J) 500 CONTINUE RETURN C C HERE NUMIN=0 AND ITERSL .GE. 2, SO WE HAVE HAD AT LEAST 2 SUCCESSES, C WITH THE COEFFICIENTS AND BOUNDS FOR THE LAST ONE IN X AND C COFBND RESPECTIVELY, AND THE COEFFICIENTS AND BOUNDS FOR THE C PREVIOUS ONE IN XKEEP AND BNDKP RESPECTIVELY. WE WILL FORM A C NEW COFBND, AND SHIFT THE OLD COFBND INTO BNDKP AND X INTO XKEEP. 600 DO 1900 J=1,NPARM C SAVE THE OLD COFBND(J) IN BSAVE. BSAVE=COFBND(J) C IF AT BOTH THE LAST AND PREVIOUS SUCCESSFUL ITERATION THE CHANGES C IN A COEFFICIENT RELATIVE TO ITS BOUND WERE .GE. CHLM2 IN ABSOLUTE C VALUE AND IN THE SAME DIRECTION, WE LOOSEN THE BOUND BY A FACTOR C OF FACT2. IF THE RELATIVE CHANGES WERE .GE. CHLM1 IN ABSOLUTE C VALUE AND IN OPPOSITE DIRECTIONS, WE TIGHTEN THE BOUND BY A FACTOR C OF FACT1 BECAUSE OF SUSPECTED OSCILLATION. WE ALSO TIGHTEN THE C BOUND IF BOTH RELATIVE CHANGES WERE LESS THAN CHLM1 IN ABSOLUTE C VALUE IN ORDER TO PREVENT A LONG SEQUENCE OF OSCILLATIONS OF THE C SAME SMALL ORDER. OTHERWISE WE LEAVE THE BOUND ALONE. C THE NEXT FOUR IF STATEMENTS CHECK TO SEE IF THE BOUND SHOULD BE C LOOSENED. IF(X(J)-CHLM2*COFBND(J))800,700,700 700 IF(XKEEP(J)-CHLM2*BNDKP(J))1100,1000,1000 800 IF(X(J)+CHLM2*COFBND(J))900,900,1100 900 IF(XKEEP(J)+CHLM2*BNDKP(J))1000,1000,1100 C LOOSEN THE BOUND. 1000 COFBND(J)=FACT2*COFBND(J) GO TO 1800 C C HERE THE BOUND SHOULD NOT BE LOOSENED. THE NEXT FIVE IF C STATEMTENTS CHECK TO SEE IF IT SHOULD BE TIGHTENED. 1100 IF(X(J)-CHLM1*COFBND(J))1300,1200,1200 1200 IF(XKEEP(J)+CHLM1*BNDKP(J))1600,1600,1800 1300 IF(X(J)+CHLM1*COFBND(J))1400,1400,1500 1400 IF(XKEEP(J)-CHLM1*BNDKP(J))1800,1600,1600 C HERE WE HAVE ABS(X(J)) .LT. CHLM1*COFBND(J). 1500 IF(ABS(XKEEP(J))-CHLM1*BNDKP(J))1600,1800,1800 C TIGHTEN THE BOUND. 1600 COFBND(J)=FACT1*COFBND(J) C DO NOT ALLOW THE BOUND TO DROP BELOW EPSIL. IF(COFBND(J)-EPSIL)1700,1800,1800 1700 COFBND(J)=EPSIL C C SAVE X(J) AND THE OLD COFBND(J). 1800 BNDKP(J)=BSAVE XKEEP(J)=X(J) 1900 CONTINUE C C IF THE LAST PROJECTION FACTOR IS SMALLER THAN .499, WE TIGHTEN THE C BOUNDS BY A FACTOR OF 0.2, WITH THE RESTRICTION THAT WE DO NOT C ALLOW THE BOUNDS TO DROP BELOW EPSIL. IF(PRJSLP-TSTPRJ)1920,1980,1980 1920 DO 1960 J=1,NPARM COFBND(J)=FACT4*COFBND(J) IF(COFBND(J)-EPSIL)1940,1960,1960 1940 COFBND(J)=EPSIL 1960 CONTINUE 1980 RETURN C C HERE NUMIN=1 SO THE LAST BNDSET CALL RESULTED IN FAILURE TO C IMPROVE THE PRINCIPAL ERROR NORM, AND WE SET FACT3= C FACT3A AND TIGHTEN THE BOUNDS. 2000 FACT3=FACT3A GO TO 2200 C C HERE NUMIN .GT. 1 SO THERE HAVE BEEN AT LEAST 2 SUCCESSIVE C FAILURES, AND WE SET FACT3=FACT3B AND TIGHTEN THE BOUNDS. 2100 FACT3=FACT3B C C TIGHTEN THE BOUNDS BY A FACTOR OF FACT3. 2200 ITIGHT=1 DO 2700 J=1,NPARM BSAVE=COFBND(J) COFBND(J)=FACT3*BSAVE C WE DO NOT ALLOW A BOUND TO DROP BELOW EPSIL. IF(COFBND(J)-EPSIL)2300,2600,2600 C IF THE BOUND WAS ALREADY (ESSENTIALLY) AT EPSIL, KEEP TRACK OF C THIS BY NOT SETTING ITIGHT=0. 2300 IF(BSAVE-EPSIL1)2500,2500,2400 2400 ITIGHT=0 2500 COFBND(J)=EPSIL GO TO 2700 2600 ITIGHT=0 2700 CONTINUE C IF ALL THE BOUNDS WERE ALREADY (ESSENTIALLY) AT EPSIL, WE TRY C RESETTING THE BOUNDS TO THEIR ORIGINAL VALUES. IF(ITIGHT)1980,1980,2800 2800 CONTINUE C2800 WRITE(NWRIT,2900) C2900 FORMAT(/52H *****RESETTING BOUNDS TO THEIR ORIGINAL VALUES*****) GO TO 200 END SUBROUTINE SETU1(IOPTN,NUMGR,NPARM,NUMIN,RCHIN,PTTBL,IPTB, *INDM,FUN,IFUN,FUNTBL,COFBND,PARAM,ICNTYP,RCHDWN,ERROR,MACT1, *IACT1,BNDLGT,IYRCT,IPHSE,IWORK,LIWRK,WORK,LWRK,CONFUN,IACT,V,M) C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION PTTBL(IPTB,INDM),FUN(IFUN),FUNTBL(NUMGR,NPARM+1), *COFBND(NPARM),PARAM(NPARM),ERROR(NUMGR+3), *V(NUMGR+2*NPARM+1,NPARM+2),IACT(NUMGR),IACT1(NUMGR), *IYRCT(NUMGR+2*NPARM),ICNTYP(NUMGR),CONFUN(NUMGR,NPARM+1), *IWORK(LIWRK),WORK(LWRK) C C THIS SUBROUTINE SETS UP V FOR SLNPRO TO SOLVE A MODIFIED LINEARIZED C (ABOUT THE OLD PARAMETERS IN PARAM) VERSION OF OUR PROBLEM. C C SET MACHINE AND PRECISION CONSTANTS FOR SETU1. C NWRIT=I1MACH(2) ONE=1.0D0 ZERO=ONE-ONE TWO=ONE+ONE FOUR=TWO+TWO TEN=FOUR+FOUR+TWO C END OF SETTING MACHINE AND PRECISION DEPENDENT CONSTANTS FOR SETU1. C ILC22=ILOC(22,NPARM,NUMGR) ILC24=ILOC(24,NPARM,NUMGR) NPAR1=NPARM+1 NPAR2=NPARM+2 IOPTTH=(IOPTN-(IOPTN/100000)*100000)/10000 C C THE LINEARIZED PROBLEM REPLACES THE APPROXIMATING FUNCTION BY ITS C FIRST ORDER TAYLOR SERIES, SO FUN(I)-(APPROXIMATING FUNCTION)(I) IS C REPLACED BY ERROR(I)-(SUMMATION OF COEFFICIENT CHANGES TIMES PARTIAL C DERIVATIVES OF THE APPROXIMATING FUNCTION WITH RESPECT TO THOSE C COEFFICIENTS) IF ICNTYP(I)=2, AND IF ICNTYP(I)=1 WE REPLACE THE LEFT C SIDE OF CONSTRAINT I BY ERROR(I)+(SUMMATION OF COEFFICIENT CHANGES TIMES C PARTIAL DERIVATIVES OF THE LEFT SIDE OF CONSTRAINT I). C V AND M ARE THE OUTPUT QUANTITIES. M WILL KEEP TRACK OF THE NUMBER C OF CONSTRAINTS IN THE LP PROBLEM TO BE SOLVED BY SLNPRO. M=0 ENORM=ERROR(NUMGR+1) STFUDG=ONE/TEN C C COMPUTE THE LENGTH OF THE LONGEST X VECTOR SATISFYING THE COEFFICIENT C CHANGE BOUNDS. SUM=ZERO DO 20 J=1,NPARM SUM=SUM+(COFBND(J))**2 20 CONTINUE BNDLGT=SQRT(SUM) BNDFUD=STFUDG*BNDLGT C C WE WILL SAY A PRIMARY CONSTRAINT IS ACTIVE IF ERROR(I) (OR ABS(ERROR(I C IF ICNTYP(I)=2) .GE. ENORM-RCHDWN*BNDLGT. ACTLIM=ENORM-RCHDWN*BNDLGT C C WE WILL SAY A TYPE -2 CONSTRAINT IS ACTIVE IF ERROR(I) .GE. -RCHIND. RCHIND=RCHIN*BNDLGT C IF(NUMIN)80,80,40 C HERE NUMIN IS NOT 0, AND WE WILL KEEP THE OLD ACTIVE CONSTRAINT SET C AND FOREGO RECOMPUTING FUNCTION VALUES AND GRADIENTS. 40 MACT=MACT1 M=MACT DO 60 L=1,MACT IACT(L)=IACT1(L) 60 CONTINUE GO TO 440 C C HERE NUMIN=0, SO WE WILL FIRST COMPUTE A NEW SET OF ACTIVE INDICES, C THEN PUT THE FUNCTION VALUES AND GRADIENTS FOR THESE INDICES IN C FUNTBL, WHERE THEY WILL REMAIN THROUGHOUT THIS CALL TO SLPCON. 80 DO 240 I=1,NUMGR IF(ICNTYP(I))220,240,100 100 IF(ICNTYP(I)-1)120,120,160 C C HERE ICNTYP(I)=1 AND WE WILL DECLARE THE CONSTRAINT TO BE +ACTIVE IF AND C ONLY IF ERROR(I) .GE. ACTLIM. 120 IF(ERROR(I)-ACTLIM)240,140,140 C C DECLARE CONSTRAINT I TO BE (+)ACTIVE. 140 M=M+1 IACT(M)=I GO TO 240 C C HERE ICNTYP(I)=2 AND WE WILL DECLARE THE CONSTRAINT TO BE +ACTIVE IF AND C ONLY IF ERROR(I) .GE. ACTLIM OR -ACTIVE IF AND ONLY IF ERROR(I) .LE. C -ACTLIM. 160 IF(ERROR(I)-ACTLIM)180,140,140 180 IF(ERROR(I)+ACTLIM)200,200,240 C C DECLARE CONSTRAINT I TO BE -ACTIVE. 200 M=M+1 IACT(M)=-I GO TO 240 C C HERE ICNTYP(I) .LT. 0 AND WE WILL DECLARE THE CONSTRAINT TO BE ACTIVE IF C AND ONLY IF ICNTYP(I)=-1, OR ICNTYP(I)=-2 AND ERROR(I) .GE. -RCHIND. 220 IF(ICNTYP(I)+1)230,140,140 230 IF(ERROR(I)+RCHIND)240,140,140 240 CONTINUE MACT=M C C NOW PUT ACTIVE VALUES AND GRADIENTS IN FUNTBL. IF(IOPTTH)260,260,380 C HERE IOPTTH=0 AND WE CALL DERST FOR EACH ACTIVE CONSTRAINT. 260 DO 360 L=1,MACT I=IABS(IACT(L)) IPT=I C CALL DERST TO COMPUTE BOTH FUNCTION AND GRADIENT VALUES. CALL DERST(IOPTN,NPARM,NUMGR,PTTBL,IPTB,INDM,PARAM,IPT, * WORK(ILC24),V,IWORK(ILC22),CONFUN) C COPY THE VALUES FOR CONSTRAINT I INTO FUNTBL. DO 340 J=1,NPAR1 FUNTBL(I,J)=CONFUN(I,J) 340 CONTINUE 360 CONTINUE GO TO 440 C C HERE IOPTTH=1 AND ONLY ONE DERST CALL IS NEEDED. C IF IPHSE .LT. 0 OR NO ICNTYP(L) IS POSITIVE, SET IPT=-1 TO TELL DERST C TO COMPUTE STANDARD CONSTRAINTS ONLY, WHILE OTHERWISE SET IPT=0 TO C TELL DERST TO COMPUTE ALL CONSTRAINTS. 380 IF(IPHSE)389,383,383 383 DO 386 L=1,NUMGR IF(ICNTYP(L))386,386,392 386 CONTINUE 389 IPT=-1 GO TO 395 392 IPT=0 395 CALL DERST(IOPTN,NPARM,NUMGR,PTTBL,IPTB,INDM,PARAM,IPT, *WORK(ILC24),V,IWORK(ILC22),CONFUN) C COPY THE ACTIVE FUNCTION AND GRADIENT VALUES INTO FUNTBL. DO 420 L=1,MACT I=IABS(IACT(L)) DO 400 J=1,NPAR1 FUNTBL(I,J)=CONFUN(I,J) 400 CONTINUE 420 CONTINUE C C NOW SET UP THE ACTIVE CONSTRAINTS IN V FOR SLNPRO. 440 DO 680 L=1,MACT I=IABS(IACT(L)) IF(ICNTYP(I))610,680,460 460 IF(ICNTYP(I)-1)480,480,520 C C HERE ICNTYP(I)=1 AND WE SET UP A CONSTRAINT OF THE FORM C GRADIENT.CHANGE - W .LE. -CONSTRAINT VALUE. 480 DO 500 J=1,NPARM V(L,J)=FUNTBL(I,J+1) 500 CONTINUE V(L,NPAR1)=-ONE V(L,NPAR2)=-ERROR(I) GO TO 680 520 IF(IACT(L))580,580,540 C C HERE ICNTYP(I)=2 AND IACT(L) .GT. 0, AND WE SET UP A CONSTRAINT OF THE C FORM -GRADIENT.CHANGE - W .LE. -(FUN - CONSTRAINT VALUE). 540 DO 560 J=1,NPARM V(L,J)=-FUNTBL(I,J+1) 560 CONTINUE V(L,NPAR1)=-ONE V(L,NPAR2)=-ERROR(I) GO TO 680 C C HERE ICNTYP(I)=2 AND IACT(L) .LT. 0, AND WE SET UP A CONSTRAINT OF THE C FORM GRADIENT.CHANGE - W .LE. FUN - CONSTRAINT VALUE. 580 DO 600 J=1,NPARM V(L,J)=FUNTBL(I,J+1) 600 CONTINUE V(L,NPAR1)=-ONE V(L,NPAR2)=ERROR(I) GO TO 680 C 610 IF(ICNTYP(I)+1)630,620,620 C C HERE ICNTYP(I)=-1 AND WE SET UP A CONSTRAINT OF THE FORM C GRADIENT.CHANGE .LE. -CONSTRAINT VALUE. 620 DO 625 J=1,NPARM V(L,J)=FUNTBL(I,J+1) 625 CONTINUE V(L,NPAR1)=ZERO V(L,NPAR2)=-ERROR(I) GO TO 680 C C HERE ICNTYP(I)=-2 AND WE FIRST COMPUTE THE LENGTH OF THE GRADIENT. 630 SUM=ZERO DO 640 J=1,NPARM SUM=SUM+(FUNTBL(I,J+1))**2 640 CONTINUE GRDLGT=SQRT(SUM) C NOW SET UP A CONSTRAINT OF THE FORM GRADIENT.CHANGE .LE. C -MIN(1.0,CONSTRAINT VALUE)*BNDFUD*GRDLGT, SO IF GRDLGT .GT. 0.0 WE C HAVE (-GRADIENT/GRDLGT).(CHANGE/BNDLGT) .GE. STFUDG*MIN(1.0, C CONSTRAINT VALUE). DO 660 J=1,NPARM V(L,J)=FUNTBL(I,J+1) 660 CONTINUE V(L,NPAR1)=ZERO RT=ERROR(I) IF(RT-ONE)675,675,665 665 RT=ONE 675 V(L,NPAR2)=-RT*BNDFUD*GRDLGT 680 CONTINUE C C SET THE CONSTRAINTS OF THE FORM -X(J) .LE. COFBND(J) AND C X(J) .LE. COFBND(J). DO 800 J=1,NPARM M=M+2 MM1=M-1 DO 700 K=1,NPAR1 V(MM1,K)=ZERO V(M,K)=ZERO 700 CONTINUE V(MM1,J)=-ONE V(M,J)=ONE V(MM1,NPAR2)=COFBND(J) V(M,NPAR2)=COFBND(J) 800 CONTINUE C C NOW SET THE BOTTOM ROW. TO MINIMIZE W = X(NPARM+1) WE MAXIMIZE -W. MP1=M+1 DO 900 J=1,NPAR2 V(MP1,J)=ZERO 900 CONTINUE V(MP1,NPAR1)=ONE C C THIS SECTION ADJUSTS IYRCT TO EITHER TELL SLNPRO TO DO THE INITIAL C EXCHANGES STRICTLY ACCORDING TO A PIVOTING STRATEGY (BY SETTING C IYRCT(1)=-1) OR TO SPECIFY AN INITIAL VERTEX FOR SLNPRO, NAMELY THE C VERTEX CORRESPONDING TO THE LAST LINEAR PROGRAMMING SOLUTION. C IF IYRCT(1) IS -1 ALREADY WE DO NOT ATTEMPT TO SPECIFY A VERTEX, BUT C WE STORE MACT IN MACT1 AND IACT IN IACT1 FOR POSSIBLE LATER USE. IF(IYRCT(1))1700,1100,1100 C HERE IYRCT(1) .NE. -1, AND WE CONSIDER THE PRESENT ENTRIES IN IYRCT C ONE BY ONE. 1100 DO 1600 J=1,NPAR1 JJ=IYRCT(J) IF(JJ-MACT1)1200,1200,1500 C HERE ENTRY J OF IYRCT CORRESPONDS TO A FORMER ACTIVE CONSTRAINT AT C SOME POINT IABS(KK), WHERE THE SIGN OF KK WILL INDICATE WHETHER THE C CONSTRAINT WAS +ACTIVE OR -ACTIVE. 1200 KK=IACT1(JJ) C WE NOW CHECK TO SEE IF THIS FORMER ACTIVE CONSTRAINT IS STILL C ACTIVE WITH THE SAME SIGN. IF SO, WE RESET IYRCT(J) TO THE PRESENT C NUMBER OF THIS CONSTRAINT, AND IF NOT (WHICH WILL OCCUR IFF THE K C LOOP BELOW IS COMPLETED), WE WILL NOT TRY TO DETERMINE A VERTEX, SO C WE WILL SET IYRCT(1)=-1 AND LEAVE THE J LOOP. DO 1400 K=1,MACT IF(KK-IACT(K))1400,1300,1400 1300 IYRCT(J)=K GO TO 1600 1400 CONTINUE IYRCT(1)=-1 GO TO 1700 C HERE ENTRY J OF IYRCT CORRESPONDS TO A CONSTRAINT BEYOND THE ACTIVE C POINT CONSTRAINTS, AND WE ADJUST IYRCT(J) BY THE DIFFERENCE OF THE C PRESENT AND FORMER NUMBER OF ACTIVE CONSTRAINTS. 1500 IYRCT(J)=IYRCT(J)+MACT-MACT1 1600 CONTINUE C WE HAVE NOW FILLED IN IYRCT(1),...,IYRCT(NPARM+1) WITH DISTINCT C POSITIVE INTEGERS BETWEEN 1 AND M, AND WE FILL IN THE REST OF IYRCT C SO THAT IYRCT WILL CONTAIN A PERMUTATION OF 1,...,M. TO BE CONSISTENT C WITH SLNPRO WE PUT IYRCT(NPARM+2),...,IYRCT(M) IN DECREASING ORDER. L=NPAR1 DO 1660 I=1,M II=M-I+1 C SKIP II IF IT IS ALREADY IN IYRCT. DO 1640 J=1,NPAR1 IF(II-IYRCT(J))1640,1660,1640 1640 CONTINUE L=L+1 IYRCT(L)=II 1660 CONTINUE C C SAVE MACT IN MACT1 AND IACT IN IACT1 AND RETURN. 1700 MACT1=MACT DO 1800 J=1,MACT IACT1(J)=IACT(J) 1800 CONTINUE RETURN END SUBROUTINE SLNPRO(V,M,N,IYRCT,Y,IXRCT,IYCCT,NPARM,NUMGR,X, *INDIC) C***BEGIN PROLOGUE SLNPRO C***ROUTINES CALLED SJELIM C***PURPOSE THIS SUBROUTINE SOLVES THE LINEAR PROGRAMMING PROBLEM C MAXIMIZE Z = -V(M+1,1)*X(1)-...-V(M+1,N)*X(N) C WHERE X(1),...,X(N) ARE FREE VARIABLES, SUBJECT TO C V(I,1)*X(1)+...+V(I,N)*X(N) .LE. V(I,N+1), FOR I=1,..,M, C WHERE M .GE. N. C (INFORMATION CONCERNING TOLERANCES AND BASIC VARIABLES C IS ALSO TRANSMITTED USING M, N, AND IYRCT.) C***REFERENCES AVDEYEVA, L. I. AND ZUKHOVITSKIY, S. I., C LINEAR AND CONVEX PROGRAMMING, C SAUNDERS, PHILADELPHIA, 1966. C***END PROLOGUE SLNPRO C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION V(NUMGR+2*NPARM+1,NPARM+2),IYRCT(NUMGR+2*NPARM), *X(NPARM+1),Y(NUMGR+2*NPARM),IXRCT(NUMGR+2*NPARM),IYCCT(NPARM+1) C C GIVEN INTEGERS M AND N (WITH M .GE. N) AND A MATRIX V, C THIS SUBROUTINE SOLVES THE LINEAR PROGRAMMING PROBLEM C MAXIMIZE Z=-V(M+1,1)X(1)-...-V(M+1,N)X(N)+V(M+1,N+1) C SUBJECT TO THE CONSTRAINTS C V(I,1)X(1)+...+V(I,N)X(N) .LE. V(I,N+1), I=1,...,M C USING ESSENTIALLY THE METHOD IN THE BOOK BY AVDEYEVA AND C ZUKHOVITSKIY. Y(I)=-V(I,1)X(1)-...-V(I,N)X(N)+V(I,N+1), C I=1,...,M ARE SLACK VARIABLES. THE METHOD HAS 4 PHASES. C C FIRST, XS ARE EXCHANGED FOR YS TO GET A PROBLEM C INVOLVING ONLY NONNEGATIVE VARIABLES, THE YS BEING C SELECTED IN THE ORDER DETERMINED BY IYRCT AND A PIVOTING C STRATEGY. AT THE BEGINNING OF THIS ROUTINE WE MUST HAVE C IYRCT(1) NONPOSITIVE, OR IYRCT MUST CONTAIN SOME C PERMUTATION OF THE INTEGERS 1,...,M (SEE BELOW). C SECOND, THE SLACK CONSTANTS OF THE DUAL PROBLEM ARE MADE C (SIGNIFICANTLY) NONNEGATIVE. C THIRD, THE COST COEFFICIENTS OF THE DUAL PROBLEM ARE MADE C (SIGNIFICANTLY) NONNEGATIVE. C FINALLY, THE SOLUTION VECTOR IS COMPUTED. C C THE VARIABLE INDIC WILL BE GIVEN VALUE C 0, IF A SOLUTION WAS FOUND NORMALLY C 1, IF THERE WAS TROUBLE IN PHASE 1 C 2, IF THERE WAS TROUBLE IN PHASE 2 (EITHER ROUND OFF C ERROR, OR THE ORIGINAL PROBLEM WAS INCONSISTENT OR C UNBOUNDED) C 3, IF THERE WAS TROUBLE IN PHASE 3 (EITHER ROUND OFF C ERROR, OR THE ORIGINAL CONSTRAINTS WERE INCONSISTENT) C 4, IF LIMJOR MODIFIED JORDAN ELIMINATIONS WERE USED IN C PHASES 2 AND 3 C -1, IF A SOLUTION WAS FOUND BUT IN ORDER TO OVERCOME C TROUBLE IN PHASE 2 OR 3 IT WAS NECESSARY TO TEMPORARILY C RELAX THE RESTRICTION ON DIVISION BY NUMBERS WITH SMALL C ABSOLUTE VALUE. THUS THERE IS AN INCREASED CHANCE OF C SERIOUS ROUNDOFF ERROR IN THE RESULTS. C -2, IF A SOLUTION WAS FOUND NORMALLY, EXCEPT THAT C THE PARAMETERS REA AND REA1 WERE INCREASED BY A SIGNAL C FROM THE CALLING PROGRAM (NAMELY, M=-M). THE INCREASED C TOLERANCES MAY HAVE ALLOWED MORE ERROR THAN USUAL. C -3, IF IN ORDER TO COMPLETE PHASE 1 IT WAS NECESSARY TO C TEMPORARILY RELAX THE RESTRICTION ON DIVISION BY NUMBERS C WITH SMALL ABSOLUTE VALUE. THUS THERE IS AN INCREASED C CHANCE OF SERIOUS ROUNDOFF ERROR IN THE RESULTS. C -4, IF A SOLUTION WAS FOUND NORMALLY, EXCEPT THAT REA AND REA1 C WERE DECREASED BY A SIGNAL FROM THE CALLING PROGRAM (NAMELY C N=-N) IN ORDER TO TRY FOR MORE ACCURACY. THIS INCREASES THE C CHANCES OF SERIOUS ROUNDOFF ERROR IN THE RESULTS. C NOTE THAT INDIC=-3 WILL OVERWRITE (AND THUS CONCEAL) INDIC=-2 C OR INDIC=-4, AND INDIC=-1 WILL OVERWRITE INDIC=-2, -3, OR -4 C C SET MACHINE DEPENDENT PARAMETERS FOR SUBROUTINE SLNPRO. C SET SPCMN=B**(-ITT), WHERE B IS THE BASE AND ITT IS THE NUMBER C OF BASE B DIGITS IN THE MANTISSA. SPCMN IS THE MINIMUM C RELATIVE SPACING ABS((X1-X2)/X2) BETWEEN TWO SUCCESSIVE C FLOATING POINT NUMBERS, SO IT IS THE SPACING BETWEEN TWO C SUCCESSIVE FLOATING POINT NUMBERS IN THE CLOSED INTERVAL C (0.1,1.0). WE ALSO HAVE SPCMN=10.0**(-ITT*(LOG10)(B))= C 10.0**(-TNMAN), WHERE TNMAN IS THE BASE 10 EQUIVALENT OF C THE LENGTH OF THE MANTISSA. C C***FIRST EXECUTABLE STATEMENT SLNPRO SPCMN=D1MACH(3) C SET PRECISION DEPENDENT CONSTANTS FOR SLNPRO. ONE=1.0D0 ZERO=ONE-ONE TWO=ONE+ONE FOUR=TWO+TWO TEN=FOUR+FOUR+TWO C SET REA (ROUND OFF ERROR ADJUSTMENT) = C MAX(10.0**(-8),100.0*SPCMN). THUS REA=10.0**(-EXREA), C WHERE EXREA=MIN(8,TNMAN-2). C DIVISION BY NUMBERS .LE. REA IN ABSOLUTE VALUE WILL NOT BE C PERMITTED. REA=TEN*TEN*SPCMN IF(REA-TEN**(-8))10000,10010,10010 10000 REA=TEN**(-8) C SET REA1=10.0*SPCMN. THUS REA1=10.0**(-(TNMAN-1)). C NUMBERS IN ROW M+1 OR COLUMN N+1 WHICH ARE .LE. REA1 IN C ABSOLUTE VALUE WILL BE TREATED AS ZEROES. SLNPRO ASSUMES C THAT 0.0 .LT. REA1 .LE. REA. 10010 REA1=TEN*SPCMN C END OF INITIAL SETTING OF MACHINE DEPENDENT PARAMETERS FOR C SLNPRO. THESE PARAMETERS MAY BE ADJUSTED BY A COMMAND FROM C THE CALLING PROGRAM. C INDIC=0 LIMJOR=300 C M BEING NEGATIVE IS A SIGNAL TO INCREASE REA AND REA1, C THUS TREATING MORE NUMBERS WITH SMALL ABSOLUTE VALUES AS C ZEROES. THIS MAY GIVE THIS ROUTINE A BETTER CHANCE TO C SUCCEED, BUT MAY ALSO CAUSE LARGER ERRORS. IF(M)1001,10001,10001 C RESET M. 1001 M=-M REA=SQRT(REA) REA1=SQRT(REA1) INDIC=-2 C N BEING NEGATIVE IS A SIGNAL TO DECREASE REA AND REA1 TO TRY C FOR MORE ACCURACY. AMONG OTHER THINGS, THIS MAKES IT MORE C LIKELY THAT THE PREVIOUS VERTEX WILL BE RETAINED IN PHASE 1 C BELOW, BUT IT ALSO COULD INCREASE ROUND OFF ERROR. 10001 IF(N)10002,1002,1002 C RESET N. 10002 N=-N REA=REA1 REA1=REA1/(TEN*TEN) INDIC=-4 C PRESERVE REA IN CASE IT MUST BE TEMPORARILY RELAXED. C IRLAX=0 INDICATES REA IS NOT RELAXED AT THIS STAGE. 1002 REAKP=REA IRLAX=0 C IN COLUMN N+1, NUMBERS .LE. REA2 IN ABSOLUTE VALUE WILL BE C TREATED AS ZEROES. REA2=REA1 NP1=N+1 MP1=M+1 KTJOR=0 IBACK=0 C SET V(MP1,NP1)=0.0 SO THE DESCRIPTIONS IN AND FOLLOWING THE C PROLOGUE WILL AGREE. V(MP1,NP1)=ZERO C THE ONLY REASON FOR THE FOLLOWING THREE STATEMENTS IS TO C AVOID THE ERROR MESSAGE ON SOME MACHINES THAT THESE C VARIABLES HAVE NOT BEEN ASSIGNED A VALUE. THEY WILL BE C REASSIGNED A VALUE BEFORE THE PROGRAM REACHES A SPOT WHERE C THEY WILL ACTUALLY BE USED. DIST=ONE AMPRV=ONE AMPR2=ONE C SET IXRCT. IXRCT(I)=0 MEANS SOME Y IS IN ROW I, WHILE C IXRCT(I)=K.NE.0 MEANS X(K) IS IN ROW I. DO 1 I=1,M IXRCT(I)=0 1 CONTINUE C C EXCHANGE THE XS AT THE TOP OF THE TABLE FOR YS. C IF IYRCT(1) IS NONPOSITIVE, WE SET IYRCT AND CHOOSE THE C LARGEST POSSIBLE RESOLVENTS FOR THE EXCHANGES. IF C IYRCT(1) IS POSITIVE, IYRCT WILL HAVE BEEN PREVIOUSLY SET C AND WE TRY TO EXCHANGE IN ROWS IYRCT(1),...,IYRCT(N), C STILL EMPLOYING A PIVOTING STRATEGY, BUT IF WE CANNOT, WE C EXCHANGE IN ROWS IYRCT(N+1),...,IYRCT(M). IF(IYRCT(1))1003,1003,1005 1003 I10=1 I20=M C IF WE HAVE NO INFORMATION FROM A PREVIOUS VERTEX, WE GIVE C UP A LITTLE ACCURACY IN COLUMN N+1 TO HAVE A BETTER CHANCE C OF SUCCESS. REA2=REA C THIS ROUTINE HAS A BACKTRACKING OPTION WHICH SOMETIMES C INCREASES ACCURACY BUT SOMETIMES LEADS TO FAILURE DUE TO C CYCLING. IT IS SUGGESTED THAT THIS OPTION BE EMPLOYED IF C INFORMATION ABOUT A STARTING VERTEX IS AVAILABLE, AND C OTHERWISE BE DISABLED BY SETTING IBACK=1. IBACK=1 DO 1004 I=1,M IYRCT(I)=I 1004 CONTINUE GO TO 1006 1005 I10=1 I20=N 1006 J=0 C SET THE LOWER BOUND ON THE ABSOLUTE VALUE OF A RESOLVENT IN C PHASE 1. ALSO SET IFAIL=0 TO INDICATE THE RESOLVENT SEARCH C IN THIS COLUMN HAS NOT FAILED. REA3=REA IFAIL=0 2 J=J+1 IF(J-N)1007,1007,9 C SET I1, I2 ACCORDING TO THE STRATEGY WE ARE USING. 1007 I1=I10 I2=I20 AMAX=ZERO C SEARCH FOR A RESOLVENT IN ROWS IYRCT(I1),...,IYRCT(I2). 10003 DO 1012 I=I1,I2 IYTMP=IYRCT(I) IF(IXRCT(IYTMP))1012,1009,1012 1009 ABSV=ABS(V(IYTMP,J)) IF(ABSV-AMAX)1012,1012,1011 1011 IYRI=IYTMP AMAX=ABSV 1012 CONTINUE C CHECK TO SEE IF THE PROSPECTIVE RESOLVENT IS LARGE ENOUGH C IN ABSOLUTE VALUE. IF(AMAX-REA3)1013,1013,7 C EXCHANGE X(J) FOR Y(IYRI). 7 CALL SJELIM(MP1,1,NP1,IYRI,J,NPARM,NUMGR,V) IXRCT(IYRI)=J IYCCT(J)=IYRI C IYCCT(J)=IYRI MEANS Y(IYRI) IS IN COLUMN J. C RESET REA3 AND IFAIL SINCE WE SUCCESSFULLY FOUND A RESOLVENT IN C THIS COLUMN, AND THE RESOLVENT SEARCH IN THE NEXT COLUMN HAS C NOT FAILED. REA3=REA IFAIL=0 GO TO 2 C WE HAVE NOT FOUND A SUITABLE RESOLVENT IN ROWS IYRCT(I1), C ...IYRCT(I2). IF I2 .LT. M WE SEARCH THE REST OF COLUMN J. 1013 IF(I2-M)1014,10004,10004 1014 I1=I2+1 I2=M GO TO 10003 C HERE WE FAILED TO FIND A RESOLVENT IN COLUMN J WITH ABSOLUTE C VALUE .GT. REA3. IF IFAIL=0 WE SET INDIC=-3 AND TRY AGAIN C WITH REA3 REDUCED. IF THIS HAS ALREADY BEEN TRIED WE SET C INDIC=1 AND RETURN. 10004 IF(IFAIL)10005,10005,8 10005 IFAIL=1 INDIC=-3 REA3=REA1 GO TO 1007 C 8 INDIC=1 RETURN C C REARRANGE THE ROWS OF V SO THAT X(1),...,X(N) COME FIRST C IN THAT ORDER. REDEFINE IYRCT SO THAT AFTER THE C REARRANGEMENT IS DONE, IYRCT(I)=K WILL MEAN Y(K) IS IN C ROW I (FOR I GREATER THAN N). 9 DO 10 I=1,M IYRCT(I)=I 10 CONTINUE IROW=0 11 IROW=IROW+1 IF(IROW-M)12,12,20 12 IF(IXRCT(IROW))13,11,13 13 IF(IXRCT(IROW)-IROW)14,11,14 C NOW X(L) IS IN ROW IROW, BUT WE WANT IT IN ROW L. 14 L=IXRCT(IROW) LL=IXRCT(L) IF(LL)15,16,15 C X(L) IS IN ROW IROW, WHILE X(LL) IS IN ROW L. 15 IXRCT(IROW)=LL IXRCT(L)=L GO TO 17 C X(L) IS IN ROW IROW, WHILE Y(IYRCT(L)) IS IN ROW L. 16 IXRCT(IROW)=0 IYRCT(IROW)=IYRCT(L) IXRCT(L)=L C NOW EXCHANGE THE CONTENTS OF ROWS IROW AND L. 17 DO 18 J=1,NP1 TEMP=V(IROW,J) V(IROW,J)=V(L,J) V(L,J)=TEMP 18 CONTINUE IF(IXRCT(IROW))19,11,19 19 IF(IXRCT(IROW)-IROW)14,11,14 C NOW IXRCT IS NO LONGER NEEDED, SO STORE THE PRESENT IYCCT C IN IT. 20 DO 21 I=1,N IXRCT(I)=IYCCT(I) 21 CONTINUE C END OF PHASE 1. C C THE FIRST N ROWS OF V GIVE THE XS IN TERMS OF CERTAIN C YS. THESE ROWS WILL NOT BE CHANGED BY LATER OPERATIONS. C C WE NOW ATTACK THE MAXIMIZATION PROBLEM BY LOOKING AT THE C DUAL PROBLEM OF MINIMIZING A FORM GIVEN BY THE C COEFFICIENTS IN V(N+1,N+1) THROUGH V(M,N+1) WITH SLACK C TERMS IN THE BOTTOM ROW OF V. C SEARCH ROW M+1 FOR A SIGNIFICANTLY NEGATIVE ELEMENT. IF C THERE ARE NONE, PROCEED TO THE ACTUAL MINIMIZATION C PROBLEM. STICK WITH COLUMN JJ UNTIL V(M+1,JJ) .GE. -REA1. JJ=0 22 JJ=JJ+1 IF(JJ-N)1015,1015,1035 1015 IF(V(MP1,JJ)+REA1)24,22,22 C C WE HAVE V(M+1,JJ) SIGNIFICANTLY NEGATIVE. SEARCH COLUMN C JJ FOR A POSITIVE ELEMENT, TREATING A VERY SMALL V(I,J) C AS A ZERO. IF THERE ARE NO POSITIVE ELEMENTS THE DUAL C CONSTRAINTS WERE INCONSISTENT, SO THE ORIGINAL PROBLEM WAS C INCONSISTENT OR UNBOUNDED. 24 I=N INAMP=0 25 I=I+1 IF(I-M)1016,1016,1020 1016 IF(V(I,JJ)-REA)25,25,1017 C C NOW V(I,JJ) .GT. REA. WE SEARCH ROW I FOR INDICES K SUCH C THAT V(M+1,K) .GE. 0.0.OR.K .LT. JJ, AND V(I,K) .LT. -REA, AND C FIND THE MAXIMUM RATIO (I.E. THE RATIO WITH SMALLEST C ABSOLUTE VALUE, IF V(M+1,K) .GE. 0.0) V(M+1,K)/V(I,K). IF C THERE IS NO SUCH K WE LOOK AT POSITIVE V(I,K) BELOW. 1017 INDST=0 DO 32 J=1,N IF(V(MP1,J))1018,28,28 1018 IF(J-JJ)28,32,32 28 IF(V(I,J)+REA)29,32,32 29 DIST1=V(MP1,J)/V(I,J) IF(INDST)31,31,30 30 IF(DIST1-DIST)32,32,31 31 DIST=DIST1 INDST=1 K=J 32 CONTINUE IF(INDST)35,35,33 C C WE NOW COMPUTE V(I,JJ)*DIST AND GO ON TO LOOK AT OTHER C ROWS TO MINIMIZE THIS QUANTITY (I.E. TO MAXIMIZE ITS C ABSOLUTE VALUE, IF V(M+1,K) .GE. 0.0). THIS IS THE NEGATIVE C OF THE CHANGE IN V(M+1,JJ). 33 BMPRV=V(I,JJ)*DIST IF(INAMP)34,34,1019 1019 IF(BMPRV-AMPRV)34,25,25 34 AMPRV=BMPRV INAMP=1 KPMP1=I KPMP2=K C (KPMP1,KPMP2) GIVES THE POSITION OF THE BEST PROSPECTIVE C RESOLVENT FOUND SO FAR. GO TO 25 C C IF THERE WAS NO INDEX K SUCH THAT V(M+1,K) .GE. 0.0.OR.K .LT. C JJ, AND V(I,K) .LT. -REA, WE LOOK FOR THE SMALLEST (I.E. C LARGEST IN ABSOLUTE VALUE) RATIO V(M+1,K)/V(I,K) FOR C V(I,K) .GT. REA AND V(M+1,K) .LT. 0.0, AND PERFORM AN C ELIMINATION WITH RESOLVENT V(I,K). THERE IS AT LEAST ONE C SUCH K, NAMELY JJ. C THIS WILL FINISH PHASE 2 UNLESS BACKTRACKING IS NECESSARY. 35 DIST=ONE DO 39 J=1,N IF(V(MP1,J))36,39,39 36 IF(V(I,J)-REA)39,39,37 37 DIST1=V(MP1,J)/V(I,J) IF(DIST1-DIST)38,39,39 38 DIST=DIST1 K=J 39 CONTINUE GO TO 49 C 1020 IF(INAMP)1021,1021,1023 C AT THIS POINT INAMP IS POSITIVE IFF THERE WAS AT LEAST ONE C ELEMENT .GT. REA IN COLUMN JJ. IF THERE WERE NONE, WE C TEMPORARILY RELAX REA AND TRY AGAIN. 1021 IF(IRLAX)1022,1022,41 1022 IRLAX=1 INDIC=-1 REA=REA1 GO TO 24 C 41 INDIC=2 RETURN C C CHECK TO SEE IF V(MP1,KPMP2) IS VERY SMALL IN ABSOLUTE C VALUE OR NEGATIVE. THIS INDICATES DEGENERACY. 1023 IF(V(MP1,KPMP2)-REA)1024,1024,43 C DO AN ELIMINATION WITH RESOLVENT V(KPMP1,KPMP2). 43 I=KPMP1 K=KPMP2 GO TO 49 C C WE ARE NOW STUCK IN DEGENERATE COLUMN KPMP2. WE SEARCH C EACH DEGENERATE COLUMN IN WHICH WE ARE STUCK FOR A C RESOLVENT WHICH WILL KEEP US FROM GETTING STUCK IN THIS C COLUMN NEXT TIME, AND TO REDUCE THE ROUND-OFF ERROR WE C TAKE THE SMALLEST OF THESE (I.E. LARGEST IN ABSOLUTE C VALUE) AS OUR ACTUAL RESOLVENT. 1024 AMIN=ONE DO 1034 J=1,N C COLUMN J MAY BE DEGENERATE IF 0.0 .LE. V(M+1,J) .LE. REA, C OR V(M+1,J) .LT. 0.0.AND.J .LT. JJ. IF(V(MP1,J))1025,1026,1026 1025 IF(J-JJ)1027,1034,1034 1026 IF(V(MP1,J)-REA)1027,1027,1034 C WE WILL BE STUCK IN COLUMN J IFF THERE IS AN INDEX ID FOR C WHICH V(ID,JJ) .GT. REA AND V(ID,J) .LT. -REA. IF THIS IS THE C CASE, CHOOSING SUCH AN ID SO THAT V(ID,JJ)/V(ID,J) IS C MINIMIZED (I.E. MAXIMIZED IN ABSOLUTE VALUE) AND TAKING C V(ID,J) AS THE RESOLVENT WILL INSURE THAT WE DONT GET C STUCK IN COLUMN J NEXT TIME. 1027 DIST=ONE DO 1031 I=NP1,M IF(V(I,JJ)-REA)1031,1031,1028 1028 IF(V(I,J)+REA)1029,1031,1031 1029 DIST1=V(I,JJ)/V(I,J) IF(DIST1-DIST)1030,1031,1031 1030 DIST=DIST1 ID=I 1031 CONTINUE IF(DIST-ONE/TWO)1032,1034,1034 C WE HAVE NOW DETERMINED THAT WE ARE STUCK IN COLUMN J. C IF V(ID,J) .LT. AMIN THEN V(ID,J) IS THE BEST RESOLVENT C FOUND SO FAR. 1032 IF(V(ID,J)-AMIN)1033,1034,1034 1033 AMIN=V(ID,J) KPMP1=ID KPMP2=J 1034 CONTINUE C THE BEST RESOLVENT IS V(KPMP1,KPMP2), SO WE DO AN C ELIMINATION. GO TO 43 C 49 KTJOR=KTJOR+1 IF(KTJOR-LIMJOR)50,50,73 50 CALL SJELIM(MP1,NP1,NP1,I,K,NPARM,NUMGR,V) ITEMP=IYRCT(I) IYRCT(I)=IYCCT(K) IYCCT(K)=ITEMP C RESET REA AND IRLAX. REA=REAKP IRLAX=0 C IF NOW V(M+1,JJ) HAS BEEN MADE NOT SIGNIFICANTLY NEGATIVE, C WE GO TO THE NEXT COLUMN. OTHERWISE WE TRY AGAIN IN C COLUMN JJ. IF(V(MP1,JJ)+REA1)24,22,22 C C IN THE UNLIKELY EVENT THAT SOME V(M+1,J) IS STILL VERY C SIGNIFICANTLY NEGATIVE WE BACKTRACK TO COLUMN J. THIS C COULD NOT HAPPEN IF THERE WERE NO ROUNDOFF ERROR AND WE C COULD ALLOW DIVISION BY NUMBERS WITH VERY SMALL ABSOLUTE C VALUE. OMIT BACKTRACKING IF IBACK=1. 1035 IF(IBACK)1036,1036,51 1036 DO 1038 J=1,N IF(V(MP1,J)+REA)1037,1037,1038 1037 JJ=J GO TO 24 1038 CONTINUE C END OF PHASE 2. C 51 I=N KKK=0 C C SEARCH FOR A SIGNIFICANTLY NEGATIVE ELEMENT BETWEEN C V(N+1,N+1) AND V(N+1,M). IF THERE ARE NONE WE HAVE THE C MINIMAL POINT OF THE DUAL PROBLEM (AND THUS THE MAXIMAL C POINT OF THE DIRECT PROBLEM) ALREADY. 52 I=I+1 IF(I-M)1039,1039,1043 1039 IF(V(I,NP1)+REA2)1040,52,52 C C SEARCH FOR A NEGATIVE ELEMENT IN ROW I, TREATING A NUMBER C WHICH IS VERY SMALL IN ABSOLUTE VALUE AS A ZERO. IF THERE C ARE NO NEGATIVE ELEMENTS THE DUAL PROBLEM WAS UNBOUNDED C BELOW, SO THE ORIGINAL CONSTRAINTS WERE INCONSISTENT. C FIND THE INDEX K OF THE LARGEST (I.E. SMALLEST IN ABSOLUTE C VALUE, IF V(M+1,K) .GE. 0.0) RATIO V(M+1,K)/V(I,K) WITH C V(I,K) .LT. -REA. 1040 INDST=0 DO 58 J=1,N IF(V(I,J)+REA)55,58,58 55 DIST1=V(MP1,J)/V(I,J) IF(INDST)57,57,56 56 IF(DIST1-DIST)58,58,57 57 K=J INDST=1 DIST=DIST1 58 CONTINUE IF(INDST)1041,1041,60 C RELAX REA AND LOOK FOR NEGATIVE ELEMENTS WITH SMALLER C ABSOLUTE VALUE. 1041 IF(IRLAX)1042,1042,59 1042 IRLAX=1 INDIC=-1 REA=REA1 GO TO 1040 C 59 INDIC=3 RETURN C C COMPUTE THE IMPROVEMENT DIST*V(I,N+1) IN THE VALUE OF THE C FORM USING V(I,K) AS THE RESOLVENT. SET KKK=1 TO INDICATE C A SIGNIFICANTLY NEGATIVE V(I,N+1) WAS FOUND, AND LOOK AT C THE OTHER ROWS TO FIND THE RESOLVENT GIVING THE LARGEST C IMPROVEMENT. 60 BMPR2=DIST*V(I,NP1) C RESET IRLAX SO THAT THE NEXT ROW WHICH NEEDS RELAXING DOES C NOT TERMINATE THE ROUTINE. REA WILL REMAIN RELAXED UNTIL C AFTER THE NEXT ELIMINATION. IRLAX=0 IF(KKK)62,62,61 61 IF(BMPR2-AMPR2)52,52,62 62 KKK=1 KEEP=I KEEP1=K AMPR2=BMPR2 GO TO 52 C KKK=0 HERE IFF NONE OF THE COST COEFFICIENTS ARE C SIGNIFICANTLY NEGATIVE. 1043 IF(KKK)1048,1044,1048 C CHECK TO SEE IF ANY OF THE NUMBERS IN THE BOTTOM ROW HAVE C BECOME VERY SIGNIFICANTLY NEGATIVE. IF SO, WE MUST C BACKTRACK TO PHASE 2 (SEE COMMENT ABOVE STATEMENT 1035). C OMIT BACKTRACKING IF IBACK=1. 1044 IF(IBACK)1045,1045,74 1045 DO 1047 J=1,N IF(V(MP1,J)+REA)1046,1046,1047 1046 JJ=J GO TO 24 1047 CONTINUE GO TO 74 C CHECK TO SEE IF V(MP1,KEEP1) IS VERY SMALL IN ABSOLUTE C VALUE OR NEGATIVE. THIS INDICATES DEGENERACY. 1048 IF(V(MP1,KEEP1)-REA)1049,1049,65 C DO AN ELIMINATION WITH RESOLVENT V(KEEP,KEEP1). 65 I=KEEP K=KEEP1 GO TO 71 C C WE ARE NOW STUCK IN DEGENERATE COLUMN KEEP1. WE SEARCH C EACH DEGENERATE COLUMN IN WHICH WE ARE STUCK FOR A C RESOLVENT WHICH WILL KEEP US FROM GETTING STUCK IN THIS C COLUMN NEXT TIME. IF WE ARE NOT USING THE OPTION C DESCRIBED IN THE COMMENTS PRECEDING STATEMENT 1055, WE C TAKE THE SMALLEST OF THESE (I.E. THE LARGEST IN ABSOLUTE C VALUE) AS OUR ACTUAL RESOLVENT IN ORDER TO REDUCE THE C GROWTH OF ROUND-OFF ERROR. 1049 AMIN=ONE MXRKN=NP1 DO 1072 J=1,N C COLUMN J MAY BE DEGENERATE IF V(M+1,J) .LE. REA. IF(V(MP1,J)-REA)1050,1050,1072 C WE WILL BE STUCK IN COLUMN J IFF THERE IS AN INDEX ID FOR C WHICH V(ID,N+1) .LT. -REA2 AND V(ID,J) .LT. -REA. IF THIS C IS THE CASE, CHOOSING SUCH AN ID SO THAT V(ID,N+1)/V(ID,J) C IS MAXIMIZED AND TAKING V(ID,J) AS THE RESOLVENT WILL C INSURE THAT WE DONT GET STUCK IN COLUMN J NEXT TIME. 1050 DIST=-ONE DO 1054 I=NP1,M IF(V(I,NP1)+REA2)1051,1054,1054 1051 IF(V(I,J)+REA)1052,1054,1054 1052 DIST1=V(I,NP1)/V(I,J) IF(DIST1-DIST)1054,1054,1053 1053 DIST=DIST1 ID=I 1054 CONTINUE IF(DIST+ONE/TWO)1072,1072,1055 C C WE HAVE NOW DETERMINED THAT WE ARE STUCK IN COLUMN J. C THE FOLLOWING STATEMENTS ATTEMPT TO BREAK DEGENERACY C FASTER BY LOOKING ONE ITERATION INTO THE FUTURE, THAT IS, C BY CHOOSING FROM THE PROSPECTIVE RESOLVENTS FOUND ABOVE C THAT ONE WHICH MINIMIZES THE MINIMUM NUMBER OF STICKING C PLACES IN ANY ROW AT THE NEXT STAGE. C BECAUSE OF COMPUTER TIME AND THE POSSIBLE LOSS OF ACCURACY C DUE TO LESSENED PIVOTING (EVEN THOUGH TIES ARE ALWAYS C BROKEN IN FAVOR OF THE RESOLVENT WITH GREATEST ABSOLUTE C VALUE), IT IS SUGGESTED THAT THIS OPTION BE OMITTED IF C INFORMATION WAS AVAILABLE FROM A PREVIOUS VERTEX. THIS C WILL BE THE CASE IFF THE BACKTRACKING OPTION WAS USED, C THAT IS, IFF IBACK=0. 1055 IF(IBACK)1070,1070,1056 C COMPUTE WHAT THE NEW BOTTOM ROW WOULD BE (EXCEPT FOR C POSITION J) IF V(ID,J) WERE USED AS THE RESOLVENT, AND C PUT THE RESULTS INTO Y. 1056 ROWQ=V(MP1,J)/V(ID,J) DO 1058 L=1,N IF(L-J)1057,1058,1057 1057 Y(L)=V(MP1,L)-V(ID,L)*ROWQ 1058 CONTINUE LRKNT=-1 C WE LOOK FOR A ROW WHICH WILL HAVE A SIGNIFICANTLY NEGATIVE C LAST ELEMENT BUT A MINIMUM NUMBER OF PLACES WHERE WE WILL C BE STUCK IN DEGENERATE COLUMNS. LRKNT=-1 MEANS WE HAVE C NOT YET FOUND A ROW WHICH WILL HAVE A SIGNIFICANTLY C NEGATIVE LAST ELEMENT. DO 1068 II=NP1,M IF(II-ID)1059,1068,1059 1059 ROWQ=V(II,J)/V(ID,J) RTCOL=V(II,NP1)-V(ID,NP1)*ROWQ IF(RTCOL+REA2)1060,1068,1068 C IF WE HAVE ALREADY LOCATED A RESOLVENT WHICH WILL FINISH C THE ROUTINE, BUT THE PRESENT PROSPECTIVE RESOLVENT WOULD C GIVE A ROW WITH A SIGNIFICANTLY NEGATIVE LAST ELEMENT, WE C LOOK AT THE NEXT PROSPECTIVE RESOLVENT FOR PIVOTING C PURPOSES. 1060 IF(MXRKN+1)1061,1072,1061 1061 LRKNT=0 C NOW COUNT THE NUMBER (LRKNT) OF STICKING PLACES IN ROW II C AT THE NEXT ITERATION. DO 1065 JJ=1,N IF(JJ-J)1062,1065,1062 1062 IF(Y(JJ)-REA)1063,1063,1065 1063 IF(V(II,JJ)-V(ID,JJ)*ROWQ+REA)1064,1065,1065 1064 LRKNT=LRKNT+1 IF(LRKNT-MXRKN)1065,1065,1068 1065 CONTINUE IF(LRKNT-MXRKN)1067,1066,1068 1066 IF(V(ID,J)-AMIN)1067,1068,1068 1067 MXRKN=LRKNT AMIN=V(ID,J) KEEP=ID KEEP1=J 1068 CONTINUE C LRKNT=-1 HERE WOULD MEAN THIS RESOLVENT WOULD FINISH THE C ROUTINE. IF LRKNT .GE. 0 THEN MXRKN .GE. 0 ALSO, SO WE WILL C NOT HAVE EARLIER FOUND A RESOLVENT WHICH WILL FINISH THE C ROUTINE. IF(LRKNT+1)1072,1069,1072 1069 IF(MXRKN+1)1071,1070,1071 1070 IF(V(ID,J)-AMIN)1071,1072,1072 1071 MXRKN=-1 AMIN=V(ID,J) KEEP=ID KEEP1=J 1072 CONTINUE C THE BEST RESOLVENT IS V(KEEP,KEEP1), SO WE DO AN C ELIMINATION. GO TO 65 C 71 KTJOR=KTJOR+1 IF(KTJOR-LIMJOR)72,72,73 72 CALL SJELIM(MP1,NP1,NP1,I,K,NPARM,NUMGR,V) ITEMP=IYRCT(I) IYRCT(I)=IYCCT(K) IYCCT(K)=ITEMP C RESET REA AND IRLAX. REA=REAKP IRLAX=0 GO TO 51 C 73 INDIC=4 RETURN C END OF PHASE 3. WE NOW HAVE THE VERTEX WE ARE SEEKING. C C READ OFF THE Y VALUES FOR THIS VERTEX. 74 DO 75 J=1,N IYCJ=IYCCT(J) Y(IYCJ)=ZERO 75 CONTINUE DO 76 I=NP1,M IYRI=IYRCT(I) Y(IYRI)=V(I,NP1) 76 CONTINUE C COMPUTE THE XS FROM THE YS. RECALL THAT IXRCT CONTAINS C THE FORMER IYCCT. DO 78 I=1,N X(I)=V(I,NP1) DO 77 J=1,N IXRJ=IXRCT(J) X(I)=X(I)-V(I,J)*Y(IXRJ) 77 CONTINUE 78 CONTINUE C C NOW PUT THE VALUES IN IYCCT INTO THE FIRST N POSITIONS OF C IYRCT IN DECREASING ORDER. C TO ACCOMPLISH THIS, MAKE IXRCT(I)=-1 IF IYCCT(J)=I FOR C SOME J, THEN SCAN IXRCT BACKWARDS. DO 79 J=1,N IYCJ=IYCCT(J) IXRCT(IYCJ)=-1 79 CONTINUE K=1 I=MP1 80 I=I-1 IF(I)83,83,81 81 IF(IXRCT(I)+1)80,82,80 82 IYRCT(K)=I K=K+1 GO TO 80 C NOW FILL IN THE REST OF IYRCT BY SCANNING IXRCT AGAIN. 83 I=MP1 84 I=I-1 IF(I)87,87,85 85 IF(IXRCT(I))84,86,86 86 IYRCT(K)=I K=K+1 GO TO 84 87 RETURN END SUBROUTINE SJELIM(L,LL,K,IR,IS,NPARM,NUMGR,V) C***BEGIN PROLOGUE SJELIM C***REFER TO SLNPRO C***ROUTINES CALLED (NONE) C***PURPOSE THIS SUBROUTINE PERFORMS A MODIFIED JORDAN C ELIMINATION ON THE L-LL+1 BY K MATRIX C CONSISTING OF ROWS LL THROUGH L OF V AND C COLUMNS 1 THROUGH K OF V. THE RESOLVENT C IS V(IR,IS). C***END PROLOGUE SJELIM C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION V(NUMGR+2*NPARM+1,NPARM+2) C C SET PRECISION DEPENDENT CONSTANTS FOR SJELIM. C***FIRST EXECUTABLE STATEMENT SJELIM ONE=1.0D0 C EMD OF SETTING PRECISION DEPENDENT CONSTANTS FOR SJELIM. C C DIVIDE THE ENTRIES IN THE RESOLVENT ROW (EXCEPT FOR THE C RESOLVENT) BY THE RESOLVENT. RESOL=V(IR,IS) DO 2 J=1,K IF(J-IS)1001,2,1001 1001 V(IR,J)=V(IR,J)/RESOL 2 CONTINUE C SWEEP OUT IN ALL BUT ROW IR AND COLUMN IS. DO 6 I=LL,L IF(I-IR)1002,6,1002 1002 FACT=-V(I,IS) DO 5 J=1,K IF(J-IS)1003,5,1003 1003 V(I,J)=V(I,J)+V(IR,J)*FACT 5 CONTINUE 6 CONTINUE C DIVIDE THE ENTRIES IN THE RESOLVENT COLUMN (EXCEPT FOR THE C RESOLVENT) BY THE NEGATIVE OF THE RESOLVENT. DO 8 I=LL,L IF(I-IR)1004,8,1004 1004 V(I,IS)=-V(I,IS)/RESOL 8 CONTINUE C REPLACE THE RESOLVENT BY ITS RECIPROCAL. V(IR,IS)=ONE/RESOL RETURN END SUBROUTINE SEARSL(IOPTN,NUMGR,NPARM,PRJLIM,TOL1,X,FUN,IFUN, *PTTBL,IPTB,INDM,PARAM,ERROR,RCHDWN,MACT,IACT,IPHSE,UNIT, *TOLCON,RCHIN,ITYPM1,ITYPM2,IWORK,LIWRK,WORK,LWRK,ERR1,PARPRJ, *PROJCT,EMIN,EMIN1,PARSER,NSRCH) C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION FUN(IFUN),PTTBL(IPTB,INDM),PARAM(NPARM),ERR1(NUMGR+3), *PARPRJ(NPARM),X(NPARM+1),ERROR(NUMGR+3), *IACT(NUMGR),PARSER(NPARM),IWORK(LIWRK),WORK(LWRK) C C THIS SUBROUTINE USES A MODIFIED QUADRATIC FITTING PROCESS TO C SEARCH FOR THE MINIMUM OF A FUNCTION F. IT REQURES AN INITIAL C GUESS IN PROJCT, A TOLERANCE TOL1 ON THE SEARCH INTERVAL LENGTH, C AN UPPER BOUND PRJLIM ON THE MINIMIZING POINT (WHICH SHOULD BE SET C VERY LARGE IF NO UPPER BOUND IS DESIRED), AND A WAY TO COMPUTE F(X) C FOR A GIVEN X. THE SUBROUTINE WILL PRINT A WARNING AND RETURN IF C IT WOULD NEED TO COMPUTE F MORE THAN INITLM TIMES IN THE INITIALIZATION C OR MORE THAN NADD ADDITIONAL TIMES IN THE MAIN PART OF THE PROGRAM. C WHEN THE SUBROUTINE RETURNS, IT WILL HAVE PUT THE MINIMUM LOCATION IN C PROJCT, THE MINIMUM F VALUE IN EMIN, THE F VALUE FOR THE INITIAL C PROJCT IN EMIN1, AND THE NUMBER OF TIMES IT COMPUTED F IN NSRCH. C C SET MACHINE AND PRECISION DEPENDENT CONSTANTS FOR SEARSL. ONE=1.0D0 TWO=ONE+ONE FOUR=TWO+TWO TEN=FOUR+FOUR+TWO C NWRIT=I1MACH(2) SPCMN=D1MACH(3) BIG=ONE/SPCMN TOLDEN=TEN*SPCMN TOL4=TOL1/FOUR BALFCT=TEN BALADJ=(TEN-ONE)/TEN ILC08=ILOC(8,NPARM,NUMGR) ILC10=ILOC(10,NPARM,NUMGR) ILC17=ILOC(17,NPARM,NUMGR) ILC21=ILOC(21,NPARM,NUMGR) ILC27=ILOC(27,NPARM,NUMGR) ILC29=ILOC(29,NPARM,NUMGR) ILC48=ILOC(48,NPARM,NUMGR) C C THE INITIAL PROJCT CAN BE INCREASED (OR DECREASED) BY A FACTOR OF C 2.0**((INITLM-1)*INITLM-2)/2) (ASSUMING WE TAKE INITLM .GE. 3, AS C WE SHOULD). WE TAKE INITLM=6 SINCE A FACTOR OF 1024 SEEMS SUFFICIENT. INITLM=6 C NADD=4 SEEMS TO BE SUFFICIENT SINCE THIS NUMBER OF ITERATIONS PAST THE C INITIALIZATION SEEMS TO ONLY RARELY BE EXCEEDED. NADD=4 NSRCH=0 ILF=0 IRT=0 IUPBAR=0 ISAVE=0 C INITIALLY PUT PARAM IN PARSER SO THERE WILL BE SOMETHING THERE IF C WE NEVER GET A CORRECTIBLE PARPRJ. DO 55 J=1,NPARM PARSER(J)=PARAM(J) 55 CONTINUE C WE NOW TRY TO COMPUTE VALUES AT POINTS P2=PROJCT, P1=P2/2.0, AND C P3=2.0*P2 (BUT P3 CANNOT EXCEED PRJLIM). P2=PROJCT C SET LLL=2 AS THE THREAD THROUGH THE MINOTAURS CAVERN AND JUMP C DOWN TO PUT F(P2) IN F2. WE WILL JUMP BACK AFTER ALL SUCH JUMPS LLL=2 PVAL=P2 GO TO 3500 C 77 F2=FVAL C SET EMIN1 = THE VALUE OF F USING THE GIVEN PROJECTION FACTOR PROJCT. EMIN1=FVAL P1=P2/TWO C SET LLL=1 AND PUT F(P1) IN F1. LLL=1 PVAL=P1 GO TO 3500 C 97 F1=FVAL P3=TWO*P2 C IF P3 .GT. PRJLIM, SET IUPBAR=1 AS AN INDICATOR WE CANNOT LATER C EXPAND THE INTERVAL TO THE RIGHT. THEN IF PRJLIM .GE. P2+TOL4 C REPLACE P3 BY PRJLIM, AND OTHERWISE EXPAND THE INTERVAL TO THE C LEFT TO GET THE DESIRED THIRD POINT. IF(P3-PRJLIM)160,160,120 120 IUPBAR=1 IF(PRJLIM-P2-TOL4)220,140,140 140 P3=PRJLIM C HERE SET LLL=3 AND PUT F(P3) IN F3. 160 LLL=3 PVAL=P3 GO TO 3500 C 187 F3=FVAL GO TO 280 C C EXPAND LEFT TO GET THE INITIAL THIRD POINT SINCE THERE IS NO ROOM C TO EXPAND RIGHT. 220 P3=P2 F3=F2 P2=P1 F2=F1 P1=P1/TWO C SET LLL=4 AND PUT F(P1) IN F1. LLL=4 PVAL=P1 GO TO 3500 C 247 F1=FVAL C C WE NOW HAVE FOUND P1, P2, AND P3 WITH CORRESPONDING VALUES C F1, F2, AND F3. WE EXPAND THE INTERVAL IF NECESSARY TO TRY C TO FIND NEW VALUES WITH F2 .LE. MIN(F1,F3). 280 IF(F2-F1)500,500,300 300 IF(F1-F3)320,320,520 C C HERE WE WILL EXPAND THE INTERVAL TO THE LEFT, PROVIDING THAT C NSRCH .LT. INITLM AND P1-P1/2.0**(NSRCH-1) .GE. TOL4. 320 IF(NSRCH-INITLM)340,360,360 340 IF(P1-P1/TWO**(NSRCH-1)-TOL4)360,380,380 C C HERE WE CANNOT EXPAND LEFT AND WE RETURN WITH THE BEST VALUES C FOUND SO FAR. 360 PROJCT=P1 EMIN=F1 RETURN C C EXPAND LEFT. 380 P3=P2 F3=F2 P2=P1 F2=F1 P1=P1/TWO**(NSRCH-1) C SET LLL=5 AND PUT F(P1) IN F1. LLL=5 PVAL=P1 GO TO 3500 407 F1=FVAL C C HERE F2 .LE. F3 AND WE HAVE JUST EXPANDED LEFT. IF F2 .GT. F1 WE C TRY TO EXPAND LEFT AGAIN, OTHERWISE WE CHECK TO SEE IF WE ARE DONE C INITIALIZING. IF(F2-F1)440,440,320 C C HERE WE CHECK TO SEE IF THE F COMPUTATION HAS FAILED EVERY TIME C (INDICATED BY F1=F2=F3=BIG), AND IF SO WE TRY TO EXPAND LEFT. C IF NOT, WE ARE DONE WITH THE INITIALIZATION. 440 IF(F1-BIG)1100,460,460 460 IF(F2-BIG)1100,480,480 480 IF(F3-BIG)1100,320,320 C C HERE F2 .LE. F1. IF F2 .LE. F3 AND WE HAVE NOT HAD ALL FAILURES OF C THE F COMPUTATION, WE ARE DONE INITIALIZING. 500 IF(F2-F3)440,440,520 C C HERE F3 .LT. MIN(F1,F2) AND WE EXPAND THE INTERVAL TO THE RIGHT IF C NSRCH .LT. INITLM AND IUPBAR=0. 520 IF(NSRCH-INITLM)540,560,560 540 IF(IUPBAR)580,580,560 C C HERE WE CANNOT EXPAND RIGHT AND WE RETURN WITH THE BEST VALUES C FOUND SO FAR. 560 PROJCT=P3 EMIN=F3 RETURN C C EXPAND RIGHT. 580 P1=P2 F1=F2 P2=P3 F2=F3 P3=TWO**(NSRCH-1)*P2 C IF P3 .GT. PRJLIM, SET IUPBAR=1 AS AN INDICATOR WE CANNOT LATER C EXPAND THE INTERVAL TO THE RIGHT. THEN IF PRJLIM .GE. P2+TOL4 C REPLACE P3 BY PRJLIM, AND OTHERWISE RETURN WITH THE BEST VALUES C FOUND SO FAR. IF(P3-PRJLIM)660,660,600 600 IUPBAR=1 IF(PRJLIM-P2-TOL4)620,640,640 620 PROJCT=P2 EMIN=F2 RETURN C 640 P3=PRJLIM C C SET LLL=6 AND PUT F(P3) IN F3. 660 LLL=6 PVAL=P3 GO TO 3500 687 F3=FVAL C C HERE F2 .LT. F1 AND WE HAVE JUST EXPANDED RIGHT. IF F2 .LE. F3 C WE ARE DONE INITIALIZING, OTHERWISE WE TRY TO EXPAND RIGHT AGAIN. IF(F2-F3)1100,1100,520 C END OF INITIALIZATION. C C ASSUMING THAT P3-P1 .GE. TOL1, WE NOW HAVE POINTS P1, P2, P3 WITH C P1 .LE. P2-TOL4, P2 .LE. P3-TOL4, F1=F(P1) .GE. F2=F(P2), AND F3=F(P3) C .GE. F2. THESE CONDITIONS WILL BE MAINTAINED THROUGHOUT THE PROGRAM. C SET LLL=7, WHERE IT WILL REMAIN FROM NOW ON. 1100 LLL=7 C C RESET LIMS1 SO THAT AT MOST NADD MORE COMPUTATIONS OF F WILL BE DONE. LIMS1=NSRCH+NADD C C IF WE HAVE COMPUTED F LIMS1 TIMES, WE PUT P2 IN PROJCT, PUT F2 IN C EMIN, AND RETURN. 1200 IF(NSRCH-LIMS1)1250,1300,1300 C C IF THE SEARCH INTERVAL LENGTH IS LESS THAN TOL1 WE PUT P2 IN C PROJCT, PUT F2 IN EMIN, AND RETURN. 1250 IF(P3-P1-TOL1)1300,1400,1400 C 1300 PROJCT=P2 EMIN=F2 RETURN C C COMPUTE S1 = THE ABSOLUTE VALUE OF THE SLOPE OF THE LINE THROUGH C (P1,F1) AND (P2,F2), AND S2 = THE (ABSOLUTE VALUE OF THE) SLOPE C OF THE LINE THROUGH (P2,F2) AND (P3,F3). C***MOD CONSIDER INCREASING TOL1 TO 10**4*SPCMN 1400 S1=(F1-F2)/(P2-P1) S2=(F3-F2)/(P3-P2) C IF S1+S2 IS VERY SMALL WE RETURN WITH THE BEST VALUES FOUND SO FAR. IF(S1+S2-TOLDEN)1300,1600,1600 C 1600 RLF=S2/(S1+S2) RRT=ONE-RLF C THE MINIMUM OF THE QUADRATIC POLYNOMIAL PASSING THROUGH C (P1,F1), (P2,F2), AND (P3,F3) WILL OCCUR AT (RLF*P1+ C RRT*P3+P2)/2.0. NOTE THAT THE THREE POINTS CANNOT BE C COLLNEAR, ELSE WE WOULD HAVE TERMINATED ABOVE. SINCE THE C MINIMUM OCCURS AT THE AVERAGE OF P2 AND A CONVEX COMBINATION C OF P1 AND P3, IT WILL BE AT LEAST AS CLOSE TO P2 AS TO THE C ENDPOINT ON THE SAME SIDE. IF(ILF-1)1800,1800,1700 C C HERE THE LEFT ENDPOINT WAS DROPPED AT THE LAST ILF .GT. 1 C ITERATIONS, SO TO PREVENT A LONG STRING OF SUCH OCCURRENCES C WITH LITTLE REDUCTION OF P3-P1 WE WILL SHIFT THE NEW POINT C TO THE RIGHT BY DECREASING RLF RELATIVE TO RRT. 1700 RLF=RLF/TWO**(ILF-1) RRT=ONE-RLF GO TO 2400 1800 IF(IRT-1)2000,2000,1900 C C HERE THE RIGHT ENDPOINT WAS DROPPED AT THE LAST IRT .GT. 1 C ITERATIONS, AND WE WILL SHIFT THE NEW POINT TO THE LEFT. 1900 RRT=RRT/TWO**(IRT-1) RLF=ONE-RRT GO TO 2400 C C HERE WE HAVE NOT JUST HAD A STRING OF TWO OR MORE MOVES IN C THE SAME DIRECTION, BUT IF THE SUBINTERVALS ARE OUT OF C BALANCE BY MORE THAN A FACTOR OF BALFCT, WE SHIFT THE NEW C POINT SLIGHTLY IN THE DIRECTION OF THE LONGER INTERVAL. THE C IDEA HERE IS THAT THE TWO CLOSE POINTS ARE PROBABLY NEAR THE C SOLUTION, AND IF WE CAN BRACKET THE SOLUTION WE MAY BE ABLE TO C CUT OFF THE MAJOR PORTION OF THE LONGER SUBINTERVAL. 2000 IF(P2-P1-BALFCT*(P3-P2))2200,2200,2100 C C HERE THE LEFT SUBINTERVAL IS MORE THAN BALFCT TIMES LONGER THAN C THE RIGHT SUBINTERVAL, SO WE DECREASE RRT RRELATIVE TO RLF. 2100 RRT=BALADJ*RRT RLF=ONE-RRT GO TO 2400 2200 IF(P3-P2-BALFCT*(P2-P1))2400,2400,2300 C C HERE THE RIGHT SUBINTERVAL IS MORE THAN BALFCT TIMES LONGER C THAN THE LEFT SUBINTERVAL, SO WE DECREASE RLF RELATIVE TO RRT. 2300 RLF=BALADJ*RLF RRT=ONE-RLF C C COMPUTE THE (POSSIBLY MODIFIED) MINIMUM OF THE QUADRATIC FIT. 2400 P4=(RLF*P1+RRT*P3+P2)/TWO C C THE NEXT SECTION (FROM HERE TO STATEMENT 2800) MODIFIES P4 IF NECESSARY C TO GET P1+TOL4 .LE. P2,P4 .LE. P3-TOL4 AND ABS(P4-P2) .GE. TOL4. IN C THE UNLIKELY EVENT THIS IS NOT POSSIBLE WE SET PROJCT=P2, EMIN=F2 C AND RETURN. C C IF ABS(P4-P2) .LT. TOL4 WE REDEFINE P4 BY MOVING TOL4 FROM C P2 INTO THE LONGER SUBINTERVAL. NOTE THAT THE LENGTH OF THIS C SUBINTERVAL MUST BE AT LEAST TOL1/2.0 = 2.0*TOL4, ELSE WE C WOULD HAVE TERMINATED EARLIER. IF(ABS(P4-P2)-TOL4)2500,2710,2710 2500 IF(P3-P2-(P2-P1))2700,2700,2600 2600 P4=P2+TOL4 C IF TOL4 WAS SMALL ENOUGH RELATIVE TO P2 THAT THE MACHINE THINKS P4 C STILL EQUALS P2, WHICH IS MORE LIKELY IF P2 IS LARGE, THIS COULD RESULT C IN A DIVIDE FAULT LATER. TO AVOID THIS, WE REDEFINE P4 AS THE AVERAGE C OF P2 AND P3 I