1 /*******************************************************************
2 Copyright (C) 2017 AMPL Optimization, Inc.; written by David M. Gay.
3
4 Permission to use, copy, modify, and distribute this software and its
5 documentation for any purpose and without fee is hereby granted,
6 provided that the above copyright notice appear in all copies and that
7 both that the copyright notice and this permission notice and warranty
8 disclaimer appear in supporting documentation.
9
10 The author and AMPL Optimization, Inc. disclaim all warranties with
11 regard to this software, including all implied warranties of
12 merchantability and fitness. In no event shall the author be liable
13 for any special, indirect or consequential damages or any damages
14 whatsoever resulting from loss of use, data or profits, whether in an
15 action of contract, negligence or other tortious action, arising out
16 of or in connection with the use or performance of this software.
17 *******************************************************************/
18
19 #include "jacpdim.h"
20
21 #ifdef __cplusplus
22 extern "C" {
23 #endif
24
25 static void
hv_fwd(expr * e)26 hv_fwd(expr *e)
27 {
28 argpair *da, *dae;
29 expr *e1, **ep;
30 real dO;
31
32 for(; e; e = e->fwd) {
33 e->aO = e->adO = 0;
34 switch(e->a) {
35
36 case Hv_timesR:
37 case Hv_binaryR:
38 e->dO.r = e->R.e->dO.r * e->dR;
39 break;
40
41 case Hv_timesLR:
42 case Hv_binaryLR:
43 case Hv_divLR:
44 e->dO.r = e->L.e->dO.r*e->dL + e->R.e->dO.r*e->dR;
45 break;
46
47 case Hv_timesL:
48 case Hv_unary:
49 e->dO.r = e->L.e->dO.r * e->dL;
50 break;
51
52 case Hv_vararg:
53 case Hv_if:
54 if ((e1 = ((expr_va *)e)->valf)) {
55 hv_fwd(e1);
56 e->dO.r = ((expr_va *)e)->vale->dO.r;
57 }
58 else if ((e1 = ((expr_va *)e)->val) && e1->op != f_OPNUM)
59 e->dO.r = e1->dO.r;
60 else
61 e->dO.r = 0;
62 break;
63
64 case Hv_plterm:
65 e->dO.r = e->dL * e->R.e->dO.r;
66 break;
67
68 case Hv_sumlist:
69 ep = e->R.ep;
70 for(dO = 0; (e1 = *ep); ep++)
71 dO += e1->dO.r;
72 e->dO.r = dO;
73 break;
74
75 case Hv_func:
76 da = ((expr_f *)e)->da;
77 dae = ((expr_f *)e)->dae;
78 for(dO = 0; da < dae; da++)
79 dO += da->e->dO.r * *da->u.v;
80 e->dO.r = dO;
81 break;
82
83 case Hv_negate:
84 e->dO.r = -e->L.e->dO.r;
85 break;
86
87 case Hv_plusR:
88 e->dO.r = e->R.e->dO.r;
89 break;
90
91 case Hv_plusL:
92 e->dO.r = e->L.e->dO.r;
93 break;
94
95 case Hv_plusLR:
96 e->dO.r = e->L.e->dO.r + e->R.e->dO.r;
97 break;
98
99 case Hv_minusR:
100 e->dO.r = -e->R.e->dO.r;
101 break;
102
103 case Hv_minusLR:
104 e->dO.r = e->L.e->dO.r - e->R.e->dO.r;
105 break;
106
107 default:/*DEBUG*/
108 fprintf(Stderr, "bad e->a = %d in hv_fwd\n", e->a);
109 exit(1);
110 }
111 }
112 }
113
114 static void
func_back(expr_f * f)115 func_back(expr_f *f)
116 {
117 argpair *da, *da1, *dae;
118 expr *e;
119 real **fh;
120 real aO, adO, t;
121
122 fh = f->fh;
123 aO = f->aO;
124 adO = f->adO;
125 da = f->da;
126 for(dae = f->dae; da < dae; da++) {
127 e = da->e;
128 e->adO += (t = *da->u.v) * adO;
129 e->aO += t * aO;
130 t = adO*e->dO.r;
131 for(da1 = f->da; da1 < dae; da1++) {
132 e = da1->e;
133 e->aO += t * **fh++;
134 }
135 }
136 }
137
138 static void
funnel_back(ASL_pfgh * asl,cexp * c,expr_v * v,real t)139 funnel_back(ASL_pfgh *asl, cexp *c, expr_v *v, real t)
140 {
141 expr_v **vp, **vp1, **vpe;
142 hes_fun *hf;
143 ograd *og;
144 real *g, *h;
145 real aO, adO;
146
147 aO = v->aO = t;
148 adO = v->adO;
149 hf = c->hfun;
150 if ((og = hf->og)) {
151 do {
152 v = var_e + og->varno;
153 v->adO += (t = og->coef) * adO;
154 v->aO += t*aO;
155 }
156 while((og = og->next));
157 return;
158 }
159 g = hf->grdhes;
160 h = g + hf->n;
161 vp = hf->vp;
162 vpe = vp + hf->n;
163 do {
164 v = *vp++;
165 v->adO += (t = *g++) * adO;
166 v->aO += t*aO;
167 t = adO * v->dO.r;
168 vp1 = hf->vp;
169 do (*vp1++)->aO += t * *h++;
170 while(vp1 < vpe);
171 }
172 while(vp < vpe);
173 }
174
175 static void
hv_back(expr * e)176 hv_back(expr *e)
177 {
178 expr *e1, **ep, *e2;
179 real adO, t1, t2;
180
181 if (!e || (!e->aO && !e->adO))
182 return;
183 for(; e; e = e->bak)
184 switch(e->a) {
185 case Hv_binaryR:
186 e1 = e->R.e;
187 e1->adO += e->adO * e->dR;
188 e1->aO += e->aO * e->dR + e->adO * e1->dO.r * e->dR2;
189 break;
190
191 case Hv_binaryLR:
192 e1 = e->L.e;
193 e2 = e->R.e;
194 adO = e->adO;
195 t1 = adO * e1->dO.r;
196 t2 = adO * e2->dO.r;
197 e1->aO += e->aO*e->dL + t1*e->dL2 + t2*e->dLR;
198 e2->aO += e->aO*e->dR + t1*e->dLR + t2*e->dR2;
199 e1->adO += adO * e->dL;
200 e2->adO += adO * e->dR;
201 break;
202
203 case Hv_divLR:
204 e1 = e->L.e;
205 e2 = e->R.e;
206 adO = e->adO;
207 t1 = adO * e1->dO.r;
208 t2 = adO * e2->dO.r;
209 e1->aO += e->aO*e->dL + t2*e->dLR;
210 e2->aO += e->aO*e->dR + t1*e->dLR + t2*e->dR2;
211 e1->adO += adO * e->dL;
212 e2->adO += adO * e->dR;
213 break;
214
215 case Hv_unary:
216 e1 = e->L.e;
217 e1->adO += e->adO * e->dL;
218 e1->aO += e->aO * e->dL + e->adO * e1->dO.r * e->dL2;
219 break;
220
221 case Hv_vararg:
222 case Hv_if:
223 if ((e1 = ((expr_va *)e)->vale)) {
224 e1->aO = e->aO;
225 e1->adO = e->adO;
226 hv_back(e1);
227 }
228 else {
229 e1 = ((expr_va *)e)->val;
230 if (e1->op != f_OPNUM) {
231 e1->aO += e->aO;
232 e1->adO += e->adO;
233 }
234 }
235 break;
236
237 case Hv_plterm:
238 e->R.e->aO += e->dL * e->aO;
239 break;
240
241 case Hv_sumlist:
242 ep = e->R.ep;
243 t1 = e->aO;
244 t2 = e->adO;
245 while((e1 = *ep++)) {
246 e1->aO += t1;
247 e1->adO += t2;
248 }
249 break;
250
251 case Hv_func:
252 func_back((expr_f *)e);
253 break;
254
255 case Hv_negate:
256 e1 = e->L.e;
257 neg_end:
258 e1->aO -= e->aO;
259 e1->adO -= e->adO;
260 break;
261
262 case Hv_plusR:
263 e1 = e->R.e;
264 goto plus_end;
265
266 case Hv_plusL:
267 e1 = e->L.e;
268 plus_end:
269 e1->aO += e->aO;
270 e1->adO += e->adO;
271 break;
272
273 case Hv_plusLR:
274 e1 = e->L.e;
275 e1->aO += t1 = e->aO;
276 e1->adO += t2 = e->adO;
277 e1 = e->R.e;
278 e1->aO += t1;
279 e1->adO += t2;
280 break;
281
282 case Hv_minusR:
283 e1 = e->R.e;
284 goto neg_end;
285
286 case Hv_minusLR:
287 e1 = e->L.e;
288 e1->aO += t1 = e->aO;
289 e1->adO += t2 = e->adO;
290 e1 = e->R.e;
291 e1->aO -= t1;
292 e1->adO -= t2;
293 break;
294
295 case Hv_timesR:
296 e1 = e->R.e;
297 e1->adO += e->adO * (t1 = e->dR);
298 e1->aO += e->aO * t1;
299 break;
300
301 case Hv_timesL:
302 e1 = e->L.e;
303 e1->adO += e->adO * (t1 = e->dL);
304 e1->aO += e->aO * t1;
305 break;
306
307 case Hv_timesLR:
308 e1 = e->L.e;
309 e2 = e->R.e;
310 adO = e->adO;
311 e1->aO += e->aO*e->dL + adO * e2->dO.r;
312 e2->aO += e->aO*e->dR + adO * e1->dO.r;
313 e1->adO += adO * e->dL;
314 e2->adO += adO * e->dR;
315 break;
316
317 default:/*DEBUG*/
318 fprintf(Stderr, "bad e->a = %d in hv_back\n", e->a);
319 exit(1);
320 }
321 }
322
323 static void
hv_fwd0(ASL_pfgh * asl,cexp * c,expr_v * v)324 hv_fwd0(ASL_pfgh *asl, cexp *c, expr_v *v)
325 {
326 expr_v **vp, **vpe;
327 hes_fun *hf;
328 int i;
329 linarg *la;
330 linpart *L, *Le;
331 ograd *og;
332 real *g, x;
333
334 v->aO = v->adO = 0;
335 if ((hf = c->hfun)) {
336 x = 0;
337 if ((og = hf->og))
338 do x += og->coef * var_e[og->varno].dO.r;
339 while((og = og->next));
340 else {
341 g = hf->grdhes;
342 vp = hf->vp;
343 vpe = vp + hf->n;
344 do x += *g++ * (*vp++)->dO.r;
345 while(vp < vpe);
346 }
347 }
348 else if (c->ef) {
349 hv_fwd(c->ef);
350 x = c->ee->dO.r;
351 }
352 else if (c->e->op != f_OPNUM)
353 x = c->e->dO.r;
354 else
355 x = 0;
356 if ((la = c->la)) {
357 i = c - asl->I.cexps2_;
358 x += asl->P.dv[i].scale * la->v->dO.r;
359 }
360
361 else if ((L = c->L)) {
362 for(Le = L + c->nlin; L < Le; L++)
363 x += L->fac * ((expr_v*)L->v.vp)->dO.r;
364 }
365 v->dO.r = x;
366 }
367
368 static void
hfg_fwd(expr * e)369 hfg_fwd(expr *e)
370 {
371 expr *e1;
372
373 for(; e; e = e->fwd) {
374 e->aO = 0;
375 switch(e->a) {
376 case Hv_vararg:
377 case Hv_if:
378 if ((e1 = ((expr_va *)e)->valf))
379 hfg_fwd(e1);
380 }
381 }
382 }
383
384 static void
hfg_back(expr * e)385 hfg_back(expr *e)
386 {
387 expr *e1, **ep;
388 real aO;
389
390 if (!e || !e->aO)
391 return;
392 for(; e; e = e->bak)
393 switch(e->a) {
394 case Hv_timesR:
395 case Hv_binaryR:
396 e->R.e->aO += e->aO * e->dR;
397 break;
398
399 case Hv_binaryLR:
400 case Hv_timesLR:
401 case Hv_divLR:
402 aO = e->aO;
403 e->L.e->aO += aO * e->dL;
404 e->R.e->aO += aO * e->dR;
405 break;
406
407 case Hv_unary:
408 case Hv_timesL:
409 e->L.e->aO += e->aO * e->dL;
410 break;
411
412 case Hv_vararg:
413 case Hv_if:
414 if ((e1 = ((expr_va *)e)->vale)) {
415 e1->aO = e->aO;
416 hfg_back(e1);
417 }
418 else {
419 e1 = ((expr_va *)e)->val;
420 if (e1->op != f_OPNUM)
421 e1->aO += e->aO;
422 }
423 break;
424
425 case Hv_plterm:
426 e->R.e->aO += e->dL * e->aO;
427 break;
428
429 case Hv_sumlist:
430 ep = e->R.ep;
431 aO = e->aO;
432 while((e1 = *ep++))
433 e1->aO += aO;
434 break;
435
436 case Hv_func:
437 func_back((expr_f *)e);
438 break;
439
440 case Hv_negate:
441 e->L.e->aO -= e->aO;
442 break;
443
444 case Hv_plusR:
445 e->R.e->aO += e->aO;
446 break;
447
448 case Hv_plusL:
449 e->L.e->aO += e->aO;
450 break;
451
452 case Hv_plusLR:
453 e->L.e->aO += aO = e->aO;
454 e->R.e->aO += aO;
455 break;
456
457 case Hv_minusR:
458 e->R.e->aO -= e->aO;
459 break;
460
461 case Hv_minusLR:
462 e->L.e->aO += aO = e->aO;
463 e->R.e->aO -= aO;
464 break;
465
466 default:/*DEBUG*/
467 fprintf(Stderr, "bad e->a = %d in hfg_back\n", e->a);
468 exit(1);
469 }
470 }
471
472 static void
funnelhes(ASL_pfgh * asl)473 funnelhes(ASL_pfgh *asl)
474 {
475 cexp *c;
476 expr *e;
477 expr_v *v, **vp, **vp1, **vpe;
478 hes_fun *hf;
479 int n;
480 real *g, *h;
481
482 x0kind &= ~ASL_need_funnel;
483 for(hf = asl->I.hesthread; hf; hf = hf->hfthread) {
484 if (hf->og)
485 continue;
486 n = hf->n;
487 g = hf->grdhes;
488 h = g + n;
489 c = hf->c;
490 vp = hf->vp;
491 vpe = vp + n;
492
493 do (*vp++)->aO = 0;
494 while(vp < vpe);
495
496 hfg_fwd(c->ef);
497 e = c->ee;
498 e->aO = 1;
499 hfg_back(e);
500
501 vp = hf->vp;
502 do {
503 v = *vp++;
504 *g++ = v->aO;
505 v->dO.r = v->aO = v->adO = 0;
506 }
507 while(vp < vpe);
508
509 vp = hf->vp;
510 do {
511 v = *vp++;
512 v->dO.r = 1;
513 if ((e = c->ef))
514 hv_fwd(e);
515 if ((e = c->ee)) {
516 e->aO = 0;
517 e->adO = 1;
518 hv_back(e);
519 }
520 else if ((e = c->e)->op != f_OPNUM) {
521 e->aO = 0;
522 e->adO = 1;
523 }
524 v->dO.r = 0;
525 vp1 = hf->vp;
526 do {
527 v = *vp1++;
528 *h++ = v->aO;
529 v->aO = v->adO = 0;
530 }
531 while(vp1 < vpe);
532 }
533 while(vp < vpe);
534 }
535 }
536
537 static void
hvp0comp_ASL(ASL_pfgh * asl,real * hv,real * p,int nobj,real * ow,real * y)538 hvp0comp_ASL(ASL_pfgh *asl, real *hv, real *p, int nobj, real *ow, real *y)
539 /* p = direction */
540 /* y = Lagrange multipliers */
541 /* hv = result */
542 {
543 cexp *c, *c1, *ce;
544 expr *e;
545 expr_v *v, *x, *x0, *xe;
546 int *dvsp0, i, i0, i1, n, nc, no, noe;
547 linarg *la;
548 linpart *L, *Le;
549 ograd *og;
550 ps_func *f, *f0;
551 psb_elem *b, *be;
552 psg_elem *g, *ge;
553 real *cscale, t, t2, *p1;
554
555 #ifdef IGNORE_BOGUS_WARNINGS
556 c1 = ce = 0;
557 dvsp0 = 0;
558 i1 = 0;
559 #endif
560
561 if (x0kind & ASL_need_funnel)
562 funnelhes(asl);
563 if (nobj >= 0 && nobj < n_obj) {
564 ow = ow ? ow + nobj : &edag_one_ASL;
565 no = nobj;
566 noe = no + 1;
567 }
568 else {
569 no = noe = 0;
570 if (ow)
571 noe = n_obj;
572 }
573 n = c_vars >= o_vars ? c_vars : o_vars;
574 for(la = asl->P.lalist; la; la = la->lnext) {
575 og = la->nz;
576 t = p[og->varno]*og->coef;
577 while((og = og->next))
578 t += p[og->varno]*og->coef;
579 x = la->v;
580 x->dO.r = t;
581 x->aO = x->adO = 0;
582 }
583 p1 = p;
584 x = x0 = var_e;
585 for(xe = x + n; x < xe; x++) {
586 x->dO.r = *p1++;
587 x->aO = x->adO = 0.;
588 }
589 if (asl->P.ncom) {
590 x = var_ex;
591 dvsp0 = asl->P.dvsp0;
592 c = cexps;
593 i0 = *dvsp0;
594 for(ce = c1 = c + asl->P.ncom; c < ce; c++) {
595 for(i1 = *++dvsp0; i0 < i1; i0++)
596 hv_fwd0(asl, c1++, asl->P.vp[i0]);
597 hv_fwd0(asl, c, x++);
598 }
599 }
600 if (!y || (nc = n_con) <= 0)
601 goto no_y;
602 cscale = asl->i.lscale;
603 f0 = asl->P.cps;
604 for(i = 0; i < nc; ++i) {
605 if ((t2 = y[i])) {
606 if (cscale)
607 t2 *= cscale[i];
608 f = f0 + i;
609 for(b = f->b, be = b + f->nb; b < be; b++) {
610 if ((e = b->D.ef)) {
611 hv_fwd(e);
612 e = b->D.ee;
613 e->aO = 0;
614 e->adO = t2;
615 hv_back(e);
616 }
617 else if ((e = b->D.e)->op != f_OPNUM) {
618 e->aO = 0;
619 e->adO = t2;
620 }
621 }
622 for(g = f->g, ge = g + f->ng; g < ge; g++) {
623 for(b = g->E, be = b + g->ns; b < be; b++) {
624 if ((e = b->D.ef)) {
625 hv_fwd(e);
626 e = b->D.ee;
627 e->aO = 0;
628 e->adO = t2*g->g1;
629 hv_back(e);
630 }
631 else if ((e = b->D.e)->op != f_OPNUM) {
632 e->aO = 0;
633 e->adO = t2*g->g1;
634 }
635 }
636 if (g->g2) {
637 t = 0.;
638 for(og = g->og; og; og = og->next)
639 t += og->coef * p[og->varno];
640 t *= t2*g->g2;
641 for(og = g->og; og; og = og->next)
642 x0[og->varno].aO += t*og->coef;
643 }
644 }
645 }
646 }
647 no_y:
648 for(; no < noe; no++) {
649 if ((t2 = *ow++)) {
650 f = asl->P.ops + no;
651 for(b = f->b, be = b + f->nb; b < be; b++) {
652 if ((e = b->D.ef)) {
653 hv_fwd(e);
654 e = b->D.ee;
655 e->aO = 0;
656 e->adO = t2;
657 hv_back(e);
658 }
659 else if ((e = b->D.e)->op != f_OPNUM) {
660 e->aO = 0;
661 e->adO = t2;
662 }
663 }
664 for(g = f->g, ge = g + f->ng; g < ge; g++) {
665 for(b = g->E, be = b + g->ns; b < be; b++) {
666 if ((e = b->D.ef)) {
667 hv_fwd(e);
668 e = b->D.ee;
669 e->aO = 0;
670 e->adO = t2*g->g1;
671 hv_back(e);
672 }
673 else if ((e = b->D.e)->op != f_OPNUM) {
674 e->aO = 0;
675 e->adO = t2*g->g1;
676 }
677 }
678 if (g->g2) {
679 t = 0.;
680 for(og = g->og; og; og = og->next)
681 t += og->coef * p[og->varno];
682 t *= t2*g->g2;
683 for(og = g->og; og; og = og->next)
684 x0[og->varno].aO += t*og->coef;
685 }
686 }
687 }
688 }
689 if (asl->P.ncom) {
690 for(c = cexps; c < ce--; ) {
691 for(i0 = *--dvsp0; i0 < i1; ) {
692 v = asl->P.vp[--i1];
693 --c1;
694 if ((t = v->aO) && (L = c1->L))
695 for(Le = L + c1->nlin; L < Le; L++)
696 ((expr_v*)L->v.vp)->aO += t * L->fac;
697 if (c1->hfun)
698 funnel_back(asl, c1, v, t);
699 else if ((e = c1->ee)) {
700 e->aO = t;
701 e->adO = v->adO;
702 hv_back(e);
703 }
704 else if ((e = c1->e)->op != f_OPNUM) {
705 e->aO = t;
706 e->adO = v->adO;
707 }
708 }
709 if ((t = (--x)->aO) && (L = ce->L))
710 for(Le = L + ce->nlin; L < Le; L++)
711 ((expr_v*)L->v.vp)->aO += t * L->fac;
712 if (ce->hfun)
713 funnel_back(asl, ce, x, t);
714 else if ((e = ce->ee)) {
715 e->aO = t;
716 e->adO = x->adO;
717 hv_back(e);
718 }
719 else if ((e = ce->e)->op != f_OPNUM) {
720 e->aO = t;
721 e->adO = x->adO;
722 }
723 }
724 }
725 x = var_e;
726 for(la = asl->P.lalist; la; la = la->lnext)
727 if ((t = la->v->aO)) {
728 og = la->nz;
729 do x[og->varno].aO += t*og->coef;
730 while((og = og->next));
731 }
732 while(x < xe)
733 *hv++ = (x++)->aO;
734 }
735
736 static real * /* Compute vector x0 = mat(h)*y0, */
737 /* where h = upper triang of mat(h). */
dtmul(int n,real * x0,real * h,real * y0)738 dtmul(int n, real *x0, real *h, real *y0)
739 {
740 int i;
741 real *hi, t, *x, *y, *y1, yi;
742
743 y1 = y0;
744 --h;
745 for(i = 0; i < n; i++) {
746 hi = ++h + i;
747 yi = *y1++;
748 t = yi**hi;
749 x = x0;
750 y = y0;
751 while(h < hi) {
752 t += *y++**h;
753 *x++ += yi**h++;
754 }
755 *x = t;
756 }
757 return x0;
758 }
759
760 void
hvpcomp_ASL(ASL * a,real * hv,real * p,int nobj,real * ow,real * y)761 hvpcomp_ASL(ASL *a, real *hv, real *p, int nobj, real *ow, real *y)
762 /* p = direction */
763 /* y = Lagrange multipliers */
764 /* hv = result */
765 {
766 ASL_pfgh *asl;
767 Ihinfo *ihi;
768 expr_v *v;
769 int kp, kw, n, no, noe, ns, nv, *ui, *uie;
770 linarg *la, **lap, **lape;
771 ograd *og;
772 ps_func *ps, *pe;
773 psg_elem *g, *ge;
774 range *r;
775 real *cscale, *owi, t, t1, t2, *p0, *s, *w, *wi, *x;
776
777 ASL_CHECK(a, ASL_read_pfgh, "hvpcomp");
778 asl = (ASL_pfgh*)a;
779 xpsg_check_ASL(asl, nobj, ow, y);
780 nv = n_var;
781 kp = htcl(nv*sizeof(real));
782 p0 = 0;
783 if ((s = asl->i.vscale)) {
784 p0 = (real*)new_mblk(kp);
785 for(n = 0; n < nv; n++)
786 p0[n] = s[n] * p[n];
787 p = p0;
788 }
789 if (!asl->P.ihdcur) {
790 if (asl->P.ndhmax <= 0) {
791 hvp0comp_ASL(asl,hv,p,nobj,ow,y);
792 goto done;
793 }
794 if (!(n = asl->P.nhvprod))
795 n = asl->P.ihdmin;
796 if (n >= asl->P.ihdmin)
797 hvpinit_ASL(a, ihd_limit, nobj, ow, y);
798 }
799 asl->P.nhvprod++;
800 memset(hv, 0, nv*sizeof(real));
801 for(la = asl->P.lalist; la; la = la->lnext) {
802 og = la->nz;
803 t = p[og->varno]*og->coef;
804 while((og = og->next))
805 t += p[og->varno]*og->coef;
806 v = la->v;
807 v->dO.r = t;
808 v->aO = v->adO = 0;
809 }
810 kw = kp + 1;
811 w = (real*)new_mblk(kw);
812 x = w + n_var;
813 s = asl->P.dOscratch;
814 ns = 0;
815 for(ihi = asl->P.ihi1; ihi; ihi = ihi->next) {
816 r = ihi->r;
817 if (ihi->hest)
818 for(; r; r = r->rlist.prev) {
819 n = r->n;
820 nv = r->nv;
821 wi = w;
822 if (n < nv) {
823 lap = r->lap;
824 lape = lap + n;
825 do {
826 og = (*lap++)->nz;
827 t = p[og->varno]*og->coef;
828 while((og = og->next))
829 t += p[og->varno]*og->coef;
830 *wi++ = t;
831 }
832 while(lap < lape);
833 wi = dtmul(n, x, r->hest, w);
834 lap = r->lap;
835 do if ((t = *wi++)) {
836 og = (*lap)->nz;
837 do hv[og->varno] +=
838 t*og->coef;
839 while((og = og->next));
840 }
841 while(++lap < lape);
842 }
843 else {
844 ui = r->ui;
845 uie = ui + nv;
846 do *wi++ = p[*ui++];
847 while(ui < uie);
848 wi = dtmul(nv, x, r->hest, w);
849 ui = r->ui;
850 do hv[*ui++] += *wi++;
851 while(ui < uie);
852 }
853 }
854 else
855 for(; r; r = r->rlist.prev) {
856 n = r->n;
857 if (ns < n)
858 ns = n;
859 wi = s;
860 lap = r->lap;
861 lape = lap + n;
862 do {
863 og = (*lap++)->nz;
864 t = p[og->varno]*og->coef;
865 while((og = og->next))
866 t += p[og->varno]*og->coef;
867 *wi++ = t;
868 }
869 while(lap < lape);
870 pshv_prod_ASL(asl, r, nobj, ow, y);
871 lap = r->lap;
872 do {
873 la = *lap++;
874 if ((t = la->v->aO)) {
875 og = la->nz;
876 do hv[og->varno] += t*og->coef;
877 while((og = og->next));
878 }
879 }
880 while(lap < lape);
881 }
882 }
883 del_mblk(kw,w);
884 wi = s + ns;
885 while(wi > s)
886 *--wi = 0.;
887 if (asl->P.nobjgroups) {
888 if (nobj >= 0 && nobj < n_obj) {
889 owi = ow ? ow + nobj : &edag_one_ASL;
890 no = nobj;
891 noe = no + 1;
892 }
893 else {
894 nobj = -1;
895 no = noe = 0;
896 if ((owi = ow))
897 noe = n_obj;
898 }
899 for(; no < noe; no++)
900 if ((t = *owi++)) {
901 ps = asl->P.ops + no;
902 g = ps->g;
903 for(ge = g + ps->ng; g < ge; g++)
904 if ((t2 = g->g2) && (og = g->og)) {
905 t1 = p[og->varno]*og->coef;
906 while((og = og->next))
907 t1 += p[og->varno]*og->coef;
908 t2 *= t*t1;
909 og = g->og;
910 do hv[og->varno] += t2*og->coef;
911 while((og = og->next));
912 }
913 }
914 }
915 if (asl->P.ncongroups && y) {
916 cscale = a->i.lscale;
917 ps = asl->P.cps;
918 for(pe = ps + n_con; ps < pe; ps++, y++)
919 if ((t = cscale ? *cscale++ * *y : *y))
920 for(g = ps->g, ge = g + ps->ng; g < ge; g++)
921 if ((t2 = g->g2) && (og = g->og)) {
922 t1 = p[og->varno]*og->coef;
923 while((og = og->next))
924 t1 += p[og->varno]*og->coef;
925 t2 *= t*t1;
926 og = g->og;
927 do hv[og->varno] += t2*og->coef;
928 while((og = og->next));
929 }
930 }
931 v = var_e;
932 for(la = asl->P.lalist; la; la = la->lnext)
933 if ((t = la->v->aO)) {
934 og = la->nz;
935 do v[og->varno].aO += t*og->coef;
936 while((og = og->next));
937 }
938 done:
939 if (p0) {
940 del_mblk(kp, p0);
941 s = asl->i.vscale;
942 w = hv + n_var;
943 while(hv < w)
944 *hv++ *= *s++;
945 }
946 }
947
948 void
pshv_prod_ASL(ASL_pfgh * asl,range * r,int nobj,real * ow,real * y)949 pshv_prod_ASL(ASL_pfgh *asl, range *r, int nobj, real *ow, real *y)
950 {
951 cexp *c;
952 expr *e;
953 expr_v *v;
954 int *cei, *cei0, *ceie, i;
955 linarg *la, **lap, **lape;
956 linpart *L, *Le;
957 ps_func *p;
958 psb_elem *b;
959 psg_elem *g;
960 real *cscale, *s, owi, t;
961
962 cscale = asl->i.lscale;
963 owi = 1.;
964 if (nobj >= 0 && nobj < n_obj) {
965 if (ow) {
966 if ((owi = ow[nobj]) == 0.)
967 nobj = -1;
968 ow = 0;
969 }
970 }
971 if (x0kind & ASL_need_funnel)
972 funnelhes(asl);
973 s = asl->P.dOscratch;
974 lap = r->lap;
975 lape = lap + r->n;
976 while(lap < lape) {
977 la = *lap++;
978 v = la->v;
979 v->dO.r = *s++;
980 v->adO = v->aO = 0.;
981 }
982 if ((cei = cei0 = r->cei)) {
983 i = *cei0++;
984 ceie = (cei = cei0) + i;
985 do {
986 i = *cei++;
987 hv_fwd0(asl, cexps + i, asl->P.vp[i]);
988 }
989 while(cei < ceie);
990 }
991 for(b = r->refs; b; b = b->next) {
992 if ((i = b->conno) < 0) {
993 i = -2 - i;
994 if (i == nobj)
995 t = owi;
996 else if (ow) {
997 if (!(t = ow[i]))
998 continue;
999 }
1000 else
1001 continue;
1002 p = asl->P.ops;
1003 }
1004 else {
1005 if (!y || !(t = y[i]))
1006 continue;
1007 if (cscale)
1008 t *= cscale[i];
1009 p = asl->P.cps;
1010 }
1011 if (b->groupno) {
1012 p += i;
1013 g = p->g + (b->groupno - 1);
1014 if (asl->P.pshv_g1)
1015 t *= g->g1;
1016 }
1017 if ((e = b->D.ef)) {
1018 hv_fwd(e);
1019 e = b->D.ee;
1020 e->aO = 0;
1021 e->adO = t;
1022 hv_back(e);
1023 }
1024 else if ((e = b->D.e)->op != f_OPNUM)
1025 e->adO += t;
1026 }
1027 while(cei > cei0) {
1028 i = *--cei;
1029 c = cexps + i;
1030 v = asl->P.vp[i];
1031 if ((t = v->aO) && (L = c->L)) {
1032 if ((la = c->la))
1033 la->v->aO += t * asl->P.dv[i].scale;
1034 else {
1035 for(Le = L + c->nlin; L < Le; L++)
1036 ((expr_v*)L->v.vp)->aO += t * L->fac;
1037 }
1038 }
1039 if (c->hfun)
1040 funnel_back(asl, c, v, t);
1041 else if ((e = c->ee)) {
1042 e->aO = t;
1043 e->adO = v->adO;
1044 hv_back(e);
1045 }
1046 else if ((e = c->e)->op != f_OPNUM) {
1047 e->aO += t;
1048 e->adO += v->adO;
1049 }
1050 }
1051 }
1052
1053 void
funpset_ASL(ASL_pfgh * asl,funnel * f)1054 funpset_ASL(ASL_pfgh *asl, funnel *f)
1055 {
1056 cplist *cl;
1057 derp *d;
1058
1059 for(; f; f = f->next) {
1060 memset(adjoints_nv1, 0, f->fcde.zaplen);
1061 cl = f->cl;
1062 do *cl->ca.rp = 0;
1063 while((cl = cl->next));
1064 d = f->fcde.d;
1065 *d->b.rp = 1.;
1066 do *d->a.rp += *d->b.rp * *d->c.rp;
1067 while((d = d->next));
1068 cl = f->cl;
1069 do *cl->cfa = *cl->ca.rp;
1070 while((cl = cl->next));
1071 }
1072 }
1073
1074 void
hvpcompd_ASL(ASL * a,real * hv,real * p,int co)1075 hvpcompd_ASL(ASL *a, real *hv, real *p, int co)
1076 /* p = direction */
1077 /* hv = result */
1078 /* co >= 0: behave like hvpcomp_ASL with nobj = -1, ow = 0, y[i] = i == co ? 1. : 0. */
1079 /* co < 0: behave like hvpcomp_ASL with nobj = -1 - co, ow = 0, y = 0 */
1080 {
1081 ASL_pfgh *asl;
1082 cexp *c, *c1, *ce;
1083 cgrad *cg, *cg0;
1084 expr *e;
1085 expr_v *v, *x, *x0;
1086 int *dvsp0, i0, i1, kp, n, no, nx, oxk;
1087 linarg *la;
1088 linpart *L, *Le;
1089 ograd *og, *og0;
1090 ps_func *f;
1091 psb_elem *b, *be;
1092 psg_elem *g, *ge;
1093 real *s, t, t2, *p0;
1094 varno_t i;
1095
1096 #ifdef IGNORE_BOGUS_WARNINGS
1097 c1 = ce = 0;
1098 dvsp0 = 0;
1099 i1 = 0;
1100 v = 0;
1101 x = 0;
1102 #endif
1103 ASL_CHECK(a, ASL_read_pfgh, "hvpcompi");
1104 asl = (ASL_pfgh*)a;
1105 if (x0kind == ASL_first_x) {
1106 if (!(s = X0))
1107 memset(s = Lastx, 0, n_var*sizeof(real));
1108 co_index = co;
1109 xp_check_ASL(asl, s);
1110 }
1111 nx = asl->i.nxval;
1112 oxk = asl->i.x_known;
1113 asl->i.x_known = 1;
1114
1115 p0 = 0;
1116 cg0 = 0;
1117 og0 = 0;
1118 x0 = var_e;
1119 n = c_vars >= o_vars ? c_vars : o_vars;
1120 t2 = 1.;
1121 memset(hv, 0, n_var*sizeof(real));
1122 for(la = asl->P.lalist; la; la = la->lnext) {
1123 og = la->nz;
1124 t = p[og->varno]*og->coef;
1125 while((og = og->next))
1126 t += p[og->varno]*og->coef;
1127 x = la->v;
1128 x->dO.r = t;
1129 x->aO = x->adO = 0;
1130 }
1131 if (co >= 0) {
1132 if (co >= nlc)
1133 return;
1134 f = asl->P.cps + co;
1135 if (asl->i.ncxval[co] != nx)
1136 conpival_ASL(a, co, Lastx, 0);
1137 if (f->ng && f->nxval != nx)
1138 conpgrd_ASL(a, co, Lastx, 0, 0);
1139 if ((s = asl->i.lscale))
1140 t2 = s[co];
1141 cg = cg0 = Cgrad[co];
1142 if ((s = asl->i.vscale)) {
1143 kp = htcl(n*sizeof(real));
1144 p0 = (real*)new_mblk(kp);
1145 for(; cg; cg = cg->next) {
1146 i = cg->varno;
1147 x = x0 + i;
1148 x->dO.r = p0[i] = p[i]*s[i];
1149 x->aO = x->adO = 0.;
1150 }
1151 p = p0;
1152 }
1153 else {
1154 for(; cg; cg = cg->next) {
1155 i = cg->varno;
1156 x = x0 + i;
1157 x->dO.r = p[i];
1158 x->aO = x->adO = 0.;
1159 }
1160 }
1161 }
1162 else {
1163 no = -1 - co;
1164 if (no >= nlo)
1165 return;
1166 f = asl->P.ops + no;
1167 if (asl->i.ncxval[no] != nx)
1168 objpval_ASL(a, no, Lastx, 0);
1169 if (f->ng && f->nxval != nx)
1170 objpgrd_ASL(a, no, Lastx, 0, 0);
1171 og = og0 = Ograd[no];
1172 if ((s = asl->i.vscale)) {
1173 kp = htcl(n*sizeof(real));
1174 p0 = (real*)new_mblk(kp);
1175 for(; og; og = og->next) {
1176 i = og->varno;
1177 x = x0 + i;
1178 x->dO.r = p0[i] = p[i]*s[i];
1179 x->aO = x->adO = 0.;
1180 }
1181 p = p0;
1182 }
1183 else {
1184 for(; og; og = og->next) {
1185 i = og->varno;
1186 x = x0 + i;
1187 x->dO.r = p[i];
1188 x->aO = x->adO = 0.;
1189 }
1190 }
1191 }
1192 if (asl->i.Derrs) {
1193 asl->i.x_known = oxk;
1194 deriv_errchk_ASL(a, 0, co, 1);
1195 asl->i.x_known = 1;
1196 }
1197 if (asl->P.ncom) {
1198 x = var_ex;
1199 dvsp0 = asl->P.dvsp0;
1200 c = cexps;
1201 i0 = *dvsp0;
1202 for(ce = c1 = c + asl->P.ncom; c < ce; c++) {
1203 for(i1 = *++dvsp0; i0 < i1; i0++)
1204 hv_fwd0(asl, c1++, asl->P.vp[i0]);
1205 hv_fwd0(asl, c, x++);
1206 }
1207 }
1208 for(b = f->b, be = b + f->nb; b < be; b++) {
1209 if ((e = b->D.ef)) {
1210 hv_fwd(e);
1211 e = b->D.ee;
1212 e->aO = 0;
1213 e->adO = t2;
1214 hv_back(e);
1215 }
1216 else if ((e = b->D.e)->op != f_OPNUM) {
1217 e->aO = 0;
1218 e->adO = t2;
1219 }
1220 }
1221 for(g = f->g, ge = g + f->ng; g < ge; g++) {
1222 for(b = g->E, be = b + g->ns; b < be; b++) {
1223 if ((e = b->D.ef)) {
1224 hv_fwd(e);
1225 e = b->D.ee;
1226 e->aO = 0;
1227 e->adO = t2*g->g1;
1228 hv_back(e);
1229 }
1230 else if ((e = b->D.e)->op != f_OPNUM) {
1231 e->aO = 0;
1232 e->adO = t2*g->g1;
1233 }
1234 }
1235 if (g->g2) {
1236 t = 0.;
1237 for(og = g->og; og; og = og->next)
1238 t += og->coef * p[og->varno];
1239 t *= t2*g->g2;
1240 for(og = g->og; og; og = og->next)
1241 x0[og->varno].aO += t*og->coef;
1242 }
1243 }
1244 if (asl->P.ncom) {
1245 for(c = cexps; c < ce--; ) {
1246 for(i0 = *--dvsp0; i0 < i1; ) {
1247 v = asl->P.vp[--i1];
1248 --c1;
1249 if ((t = v->aO) && (L = c1->L))
1250 for(Le = L + c1->nlin; L < Le; L++)
1251 ((expr_v*)L->v.vp)->aO += t * L->fac;
1252 if (c1->hfun)
1253 funnel_back(asl, c1, v, t);
1254 else if ((e = c1->ee)) {
1255 e->aO = t;
1256 e->adO = v->adO;
1257 hv_back(e);
1258 }
1259 else if ((e = c1->e)->op != f_OPNUM) {
1260 e->aO = t;
1261 e->adO = v->adO;
1262 }
1263 }
1264 if ((t = (--x)->aO) && (L = ce->L))
1265 for(Le = L + ce->nlin; L < Le; L++)
1266 ((expr_v*)L->v.vp)->aO += t * L->fac;
1267 if (ce->hfun)
1268 funnel_back(asl, ce, x, t);
1269 else if ((e = ce->ee)) {
1270 e->aO = t;
1271 e->adO = x->adO;
1272 hv_back(e);
1273 }
1274 else if ((e = ce->e)->op != f_OPNUM) {
1275 e->aO = t;
1276 e->adO = x->adO;
1277 }
1278 }
1279 }
1280 x = var_e;
1281 for(la = asl->P.lalist; la; la = la->lnext)
1282 if ((t = la->v->aO)) {
1283 og = la->nz;
1284 do x[og->varno].aO += t*og->coef;
1285 while((og = og->next));
1286 }
1287 if ((cg = cg0)) {
1288 if (s) {
1289 while(cg) {
1290 i = cg->varno;
1291 hv[i] = s[i]*x0[i].aO;
1292 cg = cg->next;
1293 }
1294 }
1295 else {
1296 while(cg) {
1297 i = cg->varno;
1298 hv[i] = x0[i].aO;
1299 cg = cg->next;
1300 }
1301 }
1302 }
1303 else {
1304 og = og0;
1305 if (s) {
1306 while(og) {
1307 i = og->varno;
1308 hv[i] = s[i]*x0[i].aO;
1309 og = og->next;
1310 }
1311 }
1312 else {
1313 while(og) {
1314 i = og->varno;
1315 hv[i] = x0[i].aO;
1316 og = og->next;
1317 }
1318 }
1319 }
1320 if (p0)
1321 del_mblk(kp, p0);
1322 }
1323
1324 varno_t
hvpcomps_ASL(ASL * a,real * hv,real * p,int co,varno_t nz,varno_t * z)1325 hvpcomps_ASL(ASL *a, real *hv, real *p, int co, varno_t nz, varno_t *z)
1326 /* p = direction */
1327 /* hv = result */
1328 /* co >= 0: behave like hvpcomp_ASL with nobj = -1, ow = 0, y[i] = i == co ? 1. : 0. */
1329 /* co < 0: behave like hvpcomp_ASL with nobj = -1 - co, ow = 0, y = 0 */
1330 {
1331 ASL_pfgh *asl;
1332 cexp *c, *c1, *ce;
1333 cgrad *cg, *cg0;
1334 expr *e;
1335 expr_v *v, *x, *x0;
1336 int *dvsp0, i0, i1, kp, n, no, nx, oxk;
1337 linarg *la;
1338 linpart *L, *Le;
1339 ograd *og, *og0;
1340 ps_func *f;
1341 psb_elem *b, *be;
1342 psg_elem *g, *ge;
1343 real *hve, *p0, *s, t, t2, *vscale;
1344 varno_t F, i, rv, *ze;
1345
1346 #ifdef IGNORE_BOGUS_WARNINGS
1347 c1 = ce = 0;
1348 dvsp0 = 0;
1349 i1 = 0;
1350 v = 0;
1351 x = 0;
1352 #endif
1353 ASL_CHECK(a, ASL_read_pfgh, "hvpcompi");
1354 asl = (ASL_pfgh*)a;
1355 if (x0kind == ASL_first_x) {
1356 if (!(s = X0))
1357 memset(s = Lastx, 0, n_var*sizeof(real));
1358 co_index = co;
1359 xp_check_ASL(asl, s);
1360 }
1361 nx = asl->i.nxval;
1362 oxk = asl->i.x_known;
1363 asl->i.x_known = 1;
1364
1365 p0 = 0;
1366 cg0 = 0;
1367 og0 = 0;
1368 x0 = var_e;
1369 n = c_vars >= o_vars ? c_vars : o_vars;
1370 t2 = 1.;
1371 memset(hv, 0, n_var*sizeof(real));
1372 for(la = asl->P.lalist; la; la = la->lnext) {
1373 og = la->nz;
1374 t = p[og->varno]*og->coef;
1375 while((og = og->next))
1376 t += p[og->varno]*og->coef;
1377 x = la->v;
1378 x->dO.r = t;
1379 x->aO = x->adO = 0;
1380 }
1381 no = -1 - co;
1382 if (co >= 0) {
1383 if (co >= nlc)
1384 return 0;
1385 f = asl->P.cps + co;
1386 if (asl->i.ncxval[co] != nx)
1387 conpival_ASL(a, co, Lastx, 0);
1388 if (f->ng && f->nxval != nx)
1389 conpgrd_ASL(a, co, Lastx, 0, 0);
1390 if ((s = asl->i.lscale))
1391 t2 = s[co];
1392 cg = cg0 = Cgrad[co];
1393 if ((vscale = asl->i.vscale)) {
1394 kp = htcl(n*sizeof(real));
1395 p0 = (real*)new_mblk(kp);
1396 for(; cg; cg = cg->next) {
1397 i = cg->varno;
1398 x = x0 + i;
1399 x->dO.r = p0[i] = p[i]*vscale[i];
1400 x->aO = x->adO = 0.;
1401 }
1402 p = p0;
1403 }
1404 else {
1405 for(; cg; cg = cg->next) {
1406 i = cg->varno;
1407 x = x0 + i;
1408 x->dO.r = p[i];
1409 x->aO = x->adO = 0.;
1410 }
1411 }
1412 }
1413 else {
1414 if (no >= nlo)
1415 return 0;
1416 f = asl->P.ops + no;
1417 if (asl->i.ncxval[no] != nx)
1418 objpval_ASL(a, no, Lastx, 0);
1419 if (f->ng && f->nxval != nx)
1420 objpgrd_ASL(a, no, Lastx, 0, 0);
1421 og = og0 = Ograd[no];
1422 if ((vscale = asl->i.vscale)) {
1423 kp = htcl(n*sizeof(real));
1424 p0 = (real*)new_mblk(kp);
1425 for(; og; og = og->next) {
1426 i = og->varno;
1427 x = x0 + i;
1428 x->dO.r = p0[i] = p[i]*vscale[i];
1429 x->aO = x->adO = 0.;
1430 }
1431 p = p0;
1432 }
1433 else {
1434 for(; og; og = og->next) {
1435 i = og->varno;
1436 x = x0 + i;
1437 x->dO.r = p[i];
1438 x->aO = x->adO = 0.;
1439 }
1440 }
1441 }
1442 if (asl->i.Derrs) {
1443 asl->i.x_known = oxk;
1444 deriv_errchk_ASL(a, 0, co, 1);
1445 asl->i.x_known = 1;
1446 }
1447 if (asl->P.ncom) {
1448 x = var_ex;
1449 dvsp0 = asl->P.dvsp0;
1450 c = cexps;
1451 i0 = *dvsp0;
1452 for(ce = c1 = c + asl->P.ncom; c < ce; c++) {
1453 for(i1 = *++dvsp0; i0 < i1; i0++)
1454 hv_fwd0(asl, c1++, asl->P.vp[i0]);
1455 hv_fwd0(asl, c, x++);
1456 }
1457 }
1458 for(b = f->b, be = b + f->nb; b < be; b++) {
1459 if ((e = b->D.ef)) {
1460 hv_fwd(e);
1461 e = b->D.ee;
1462 e->aO = 0;
1463 e->adO = t2;
1464 hv_back(e);
1465 }
1466 else if ((e = b->D.e)->op != f_OPNUM) {
1467 e->aO = 0;
1468 e->adO = t2;
1469 }
1470 }
1471 for(g = f->g, ge = g + f->ng; g < ge; g++) {
1472 for(b = g->E, be = b + g->ns; b < be; b++) {
1473 if ((e = b->D.ef)) {
1474 hv_fwd(e);
1475 e = b->D.ee;
1476 e->aO = 0;
1477 e->adO = t2*g->g1;
1478 hv_back(e);
1479 }
1480 else if ((e = b->D.e)->op != f_OPNUM) {
1481 e->aO = 0;
1482 e->adO = t2*g->g1;
1483 }
1484 }
1485 if (g->g2) {
1486 t = 0.;
1487 for(og = g->og; og; og = og->next)
1488 t += og->coef * p[og->varno];
1489 t *= t2*g->g2;
1490 for(og = g->og; og; og = og->next)
1491 x0[og->varno].aO += t*og->coef;
1492 }
1493 }
1494 if (asl->P.ncom) {
1495 for(c = cexps; c < ce--; ) {
1496 for(i0 = *--dvsp0; i0 < i1; ) {
1497 v = asl->P.vp[--i1];
1498 --c1;
1499 if ((t = v->aO) && (L = c1->L))
1500 for(Le = L + c1->nlin; L < Le; L++)
1501 ((expr_v*)L->v.vp)->aO += t * L->fac;
1502 if (c1->hfun)
1503 funnel_back(asl, c1, v, t);
1504 else if ((e = c1->ee)) {
1505 e->aO = t;
1506 e->adO = v->adO;
1507 hv_back(e);
1508 }
1509 else if ((e = c1->e)->op != f_OPNUM) {
1510 e->aO = t;
1511 e->adO = v->adO;
1512 }
1513 }
1514 if ((t = (--x)->aO) && (L = ce->L))
1515 for(Le = L + ce->nlin; L < Le; L++)
1516 ((expr_v*)L->v.vp)->aO += t * L->fac;
1517 if (ce->hfun)
1518 funnel_back(asl, ce, x, t);
1519 else if ((e = ce->ee)) {
1520 e->aO = t;
1521 e->adO = x->adO;
1522 hv_back(e);
1523 }
1524 else if ((e = ce->e)->op != f_OPNUM) {
1525 e->aO = t;
1526 e->adO = x->adO;
1527 }
1528 }
1529 }
1530 x = var_e;
1531 for(la = asl->P.lalist; la; la = la->lnext)
1532 if ((t = la->v->aO)) {
1533 og = la->nz;
1534 do x[og->varno].aO += t*og->coef;
1535 while((og = og->next));
1536 }
1537 rv = 0;
1538 if ((ze = z))
1539 ze += nz;
1540 if ((hve = hv))
1541 hve += nz;
1542 F = Fortran;
1543 if ((cg = cg0)) {
1544 if (!hv) {
1545 while(cg) {
1546 ++rv;
1547 if (z < ze)
1548 *z++ = cg->varno;
1549 cg = cg->next;
1550 }
1551 }
1552 else if (vscale) {
1553 while(cg) {
1554 ++rv;
1555 i = cg->varno;
1556 if (z < ze)
1557 *z++ = F + i;
1558 if (hv < hve)
1559 *hv++ = vscale[i]*x0[i].aO;
1560 cg = cg->next;
1561 }
1562 }
1563 else {
1564 while(cg) {
1565 ++rv;
1566 i = cg->varno;
1567 if (z < ze)
1568 *z++ = F + i;
1569 if (hv < hve)
1570 *hv++ = x0[i].aO;
1571 cg = cg->next;
1572 }
1573 }
1574 }
1575 else {
1576 og = Ograd[no];
1577 if (!hv) {
1578 while(og) {
1579 ++rv;
1580 if (z < ze)
1581 *z++ = og->varno;
1582 og = og->next;
1583 }
1584 }
1585 else if (vscale) {
1586 while(og) {
1587 ++rv;
1588 i = og->varno;
1589 if (z < ze)
1590 *z++ = F + i;
1591 if (hv < hve)
1592 *hv++ = vscale[i]*x0[i].aO;
1593 og = og->next;
1594 }
1595 }
1596 else {
1597 while(og) {
1598 ++rv;
1599 i = og->varno;
1600 if (z < ze)
1601 *z++ = F + i;
1602 if (hv < hve)
1603 *hv++ = x0[i].aO;
1604 og = og->next;
1605 }
1606 }
1607 }
1608 if (p0)
1609 del_mblk(kp, p0);
1610 return rv;
1611 }
1612
1613 #ifdef __cplusplus
1614 }
1615 #endif
1616