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