1 /****************************************************************
2 Copyright (C) 1997-2001 Lucent Technologies
3 All Rights Reserved
4
5 Permission to use, copy, modify, and distribute this software and
6 its documentation for any purpose and without fee is hereby
7 granted, provided that the above copyright notice appear in all
8 copies and that both that the copyright notice and this
9 permission notice and warranty disclaimer appear in supporting
10 documentation, and that the name of Lucent or any of its entities
11 not be used in advertising or publicity pertaining to
12 distribution of the software without specific, written prior
13 permission.
14
15 LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
16 INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
17 IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
18 SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
19 WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
20 IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
21 ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
22 THIS SOFTWARE.
23 ****************************************************************/
24
25 #include "nlp.h"
26
27 #define Egulp 400
28 #define GAP_MAX 10
29
30 #ifdef Just_Linear
31 #define who "f_read"
32 #define fg_read_ASL f_read_ASL
33 #else
34 #define who "fg_read"
35 #ifdef __cplusplus
36 extern "C" {
37 static real Missing_func(arglist*);
38 static int compar(const void*, const void*, void*);
39 }
40 #endif /* __cplusplus */
41 #endif /* Just_Linear */
42
43 #undef nzc
44 #undef r_ops
45
46 typedef struct
47 Static {
48 int _k_seen, _nv0;
49 ASL *a;
50 ASL_fg *asl;
51 efunc **_r_ops;
52 #ifndef Just_Linear
53 derp *_last_d;
54 expr *(*_holread) ANSI((EdRead*));
55 expr_if *_iflist, *_if2list, *_if2list_end;
56 expr_va *_varglist, *_varg2list, *_varg2list_end;
57 relo *_relolist, *_relo2list;
58 int *_imap, *_vrefnext, *_vrefx, *_zc, *_zci;
59 int _amax1, _co_first, _firstc1, _imap_len;
60 int _last_cex, _lasta, _lasta0, _lasta00, _lastc1, _lastj;
61 int _max_var, _ncom_togo, _nderp, _nocopy;
62 int _nv01, _nv011, _nv0b, _nv0c, _nv1, _nvref, _nzc, _nzclim;
63 int nvar0, nvinc;
64 #endif /* Just_Linear */
65 } Static;
66
67 #define amax1 S->_amax1
68 #define co_first S->_co_first
69 #define firstc1 S->_firstc1
70 #define holread S->_holread
71 #define if2list S->_if2list
72 #define if2list_end S->_if2list_end
73 #define iflist S->_iflist
74 #define imap S->_imap
75 #define imap_len S->_imap_len
76 #define k_seen S->_k_seen
77 #define last_cex S->_last_cex
78 #define last_d S->_last_d
79 #define lasta S->_lasta
80 #define lasta0 S->_lasta0
81 #define lasta00 S->_lasta00
82 #define lastc1 S->_lastc1
83 #define lastj S->_lastj
84 #define max_var S->_max_var
85 #define ncom_togo S->_ncom_togo
86 #define nderp S->_nderp
87 #define nocopy S->_nocopy
88 #define nv0 S->_nv0
89 #define nv01 S->_nv01
90 #define nv011 S->_nv011
91 #define nv0b S->_nv0b
92 #define nv0c S->_nv0c
93 #define nv1 S->_nv1
94 #define nvref S->_nvref
95 #define nzc S->_nzc
96 #define nzclim S->_nzclim
97 #define r_ops S->_r_ops
98 #define relo2list S->_relo2list
99 #define relolist S->_relolist
100 #define varg2list S->_varg2list
101 #define varg2list_end S->_varg2list_end
102 #define varglist S->_varglist
103 #define vrefnext S->_vrefnext
104 #define vrefx S->_vrefx
105 #define zc S->_zc
106 #define zci S->_zci
107
108 #ifdef Just_Linear
109
110 static void
111 #ifdef KR_headers
sorry_nonlin(R)112 sorry_nonlin(R) EdRead *R;
113 #else
114 sorry_nonlin(EdRead *R)
115 #endif
116 {
117 fprintf(Stderr,
118 "Sorry, %s can't handle nonlinearities.\n",
119 progname ? progname : "");
120 exit_ASL(R,ASL_readerr_nonlin);
121 }
122 #else /* Just_Linear */
123
124 #include "r_opn.hd"
125 #endif /* Just_Linear */
126
127 static void
128 #ifdef KR_headers
sorry_CLP(R,what)129 sorry_CLP(R, what) EdRead *R; char *what;
130 #else
131 sorry_CLP(EdRead *R, char *what)
132 #endif
133 {
134 fprintf(Stderr,
135 "Sorry, %s cannot handle %s.\n",
136 progname ? progname : "", what);
137 exit_ASL(R,ASL_readerr_CLP);
138 }
139
140 static Static *
141 #ifdef KR_headers
ed_reset(S,a)142 ed_reset(S, a) Static *S; ASL *a;
143 #else
144 ed_reset(Static *S, ASL *a)
145 #endif
146 {
147 memset(S, 0, sizeof(Static));
148 S->asl = (ASL_fg*)a;
149 S->a = a;
150 a->i.memLast = a->i.memNext = 0;
151 #ifndef Just_Linear
152 co_first = 1;
153 #endif
154 return S;
155 }
156
157 #ifdef Double_Align
158 #define memadj(x) x
159 #else
160 #define memadj(x) (((x) + (sizeof(long)-1)) & ~(sizeof(long)-1))
161 #endif
162
163 #ifndef Just_Linear
164
165 static void
166 #ifdef KR_headers
fscream(R,name,nargs,kind)167 fscream(R, name, nargs, kind) EdRead *R; char *name; char *kind;
168 #else
169 fscream(EdRead *R, const char *name, int nargs, char *kind)
170 #endif
171 {
172 scream(R, ASL_readerr_argerr,
173 "line %ld: attempt to call %s with %d %sargs\n",
174 R->Line, name, nargs, kind);
175 }
176
177 static void
178 #ifdef KR_headers
new_derp(S,a,b,c)179 new_derp(S, a, b, c) Static *S; real *c;
180 #else
181 new_derp(Static *S, int a, int b, real *c)
182 #endif
183 {
184 derp *d;
185 if (a == nv1)
186 return;
187 nderp++;
188 d = (derp *)mem_ASL(S->a, sizeof(derp));
189 d->next = last_d;
190 last_d = d;
191 d->a.i = a;
192 d->b.i = b;
193 d->c.rp = c;
194 }
195
196 static derp *
197 #ifdef KR_headers
new_relo(S,e,Dnext,ap)198 new_relo(S, e, Dnext, ap) Static *S; expr *e; derp *Dnext; int *ap;
199 #else
200 new_relo(Static *S, expr *e, derp *Dnext, int *ap)
201 #endif
202 {
203 relo *r;
204 derp *d;
205
206 r = (relo *)mem_ASL(S->a, sizeof(relo));
207 r->next = relolist;
208 r->next2 = relo2list;
209 relo2list = relolist = r;
210 if (last_d != Dnext) {
211 *ap = e->a;
212 for(d = last_d; d->next != Dnext; d = d->next);
213 d->next = 0;
214 }
215 else {
216 last_d = 0;
217 new_derp(S, e->a, *ap = lasta++, &edagread_one);
218 }
219 r->D = r->Dcond = last_d;
220 r->Dnext = Dnext;
221 return r->D;
222 }
223
224 static relo *
225 #ifdef KR_headers
new_relo1(S,Dnext)226 new_relo1(S, Dnext) Static *S; derp *Dnext;
227 #else
228 new_relo1(Static *S, derp *Dnext)
229 #endif
230 {
231 relo *r;
232
233 r = (relo *)mem_ASL(S->a, sizeof(relo));
234 r->next = relolist;
235 relolist = r;
236 r->D = 0;
237 r->Dnext = Dnext;
238 return r;
239 }
240
241 static expr *
242 #ifdef KR_headers
new_expr(S,opcode,L,R,deriv)243 new_expr(S, opcode, L, R, deriv) Static *S; expr *L, *R;
244 #else
245 new_expr(Static *S, int opcode, expr *L, expr *R, int deriv)
246 #endif
247 {
248 expr *rv;
249 efunc *o;
250 int L1, R1;
251
252 o = r_ops[opcode];
253 if (o == f_OPPOW)
254 if (R->op == f_OPNUM)
255 if (((expr_n *)R)->v == 2.) {
256 o = f_OP2POW;
257 R = 0;
258 }
259 else
260 o = f_OP1POW;
261 else if (L->op == f_OPNUM)
262 o = f_OPCPOW;
263 rv = (expr *)mem_ASL(S->a, sizeof(expr));
264 rv->op = o;
265 rv->L.e = L;
266 rv->R.e = R;
267 rv->a = nv1;
268 if (deriv) {
269 L1 = L && L->op != f_OPNUM && L->a != nv1;
270 R1 = R && R->op != f_OPNUM && R->a != nv1;
271 if (L1 | R1) {
272 rv->a = lasta++;
273 if (L1)
274 new_derp(S, L->a, rv->a, &rv->dL);
275 if (R1)
276 new_derp(S, R->a, rv->a, &rv->dR);
277 }
278 }
279 return rv;
280 }
281
282 #endif /* Just_Linear */
283
284 static expr *
285 #ifdef KR_headers
eread(R,deriv)286 eread(R, deriv) EdRead *R; int deriv;
287 #else
288 eread(EdRead *R, int deriv)
289 #endif
290 {
291 expr_n *rvn;
292 short sh;
293 fint L1;
294 real r;
295 Static *S = (Static *)R->S;
296 ASL_fg *asl = S->asl;
297 #ifndef Just_Linear
298 char **sa;
299 int a0, a1, *at, i, j, k, kd, ks, numargs, symargs;
300 real *b, *ra;
301 expr *L, *arg, **args, **argse, *rv;
302 expr_va *rva;
303 plterm *p;
304 de *d;
305 derp *dsave;
306 efunc *op;
307 expr_if *rvif;
308 expr_f *rvf;
309 func_info *fi;
310 arglist *al;
311 argpair *ap, *sap;
312 char *dig;
313 static real dvalue[] = {
314 #include "dvalue.hd"
315 };
316 #else /* Just_Linear */
317 Not_Used(deriv);
318 #endif /* Just_Linear */
319
320 switch(edag_peek(R)) {
321 #ifdef Just_Linear
322 case 'f':
323 case 'h':
324 case 'o':
325 case 'v':
326 sorry_nonlin(R);
327 #else
328 case 'f':
329 if (xscanf(R, "%d %d", &i, &j) != 2)
330 badline(R);
331 fi = funcs[i];
332 if (fi->nargs >= 0) {
333 if (fi->nargs != j) {
334 bad_nargs:
335 fscream(R, fi->name, j, "");
336 }
337 }
338 else if (-(1+fi->nargs) > j)
339 goto bad_nargs;
340 rvf = (expr_f *)mem(sizeof(expr_f)
341 + (j-1)*sizeof(expr *));
342 rvf->op = f_OPFUNCALL;
343 rvf->fi = fi;
344 args = rvf->args;
345 argse = args + j;
346 k = ks = symargs = numargs = 0;
347 while(args < argse) {
348 arg = *args++ = eread(R, deriv);
349 if ((op = arg->op) == f_OPHOL)
350 symargs++;
351 else if (op == f_OPIFSYM)
352 ks++;
353 else {
354 numargs++;
355 if (op != f_OPNUM)
356 k++;
357 }
358 }
359 symargs += ks;
360 if (symargs && !(fi->ftype & 1))
361 fscream(R, fi->name, symargs, "symbolic ");
362 if (deriv && k) {
363 kd = numargs;
364 rvf->a = lasta++;
365 }
366 else {
367 kd = 0;
368 rvf->a = nv1;
369 }
370 ra = (real *)mem(sizeof(arglist)
371 + (k+ks)*sizeof(argpair)
372 + (numargs+kd)*sizeof(real)
373 + symargs*sizeof(char *)
374 + j*sizeof(int));
375 dig = 0;
376 if (k < kd)
377 dig = (char*)mem(numargs);
378 b = ra + kd;
379 al = rvf->al = (arglist *)(b + numargs);
380 al->n = numargs + symargs;
381 al->nr = numargs;
382 al->ra = ra;
383 if (kd)
384 memset(al->derivs = b, 0, kd*sizeof(real));
385 else
386 al->derivs = 0;
387 al->hes = 0;
388 al->dig = dig;
389 al->funcinfo = fi->funcinfo;
390 al->AE = asl->i.ae;
391 al->sa = (Const char**)(sa = (char **)(al + 1));
392 ap = rvf->ap = (argpair *)(sa + symargs);
393 sap = rvf->sap = ap + k;
394 at = al->at = (int *)(sap + ks);
395 symargs = numargs = 0;
396 for(args = rvf->args; args < argse; at++) {
397 arg = *args++;
398 if ((op = arg->op) == f_OPHOL) {
399 *at = --symargs;
400 *sa++ = ((expr_h *)arg)->sym;
401 }
402 else if (op == f_OPIFSYM) {
403 *at = --symargs;
404 sap->e = arg;
405 (sap++)->u.s = sa++;
406 }
407 else {
408 *at = numargs++;
409 j = 1;
410 if (op == f_OPNUM)
411 *ra = ((expr_n *)arg)->v;
412 else {
413 ap->e = arg;
414 (ap++)->u.v = ra;
415 if (kd) {
416 j = 0;
417 new_derp(S, arg->a,
418 rvf->a, b);
419 *b = 0;
420 }
421 }
422 if (dig)
423 *dig++ = j;
424 b++;
425 ra++;
426 }
427 }
428 rvf->ape = ap;
429 rvf->sape = sap;
430 return (expr *)rvf;
431
432 case 'h':
433 return holread(R);
434 #endif /* Just_Linear */
435 case 's':
436 if (xscanf(R, "%hd", &sh) != 1)
437 badline(R);
438 r = sh;
439 goto have_r;
440
441 case 'l':
442 if (xscanf(R, "%ld", &L1) != 1)
443 badline(R);
444 r = L1;
445 goto have_r;
446
447 case 'n':
448 if (xscanf(R, "%lf", &r) != 1)
449 badline(R);
450 have_r:
451 rvn = (expr_n *)mem(size_expr_n);
452 rvn->op = (efunc_n*)f_OPNUM;
453 rvn->v = r;
454 return (expr *)rvn;
455 #ifndef Just_Linear
456
457 case 'o':
458 break;
459
460 case 'v':
461 if (xscanf(R,"%d",&k) != 1 || k < 0)
462 badline(R);
463 if (k >= S->nvar0)
464 k += S->nvinc;
465 if (k > max_var)
466 badline(R);
467 if (k < nv01 && deriv && !zc[k]++)
468 zci[nzc++] = k;
469 return (expr *)(var_e + k);
470
471 #endif /* Just_Linear */
472 default:
473 badline(R);
474 }
475
476 #ifndef Just_Linear
477
478 if (xscanf(R, "%d", &k) != 1 || k < 0 || k >= N_OPS)
479 badline(R);
480 switch(optype[k]) {
481
482 case 1: /* unary */
483 rv = new_expr(S, k, eread(R, deriv), 0, deriv);
484 rv->dL = dvalue[k]; /* for UMINUS, FLOOR, CEIL */
485 return rv;
486
487 case 2: /* binary */
488 if (dvalue[k] == 11)
489 deriv = 0;
490 L = eread(R, deriv);
491 rv = new_expr(S, k, L, eread(R, deriv), deriv);
492 rv->dL = 1.;
493 rv->dR = dvalue[k]; /* for PLUS, MINUS, REM */
494 return rv;
495
496 case 3: /* vararg (min, max) */
497 i = -1;
498 xscanf(R, "%d", &i);
499 if (i <= 0)
500 badline(R);
501 rva = (expr_va *)mem(sizeof(expr_va));
502 rva->op = r_ops[k];
503 rva->L.d = d = (de *)mem(i*sizeof(de) + sizeof(expr *));
504 rva->next = varglist;
505 varglist = varg2list = rva;
506 if (!last_d && deriv) {
507 new_derp(S, lasta, lasta, &edagread_one);
508 lasta++;
509 }
510 rva->d0 = dsave = last_d;
511 a0 = a1 = lasta;
512 for(j = 0; i > 0; i--, d++) {
513 last_d = dsave;
514 d->e = L = eread(R, deriv);
515 if (L->op == f_OPNUM || L->a == nv1 || !deriv) {
516 d->d = dsave;
517 d->dv.i = nv1;
518 }
519 else {
520 d->d = new_relo(S, L, dsave, &d->dv.i);
521 j++;
522 if (a1 < lasta)
523 a1 = lasta;
524 lasta = a0;
525 }
526 }
527 d->e = 0; /* sentinnel expr * */
528 last_d = dsave;
529 if (j) {
530 rva->a = lasta = a1;
531 new_derp(S, 0, lasta++, &edagread_one);
532 /* f_MINLIST or f_MAXLIST will replace the 0 */
533 rva->R.D = last_d;
534 nocopy = 1;
535 }
536 else {
537 rva->a = nv1;
538 rva->R.D = 0;
539 }
540 return (expr *)rva;
541
542 case 4: /* piece-wise linear */
543 i = -1;
544 xscanf(R, "%d", &i);
545 if (i <= 1)
546 badline(R);
547 plterms++;
548 j = 2*i - 1;
549 p = (plterm *)mem(sizeof(plterm) + (j-1)*sizeof(real));
550 p->n = i;
551 b = p->bs;
552 do {
553 switch(edag_peek(R)) {
554 case 's':
555 if (xscanf(R,"%hd",&sh) != 1)
556 badline(R);
557 r = sh;
558 break;
559 case 'l':
560 if (xscanf(R,"%ld",&L1) != 1)
561 badline(R);
562 r = L1;
563 break;
564 case 'n':
565 if (xscanf(R,"%lf",&r) == 1)
566 break;
567 default:
568 badline(R);
569 }
570 *b++ = r;
571 }
572 while(--j > 0);
573 rv = (expr *)mem(sizeof(expr));
574 rv->op = f_OPPLTERM;
575 rv->L.p = p;
576 rv->R.e = L = eread(R, deriv);
577 if (deriv)
578 new_derp(S, L->a, rv->a = lasta++, &rv->dL);
579 return rv;
580
581 case 5: /* if */
582 rvif = (expr_if *)mem(sizeof(expr_if));
583 rvif->op = r_ops[k];
584 rvif->next = iflist;
585 iflist = if2list = rvif;
586 if (!last_d && deriv) {
587 new_derp(S, lasta, lasta, &edagread_one);
588 lasta++;
589 }
590 rvif->d0 = dsave = last_d;
591 rvif->e = eread(R, 0);
592 a0 = lasta;
593 rvif->T = L = eread(R, deriv);
594 j = 0;
595 if (L->op == f_OPNUM || L->a == nv1 || !(j = deriv)) {
596 rvif->dT = dsave;
597 rvif->Tv.i = nv1;
598 }
599 else
600 rvif->dT = new_relo(S, L, dsave, &rvif->Tv.i);
601 a1 = lasta;
602 lasta = a0;
603 last_d = dsave;
604 rvif->F = L = eread(R, deriv);
605 if (L->op == f_OPNUM || L->a == nv1 || !(j = deriv)) {
606 rvif->dF = dsave;
607 rvif->Fv.i = nv1;
608 }
609 else
610 rvif->dF = new_relo(S, L, dsave, &rvif->Fv.i);
611 if (lasta < a1)
612 lasta = a1;
613 last_d = dsave;
614 if (j) {
615 new_derp(S, 0, rvif->a = lasta++, &edagread_one);
616 rvif->D = last_d;
617 nocopy = 1;
618 }
619 else {
620 rvif->a = nv1;
621 rvif->D = 0;
622 }
623 return (expr *)rvif;
624
625 case 11: /* OPCOUNT */
626 deriv = 0;
627 /* no break */
628 case 6: /* sumlist */
629 i = j = 0;
630 xscanf(R, "%d", &i);
631 if (i <= 2 && (optype[k] == 6 || i < 1))
632 badline(R);
633 rv = (expr *)mem(sizeof(expr) - sizeof(real)
634 + i*sizeof(expr *));
635 rv->op = r_ops[k];
636 a0 = lasta;
637 rv->L.ep = args = (expr **)&rv->dR;
638 if (deriv) {
639 rv->a = lasta++;
640 do {
641 *args++ = L = eread(R, deriv);
642 if (L->op != f_OPNUM && L->a != nv1) {
643 new_derp(S, L->a, rv->a,
644 &edagread_one);
645 j++;
646 }
647 }
648 while(--i > 0);
649 }
650 else do
651 *args++ = eread(R, deriv);
652 while(--i > 0);
653 rv->R.ep = args;
654 if (!j) {
655 rv->a = nv1;
656 lasta = a0;
657 }
658 return rv;
659 }
660 #endif /* Just_Linear */
661 badline(R);
662 return 0;
663 }
664
665 #ifndef Just_Linear
666
667 static list *
668 #ifdef KR_headers
new_list(asl,nxt)669 new_list(asl, nxt) ASL *asl; list *nxt;
670 #else
671 new_list(ASL *asl, list *nxt)
672 #endif
673 {
674 list *rv = (list *)mem(sizeof(list));
675 rv->next = nxt;
676 return rv;
677 }
678
679 static list *
680 #ifdef KR_headers
crefs(S)681 crefs(S) Static *S;
682 #else
683 crefs(Static *S)
684 #endif
685 {
686 int i;
687 list *rv = 0;
688
689 while(nzc > 0) {
690 if ((i = zci[--nzc]) >= nv0) {
691 rv = new_list(S->a, rv);
692 rv->item.i = i;
693 }
694 zc[i] = 0;
695 }
696 return rv;
697 }
698
699 static funnel *
700 #ifdef KR_headers
funnelfix(f)701 funnelfix(f) funnel *f;
702 #else
703 funnelfix(funnel *f)
704 #endif
705 {
706 cexp *ce;
707 funnel *fnext, *fprev;
708
709 for(fprev = 0; f; f = fnext) {
710 fnext = f->next;
711 f->next = fprev;
712 fprev = f;
713 ce = f->ce;
714 ce->z.i = ce->d->b.i;
715 }
716 return fprev;
717 }
718
719 static derp *
720 #ifdef KR_headers
derpadjust(S,d0,a,dnext)721 derpadjust(S, d0, a, dnext) Static *S; derp *d0; int a; derp *dnext;
722 #else
723 derpadjust(Static *S, derp *d0, int a, derp *dnext)
724 #endif
725 {
726 derp *d, *d1;
727 int *r, *re;
728 relo *rl;
729 expr_if *il, *ile;
730 expr_va *vl, *vle;
731 de *de1;
732 ASL_fg *asl;
733
734 if (!(d = d0))
735 return dnext;
736 asl = S->asl;
737 r = imap + lasta0;
738 re = imap + lasta;
739 while(r < re)
740 *r++ = a++;
741 if (amax < a)
742 amax = a;
743 r = imap;
744 for(;; d = d1) {
745 d->a.i = r[d->a.i];
746 d->b.i = r[d->b.i];
747 if (!(d1 = d->next))
748 break;
749 }
750 d->next = dnext;
751 if (rl = relo2list) {
752 relo2list = 0;
753 do {
754 d = rl->Dcond;
755 do {
756 d->a.i = r[d->a.i];
757 d->b.i = r[d->b.i];
758 }
759 while(d = d->next);
760 }
761 while(rl = rl->next2);
762 }
763 if (if2list != if2list_end) {
764 ile = if2list_end;
765 if2list_end = il = if2list;
766 do {
767 il->Tv.i = r[il->Tv.i];
768 il->Fv.i = r[il->Fv.i];
769 }
770 while((il = il->next) != ile);
771 }
772 if (varg2list != varg2list_end) {
773 vle = varg2list_end;
774 varg2list_end = vl = varg2list;
775 do {
776 for(de1 = vl->L.d; de1->e; de1++)
777 de1->dv.i = r[de1->dv.i];
778 }
779 while((vl = vl->next) != vle);
780 }
781 return d0;
782 }
783
784 static derp *
785 #ifdef KR_headers
derpcopy(S,ce,dnext)786 derpcopy(S, ce, dnext) Static *S; cexp *ce; derp *dnext;
787 #else
788 derpcopy(Static *S, cexp *ce, derp *dnext)
789 #endif
790 {
791 derp *d, *dprev;
792 int *map;
793 derp d00;
794
795 if (!(d = ce->d))
796 return dnext;
797 map = imap;
798 for(dprev = &d00; d; d = d->next) {
799 new_derp(S, map[d->a.i], map[d->b.i], d->c.rp);
800 dprev = dprev->next = last_d;
801 }
802 dprev->next = dnext;
803 return d00.next;
804 }
805
806 static void
807 #ifdef KR_headers
imap_alloc(S)808 imap_alloc(S) Static *S;
809 #else
810 imap_alloc(Static *S)
811 #endif
812 {
813 int i, *r, *re;
814
815 if (imap) {
816 imap_len += lasta;
817 imap = (int *)Realloc(imap, imap_len*Sizeof(int));
818 return;
819 }
820 imap_len = amax1 > lasta ? amax1 : lasta;
821 imap_len += 100;
822 r = imap = (int *)Malloc(imap_len*Sizeof(int));
823 for(i = 0, re = r + nv1+1; r < re;)
824 *r++ = i++;
825 }
826
827 static int
828 #ifdef KR_headers
compar(a,b,v)829 compar(a, b, v) char *a, *b, *v;
830 #else
831 compar(const void *a, const void *b, void *v)
832 #endif
833 {
834 Not_Used(v);
835 return *(int*)a - *(int*)b;
836 }
837
838 static void
839 #ifdef KR_headers
comsubs(S,alen,d,z)840 comsubs(S, alen, d, z) Static *S; cde *d; int alen, **z;
841 #else
842 comsubs(Static *S, int alen, cde *d, int **z)
843 #endif
844 {
845 list *L;
846 int a, i, j, k;
847 int *r, *re, *z1;
848 cexp *ce;
849 derp *D, *dnext;
850 relo *R;
851 ASL_fg *asl = S->asl;
852
853 D = last_d;
854 a = lasta00;
855 dnext = 0;
856 R = 0;
857 for(i = j = 0; i < nzc; i++)
858 if ((k = zci[i]) >= nv0)
859 zci[j++] = k;
860 else
861 zc[k] = 0;
862 if (nzc = j) {
863 for(i = 0; i < nzc; i++)
864 for(L = cexps[zci[i]-nv0].cref; L; L = L->next)
865 if (!zc[L->item.i]++)
866 zci[nzc++] = L->item.i;
867 if (nzc > 1)
868 if (nzc < nzclim)
869 qsortv(zci, nzc, sizeof(int), compar, NULL);
870 else for(i = nv0, j = 0; i < max_var; i++)
871 if (zc[i])
872 zci[j++] = i;
873 }
874 if (z && (k = lastc1 - firstc1 + nzc)) {
875 i = (2*k + 1)*sizeof(int);
876 *z = z1 = k > 20 ? (int *)M1alloc(i) : (int *)mem(i);
877 *z1++ = k;
878 }
879 if (nzc > 0) {
880 R = new_relo1(S, dnext);
881 i = 0;
882 do {
883 j = zci[i];
884 zc[j] = 0;
885 ce = &cexps[j - nv0];
886 if (ce->funneled)
887 imap[var_e[j].a] = a++;
888 else {
889 r = imap + ce->z.i;
890 re = r + ce->zlen;
891 while(r < re)
892 *r++ = a++;
893 }
894 if (z) {
895 *z1++ = j;
896 *z1++ = a - 1;
897 }
898 dnext = R->D = derpcopy(S, ce, R->D);
899 }
900 while(++i < nzc);
901 nzc = 0;
902 }
903 if (D || R) {
904 if (!R)
905 R = new_relo1(S, dnext);
906 D = R->D = derpadjust(S, D, a, R->D);
907 if (d->e->op != f_OPVARVAL)
908 d->e->a = imap[d->e->a];
909 }
910 if (z)
911 for(i = firstc1 + nv01, j = lastc1 + nv01; i < j; i++) {
912 *z1++ = i;
913 *z1++ = imap[var_e[i].a];
914 }
915 d->d = D;
916 a += alen;
917 d->zaplen = (a > lasta00 ? a - nv1 : 0)*sizeof(real);
918 if (amax < a)
919 amax = a;
920 }
921 #endif /* Just_Linear */
922
923 static void
924 #ifdef KR_headers
co_read(R,d,cexp1_end,k,z,wd)925 co_read(R, d, cexp1_end, k, z, wd) EdRead *R; cde *d; int *cexp1_end, k, wd, **z;
926 #else
927 co_read(EdRead *R, cde *d, int *cexp1_end, int k, int **z, int wd)
928 #endif
929 {
930 #ifdef Just_Linear
931 Not_Used(cexp1_end);
932 Not_Used(z);
933 #else
934 Static *S = (Static *)R->S;
935 ASL_fg *asl = S->asl;
936 lastc1 = last_cex - nv011;
937 if (cexp1_end)
938 cexp1_end[k+1] = lastc1;
939 if (amax1 < lasta)
940 amax1 = lasta;
941 if (co_first) {
942 co_first = 0;
943 if (imap_len < lasta)
944 imap_alloc(S);
945 f_b = funnelfix(f_b);
946 f_c = funnelfix(f_c);
947 f_o = funnelfix(f_o);
948 }
949 if (!lastj) {
950 lasta = lasta0;
951 last_d = 0;
952 }
953 lastj = 0;
954 #endif /* Just_Linear */
955 d += k;
956 d->e = eread(R, wd);
957 #ifndef Just_Linear
958 { int alen;
959 alen = lasta - lasta0;
960 if (imap_len < lasta)
961 imap_alloc(S);
962 if (z) {
963 z += k;
964 *z = 0;
965 }
966 comsubs(S, alen, d, z);
967 firstc1 = lastc1;
968 }
969 #endif /* Just_Linear */
970 }
971
972 #ifndef Just_Linear
973 static linpart *
974 #ifdef KR_headers
linpt_read(R,nlin)975 linpt_read(R, nlin) EdRead *R;
976 #else
977 linpt_read(EdRead *R, int nlin)
978 #endif
979 {
980 linpart *L, *rv;
981 ASL *asl = R->asl;
982
983 if (nlin <= 0)
984 return 0;
985 L = rv = (linpart *)mem(nlin*sizeof(linpart));
986 do {
987 if (xscanf(R, "%d %lf", &L->v.i, &L->fac) != 2)
988 badline(R);
989 L++;
990 }
991 while(--nlin > 0);
992 return rv;
993 }
994
995 static int
996 #ifdef KR_headers
funnelkind(S,ce,ip)997 funnelkind(S, ce, ip) Static *S; cexp *ce; int *ip;
998 #else
999 funnelkind(Static *S, cexp *ce, int *ip)
1000 #endif
1001 {
1002 int i, j, k, nzc0, rv;
1003 int *vr, *vre;
1004 ASL_fg *asl = S->asl;
1005
1006 ce->vref = 0;
1007 if (!(nzc0 = nzc) || maxfwd <= 0)
1008 return 0;
1009 for(i = k = rv = 0; i < nzc; i++)
1010 if ((j = zci[i]) < nv0) {
1011 if (k >= maxfwd)
1012 goto done;
1013 vrefx[k++] = j;
1014 }
1015 else {
1016 if (!(vr = cexps[j-nv0].vref))
1017 goto done;
1018 vre = vr + *vr;
1019 while(++vr <= vre)
1020 if (!zc[*vr]++)
1021 zci[nzc++] = *vr;
1022 }
1023 if (k >= nvref) {
1024 nvref = (maxfwd + 1)*(ncom_togo < vrefGulp
1025 ? ncom_togo : vrefGulp);
1026 vrefnext = (int *)M1alloc(nvref*Sizeof(int));
1027 }
1028 vr = ce->vref = vrefnext;
1029 *vr = k;
1030 vrefnext += i = k + 1;
1031 nvref -= i;
1032 for(i = 0; i < k; i++)
1033 *++vr = vrefx[i];
1034 if (nderp > 3*k && !nocopy) {
1035 *ip = k;
1036 return 2;
1037 }
1038 else {
1039 done:
1040 if (nocopy || nderp > 3*nzc0)
1041 rv = 1;
1042 }
1043 while(nzc > nzc0)
1044 zc[zci[--nzc]] = 0;
1045 return rv;
1046 }
1047
1048 static void
1049 #ifdef KR_headers
cexp_read(R,k,nlin)1050 cexp_read(R, k, nlin) EdRead *R; int k, nlin;
1051 #else
1052 cexp_read(EdRead *R, int k, int nlin)
1053 #endif
1054 {
1055 int a, fk, i, j, la0, na;
1056 int *z1, **zp;
1057 funnel *f, **fp;
1058 linpart *L, *Le;
1059 expr *e;
1060 cplist *cl, *cl0;
1061 cexp *ce;
1062 Static *S = (Static *)R->S;
1063 ASL_fg *asl = S->asl;
1064
1065 ce = cexps + k - nv0;
1066 L = ce->L = linpt_read(R, ce->nlin = nlin);
1067 nocopy = 0;
1068 last_d = 0;
1069 ce->z.i = la0 = lasta;
1070 nderps += nderp;
1071 nderp = 0;
1072 e = ce->e = eread(R, want_derivs);
1073 if (la0 == lasta) {
1074 a = lasta++;
1075 if (e->op != f_OPNUM)
1076 new_derp(S, e->a, a, &edagread_one);
1077 }
1078 else
1079 a = e->a;
1080 var_e[k].a = a;
1081 ce->zlen = lasta - la0;
1082 for(Le = L + nlin; L < Le; L++) {
1083 new_derp(S, i = L->v.i, a, &L->fac);
1084 if (!zc[i]++)
1085 zci[nzc++] = i;
1086 }
1087 if (zp = zaC)
1088 *(zp += k - nv0) = 0;
1089 if (fk = funnelkind(S, ce, &i)) {
1090 /* arrange to funnel */
1091 fp = k < nv0b ? &f_b : k < nv0c ? &f_c : &f_o;
1092 ce->funneled = f = (funnel *)mem(sizeof(funnel));
1093 f->next = *fp;
1094 *fp = f;
1095 f->ce = ce;
1096 if (imap_len < lasta)
1097 imap_alloc(S);
1098 if (fk == 1) {
1099 f->fulld = last_d;
1100 a = lasta00;
1101 z1 = 0;
1102 if (zp) {
1103 for(i = j = 0; i < nzc; i++)
1104 if (zci[i] >= nv0)
1105 j++;
1106 if (j) {
1107 i = (2*j + 1)*sizeof(int);
1108 *zp = z1 = j > 20
1109 ? (int *)M1alloc(i)
1110 : (int *)mem(i);
1111 *z1++ = j;
1112 }
1113 }
1114 for(i = nzc; --i >= 0; )
1115 if ((j = zci[i]) >= nv0) {
1116 if (z1) {
1117 *z1++ = j;
1118 *z1++ = a;
1119 }
1120 imap[var_e[j].a] = a++;
1121 }
1122 if ((na = ce->zlen) || a > lasta00)
1123 na += a - nv1;
1124 f->fcde.zaplen = na*sizeof(real);
1125 i = nzc;
1126 derpadjust(S, last_d, a, 0);
1127 }
1128 else {
1129 f->fulld = 0;
1130 f->fcde.e = e;
1131 comsubs(S, ce->zlen, &f->fcde, zp);
1132 memcpy(zci, vrefx, i*sizeof(int));
1133 }
1134 last_d = 0;
1135 cl0 = 0;
1136 while(i > 0)
1137 if ((a = var_e[zci[--i]].a) != nv1) {
1138 new_derp(S, a, lasta0, 0);
1139 cl = (cplist *)mem(sizeof(cplist));
1140 cl->next = cl0;
1141 cl0 = cl;
1142 cl->ca.i = imap[last_d->a.i];
1143 cl->cfa = last_d->c.rp =
1144 (real *)mem(sizeof(real));
1145 }
1146 f->cl = cl0;
1147 var_e[k].a = lasta0++;
1148 lasta = lasta0;
1149 }
1150 lasta0 = lasta;
1151 ce->d = last_d;
1152 ce->cref = crefs(S);
1153 --ncom_togo;
1154 }
1155
1156 static void
1157 #ifdef KR_headers
cexp1_read(R,j,k,nlin)1158 cexp1_read(R, j, k, nlin) EdRead *R;
1159 #else
1160 cexp1_read(EdRead *R, int j, int k, int nlin)
1161 #endif
1162 {
1163 Static *S = (Static *)R->S;
1164 ASL_fg *asl = S->asl;
1165 linpart *L, *Le;
1166 cexp1 *ce = cexps1 + (k - nv01);
1167 expr *e;
1168 int la0;
1169
1170 L = ce->L = linpt_read(R, ce->nlin = nlin);
1171
1172 if (!lastj) {
1173 last_d = 0;
1174 if (amax1 < lasta)
1175 amax1 = lasta;
1176 lasta = lasta0;
1177 lastj = j;
1178 }
1179 la0 = lasta;
1180 e = ce->e = eread(R, want_derivs);
1181 if (la0 == lasta) {
1182 j = lasta++;
1183 if (e->op != f_OPNUM)
1184 new_derp(S, e->a, j, &edagread_one);
1185 }
1186 else
1187 j = e->a;
1188 var_e[k].a = j;
1189 for(Le = L + nlin; L < Le; L++)
1190 new_derp(S, L->v.i, j, &L->fac);
1191 last_cex = k;
1192 }
1193
1194 static void
1195 #ifdef KR_headers
ifadjust(a,e)1196 ifadjust(a, e) real *a; expr_if *e;
1197 #else
1198 ifadjust(real *a, expr_if *e)
1199 #endif
1200 {
1201 for(; e; e = e->next) {
1202 e->Tv.rp = &a[e->Tv.i];
1203 e->Fv.rp = &a[e->Fv.i];
1204 }
1205 }
1206
1207 static void
1208 #ifdef KR_headers
vargadjust(a,e)1209 vargadjust(a, e) real *a; expr_va *e;
1210 #else
1211 vargadjust(real *a, expr_va *e)
1212 #endif
1213 {
1214 de *d;
1215
1216 for(; e; e = e->next) {
1217 for(d = e->L.d; d->e; d++)
1218 d->dv.rp = &a[d->dv.i];
1219 }
1220 }
1221
1222 static void
1223 #ifdef KR_headers
funneladj1(a,f)1224 funneladj1(a, f) real *a; funnel *f;
1225 #else
1226 funneladj1(real *a, funnel *f)
1227 #endif
1228 {
1229 derp *d;
1230 cplist *cl;
1231
1232 for(; f; f = f->next) {
1233 if (d = f->fulld) {
1234 f->fcde.d = d;
1235 do {
1236 d->a.rp = &a[d->a.i];
1237 d->b.rp = &a[d->b.i];
1238 }
1239 while(d = d->next);
1240 }
1241 for(cl = f->cl; cl; cl = cl->next)
1242 cl->ca.rp = &a[cl->ca.i];
1243 }
1244 }
1245
1246 static void
1247 #ifdef KR_headers
funneladjust(asl)1248 funneladjust(asl) ASL_fg *asl;
1249 #else
1250 funneladjust(ASL_fg *asl)
1251 #endif
1252 {
1253 cexp *c, *ce;
1254 linpart *L, *Le;
1255 real *a = adjoints;
1256 c = cexps;
1257 for(ce = c + ncom0; c < ce; c++)
1258 if (L = c->L)
1259 for(Le = L + c->nlin; L < Le; L++)
1260 L->v.rp = &var_e[L->v.i].v;
1261
1262 funneladj1(a, f_b);
1263 funneladj1(a, f_c);
1264 funneladj1(a, f_o);
1265 }
1266
1267 static void
1268 #ifdef KR_headers
com1adjust(asl)1269 com1adjust(asl) ASL_fg *asl;
1270 #else
1271 com1adjust(ASL_fg *asl)
1272 #endif
1273 {
1274 cexp1 *c, *ce;
1275 linpart *L, *Le;
1276 expr_v *v = var_e;
1277
1278 for(c = cexps1, ce = c + ncom1; c < ce; c++)
1279 for(L = c->L, Le = L + c->nlin; L < Le; L++)
1280 L->v.rp = &v[L->v.i].v;
1281 }
1282 #endif /* Just_Linear */
1283
1284 static void
1285 #ifdef KR_headers
goff_comp(asl)1286 goff_comp(asl) ASL_fg *asl;
1287 #else
1288 goff_comp(ASL_fg *asl)
1289 #endif
1290 {
1291 int *ka = A_colstarts + 1;
1292 cgrad **cgx, **cgxe;
1293 cgrad *cg;
1294
1295 cgx = Cgrad;
1296 cgxe = cgx + asl->i.n_con0;
1297 while(cgx < cgxe)
1298 for(cg = *cgx++; cg; cg = cg->next)
1299 cg->goff = ka[cg->varno]++;
1300 }
1301
1302 static void
1303 #ifdef KR_headers
colstart_inc(S)1304 colstart_inc(S) Static *S;
1305 #else
1306 colstart_inc(Static *S)
1307 #endif
1308 {
1309 int *ka, *kae;
1310 ASL_fg *asl = S->asl;
1311
1312 ka = A_colstarts;
1313 kae = ka + asl->i.n_var0;
1314 while(ka <= kae)
1315 ++*ka++;
1316 }
1317
1318 #ifndef Just_Linear
1319
1320 static void
1321 #ifdef KR_headers
zerograd_chk(S)1322 zerograd_chk(S) Static *S;
1323 #else
1324 zerograd_chk(Static *S)
1325 #endif
1326 {
1327 ASL_fg *asl = S->asl;
1328 int j, n, nv, *z, **zg;
1329 ograd *og, **ogp, **ogpe;
1330
1331 if (!(nv = asl->i.nlvog))
1332 nv = nv0;
1333 zerograds = 0;
1334 ogp = Ograd;
1335 ogpe = ogp + (j = n_obj);
1336 while(ogp < ogpe) {
1337 og = *ogp++;
1338 n = 0;
1339 while(og) {
1340 j += og->varno - n;
1341 n = og->varno + 1;
1342 if (n >= nv)
1343 break;
1344 og = og->next;
1345 }
1346 if (n < nv)
1347 j += nv - n;
1348 }
1349 if (j == n_obj)
1350 return;
1351 zerograds = zg = (int **)mem(n_obj*sizeof(int*)+j*sizeof(int));
1352 z = (int*)(zg + n_obj);
1353 ogp = Ograd;
1354 while(ogp < ogpe) {
1355 *zg++ = z;
1356 og = *ogp++;
1357 n = 0;
1358 while(og) {
1359 while(n < og->varno)
1360 *z++ = n++;
1361 og = og->next;
1362 if (++n >= nv)
1363 break;
1364 }
1365 while(n < nv)
1366 *z++ = n++;
1367 *z++ = -1;
1368 }
1369 }
1370 #endif /* Just_Linear */
1371
1372 static void
1373 #ifdef KR_headers
adjust_compl_rhs(asl,opnum)1374 adjust_compl_rhs(asl, opnum) ASL_fg *asl; efunc *opnum;
1375 #else
1376 adjust_compl_rhs(ASL_fg *asl, efunc *opnum)
1377 #endif
1378 {
1379 cde *C;
1380 expr *e;
1381 int *Cvar, i, j, nc, stride;
1382 real *L, *U, t;
1383
1384 L = LUrhs;
1385 if (U = Urhsx)
1386 stride = 1;
1387 else {
1388 U = L + 1;
1389 stride = 2;
1390 }
1391 C = con_de;
1392 Cvar = cvar;
1393 nc = n_con;
1394 for(i = nlc; i < nc; i++)
1395 if (Cvar[i] && (e = C[i].e) && e->op == opnum
1396 && (t = ((expr_n*)e)->v) != 0.) {
1397 ((expr_n*)e)->v = 0.;
1398 if (L[j = stride*i] > negInfinity)
1399 L[j] -= t;
1400 if (U[j] < Infinity)
1401 U[j] -= t;
1402 }
1403 }
1404
1405 static void
1406 #ifdef KR_headers
adjust(S,flags)1407 adjust(S, flags) Static *S; int flags;
1408 #else
1409 adjust(Static *S, int flags)
1410 #endif
1411 {
1412 ASL_fg *asl = S->asl;
1413 #ifndef Just_Linear
1414 derp *d, **dp;
1415 real *a = adjoints;
1416 relo *r;
1417
1418 for(r = relolist; r; r = r->next) {
1419 dp = &r->D;
1420 while(d = *dp) {
1421 d->a.rp = &a[d->a.i];
1422 d->b.rp = &a[d->b.i];
1423 dp = &d->next;
1424 }
1425 *dp = r->Dnext;
1426 }
1427 ifadjust(adjoints, iflist);
1428 vargadjust(adjoints, varglist);
1429 if (ncom0)
1430 funneladjust(asl);
1431 com1adjust(asl);
1432 if (n_obj)
1433 zerograd_chk(S);
1434 #endif /* Just_Linear */
1435 if (k_seen) {
1436 if (!A_vals)
1437 goff_comp(asl);
1438 else if (Fortran)
1439 colstart_inc(S);
1440 }
1441 if (n_cc > nlcc && nlc < n_con
1442 && !(flags & ASL_no_linear_cc_rhs_adjust))
1443 adjust_compl_rhs(asl, f_OPNUM);
1444 }
1445
1446 static void
1447 #ifdef KR_headers
br_read(R,nc,nc1,Lp,U,Cvar,nv)1448 br_read(R, nc, nc1, Lp, U, Cvar, nv)
1449 EdRead *R; real **Lp, *U; int nc, nc1, *Cvar, nv;
1450 #else
1451 br_read(EdRead *R, int nc, int nc1, real **Lp, real *U, int *Cvar, int nv)
1452 #endif
1453 {
1454 int i, inc, j, k;
1455 real *L;
1456 ASL *asl = R->asl;
1457
1458 if (!(L = *Lp)) {
1459 if (nc1 < nc)
1460 nc1 = nc;
1461 L = *Lp = (real *)M1alloc(2*sizeof(real)*nc1);
1462 }
1463 if (U)
1464 inc = 1;
1465 else {
1466 U = L + 1;
1467 inc = 2;
1468 }
1469 xscanf(R, ""); /* purge line */
1470 for(i = 0; i < nc; i++, L += inc, U += inc) {
1471 switch(edag_peek(R) - '0') {
1472 case 0:
1473 if (xscanf(R,"%lf %lf",L,U)!= 2)
1474 badline(R);
1475 break;
1476 case 1:
1477 if (xscanf(R, "%lf", U) != 1)
1478 badline(R);
1479 *L = negInfinity;
1480 break;
1481 case 2:
1482 if (xscanf(R, "%lf", L) != 1)
1483 badline(R);
1484 *U = Infinity;
1485 break;
1486 case 3:
1487 *L = negInfinity;
1488 *U = Infinity;
1489 xscanf(R, ""); /* purge line */
1490 break;
1491 case 4:
1492 if (xscanf(R, "%lf", L) != 1)
1493 badline(R);
1494 *U = *L;
1495 break;
1496 case 5:
1497 if (Cvar) {
1498 if (xscanf(R, "%d %d", &k, &j) == 2
1499 && j > 0 && j <= nv) {
1500 Cvar[i] = j;
1501 *L = k & 2 ? negInfinity : 0.;
1502 *U = k & 1 ? Infinity : 0.;
1503 break;
1504 }
1505 }
1506 default:
1507 badline(R);
1508 }
1509 }
1510 }
1511
1512 #ifndef Just_Linear
1513
1514 static expr *
1515 #ifdef KR_headers
aholread(R)1516 aholread(R) EdRead *R;
1517 #else
1518 aholread(EdRead *R)
1519 #endif
1520 {
1521 int i, k;
1522 expr_h *rvh;
1523 char *s1;
1524 FILE *nl = R->nl;
1525 Static *S = (Static *)R->S;
1526
1527 k = getc(nl);
1528 if (k < '1' || k > '9')
1529 badline(R);
1530 i = k - '0';
1531 while((k = getc(nl)) != ':') {
1532 if (k < '0' || k > '9')
1533 badline(R);
1534 i = 10*i + k - '0';
1535 }
1536 rvh = (expr_h *)mem_ASL(R->asl, memadj(sizeof(expr_h) + i));
1537 for(s1 = rvh->sym;;) {
1538 if ((k = getc(nl)) < 0) {
1539 fprintf(Stderr,
1540 "Premature end of file in aholread, line %ld of %s\n",
1541 R->Line, R->asl->i.filename_);
1542 exit_ASL(R,1);
1543 }
1544 if (k == '\n') {
1545 R->Line++;
1546 if (!i)
1547 break;
1548 }
1549 if (--i < 0)
1550 badline(R);
1551 *s1++ = k;
1552 }
1553 *s1 = 0;
1554 rvh->op = f_OPHOL;
1555 rvh->a = nv1;
1556 return (expr *)rvh;
1557 }
1558
1559 static expr *
1560 #ifdef KR_headers
bholread(R)1561 bholread(R) EdRead *R;
1562 #else
1563 bholread(EdRead *R)
1564 #endif
1565 {
1566 int i;
1567 expr_h *rvh;
1568 char *s;
1569 Static *S = (Static *)R->S;
1570 ASL_fg *asl = S->asl;
1571
1572 if (xscanf(R, "%d", &i) != 1)
1573 badline(R);
1574 rvh = (expr_h *)mem(memadj(sizeof(expr_h) + i));
1575 s = rvh->sym;
1576 if (fread(s, i, 1, R->nl) != 1)
1577 badline(R);
1578 s[i] = 0;
1579 rvh->op = f_OPHOL;
1580 rvh->a = nv1;
1581 for(;;) switch(*s++) {
1582 case 0: goto break2; /* so we return at end of fcn */
1583 case '\n': R->Line++;
1584 }
1585 break2:
1586 return (expr *)rvh;
1587 }
1588
1589 static void
1590 #ifdef KR_headers
nlvzap(S,i,j)1591 nlvzap(S, i, j) Static *S; int i, j;
1592 #else
1593 nlvzap(Static *S, int i, int j)
1594 #endif
1595 {
1596 int n = nv1;
1597 expr_v *v = S->var_e;
1598
1599 i -= j;
1600 while(--j >= 0)
1601 v[i+j].a = n;
1602 }
1603
1604 static real
1605 #ifdef KR_headers
Missing_func(al)1606 Missing_func(al) arglist *al;
1607 #else
1608 Missing_func(arglist *al)
1609 #endif
1610 {
1611 AmplExports *ae = al->AE;
1612 func_info *fi = (func_info*)al->funcinfo;
1613
1614 char *s = (char*)(*ae->Tempmem)(al->TMI, strlen(fi->name) + 64);
1615 (*ae->SprintF)(al->Errmsg = s,
1616 "Attempt to call unavailable function %s.",
1617 fi->name);
1618 return 0.;
1619 }
1620 #endif /* Just_Linear */
1621
1622 int
1623 #ifdef KR_headers
fg_read_ASL(a,nl,flags)1624 fg_read_ASL(a, nl, flags) ASL *a; FILE *nl; int flags;
1625 #else
1626 fg_read_ASL(ASL *a, FILE *nl, int flags)
1627 #endif
1628 {
1629 cgrad *cg, **cgp;
1630 expr_v *e;
1631 int *ka;
1632 int i, i1, j, k, nc, nco, nlcon, no, nv, nvc, nvo, nz, readall;
1633 ograd *og, **ogp;
1634 real t;
1635 unsigned x;
1636 #ifdef Just_Linear
1637 #define ASL_readtype ASL_read_f
1638 #else /* Just_Linear */
1639 #define ASL_readtype ASL_read_fg
1640 int maxfwd1, ncom, nlin;
1641 func_info *fi;
1642 char fname[128];
1643 #endif /* Just_Linear */
1644 EdRead ER, *R;
1645 Static SS, *S;
1646 Jmp_buf JB;
1647
1648 ASL_CHECK(a, ASL_readtype, who);
1649 S = ed_reset(&SS, a);
1650 R = EdReadInit_ASL(&ER, a, nl, S);
1651 if (flags & ASL_return_read_err) {
1652 a->i.err_jmp_ = &JB;
1653 i = setjmp(JB.jb);
1654 if (i) {
1655 a->i.err_jmp_ = 0;
1656 return i;
1657 }
1658 }
1659 nlcon = a->i.n_lcon_;
1660 if (nlcon && !(flags & ASL_allow_CLP)) {
1661 if (a->i.err_jmp_)
1662 return ASL_readerr_CLP;
1663 sorry_CLP(R, "logical constraints");
1664 }
1665 if ((readall = flags & ASL_keep_all_suffixes)
1666 && a->i.nsuffixes)
1667 readall |= 1;
1668 #define asl ((ASL_fg*)a)
1669 ER.lineinc = 1;
1670 if (!size_expr_n) size_expr_n = sizeof(expr_n);
1671 #ifdef Just_Linear
1672 /* Tried "xscanf = binary ? bscanf : ascanf;", but this
1673 * encounters a bug in SGI's C compiler.
1674 */
1675 if (binary_nl)
1676 xscanf = bscanf;
1677 else
1678 xscanf = ascanf;
1679 #else
1680 if (!(SS._r_ops = asl->I.r_ops_))
1681 SS._r_ops = r_ops_ASL;
1682 if (c_cexp1st)
1683 *c_cexp1st = 0;
1684 if (o_cexp1st)
1685 *o_cexp1st = comc1;
1686 if (nfunc)
1687 func_add(a);
1688 if (binary_nl) {
1689 holread = bholread;
1690 xscanf = bscanf;
1691 }
1692 else {
1693 holread = aholread;
1694 xscanf = ascanf;
1695 }
1696
1697 ncom = comb + comc + como + comc1 + como1;
1698 #endif /* Just_Linear */
1699 nc = n_con;
1700 no = n_obj;
1701 nvc = c_vars;
1702 nvo = o_vars;
1703 if (no < 0 || (nco = nc + no + nlcon) <= 0)
1704 scream(R, ASL_readerr_corrupt,
1705 "edagread: nc = %d, no = %d, nlcon = %d\n",
1706 nc, no, nlcon);
1707 if (pi0) {
1708 memset(pi0, 0, nc*sizeof(real));
1709 if (havepi0)
1710 memset(havepi0, 0, nc);
1711 }
1712 #ifdef Just_Linear
1713 nv = nv0 = nvc > nvo ? nvc : nvo;
1714 x = nco*sizeof(cde) + no*sizeof(ograd *) + nv*sizeof(expr_v) + no;
1715 #else
1716 nv1 = nv0 = nvc > nvo ? nvc : nvo;
1717 max_var = nv = nv0 + ncom;
1718 combc = comb + comc;
1719 ncom0 = ncom_togo = combc + como;
1720 nzclim = ncom0 >> 3;
1721 ncom1 = comc1 + como1;
1722 nv0b = nv0 + comb;
1723 nv0c = nv0b + comc;
1724 nv01 = nv0 + ncom0;
1725 last_cex = nv011 = nv01 - 1;
1726 amax = lasta = lasta0 = lasta00 = nv1 + 1;
1727 if ((maxfwd1 = maxfwd + 1) > 1)
1728 nvref = maxfwd1*((ncom0 < vrefGulp ? ncom0 : vrefGulp) + 1);
1729 x = nco*sizeof(cde) + no*sizeof(ograd *)
1730 + nv*(sizeof(expr_v) + 2*sizeof(int))
1731 + ncom0*sizeof(cexp)
1732 + ncom1*sizeof(cexp1)
1733 + nfunc*sizeof(func_info *)
1734 + nvref*sizeof(int)
1735 + no;
1736 SS.nvar0 = a->i.n_var0;
1737 if (!(SS.nvinc = a->i.n_var_ - SS.nvar0))
1738 SS.nvar0 += ncom0 + ncom1;
1739 #endif /* Just_Linear */
1740 if (X0)
1741 memset(X0, 0, nv0*sizeof(real));
1742 if (havex0)
1743 memset(havex0, 0, nv0);
1744 e = var_e = (expr_v *)M1zapalloc(x);
1745 con_de = (cde *)(e + nv);
1746 lcon_de = con_de + nc;
1747 obj_de = lcon_de + nlcon;
1748 Ograd = (ograd **)(obj_de + no);
1749 if (A_vals) {
1750 if (!A_rownos)
1751 A_rownos = (int *)M1alloc(asl->i.nzc_*sizeof(int));
1752 }
1753 else if (nc)
1754 Cgrad = (cgrad **)M1zapalloc(nc*sizeof(cgrad *));
1755 #ifdef Just_Linear
1756 objtype = (char *)(Ograd + no);
1757 #else
1758 var_ex = e + nv0;
1759 var_ex1 = var_ex + ncom0;
1760 for(k = 0; k < nv; e++) {
1761 e->op = f_OPVARVAL;
1762 e->a = k++;
1763 }
1764 if (skip_int_derivs) {
1765 if (nlvbi)
1766 nlvzap(S, nlvb, nlvbi);
1767 if (nlvci)
1768 nlvzap(S, nlvb+nlvc, nlvci);
1769 if (nlvoi)
1770 nlvzap(S, nlvb+nlvc+nlvo, nlvoi);
1771 }
1772 cexps = (cexp *)(Ograd + no);
1773 cexps1 = (cexp1 *)(cexps + ncom0);
1774 funcs = (func_info **)(cexps1 + ncom1);
1775 zc = (int*)(funcs + nfunc);
1776 zci = zc + nv;
1777 vrefx = zci + nv;
1778 objtype = (char *)(vrefx + nvref);
1779 if (nvref) {
1780 vrefnext = vrefx + maxfwd1;
1781 nvref -= maxfwd1;
1782 }
1783 last_d = 0;
1784 #endif
1785 if (n_cc && !cvar)
1786 cvar = (int*)M1alloc(nc*sizeof(int));
1787 if (cvar)
1788 memset(cvar, 0, nc*sizeof(int));
1789 ka = 0;
1790 nz = 0;
1791 for(;;) {
1792 ER.can_end = 1;
1793 i = edag_peek(R);
1794 if (i == EOF) {
1795 #ifndef Just_Linear
1796 free(imap);
1797 /* Make amax long enough for nlc to handle */
1798 /* var_e[i].a for common variables i. */
1799 if (ncom0) {
1800 i = comb + como;
1801 if (i < combc)
1802 i = combc;
1803 if ((i += nv1 + 1) > amax)
1804 amax = i;
1805 }
1806 adjoints = (real *)M1zapalloc(amax*Sizeof(real));
1807 adjoints_nv1 = &adjoints[nv1];
1808 nderps += nderp;
1809 #endif /* Just_Linear */
1810 adjust(S, flags);
1811 nzjac = nz;
1812 if (!Lastx)
1813 Lastx = (real *)M1alloc(nv0*sizeof(real));
1814 fclose(nl);
1815 #ifndef Just_Linear
1816 a->p.Objval = obj1val_ASL;
1817 a->p.Objgrd = obj1grd_ASL;
1818 a->p.Conval = con1val_ASL;
1819 a->p.Jacval = jac1val_ASL;
1820 a->p.Conival = con1ival_ASL;
1821 a->p.Congrd = con1grd_ASL;
1822 a->p.Lconval = lcon1val_ASL;
1823 a->p.Xknown = x1known_ASL;
1824 #endif /* Just_Linear */
1825 a->i.err_jmp_ = 0;
1826 return 0;
1827 }
1828 ER.can_end = 0;
1829 k = -1;
1830 switch(i) {
1831 case 'C':
1832 xscanf(R, "%d", &k);
1833 if (k < 0 || k >= nc)
1834 badline(R);
1835 co_read(R,con_de,c_cexp1st,k,zac,want_derivs);
1836 break;
1837 #ifdef Just_Linear
1838 case 'F':
1839 case 'V':
1840 sorry_nonlin(R);
1841 #else
1842 case 'F':
1843 if (xscanf(R, "%d %d %d %127s",
1844 &i, &j, &k, fname) != 4
1845 || i < 0 || i >= nfunc)
1846 badline(R);
1847 if (fi = func_lookup(a, fname,0)) {
1848 if (fi->nargs != k && fi->nargs >= 0
1849 && (k >= 0 || fi->nargs < -(k+1)))
1850 scream(R, ASL_readerr_argerr,
1851 "function %s: disagreement of nargs: %d and %d\n",
1852 fname,fi->nargs, k);
1853 }
1854 else {
1855 fi = (func_info *)mem(sizeof(func_info));
1856 fi->ftype = j;
1857 fi->nargs = k;
1858 fi->funcp = 0;
1859 fi->name = (Const char *)strcpy((char*)
1860 mem(memadj(strlen(fname)+1)),
1861 fname);
1862 }
1863 if (!fi->funcp && !(fi->funcp = dynlink(fname))){
1864 if (!(flags & ASL_allow_missing_funcs))
1865 scream(R, ASL_readerr_unavail,
1866 "function %s not available\n",
1867 fname);
1868 fi->funcp = Missing_func;
1869 fi->funcinfo = (Char*)fi;
1870 }
1871 funcs[i] = fi;
1872 break;
1873 case 'L':
1874 xscanf(R, "%d", &k);
1875 if (k < 0 || k >= nlcon)
1876 badline(R);
1877 co_read(R, lcon_de, 0, k, 0, 0);
1878 break;
1879 case 'V':
1880 if (xscanf(R, "%d %d %d", &k, &nlin, &j) != 3)
1881 badline(R);
1882 if (k >= SS.nvar0)
1883 k += SS.nvinc;
1884 if (k < nv0 || k >= nv)
1885 badline(R);
1886 if (j)
1887 cexp1_read(R, j, k, nlin);
1888 else
1889 cexp_read(R, k, nlin);
1890 break;
1891 #endif /* Just_Linear */
1892 case 'G':
1893 if (xscanf(R, "%d %d", &j, &k) != 2
1894 || j < 0 || j >= no || k <= 0 || k > nvo)
1895 badline(R);
1896 ogp = Ograd + j;
1897 while(k--) {
1898 *ogp = og = (ograd *)mem(sizeof(ograd));
1899 ogp = &og->next;
1900 if (xscanf(R, "%d %lf", &og->varno,
1901 &og->coef) != 2)
1902 badline(R);
1903 }
1904 *ogp = 0;
1905 break;
1906 case 'J':
1907 if (xscanf(R, "%d %d", &j, &k) != 2
1908 || j < 0 || j >= nc || k <= 0 || k > nvc)
1909 badline(R);
1910 nz += k;
1911 if (ka) {
1912 if (!A_vals)
1913 goto cg_read;
1914 j += Fortran;
1915 while(k--) {
1916 if (xscanf(R, "%d %lf",
1917 &i, &t) != 2)
1918 badline(R);
1919 i1 = ka[i]++;
1920 A_vals[i1] = t;
1921 A_rownos[i1] = j;
1922 }
1923 break;
1924 }
1925 cg_read:
1926 cgp = Cgrad + j;
1927 j = 0;
1928 while(k--) {
1929 *cgp = cg = (cgrad *)mem(sizeof(cgrad));
1930 cgp = &cg->next;
1931 if (ka) {
1932 if (xscanf(R, "%d %lf",
1933 &cg->varno,
1934 &cg->coef) != 2)
1935 badline(R);
1936 }
1937 else
1938 if (xscanf(R, "%d %d %lf",
1939 &cg->varno, &j,
1940 &cg->coef) != 3)
1941 badline(R);
1942 cg->goff = j;
1943 }
1944 *cgp = 0;
1945 break;
1946 case 'O':
1947 if (xscanf(R, "%d %d", &k, &j) != 2
1948 || k < 0 || k >= no)
1949 badline(R);
1950 objtype[k] = j;
1951 co_read(R,obj_de,o_cexp1st,k,zao,want_derivs);
1952 break;
1953 case 'S':
1954 Suf_read_ASL(R, readall);
1955 break;
1956 case 'r':
1957 br_read(R, asl->i.n_con0, nc, &LUrhs,
1958 Urhsx, cvar, nv0);
1959 break;
1960 case 'b':
1961 br_read(R, asl->i.n_var0, nv0, &LUv,
1962 Uvx, 0, 0);
1963 break;
1964 case 'k':
1965 k_seen++;
1966 k = asl->i.n_var0;
1967 if (!xscanf(R,"%d",&j) || j != k - 1)
1968 badline(R);
1969 if (!(ka = A_colstarts)) {
1970 if ((i = k) < n_var)
1971 i = n_var;
1972 ka = A_colstarts = (int *)
1973 M1alloc((i+1)*Sizeof(int));
1974 }
1975 *ka++ = 0;
1976 *ka++ = 0; /* sic */
1977 while(--k > 0)
1978 if (!xscanf(R, "%d", ka++))
1979 badline(R);
1980 ka = A_colstarts + 1;
1981 break;
1982 case 'x':
1983 if (!xscanf(R,"%d",&k)
1984 || k < 0 || k > nv0)
1985 badline(R);
1986 if (!X0 && want_xpi0 & 1) {
1987 x = nv0*sizeof(real);
1988 if (want_xpi0 & 4)
1989 x += nv0;
1990 X0 = (real *)M1zapalloc(x);
1991 if (want_xpi0 & 4)
1992 havex0 = (char*)(X0 + nv0);
1993 }
1994 while(k--) {
1995 if (xscanf(R, "%d %lf", &j, &t) != 2
1996 || j < 0 || j >= nv)
1997 badline(R);
1998 if (X0) {
1999 X0[j] = t;
2000 if (havex0)
2001 havex0[j] = 1;
2002 }
2003 }
2004 break;
2005 case 'd':
2006 if (!xscanf(R,"%d",&k)
2007 || k < 0 || k > nc)
2008 badline(R);
2009 if (!pi0 && want_xpi0 & 2) {
2010 x = nc*sizeof(real);
2011 if (want_xpi0 & 4)
2012 x += nc;
2013 pi0 = (real *)M1zapalloc(x);
2014 if (want_xpi0 & 4)
2015 havepi0 = (char*)(pi0 + nc);
2016 }
2017 while(k--) {
2018 if (xscanf(R, "%d %lf", &j, &t) != 2
2019 || j < 0 || j >= nc)
2020 badline(R);
2021 if (pi0) {
2022 pi0[j] = t;
2023 if (havepi0)
2024 havepi0[j] = 1;
2025 }
2026 }
2027 break;
2028 default:
2029 badline(R);
2030 }
2031 }
2032 }
2033