Initial commit from Polytechnique Montreal

1 files changed, 341 insertions, 0 deletions

diff --git a/Donjon/src/PLMAP1.f b/Donjon/src/PLMAP1.f
new file mode 100644
index 0000000..3f4d677
--- /dev/null
+++ b/Donjon/src/PLMAP1.f
@@ -0,0 +1,341 @@
+*DECK PLMAP1
+      SUBROUTINE PLMAP1(N0,M0,APLUS,PDG,BPLUS,INPLUS,XDROIT,COUT,OBJ,
+     > XOBJ,IMTHD,IMPR,IERR,BINF,BSUP,SCALE,PX,RX,DELTA,BGAR)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Solves a linear optimization problem with quadratic constraint using
+* the method of approximation programming (MAP). 
+* PLMAP1 = Linear Programmation MAP1
+*
+*Reference:
+*  R.E. Griffith and R.A. Stewart, 'A non-linear programming technique
+*  for the optimization of continuous processing systems', Management
+*  Science, Vol. 7, NO. 4, 379 (1961).
+*
+*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
+* 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.
+* COUT    costs of control variables.
+* OBJ     objective function.
+* XOBJ    control variables.
+* IMTHD   type of solution (=1: SIMPLEX/LEMKE; =3: MAP).
+* IMPR    print flag.
+*
+*Parameters: ouput
+* IERR    return code (=0: normal completion).
+*
+*Parameters: scratch
+* BINF
+* BSUP
+* SCALE
+* PX
+* RX
+* DELTA
+* BGAR
+*
+*-----------------------------------------------------------------------
+*
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER N0,M0,INPLUS(M0+1),IMTHD,IMPR,IERR
+      DOUBLE PRECISION PDG(N0),BPLUS(M0+2),XDROIT,XOBJ(N0),BINF(N0),
+     > BSUP(N0),SCALE(N0),PX(N0),RX(N0),DELTA(N0),BGAR(M0+1),
+     > APLUS(M0+2,N0),COUT(N0),OBJ
+*----
+*  LOCAL VARIABLES
+*----
+      DOUBLE PRECISION  DELF,X,DUMY,ZMAX,XVAL,SCAL,DELX,TEMP,ERR,
+     > CONT,EPSIR,EPS,EPSS
+      INTEGER     ITER,ITMAX,I,J,M
+      CHARACTER   CLNAME*6
+      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: ZMAXC,BGAR0
+*----
+*  DATA STATEMENTS
+*----
+      DATA     CLNAME /'PLMAP'/
+*
+      EPS   = 0.0001D0
+      EPSIR =EPS
+      ITMAX = 100
+*----
+*  CONTROL-VARIABLE SCALING
+*----
+      DO 10 J=1,N0
+         SCAL = SQRT(XDROIT/PDG(J))
+         SCALE(J) = SCAL
+         COUT(J) = COUT(J)*SCAL
+*
+         PDG(J)  = 1.0D0
+         BINF(J) = 0.0D0
+         BSUP(J) = 0.0D0
+*
+         DO 20 I=1,M0
+            APLUS(I,J) = APLUS(I,J)*SCAL
+   20    CONTINUE
+*
+  10  CONTINUE
+*----
+*  PRINT TABLES AFTER SCALING OF CONTROL VARIABLES
+*----
+      IF(IMPR.GE.5) THEN
+        CALL PLNTAB(COUT,APLUS,INPLUS,BPLUS,PDG,BINF,BSUP,
+     >              N0,M0,CLNAME//' AFTER SCALING OF CONTROL VARIABLES')
+      ENDIF
+*
+      XDROIT=1.0D0
+*----
+*  CONSTRAINT SCALING
+*----
+      ALLOCATE(ZMAXC(M0))
+      DO 30 I=1,M0
+         ZMAX = ABS(BPLUS(I))
+*
+         DO 40 J=1,N0
+            ZMAX = MAX(ZMAX,ABS(APLUS(I,J)))
+   40    CONTINUE
+         BGAR(I) = BPLUS(I)/ZMAX
+*
+         DO 42 J=1,N0
+            APLUS(I,J) = APLUS(I,J)/ZMAX
+   42    CONTINUE
+         ZMAXC(I) = ZMAX
+  30  CONTINUE
+*----
+*  COST SCALING
+*----
+      ZMAX = 0.0D0
+      DO 45 J=1,N0
+         ZMAX = MAX(ZMAX,ABS(COUT(J)))
+   45 CONTINUE
+      DO 50 J=1,N0
+         COUT(J) = COUT(J)/ZMAX
+   50 CONTINUE
+*----
+*  PRINT TABLES AFTER SCALING OF COSTS AND CONSTRAINTS
+*----
+      IF(IMPR.GE.5) THEN
+        CALL PLNTAB(COUT,APLUS,INPLUS,BGAR,PDG,BINF,BSUP,N0,M0,
+     >              CLNAME//' AFTER SCALING OF COSTS AND CONSTRAINTS')
+      ENDIF
+      ALLOCATE(BGAR0(M0))
+      DO 52 I=1,M0
+         BGAR0(I) = BGAR(I)
+   52 CONTINUE
+*----
+*  INITIAL ESTIMATES
+*----
+      DELX = SQRT(XDROIT)
+      EPSS = EPS*DELX
+*
+      DO 55 I=1,N0
+         DELTA(I) = DELX/10.0
+         RX(I)    = 0.0
+   55 CONTINUE
+      TEMP = DELX/SQRT(REAL(N0))/10
+*----
+*  MAP ITERATIONS
+*----
+      ITER = 0
+   60 ITER = ITER + 1
+      CONT = 0.0
+      DO 70 I=1,M0+1
+         BGAR(I) = BGAR0(I)
+   70 CONTINUE
+*----
+*  CONTROL VARIABLE BOUNDS
+*----
+      DO 90 I=1,N0
+         IF(ITER.EQ.1) THEN
+!            XOBJ(I) = EPSIR*10.0
+            XOBJ(I) = 0.0
+            BINF(I) = -TEMP
+            BSUP(I) = TEMP
+         ELSE
+            BINF(I) = -DELTA(I)
+            BSUP(I) = DELTA(I)
+         ENDIF
+*----
+*  LINEARIZATION OF THE QUADRATIC CONSTRAINT
+*----
+         CONT = CONT + XOBJ(I)**2
+         APLUS(M0+1,I) = 2.0*XOBJ(I)
+         DO 95 J=1,M0
+            BGAR(J) = BGAR(J) - APLUS(J,I)*XOBJ(I)
+   95    CONTINUE
+  90  CONTINUE
+*
+      INPLUS(M0+1) = 1
+      BGAR(M0+1)   = XDROIT - CONT
+      M = M0 + 1
+*----
+*  REORGANIZE TABLES FOR SIMPLEX
+*----
+      DO 120 I=1,M
+         DUMY = 0.0D0
+*
+         DO 100 J=1,N0
+            DUMY = DUMY + APLUS(I,J)*BINF(J)
+  100    CONTINUE
+*
+         BGAR(I) = BGAR(I) - DUMY
+         IF(BGAR(I).GE.0.0) GOTO 120
+*
+         DO 110 J=1,N0
+            APLUS(I,J) = -APLUS(I,J)
+  110    CONTINUE
+*
+         BGAR(I)   = -BGAR(I)
+         BGAR0(I)= -BGAR0(I)
+         BPLUS(I)  = -BPLUS(I)
+         INPLUS(I) = -INPLUS(I)
+*
+  120 CONTINUE
+*
+      DO 130 J=1,N0
+         BSUP(J) = BSUP(J) - BINF(J)
+         BINF(J) = 0.0
+  130 CONTINUE
+*----
+*  PRINT SIMPLEX TABLES
+*----
+      IF(IMPR.GE.5) THEN
+        CALL PLNTAB(COUT,APLUS,INPLUS,BGAR ,XOBJ  ,BINF  ,BSUP,N0,M0,
+     >              CLNAME//' AFTER REORGANIZATION FOR SIMPLEX')
+      ENDIF
+*----
+*  SOLUTION OF A LINEAR PROGRAMMING PROBLEM USING THE SIMPLEX
+*----
+      CALL PLSPLX(N0,M,M0+2,1,COUT,APLUS,BGAR,INPLUS,BINF,BSUP,PX,
+     >          DELF,EPSS,IMTHD,IMPR,IERR)
+*
+      DO 140 I=1,N0
+         IF(ITER.EQ.1) THEN
+            PX(I) = PX(I) - TEMP
+         ELSE
+            PX(I) = PX(I) - DELTA(I)
+         ENDIF
+  140 CONTINUE
+*----
+*  SOLUTION OF CURRENT ITERATION
+*----
+      IF(IMPR.GE.2) THEN
+         IF(((ITER.GE.1).AND.(IMPR.LE.2)).OR.(IMPR.GE.3)) THEN
+            WRITE (6,1000)
+         ENDIF
+         WRITE (6,2000) ITER,DELF,(PX(I),I=1,N0)
+      ENDIF
+*----
+*  DEGENERESCENCE OR EPS TOO SMALL
+*----
+      IF(IERR.EQ.1) THEN
+         WRITE(6,3000) ITER
+         IERR = 3
+         RETURN
+*----
+*  NO SOLUTION IF ITER=1
+*----
+      ELSE IF(IERR.EQ.2) THEN
+         IF(IMPR.GE.1) WRITE(6,4000) ITER
+         IF(ITER.GE.ITMAX) RETURN
+      ENDIF
+*
+      ERR = 0.0
+      DO 160 I=1,N0
+*
+         IF((RX(I)*PX(I).LT.0.0).AND.(IERR.EQ.0)) THEN
+            DELTA(I) = DELTA(I)*0.5
+         ENDIF
+*
+         RX(I)   = PX(I)
+         XOBJ(I) = XOBJ(I) + PX(I)
+         ERR = ERR + PX(I)**2
+  160 CONTINUE
+*
+      ERR  = SQRT(ERR)
+      EPSS = EPS*DELX/10.0
+*
+      IF(IMPR.GE.1) THEN
+         WRITE(6,2000) ITER,DELF,(XOBJ(I),I=1,N0)
+         WRITE(6,2000) ITER,0.0,(DELTA(I),I=1,N0)
+      ENDIF
+*
+      IF(ERR.LE.EPSS) THEN
+         IERR = 0
+         GOTO 170
+      ENDIF
+*
+      IF(ITER.GE.ITMAX) THEN
+         IERR = 5
+         WRITE (6,5000) ITER
+         RETURN
+      ENDIF
+      GO TO 60
+*----
+*  RESCALE BACK AND PRINT THE SOLUTION
+*----
+  170 DO 175 J=1,N0
+         SCAL    = SCALE(J)
+         COUT(J) = COUT(J)*ZMAX/SCAL
+         XOBJ(J) = XOBJ(J)*SCAL
+         PDG(J)  = XDROIT/SCAL**2
+*
+         DO 177 I=1,M0
+            APLUS(I,J) = APLUS(I,J)/SCAL
+  177    CONTINUE
+  175 CONTINUE
+*
+      X   = 0.0D0
+      OBJ = 0.0D0
+      DO 180 J=1,N0
+         X   = X   + PDG(J)*XOBJ(J)*XOBJ(J)
+         OBJ = OBJ + XOBJ(J)*COUT(J)
+  180 CONTINUE
+*
+      IF(IMPR.GE.1) THEN
+         WRITE (6,6000) OBJ,X,(XOBJ(J),J=1,N0)
+         IF(M0.GT.0) WRITE (6,7000)
+*
+         DO 190 I=1,M0
+            XVAL = BPLUS(I)
+            DO 185 J=1,N0
+               XVAL = XVAL - APLUS(I,J)*XOBJ(J)*ZMAXC(I)/SCALE(J)
+  185       CONTINUE
+            WRITE (6,8000) I,XVAL
+  190    CONTINUE
+      ENDIF
+      DEALLOCATE(ZMAXC,BGAR0)
+      RETURN
+*
+1000  FORMAT(/,5X,'ITERATION',8X,'DELF',5X,'CONTROL VARIABLES')
+2000  FORMAT(5X,I6,5X,8E12.4,/,(28X,5E12.4))
+3000  FORMAT(5X,I6,5X,'DEGENERESCENCE OR EPS TOO SMALL')
+4000  FORMAT(5X,I6,5X,'NO SOLUTION')
+5000  FORMAT(5X,I6,5X,'MAXIMUM ITERATION REACHED')
+6000  FORMAT(//,5X,'FINAL SOLUTION (MAP1-SIMPLEX) ',
+     >        /,5X,'------------------------',
+     >        /,5X,'OBJECTIVE FUNCTION    : ',1P,E12.5,
+     >        /,5X,'QUADRATIC CONSTRAINT  : ',1P,E12.5,
+     >        /,5X,'CONTROL VARIABLES     : ',/,(10X,10E12.4))
+7000  FORMAT(//,5X,'CONSTRAINT DEVIATIONS : ',/)
+8000  FORMAT(2X,I3,'...',2X,1P,D12.4)
+      END