1 /****************************************************************
2 Copyright (C) 2002 AMPL Optimization LLC
3 All Rights Reserved
4 Based largely on suf_sos.c, which bears the following Copyright
5 notice and disclaimer. AMPL Optimization LLC similarly disclaims
6 all warranties with regard to this software and grants permission
7 to use, copy, modify, and distribute it and its documentation. */
8
9 /****************************************************************
10 Copyright (C) 1999-2001 Lucent Technologies
11 All Rights Reserved
12
13 Permission to use, copy, modify, and distribute this software and
14 its documentation for any purpose and without fee is hereby
15 granted, provided that the above copyright notice appear in all
16 copies and that both that the copyright notice and this
17 permission notice and warranty disclaimer appear in supporting
18 documentation, and that the name of Lucent or any of its entities
19 not be used in advertising or publicity pertaining to
20 distribution of the software without specific, written prior
21 permission.
22
23 LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
24 INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
25 IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
26 SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
27 WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
28 IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
29 ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
30 THIS SOFTWARE.
31 ****************************************************************/
32
33 /* This is a replacement for suf_sos() for solvers that know about
34 SOS1 sets but do not supply the requisite convexity constraints
35 (or do not know about SOS2 sets), or that know nothing of SOS
36 sets but do handle integer variables (in which case sos_finish()
37 should be called with null values for the arguments in which
38 SOS information is returned, such as sosref_p).
39
40 Before calling the .nl reader, one invokes
41
42 void *SI = sos_add(nl,flags);
43
44 (where nl is the open FILE* that will be passed to the .nl reader)
45 to prepare for sos_finish to insert convexity constraints and adjust
46 n_con, n_var, nzc, and nbv, etc. After calling the .nl reader,
47 one invokes
48
49 if (SI)
50 sos_finish(SI,...);
51
52 to finish modifying the problem -- adding the convexity constraints,
53 etc., for SOS sets implied by suffixes .sosno and .ref, and retaining
54 SOS1 and SOS2 sets added when AMPL linearizes nonconvex piecewise-
55 linear terms.
56 */
57
58 #include "asl.h"
59 #define SKIP_NL2_DEFINES
60 #include "nlp.h"
61 #include "nlp2.h"
62 #include "asl_pfg.h"
63 #include "asl_pfgh.h"
64 #undef cde
65
66 #ifdef __cplusplus
67 extern "C" {
68 static int compar(const void*, const void*, void*);
69 static int rcompar(const void*, const void*, void*);
70 }
71 #endif
72
73 #undef nzc
74
75 typedef struct
76 SOSinfo {
77 SufDesc mysd[2];
78 Char **zap[5];
79 ASL *asl;
80 SufDesc *grefd, *refd;
81 int *col1, *g0, **gp0, **gp1, **gpe, *v0, *ve, *z;
82 int nnc, nnv, nnz, nsos, nsos1, nsos2, nsosnz, nsosnz1, nzc;
83 }
84 SOSinfo;
85
86 typedef struct
87 Coninfo {
88 cgrad **cgp, *cg;
89 int *z;
90 real *clu, *cu, *lu, *u;
91 int m, nbin;
92 }
93 Coninfo;
94
95 static real pl_bigM, mpl_bigM;
96
97 static SufDesc*
refd_copy(SufDesc * nu,SufDesc * ol,int n)98 refd_copy(SufDesc *nu, SufDesc *ol, int n)
99 {
100 if (!ol)
101 return ol;
102 memcpy(nu, ol, sizeof(SufDesc));
103 memcpy(nu->u.r = (real*)Malloc(n*sizeof(real)), ol->u.r, n*sizeof(real));
104 return nu;
105 }
106
107 static int
compar(const void * a,const void * b,void * v)108 compar(const void *a, const void *b, void *v)
109 {
110 int k;
111 Not_Used(v);
112 if ((k = *(int*)a - *(int*)b))
113 return k;
114 return ((int*)a)[1] - ((int*)b)[1];
115 }
116
117 static int
rcompar(const void * a,const void * b,void * v)118 rcompar(const void *a, const void *b, void *v)
119 {
120 real *r = (real *)v, t;
121 t = r[*(int*)a] - r[*(int*)b];
122 if (t == 0.)
123 return 0;
124 return t < 0. ? -1 : 1;
125 }
126
127 static int
refcomp(const void * a,const void * b,void * v)128 refcomp(const void *a, const void *b, void *v)
129 {
130 real *r = (real *)v, t;
131 t = r[((int*)a)[1]] - r[((int*)b)[1]];
132 if (t == 0.)
133 return 0;
134 return t < 0. ? -1 : 1;
135 }
136
137 static void
reorder(int * ind,real * ref,int j0,int k,int * p)138 reorder(int *ind, real *ref, int j0, int k, int *p)
139 {
140 int i, j, n, ti;
141 real tr;
142
143 ref += j0;
144 ind += j0;
145 n = k - j0;
146 for(i = 0; i < n; i++)
147 p[i] = i;
148 qsortv(p, n, sizeof(int), rcompar, ref);
149 for(i = 0; i < n; i++) {
150 if ((j = p[i]) > i) {
151 ti = ind[i];
152 tr = ref[i];
153 for(k = i;;) {
154 ind[k] = ind[j];
155 ref[k] = ref[j];
156 j = p[k = j];
157 p[k] = -1 - j;
158 if (j == i) {
159 ind[k] = ti;
160 ref[k] = tr;
161 break;
162 }
163 }
164 }
165 }
166 }
167
168 static void
turnon(int * z,int n)169 turnon(int *z, int n)
170 {
171 while(n-- > 0)
172 *z++ = 1;
173 }
174
175 static int
nonbinary(int k,Coninfo * CI,real * LU)176 nonbinary(int k, Coninfo *CI, real *LU)
177 {
178 int *z;
179 real *lu, *u;
180
181 if ((u = CI->u)) {
182 lu = CI->lu + k;
183 u += k;
184 }
185 else {
186 lu = CI->lu + 2*k;
187 u = lu + 1;
188 }
189 if ((LU[1] = *u) > pl_bigM)
190 LU[1] = pl_bigM;
191 if ((LU[0] = *lu) < mpl_bigM) {
192 LU[0] = mpl_bigM;
193 return 1;
194 }
195 if (LU[0] != 0.)
196 return 1;
197 if (!(z = CI->z) || !z[k] || *u != 1.)
198 return 1;
199 CI->nbin++;
200 return 0;
201 }
202
203 /* Call sos_add() prior to the .nl reader: sos_add scans incoming
204 suffixes for sosno, ref, sos and sosref and if found, increments
205 n_var and n_con in preparation for sos_finish(), which must be
206 called after the .nl reader to add convexity constraints, etc.
207 if sos_add has provided a nonzero return.
208
209 We use a private suftab. Users can omit the items in our suftab
210 from theirs, but should provide "priority" if desired.
211 */
212
213 static SufDecl
214 suftab[] = {
215 { (char*)"ref", 0, ASL_Sufkind_var | ASL_Sufkind_real },
216 { (char*)"sos", 0, ASL_Sufkind_var },
217 { (char*)"sosno", 0, ASL_Sufkind_var | ASL_Sufkind_real },
218 { (char*)"sosref", 0, ASL_Sufkind_var | ASL_Sufkind_real }
219 };
220
221 Char*
sos_add_ASL(ASL * asl,FILE * f,int flags)222 sos_add_ASL(ASL *asl, FILE *f, int flags)
223 {
224 Char **M1state1, **M1state2;
225 EdRead ER, *R;
226 SOSinfo *SI;
227 SufDesc *gd, *grefd, *refd, *sds[4], *vd;
228 int i, j, k, n, ng, nnc, nnv, nnz;
229 int nsos, nsos1, nsos2, nsosnz, nsosnz1, tg;
230 int *col1, *g0, *g, *g1, *ge, **gp, **gp0, **gp1, **gpe;
231 int nss[5], *v, *v0, *ve, *z, *zg;
232 long ft0;
233 real *gn;
234 size_t L;
235
236 SI = 0;
237 M1state1 = asl->i.Mbnext;
238 M1state2 = asl->i.Mblast;
239
240 /* Save and restore asl->i.suff stuff in case caller */
241 /* has already called suf_declare(). */
242
243 nss[0] = asl->i.nsuffixes;
244 memcpy(nss+1, asl->i.nsuff, 4*sizeof(int));
245 memcpy(sds, asl->i.suffixes, 4*sizeof(SufDesc*));
246 suf_declare(suftab, sizeof(suftab)/sizeof(SufDecl));
247
248 R = EdReadInit_ASL(&ER, asl, f, 0);
249 if (binary_nl)
250 xscanf = bscanf;
251 else
252 xscanf = ascanf;
253 ft0 = ftell(f);
254 while((i = edag_peek(R)) == 'S')
255 Suf_read_ASL(R, 0);
256 fseek(f, ft0, SEEK_SET);
257
258 nsos1 = nsos2 = 0;
259 gd = refd = grefd = 0;
260 z = 0;
261 if (!(flags & ASL_suf_sos_ignore_sosno)
262 && (gd = suf_get("sosno", ASL_Sufkind_var | ASL_Sufkind_input))
263 && (grefd = suf_get("ref", ASL_Sufkind_var | ASL_Sufkind_input)))
264 nsos1 = 1;
265 if (!(flags & ASL_suf_sos_ignore_amplsos)
266 && (vd = suf_get("sos", ASL_Sufkind_var | ASL_Sufkind_input))
267 && (refd = suf_get("sosref", ASL_Sufkind_var|ASL_Sufkind_input)))
268 nsos2 = 1;
269 else if (!nsos1)
270 goto ret;
271 nnc = nnv = nnz = nsos = nsosnz = nsosnz1 = 0;
272 n = n_var;
273 col1 = v0 = ve = 0;
274 if (nsos2) {
275 v = v0 = vd->u.i;
276 ve = v + n;
277 for(j = 0; j < n && v[j] <= 0; j++);
278 col1 = v + j;
279 for(; j < n; j++) {
280 if ((i = v[j]) & 2 && !(i&1)) {
281 nsosnz++;
282 if (nsos < (i >>= 4))
283 nsos = i;
284 }
285 }
286 }
287 gp0 = gp1 = gpe = 0;
288 g0 = 0;
289 if (nsos1) {
290 nsos1 = 0;
291 gn = gd->u.r;
292 for(i = ng = 0; i < n; i++) {
293 if ((j = gn[i]))
294 ng++;
295 }
296 if (!ng) {
297 /* should not happen */
298 if (nsos2)
299 goto havenz;
300 goto ret;
301 }
302 j = (niv | nbv | nlvbi | nlvci | nlvoi) ? n : 0;
303 L = (ng+1)*sizeof(int*) + (2*ng + j)*sizeof(int);
304 gp = gp0 = (int**)Malloc(L);
305 zg = g0 = (int*)(gpe = gp + ng + 1);
306 for(i = 0; i < n; i++)
307 if ((k = gn[i])) {
308 *zg++ = k;
309 *zg++ = i;
310 }
311 ge = zg;
312 qsortv(g0, ng, 2*sizeof(int), compar, 0);
313 z = 0;
314 if (j) {
315 memset(z = zg, 0, n*sizeof(int));
316 if (nlvbi)
317 turnon(z + (nlvb - nlvbi), nlvbi);
318 if (nlvci)
319 turnon(z + (nlvc - nlvci), nlvci);
320 if (nlvoi)
321 turnon(z + (nlvo - nlvoi), nlvoi);
322 if ((j = niv + nbv))
323 turnon(z + (n-j), j);
324 }
325 *gp++ = g = g1 = zg = g0;
326 j = -((tg = *g) < 0);
327 for(;;) {
328 if ((g += 2) >= ge || *g != tg) {
329 /* Ignore SOS1 sets of 1 element */
330 /* and SOS2 sets of <= 2 elements. */
331 if ((j += ((g-zg)>>1)) >= 2) {
332 nsosnz1 += j;
333 nsos1++;
334 if (!gp1 && tg > 0)
335 gp1 = gp - 1;
336 if (g1 == zg)
337 g1 = g;
338 else while(zg < g)
339 *g1++ = *zg++;
340 *gp++ = g1;
341 }
342 if (g >= ge)
343 break;
344 j = -((tg = *(zg = g)) < 0);
345 }
346 }
347 nsosnz += nsosnz1;
348 ge = g1;
349 gpe = gp - 1;
350 nsos += nsos1;
351 }
352 havenz:
353 if (!nsos) {
354 if (gp0)
355 free(gp0);
356 goto ret;
357 }
358 if (nsos1) {
359 if (!gp1)
360 gp1 = gpe;
361 /* scan SOS2s */
362 /* Here we must assume the worst case for nonbinary
363 * variables (since their bounds are not yet available):
364 * both lower and upper bounds nonzero. In sos_finish(),
365 * we may discover that we need fewer new constraints
366 * if some bounds turn out to be zero.
367 */
368 for(gp = gp0; gp < gp1; gp++) {
369 g = gp[0];
370 ge = gp[1] - 2;
371 j = (ge - g) >> 1;
372 nnv += j - 2;
373 nnc += j;
374 nnz += 4*j - 3;
375 for(i = 0; i < 2; i++) {
376 k = i ? ge[1] : g[1];
377 if (!z || !z[k]) {
378 nnv++;
379 nnc += 2;
380 nnz += 4;
381 }
382 }
383 nnc += k = (ge - g) - 2;
384 nnz += 3*k;
385 }
386 /* scan SOS1s */
387 for(; gp < gpe; gp++) {
388 g = gp[0];
389 ge = gp[1];
390 j = (ge - g) >> 1;
391 nnc++;
392 nnz += j;
393 nnc += k = ge - g;
394 nnv += k >> 1;
395 nnz += k << 1;
396 }
397 }
398 SI = (SOSinfo*)Malloc(sizeof(SOSinfo));
399 memset(SI->zap, 0, sizeof(SI->zap));
400 SI->asl = asl;
401 if (v0) {
402 L = n*sizeof(int);
403 memcpy(SI->v0 = (int*)Malloc(L), v0, L);
404 col1 = SI->v0 + (col1 - v0);
405 ve = (v0 = SI->v0) + n;
406 }
407 SI->col1 = col1;
408 SI->g0 = g0;
409 SI->gp0 = gp0;
410 SI->gp1 = gp1;
411 SI->gpe = gpe;
412 SI->grefd = grefd = refd_copy(&SI->mysd[0], grefd, n);
413 asl->i.nsufext[ASL_Sufkind_con] += SI->nnc = nnc;
414 asl->i.nsufext[ASL_Sufkind_var] += SI->nnv = nnv;
415 asl->i.nsufext[ASL_Sufkind_prob] += SI->nnz = nnz;
416 SI->nsos = nsos;
417 SI->nsos1 = nsos1;
418 SI->nsos2 = nsos2;
419 SI->nsosnz = nsosnz;
420 SI->nsosnz1 = nsosnz1;
421 SI->refd = refd = refd_copy(&SI->mysd[1], refd, n);
422 SI->v0 = v0;
423 SI->ve = ve;
424 SI->z = z;
425 SI->nzc = asl->i.nzc_;
426 asl->i.nzc_ += nnz;
427 ret:
428 M1free_ASL(&asl->i, M1state1, M1state2);
429 if (SI) {
430 SI->zap[0] = M1record(SI);
431 if (gp0)
432 SI->zap[1] = M1record(gp0);
433 if (v0)
434 SI->zap[2] = M1record(v0);
435 if (grefd)
436 SI->zap[3] = M1record(grefd->u.r);
437 if (refd)
438 SI->zap[4] = M1record(refd->u.r);
439 }
440 asl->i.nsuffixes = nss[0];
441 memcpy(asl->i.nsuff, nss+1, 4*sizeof(int));
442 memcpy(asl->i.suffixes, sds, 4*sizeof(SufDesc*));
443
444 return (Char*)SI;
445 }
446
447 static real LUge[2];
448
449 static cgrad**
newcon(Coninfo * CI,int ge)450 newcon(Coninfo *CI, int ge)
451 {
452 int m = CI->m++;
453 real *lu;
454 static real LU1[2] = { 0., 1. };
455
456 lu = ge ? LUge : LU1;
457 if (CI->cu) {
458 CI->clu[m] = lu[0];
459 CI-> cu[m] = lu[1];
460 }
461 else {
462 m <<= 1;
463 CI->clu[m] = lu[0];
464 CI->clu[m+1] = lu[1];
465 }
466 return CI->cgp++;
467 }
468
469 static void
newcoef(Coninfo * CI,cgrad *** p,int k,real t)470 newcoef(Coninfo *CI, cgrad ***p, int k, real t)
471 {
472 cgrad *cg = CI->cg++;
473 **p = cg;
474 *p = &cg->next;
475 cg->varno = k;
476 cg->coef = t;
477 }
478
479 static void
Bound(Coninfo * CI,int j,int k,real * LU)480 Bound(Coninfo *CI, int j, int k, real *LU)
481 {
482 cgrad **cgb;
483 if (LU[1]) {
484 cgb = newcon(CI, 1);
485 if (j < k) {
486 newcoef(CI, &cgb, j, -1.);
487 newcoef(CI, &cgb, k, LU[1]);
488 }
489 else {
490 newcoef(CI, &cgb, k, LU[1]);
491 newcoef(CI, &cgb, j, -1.);
492 }
493 *cgb = 0;
494 }
495 if (LU[0]) {
496 cgb = newcon(CI, 1);
497 if (j < k) {
498 newcoef(CI, &cgb, j, 1.);
499 newcoef(CI, &cgb, k, -LU[0]);
500 }
501 else {
502 newcoef(CI, &cgb, k, -LU[0]);
503 newcoef(CI, &cgb, j, 1.);
504 }
505 *cgb = 0;
506 }
507 }
508
509 static void
Bound2(Coninfo * CI,int j,int k0,int k,real * LU)510 Bound2(Coninfo *CI, int j, int k0, int k, real *LU)
511 {
512 cgrad **cgb;
513 if (LU[1]) {
514 cgb = newcon(CI, 1);
515 if (j < k0) {
516 newcoef(CI, &cgb, j, -1.);
517 newcoef(CI, &cgb, k0, LU[1]);
518 }
519 else {
520 newcoef(CI, &cgb, k0, LU[1]);
521 newcoef(CI, &cgb, j, -1.);
522 }
523 newcoef(CI, &cgb, k, LU[1]);
524 *cgb = 0;
525 }
526 if (LU[0]) {
527 cgb = newcon(CI, 1);
528 if (j < k0) {
529 newcoef(CI, &cgb, j, 1.);
530 newcoef(CI, &cgb, k0, -LU[0]);
531 }
532 else {
533 newcoef(CI, &cgb, k0, -LU[0]);
534 newcoef(CI, &cgb, j, 1.);
535 }
536 newcoef(CI, &cgb, k, -LU[0]);
537 *cgb = 0;
538 }
539 }
540
541 static void
debugchk(const char * what,int expected,int got,int exact)542 debugchk(const char *what, int expected, int got, int exact)
543 {
544 if (got != expected && (exact || got > expected)) {
545 fprintf(Stderr, "sos_finish: expected %s = %d, got %d\n",
546 what, expected, got);
547 mainexit_ASL(2);
548 }
549 }
550
551 int
sos_finish_ASL(ASL * asl,void ** VP,int flags,int * nsosnz_p,int ** sospri_p,int * copri,int ** sosbeg_p,int ** sosind_p,real ** sosref_p)552 sos_finish_ASL(ASL *asl, void **VP, int flags, int *nsosnz_p, int **sospri_p,
553 int *copri, int **sosbeg_p, int **sosind_p, real **sosref_p)
554 {
555 Char **vp;
556 Coninfo CI;
557 SOSinfo *SI;
558 SufDesc *grefd, *pd, *refd;
559 cde *Cde;
560 cde2 *Cde2;
561 cgrad **Cgrd, *cg, *cg0, *cg1, **cgb, **cgp0, **cgx;
562 const char *s;
563 expr_n *en;
564 int f, i, j, j0, k, k0, m, m0, n, niv0, nnc, nnv, nnv0, nnz, ns;
565 int nsos, nsos1, nsos2, nsosnz, nsosnz1, p1;
566 int *col1, *cs, *g, *ge, **gp, **gp0, **gp1, **gpe, *myp[2], *p;
567 int *pri, *sospri, *sosbeg, *sosbeg1, *sosind, *v, *v0, *ve, *z, **zg;
568 ograd *og;
569 real LU[2], LU1[2], *a, *lu, *mysr, *sosref, *sufref, *u, t, t1;
570 size_t L;
571
572 if (!VP || !(SI = (SOSinfo*)*VP))
573 return 0;
574 /* consistency check... */
575 if (asl != SI->asl) {
576 fprintf(Stderr, "Botched VP argument to sos_finish\n");
577 mainexit_ASL(1);
578 }
579 *VP = 0;
580 col1 = SI->col1;
581 gp0 = SI->gp0;
582 gp1 = SI->gp1;
583 gpe = SI->gpe;
584 grefd = SI->grefd;
585 nnc = SI->nnc;
586 nnv = SI->nnv;
587 nnz = SI->nnz;
588 nsos = SI->nsos;
589 nsos1 = SI->nsos1;
590 nsos2 = SI->nsos2;
591 nsosnz = SI->nsosnz;
592 nsosnz1 = SI->nsosnz1;
593 refd = SI->refd;
594 v0 = SI->v0;
595 ve = SI->ve;
596 z = SI->z;
597
598 if (pl_bigM <= 0.) {
599 if ((s = getenv("pl_bigM")))
600 pl_bigM = strtod(s,0);
601 if (pl_bigM <= 0.)
602 pl_bigM = 1e6;
603 mpl_bigM = -pl_bigM;
604 }
605
606 if (nsosnz_p)
607 *nsosnz_p = nsosnz;
608 L = nsos + 2*(nsos + 1)*sizeof(int)
609 + nsosnz*(sizeof(int) + sizeof(double));
610 if (!copri) {
611 if (sospri_p)
612 *sospri_p = 0;
613 sospri_p = 0;
614 }
615 else if (sospri_p)
616 L += nsos*sizeof(int);
617 mysr = 0;
618 Cgrd = Cgrad;
619 if (!sosref_p || !sosind_p || !sosbeg_p) {
620 if (sosref_p)
621 *sosref_p = 0;
622 if (sosind_p)
623 *sosind_p = 0;
624 if (sosbeg_p)
625 *sosbeg_p = 0;
626 if (sospri_p) {
627 *sospri_p = 0;
628 sospri_p = 0;
629 }
630 flags |= ASL_suf_sos_explict_free;
631 sosref_p = &mysr;
632 sosind_p = &myp[0];
633 sosbeg_p = &myp[1];
634 }
635 sosref = *sosref_p = flags & ASL_suf_sos_explict_free
636 ? (real*)Malloc(L)
637 : (real*)M1alloc(L);
638 sosind = *sosind_p = (int *)(sosref + nsosnz);
639 sosbeg = *sosbeg_p = sosind + nsosnz;
640 sosbeg1 = sosbeg + nsos + 1; /* scratch */
641 sospri = 0;
642 if (sospri_p)
643 sospri = *sospri_p = (int*)(sosbeg1 + nsos + 1);
644 memset(sosbeg, 0, (nsos+1)*sizeof(int));
645 f = Fortran;
646 cg0 = 0; /* silence buggy "not-initialized" warning */
647 cgp0 = 0; /* ditto */
648 n = n_var; /* ditto */
649 k0 = nnv0 = 0; /* ditto */
650 m0 = n_con;
651 if (nsos1) {
652 LUge[0] = 0.;
653 LUge[1] = Infinity;
654 CI.z = z;
655 CI.lu = LUv;
656 CI.u = Uvx;
657 CI.clu = LUrhs;
658 CI.cu = Urhsx;
659 CI.m = m0;
660 CI.nbin = 0;
661 if (Cgrd) {
662 CI.cgp = Cgrd + m0;
663 CI.cg = (cgrad*)M1alloc(nnz*sizeof(cgrad));
664 }
665 else {
666 CI.cgp = (cgrad**)Malloc(nnc*sizeof(cgrad*));
667 CI.cg = (cgrad*)Malloc(nnz*sizeof(cgrad));
668 }
669 cg0 = CI.cg;
670 cgp0 = CI.cgp;
671 pri = 0;
672 if (sospri) {
673 memset(sospri, 0, nsos1*sizeof(int));
674 if ((pd = suf_get("priority",
675 ASL_Sufkind_var | ASL_Sufkind_input)))
676 pri = pd->u.i;
677 }
678 i = 0;
679 for(gp = gp0; gp < gp1; gp++)
680 sosbeg[i++] = ((gp[1] - gp[0]) >> 1) - 1;
681 for(; gp < gpe; gp++)
682 sosbeg[i++] = (gp[1] - gp[0]) >> 1;
683 for(i = j = k = 0; i < nsos1; i++) { /* "k = 0" to omit an erroneous warning */
684 k = sosbeg[i] + j;
685 sosbeg[i] = j;
686 j = k;
687 }
688 sosbeg[nsos1] = k;
689 sufref = grefd->u.r;
690 if (A_vals) {
691 cs = A_colstarts;
692 k = cs[i = n];
693 j = n_var;
694 while(i < j)
695 cs[++i] = k;
696 }
697 p = get_vcmap_ASL(asl, ASL_Sufkind_var);
698 for(i = n, k = n + nnv; i < k; ++i)
699 p[i] = -1;
700 if (niv) {
701 /* tell write_sol about moved integer variables */
702
703 k = n + nnv;
704 p1 = n;
705 n -= niv;
706 while(--k > p1)
707 p[k] = p[k-nnv];
708 for(i = nnv; i > 0; --i)
709 p[k--] = -1;
710
711 /* copy integer variables up */
712
713 nnv0 = n + nnv;
714 for(i = n_obj; --i >= 0; )
715 for(og = Ograd[i]; og; og = og->next) {
716 if (og->varno >= n)
717 og->varno += nnv;
718 }
719 if (A_vals) {
720 cs = A_colstarts;
721 for(i = n + niv; --i >= n; )
722 cs[i+nnv] = cs[i];
723 j = cs[n];
724 for(i = n + nnv; i > n;)
725 cs[--i] = j;
726 }
727 else for(i = n_con - nnc; --i >= 0; )
728 for(cg = Cgrd[i]; cg; cg = cg->next) {
729 if (cg->varno >= n)
730 cg->varno += nnv;
731 }
732 }
733 niv0 = n;
734 /* add SOS2 stuff */
735 ns = 0;
736 for(gp = gp0; gp < gp1; gp++) {
737 cgx = newcon(&CI, 0);
738 g = gp[0];
739 ge = gp[1] - 4;
740 /* Order on ref row values, then */
741 /* average those values. */
742 qsortv(g, ((int)(ge-g)>>1) + 2, 2*sizeof(int),
743 refcomp, sufref);
744 for(; g <= ge; g += 2) {
745 j = g[1];
746 sufref[j] += 0.5*(sufref[g[3]] + sufref[j]);
747 }
748 for(g = gp[0]; g <= ge; g += 2) {
749 if ((i = g[1]) >= niv0)
750 i += nnv;
751 j = nonbinary(i, &CI, LU);
752 if (g == gp[0]) {
753 if (!j)
754 k = i;
755 else
756 Bound(&CI, i, k = n++, LU);
757 }
758 else if (g == ge) {
759 if (g[3] >= niv0)
760 g[3] += nnv;
761 if (nonbinary(g[3], &CI, LU1)) {
762 Bound2(&CI, i,k0,k=n++,LU);
763 Bound(&CI, g[3], k, LU1);
764 }
765 else
766 Bound(&CI, i, g[3], LU);
767 }
768 else
769 Bound2(&CI, i, k0, k = n++, LU);
770 newcoef(&CI, &cgx, k, 1.);
771 sosref[ns] = sufref[g[1]];
772 sosind[ns++] = k + f;
773 k0 = k;
774 }
775 *cgx = 0;
776 }
777 /* add SOS1 stuff */
778 for(; gp < gpe; gp++) {
779 cgx = newcon(&CI, 0);
780 for(g = gp[0], ge = gp[1]; g < ge; g += 2) {
781 if (g[1] >= niv0)
782 g[1] += nnv;
783 if (nonbinary(k = g[1], &CI, LU))
784 Bound(&CI, g[1], k = n++, LU);
785 newcoef(&CI, &cgx, k, 1.);
786 sosref[ns] = sufref[g[1]];
787 sosind[ns++] = k + f;
788 }
789 *cgx = 0;
790 }
791 if (pri)
792 for(i = 0, gp = gp0; gp < gpe; gp++) {
793 j = 0;
794 for(g = gp[0], ge = gp[1]; g < ge; g += 2) {
795 if (j < pri[g[1]])
796 j = pri[g[1]];
797 }
798 sospri[i++] = j;
799 }
800
801 /* reorder each SOS1 set by sosref */
802 p = SI->g0;
803 for(i = k = 0; i < nsos1; i++) {
804 j = j0 = k;
805 k = sosbeg[i];
806 t = sosref[j++];
807 while(j < k) {
808 t1 = sosref[j++];
809 if (t1 < t) {
810 reorder(sosind, sosref, j0, k, p);
811 break;
812 }
813 t = t1;
814 }
815 }
816
817 if (!nsos2)
818 goto sosbeg_adjust;
819 nsos -= nsos1;
820 sosbeg += nsos1;
821 *sosbeg = 0;
822 sosind += nsosnz1;
823 sosref += nsosnz1;
824 if (sospri)
825 sospri += nsos1;
826 }
827 for(v = col1; v < ve; )
828 if (((i = *v++) & 3) == 2 && !sosbeg[j = i>>4]++) {
829 if (sospri)
830 sospri[j-1] = copri[(i & 4) >> 2];
831 }
832 for(j = 0, i = 1; i <= nsos; i++) {
833 k = sosbeg[i] + j;
834 sosbeg1[i] = sosbeg[i] = j;
835 j = k;
836 }
837 sufref = refd->u.r;
838 for(v = col1; v < ve; v++)
839 if (((i = *v) & 3) == 2)
840 sosind[sosbeg[i>>4]++] = (v - v0) + f;
841 for(v = col1; v < ve; v++)
842 if ((i = *v) && !(i & 2)) {
843 j = i >> 4;
844 if ((k = sosbeg1[j]++) < sosbeg[j]) {
845 k0 = v - v0;
846 sosref[k] = sufref[k0] +
847 0.5*(sufref[k0+1] - sufref[k0]);
848 }
849 }
850 if (nsos1)
851 for(i = 0; i <= nsos; i++)
852 sosbeg[i] += nsosnz1;
853 nsos += nsos1;
854
855 sosbeg_adjust:
856 if (nsos1) {
857 k = (int)(CI.cg - cg0);
858 debugchk("n", n_var - niv - CI.nbin, n, 1);
859 debugchk("m", n_con, CI.m, 0);
860 debugchk("nnz", nnz, k, 0);
861
862 if (niv && CI.nbin) {
863
864 /* Correct overestimation of number */
865 /* of new binary variables. */
866
867 p = asl->i.vmap;
868 for(i = n_var - niv - CI.nbin, k = i + niv; i < k; ++i)
869 p[i] = p[i + CI.nbin];
870 for(i = n_obj; --i >= 0; )
871 for(og = Ograd[i]; og; og = og->next) {
872 if (og->varno >= nnv0)
873 og->varno -= CI.nbin;
874 }
875 if (A_vals) {
876 cs = A_colstarts;
877 for(i = n, j = n + niv; i <= j; i++)
878 cs[i] = cs[i+CI.nbin];
879 }
880 else for(i = n_con - nnc; --i >= 0; )
881 for(cg = Cgrd[i]; cg; cg = cg->next) {
882 if (cg->varno >= nnv0)
883 cg->varno -= CI.nbin;
884 }
885 }
886 if (CI.nbin) {
887 n_var -= CI.nbin;
888 c_vars -= CI.nbin;
889 o_vars -= CI.nbin;
890 nbv -= CI.nbin;
891 nnv -= CI.nbin;
892 if ((zg = zerograds)) {
893 j = n_var;
894 for(i = n_obj; i > 0; --i) {
895 for(p = *zg++; *p >= 0 && *p < j; p++);
896 *p = -1;
897 }
898 }
899 }
900
901 nnc = CI.m - m0;
902 n_conjac[1] = n_con = CI.m;
903 asl->i.nzc_ -= nnz - k;
904 nnz = k;
905
906 en = (expr_n *)mem_ASL(asl, sizeof(expr_n));
907 en->v = 0.;
908 en->op = (efunc_n *)f_OPNUM;
909 i = m0;
910 switch(asl->i.ASLtype) {
911 case ASL_read_f:
912 case ASL_read_fg:
913 Cde = ((ASL_fg*)asl)->I.con_de_;
914 goto more_Cde;
915 case ASL_read_pfg:
916 Cde = ((ASL_pfg*)asl)->I.con_de_;
917 more_Cde:
918 while(i < CI.m)
919 Cde[i++].e = (expr*)en;
920 break;
921 case ASL_read_fgh:
922 Cde2 = ((ASL_fgh*)asl)->I.con2_de_;
923 goto more_Cde2;
924 case ASL_read_pfgh:
925 Cde2 = ((ASL_pfgh*)asl)->I.con2_de_;
926 more_Cde2:
927 while(i < CI.m)
928 Cde2[i++].e = (expr2*)en;
929 }
930 j = asl->i.n_con1;
931 asl->i.n_con1 = j + nnc;
932 if ((p = asl->i.cmap))
933 for(i = m0; i < CI.m; ++i)
934 p[i] = j++;
935 }
936 if (f)
937 for(sosbeg = *sosbeg_p, i = 0; i <= nsos; i++)
938 sosbeg[i] += f;
939
940 if (!nnc) /* no new constraints */
941 goto freeup;
942
943 if (nnv) {
944 j = asl->i.n_var1;
945 asl->i.n_var1 = j + nnv;
946 if ((p = asl->i.vminv)) {
947 i = n_var;
948 k = i + nnv;
949 do p[i++] = j++; while(j < k);
950 }
951 lu = LUv;
952 k = n; /* n == n_var - niv */
953 i = n + niv;
954 if ((u = Uvx)) {
955 while(i > k) {
956 --i;
957 lu[i] = lu[i-nnv];
958 u[i] = u[i-nnv];
959 }
960 for(k -= nnv; --i >= k; ) {
961 lu[i] = 0.;
962 u[i] = 1.;
963 }
964 }
965 else {
966 for(i <<= 1, j = nnv << 1, k <<= 1; --i >= k;)
967 lu[i] = lu[i-j];
968 k -= 2*nnv;
969 do {
970 lu[i--] = 1.;
971 lu[i--] = 0.;
972 }
973 while(i >= k);
974 }
975 }
976 n = n_var;
977 if (Cgrd) {
978
979 /* (re)compute goff fields */
980
981 m = m0 + nnc;
982 n += nnv;
983 p = (int*)Malloc(L = n*sizeof(int));
984 memset(p, 0, L);
985 for(i = 0; i < m; i++)
986 for(cg = Cgrd[i]; cg; cg = cg->next)
987 p[cg->varno]++;
988 for(i = k = 0; i < n; i++) {
989 j = p[i];
990 p[i] = k;
991 k += j;
992 }
993 for(i = 0; i < m; i++)
994 for(cg = Cgrd[i]; cg; cg = cg->next)
995 cg->goff = p[cg->varno]++;
996 free(p);
997 }
998 else {
999 /* insert new nonzeros into A_val, A_rownos, A_colstarts */
1000
1001 cgb = (cgrad**)Malloc(L = n*sizeof(cgrad*));
1002 memset(cgb, 0, L);
1003 for(i = 0; i < nnc; i++)
1004 for(cg = cgp0[i]; cg; cg = cg1) {
1005 cg1 = cg->next;
1006 j = cg->varno;
1007 cg->varno = i + m0;
1008 cg->next = cgb[j];
1009 cgb[j] = cg;
1010 }
1011 cs = A_colstarts;
1012 p = A_rownos;
1013 a = A_vals;
1014 k = asl->i.nzc_;
1015 for(i = n; i > 0; ) {
1016 j = cs[i];
1017 cs[i--] = k + f;
1018 for(cg = cgb[i]; cg; cg = cg->next) {
1019 a[--k] = cg->coef;
1020 p[k] = cg->varno + f;
1021 --nnz;
1022 }
1023 if (nnz <= 0)
1024 break;
1025 j -= cs[i];
1026 while(j-- > 0) {
1027 j0 = --k - nnz;
1028 a[k] = a[j0];
1029 p[k] = p[j0];
1030 }
1031 }
1032 free(cgb);
1033 free(cg0);
1034 free(cgp0);
1035 }
1036 freeup:
1037 if (mysr)
1038 free(mysr);
1039 for(i = 5; i--; )
1040 if ((vp = SI->zap[i])) {
1041 free(*vp);
1042 *vp = 0;
1043 }
1044 n_var += nnv;
1045 c_vars += nnv;
1046 o_vars += nnv;
1047 nbv += nnv;
1048 n_con += nnc;
1049 return nsos;
1050 }
1051