Initial commit from Polytechnique Montreal

1 files changed, 190 insertions, 0 deletions

diff --git a/Donjon/src/PLDRV.f b/Donjon/src/PLDRV.f
new file mode 100644
index 0000000..7b75d2a
--- /dev/null
+++ b/Donjon/src/PLDRV.f
@@ -0,0 +1,190 @@
+*DECK PLDRV
+      SUBROUTINE PLDRV(IPOPT,N0,NCST,M0,MINMAX,IMTHD,FCOST,XOBJ,PDG,
+     >           GRAD,INEGAL,CONTR,DINF,DSUP,XDROIT,EPSIM,IMPR,IERR)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Prepares the different matrices for the resolution of a linear
+* optimisation problem with a quadratic constraint. The different
+* available methods are the MAP technique, the lemke, the augmented-
+* Lagrangian and the penalty function.
+* PLDRV = Linear Programmation DRiVer for options
+*
+*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): 
+* A. Hebert and R. Chambon
+*
+*Parameters: input
+* IPOPT   pointer to the L_OPTIMIZE object.
+* N0      number of control variables.
+* NCST    number of constraints.
+* M0      number of constraints plus the number of lower/upper bounds
+*         intercepting the quadratic constraint.
+* MINMAX  type of optimization (=-1: minimize; =1: maximize).
+* IMTHD   type of solution (=1: SIMPLEX/LEMKE; =2: LEMKE/LEMKE;
+*         =3: MAP; =4: Augmented Lagragian; =5: External penalty
+*         funnction).
+* FCOST   objective function.
+* XOBJ    control variables.
+* PDG     weights assigned to control variables in the quadratic
+*         constraint.
+* GRAD    linearized gradients (GRAD(:,1) are control variable costs
+*         and GRAD(:,2:NCST+1) are linear constraint coefficients).
+* INEGAL  constraint relations (=-1 for .GE.; =0 for .EQ.; =1 for .LE.).
+* CONTR   constraint right hand sides.
+* DINF    lower bounds of control variables.
+* DSUP    upper bounds of control variables.
+* XDROIT  quadratic constraint radius squared.
+* EPSIM   tolerence used for inner linear LEMKE or SIMPLEX calculation.
+* IMPR    print flag.
+*
+*Parameters: ouput
+* IERR    return code (=0: normal completion).
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+      IMPLICIT NONE
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      TYPE(C_PTR) IPOPT
+      INTEGER N0,NCST,M0,MINMAX,IMTHD,INEGAL(NCST),IMPR,IERR
+      DOUBLE PRECISION FCOST,XOBJ(N0),PDG(N0),GRAD(N0,NCST+1),
+     >     CONTR(NCST),DINF(N0),DSUP(N0),XDROIT,EPSIM
+*----
+*  LOCAL VARIABLES
+*----
+      INTEGER      ME,MI,I,J,NPM
+      DOUBLE PRECISION  XX
+      CHARACTER    CLNAME*6
+      INTEGER, ALLOCATABLE, DIMENSION(:) :: INPLUS
+      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: BPLUS,GF
+      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: APLUS
+*----
+*  DATA STATEMENTS
+*----
+      DATA         CLNAME /'PLDRV'/
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(INPLUS(M0+1))
+      ALLOCATE(BPLUS(M0+2),GF(N0),APLUS(M0+2,(M0+1)+N0))
+*----
+*  SET COSTS OF CONTROL VARIABLES
+*----
+      DO 10 I=1,N0
+         GF(I)=GRAD(I,1)*REAL(MINMAX)
+   10 CONTINUE
+*----
+*  ORGANIZE TABLES APLUS AND BPLUS FOR EQUALITY CONSTRAINTS
+*----
+      ME=0
+      DO 30 I=1,NCST
+         IF(INEGAL(I).EQ.0) THEN
+            ME = ME + 1
+            DO 20 J=1,N0
+               APLUS(ME,J) = GRAD(J,I+1)
+   20       CONTINUE
+            BPLUS(ME)  = CONTR(I)
+            INPLUS(ME) = 0
+         ENDIF
+   30 CONTINUE
+*----
+*  ORGANIZE TABLES APLUS AND BPLUS FOR INEQUALITY CONSTRAINTS
+*----
+      MI=0
+      DO 50 I=1,NCST
+         IF(INEGAL(I).NE.0) THEN
+            MI = MI + 1
+            DO 40 J=1,N0
+               APLUS(ME+MI,J) = GRAD(J,I+1)
+   40       CONTINUE
+            BPLUS(ME+MI)  = CONTR(I)
+            INPLUS(ME+MI) = INEGAL(I)
+         ENDIF
+   50 CONTINUE
+*----
+*  ORGANIZE TABLES APLUS AND BPLUS FOR CONTROL-VARIABLE BOUNDS
+*----
+      DO 80 I=1,N0
+         XX = SQRT(XDROIT/PDG(I))
+         IF(DINF(I).GT.-XX) THEN
+            MI = MI + 1
+            DO 60 J=1,N0
+               APLUS(ME+MI,J) = 0.0D0
+   60       CONTINUE
+            APLUS(ME+MI,I) = 1.0D0
+            BPLUS(ME+MI)   = DINF(I)
+            INPLUS(ME+MI)  = -1
+         ENDIF
+         IF(DSUP(I).LT.XX) THEN
+            MI = MI + 1
+            DO 70 J=1,N0
+               APLUS(ME+MI,J) = 0.0D0
+   70       CONTINUE
+            APLUS(ME+MI,I) = 1.0D0
+            BPLUS(ME+MI)   = DSUP(I)
+            INPLUS(ME+MI)  = 1
+         ENDIF
+   80 CONTINUE
+*
+      DO 90 J=1,N0
+         APLUS(M0+1,J) = 0.0D0
+   90 CONTINUE
+      BPLUS(M0+1)   = 0.0D0
+      INPLUS(M0+1)  = 0
+*
+      IF(M0.NE.ME+MI) THEN
+         WRITE (6,1000) M0,ME,MI
+         CALL XABORT('PLDRV: M0 AND ME+MI ARE NOT THE SAME')
+      ENDIF
+*----
+*  PRINT THE QUASILINEAR PROBLEM
+*----
+      IF(IMPR.GE.5) THEN
+         CALL PLNTAB(GF,APLUS,INPLUS,BPLUS,PDG,DINF,DSUP,N0,M0,
+     >             CLNAME)
+      ENDIF
+*----
+*  LEMKE METHOD
+*----
+      NPM=(M0+1)+N0
+      IF(IMTHD.LE.2) THEN
+         CALL PLMAP2(N0,M0,APLUS,PDG,BPLUS,INPLUS,XDROIT,GF,FCOST,XOBJ,
+     1   IMTHD,EPSIM,IMPR,IERR)
+*----
+*  MAP
+*----
+      ELSE IF(IMTHD.EQ.3) THEN
+         CALL PLMAP1(N0,M0,APLUS,PDG,BPLUS,INPLUS,XDROIT,GF,FCOST,
+     1   XOBJ,IMTHD,IMPR,IERR)
+*----
+*  AUGMENTED LAGRANGIAN
+*----
+      ELSE IF(IMTHD.EQ.4) THEN
+         CALL PLLA(IPOPT,N0,M0,APLUS,PDG,BPLUS,INPLUS,XDROIT,
+     1   FCOST,GF,XOBJ,IMPR,EPSIM,NCST,GRAD,CONTR,INEGAL,IERR)
+*----
+*  EXTERNAL PENALTY METHOD
+*----
+      ELSE IF(IMTHD.EQ.5) THEN
+         CALL PLPNLT(IPOPT,N0,M0,APLUS,PDG,BPLUS,INPLUS,XDROIT,
+     1   FCOST,GF,XOBJ,IMPR,EPSIM,NCST,GRAD,CONTR,INEGAL,IERR)
+      ENDIF
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(APLUS,GF,BPLUS)
+      DEALLOCATE(INPLUS)
+      RETURN
+*
+1000  FORMAT(/' PLDRV: INCONSISTENCY BETWEEN M0 AND ME+MI ',3I5)
+      END