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