CONMAX consists of two programs for solving the problem
minimize w
subject to
abs(fun(i) - confun(i,1)) <= w if icntyp(i)=2,
confun(i,1) <= w if icntyp(i)=1,
confun(i,1) <= 0.0 if icntyp(i)=-1 or -2,
Where:
icntyp(i)=-1 the constraint is linear (i.e. the left side
consists of a linear combination of the parameters in the vector array
param plus a constant)icntyp(i)=-2 the constraint may be nonlinear.Thus we are choosing the parameters to minimize the left sides of the type 2 and type 1 (i.e. primary) constraints subject to the type -1 and type -2 (i.e. standard) constraints.
If there are no primary constraints the program will attempt to find a feasible point (that is, a vector param for which the type -1 and type -2 constraints are satisfied within a tolerance tolcon described below) which is relatively close to the user supplied initial approximation, then return.
iwork and work are changed by the program,
the user will need to reset these values each time conmax is called.| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(conmax_solver), | intent(inout) | :: | me | |||
| integer, | intent(in) | :: | ioptn |
THIS IS THE OPTION SWITCH, WHICH SHOULD BE SET TO 0 UNLESS ONE OR MORE OF THE EXTRA OPTIONS DESCRIBED BELOW IS USED. THE USER HAS SEVERAL EXTRA OPTIONS WHICH ARE ACTIVATED BY SETTING IOPTN TO A VALUE OTHER THAN 0; MORE THAN ONE AT A TIME CAN BE USED. IN PARTICULAR:
|
||
| integer, | intent(in) | :: | nparm |
THIS IS THE NUMBER OF PARAMETERS IN THE PROBLEM.
(THEY ARE STORED IN |
||
| integer, | intent(in) | :: | numgr |
THIS IS THE TOTAL NUMBER OF CONSTRAINTS. |
||
| integer, | intent(in) | :: | itlim |
THIS IS THE LIMIT ON THE NUMBER OF ITERATIONS, I.E. THE LIMIT ON THE NUMBER OF TIMES THE PROGRAM REDUCES W. IF ITLIM IS SET TO 0 THE PROGRAM WILL COMPUTE THE ERRORS FOR THE INITIAL APPROXIMATION AND STOP WITHOUT CHECKING FEASIBILITY. |
||
| real(kind=wp), | intent(in) | :: | fun(Ifun) |
(VECTOR ARRAY OF DIMENSION IFUN) THIS IS A VECTOR ARRAY IN WHICH DATA OR FUNCTION VALUES IN TYPE 2 CONSTRAINTS (SEE ABOVE) CAN BE STORED. FUN(I) NEED NOT BE ASSIGNED A VALUE IF IT IS NOT GOING TO BE USED. |
||
| integer, | intent(in) | :: | ifun |
THIS IS THE DIMENSION OF FUN IN THE DRIVER PROGRAM. IT MUST BE >= THE LARGEST INDEX I FOR WHICH FUN(I) IS USED UNLESS NO FUN(I) IS USED, IN WHICH CASE IT MUST BE >= 1. |
||
| real(kind=wp), | intent(in) | :: | pttbl(Iptb,Indm) |
(MATRIX ARRAY OF DIMENSION (IPTB,INDM)) ROW I OF THIS ARRAY NORMALLY CONTAINS A POINT USED IN THE ITH CONSTRAINT. THE ENTRIES IN ROW I NEED NOT BE ASSIGNED VALUES IF SUCH A POINT IS NOT USED IN THE ITH CONSTRAINT. (EXAMPLE: IF THE LEFT SIDE OF CONSTRAINT I IS A POLYNOMIAL IN ONE INDEPENDENT VARIABLE, THEN THE VALUE OF THE INDEPENDENT VARIABLE SHOULD BE IN PTTBL(I,1), AND THE COEFFICIENTS SHOULD BE IN PARAM.) PTTBL CAN ALSO BE USED TO PASS OTHER INFORMATION FROM THE DRIVER PROGRAM TO SUBROUTINE FNSET. |
||
| integer, | intent(in) | :: | iptb |
THIS IS THE FIRST DIMENSION OF PTTBL IN THE DRIVER PROGRAM. IT MUST BE >= THE LARGEST SUBSCRIPT I FOR WHICH A VALUE PTTBL(I,J) IS USED, AND MUST BE >= 1 IF NO SUCH VALUES ARE USED. |
||
| integer, | intent(in) | :: | Indm |
THIS IS THE SECOND DIMENSION OF PTTBL IN THE DRIVER PROGRAM. IT MUST BE >= THE LARGEST SUBSCRIPT J FOR WHICH A VALUE PTTBL(I,J) IS USED, AND MUST BE >= 1 IF NO SUCH VALUES ARE USED. |
||
| integer, | intent(inout) | :: | iwork(Liwrk) |
(VECTOR ARRAY OF DIMENSION LIWRK) THIS IS AN INTEGER WORK ARRAY. THE USER NEED NOT PLACE ANY VALUES IN IT, EXCEPT POSSIBLY CERTAIN OPTIONAL INFORMATION AS DESCRIBED BELOW. |
||
| integer, | intent(in) | :: | liwrk |
THIS IS THE DIMENSION OF IWORK IN THE DRIVER PROGRAM. IT MUST BE AT LEAST 7NUMGR + 7NPARM + 3. IF NOT, CONMAX WILL RETURN WITH THIS MINIMUM VALUE MULTIPLIED BY -1 AS A WARNING. |
||
| real(kind=wp), | intent(inout) | :: | work(Lwrk) |
(VECTOR ARRAY OF DIMENSION LWRK) THIS IS A FLOATING POINT WORK ARRAY. THE USER NEED NOT PLACE ANY VALUES IN IT, EXCEPT POSSIBLY CERTAIN OPTIONAL INFORMATION AS DESCRIBED BELOW. |
||
| integer, | intent(in) | :: | lwrk |
THIS IS THE DIMENSION OF WORK IN THE DRIVER PROGRAM. IT MUST BE AT LEAST 2NPARM2 + 4NUMGRNPARM + 11NUMGR + 27*NPARM + 13. IF NOT, CONMAX WILL RETURN WITH THIS MINIMUM VALUE MULTIPLIED BY -1 AS A WARNING. |
||
| integer, | intent(out) | :: | iter |
THIS IS THE NUMBER OF ITERATIONS PERFORMED BY CONMAX, INCLUDING THOSE USED IN ATTEMPTING TO GAIN FEASIBILITY, UNTIL EITHER IT CAN NO LONGER IMPROVE THE SITUATION OR THE ITERATION LIMIT IS REACHED. IF ITER=ITLIM IT IS POSSIBLE THAT THE PROGRAM COULD FURTHER REDUCE W IF RESTARTED (POSSIBLY WITH THE NEW PARAMETERS).
|
||
| real(kind=wp), | intent(inout) | :: | param(Nparm) |
(VECTOR ARRAY OF DIMENSION AT LEAST NPARM IN THE DRIVER PROGRAM) THE USER SHOULD PLACE AN INITIAL GUESS FOR THE PARAMETERS IN PARAM, AND ON OUTPUT PARAM WILL CONTAIN THE BEST PARAMETERS CONMAX HAS BEEN ABLE TO FIND. IF THE INITIAL PARAM IS NOT FEASIBLE THE PROGRAM WILL FIRST TRY TO FIND A FEASIBLE PARAM. |
||
| real(kind=wp), | intent(out) | :: | error(Numgr+3) |
(VECTOR ARRAY OF DIMENSION AT LEAST NUMGR + 3 IN THE DRIVER PROGRAM) FOR I=1,...,NUMGR, CONMAX WILL PLACE IN ERROR(I) THE ERROR IN CONSTRAINT I (DEFINED TO BE THE VALUE OF THE LEFT SIDE OF CONSTRAINT I, EXCEPT WITHOUT THE ABSOLUTE VALUE IN TYPE 2 CONSTRAINTS). FURTHER,
|
subroutine conmax(me, Ioptn, Nparm, Numgr, Itlim, Fun, Ifun, Pttbl, Iptb, & Indm, Iwork, Liwrk, Work, Lwrk, Iter, Param, Error) implicit none class(conmax_solver), intent(inout) :: me integer, intent(in) :: ifun !! THIS IS THE DIMENSION OF FUN IN THE DRIVER PROGRAM. !! IT MUST BE >= THE LARGEST INDEX I FOR WHICH FUN(I) IS USED !! UNLESS NO FUN(I) IS USED, IN WHICH CASE IT MUST BE >= 1. integer, intent(in) :: Indm !! THIS IS THE SECOND DIMENSION OF PTTBL IN THE DRIVER !! PROGRAM. IT MUST BE >= THE LARGEST SUBSCRIPT J FOR WHICH A !! VALUE PTTBL(I,J) IS USED, AND MUST BE >= 1 IF NO SUCH VALUES !! ARE USED. integer, intent(in) :: ioptn !! THIS IS THE OPTION SWITCH, WHICH SHOULD BE SET TO !! 0 UNLESS ONE OR MORE OF THE EXTRA OPTIONS DESCRIBED BELOW !! IS USED. THE USER HAS SEVERAL EXTRA OPTIONS WHICH ARE ACTIVATED BY SETTING !! IOPTN TO A VALUE OTHER THAN 0; MORE THAN ONE AT A TIME CAN BE USED. !! IN PARTICULAR: !! !! * IF THE ONES DIGIT OF IOPTN IS 0, THEN THE USER SHOULD GIVE FORMULAS !! IN FNSET FOR COMPUTING THE PARTIAL DERIVATIVES WHEN INDFN=1 AS DESCRIBED !! ABOVE. !! !! * IF THE USER SETS THE ONES DIGIT OF IOPTN TO 1, THEN INDFN WILL ALWAYS !! BE 0 WHEN FNSET IS CALLED, AND THE PROGRAM WILL AUTOMATICALLY !! APPROXIMATE THE PARTIAL DERIVATIVES WHEN REQUIRED USING THE FOLLOWING !! CENTERED DIFFERENCE FORMULA: !! PARTIAL DERIVATIVE WITH RESPECT TO THE JTH VARIABLE (WHERE 1 <= J !! <= NPARM) OF THE FUNCTION WHOSE VALUE IS COMPUTED IN CONFUN(IPT,1) !! (WHERE 1 <= IPT <= NUMGR) IS APPROXIMATELY THE VALUE OF THIS !! FUNCTION WHEN THE JTH COMPONENT OF PARAM IS INCREASED BY DELT, MINUS !! THE VALUE OF THIS FUNCTION WHEN THE JTH COMPONENT OF PARAM IS !! DECREASED BY DELT, ALL DIVIDED BY 2.0*DELT, WHERE DELT = SQRT(B**(-ITT)), !! WHERE B IS THE BASE FOR FOR FLOATING POINT NUMBERS AND ITT IS THE NUMBER !! OF BASE B DIGITS IN THE MANTISSA. !! !! * IF THE TENS DIGIT OF IOPTN IS 0, THE PROGRAM WILL NOT GIVE UP !! UNTIL EITHER AN ITERATION FAILS TO PRODUCE A REDUCTION ABS(ENCHG) OF !! MORE THAN 100.0*B**(-ITT) IN THE PRINCIPAL ERROR NORM, OR ITLIM !! ITERATIONS HAVE BEEN USED. !! !! * IF THE TENS DIGIT OF IOPTN IS 1, THE PROGRAM WILL ALSO GIVE UP IF !! ABS(ENCHG) < ENCSM FOR LIMSM CONSECUTIVE STEPS IN THE MAIN PART OF !! THE PROGRAM (I.E. AFTER TYPE -1 AND TYPE -2 FEASIBILITY HAVE BOTH BEEN !! ACHIEVED). THIS OPTION MAY BE USEFUL IN SAVING SOME TIME BY !! FORESTALLING A LONG STRING OF ITERATIONS AT THE END OF A RUN WITH ONLY !! A TINY IMPROVEMENT IN EACH ONE. ENCSM AND LIMSM ARE TRANSMITTED TO !! CONMAX IN WORK(1) AND IWORK(1) RESPECTIVELY. WORK(1) AND IWORK(1) DO !! NOT NEED TO BE ASSIGNED VALUES IF THE TENS DIGIT OF IOPTN IS 0. !! !! * IF THE HUNDREDS DIGIT OF IOPTN IS 0 OR 2, THEN THE INTERNAL PARAMETER !! NSTEP WILL BE GIVEN THE DEFAULT VALUE 1. NSTEP IS THE NUMBER OF !! RUNGE-KUTTA STEPS USED IN EACH RK ITERATION. !! !! * IF THE HUNDREDS DIGIT OF IOPTN IS 1 OR 3, THEN THE VALUE OF NSTEP USED !! WILL BE THE POSITIVE INTEGER VALUE PLACED IN IWORK(2) BY THE USER IN THE !! DRIVER PROGRAM. SETTING NSTEP LARGER THAN 1 MAY ALLOW THE PROGRAM TO !! SUCCEED ON DIFFICULT PROBLEMS WHERE THE CONVERGENCE WOULD BE EXTREMELY !! SLOW WITH NSTEP=1, BUT IT WILL GENERALLY CAUSE MORE COMPUTER TIME TO BE !! USED IN EACH RK ITERATION. IWORK(2) DOES NOT NEED TO BE ASSIGNED A !! VALUE IF THE HUNDREDS DIGIT OF IOPTN IS 0 OR 2. (NSTEP IS SOMETIMES !! CALLED THE OVERDRIVE PARAMETER.) !! !! * IF THE HUNDREDS DIGIT OF IOPTN IS 0 OR 1, THEN THE INTERNAL PARAMETER !! TOLCON WILL BE GIVEN THE DEFAULT VALUE SQRT(B**(-ITT)), WHERE B IS THE !! BASE FOR FLOATING POINT NUMBERS AND ITT IS THE NUMBER OF BASE B DIGITS !! IN THE MANTISSA. !! !! * IF THE HUNDREDS DIGIT OF IOPTN IS 2 OR 3, THEN THE VALUE OF TOLCON USED !! WILL BE THE VALUE PLACED IN WORK(2) BY THE USER IN THE DRIVER PROGRAM. !! THIS VALUE SHOULD ALWAYS BE NONNEGATIVE. IF THERE ARE NO TYPE -2 OR -1 !! CONSTRAINTS THEN THE SETTING OF TOLCON WILL HAVE NO EFFECT, BUT IF !! THERE ARE TYPE -2 OR -1 CONSTRAINTS THEN IN GENERAL SMALLER VALUES OF !! TOLCON MAY GIVE MORE ACCURACY IN THE FINAL ANSWER, BUT MAY SLOW DOWN !! OR PREVENT CONVERGENCE, WHILE LARGER VALUES OF TOLCON MAY HAVE THE !! REVERSE EFFECT. IN PARTICULAR, IF THE TYPE -2 AND -1 CONSTRAINTS !! CANNOT BE SATISFIED SUMULTANEOUSLY WITH STRICT INEQUALITY (THIS CASE !! COULD OCCUR, FOR EXAMPLE, IF AN EQUALITY CONSTRAINT G = 0.0 WAS !! ENTERED AS THE TWO INEQUALITY CONSTRAINTS G <= 0.0 AND !! -G <= 0.0), THEN SETTING TOLCON=0.0 WILL ALMOST CERTAINLY CAUSE THE !! PROGRAM TO FAIL, SINCE AT EACH ITERATION THE PROGRAM WILL NOT ACCEPT !! PROSPECTIVE NEW VALUES OF THE PARAMETERS UNLESS IT CAN CORRECT THEM !! BACK INTO THE RELAXED FEASIBLE REGION WHERE CONFUN(IPT,1) <= TOLCON !! FOR ALL THE TYPE -2 AND -1 CONSTRAINTS. !! !! * IF THE THOUSANDS DIGIT OF IOPTN IS 0, THE PROGRAM WILL USE THE RK METHOD !! (WHICH INVOLVES APPLYING FOURTH ORDER RUNGE-KUTTA TO A SYSTEM OF !! DIFFERENTIAL EQUATIONS), THEN IF THIS FAILS IT WILL TRY TO REDUCE !! W WITH AN SLP STEP (I.E. SUCCESSIVE LINEAR PROGRAMMING WITH A TRUST !! REGION), THEN GO BACK TO RK IF THE SLP STEP REDUCES W. !! !! * IF THE THOUSANDS DIGIT OF IOPTN IS 1, THE PROGRAM WILL USE SLP STEPS ONLY. !! IN GENERAL, IN SOME PROBLEMS SLP WORKS VERY WELL, AND IN THOSE !! PROBLEMS IT WILL USUALLY BE FASTER THAN RK BECAUSE IT REQUIRES FEWER !! GRADIENT EVALUATIONS THAN RK, BUT IN OTHER PROBLEMS THE CONVERGENCE !! OF SLP MAY BE EXCRUCIATINGLY SLOW, AND THE MORE POWERFUL RK METHOD !! MAY BE MUCH FASTER. !! !! * IF THE THOUSANDS DIGIT OF IOPTN IS 2 THE PROGRAM WILL USE THE RK METHOD !! ONLY, QUITTING WHEN RK CAN NO LONGER PRODUCE AN IMPROVEMENT. THIS !! MAY GIVE A LITTLE LESS ACCURACY THAN SETTING THE THOUSANDS DIGIT TO 0, !! BUT MAY SAVE SIGNIFICANT COMPUTER TIME IN SOME CASES. !! !! * IF THE TEN THOUSANDS DIGIT OF IOPTN IS 0, THEN FNSET SHOULD BE WRITTEN AS !! DESCRIBED ABOVE. !! !! * IF THE USER SETS THE TEN THOUSANDS DIGIT OF IOPTN TO 1, THEN FNSET SHOULD BE !! WRITTEN AS DESCRIBED ABOVE EXCEPT THAT THE COMPUTATIONS SHOULD BE DONE !! FOR ALL IPT=1,..,NUMGR INSTEAD OF FOR A SINGLE GIVEN VALUE OF IPT. !! THIS OPTION MAY SAVE COMPUTER TIME IN PROBLEMS WHERE A LARGE PART OF !! THE COMPUTATION IS THE SAME FOR DIFFERENT VALUES OF IPT, SINCE IT !! AVOIDS UNNECESSARY REPITITION OF THIS COMMON COMPUTATION BY HAVING !! THE LOOP OVER IPT IN FNSET INSTEAD OF OUTSIDE FNSET. !! IF THE TEN THOUSANDS DIGIT OF IOPTN IS 1, EVEN MORE TIME CAN OFTEN BE !! SAVED IF FNSET IS WRITTEN SO THAT ALL CONSTRAINTS ARE COMPUTED IF IPT=0, !! BUT ONLY THE STANDARD (I.E. TYPE -1 OR -2) CONSTRAINTS ARE COMPUTED IF !! IPT=-1. NOTE THAT IF THE TEN THOUSANDS DIGIT OF IOPTN IS 0, THEN IPT !! WILL BE POSITIVE WHENEVER FNSET IS CALLED, INDICATING THAT ONLY !! CONSTRAINT IPT SHOULD BE COMPUTED. !! THE DRAWBACK OF USING THIS OPTION IS THAT IN GENERAL SOME CONSTRAINT !! VALUES AND DERIVATIVES WILL BE COMPUTED UNNECESSARILY, WHICH COULD COST !! TIME IN SOME PROBLEMS. integer, intent(in) :: iptb !! THIS IS THE FIRST DIMENSION OF PTTBL IN THE DRIVER !! PROGRAM. IT MUST BE >= THE LARGEST SUBSCRIPT I FOR WHICH A !! VALUE PTTBL(I,J) IS USED, AND MUST BE >= 1 IF NO SUCH VALUES !! ARE USED. integer, intent(in) :: itlim !! THIS IS THE LIMIT ON THE NUMBER OF ITERATIONS, I.E. !! THE LIMIT ON THE NUMBER OF TIMES THE PROGRAM REDUCES W. IF !! ITLIM IS SET TO 0 THE PROGRAM WILL COMPUTE THE ERRORS FOR !! THE INITIAL APPROXIMATION AND STOP WITHOUT CHECKING !! FEASIBILITY. integer, intent(in) :: liwrk !! THIS IS THE DIMENSION OF IWORK IN THE DRIVER PROGRAM. !! IT MUST BE AT LEAST 7*NUMGR + 7*NPARM + 3. IF NOT, CONMAX WILL !! RETURN WITH THIS MINIMUM VALUE MULTIPLIED BY -1 AS A WARNING. integer, intent(in) :: lwrk !! THIS IS THE DIMENSION OF WORK IN THE DRIVER PROGRAM. !! IT MUST BE AT LEAST 2*NPARM**2 + 4*NUMGR*NPARM + 11*NUMGR + !! 27*NPARM + 13. IF NOT, CONMAX WILL RETURN WITH THIS MINIMUM !! VALUE MULTIPLIED BY -1 AS A WARNING. integer, intent(in) :: nparm !! THIS IS THE NUMBER OF PARAMETERS IN THE PROBLEM. !! (THEY ARE STORED IN `PARAM` -- SEE BELOW.) integer, intent(in) :: numgr !! THIS IS THE TOTAL NUMBER OF CONSTRAINTS. real(wp), intent(in) :: fun(Ifun) !! (VECTOR ARRAY OF DIMENSION IFUN) THIS IS !! A VECTOR ARRAY IN WHICH DATA OR FUNCTION VALUES IN TYPE 2 !! CONSTRAINTS (SEE ABOVE) CAN BE STORED. FUN(I) NEED NOT BE !! ASSIGNED A VALUE IF IT IS NOT GOING TO BE USED. real(wp), intent(in) :: pttbl(Iptb, Indm) !! (MATRIX ARRAY OF DIMENSION (IPTB,INDM)) !! ROW I OF THIS ARRAY NORMALLY CONTAINS A POINT USED IN THE ITH !! CONSTRAINT. THE ENTRIES IN ROW I NEED NOT BE ASSIGNED VALUES IF !! SUCH A POINT IS NOT USED IN THE ITH CONSTRAINT. !! (EXAMPLE: IF THE LEFT SIDE OF CONSTRAINT I IS A POLYNOMIAL IN !! ONE INDEPENDENT VARIABLE, THEN THE VALUE OF THE INDEPENDENT !! VARIABLE SHOULD BE IN PTTBL(I,1), AND THE COEFFICIENTS SHOULD BE !! IN PARAM.) !! PTTBL CAN ALSO BE USED TO PASS OTHER INFORMATION FROM THE DRIVER !! PROGRAM TO SUBROUTINE FNSET. integer, intent(inout) :: iwork(Liwrk) !! (VECTOR ARRAY OF DIMENSION LIWRK) THIS IS AN INTEGER !! WORK ARRAY. THE USER NEED NOT PLACE ANY VALUES IN IT, EXCEPT !! POSSIBLY CERTAIN OPTIONAL INFORMATION AS DESCRIBED BELOW. real(wp), intent(inout) :: work(Lwrk) !! (VECTOR ARRAY OF DIMENSION LWRK) THIS IS A FLOATING !! POINT WORK ARRAY. THE USER NEED NOT PLACE ANY VALUES IN IT, !! EXCEPT POSSIBLY CERTAIN OPTIONAL INFORMATION AS DESCRIBED BELOW. real(wp), intent(inout) :: param(Nparm) !! (VECTOR ARRAY OF DIMENSION AT LEAST NPARM !! IN THE DRIVER PROGRAM) THE USER SHOULD PLACE AN INITIAL GUESS !! FOR THE PARAMETERS IN PARAM, AND ON OUTPUT PARAM WILL CONTAIN !! THE BEST PARAMETERS CONMAX HAS BEEN ABLE TO FIND. IF THE !! INITIAL PARAM IS NOT FEASIBLE THE PROGRAM WILL FIRST TRY TO !! FIND A FEASIBLE PARAM. real(wp), intent(out) :: error(Numgr + 3) !! (VECTOR ARRAY OF DIMENSION AT LEAST NUMGR + 3 IN THE !! DRIVER PROGRAM) FOR I=1,...,NUMGR, CONMAX WILL PLACE IN !! ERROR(I) THE ERROR IN CONSTRAINT I (DEFINED TO BE THE VALUE !! OF THE LEFT SIDE OF CONSTRAINT I, EXCEPT WITHOUT THE ABSOLUTE !! VALUE IN TYPE 2 CONSTRAINTS). FURTHER, !! !! * ERROR(NUMGR+1) WILL BE THE (PRINCIPAL) ERROR NORM, THAT IS, THE !! MAXIMUM VALUE OF THE LEFT SIDES OF THE TYPE 2 (INCLUDING THE !! ABSOLUTE VALUE) AND TYPE 1 CONSTRAINTS. !! !! * ERROR(NUMGR+2) WILL BE THE MAXIMUM VALUE OF THE LEFT SIDES OF THE !! TYPE -1 CONSTRAINTS, OR 0.0 IF THERE ARE NO TYPE -1 !! CONSTRAINTS. EXCEPT FOR ROUNDOFF ERROR AND SMALL TOLERANCES !! IN SOME SUBROUTINES THIS VALUE WILL NORMALLY BE <= 0.0, AND !! IT WILL NOT BE ALLOWED TO BE > TOLCON IN THE MAIN PART OF !! THE PROGRAM. !! !! * ERROR(NUMGR+3) WILL BE THE MAXIMUM VALUE OF THE LEFT SIDES OF THE !! TYPE -2 CONSTRAINTS, OR 0.0 IF THERE ARE NO TYPE -2 !! CONSTRAINTS. THIS VALUE SHOULD BE <= TOLCON, SINCE THE !! PROGRAM WILL NOT EVEN ATTEMPT TO COMPUTE VALUES FOR THE !! TYPE 2 AND TYPE 1 CONSTRAINTS OTHERWISE (EXCEPT FOR VALUES !! CORRESPONDING TO THE INITIAL PARAMETERS PLACED IN PARAM BY !! THE USER). THE USER CAN USE THIS FEATURE TO INSERT TYPE -2 !! OR -1 CONSTRAINTS TO KEEP THE PARAMETERS AWAY FROM VALUES !! WHERE A TYPE 2 OR TYPE 1 CONSTRAINT IS UNDEFINED. integer, intent(out) :: iter !! THIS IS THE NUMBER OF ITERATIONS PERFORMED BY CONMAX, !! INCLUDING THOSE USED IN ATTEMPTING TO GAIN FEASIBILITY, !! UNTIL EITHER IT CAN NO LONGER IMPROVE THE SITUATION OR THE !! ITERATION LIMIT IS REACHED. IF ITER=ITLIM IT IS POSSIBLE !! THAT THE PROGRAM COULD FURTHER REDUCE W IF RESTARTED !! (POSSIBLY WITH THE NEW PARAMETERS). !! !! * ITER=-1 IS A SIGNAL THAT TYPE -1 FEASIBILITY COULD NOT BE !! ACHIEVED, IN THIS CASE ERROR WILL CONTAIN THE VALUES COMPUTED !! USING THE USER SUPPLIED INITIAL APPROXIMATION. !! * ITER=-2 IS A SIGNAL THAT TYPE -1 FEASIBILITY WAS ACHIEVED !! BUT TYPE -2 FEASIBILITY COULD NOT BE ACHIEVED, !! IN THIS CASE VALUES IN ERROR CORRESPONDING TO TYPE 1 OR !! TYPE 2 CONSTRAINTS WILL BE ZERO. real(wp) :: enc1, enchg, encsm, enor2, enor3, enorm, & prjslp, projct, & rchdnk, rchdwn, rchin, s real(wp) :: tolcon, tollin, wdist integer :: i, i1, ii, ilc02, ilc06, ilc08, ilc11, & ilc12, ilc13, ilc14, ilc15, ilc17, ilc20, ilc21, & ilc22, ilc24, ilc25, ilc26, ilc27 integer :: ilc29, ilc30, ilc31, ilc33, ilc35, ilc40, ilc42, & ilc44, ilc46, iophun, iopten, ioptho, & ioptth, iphse, ipmax, ipt integer :: irk, ismax, isucc, itersl, itlim1, & ityp1, ityp1k, ityp2, ityp2k, itypm1, itypm2, & j, j1, j2, jflag, jiwrk, jj integer :: jwrk, kntsm, l, l1, l2, limsm, m, & mact1, ncor, nmaj, nmin, npar1, nstep, & numlim ! check to see if the dimensions liwrk and lwrk are large enough. jiwrk = 7*Numgr + 7*Nparm + 3 jwrk = 2*Nparm**2 + 4*Numgr*Nparm + 11*Numgr + 27*Nparm + 13 if (Liwrk < jiwrk .or. Lwrk < jwrk) then iter = -999 ! error flag return else ! 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) ! IF THE TENS DIGIT OF IOPTN IS 1, SET KNTSM TO 0 AND GET ENCSM ! FROM WORK(1) AND LIMSM FROM IWORK(1). iopten = (Ioptn - (Ioptn/100)*100)/10 if (iopten > 0) then kntsm = 0 encsm = Work(1) limsm = Iwork(1) end if ! IF THE HUNDREDS DIGIT OF IOPTN IS 1 OR 3, SET NSTEP = IWORK(2), ! AND OTHERWISE SET NSTEP TO ITS DEFAULT VALUE OF 1. iophun = (Ioptn - (Ioptn/1000)*1000)/100 if (iophun <= (iophun/2)*2) then nstep = 1 else nstep = Iwork(2) end if ! IF THE HUNDREDS DIGIT OF IOPTN IS 2 OR 3, SET TOLCON = WORK(2), ! AND OTHERWISE SET TOLCON TO ITS DEFAULT VALUE OF SQRT(SPCMN). if (iophun < 2) then tolcon = sqrt(spcmn) else tolcon = Work(2) end if ! IN THIS VERSION OF CONMAX WE SET THE LINEAR CONSTRAINT TOLERANCE ! EQUAL TO THE NONLINEAR CONSTRAINT TOLERANCE. tollin = tolcon ! SET IRK=1 IF THE THOUSANDS DIGIT OF IOPTN IS 0 AND OTHERWISE SET IRK=0. ioptho = (Ioptn - (Ioptn/10000)*10000)/1000 if (ioptho <= 0) then irk = 1 else irk = 0 end if ! COMPUTE THE TEN THOUSANDS DIGIT OF IOPTN FOR LATER USE. ioptth = (Ioptn - (Ioptn/100000)*100000)/10000 ! SET IPHSE=-1 TO INDICATE WE HAVE NOT CHECKED TYPE -1 FEASIBILITY YET. iphse = -1 ! SET RCHDWN = THE NUMBER OF LENGTHS OF PROJCT IN RKSACT (OR NUMBER OF ! LENGTHS OF BNDLGT IN SETU1) WE WILL GO BELOW ERROR(NUMGR+1) TO DECLARE ! A PRIMARY CONSTRAINT TO BE ACTIVE. rchdwn = two rchdnk = rchdwn ! SET RCHIN = THE NUMBER OF LENGTHS OF PROJCT (OR BNDLGT) WE WILL GO ! BELOW 0.0 TO DECLARE A TYPE -2 CONSTRAINT TO BE ACTIVE. rchin = two ! SET A NORMAL VALUE FOR NUMLIM FOR USE IN SLPCON. numlim = 11 end if ! END OF PRELIMINARY SECTION. THE STATEMENTS ABOVE THIS POINT WILL NOT ! BE EXECUTED AGAIN IN THIS CALL TO CONMAX. ! CALL ERCMP1 WITH ICNUSE=0 TO COMPUTE THE ERRORS, ERROR NORMS, AND ICNTYP. ! WE TAKE IPHSE AS 0 SO ALL CONSTRAINTS WILL BE COMPUTED BY FNSET IN CASE ! THE TEN THOUSANDS DIGIT OF IOPTN IS 1. ! THIS IS ONE OF ONLY TWO PLACES IN THE PROGRAM WHERE WE CALL ERCMP1 WITH ! ICNUSE=0, THE OTHER BEING STATEMENT 1415 BELOW.. 100 call me%ercmp1(Ioptn, Nparm, Numgr, Fun, Ifun, Pttbl, Iptb, Indm, Param, 0, 0, & Iwork, Liwrk, Work(ilc08), Iwork(ilc17), ipmax, ismax, & Error) ! IF ITLIM=0 WE RETURN. if (Itlim <= 0) then return else ! COMPUTE ITYP2, ITYP1, ITYPM1, AND ITYPM2 AS THE NUMBER OF CONSTRAINTS OF ! TYPE 2 (I.E. PRIMARY, ABS(FUN(I)-CONFUN(I,1)) <= W) OR 1 (I.E. PRIMARY, ! CONFUN(I,1) <= W) OR -1 (I.E. STANDARD LINEAR, CONFUN(I,1) <= 0.0) ! OR -2 (I.E. STANDARD NONLINEAR) RESPECTIVELY. ityp2 = 0 ityp1 = 0 itypm1 = 0 itypm2 = 0 ! NOTE THAT ARRAYS NOT IN THE CALLING SEQUENCE FOR CONMAX ARE ACCESSED ! THROUGH THEIR LOCATION IN IWORK OR WORK. CONMAX IS THE ONLY ! SUBROUTINE IN WHICH THIS IS NECESSARY. do i = 1, Numgr ii = ilc17 - 1 + i ! HERE IWORK(II)=ICNTYP(I). if (Iwork(ii) < 0) then if (Iwork(ii) + 1 < 0) then itypm2 = itypm2 + 1 else itypm1 = itypm1 + 1 end if else if (Iwork(ii) /= 0) then if (Iwork(ii) <= 1) then ityp1 = ityp1 + 1 else ityp2 = ityp2 + 1 end if end if end do end if ! COMPUTE THE ERROR NORMS. ENORM IS THE PRINCIPAL ERROR NORM. 200 enorm = Error(Numgr + 1) enor2 = Error(Numgr + 2) enor3 = Error(Numgr + 3) ! WRITE ITER, ISUCC, IRK, ENCHG, AND THE ERROR NORMS. !1050 WRITE(NWRIT,1100) ITER,ISUCC,IRK,ENCHG,ENORM,ENOR2,ENOR3 ! WRITE(9,1100) ITER,ISUCC,IRK,ENCHG,ENORM,ENOR2,ENOR3 !1100 FORMAT(/8H ITER IS,I5,10H ISUCC IS,I4,8H IRK IS,I4, ! *10H ENCHG IS,E24.14/9H ENORM IS,E24.14,10H ENOR2 IS,E24.14/ ! *9H ENOR3 IS,E24.14) ! THE NEXT SECTION DETERMINES WHETHER WE WILL TERMINATE DUE TO ITERATION ! COUNT, AND IF SO FOR OUTPUT PURPOSES IT MODIFIES ITER (OR TWO OF THE ! ERROR NORMS IF THE FAILURE IS DUE TO INABILITY TO GAIN TYPE -2 ! FEASIBILITY). ! IF IOPTEN=1 AND WE HAVE DONE AT LEAST ONE ITERATION IN THE MAIN PART ! OF CONMAX, WE WILL GIVE UP IF ABS(ENCHG) HAS BEEN LESS THAN ENCSM FOR ! LIMSM CONSECUTIVE MAIN ITERATIONS (INCLUDING THIS ONE). 300 if (iopten == 1) then if (iphse == 0) then if (Iter > 0) then if (-enchg < encsm) then kntsm = kntsm + 1 if (kntsm >= limsm) then ! FOR OUTPUT PURPOSES REPLACE ITER BY ITER + ITLIM - ITLIM1, THE TRUE ! NUMBER OF ITERATIONS COUNTING INITIALIZATION. ITLIM - ITLIM1 WILL BE ! THE NUMBER OF ITERATIONS NEEDED TO GAIN TYPE -2 FEASIBILITY. WORK ! DONE TO GAIN TYPE -1 FEASIBILITY IS NOT COUNTED AS AN ITERATION. Iter = Iter + Itlim - itlim1 return end if else kntsm = 0 end if end if end if end if if (Iter < itlim1) then ! HERE ITER < ITLIM1. IF IPHSE = 0 OR -2 HERE WE GO INTO THE ! ITERATIVE PHASE OF CONMAX. if (iphse + 1 /= 0) goto 900 ! HERE IPHSE=-1 AND WE CHECK TYPE -1 FEASIBILITY, TRY TO REGAIN IT IF ! WE DONT HAVE IT, CHECK TYPE -2 FEASIBILITY, AND SET UP FOR TYPE -2 ! FEASIBILITY ITERATIONS IF WE DONT HAVE IT. THE STATEMENTS FROM HERE ! DOWN TO THE TRIPLE BLANK LINE WILL BE EXECUTED AT MOST ONCE. ! NOTE THAT ENOR2=0.0 IF THERE ARE NO TYPE -1 CONSTRAINTS. if (enor2 <= tollin) then ! HERE WE HAD TYPE -1 FEASIBILITY INITIALLY. if (enor3 > tolcon) goto 700 goto 800 else ! HERE WE DO NOT HAVE TYPE -1 FEASIBILITY SO WE TRY TO GET IT. ! WE WILL NEED TO TELL DERST TO COMPUTE THE VALUES OF THE LEFT SIDES ! OF THE TYPE -1 CONSTRAINTS WITH THE VARIABLES EQUAL TO ZERO (I.E. ! THE CONSTANT TERMS IN THE CONSTRAINTS), SO WE SET PARWRK TO THE ! ZERO VECTOR TO CARRY THE MESSAGE. do j = 1, Nparm jj = ilc27 - 1 + j ! HERE WORK(JJ) = PARWRK(J). Work(jj) = zero end do if (ioptth > 0) then ! HERE IOPTTH=1 AND WE CALL DERST WITH IPT=-1 TO PUT ALL THE STANDARD ! CONSTRAINT AND DERIVATIVE VALUES IN CONFUN. ! WE SET IPT=-1 TO TELL DERST IT NEED ONLY COMPUTE STANDARD CONSTRAINTS. ipt = -1 call me%derst(Ioptn, Nparm, Numgr, Pttbl, Iptb, Indm, Work(ilc27), & ipt, Work(ilc24), Work(ilc35), Iwork(ilc22), & Work(ilc08)) end if m = 0 do i = 1, Numgr ii = ilc17 - 1 + i ! HERE WE CONSIDER ONLY TYPE -1 CONSTRAINTS. THERE MUST BE AT LEAST ! ONE OF THESE, SINCE OTHERWISE WE WOULD NOT BE HERE ATTEMPTING TO ! GAIN TYPE -1 FEASIBILITY. ! HERE IWORK(II)=ICNTYP(I). if (Iwork(ii) + 1 == 0) then m = m + 1 if (ioptth <= 0) then ! HERE IOPTTH=0 AND WE HAVE NOT YET CALLED DERST TO PUT CONSTRAINT I ! AND ITS DERIVATIVES IN CONFUN, SO WE DO IT NOW. ipt = i call me%derst(Ioptn, Nparm, Numgr, Pttbl, Iptb, Indm, & Work(ilc27), ipt, Work(ilc24), Work(ilc35), & Iwork(ilc22), Work(ilc08)) end if ! COPY THE DERIVATIVES INTO PMAT FOR USE BY WOLFE. do l = 1, Nparm l1 = ilc29 - 1 + l + (m - 1)*npar1 l2 = ilc08 - 1 + i + l*Numgr ! HERE WORK(L1)=PMAT(L,M) AND WORK(L2)=CONFUN(I,L+1). Work(l1) = Work(l2) end do ! NOW THE ITH CONSTRAINT (WHICH IS ALSO THE MTH TYPE -1 CONSTRAINT) HAS ! THE FORM PMAT(1,M)*Z1+...+PMAT(NPARM,M)*ZNPARM + CONFUN(I,1) <= ! 0.0. WE MAKE THE CHANGE OF VARIABLES ZZ = Z - PARAM TO TRANSLATE THE ! ORIGIN TO PARAM. THE ITH CONSTRAINT WILL THEN HAVE THE FORM ! PMAT(1,M)*ZZ1+...+PMAT(NPARM,M)*ZZNPARM + (CONFUN(I,1) + PMAT(1,M)* ! PARAM(1)+...+PMAT(NPARM,M)*PARAM(NPARM)) <= 0.0. AFTER WOLFE FINDS ! THE CLOSEST POINT TO THE ORIGIN IN THE POLYHEDRON DEFINED BY THE NEW ! CONSTRAINTS, WE WILL ADD PARAM TO TRANSLATE BACK TO THE POINT WE WANT. l1 = ilc29 - 1 + npar1 + (m - 1)*npar1 l2 = ilc08 - 1 + i ! HERE WORK(L1)=PMAT(NPAR1,1) AND WORK(L2)=CONFUN(I,1). Work(l1) = Work(l2) do l = 1, Nparm l2 = ilc29 - 1 + l + (m - 1)*npar1 ! HERE WORK(L1)=PMAT(NPAR1,1) AND WORK(L2)=PMAT(L,M). Work(l1) = Work(l1) + Work(l2)*Param(l) end do end if end do ! CALL WOLFE WITH ISTRT=0 TO COMPUTE THE SOLUTION IN THE ZZ COORDINATE ! 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 > 0) goto 600 ! HERE JFLAG <= 0 AND WE PUT PARAM+WPT IN PARWRK TO CHECK WHETHER ! THE TYPE -1 CONSTRAINTS ARE NOW FEASIBLE WITHIN TOLLIN. do j = 1, Nparm j1 = ilc27 - 1 + j j2 = ilc42 - 1 + j ! HERE WORK(J1)=PARWRK(J) AND WORK(J2)=WPT(J). Work(j1) = Param(j) + Work(j2) end do ! FOR USE IN ERCMP1 WE SET JCNTYP(I)=-1 IF ICNTYP(I)=-1 AND SET ! JCNTYP(I)=0 OTHERWISE. do i = 1, Numgr ii = ilc17 - 1 + i jj = ilc21 - 1 + i ! HERE IWORK(II)=ICNTYP(I) AND IWORK(JJ)=JCNTYP(I). if (Iwork(ii) + 1 /= 0) then Iwork(jj) = 0 else Iwork(jj) = -1 end if end do ! CALL ERCMP1 WITH ICNUSE=1. call me%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) ! HERE WORK(I1)=ERR1(NUMGR+2). if (Work(i1) > tollin) goto 600 ! HERE WE HAVE ACHIEVED TYPE -1 FEASIBILITY. WE REPLACE PARAM WITH ! PARWRK. do j = 1, Nparm jj = ilc27 - 1 + j ! HERE WORK(JJ)=PARWRK(J). Param(j) = Work(jj) end do ii = ilc11 - 1 + Numgr + 2 ! HERE WORK(II)=ERR1(NUMGR+2). ! WRITE(NWRIT,1397)WORK(II),(PARAM(J),J=1,NPARM) !1397 FORMAT(48H TYPE -1 FEASIBILITY ACHIEVED. ERR1(NUMGR+2) IS, ! *E15.5,10H PARAM IS/(4E20.12)) ! IF THERE ARE TYPE -2 CONSTRAINTS, SET JCNTYP AS ICNTYP WITH ALL BUT -2 ! VALUES ZEROED OUT AND CALL ERCMP1 WITH ICNUSE=1 TO CHECK TYPE -2 ! FEASIBILITY. WE CANNOT SIMPLY CHECK THE OLD ENOR3 HERE SINCE PARAM HAS ! BEEN CHANGED. IF THERE ARE NO TYPE -2 CONSTRAINTS WE WILL AUTOMATICALLY ! HAVE TYPE -2 FEASIBILITY. if (itypm2 <= 0) then ! HERE WE HAVE BOTH TYPE -1 AND TYPE -2 FEASIBILITY, BUT PARAM WAS ! CHANGED IN GETTING TYPE -1 FEASIBILITY, SO WE CALL ERCMP1 ! WITH ICNUSE=0 (ICNUSE=1 WOULD WORK ALSO SINCE ICNTYP HAS NOT BEEN ! CHANGED HERE) TO GET THE NEW ERROR VECTOR. call me%ercmp1(Ioptn, Nparm, Numgr, Fun, Ifun, Pttbl, Iptb, Indm, & Param, 0, iphse, Iwork, Liwrk, Work(ilc08), & Iwork(ilc17), ipmax, ismax, Error) goto 800 else do i = 1, Numgr ii = ilc17 - 1 + i jj = ilc21 - 1 + i ! HERE IWORK(II)=ICNTYP(I) AND IWORK(JJ)=JCNTYP(I). if (Iwork(ii) + 1 < 0) then Iwork(jj) = -2 else Iwork(jj) = 0 end if end do call me%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 ! HERE WORK(II)=ERR1(NUMGR+3). if (Work(ii) > tolcon) goto 700 call me%ercmp1(Ioptn, Nparm, Numgr, Fun, Ifun, Pttbl, Iptb, Indm, & Param, 0, iphse, Iwork, Liwrk, Work(ilc08), & Iwork(ilc17), ipmax, ismax, Error) goto 800 end if end if ! HERE ITER = ITLIM1, SO WE RETURN. else if (iphse >= 0) then Iter = Iter + Itlim - itlim1 return end if ! HERE WE HAVE FAILED TO ACHIEVE TYPE -2 FEASIBILITY AND WE SET ITER=-2 ! AS A WARNING, PUT ERROR(NUMGR+1) IN ITS PROPER LOCATION, SET ! ERROR(NUMGR+1) = 0.0 SINCE THE PRIMARY CONSTRAINTS WERE NOT COMPUTED, ! AND RETURN. NOTE THAT WE CANNOT HAVE IPHSE=-1 HERE SINCE THAT WOULD ! IMPLY ITER=0, THUS ITLIM=ITLIM1=0, IN WHICH CASE WE WOULD HAVE ! TERMINATED EARLIER. 400 Iter = -2 Error(Numgr + 3) = Error(Numgr + 1) Error(Numgr + 1) = zero ! write(nwrit,'(A)') '***WARNING NONLINEAR STANDARD FEASIBILITY NOT ACHIEVED***' return ! HERE WE HAVE FAILED TO ACHIEVE TYPE -1 FEASIBILITY. WE SET ITER=-1 ! AS A WARNING AND RETURN. 600 Iter = -1 ! WRITE(NWRIT,'(A)') '***WARNING LINEAR STANDARD FEASIBILITY NOT ACHIEVED***' return ! HERE WE HAVE TYPE -1 FEASIBILITY BUT NOT TYPE -2 FEASIBILITY. WE SET ! UP FOR THE TYPE -2 FEASIBILITY ITERATIONS, IN WHICH TYPE 1 AND TYPE ! 2 CONSTRAINTS ARE IGNORED AND TYPE -2 CONSTRAINTS ARE TREATED AS ! TYPE 1 CONSTRAINTS, EXCEPT WE WILL SWITCH OVER TO NORMAL ITERATIONS ! ONCE WE CAN FORCE W <= TOLCON. THUS WE SET THE INDICATOR IPHSE TO ! -2, RESET ICNTYP(I) TO 1 IF IT WAS -2, LEAVE IT AT -1 IF IT WAS -1, ! AND SET IT TO 0 OTHERWISE, RESET ITYP2, ITYP1, AND ITYPM2, AND CALL ! ERCMP1 WITH ICNUSE=1 TO PUT THE PROPER VALUES IN ERROR. 700 iphse = -2 do i = 1, Numgr ii = ilc17 - 1 + i ! HERE IWORK(II)=ICNTYP(I). if (Iwork(ii) + 1 < 0) then Iwork(ii) = 1 else if (Iwork(ii) + 1 /= 0) then Iwork(ii) = 0 end if end do ! SAVE ITYP2 AND ITYP1. ityp2k = ityp2 ityp1k = ityp1 ityp2 = 0 ityp1 = itypm2 itypm2 = 0 call me%ercmp1(Ioptn, Nparm, Numgr, Fun, Ifun, Pttbl, Iptb, Indm, Param, 1, & iphse, Iwork, Liwrk, Work(ilc08), Iwork(ilc17), ipmax, & ismax, Error) goto 900 ! HERE WE HAVE BOTH TYPE -1 AND TYPE -2 FEASIBILITY, AND WE ! SET IPHSE=0 AND GO INTO THE MAIN PART OF CONMAX (UNLESS THERE WERE ! NO TYPE 1 OR TYPE 2 CONSTRAINTS, IN WHICH CASE WE RETURN). 800 iphse = 0 if (ityp1 + ityp2 <= 0) return ! END OF INITIAL FEASIBILITY CHECKING, TYPE -1 FEASIBILITY WORK, AND ! TYPE -2 SETUP. THE BLOCK OF STATEMENTS FROM HERE UP TO THE ! PRECEDING DOUBLE BLANK LINE WILL NOT BE EXECUTED AGAIN. 900 if (irk <= 0) then ! HERE IRK IS 0 OR -1 AND WE DO AN SLP STEP. IF SLPCON CANNOT REDUCE THE ! PRINCIPAL ERROR NORM ENORM = ERROR(NUMGR+1) BY MORE THAN 100.0*B**(-ITT) ! THEN IT WILL LEAVE PARAM AND ERROR UNCHANGED. call me%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) else ! HERE IRK IS 1 OR 2 AND WE DO AN RK STEP. IF RKCON CANNOT REDUCE THE ! PRINCIPAL ERROR NORM ENORM = ERROR(NUMGR+1) BY MORE THAN 100.0*B**(-ITT) ! THEN IT WILL LEAVE PARAM AND ERROR UNCHANGED. call me%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) end if if (isucc <= 0) then ! HERE THE RK OR SLP STEP REDUCED ERROR(NUMGR+1) BY MORE THAN ! 100.0*B**(-ITT), AND WE INCREMENT ITER. Iter = Iter + 1 ! IF EITHER IPHSE=0, OR IPHSE=-2 AND ERROR(NUMGR+1) > TOLCON, WE GO ! ON AS USUAL TO SET UP ANOTHER STEP WITH THE SAME IPHSE. if (iphse < 0) then if (Error(Numgr + 1) <= tolcon) then ! HERE IPHSE=-2 AND ERROR(NUMGR+1) <= TOLCON, SO WE HAVE JUST ACHIEVED ! TYPE -2 FEASIBILITY. WE WILL SET IPHSE=0, AND IF THERE ARE ANY ! PRIMARY CONSTRAINTS WE WILL RESET ITER, ITERSL, AND ITLIM1 (SINCE ! ITER=0 AND ITERSL=0 HAVE MEANINGS TO RKCON AND SLPCON RESPECTIVELY), ! RESET RCHIN AND RCHDWN, AND GO BACK TO THE FIRST ERCMP1 CALL TO ! RESTORE ERROR AND ICNTYP (ITYP1, ITYP2, ITYPM1, AND ITYPM2 WILL ALSO ! BE RESTORED). iphse = 0 if (ityp1k + ityp2k <= 0) return itlim1 = Itlim - Iter Iter = 0 itersl = 0 rchin = rchdwn rchdwn = rchdnk goto 100 end if end if if (irk < 0) then ! HERE WE HAD AN SLP SUCCESS AND WE ARE GOING TO TRY RK AGAIN, SO WE SET ! IRK=2 TO WARN RKCON THAT THE SUCCESS CAME FROM SLP. irk = 2 else if (irk /= 0) then ! HERE IRK IS 1 OR 2, SO WE JUST HAD AN RK SUCCESS. WE RESET IRK AND ! ITERSL. irk = 1 itersl = 0 goto 200 end if ! HERE WE HAD AN SLP SUCCESS AND WE INCREMENT ITERSL = THE NUMBER OF SLP ! SUCCESSES SINCE THE LAST SUCCESSFUL RK STEP (IF ANY). ITERSL IS NEEDED ! IN SUBROUTINE BNDSET (CALLED BY SLPCON). itersl = itersl + 1 goto 200 else ! HERE RKCON OR SLPCON FAILED TO SIGNIFICANTLY REDUCE THE PRINCIPAL ERROR ! NORM. IF WE JUST TRIED SLP WE QUIT, AND IF WE JUST TRIED RK WE ATTEMPT ! AN SLP STEP UNLESS IOPTHO = 2, IN WHICH CASE WE QUIT. if (irk > 0) then if (ioptho /= 2) then irk = -1 goto 300 end if end if ! IF IPHSE=-2 HERE WE WILL SET ITER=-2 AS A WARNING AND CHANGE ! ERROR(NUMGR+1) AND ERROR(NUMGR+3) BEFORE RETURNING. OTHERWISE WE WILL ! HAVE IPHSE=0 AND WE WILL ADJUST ITER BEFORE RETURNING. if (iphse < 0) goto 400 Iter = Iter + Itlim - itlim1 end if end subroutine conmax