1 /*******************************************************************
2 Copyright (C) 2017, 2018 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 "nlp.h"
20
21 #ifdef __cplusplus
22 extern "C" {
23 extern real con1ival_nomap_ASL(ASL *a, int i, real *X, fint *nerror);
24 extern void con1grd_nomap_ASL(ASL *a, int i, real *X, real *G, fint *nerror);
25 }
26 #endif
27
28 static void
INchk(ASL * asl,const char * who,int i,int ix)29 INchk(ASL *asl, const char *who, int i, int ix)
30 {
31 ASL_CHECK(asl, ASL_read_fg, who);
32 if (i < 0 || i >= ix) {
33 fprintf(Stderr, "%s: got I = %d; expected 0 <= I < %d\n",
34 who, i, ix);
35 exit(1);
36 }
37 }
38
39 static real
cival(ASL_fg * asl,int i,real * X,fint * nerror)40 cival(ASL_fg *asl, int i, real *X, fint *nerror)
41 {
42 Jmp_buf err_jmp0;
43 expr *e;
44 int ij, nc;
45 real f;
46
47 if (nerror && *nerror >= 0) {
48 err_jmp = &err_jmp0;
49 ij = setjmp(err_jmp0.jb);
50 if ((*nerror = ij))
51 return 0.;
52 }
53 want_deriv = want_derivs;
54 errno = 0; /* in case f77 set errno opening files */
55 co_index = i;
56 if (!asl->i.x_known)
57 x0_check_ASL(asl,X);
58 if (!asl->i.ncxval)
59 asl->i.ncxval = (int*)M1zapalloc(nclcon*sizeof(int));
60 if (!(x0kind & ASL_have_concom)) {
61 if (comb < combc)
62 comeval_ASL(asl, comb, combc);
63 if (comc1)
64 com1eval_ASL(asl, 0,comc1);
65 x0kind |= ASL_have_concom;
66 }
67 if (i >= (nc = asl->i.n_con0))
68 e = lcon_de[i-nc].e;
69 else
70 e = con_de[i].e;
71 f = (*e->op)(e C_ASL);
72 asl->i.ncxval[i] = asl->i.nxval;
73 err_jmp = 0;
74 return f;
75 }
76
77 static real
Conival1(ASL_fg * asl,int i,real * X,fint * nerror)78 Conival1(ASL_fg *asl, int i, real *X, fint *nerror)
79 {
80 cgrad *gr;
81 int j1, kv, *vmi;
82 real f, *vscale;
83
84 if (i < asl->i.n_con0)
85 f = cival(asl, i, X, nerror);
86 else
87 f = 0.;
88 kv = 0;
89 vmi = 0;
90 if ((vscale = asl->i.vscale))
91 kv = 2;
92 if (asl->i.vmap) {
93 vmi = get_vminv_ASL((ASL*)asl);
94 ++kv;
95 }
96 gr = asl->i.Cgrad0[i];
97 switch(kv) {
98 case 3:
99 for(; gr; gr = gr->next) {
100 j1 = vmi[gr->varno];
101 f += X[j1] * vscale[j1] * gr->coef;
102 }
103 break;
104 case 2:
105 for(; gr; gr = gr->next) {
106 j1 = gr->varno;
107 f += X[j1] * vscale[j1] * gr->coef;
108 }
109 break;
110 case 1:
111 for(; gr; gr = gr->next)
112 f += X[vmi[gr->varno]] * gr->coef;
113 break;
114 case 0:
115 for(; gr; gr = gr->next)
116 f += X[gr->varno] * gr->coef;
117 }
118 return f;
119 }
120
121 real
con1ival_nomap_ASL(ASL * a,int i,real * X,fint * nerror)122 con1ival_nomap_ASL(ASL *a, int i, real *X, fint *nerror)
123 {
124 INchk(a, "con1ival_nomap", i, a->i.n_con0);
125 return Conival1((ASL_fg*)a, i, X, nerror);
126 }
127
128 real
con1ival(ASL * a,int i,real * X,fint * nerror)129 con1ival(ASL *a, int i, real *X, fint *nerror)
130 {
131 int *cm;
132
133 INchk(a, "con1ival", i, a->i.n_con_);
134 if ((cm = a->i.cmap))
135 i = cm[i];
136 return Conival1((ASL_fg*)a, i, X, nerror);
137 }
138
139 int
lcon1val(ASL * a,int i,real * X,fint * nerror)140 lcon1val(ASL *a, int i, real *X, fint *nerror)
141 {
142 real f;
143
144 INchk(a, "lcon1ival", i, a->i.n_lcon_);
145 f = cival((ASL_fg*)a, i + a->i.n_con0, X, nerror);
146 return f != 0.;
147 }
148
149 static void
Congrd1(ASL_fg * asl,int i,real * X,real * G,fint * nerror)150 Congrd1(ASL_fg *asl, int i, real *X, real *G, fint *nerror)
151 {
152 Jmp_buf err_jmp0;
153 cde *d;
154 cgrad *gr, *gr1;
155 int i0, ij, j, *vmi, xksave;
156 real *Adjoints, *vscale;
157 size_t L;
158
159 if (nerror && *nerror >= 0) {
160 err_jmp = &err_jmp0;
161 ij = setjmp(err_jmp0.jb);
162 if ((*nerror = ij))
163 return;
164 }
165 errno = 0; /* in case f77 set errno opening files */
166 if (!asl->i.x_known) {
167 co_index = i;
168 x0_check_ASL(asl,X);
169 }
170 if ((!asl->i.ncxval || asl->i.ncxval[i] != asl->i.nxval)
171 && (!(x0kind & ASL_have_conval)
172 || i < n_conjac[0] || i >= n_conjac[1])) {
173 xksave = asl->i.x_known;
174 asl->i.x_known = 1;
175 con1ival_nomap_ASL((ASL*)asl,i,X,nerror);
176 asl->i.x_known = xksave;
177 if (nerror && *nerror)
178 return;
179 }
180 if (asl->i.Derrs)
181 deriv_errchk_ASL((ASL*)asl, nerror, i, 1);
182 if (!(x0kind & ASL_have_funnel)) {
183 if (f_b)
184 funnelset_ASL(asl, f_b);
185 if (f_c)
186 funnelset_ASL(asl, f_c);
187 x0kind |= ASL_have_funnel;
188 }
189 Adjoints = adjoints;
190 d = &con_de[i];
191 gr1 = asl->i.Cgrad0[i];
192 for(gr = gr1; gr; gr = gr->next)
193 Adjoints[gr->varno] = gr->coef;
194 if ((L = d->zaplen)) {
195 memset(adjoints_nv1, 0, L);
196 derprop(d->d);
197 }
198 vmi = 0;
199 if (asl->i.vmap)
200 vmi = get_vminv_ASL((ASL*)asl);
201 if ((vscale = asl->i.vscale)) {
202 if (vmi)
203 for(gr = gr1; gr; gr = gr->next) {
204 i0 = gr->varno;
205 Adjoints[i0] *= vscale[vmi[i0]];
206 }
207 else
208 for(gr = gr1; gr; gr = gr->next) {
209 i0 = gr->varno;
210 Adjoints[i0] *= vscale[i0];
211 }
212 }
213 gr = gr1;
214 i0 = 0;
215 switch(asl->i.congrd_mode) {
216 case 1:
217 for(; gr; gr = gr->next)
218 G[i0++] = Adjoints[gr->varno];
219 break;
220 case 2:
221 for(; gr; gr = gr->next)
222 G[gr->goff] = Adjoints[gr->varno];
223 break;
224 default:
225 if (vmi) {
226 for(; gr; gr = gr->next) {
227 i = vmi[j = gr->varno];
228 while(i0 < i)
229 G[i0++] = 0;
230 G[i] = Adjoints[j];
231 i0 = i + 1;
232 }
233 }
234 else
235 for(; gr; gr = gr->next) {
236 i = gr->varno;
237 while(i0 < i)
238 G[i0++] = 0;
239 G[i] = Adjoints[i];
240 i0 = i + 1;
241 }
242 i = n_var;
243 while(i0 < i)
244 G[i0++] = 0;
245 }
246 err_jmp = 0;
247 }
248
249 void
con1grd_nomap_ASL(ASL * a,int i,real * X,real * G,fint * nerror)250 con1grd_nomap_ASL(ASL *a, int i, real *X, real *G, fint *nerror)
251 {
252 ASL_fg *asl;
253 static char who[] = "con1grd_nomap";
254
255 INchk(a, who, i, a->i.n_con0);
256 asl = (ASL_fg*)a;
257 if (!want_derivs)
258 No_derivs_ASL(who);
259 Congrd1(asl, i, X, G, nerror);
260 }
261
262 void
con1grd_ASL(ASL * a,int i,real * X,real * G,fint * nerror)263 con1grd_ASL(ASL *a, int i, real *X, real *G, fint *nerror)
264 {
265 ASL_fg *asl;
266 int *cm;
267 static char who[] = "con1grd";
268
269 INchk(a, who, i, a->i.n_con_);
270 asl = (ASL_fg*)a;
271 if (!want_derivs)
272 No_derivs_ASL(who);
273 if ((cm = a->i.cmap))
274 i = cm[i];
275 Congrd1(asl, i, X, G, nerror);
276 }
277