diff options
Diffstat (limited to 'Donjon/src/PLLA.f')
| -rw-r--r-- | Donjon/src/PLLA.f | 240 |
1 files changed, 240 insertions, 0 deletions
diff --git a/Donjon/src/PLLA.f b/Donjon/src/PLLA.f new file mode 100644 index 0000000..da3b884 --- /dev/null +++ b/Donjon/src/PLLA.f @@ -0,0 +1,240 @@ +*DECK PLLA + SUBROUTINE PLLA(IPOPT,N0,M0,APLUS,PDG,BPLUS,INPLUS,XDROIT,FCOST, + > GF,XOBJ,IMPR,EPSIM,NCST,GRAD,CONTR,INEGAL,IERR) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Solves the quasi-linear problem using the augmented Lagrangian. +* PLLA = Linear Programmation Augmented Lagrangian +* +*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/ouput +* IPOPT pointer to the L_OPTIMIZE object. +* N0 number of control variables. +* M0 number of constraints. +* APLUS coefficient matrix for the linear constraints. +* PDG weights assigned to control variables in the quadratic +* constraint. +* BPLUS right hand sides corresponding to the coefficient matrix. +* INPLUS constraint relations (=-1 for .GE.; =0 for .EQ.; =1 for .LE.). +* XDROIT quadratic constraint radius squared. +* FCOST costs of control variables. +* GF objective function. +* XOBJ control variables. +* IMPR print flag. +* EPSIM tolerence used for inner linear SIMPLEX calculation. +* NCST number of constraints. +* GRAD linearized gradients (GRAD(:,1) are control variable costs +* and GRAD(:,2:NCST+1) are linear constraint coefficients). +* CONTR constraint right hand sides. +* INEGAL constraint relations (=-1 for .GE.; =0 for .EQ.; =1 for .LE.). +* +*Parameters: ouput +* IERR return code (=0: normal completion). +* +*----------------------------------------------------------------------- +* + USE GANLIB + IMPLICIT NONE +*---- +* SUBROUTINE ARGUMENTS +*---- + TYPE(C_PTR) IPOPT + INTEGER N0,M0,INPLUS(M0+1),IMPR,NCST,INEGAL(NCST),IERR + DOUBLE PRECISION PDG(N0),BPLUS(M0+2),XDROIT,XOBJ(N0),EPSIM, + > GRAD(N0,NCST+1),CONTR(NCST),APLUS(M0+2,M0+N0+1),GF(N0),FCOST +*---- +* LOCAL VARIABLES +*---- + INTEGER ITERMX,LENGT,ITYP,I,J,K,ITER,NPM,M0B + PARAMETER (ITERMX=10) + LOGICAL LCST(NCST),LCST2(NCST),LTST + DOUBLE PRECISION LACOST,NORM,CRIT,LA0E + INTEGER, ALLOCATABLE, DIMENSION(:) :: INPL2 + DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: LAFAC,CSTWGT,B2, + > CONTR2,LAGF + DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: APLUS2 +*---- +* SCRATCH STORAGE ALLOCATION +*---- + ALLOCATE(INPL2(M0-NCST+1)) + ALLOCATE(LAFAC(NCST),CSTWGT(NCST),B2(M0-NCST+2),CONTR2(NCST)) + ALLOCATE(LAGF(N0),APLUS2(M0-NCST+1,N0+M0-NCST)) +*---- +* STEP 0: INITIALIZATION +* NPM: SIZE OF THE LINEARIZED SYSTEM. +* M0B: NUMBER OF LINEARIZED CONSTRAINTS FOR THE LA ALGORITHM. +* CORRESPONDS TO THE NUMBER OF POSSIBLY ACTIVE BOUNDS. +*---- + NPM=(M0+1)+N0 + M0B=M0-NCST + IF(NCST.GT.0) THEN + CALL LCMLEN(IPOPT,'F-LA-MULT',LENGT,ITYP) + IF(LENGT.EQ.0) THEN + CALL XABORT('PLLA: LAGRANGIAN FACTORS NON INITIALIZED') + ELSEIF(LENGT.EQ.NCST) THEN + CALL LCMGET(IPOPT,'F-LA-MULT',LAFAC) + ELSE + CALL XABORT('PLLA: WRONG NUMBER OF LA COEFFICIENTS') + ENDIF + CALL LCMLEN(IPOPT,'CST-WEIGHT',LENGT,ITYP) + IF(LENGT.EQ.0) THEN + CALL XABORT('PLLA: CONSTRAINTS PENALTIES NON INITIALIZED') + ELSEIF(LENGT.EQ.NCST) THEN + CALL LCMGET(IPOPT,'CST-WEIGHT',CSTWGT) + ELSE + CALL XABORT('PLLA: WRONG NUMBER OF CONSTRAINT WEIGHTS') + ENDIF + DO 10 J=1,NCST + CONTR2(J)=-CONTR(J) + 10 CONTINUE + LCST2(:NCST)=.TRUE. + LCST(:NCST)=.TRUE. + ENDIF + XOBJ(:N0)=0.0D0 +*---- +* INTERNAL ITERATIONS FOR THE LINEAR PROBLEM +*---- + ITER=0 + 99 ITER=ITER+1 + LTST=.TRUE. +*---- +* STEP 1: LA FUNCTION EVALUATION +*---- + DO 110 J=1,NCST + IF(INEGAL(J).NE.0) THEN + CRIT=CONTR2(J) + DO 100 I=1,N0 + CRIT=CRIT+GRAD(I,J+1)*XOBJ(I) + 100 CONTINUE + CRIT=CSTWGT(J)*INEGAL(J)*CRIT+LAFAC(J) + LCST(J)=(CRIT.LE.0.0) + ENDIF + 110 CONTINUE + + DO 150 I=1,N0 + LAGF(I)=GF(I) + DO 140 J=1,NCST + IF(INEGAL(J).EQ.0) THEN + LAGF(I)=LAGF(I)+GRAD(I,J+1)*(LAFAC(J)+CSTWGT(J)*CONTR2(J)) + ELSEIF(.NOT.LCST(J)) THEN + LAGF(I)=LAGF(I)+INEGAL(J)*GRAD(I,J+1) + 1 *(CSTWGT(J)*INEGAL(J)*CONTR2(J)+2*LAFAC(J)) + ENDIF + 140 CONTINUE + 150 CONTINUE + + LACOST=FCOST + DO 160 J=1,NCST + IF(INEGAL(J).EQ.0) THEN + LACOST=LACOST+(LAFAC(J)+CSTWGT(J)/2.0*CONTR2(J))*CONTR2(J) + ELSEIF(.NOT.LCST(J)) THEN + LACOST=LACOST+INEGAL(J)*CONTR2(J)*2.0*LAFAC(J) + 1 +CSTWGT(J)/2.0*CONTR2(J)**2 + ELSE + LACOST=LACOST-3.0*LAFAC(J)**2/2.0/CSTWGT(J) + ENDIF + 160 CONTINUE + IF(ITER.EQ.1) LA0E=LACOST + IF(IMPR.GE.3) THEN + WRITE(6,*) 'GF',(GF(I),I=1,N0) + WRITE(6,*) 'LAGF',(LAGF(I),I=1,N0) + WRITE(6,*) 'PDG',(PDG(I),I=1,N0) + WRITE(6,*) 'LACOST',LACOST,'M0B',M0B,'XDROIT',XDROIT + ENDIF +*---- +* STEP 2: SOLUTION +* k,l +* compute DX +* case 1 +* If there is no constraints for the LA problem (M0B=0), +* then the solution is obvious: on the hypersphere(radius XDROIT) +* in the direction LAGF +*---- + IF(M0B.EQ.0) THEN +* + NORM=0.0 + DO 200 I=1,N0 + NORM=NORM+LAGF(I)**2/PDG(I) + 200 CONTINUE + NORM=NORM**0.5 +* + IF(NORM.EQ.0.0) THEN + XOBJ(:N0)=0.0D0 + ELSE + DO 210 I=1,N0 + XOBJ(I)=-XDROIT**0.5*LAGF(I)/PDG(I)/NORM + 210 CONTINUE + ENDIF +*---- +* CASE 2 +* SOLUTION WITH THE LEMKE METHOD +*---- + ELSE +* + DO 260 K=1,M0B + DO 250 I=1,N0 + APLUS2(K,I)=APLUS(NCST+K,I) + 250 CONTINUE + B2(K)=BPLUS(NCST+K) + INPL2(K)=INPLUS(NCST+K) + 260 CONTINUE + DO 270 I=1,N0 + APLUS2(M0B+1,I) = 0.0D0 + 270 CONTINUE + BPLUS(M0B+1) = 0.0 + INPL2(M0B+1) = 0 +* + CALL PLMAP2(N0,M0B,APLUS2,PDG,B2,INPL2,XDROIT,LAGF,LACOST,XOBJ,2, + > EPSIM,IMPR,IERR) +* + ENDIF + DO 410 J=1,NCST + IF(INEGAL(J).NE.0) THEN + CRIT=CONTR2(J) + DO 400 I=1,N0 + CRIT=CRIT+GRAD(I,J+1)*XOBJ(I) + 400 CONTINUE + CRIT=CSTWGT(J)*INEGAL(J)*CRIT+LAFAC(J) + LCST2(J)=(CRIT.LE.0.0) + ENDIF + 410 CONTINUE + + IF((IMPR.GE.2).AND.(NCST.GT.0)) THEN + WRITE(6,*) (LCST(J),J=1,NCST) + WRITE(6,*) (LCST2(J),J=1,NCST) + ENDIF + DO 420 J=1,NCST + LTST=LTST.AND.(LCST(J).EQV.LCST2(J)) + 420 CONTINUE + + IF((.NOT.LTST) .AND.(ITER.LE.ITERMX)) GO TO 99 +*---- +* k,l +* STEP 3: SAVE L +* a +*---- + CALL LCMSIX(IPOPT,'OLD-VALUE',1) + CALL LCMPUT(IPOPT,'F-LA-EVAL',1,4,LA0E) + IF(IMPR.GE.1) WRITE(6,*) 'LAGF',(LAGF(I),I=1,N0) + CALL LCMPUT(IPOPT,'DF-LA-PENAL',N0,4,LAGF) + CALL LCMSIX(IPOPT,' ',0) +* + IF(IMPR.GE.1) WRITE(6,*) 'Dvar',(XOBJ(I),I=1,N0) +*---- +* SCRATCH STORAGE DEALLOCATION +*---- + DEALLOCATE(APLUS2,LAGF) + DEALLOCATE(CONTR2,B2,CSTWGT,LAFAC) + DEALLOCATE(INPL2) + RETURN + END |