1 /****************************************************************
2 Copyright (C) 1997-1998, 2000-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 "asl_pfgh.h"
26
27 static void
28 #ifdef KR_headers
hv_fwd(asl,e)29 hv_fwd(asl, e) ASL_pfgh *asl; register expr *e;
30 #else
31 hv_fwd(ASL_pfgh *asl, register expr *e)
32 #endif
33 {
34 expr *e1, **ep;
35 real dO;
36 de *d;
37 derp *D, *D1;
38
39 for(; e; e = e->fwd) {
40 e->aO = e->adO = 0;
41 switch(e->a) {
42
43 case Hv_timesR:
44 case Hv_binaryR:
45 e->dO.r = e->R.e->dO.r;
46 break;
47
48 case Hv_timesLR:
49 case Hv_binaryLR:
50 e->dO.r = e->L.e->dO.r + e->R.e->dO.r;
51 break;
52
53 case Hv_timesL:
54 case Hv_unary:
55 e->dO.r = e->L.e->dO.r;
56 break;
57
58 case Hv_vararg:
59 e->dO.r = 0;
60 for(d = ((expr_va*)e)->L.d; d->e; d++) {
61 if (e1 = d->ef) {
62 hv_fwd(asl, e1);
63 e->dO.r += d->ee->dO.r;
64 }
65 else if (e1 = d->e)
66 e->dO.r += e1->dO.r;
67 }
68 if (!((expr_va*)e)->val && (D = ((expr_va*)e)->R.D)) {
69 D->a.rp = adjoints_nv1;
70 D1 = ((expr_va*)e)->d0;
71 d = ((expr_va*)e)->L.d;
72 ((expr_va*)e)->val = d->e;
73 do {
74 D->next = d->d;
75 while(D->next != D1)
76 D = D->next;
77 d->dlast = D;
78 ++d;
79 }
80 while(d->e);
81 ((expr_va*)e)->next = asl->P.valist;
82 asl->P.valist = (expr_va*)e;
83 }
84 break;
85
86 case Hv_if:
87 e->dO.r = 0;
88 if (e1 = ((expr_if*)e)->Tf) {
89 hv_fwd(asl, e1);
90 e->dO.r = ((expr_if *)e)->Te->dO.r;
91 }
92 else if (e1 = ((expr_if*)e)->T)
93 e->dO.r = e1->dO.r;
94 if (e1 = ((expr_if*)e)->Ff) {
95 hv_fwd(asl, e1);
96 e->dO.r += ((expr_if *)e)->Fe->dO.r;
97 }
98 else if (e1 = ((expr_if*)e)->F)
99 e->dO.r += e1->dO.r;
100 if (!((expr_if*)e)->val && (D = ((expr_if*)e)->D)) {
101 ((expr_if*)e)->val = ((expr_if*)e)->T;
102 D->a.rp = adjoints_nv1;
103 D1 = ((expr_if*)e)->d0;
104 D->next = ((expr_if*)e)->dT;
105 while(D->next != D1)
106 D = D->next;
107 ((expr_if*)e)->dTlast = D;
108 D->next = ((expr_if*)e)->dF;
109 ((expr_if*)e)->next = asl->P.iflist;
110 asl->P.iflist = (expr_if*)e;
111 }
112 break;
113
114 case Hv_plterm:
115 e->dO.r = e->R.e->dO.r;
116 break;
117
118 case Hv_sumlist:
119 ep = e->R.ep;
120 for(dO = 0; e1 = *ep; ep++)
121 dO += e1->dO.r;
122 e->dO.r = dO;
123 break;
124
125 case Hv_func:
126 e->dO.r = 1.;
127 break;
128
129 case Hv_negate:
130 e->dO.r = e->L.e->dO.r;
131 break;
132
133 case Hv_plusR:
134 e->dO.r = e->R.e->dO.r;
135 break;
136
137 case Hv_plusL:
138 e->dO.r = e->L.e->dO.r;
139 break;
140
141 case Hv_plusLR:
142 e->dO.r = e->L.e->dO.r + e->R.e->dO.r;
143 break;
144
145 case Hv_minusR:
146 e->dO.r = e->R.e->dO.r;
147 break;
148
149 case Hv_minusLR:
150 e->dO.r = e->L.e->dO.r + e->R.e->dO.r;
151 break;
152
153 default:/*DEBUG*/
154 fprintf(Stderr, "bad e->a = %d in hv_fwd\n", e->a);
155 exit(1);
156 }
157 }
158 }
159
160 static void
161 #ifdef KR_headers
func_back(f)162 func_back(f) expr_f *f;
163 #else
164 func_back(expr_f *f)
165 #endif
166 {
167 argpair *da, *da1, *dae;
168 expr *e;
169 real aO, adO, t;
170
171 aO = f->aO;
172 adO = f->adO;
173 da = f->da;
174 for(dae = f->dae; da < dae; da++) {
175 e = da->e;
176 e->aO += aO;
177 e->adO += adO;
178 t = adO*e->dO.r;
179 for(da1 = f->da; da1 < dae; da1++) {
180 e = da1->e;
181 e->aO += t;
182 }
183 }
184 }
185
186 static void
187 #ifdef KR_headers
hv_back(e)188 hv_back(e) register expr *e;
189 #else
190 hv_back(register expr *e)
191 #endif
192 {
193 register expr *e1, **ep, *e2;
194 real adO, t1, t2;
195 de *d;
196
197 if (!e || !e->aO && !e->adO)
198 return;
199 for(; e; e = e->bak)
200 switch(e->a) {
201 case Hv_binaryR:
202 e1 = e->R.e;
203 e1->adO += e->adO;
204 e1->aO += e->aO + e->adO * e1->dO.r;
205 break;
206
207 case Hv_binaryLR:
208 e1 = e->L.e;
209 e2 = e->R.e;
210 adO = e->adO;
211 t1 = adO * e1->dO.r;
212 t2 = adO * e2->dO.r;
213 e1->aO += e->aO + t1 + t2;
214 e2->aO += e->aO + t1 + t2;
215 e1->adO += adO;
216 e2->adO += adO;
217 break;
218
219 case Hv_unary:
220 e1 = e->L.e;
221 e1->adO += e->adO;
222 e1->aO += e->aO + e->adO * e1->dO.r;
223 break;
224
225 case Hv_vararg:
226 for(d = ((expr_va*)e)->L.d; d->e; d++) {
227 if (e1 = d->ee) {
228 e1->aO = e->aO;
229 e1->adO = e->adO;
230 hv_back(e1);
231 }
232 else {
233 e1 = d->e;
234 if (e1->op != f_OPNUM) {
235 e1->aO = e->aO;
236 e1->adO = e->adO;
237 }
238 }
239 }
240 break;
241
242 case Hv_if:
243 if (e1 = ((expr_if *)e)->Te) {
244 e1->aO = e->aO;
245 e1->adO = e->adO;
246 hv_back(e1);
247 }
248 else {
249 e1 = ((expr_if *)e)->T;
250 if (e1->op != f_OPNUM) {
251 e1->aO = e->aO;
252 e1->adO = e->adO;
253 }
254 }
255 if (e1 = ((expr_if *)e)->Fe) {
256 e1->aO = e->aO;
257 e1->adO = e->adO;
258 hv_back(e1);
259 }
260 else {
261 e1 = ((expr_if *)e)->F;
262 if (e1->op != f_OPNUM) {
263 e1->aO = e->aO;
264 e1->adO = e->adO;
265 }
266 }
267 break;
268
269 case Hv_plterm:
270 e->R.e->aO += e->aO;
271 break;
272
273 case Hv_sumlist:
274 ep = e->R.ep;
275 t1 = e->aO;
276 t2 = e->adO;
277 while(e1 = *ep++) {
278 e1->aO += t1;
279 e1->adO += t2;
280 }
281 break;
282
283 case Hv_func:
284 func_back((expr_f *)e);
285 break;
286
287 case Hv_negate:
288 e1 = e->L.e;
289 neg_end:
290 e1->aO += e->aO;
291 e1->adO += e->adO;
292 break;
293
294 case Hv_plusR:
295 e1 = e->R.e;
296 goto plus_end;
297
298 case Hv_plusL:
299 e1 = e->L.e;
300 plus_end:
301 e1->aO += e->aO;
302 e1->adO += e->adO;
303 break;
304
305 case Hv_plusLR:
306 e1 = e->L.e;
307 e1->aO += t1 = e->aO;
308 e1->adO += t2 = e->adO;
309 e1 = e->R.e;
310 e1->aO += t1;
311 e1->adO += t2;
312 break;
313
314 case Hv_minusR:
315 e1 = e->R.e;
316 goto neg_end;
317
318 case Hv_minusLR:
319 e1 = e->L.e;
320 e1->aO += t1 = e->aO;
321 e1->adO += t2 = e->adO;
322 e1 = e->R.e;
323 e1->aO += t1;
324 e1->adO += t2;
325 break;
326
327 case Hv_timesR:
328 e1 = e->R.e;
329 e1->adO += e->adO;
330 e1->aO += e->aO;
331 break;
332
333 case Hv_timesL:
334 e1 = e->L.e;
335 e1->adO += e->adO;
336 e1->aO += e->aO;
337 break;
338
339 case Hv_timesLR:
340 e1 = e->L.e;
341 e2 = e->R.e;
342 adO = e->adO;
343 e1->aO += e->aO + adO * e2->dO.r;
344 e2->aO += e->aO + adO * e1->dO.r;
345 e1->adO += adO;
346 e2->adO += adO;
347 break;
348
349 default:/*DEBUG*/
350 fprintf(Stderr, "bad e->a = %d in hv_back\n", e->a);
351 exit(1);
352 }
353 }
354
355 static void
356 #ifdef KR_headers
hv_fwd0(asl,c,v)357 hv_fwd0(asl, c, v) ASL_pfgh *asl; register cexp *c; register expr_v *v;
358 #else
359 hv_fwd0(ASL_pfgh *asl, register cexp *c, register expr_v *v)
360 #endif
361 {
362 register linpart *L, *Le;
363 real x;
364
365 v->aO = v->adO = 0;
366 if (c->ef) {
367 hv_fwd(asl, c->ef);
368 x = c->ee->dO.r;
369 }
370 else if (c->e->op != f_OPNUM)
371 x = c->e->dO.r;
372 else
373 x = 0;
374 if (L = c->L)
375 for(Le = L + c->nlin; L < Le; L++)
376 x += ((expr_v*)L->v.vp)->dO.r;
377 v->dO.r = x;
378 }
379
380 static void
381 #ifdef KR_headers
pshv_prod1(asl,r,nobj,ow,y)382 pshv_prod1(asl, r, nobj, ow, y) ASL_pfgh *asl; range *r; int nobj, ow, y;
383 #else
384 pshv_prod1(ASL_pfgh *asl, range *r, int nobj, int ow, int y)
385 #endif
386 {
387 int *cei, *cei0, *ceie, i;
388 linarg *la, **lap, **lap1, **lape;
389 linpart *L, *Le;
390 expr_v *v;
391 cexp *c;
392 expr *e;
393 real *s;
394 psb_elem *b;
395
396 s = asl->P.dOscratch;
397 lap = lap1 = r->lap;
398 lape = lap + r->n;
399 while(lap1 < lape) {
400 la = *lap1++;
401 v = la->v;
402 v->dO.r = *s++;
403 v->adO = v->aO = 0.;
404 }
405 if (cei = cei0 = r->cei) {
406 i = *cei0++;
407 ceie = (cei = cei0) + i;
408 do {
409 i = *cei++;
410 hv_fwd0(asl, cexps + i, asl->P.vp[i]);
411 }
412 while(cei < ceie);
413 }
414 for(b = r->refs; b; b = b->next) {
415 if ((i = b->conno) < 0) {
416 i = -2 - i;
417 if (!ow && i != nobj)
418 continue;
419 }
420 else if (!y)
421 continue;
422 if (e = b->D.ef) {
423 hv_fwd(asl, e);
424 e = b->D.ee;
425 e->aO = 0;
426 e->adO = 1.;
427 hv_back(e);
428 }
429 else if ((e = b->D.e)->op != f_OPNUM) {
430 e->aO = 0;
431 e->adO = 1.;
432 }
433 }
434 while(cei > cei0) {
435 i = *--cei;
436 c = cexps + i;
437 v = asl->P.vp[i];
438 if (v->aO && (L = c->L))
439 for(Le = L + c->nlin; L < Le; L++)
440 ((expr_v*)L->v.vp)->aO++;
441 if (e = c->ee) {
442 e->aO = 1.;
443 e->adO = v->adO;
444 hv_back(e);
445 }
446 else if ((e = c->e)->op != f_OPNUM) {
447 e->aO = 1.;
448 e->adO = v->adO;
449 }
450 }
451 }
452
453 #ifdef __cplusplus
454 extern "C" {
455 static int compar(const void*, const void*);
456 }
457 #endif
458
459 static int
460 #ifdef KR_headers
compar(a,b)461 compar(a, b) char *a, *b;
462 #else
463 compar(const void *a, const void *b)
464 #endif
465 { return *(int*)a - *(int*)b; }
466
467 #undef nzc
468 #undef asl
469 #undef del_mblk
470 #define del_mblk(b,c) Del_mblk_ASL(a,b,(Char*)(c))
471
472 static void
473 #ifdef KR_headers
new_Hesoprod(asl,L,R,coef)474 new_Hesoprod(asl, L, R, coef) ASL_pfgh *asl; ograd *L, *R; real coef;
475 #else
476 new_Hesoprod(ASL_pfgh *asl, ograd *L, ograd *R, real coef)
477 #endif
478 {
479 Hesoprod *h, **hp, *h1, *h2;
480 int kh;
481 Char **mblk_free;
482
483 ACQUIRE_DTOA_LOCK(HESOPROD_LOCK);
484 if (!(h = asl->P.hop_free)) {
485 mblk_free = asl->mblk_free;
486 kh = asl->P.khesoprod;
487 while(kh < 8 && !mblk_free[kh])
488 kh++;
489 asl->P.khesoprod = kh;
490 h = h1 = (Hesoprod *)new_mblk(kh);
491 h2 = h + (sizeof(Char*) << kh)/sizeof(Hesoprod) - 1;
492 while(h1 < h2)
493 h1 = h1->next = h1 + 1;
494 h1->next = 0;
495 }
496 asl->P.hop_free = h->next;
497 FREE_DTOA_LOCK(HESOPROD_LOCK);
498 h->left = L;
499 h->right = R;
500 h->coef = coef;
501 hp = asl->P.otodo + R->varno;
502 h->next = *hp;
503 *hp = h;
504 }
505
506 static void
507 #ifdef KR_headers
del_Hesoprod(asl,x)508 del_Hesoprod(asl, x) ASL_pfgh *asl; Hesoprod *x;
509 #else
510 del_Hesoprod(ASL_pfgh *asl, Hesoprod *x)
511 #endif
512 {
513 x->next = asl->P.hop_free;
514 asl->P.hop_free = x;
515 }
516
517 static real *
518 #ifdef KR_headers
saveog(asl,no,noe,y,kp)519 saveog(asl, no, noe, y, kp) ASL_pfgh *asl; int no; int noe; int y; int *kp;
520 #else
521 saveog(ASL_pfgh *asl, int no, int noe, int y, int *kp)
522 #endif
523 {
524 real *o, *ogsave;
525 int i, k, n;
526 ps_func *p, *pe;
527 psg_elem *g, *ge;
528 ograd *og;
529
530 n = 0;
531 if (asl->P.nobjgroups)
532 for(i = no; i < noe; i++) {
533 p = asl->P.ops + i;
534 g = p->g;
535 for(ge = g + p->ng; g < ge; g++)
536 for(og = g->og; og; og = og->next)
537 n++;
538 }
539 if (asl->P.ncongroups && y) {
540 p = asl->P.cps;
541 for(pe = p + n_con; p < pe; p++)
542 for(g = p->g, ge = g + p->ng; g < ge; g++)
543 for(og = g->og; og; og = og->next)
544 n++;
545 }
546 if (!n)
547 return 0;
548 k = *kp = htcl(n*sizeof(real));
549 o = ogsave = (real*)new_mblk(k);
550 if (asl->P.nobjgroups)
551 for(i = no; i < noe; i++) {
552 p = asl->P.ops + i;
553 g = p->g;
554 for(ge = g + p->ng; g < ge; g++)
555 for(og = g->og; og; og = og->next)
556 *o++ = og->coef;
557 }
558 if (asl->P.ncongroups && y) {
559 p = asl->P.cps;
560 for(pe = p + n_con; p < pe; p++)
561 for(g = p->g, ge = g + p->ng; g < ge; g++)
562 for(og = g->og; og; og = og->next)
563 *o++ = og->coef;
564 }
565 return ogsave;
566 }
567
568 static void
569 #ifdef KR_headers
restog(asl,ogsave,no,noe,y,k)570 restog(asl, ogsave, no, noe, y, k) ASL_pfgh *asl; real *ogsave; int no; int noe; int y; int k;
571 #else
572 restog(ASL_pfgh *asl, real *ogsave, int no, int noe, int y, int k)
573 #endif
574 {
575 real *o = ogsave;
576 int i;
577 ps_func *p, *pe;
578 psg_elem *g, *ge;
579 ograd *og;
580
581 if (asl->P.nobjgroups)
582 for(i = no; i < noe; i++) {
583 p = asl->P.ops + i;
584 g = p->g;
585 for(ge = g + p->ng; g < ge; g++)
586 for(og = g->og; og; og = og->next)
587 og->coef = *o++;
588 }
589 if (asl->P.ncongroups && y) {
590 p = asl->P.cps;
591 for(pe = p + n_con; p < pe; p++)
592 for(g = p->g, ge = g + p->ng; g < ge; g++)
593 for(og = g->og; og; og = og->next)
594 og->coef = *o++;
595 }
596 Del_mblk_ASL((ASL*)asl, k, ogsave);
597 }
598
599 static fint
600 #ifdef KR_headers
bothadj(asl,spi)601 bothadj(asl, spi) ASL_pfgh *asl; SputInfo *spi;
602 #else
603 bothadj(ASL_pfgh *asl, SputInfo *spi)
604 #endif
605 {
606 /* Adjust to compute both triangles of Hessian */
607 fint i, i0, i1, j, k, k0, L, n, n1, nod, nz;
608 int kz, *z, *z0, *z1;
609 fint *hcs, *hr, *hre, *hrn, *hrn0, *ucs, *ulc, *uli;
610
611 n = n_var;
612 if ((nod = spi->nod) >= 0) {
613 if (!nod)
614 return 0;
615 goto done;
616 }
617 n1 = n + 1;
618 hcs = spi->hcolstarts;
619 nod = nz = hcs[n] - hcs[0];
620 hr = spi->hrownos - 1;
621 i = i0 = Fortran;
622 for(j = i + n; i < j; i++, hcs++) {
623 hr += k = hcs[1] - hcs[0];
624 if (k && *hr == i)
625 --nod;
626 }
627 /* nod = number of off-diagonal elements in upper triangle */
628 if (!(spi->nod = nod))
629 return 0; /* diagonal Hessian */
630 nz += nod;
631 spi->khinfob = kz = htcl((nz+2*(nod+n1))*sizeof(fint));
632 spi->ulinc0 = uli = (fint*)new_mblk(kz);
633 spi->hcs[1] = hcs = uli + n1;
634 spi->hrn[1] = hrn0 = hcs + n1;
635 spi->ulcopy0 = ulc = hrn0 + nz;
636 z = z0 = (int*)new_mblk(kz = htcl(n*sizeof(int)));
637 z1 = z - Fortran;
638 ucs = spi->hcs[0];
639 hre = spi->hrn[0];
640 for(i = i0; i < j; i++, ucs++) {
641 hr = hre;
642 hre += *z++ = ucs[1] - ucs[0];
643 while(hr < hre)
644 if ((k = *hr++) != i)
645 z1[k]++;
646 }
647 ucs = spi->hcs[0];
648 hre = spi->hrn[0];
649 *uli++ = 0;
650 for(i = k = i0; i < j; i++, ucs++) {
651 hr = hre;
652 hre += L = ucs[1] - ucs[0];
653 *hcs++ = k;
654 k0 = k - i0;
655 hrn = hrn0 + k0;
656 *uli++ = z1[i] - L;
657 k += z1[i];
658 z1[i] = k0 + L;
659 while(hr < hre)
660 if ((i1 = *hrn++ = *hr++) != i) {
661 *ulc++ = k0++;
662 hrn0[*ulc++ = z1[i1]++] = i;
663 }
664 }
665 *hcs = k;
666 spi->ulcend = ulc;
667 Del_mblk_ASL((ASL*)asl, kz, z0);
668 spi->ulinc = spi->ulinc0;
669 spi->ulcopy = spi->ulcopy0;
670 done:
671 spi->hrownos = spi->hrn[1];
672 spi->hcolstarts = spi->hcs[1];
673 return nod;
674 }
675
676 static void
677 #ifdef KR_headers
upper_to_lower(asl,spi,nz)678 upper_to_lower(asl, spi, nz) ASL_pfgh *asl; SputInfo *spi; fint nz;
679 #else
680 upper_to_lower(ASL_pfgh *asl, SputInfo *spi, fint nz)
681 #endif
682 { /* convert upper to lower triangular */
683
684 fint *cs, *hcolstarts, *hrownos, *rn;
685 int f, i, j, j1, j2, k, n;
686 int *rs, *u0, *utoL, *z;
687
688 f = Fortran;
689 n = n_var;
690 hrownos = spi->hrownos;
691 hcolstarts = spi->hcolstarts;
692 k = htcl((nz + n + 1)*sizeof(fint));
693 rn = spi->hrownos = spi->ulinc0 = (fint*)new_mblk(k);
694 spi->khinfob = k;
695 spi->hcolstarts = cs = rn + nz;
696 k = htcl((n+nz)*sizeof(int));
697 rs = (int*)new_mblk(k);
698 z = rs + n;
699 memset(rs, 0, n*sizeof(fint));
700 for(i = 0; i < nz; i++)
701 rs[hrownos[i]-f]++;
702 for(i = j = 0; i < n; i++) {
703 cs[i] = j + f;
704 j1 = rs[i];
705 rs[i] = j;
706 j += j1;
707 }
708 cs[n] = nz + f;
709 j1 = hcolstarts[1] - f;
710 for(i = j = 0; i < nz; i++) {
711 while(i >= j1)
712 j1 = hcolstarts[++j + 1] - f;
713 rn[z[i] = rs[hrownos[i]-f]++] = j + f;
714 }
715 for(i = j = 0; i < nz; i++) {
716 if ((j1 = z[i]) <= i) {
717 if (j1 < 0)
718 z[i] = -(j1 + 1);
719 continue;
720 }
721 j += 3;
722 while((j2 = z[j1]) != i) {
723 z[j1] = -(j2 + 1);
724 j++;
725 j1 = j2;
726 }
727 }
728 if (j) {
729 j += 2;
730 j1 = htcl(j*sizeof(int));
731 spi->uptolow = utoL = (int*)new_mblk(j1);
732 *utoL++ = j1;
733 for(i = 0; i < nz; i++) {
734 if ((j = z[i]) <= i)
735 continue;
736 u0 = utoL++;
737 *utoL++ = i;
738 *utoL++ = j;
739 while((j2 = z[j]) != i) {
740 z[j] = -(j2 + 1);
741 j = *utoL++ = j2;
742 }
743 *u0 = (utoL - u0) - 1;
744 }
745 *utoL = 0;
746 }
747 Del_mblk_ASL((ASL*)asl, k, rs);
748 }
749
750 fint
751 #ifdef KR_headers
sphes_setup_ASL(a,pspi,nobj,ow,y,uptri)752 sphes_setup_ASL(a, pspi, nobj, ow, y, uptri) ASL *a; SputInfo **pspi;
753 int nobj, ow, y, uptri;
754 #else
755 sphes_setup_ASL(ASL *a, SputInfo **pspi, int nobj, int ow, int y, int uptri)
756 #endif
757 {
758 int i, j, k, khinfo, kog, kz, n, n1, nhinfo, no, noe, nqslim, nzc;
759 int rfilter, robjno;
760 int *ui, *zc, *zci;
761 linarg *la, **lap, **lap1, **lape;
762 expr_v *v;
763 range *r, *r0, **rp, **rtodo;
764 real *ogsave, *s, *si, t;
765 ograd *og, *og1, **ogp, **ogpe;
766 Hesoprod *hop, *hop1, **otodo, **otodoi, **otodoj;
767 uHeswork *uhw, *uhwi, **utodo, **utodoi, **utodoj;
768 fint *hcolstarts, *hr, *hre, *hrownos, rv, *tf;
769 derp *D1;
770 de *d;
771 psb_elem *b;
772 psg_elem *g, *ge;
773 ps_func *p, *pe;
774 ASL_pfgh *asl;
775 expr_va *valist;
776 expr_if *iflist;
777 SputInfo *spi, *spi1;
778
779 asl = pscheck_ASL(a, "sphes_setup");
780 if (!pspi)
781 pspi = &asl->i.sputinfo_;
782 if (nobj >= 0 && nobj < n_obj) {
783 robjno = -2 - nobj;
784 rfilter = n_obj > 1 || !y && n_con > 0;
785 ow = 0;
786 no = nobj;
787 noe = no + 1;
788 }
789 else {
790 robjno = nobj = -1;
791 rfilter = !ow && n_obj > 0 || !y && n_con > 0;
792 no = noe = 0;
793 if (ow) {
794 noe = n_obj;
795 ow = 1;
796 }
797 }
798 if (y)
799 y = 1;
800 n = n_var;
801 if (spi = *pspi) {
802 if (spi->ow == ow && spi->y == y && spi->nobj == nobj
803 && spi->uptri == uptri)
804 goto done;
805 del_mblk(spi->khinfo, spi);
806 if (spi->ulinc0)
807 del_mblk(spi->khinfob, spi->ulinc0);
808 if (ui = spi->uptolow)
809 del_mblk(*ui, ui);
810 *pspi = 0;
811 }
812 if (!asl->P.hes_setup_called)
813 (*asl->p.Hesset)(a, 1, 0, nlo, 0, nlc);
814 asl->P.hes_setup_called = 3;
815 asl->P.iflist = 0;
816 asl->P.valist = 0;
817 otodo = otodoi = asl->P.otodo;
818 rtodo = asl->P.rtodo;
819 utodo = utodoi = asl->P.utodo;
820 s = asl->P.dOscratch;
821 nqslim = n >> 3;
822 kz = htcl(2*sizeof(int)*n);
823 zc = (int*)new_mblk_ASL(a, kz);
824 zci = zc + n;
825 memset(zc, 0, n*sizeof(int));
826 n1 = n + 1;
827 khinfo = htcl((2*n + 30)*sizeof(fint) + sizeof(SputInfo));
828 spi = (SputInfo*)new_mblk_ASL(a, khinfo);
829 hcolstarts = (fint*)(spi+1);
830 hr = hrownos = hcolstarts + n1;
831 nhinfo = ((sizeof(Char*)<<khinfo) - sizeof(SputInfo)) / sizeof(fint);
832 hre = hr + (nhinfo - n1);
833 r0 = (range*)&asl->P.rlist;
834 for(r = asl->P.rlist.next; r != r0; r = r->rlist.next) {
835 if ((j = r->n) <= 0)
836 continue;
837 if (rfilter) {
838 for(b = r->refs; b; b = b->next) {
839 if (b->conno >= 0) {
840 if (y)
841 goto keep;
842 }
843 else if (b->conno == robjno)
844 goto keep;
845 }
846 continue;
847 }
848 keep:
849 i = r->lasttermno;
850 rp = rtodo + i;
851 r->hnext = *rp;
852 *rp = r;
853 }
854 ogsave = asl->P.npsgcomp ? saveog(asl, no, noe, y, &kog) : 0;
855 if (asl->P.nobjgroups)
856 for(i = no; i < noe; i++) {
857 p = asl->P.ops + i;
858 g = p->g;
859 for(ge = g + p->ng; g < ge; g++)
860 if (og = g->og) {
861 do og->coef = 1; while(og = og->next);
862 og = g->og;
863 new_Hesoprod(asl, og, og, 1.);
864 }
865 }
866 if (asl->P.ncongroups && y) {
867 p = asl->P.cps;
868 for(pe = p + n_con; p < pe; p++)
869 for(g = p->g, ge = g + p->ng; g < ge; g++)
870 if (og = g->og) {
871 do og->coef = 1; while(og = og->next);
872 og = g->og;
873 new_Hesoprod(asl, og, og, 1.);
874 }
875 }
876 for(i = 0; i < n; i++) {
877 nzc = 0;
878 rp = rtodo;
879 uhwi = *utodoi;
880 *utodoi++ = 0;
881 while(r = *rp) {
882 rp = &r->hnext;
883 lap = r->lap;
884 lape = lap + r->n;
885 if (r->n >= r->nv) {
886 k = htcl(sizeof(uHeswork)
887 + (r->n - 1)*sizeof(ograd*));
888 uhw = (uHeswork *)new_mblk_ASL(a, k);
889 uhw->k = k;
890 uhw->next = uhwi;
891 uhwi = uhw;
892 uhw->r = r;
893 uhw->ui = ui = r->ui;
894 uhw->uie = ui + r->nv;
895 ogp = uhw->ogp;
896 while(lap < lape)
897 *ogp++ = (*lap++)->nz;
898 }
899 else {
900 si = s;
901 while(lap < lape) {
902 *si = 1;
903 pshv_prod1(asl, r, nobj, ow, y);
904 *si++ = 0;
905 lap1 = lap++;
906 la = *lap1++;
907 og = la->nz;
908 v = la->v;
909 if (t = v->aO)
910 new_Hesoprod(asl,og,og,t);
911 while(lap1 < lape) {
912 la = *lap1++;
913 v = la->v;
914 if (t = v->aO) {
915 og1 = la->nz;
916 new_Hesoprod(asl,og,og1,t);
917 new_Hesoprod(asl,og1,og,t);
918 }
919 }
920 }
921 }
922 }
923 *rtodo++ = 0; /* reset */
924 while(uhw = uhwi) {
925 uhwi = uhwi->next;
926 si = s;
927 ogp = uhw->ogp;
928 r = uhw->r;
929 ogpe = ogp + r->n;
930 si = s;
931 do {
932 if ((og = *ogp) && og->varno == i)
933 *si = 1.; /* og->coef til 20080629 */
934 si++;
935 } while(++ogp < ogpe);
936 pshv_prod1(asl, r, nobj, ow, y);
937
938 lap = r->lap;
939 lape = lap + r->n;
940 do {
941 la = *lap++;
942 if (la->v->aO)
943 for(og = la->nz; og; og = og->next)
944 if ((j = og->varno) <= i
945 && !zc[j]++)
946 zci[nzc++] = j;
947 }
948 while(lap < lape);
949
950 ogp = uhw->ogp;
951 si = s;
952 do {
953 if ((og = *ogp) && og->varno == i) {
954 *si = 0;
955 *ogp = og->next;
956 }
957 si++;
958 } while(++ogp < ogpe);
959 if ((ui = ++uhw->ui) >= uhw->uie)
960 del_mblk(uhw->k, uhw);
961 else {
962 utodoj = utodo + *ui;
963 uhw->next = *utodoj;
964 *utodoj = uhw;
965 }
966 }
967
968 hop1 = *otodoi;
969 *otodoi++ = 0;
970 while(hop = hop1) {
971 hop1 = hop->next;
972 og = hop->left;
973 og1 = hop->right;
974 while((j = og->varno) <= i) {
975 if (!zc[j]++)
976 zci[nzc++] = j;
977 if (!(og = og->next))
978 break;
979 }
980 if (og = og1->next) {
981 hop->right = og;
982 otodoj = otodo + og->varno;
983 hop->next = *otodoj;
984 *otodoj = hop;
985 }
986 else
987 del_Hesoprod(asl,hop);
988 }
989 hcolstarts[i] = hr - hrownos;
990 if (nzc > hre - hr) {
991 k = khinfo++;
992 spi1 = (SputInfo*)new_mblk_ASL(a, khinfo);
993 tf = (fint*)(spi1+1);
994 memcpy(tf, hcolstarts, (hr - hcolstarts)*sizeof(fint));
995 del_mblk(k, spi);
996 spi = spi1;
997 hcolstarts = tf;
998 hrownos = tf + n1;
999 hr = hrownos + hcolstarts[i];
1000 nhinfo = ((sizeof(Char*)<<khinfo) - sizeof(SputInfo))
1001 / sizeof(fint);
1002 hre = hrownos + (nhinfo - n1);
1003 }
1004 if (nzc > nqslim) {
1005 for(j = 0; j < n; j++)
1006 if (zc[j])
1007 zc[*hr++ = j] = 0;
1008 }
1009 else {
1010 if (nzc > 1)
1011 qsort(zci, nzc, sizeof(int), compar);
1012 for(j = 0; j < nzc; j++)
1013 zc[*hr++ = zci[j]] = 0;
1014 }
1015 }
1016 for(valist = asl->P.valist; valist; valist = valist->next) {
1017 D1 = valist->d0;
1018 for(d = valist->L.d; d->e; d++)
1019 d->dlast->next = D1;
1020 }
1021 for(iflist = asl->P.iflist; iflist; iflist = iflist->next)
1022 iflist->dTlast->next = iflist->d0;
1023 hcolstarts[n] = hr - hrownos;
1024 if (j = Fortran) {
1025 for(i = 0; i <= n; i++)
1026 hcolstarts[i] += j;
1027 i = (int)(hcolstarts[n] - j);
1028 while(i)
1029 hrownos[--i] += j;
1030 }
1031 spi->hcs[0] = hcolstarts;
1032 spi->hrn[0] = hrownos;
1033 spi->nod = -1;
1034 spi->ulcend = 0;
1035 spi->khinfo = khinfo;
1036 spi->nobj = nobj;
1037 spi->ow = ow;
1038 spi->y = y;
1039 spi->uptri = uptri;
1040 *pspi = spi;
1041 if (ogsave)
1042 restog(asl, ogsave, no, noe, y, kog);
1043 spi->ulinc0 = spi->ulinc = spi->ulcopy = 0;
1044 spi->uptolow = 0;
1045 del_mblk(kz, zc);
1046 done:
1047 spi->hrownos = spi->hrn[0];
1048 spi->hcolstarts = hcolstarts = spi->hcs[0];
1049 rv = hcolstarts[n] - hcolstarts[0];
1050 switch(uptri) {
1051 case 0:
1052 rv += bothadj(asl, spi);
1053 break;
1054 case 2:
1055 upper_to_lower(asl, spi, rv);
1056 }
1057 return rv;
1058 }
1059
1060 void
1061 #ifdef KR_headers
sphes_ASL(a,pspi,H,nobj,ow,y)1062 sphes_ASL(a, pspi, H, nobj, ow, y) ASL *a; SputInfo **pspi; real *H; int nobj; real *ow, *y;
1063 #else
1064 sphes_ASL(ASL *a, SputInfo **pspi, real *H, int nobj, real *ow, real *y)
1065 #endif
1066 {
1067 /* sparse upper triangle of Hessian */
1068
1069 int i, j, k, kh, n, no, noe, *ui;
1070 linarg *la, **lap, **lap1, **lape;
1071 expr_v *v;
1072 range *r, *r0, **rp, **rtodo;
1073 real *Hi, *H0, *H00;
1074 real *cscale, *owi, *s, *si, t, t1, *vsc0, *vsc1, *vsc, *y1;
1075 ograd *og, *og1, **ogp, **ogpe;
1076 Hesoprod *hop, *hop1, **otodo, **otodoi, **otodoj;
1077 uHeswork *uhw, *uhwi, **utodo, **utodoi, **utodoj;
1078 fint *hcs, *hr, *uli;
1079 psg_elem *g, *ge;
1080 ps_func *p, *pe;
1081 ASL_pfgh *asl;
1082 SputInfo *spi;
1083
1084 asl = pscheck_ASL(a, "sputhes");
1085 xpsg_check_ASL(asl, nobj, ow, y);
1086 if (!pspi)
1087 pspi = &a->i.sputinfo_;
1088 i = j = 0;
1089 if (y)
1090 j = 1;
1091 if (nobj >= 0 && nobj < n_obj) {
1092 no = nobj;
1093 noe = no + 1;
1094 owi = ow ? ow + no : &edag_one_ASL;
1095 ow = 0;
1096 }
1097 else {
1098 nobj = -1;
1099 no = noe = 0;
1100 if (owi = ow) {
1101 noe = n_obj;
1102 i = 1;
1103 }
1104 }
1105 if (asl->P.hes_setup_called != 3)
1106 sphes_setup_ASL(a, pspi, nobj, ow != 0, y != 0, 0);
1107 spi = *pspi;
1108 if (spi->nobj != nobj || spi->ow != i || spi->y != j) {
1109 fprintf(Stderr,
1110 "\nsphes() call inconsistent with previous sphsetup()\n");
1111 exit(1);
1112 }
1113 otodo = otodoi = asl->P.otodo;
1114 rtodo = asl->P.rtodo;
1115 utodo = utodoi = asl->P.utodo;
1116 s = asl->P.dOscratch;
1117 n = n_var;
1118 Hi = H0 = (real*)new_mblk_ASL(a, kh = htcl(n*sizeof(real)));
1119 memset(Hi, 0, n * sizeof(real));
1120 H0 -= Fortran;
1121 r0 = (range*)&asl->P.rlist;
1122 for(r = asl->P.rlist.next; r != r0; r = r->rlist.next) {
1123 if ((j = r->n) <= 0)
1124 continue;
1125 i = r->lasttermno;
1126 rp = rtodo + i;
1127 r->hnext = *rp;
1128 *rp = r;
1129 }
1130 if (asl->P.nobjgroups)
1131 for(; no < noe; no++)
1132 if (t = *owi++) {
1133 p = asl->P.ops + no;
1134 g = p->g;
1135 for(ge = g + p->ng; g < ge; g++)
1136 if (t1 = t*g->g2)
1137 for(og = g->og; og; og = og->next)
1138 if (og->coef) {
1139 new_Hesoprod(asl, og, og, t1);
1140 break;
1141 }
1142 }
1143 if (asl->P.ncongroups && y) {
1144 cscale = asl->i.lscale;
1145 p = asl->P.cps;
1146 y1 = y;
1147 for(pe = p + n_con; p < pe; p++, y1++)
1148 if (t = cscale ? *cscale++ * *y1 : *y1)
1149 for(g = p->g, ge = g + p->ng; g < ge; g++)
1150 if (t1 = t*g->g2)
1151 for(og = g->og; og; og = og->next)
1152 if (og->coef) {
1153 new_Hesoprod(asl, og, og, t1);
1154 break;
1155 }
1156 }
1157 hcs = spi->hcs[0];
1158 hr = spi->hrn[0];
1159 uli = spi->ulinc;
1160 H00 = H;
1161 if (vsc = asl->i.vscale) {
1162 vsc0 = vsc - Fortran;
1163 vsc1 = vsc;
1164 }
1165 for(i = 0; i < n; i++) {
1166 rp = rtodo;
1167 uhwi = *utodoi;
1168 *utodoi++ = 0;
1169 while(r = *rp) {
1170 rp = &r->hnext;
1171 lap = r->lap;
1172 lape = lap + r->n;
1173 if (r->n >= r->nv) {
1174 k = htcl(sizeof(uHeswork)
1175 + (r->n - 1)*sizeof(ograd*));
1176 uhw = (uHeswork *)new_mblk_ASL(a, k);
1177 uhw->k = k;
1178 uhw->next = uhwi;
1179 uhwi = uhw;
1180 uhw->r = r;
1181 uhw->ui = ui = r->ui;
1182 uhw->uie = ui + r->nv;
1183 ogp = uhw->ogp;
1184 while(lap < lape)
1185 *ogp++ = (*lap++)->nz;
1186 }
1187 else {
1188 si = s;
1189 while(lap < lape) {
1190 *si = 1;
1191 pshv_prod_ASL(asl, r, nobj, ow, y);
1192 *si++ = 0;
1193 lap1 = lap++;
1194 la = *lap1++;
1195 og = la->nz;
1196 v = la->v;
1197 if (t = v->aO)
1198 new_Hesoprod(asl,og,og,t);
1199 while(lap1 < lape) {
1200 la = *lap1++;
1201 v = la->v;
1202 if (t = v->aO) {
1203 og1 = la->nz;
1204 new_Hesoprod(asl,og,og1,t);
1205 new_Hesoprod(asl,og1,og,t);
1206 }
1207 }
1208 }
1209 }
1210 }
1211 *rtodo++ = 0; /* reset */
1212 while(uhw = uhwi) {
1213 uhwi = uhwi->next;
1214 si = s;
1215 ogp = uhw->ogp;
1216 r = uhw->r;
1217 ogpe = ogp + r->n;
1218 si = s;
1219 do {
1220 if ((og = *ogp) && og->varno == i)
1221 *si = og->coef;
1222 si++;
1223 } while(++ogp < ogpe);
1224 pshv_prod_ASL(asl, r, nobj, ow, y);
1225
1226 lap = r->lap;
1227 lape = lap + r->n;
1228 do {
1229 la = *lap++;
1230 if (t = la->v->aO)
1231 for(og = la->nz; og; og = og->next)
1232 if ((j = og->varno) <= i)
1233 Hi[j] += t*og->coef;
1234 }
1235 while(lap < lape);
1236
1237 ogp = uhw->ogp;
1238 si = s;
1239 do {
1240 if ((og = *ogp) && og->varno == i) {
1241 *si = 0;
1242 *ogp = og->next;
1243 }
1244 si++;
1245 } while(++ogp < ogpe);
1246 if ((ui = ++uhw->ui) >= uhw->uie)
1247 del_mblk(uhw->k, uhw);
1248 else {
1249 utodoj = utodo + *ui;
1250 uhw->next = *utodoj;
1251 *utodoj = uhw;
1252 }
1253 }
1254
1255 hop1 = *otodoi;
1256 *otodoi++ = 0;
1257 while(hop = hop1) {
1258 hop1 = hop->next;
1259 og = hop->left;
1260 og1 = hop->right;
1261 t = hop->coef * og1->coef;
1262 while((j = og->varno) <= i) {
1263 Hi[j] += t*og->coef;
1264 if (!(og = og->next))
1265 break;
1266 }
1267 if (og = og1->next) {
1268 hop->right = og;
1269 otodoj = otodo + og->varno;
1270 hop->next = *otodoj;
1271 *otodoj = hop;
1272 }
1273 else
1274 del_Hesoprod(asl,hop);
1275 }
1276 k = (int)(hcs[1] - hcs[0]);
1277 hcs++;
1278 if (uli)
1279 H += *uli++;
1280 if (vsc) {
1281 t = *vsc1++;
1282 while(--k >= 0) {
1283 j = (int)*hr++;
1284 *H++ = t * vsc0[j] * H0[j];
1285 H0[j] = 0;
1286 }
1287 }
1288 else
1289 while(--k >= 0) {
1290 *H++ = H0[j = (int)*hr++];
1291 H0[j] = 0;
1292 }
1293 }
1294 del_mblk(kh, Hi);
1295 H = H00;
1296 if (hr = spi->ulcopy)
1297 for(uli = spi->ulcend; hr < uli; hr += 2)
1298 H[hr[1]] = H[hr[0]];
1299 else if (ui = spi->uptolow)
1300 while(k = *++ui) {
1301 t = H[j = *++ui];
1302 while(--k) {
1303 t1 = H[i = *++ui];
1304 H[i] = t;
1305 t = t1;
1306 }
1307 H[j] = t;
1308 }
1309 }
1310