1 /****************************************************************
2 Copyright (C) 2011 AMPL Optimization LLC; written by David M. Gay.
3
4 Permission to use, copy, modify, and distribute this software and its
5 documentation for any purpose and without fee is hereby granted,
6 provided that the above copyright notice appear in all copies and that
7 both that the copyright notice and this permission notice and warranty
8 disclaimer appear in supporting documentation.
9
10 The author and AMPL Optimization LLC disclaim all warranties with
11 regard to this software, including all implied warranties of
12 merchantability and fitness. In no event shall the author be liable
13 for any special, indirect or consequential damages or any damages
14 whatsoever resulting from loss of use, data or profits, whether in an
15 action of contract, negligence or other tortious action, arising out
16 of or in connection with the use or performance of this software.
17 ****************************************************************/
18
19 #include "nlp.h"
20 #define SKIP_NL2_DEFINES
21 #undef f_OPNUM
22 #include "psinfo.h"
23 #include "jacpdim.h"
24 #undef cde
25 #undef ps_func
26
27 struct
28 MPEC_Adjust {
29 int *cc; /* indices of constraints originally involved in complementarities */
30 int *cce; /* end of cc array */
31 int *ck; /* kind of complementarity modification */
32 real *rhs1; /* start of added constraint right-hand sides */
33 cgrad **Cgrda; /* linear terms added to existing constraints */
34 int incc; /* for incrementing constraint lhs and rhs */
35 int incv; /* for incrementing variable lhs and rhs */
36 int m0; /* original number of constraints */
37 int n0; /* original number of variables */
38 };
39
40 static expr_n ZeroExpr = { f_OPNUM_ASL, 0. };
41
42 static void
reverse(int * a,int * b)43 reverse(int *a, int *b)
44 {
45 int t;
46
47 while(--b > a) {
48 t = *b;
49 *b = *a;
50 *a++ = t;
51 }
52 }
53
54 void
mpec_adjust_ASL(ASL * asl)55 mpec_adjust_ASL(ASL *asl)
56 {
57 MPEC_Adjust *mpa;
58 cde *cd;
59 cde2 *cd2;
60 cgrad **Cgrd, **Cgrd1, **Cgrda, *cg, *cg1, *ncg, **pcg;
61 char *hx0;
62 int *cc, *ck, *cv, *ind1, *ind2, *vm;
63 int i, incc, incv, j, k, m, m0, n, n0, n1, n2, nb, ncc, ncc0, nib, nib0;
64 int nnv, v1, v2, v3, v4;
65 real *Lc, *Lc0, *Lc1, *Lv, *Lv0, *Lv1, *Uc, *Uc0, *Uc1, *Uv, *Uv0, *Uv1;
66 real a, b, *x;
67 size_t nz, nz0, nznew;
68 extern void f_OPVARVAL_ASL(), f2_VARVAL_ASL();
69
70 n = n0 = n_var;
71 n2 = asl->i.n_var0;
72 nib = niv + nbv;
73 n1 = n - nib;
74 nib0 = n - nib; /* offset of first linear integer or binary variable */
75 m = m0 = n_con;
76 nz = nz0 = nZc;
77 cv = cvar;
78 Cgrd = Cgrad;
79 Cgrd1 = Cgrd + m;
80 incc = incv = 1;
81 Lc0 = LUrhs;
82 if (!(Uc0 = Urhsx)) {
83 Uc0 = Lc0 + 1;
84 incc = 2;
85 }
86 Lv0 = LUv;
87 if (!(Uv0 = Uvx)) {
88 Uv0 = Lv0 + 1;
89 incv = 2;
90 }
91 ncc = ncc0 = n_cc;
92 Lc1 = Lc0 + m*incc;
93 Uc1 = Uc0 + m*incc;
94 Lv1 = Lv0 + n*incv;
95 Uv1 = Uv0 + n*incv;
96
97 for(i = k = 0; i < m0; ++i)
98 if ((j = cv[i])) {
99 ++k;
100 Lc = Lc0 + incc*i;
101 Uc = Uc0 + incc*i;
102 nb = (*Lc > negInfinity) + (*Uc < Infinity);
103 /* nb == 0 or 1 */
104 if (!nb) {
105 m += 2;
106 n += 4;
107 nz += 6;
108 ++ncc;
109 }
110 else {
111 Lv = Lv0 + incv*(j-1);
112 if (*Lv != 0.) {
113 ++m;
114 ++n;
115 nz += 2;
116 }
117 /* Even if constraint i has the form v >= 0, */
118 /* add a new variable v1 >= 0 and change the */
119 /* constraint to v1 = v - rhs, in case v is */
120 /* involved in more than one complementarity */
121 ++n;
122 ++nz;
123 }
124 }
125 if (k != ncc0) {
126 fprintf(Stderr,
127 "\nERROR: mpec_adjust saw %d rather than %d incoming complementarities.\n",
128 k, ncc0);
129 exit(1);
130 }
131 n_var = n;
132 n_con = m;
133 asl->i.n_var1 += nnv = n - n0;
134 if (n_obj)
135 adjust_zerograds_ASL(asl, nnv);
136 if (n_conjac[1] >= m0)
137 n_conjac[1] = m;
138 nzc = nZc = nz;
139 n_cc = ncc;
140 nznew = nz - nz0;
141 ncg = (cgrad*)M1alloc(2*(ncc + ncc0)*sizeof(int) + nznew*sizeof(cgrad)
142 + ncc0*sizeof(cgrad*) + sizeof(MPEC_Adjust));
143 asl->i.mpa = mpa = (MPEC_Adjust*)(ncg + nznew);
144 Cgrda = mpa->Cgrda = (cgrad**)(mpa + 1);
145 asl->i.ccind1 = ind1 = (int*)(Cgrda + ncc0);
146 asl->i.ccind2 = ind2 = ind1 + ncc;
147 mpa->cc = cc = ind2 + ncc;
148 mpa->ck = ck = mpa->cce = cc + ncc0;
149 mpa->m0 = m0;
150 mpa->n0 = n0 - nib;
151 mpa->rhs1 = Lc1;
152 mpa->incc = incc;
153 mpa->incv = incv;
154 if (nib) {
155 vm = get_vcmap_ASL(asl, ASL_Sufkind_var);
156 /* Three reverse calls move nib values of vm up nnv places. */
157 j = n0 - nib;
158 reverse(vm+j, vm + n0 + nnv);
159 reverse(vm+j, vm + j + nnv);
160 reverse(vm + j + nnv, vm + n0 + nnv);
161 i = n0 + nnv;
162 while(--i >= n0) {
163 j = i - nnv;
164 Lv0[incv*i] = Lv0[incv*j];
165 Uv0[incv*i] = Uv0[incv*j];
166 }
167 if ((x = X0)) {
168 i = n0 + nnv;
169 while(--i >= n0)
170 x[i] = x[i-nnv];
171 for(i = n0 - nnv; i < n0; ++i)
172 x[i] = 0.;
173 if ((hx0 = havex0)) {
174 for(i = n0 + nnv; --i >= n0; )
175 hx0[i] = hx0[i-nnv];
176 for(i = n0 - nnv;i < n0; ++i)
177 hx0[i] = 0;
178 }
179 }
180 Lv1 -= j = incv*nib;
181 Uv1 -= j;
182 }
183 else if ((x = X0)) {
184 memset(x + n0, 0, nnv*sizeof(real));
185 if ((hx0 = havex0))
186 memset(hx0 + n0, 0, nnv);
187 }
188
189 #define vset(x,y) *x = y; x += incv;
190 for(i = 0; i < m0; ++i)
191 if ((j = cv[i])) {
192 if (j > nib0)
193 j += nnv;
194 *cc++ = i;
195 pcg = &Cgrd[i];
196 cg = 0;
197 while((cg1 = *pcg))
198 pcg = &(cg = cg1)->next;
199 *Cgrda++ = cg;
200 Lc = Lc0 + incv*i;
201 Uc = Uc0 + incc*i;
202 Lv = Lv0 + incv*--j;
203 Uv = Uv0 + incv*j;
204 a = *Lc;
205 b = *Uc;
206 *ck++ = nb = (a > negInfinity) + (b < Infinity);
207 if (nb == 0) {
208 /* change L <= v = _svar[j] <= U */
209 /* and -Infinity <= body <= Infinity into */
210 /* v1 = v - L >= 0, v2 = U - v >= 0, */
211 /* v3 - v4 = body, v3 >= 0, v4 >= 0, */
212 /* v1 complements v3, v2 complements v4 */
213
214 *Lc = *Uc = 0.;
215 *ind1++ = v1 = n1++;
216 *ind1++ = v2 = n1++;
217 *ind2++ = v3 = n1++;
218 *ind2++ = v4 = n1++;
219 for(k = 0; k < 4; ++k) {
220 vset(Lv1, 0.);
221 vset(Uv1, Infinity);
222 }
223 ncg[1].varno = n2+3;
224 ncg[1].coef = 1.;
225 ncg[1].next = 0;
226 ncg[0].varno = n2+2;
227 ncg[0].coef = -1.;
228 ncg[0].next = &ncg[1];
229 *pcg = ncg;
230 ncg += 2;
231 ncg[1].varno = n2;
232 ncg[1].coef = -1.;
233 ncg[1].next = 0;
234 ncg[0].varno = j;
235 ncg[0].coef = 1.;
236 ncg[0].next = &ncg[1];
237 *Lc1 = *Uc1 = *Lv;
238 Lc1 += incc;
239 Uc1 += incc;
240 *Cgrd1++ = ncg;
241 ncg += 2;
242 ncg[1].varno = n2+1;
243 ncg[1].coef = 1.;
244 ncg[1].next = 0;
245 ncg[0].varno = j;
246 ncg[0].coef = 1.;
247 ncg[0].next = &ncg[1];
248 *Lc1 = *Uc1 = *Uv;
249 Lc1 += incc;
250 Uc1 += incc;
251 *Cgrd1++ = ncg;
252 ncg += 2;
253 n2 += 4;
254 }
255 else {
256 /*nb == 1*/
257 v1 = j;
258 if (*Lv != 0.) {
259 /* For v = _svar[j], replace */
260 /* v >= a with v1 = v - a, v1 >= 0, or */
261 /* v <= b with v1 = b - v, v1 >= 0 */
262 v1 = n1++;
263 vset(Lv1, 0.);
264 vset(Uv1, Infinity);
265 ncg[1].varno = n2++;
266 ncg[1].next = 0;
267 ncg[0].varno = j;
268 ncg[0].coef = 1.;
269 ncg[0].next = &ncg[1];
270 if (*Lv > negInfinity) {
271 ncg[1].coef = -1.;
272 *Lc1 = *Uc1 = *Lv;
273 }
274 else {
275 ncg[1].coef = 1.;
276 *Lc1 = *Uc1 = *Uv;
277 }
278 Lc1 += incc;
279 Uc1 += incc;
280 *Cgrd1++ = ncg;
281 ncg += 2;
282 }
283 *ind1++ = v1;
284 *ind2++ = n1++;
285 ncg->varno = n2++;
286 ncg->next = 0;
287 vset(Lv1, 0.);
288 vset(Uv1, Infinity);
289 if (*Lv > negInfinity) {
290 ncg->coef = -1.;
291 *Uc = *Lc;
292 }
293 else {
294 ncg->coef = 1.;
295 *Lc = *Uc;
296 }
297 *pcg = ncg++;
298 }
299 }
300 #undef vset
301 i = m0;
302 asl->i.n_con1 += k = m - m0;
303 switch(asl->i.ASLtype) {
304 case ASL_read_pfg:
305 memset(((ASL_pfg*)asl)->P.cps + m0, 0, k*sizeof(ps_func));
306 cd = ((ASL_pfg*)asl)->I.con_de_;
307 goto have_cd;
308 case ASL_read_f:
309 case ASL_read_fg:
310 cd = ((ASL_fg*)asl)->I.con_de_;
311 have_cd:
312 while(i < m)
313 cd[i++].e = (expr*)&ZeroExpr;
314 break;
315 case ASL_read_fgh:
316 cd2 = ((ASL_fgh*)asl)->I.con2_de_;
317 goto have_cd2;
318 case ASL_read_pfgh:
319 memset(((ASL_pfgh*)asl)->P.cps + m0, 0, k*sizeof(ps_func2));
320 cd2 = ((ASL_pfgh*)asl)->I.con2_de_;
321 have_cd2:
322 while(i < m)
323 cd2[i++].e = (expr2*)&ZeroExpr;
324 }
325
326 }
327
328 void
mpec_auxvars_ASL(ASL * asl,real * c,real * x)329 mpec_auxvars_ASL(ASL *asl, real *c, real *x)
330 {
331 /* Adjust variables added by mpec_adjust_ASL() so the constraints */
332 /* added by mpec_adjust_ASL() are satisfied. */
333
334 MPEC_Adjust *mpa;
335 cgrad **Cg, **Cga, *cg;
336 int *cc, *cce, *ck, *cv, i, incc, incv, j, m0, n0, *vmi;
337 real *Lc, *Lc0, *Lc1, *Lv0, *ca, t;
338
339 mpa = asl->i.mpa;
340 cv = cvar;
341 cc = mpa->cc;
342 cce = mpa->cce;
343 ck = mpa->ck;
344 Cga = mpa->Cgrda;
345 m0 = mpa->m0;
346 n0 = mpa->n0;
347 Cg = Cgrad + m0;
348 ca = c + m0;
349 Lc0 = LUrhs;
350 Lc1 = mpa->rhs1;
351 Lv0 = LUv;
352 incc = mpa->incc;
353 incv = mpa->incv;
354 vmi = get_vminv_ASL(asl);
355 do {
356 t = c[i = *cc++];
357 c[i] = 0.;
358 j = cv[i] - 1;
359 for(cg = *Cga++; cg->varno < n0; cg = cg->next);
360 Lc = Lc0 + i*incc;
361 if (!*ck++) {
362 if (t >= 0.)
363 x[vmi[cg->varno]] = t;
364 else {
365 cg = cg->next;
366 x[vmi[cg->varno]] = -t;
367 }
368 cg = (*Cg++)->next;
369 x[vmi[cg->varno]] = x[j] - *Lc1;
370 *ca++ = *Lc1;
371 Lc1 += incc;
372 cg = (*Cg++)->next;
373 x[vmi[cg->varno]] = *Lc1 - x[j];
374 *ca++ = *Lc1;
375 Lc1 += incc;
376 }
377 else {
378 x[vmi[cg->varno]] = cg->coef*(*Lc - t);
379 c[i] = *Lc;
380 if (Lv0[incv*j] != 0.) {
381 cg = (*Cg++)->next;
382 x[vmi[cg->varno]] = cg->coef*(*Lc1 - x[j]);
383 *ca++ = *Lc1;
384 Lc1 += incc;
385 }
386 }
387 } while(cc < cce);
388 }
389