1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
|
*DECK PLQPRO
SUBROUTINE PLQPRO(N0,M1,MAXM,COUT,QUAD,APLUS,BPLUS,BINF,BSUP,XOBJ,
> F,EPS,IMPR,IERR)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Minimizes a quadratic programming problem using LEMKE pivot.
* PLQPRO=Quadratic PROgramming problem resolution.
* The objective function is (1/2)X^T QUAD X + COUT^T X.
*
*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
*
*Parameters: input
* N0 number of control variables.
* M1 number of constraints.
* MAXM first dimension of matrix APLUS.
* COUT costs of control variables.
* QUAD quadratic matrix for the objective function. QUAD is assumed
* to be symmetric positive definite.
* APLUS coefficient matrix for the linear constraints.
* BPLUS right hand sides corresponding to the coefficient matrix.
* BINF lower bounds of control variables.
* BSUP upper bounds of control variables.
* EPS tolerence used for pivoting.
* IMPR print flag.
*
*Parameters: ouput
* XOBJ control variables.
* F objective function.
* IERR return code (=0: normal completion).
*
*-----------------------------------------------------------------------
*
IMPLICIT NONE
*----
* SUBROUTINE ARGUMENTS
*----
INTEGER N0,M1,MAXM,IMPR,IERR
DOUBLE PRECISION BPLUS(M1),BINF(N0),BSUP(N0),XOBJ(N0),EPS,
> APLUS(MAXM,N0),COUT(N0),QUAD(N0,N0),F
*----
* LOCAL VARIABLES
*----
INTEGER N,NP1,NP2,I,J,II,IR
DOUBLE PRECISION XVAL
CHARACTER*4 ROW(7)
INTEGER, ALLOCATABLE, DIMENSION(:) :: IROW,ICOL
DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: P
*----
* SCRATCH STORAGE ALLOCATION
*----
ALLOCATE(IROW(2*N0+M1),ICOL(2*N0+M1+1))
ALLOCATE(P(2*N0+M1,2*N0+M1+2))
*----
* CHECK SYMMETRY OF MATRIX QUAD.
*----
DO I=1,N0
DO J=I+1,N0
IF(ABS(QUAD(I,J)-QUAD(J,I)).GT.EPS) THEN
IERR=10 ! matrix QUAD is not symmetric
WRITE(6,1000) IERR
GO TO 500
ENDIF
ENDDO
ENDDO
*----
* SET-UP AND SOLVE THE PARAMETRIC COMPLEMENTARITY PROBLEM.
*----
N=2*N0+M1
NP1=N+1
NP2=N+2
*
P(:N,:NP2)=0.0D0
DO I=1,N0
P(:N0,:N0)=QUAD(:N0,:N0)
XVAL=0.0D0
DO J=1,N0
XVAL=XVAL-QUAD(I,J)*BINF(J)
ENDDO
P(I,NP2)=COUT(I)-XVAL
P(N0+M1+I,NP2)=BSUP(I)-BINF(I)
P(I,N0+M1+I)=1.0D0
P(N0+M1+I,I)=-1.0D0
DO J=1,M1
P(I,N0+J)=APLUS(J,I)
ENDDO
ENDDO
*
DO I=1,M1
XVAL=0.0D0
DO J=1,N0
XVAL=XVAL+APLUS(I,J)*BINF(J)
ENDDO
P(N0+I,NP2)=BPLUS(I)-XVAL
DO J=1,N0
P(N0+I,J)=-APLUS(I,J)
ENDDO
ENDDO
*
DO I=1,N
IROW(I)= I
ICOL(I)=-I
P(I,NP1)=1.0D0
ENDDO
ICOL(NP1)=-NP1
*
CALL PLLEMK(N,NP2,EPS,IMPR,P,IROW,ICOL,IERR)
IF (IERR.GE.1) THEN
WRITE(6,1000) IERR
GO TO 500
ENDIF
*
IF(IMPR.GE.3) THEN
WRITE(6,2000)
DO I=1,N,7
II=MIN0(I+6,N)
DO J=I,II
IF (IROW(J).LT.0) THEN
WRITE (ROW(J-I+1),'(1HX,I3.3)') (-IROW(J))
ELSE
WRITE (ROW(J-I+1),'(1HY,I3.3)') IROW(J)
ENDIF
ENDDO
WRITE(6,3000) (ROW(J-I+1),P(J,NP2),J=I,II)
ENDDO
ENDIF
*
XOBJ(:N0)=BINF(:N0)
DO I=1,N
IR=-IROW(I)
IF ((IR.GT.0).AND.(IR.LE.N0)) THEN
XOBJ(IR)=BINF(IR)+P(I,NP2)
ENDIF
ENDDO
F=DOT_PRODUCT(COUT(:N0),XOBJ(:N0))
XVAL=DOT_PRODUCT(MATMUL(QUAD(:N0,:N0),XOBJ(:N0)),XOBJ(:N0))
F=F+0.5*XVAL
*----
* SCRATCH STORAGE DEALLOCATION
*----
500 DEALLOCATE(P)
DEALLOCATE(ICOL,IROW)
RETURN
*
1000 FORMAT(//,5X,'PLQPRO: FAILURE OF THE LINEAR COMPLEMENTARITY SO',
> 'LUTION (IERR=',I5,').')
2000 FORMAT(//,5X,'SOLUTION OF THE LINEAR COMPLEMENTARITY PROBLEM :',
> '*** X: KUHN-TUCKER MULTIPLIERS ;',5X,
> '*** Y: SLACK VARIABLES ',/)
3000 FORMAT(7(1X,A4,'=',E12.5),/)
END
|
