diff options
Diffstat (limited to 'Donjon/src/PLQ.f')
| -rw-r--r-- | Donjon/src/PLQ.f | 628 |
1 files changed, 628 insertions, 0 deletions
diff --git a/Donjon/src/PLQ.f b/Donjon/src/PLQ.f new file mode 100644 index 0000000..91b1392 --- /dev/null +++ b/Donjon/src/PLQ.f @@ -0,0 +1,628 @@ +*DECK PLQ + SUBROUTINE PLQ(NENTRY,HENTRY,IENTRY,JENTRY,KENTRY) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Solves a linear optimization problem with a quadratic constraint. +* PLQ = Quasi Linear Programmation (aka Optex method) +* +*Copyright: +* Copyright (C) 2002 Ecole Polytechnique de Montreal +* This library is free software; you can redistribute it and/or +* modify it under the terms of the GNU Lesser General Public +* License as published by the Free Software Foundation; either +* version 2.1 of the License, or (at your option) any later version +* +*Author(s): +* R. Chambon +* +*Parameters: input +* NENTRY number of data structures transfered to this module. +* HENTRY name of the data structures. +* IENTRY data structure type where: +* IENTRY=1 for LCM memory object; +* IENTRY=2 for XSM file; +* IENTRY=3 for sequential binary file; +* IENTRY=4 for sequential ASCII file. +* JENTRY access permission for the data structure where: +* JENTRY=0 for a data structure in creation mode; +* JENTRY=1 for a data structure in modifications mode; +* JENTRY=2 for a data structure in read-only mode. +* KENTRY data structure pointer. +* +*Comments: +* The calling specifications are: +* OPTIM := PLQ: OPTIM :: (plq\_data) ; +* where +* OPTIM : name of the \emph{optimize} object (L\_OPTIMIZE signature) +* containing the optimization informations. Object OPTIM must appear on +* both LHS and RHS to be able to update the previous values. +* (plq\_data) : structure containing the data to the module PLQ:. +* +*----------------------------------------------------------------------- +* + USE GANLIB + IMPLICIT NONE +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER NENTRY,IENTRY(NENTRY),JENTRY(NENTRY) + TYPE(C_PTR) KENTRY(NENTRY) + CHARACTER HENTRY(NENTRY)*12 +*---- +* LOCAL VARIABLES +*---- + INTEGER NSTATE + PARAMETER (NSTATE=40) + INTEGER NITMA,ITYP,ITYP1,ICONV,ICST,IEDSTP,LENGT1,LENGT2,ITYLCM, + 1 ISTEP,IVAR + REAL FLOTT,ECOSTR + CHARACTER TEXT12*12,HSIGN*12,TEXT16*16 + INTEGER OPTPRI(NSTATE) + DOUBLE PRECISION OPTPRR(NSTATE) + TYPE(C_PTR) IPOPT + INTEGER I,NVAR,NCST,LENGT,IPRINT,NSTPEX,IMTHD,M0,MINMAX,CNVTST, + 1 IERR + DOUBLE PRECISION DFLOTT,XDROIT,XS,XXS,EPS1,EPS4,EPSIM,ECOST, + 1 DELTA,SR,NORM,EPSEXT,COST,CQUAD,OBJNEW,OBJOLD, + 2 DERR,NORX,ERRX,DDX + LOGICAL LSAVE,LNORM2,LWAR,LBACK +*---- +* ALLOCATABLE ARRAYS +*---- + INTEGER, ALLOCATABLE, DIMENSION(:) :: INEGAL + DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: VARVAL,VARWGT, + > FCSTV,GRAD,ODX,ODF,DX,CONTR,VALMAX,VALMIN,DINF,DSUP,VARVL2, + > GRAD0,VARV0,WEIGH,DERIV0,CSTV0 +*---- +* CHECK THE VALIDITY OF OBJECTS +*---- + IF(NENTRY.NE.1) CALL XABORT('PLQ:ONE OBJECT EXPECTED.') + IF(JENTRY(1).NE.1) CALL XABORT('PLQ: OBJECT IN MODIFICATION ' + 1 //'MODE ONLY') + IPOPT=KENTRY(1) + IF((IENTRY(1).NE.1).AND.(IENTRY(1).NE.2))CALL XABORT('PLQ:' + 1 //' LCM OBJECT EXPECTED') + CALL LCMGTC(KENTRY(1),'SIGNATURE',12,HSIGN) + IF(HSIGN.NE.'L_OPTIMIZE') THEN + TEXT12=HENTRY(1) + CALL XABORT('PLQ: SIGNATURE OF '//TEXT12//' IS '//HSIGN// + 1 '. L_OPTIMIZE EXPECTED.') + ENDIF +*---- +* RECOVER STATE VECTOR INFORMATION +*---- + CALL LCMGET(IPOPT,'STATE-VECTOR',OPTPRI) + NVAR =OPTPRI(1) + NCST =OPTPRI(2) + MINMAX=OPTPRI(3) + ICONV =OPTPRI(4) + IF((MINMAX.NE.1).AND.(MINMAX.NE.-1)) CALL XABORT('PLQ: ' + 1 //'MINMAX not equal to 1 or -1') + NSTPEX=OPTPRI(5)+1 + IEDSTP=OPTPRI(6) + IMTHD =OPTPRI(9) + ISTEP= OPTPRI(10) + CALL LCMGET(IPOPT,'OPT-PARAM-R',OPTPRR) + SR =OPTPRR(1) + EPS1 =OPTPRR(2) + EPSEXT=OPTPRR(3) + EPSIM =OPTPRR(4) + EPS4 =OPTPRR(5) + ECOST =OPTPRR(6) +*---- +* SET CONTROL-VARIABLE VALUES +*---- + ALLOCATE(VARVAL(NVAR)) + CALL LCMLEN(IPOPT,'VAR-VALUE',LENGT,ITYP) + IF(LENGT.NE.NVAR) CALL XABORT('PLQ: WRONG NUMBER OF VARIABLE') +*---- +* SET CONTROL-VARIABLE WEIGHTS +*---- + ALLOCATE(VARWGT(NVAR)) + CALL LCMLEN(IPOPT,'VAR-WEIGHT',LENGT,ITYP) + IF(LENGT.EQ.0) THEN + VARWGT(:NVAR)=1.0D0 + ELSE IF(LENGT.EQ.NVAR) THEN + CALL LCMGET(IPOPT,'VAR-WEIGHT',VARWGT) + ELSE + CALL XABORT('PQL: NVAR - LENGT ARE NOT THE SAME') + ENDIF +*---- +* MEMORY ALLOCATION +*---- + ALLOCATE(FCSTV(NCST+1),GRAD(NVAR*(NCST+1)),ODX(NVAR),ODF(NVAR)) +*---- +* SET SYSTEM CHARACTERISTICS (THE OBJECTIVE FUNCTION IS THE FIRST ONE) +*---- + CALL LCMLEN(IPOPT,'FOBJ-CST-VAL',LENGT,ITYP) + IF(LENGT.EQ.0) CALL XABORT('PLQ: OBJECTIVE FUNCTION AND CONSTRA' + 1 //'INTS NOT YET EVALUATED') + CALL LCMGET(IPOPT,'FOBJ-CST-VAL',FCSTV) + COST=FCSTV(1) +*---- +* READ USER INPUT: +*---- + IPRINT=0 + LWAR=.FALSE. + 20 CALL REDGET(ITYP,NITMA,FLOTT,TEXT12,DFLOTT) +* Edition level + 30 IF(ITYP.NE.3) CALL XABORT('PLQ: CHARACTER DATA EXPECTED(1)') + IF(TEXT12.EQ.'EDIT') THEN + CALL REDGET(ITYP,IPRINT,FLOTT,TEXT12,DFLOTT) + IF(ITYP.NE.1) CALL XABORT('PLQ: *IPRINT* MUST BE INTEGER') + ELSE IF(TEXT12.EQ.'MINIMIZE') THEN + MINMAX=1 + ELSE IF(TEXT12.EQ.'MAXIMIZE') THEN + MINMAX=-1 + ELSE IF(TEXT12.EQ.'METHOD') THEN + CALL REDGET(ITYP,NITMA,FLOTT,TEXT12,DFLOTT) + IF(ITYP.NE.3) CALL XABORT('PLQ: CHARACTER DATA EXPECTED(2)') + IF(TEXT12.EQ.'SIMPLEX') THEN + IMTHD=1 + ELSE IF(TEXT12.EQ.'LEMKE') THEN + IMTHD=2 + ELSE IF(TEXT12.EQ.'MAP') THEN + IMTHD=3 + ELSE IF(TEXT12.EQ.'AUG-LAGRANG') THEN + IMTHD=4 + ELSE IF(TEXT12.EQ.'PENAL-METH') THEN + IMTHD=5 + ELSE + CALL XABORT('PLQ: WRONG METHOD KEYWORD') + ENDIF + ELSE IF(TEXT12.EQ.'OUT-STEP-LIM') THEN + CALL REDGET(ITYP,NITMA,FLOTT,TEXT12,DFLOTT) + IF(ITYP.EQ.2) THEN + SR=FLOTT + ELSE IF(ITYP.EQ.4) THEN + SR=DFLOTT + ELSE + CALL XABORT('PLQ: REAL OR DOUBLE PRECISION VALUE EXPECTED.') + ENDIF + ELSE IF(TEXT12.EQ.'INN-STEP-EPS') THEN +* Set the tolerence used for inner linear LEMKE or SIMPLEX +* calculation. + CALL REDGET(ITYP,NITMA,FLOTT,TEXT12,DFLOTT) + IF(ITYP.EQ.2) THEN + EPSIM=FLOTT + ELSE IF(ITYP.EQ.4) THEN + EPSIM=DFLOTT + ELSE + CALL XABORT('PLQ: REAL OR DOUBLE PRECISION VALUE EXPECTED.') + ENDIF + ELSE IF(TEXT12.EQ.'OUT-STEP-EPS') THEN +* Set the tolerence used for external iterations. + CALL REDGET(ITYP,NITMA,FLOTT,TEXT12,DFLOTT) + IF(ITYP.EQ.2) THEN + EPSEXT=FLOTT + ELSE IF(ITYP.EQ.4) THEN + EPSEXT=DFLOTT + ELSE + CALL XABORT('PLQ: REAL OR DOUBLE PRECISION VALUE EXPECTED.') + ENDIF + ELSE IF(TEXT12.EQ.'CST-QUAD-EPS') THEN + CALL REDGET(ITYP,NITMA,FLOTT,TEXT12,DFLOTT) + IF(ITYP.NE.2) CALL XABORT('PLQ: REAL DATA EXPECTED.') + EPS4=FLOTT + ELSE IF(TEXT12.EQ.'STEP-REDUCT') THEN + CALL REDGET(ITYP,NITMA,FLOTT,TEXT12,DFLOTT) + IF(ITYP.NE.3) CALL XABORT('PLQ: CHARACTER DATA EXPECTED(3).') + IF(TEXT12.EQ.'HALF') THEN + IEDSTP=1 + ELSE IF(TEXT12.EQ.'PARABOLIC') THEN + IEDSTP=2 + ELSE + CALL XABORT('PLQ: WRONG STEP REDUCTION KEYWORD.') + ENDIF + ELSE IF(TEXT12.EQ.'WARNING-ONLY')THEN +* Warning Only for failure of recovery of a valid point + LWAR=.TRUE. + ELSE IF(TEXT12.EQ.'CALCUL-DX')THEN +* Calculation of next point + GO TO 100 + ELSE IF(TEXT12.EQ.'COST-EXTRAP')THEN +* Cost extrapolation + GO TO 200 + ELSE IF(TEXT12.EQ.'OUT-CONV-TST') THEN +* Convergence test + GO TO 300 + ELSE IF( TEXT12.EQ.';' )THEN +* End of this subroutine + GO TO 1000 + ELSE + CALL XABORT('PLQ: '//TEXT12//' IS AN INVALID KEYWORD.') + ENDIF + GO TO 20 +*---- +* TEST FOR IMPROVEMENT FOR THE OBJECTIVE FUNCTION +*---- + 100 LBACK=.FALSE. + CALL LCMLEN(IPOPT,'OLD-VALUE',LENGT,ITYP) + IF((JENTRY(1).EQ.1).AND.(LENGT.NE.0)) THEN + ALLOCATE(CSTV0(NCST+1)) + OBJNEW=FCSTV(1) + CALL LCMSIX(IPOPT,'OLD-VALUE',1) + CALL LCMLEN(IPOPT,'FOBJ-CST-VAL',LENGT1,ITYP) + CALL LCMLEN(IPOPT,'VAR-VALUE',LENGT2,ITYP) + IF(LENGT1.EQ.0) THEN + CALL XABORT('PLQ: MISSING OLD OBJECTIVE FUNCTION VALUE') + ELSE IF(LENGT1.NE.NCST+1) THEN + CALL XABORT('PLQ: WRONG NUMBER OF CONSTRAINTS') + ELSE IF(LENGT2.EQ.0) THEN + CALL XABORT('PLQ: MISSING CONTROL VARIABLES RECORD') + ELSE IF(LENGT2.NE.NVAR) THEN + CALL XABORT('PLQ: WRONG NUMBER OF CONTROL VARIABLES') + ENDIF + CALL LCMGET(IPOPT,'FOBJ-CST-VAL',CSTV0) + OBJOLD=CSTV0(1) + IF(OBJNEW.GE.OBJOLD) THEN + LBACK=.TRUE. + IF(IPRINT.GT.1) WRITE(6,4005) OBJOLD,OBJNEW + ENDIF + DEALLOCATE(CSTV0) + CALL LCMSIX(IPOPT,' ',2) + ENDIF +*---- +* RECOVER OBJECTIVE FUNCTION AND GRADIENTS FROM PRECEDING ITERATION +*---- + IF(LBACK) THEN + ISTEP=0 + CALL LCMGET(IPOPT,'VAR-VALUE',VARVAL) + IF(IPRINT.GT.1) THEN + WRITE(6,4001) 'REJECTED CONTROL VARIABLES:', + 1 (VARVAL(IVAR),IVAR=1,NVAR) + ENDIF + CALL LCMSIX(IPOPT,'OLD-VALUE',1) + ALLOCATE(CSTV0(NCST+1),VARV0(NVAR),DERIV0(NVAR*(NCST+1)), + 1 WEIGH(NVAR)) + CALL LCMGET(IPOPT,'FOBJ-CST-VAL',CSTV0) + CALL LCMGET(IPOPT,'VAR-VALUE',VARV0) + CALL LCMGET(IPOPT,'GRADIENT',DERIV0) + CALL LCMSIX(IPOPT,' ',2) + CALL LCMPUT(IPOPT,'FOBJ-CST-VAL',NCST+1,4,CSTV0) + CALL LCMPUT(IPOPT,'VAR-VALUE',NVAR,4,VARV0) + CALL LCMPUT(IPOPT,'GRADIENT',NVAR*(NCST+1),4,DERIV0) + IF(IEDSTP.LE.1) THEN + SR=SR*0.5 + ELSE IF(IEDSTP.EQ.2) THEN + CALL LCMLEN(IPOPT,'VAR-WEIGHT',LENGT,ITYLCM) + IF(LENGT.EQ.NVAR) THEN + CALL LCMGET(IPOPT,'VAR-WEIGHT',WEIGH) + ELSE + WEIGH(:NVAR)=1.0D0 + ENDIF + NORX=0.0D0 + DERR=0.0D0 + DO 110 I=1,NVAR + DDX=VARVAL(I)-VARV0(I) + NORX=NORX+WEIGH(I)*DDX**2 + DERR=DERR+SQRT(WEIGH(I))*DDX*DERIV0(I) + 110 CONTINUE + NORX=NORX**0.5 + DERR=DERR/NORX + ERRX=ABS(0.5*DERR*NORX*NORX/(DERR*NORX-(OBJNEW-OBJOLD))) + SR=MAX(MIN(SR,ERRX),SR/20.0) + DEALLOCATE(WEIGH) + ENDIF + IF(IPRINT.GT.1) WRITE(6,'(/31H PLQ: REDUCES QUADRATIC CONSTRA, + 1 13HINT RADIUS TO,1P,E11.4,8H IEDSTP=,I4)') SR,IEDSTP + IF(SR.LE.EPS4) THEN + WRITE(6,4006) + ICONV=1 + ENDIF + DEALLOCATE(DERIV0,VARV0,CSTV0) +*---- +* USES NEW GRADIENTS FROM MODULE GRAD: +*---- + ELSE +* count the number of iterations without step back + ISTEP=ISTEP+1 + IF(ISTEP.GT.10) THEN + SR=2.0*SR + ISTEP=5 + IF(IPRINT.GT.1) WRITE(6,'(/29H PLQ: INCREASES QUADRATIC CON, + 1 17HSTRAINT RADIUS TO,1P,E11.4)') SR + ENDIF + CALL LCMGET(IPOPT,'VAR-VALUE',VARVAL) + CALL LCMSIX(IPOPT,'OLD-VALUE',1) + CALL LCMPUT(IPOPT,'VAR-VALUE2',NVAR,4,VARVAL) + CALL LCMSIX(IPOPT,' ',2) + ENDIF +*---- +* SET GRADIENTS +*---- + CALL LCMGET(IPOPT,'GRADIENT',GRAD) +*---- +* PRINT INFORMATION +*---- + IF(IPRINT.GT.0) THEN + WRITE(6,'(/47H PLQ: INFORMATION AT QUADRATIC CONSTRAINT ITERA, + 1 4HTION,I5)') NSTPEX + WRITE(6,3999) NSTPEX,FCSTV(1) + WRITE(6,4000) 'QUADRATIC CONSTRAINT RADIUS:',SR + IF(NCST.GT.0) WRITE(6,4001) 'CONSTRAINTS:',(FCSTV(ICST), + 1 ICST=2,NCST+1) + CALL LCMLEN(IPOPT,'VAR-VALUE',LENGT1,ITYLCM) + IF(LENGT1.GT.0) THEN + CALL LCMGET(IPOPT,'VAR-VALUE',VARVAL) + WRITE(6,4001) 'CONTROL VARIABLES:',(VARVAL(IVAR),IVAR=1,NVAR) + ENDIF + IF(IPRINT.GT.1) THEN + ALLOCATE(DERIV0(NVAR*(NCST+1))) + CALL LCMGET(IPOPT,'GRADIENT',DERIV0) + WRITE(6,'(/29H GRADIENTS-------------------)') + WRITE(6,4001) 'OBJECTIVE FUNCTION:',(DERIV0(IVAR),IVAR=1,NVAR) + IF(IPRINT.GT.2) THEN + DO 120 ICST=1,NCST + WRITE(TEXT16,'(10HCONSTRAINT,I4,1H:)') ICST + WRITE(6,4001) TEXT16,(DERIV0(ICST*NVAR+IVAR),IVAR=1,NVAR) + 120 CONTINUE + ENDIF + DEALLOCATE(DERIV0) + ENDIF + IF(LBACK) WRITE(6,'(28H *** STEP BACK ITERATION ***)') + ENDIF +*---- +* NEXT STEP CALCULATION +*---- + CALL LCMGET(IPOPT,'VAR-VALUE',VARVAL) + CALL REDGET(ITYP,NITMA,FLOTT,TEXT12,DFLOTT) + IF(TEXT12.EQ.'NO-STORE-OLD') THEN + LSAVE=.TRUE. + CALL REDGET(ITYP,NITMA,FLOTT,TEXT12,DFLOTT) + ELSE + LSAVE=.FALSE. + ENDIF + ITYP1=ITYP + ALLOCATE(DX(NVAR)) + IF(NCST.GT.0) THEN +* INEQUAL + ALLOCATE(INEGAL(NCST)) + CALL LCMLEN(IPOPT,'CST-TYPE',LENGT,ITYP) + IF(LENGT.NE.NCST) CALL XABORT('PLQ: NCST - LENGT NOT EQUAL') + CALL LCMGET(IPOPT,'CST-TYPE',INEGAL) +* +* CONTR + ALLOCATE(CONTR(NCST)) + CALL LCMLEN(IPOPT,'CST-OBJ',LENGT,ITYP) + IF(LENGT.NE.NCST) CALL XABORT('PLQ: NCST - LENGT NOT EQUAL') + CALL LCMGET(IPOPT,'CST-OBJ',CONTR) + DO 130 I=1,NCST + CONTR(I) = CONTR(I)-FCSTV(I+1) + 130 CONTINUE + ENDIF +* +* DINF AND DSUP + CALL LCMLEN(IPOPT,'VAR-VAL-MAX',LENGT,ITYP) + IF(LENGT.EQ.0) CALL XABORT('PLQ: NO MAXIMUM VALUE DEFINED') + ALLOCATE(VALMAX(NVAR),VALMIN(NVAR)) + CALL LCMGET(IPOPT,'VAR-VAL-MAX',VALMAX) + CALL LCMLEN(IPOPT,'VAR-VAL-MIN',LENGT,ITYP) + IF(LENGT.EQ.0) CALL XABORT('PLQ: NO MAXIMUM VALUE DEFINED') + CALL LCMGET(IPOPT,'VAR-VAL-MIN',VALMIN) + ALLOCATE(DINF(NVAR),DSUP(NVAR)) + DO 140 I=1,NVAR + DINF(I) = VALMIN(I) - VARVAL(I) + DSUP(I) = VALMAX(I) - VARVAL(I) + 140 CONTINUE + DEALLOCATE(VALMAX,VALMIN) +* + M0 = NCST + XDROIT = SR**2 + IF(IPRINT.GE.1) WRITE(6,4002) XDROIT,(VARWGT(I),I=1,NVAR) +*---- +* FIND ACTIVE CONSTRAINTS FOR XK(I) LIMITS +*---- + DO 150 I=1,NVAR + XS = SQRT(XDROIT/VARWGT(I)) + XXS=-XS + IF(DINF(I).GT.XXS) THEN + M0 = M0 + 1 + ENDIF + IF(DSUP(I).LT.XS) THEN + M0 = M0 + 1 + ENDIF + 150 CONTINUE +*---- +* SOLUTION OF A LINEAR OPTIMIZATION PROBLEM WITH A QUADRATIC CONSTRAINT +*---- + IERR=0 + CALL PLDRV(IPOPT,NVAR,NCST,M0,MINMAX,IMTHD,COST,DX,VARWGT,GRAD, + > INEGAL,CONTR,DINF,DSUP,XDROIT,EPSIM,IPRINT,IERR) +*---- +* STEP-BACK IN CASE OF FAILURE +*---- + IF(IERR.GE.1) THEN + OPTPRI(14)=OPTPRI(14)+1 + CALL LCMSIX(IPOPT,'OLD-VALUE',1) + CALL LCMLEN(IPOPT,'VAR-VALUE2',LENGT,ITYP) + IF(LENGT.EQ.0) THEN + IF(LWAR) THEN + WRITE(6,*) 'WARNING: UNABLE TO RECOVER A VALID POINT' + 1 //' WITH SUCCESSFUL "PLQ" RESOLUTION' + ELSE + CALL LCMLIB(IPOPT) + CALL XABORT('PLQ: UNABLE TO RECOVER A VALID POINT WITH ' + 1 //'SUCCESSFUL "PLQ" RESOLUTION') + ENDIF + ELSE + ALLOCATE(VARVL2(NVAR)) + CALL LCMGET(IPOPT,'VAR-VALUE2',VARVL2) + DO 160 I=1,NVAR + DX(I)=(VARVL2(I)-VARVAL(I))/2.0 + 160 CONTINUE + DEALLOCATE(VARVL2) + ENDIF + CALL LCMSIX(IPOPT,' ',2) + IF(IPRINT.GE.1) WRITE(6,*) 'IERR>0' + IF(IPRINT.GE.1) WRITE(6,*) 'DX=',(DX(I),I=1,NVAR) + ELSE + OPTPRI(14)=0 + ENDIF +* + DO 170 I=1,NVAR + ODX(I)=DX(I) + ODF(I)=GRAD(I) + 170 CONTINUE + DEALLOCATE(DX) + IF(NCST.GT.0) DEALLOCATE(INEGAL) + DEALLOCATE(DINF,DSUP) + IF(NCST.GT.0) DEALLOCATE(CONTR) +*---- +* BACKUP VALUES OF THE PRECEDING ITERATION +*---- + IF(.NOT.LSAVE) THEN + CALL LCMSIX(IPOPT,'OLD-VALUE',1) + CALL LCMPUT(IPOPT,'VAR-VALUE',NVAR,4,VARVAL) + CALL LCMPUT(IPOPT,'FOBJ-CST-VAL',NCST+1,4,FCSTV) + CALL LCMPUT(IPOPT,'GRADIENT',NVAR*(NCST+1),4,GRAD) + CALL LCMSIX(IPOPT,' ',2) + ENDIF +*---- +* BACKUP VALUES OF THE NEW ITERATION +*---- + DO 180 I=1,NVAR + VARVAL(I)=VARVAL(I)+ODX(I) + 180 CONTINUE + CALL LCMPUT(IPOPT,'VAR-VALUE',NVAR,4,VARVAL) + ITYP=ITYP1 +*---- +* EXTRAPOLATE OBJECTIVE FUNCTION +*---- + ECOST=COST + DO 190 I=1,NVAR + ECOST=ECOST+ODX(I)*ODF(I) + 190 CONTINUE +*---- +* REINITIALIZE GRADIENTS FOR THE NEXT ITERATION +*---- + ALLOCATE(GRAD0(NVAR*(NCST+1))) + GRAD0(:NVAR*(NCST+1))=0.0D0 + CALL LCMPUT(IPOPT,'GRADIENT',NVAR*(NCST+1),4,GRAD0) + DEALLOCATE(GRAD0) + GO TO 30 +*---- +* OUTPUT THE EXTRAPOLATED OBJECTIVE FUNCTION +*---- + 200 ECOSTR=REAL(ECOST) + CALL REDGET(ITYP,NITMA,ECOSTR,TEXT12,DFLOTT) + IF(ITYP.NE.-2) CALL XABORT('PLQ: OUTPUT REAL EXPECTED') + ITYP=2 + CALL REDPUT(ITYP,NITMA,ECOSTR,TEXT12,DFLOTT) + GO TO 20 +*---- +* TEST CONVERGENCE +*---- + 300 LNORM2=.TRUE. + CALL REDGET(ITYP,CNVTST,FLOTT,TEXT12,DFLOTT) + IF((ITYP.EQ.3).AND.(TEXT12.EQ.'NORM-INF')) THEN + LNORM2=.FALSE. + CALL REDGET(ITYP,CNVTST,FLOTT,TEXT12,DFLOTT) + ENDIF + IF(ITYP.NE.-5) CALL XABORT('PLQ: OUTPUT LOGICAL EXPECTED') + DELTA=ABS((ECOST-COST)/COST) + NORM=0.0 + CQUAD=0.0 + IF(LNORM2) THEN + DO 350 I=1,NVAR + NORM=NORM+VARWGT(I)*VARVAL(I)*VARVAL(I) + CQUAD=CQUAD+VARWGT(I)*ODX(I)*ODX(I) + 350 CONTINUE + IF(NORM.NE.0.0) THEN + CQUAD=SQRT(CQUAD/NORM) + ELSE + CQUAD=0.0 + ENDIF + ELSE + DO 360 I=1,NVAR + NORM=MAX(NORM,ABS(VARWGT(I)**0.5*VARVAL(I))) + CQUAD=MAX(CQUAD,ABS(VARWGT(I)**0.5*ODX(I))) + 360 CONTINUE + IF(NORM.NE.0.0) THEN + CQUAD=CQUAD/NORM + ELSE + CQUAD=0.0 + ENDIF + ENDIF + IF(EPSEXT.EQ.0.0) EPSEXT = 0.001D0 + IF(((DELTA.LT.EPSEXT).AND.(CQUAD.LE.EPSEXT)) .OR. + 1 (CQUAD.LE.(EPSEXT/10.0))) THEN + CNVTST=1 + ICONV =1 + ELSE + CNVTST=-1 + ICONV =0 + ENDIF + IF(IPRINT.GE.1) THEN + WRITE(6,*) 'It= convergence?', DELTA,CQUAD,EPSEXT + IF(IPRINT.GT.2) THEN + WRITE(6,*) 'DX',(ODX(I),I=1,NVAR) + WRITE(6,*) 'X',(VARVAL(I),I=1,NVAR) + ENDIF + ENDIF + ITYP=5 + CALL REDPUT(ITYP,CNVTST,FLOTT,TEXT12,DFLOTT) + GO TO 20 +*---- +* END +*---- + 1000 DEALLOCATE(VARWGT,FCSTV,GRAD,ODX,ODF,VARVAL) +*---- +* SAVE THE STATE VECTORS +*---- + OPTPRI(:NSTATE)=0 + OPTPRI(1)=NVAR + OPTPRI(2)=NCST + OPTPRI(3)=MINMAX + OPTPRI(4)=ICONV + OPTPRI(5)=NSTPEX + OPTPRI(6)=IEDSTP + OPTPRI(7)=0 + OPTPRI(8)=1 + OPTPRI(9)=IMTHD + OPTPRI(10)=ISTEP + IF(IPRINT.GT.0) WRITE(6,4003) (OPTPRI(I),I=1,10) + CALL LCMPUT(IPOPT,'STATE-VECTOR',NSTATE,1,OPTPRI) + OPTPRR(:NSTATE)=0.0D0 + OPTPRR(1)=SR + OPTPRR(2)=EPS1 + OPTPRR(3)=EPSEXT + OPTPRR(4)=EPSIM + OPTPRR(5)=EPS4 + OPTPRR(6)=ECOST + IF(IPRINT.GT.0) WRITE(6,4004) (OPTPRR(I),I=1,6) + CALL LCMPUT(IPOPT,'OPT-PARAM-R',NSTATE,4,OPTPRR) + IF(IPRINT.GT.1) CALL LCMLIB(IPOPT) + RETURN +* + 3999 FORMAT(/13H PLQ: ##ITER=,I8,20H OBJECTIVE FUNCTION=,1P,E14.6) + 4000 FORMAT(1X,A28,1P,E14.6) + 4001 FORMAT(1X,A28,1P,8E12.4/(29X,8E12.4)) + 4002 FORMAT(//,5X,'SR**2 (XDROIT) = ',1P,D13.5, + > /,5X,'FPOIDS = ',/,(11X,1P,8D13.5)) + 4003 FORMAT(/8H OPTIONS/8H -------/ + 1 7H NVAR ,I8,32H (NUMBER OF CONTROL VARIABLES)/ + 2 7H NCST ,I8,26H (NUMBER OF CONSTRAINTS)/ + 3 7H MINMAX,I8,37H (=1/-1: MINIMIZATION/MAXIMIZATION)/ + 4 7H ICONV ,I8,43H (=0/1: EXTERNAL NOT CONVERGED/CONVERGED)/ + 5 7H NSTPEX,I8,44H (ITERATION INDEX OF QUADRATIC CONSTRAINT)/ + 6 7H IEDSTP,I8,43H (=1/2: HALF REDUCTION/PARABOLIC FORMULA)/ + 7 7H IHESS ,I8,29H (=0/1/2: STEEPEST/CG/BFGS)/ + 8 7H ISEARC,I8,35H (=0/1/2: NO SEARCH/OPTEX/NEWTON)/ + 9 7H IMTHD ,I8,42H (=1/2/3: SIMPLEX-LEMKE/LEMKE-LEMKE/MAP)/ + 1 7H ISTEP ,I8,43H (NUMBER OF ITERATIONS WITHOUT STEP-BACK)) + 4004 FORMAT(/ + 1 12H REAL PARAM:,1P/12H -----------/ + 2 7H SR ,D12.4,39H (RADIUS OF THE QUADRATIC CONSTRAINT)/ + 3 7H EPS1 ,D12.4,13H (NOT USED)/ + 4 7H EPSEXT,D12.4,31H (EXTERNAL CONVERGENCE LIMIT)/ + 5 7H EPSIM ,D12.4,31H (INTERNAL CONVERGENCE LIMIT)/ + 6 7H EPS4 ,D12.4,43H (QUADRATIC CONSTRAINT CONVERGENCE LIMIT)/ + 7 7H ECOST ,D12.4,17H (UPDATED COST)) + 4005 FORMAT(/38H PLQ: OBJECTIVE FUNCTION INCREASE FROM,1P,E12.4, + 1 3H TO,E12.4/35H RETURN BACK TO PREVIOUS ITERATION.) + 4006 FORMAT(/1X,'PLQ: THE QUADRATIC CONSTRAINT RADIUS CANNOT BE FUR', + 1 'THER REDUCED') + END |