1 /* rqfnb.f -- translated by f2c (version 20050501)
2 and slightly cleaned up by Allin Cottrell
3 */
4
5 #include "libgretl.h"
6 #include "gretl_f2c.h"
7
8 /* Table of constant values */
9
10 static double c_b4 = 1.;
11 static integer one = 1;
12 static double zero = 0.;
13 static double c_b13 = -1.;
14
15 /* lapack/blas functions called below */
16
17 extern int dsyr_ (char *, integer *, double *, double *,
18 integer *, double *, integer *);
19
20 extern int dposv_ (char *, integer *, integer *, double *, integer *,
21 double *, integer *, integer *);
22
23 extern int dgemv_ (char *, integer *, integer *, double *, double *,
24 integer *, double *, integer *, double *,
25 double *, integer *);
26
27 extern int dcopy_ (integer *, double *, integer *, double *,
28 integer *);
29
30 extern int dswap_ (integer *, double *, integer *, double *,
31 integer *);
32
33 extern int daxpy_ (integer *, double *, double *, integer *,
34 double *, integer *);
35
36 extern int dpotrs_ (char *, integer *, integer *, double *, integer *,
37 double *, integer *, integer *);
38
39 extern double ddot_ (integer *, double *, integer *, double *,
40 integer *);
41
stepy_(integer * n,integer * p,double * a,double * d,double * b,double * ada,integer * info)42 static int stepy_ (integer *n, integer *p, double *a,
43 double *d, double *b, double *ada,
44 integer *info)
45 {
46 integer i, m = *p * *p;
47 int attempt = 0;
48 int err = 0;
49
50 try_again:
51
52 for (i=0; i<m; i++) {
53 ada[i] = 0.0;
54 }
55
56 for (i=0; i<*n; i++) {
57 dsyr_("U", p, &d[i], &a[i * *p], &one, ada, p);
58 }
59
60 if (attempt == 0) {
61 dposv_("U", p, &one, ada, p, b, p, info);
62 if (*info != 0) {
63 fprintf(stderr, "stepy: dposv gave info = %d\n", *info);
64 attempt = 1;
65 goto try_again;
66 }
67 } else {
68 gretl_matrix A, B;
69
70 gretl_matrix_init(&A);
71 gretl_matrix_init(&B);
72
73 A.rows = A.cols = *p;
74 A.val = ada;
75 B.rows = *p;
76 B.cols = 1;
77 B.val = b;
78
79 err = gretl_LU_solve(&A, &B);
80 if (err) {
81 fprintf(stderr, "stepy: gretl_LU_solve: err = %d\n", err);
82 }
83 }
84
85 return err;
86 }
87
88 #define ITERSTEP 5
89
lpfnb_(integer * n,integer * p,double * a,double * c__,double * b,double * d__,double * u,double * beta,double * eps,double * x,double * s,double * y,double * z__,double * w,double * dx,double * ds,double * dy,double * dz,double * dw,double * dr,double * rhs,double * ada,integer * nit,integer * info,void (* callback)(void))90 static int lpfnb_ (integer *n, integer *p, double *a, double *c__,
91 double *b, double *d__, double *u, double *beta,
92 double *eps, double *x, double *s, double *y,
93 double *z__, double *w, double *dx, double *ds,
94 double *dy, double *dz, double *dw, double *dr,
95 double *rhs, double *ada, integer *nit, integer *info,
96 void (*callback)(void))
97 {
98 integer a_dim1 = *p, ada_dim1 = *p;
99 integer a_offset = 1 + a_dim1, ada_offset = 1 + ada_dim1;
100 double d1, d2;
101 static double g;
102 static integer i;
103 static double mu, gap;
104 static double dsdw, dxdz;
105 static double deltad, deltap;
106 int main_iters = 0;
107 int err = 0;
108
109 /* Parameter adjustments */
110 --dr;
111 --dw;
112 --dz;
113 --ds;
114 --dx;
115 --w;
116 --z__;
117 --s;
118 --x;
119 --u;
120 --d__;
121 --c__;
122 ada -= ada_offset;
123 --rhs;
124 --dy;
125 --y;
126 --b;
127 a -= a_offset;
128 --nit;
129
130 /* Function Body */
131 nit[1] = 0;
132 nit[2] = 0;
133 nit[3] = *n;
134 dgemv_("N", p, n, &c_b4, &a[a_offset], p, &c__[1], &one, &zero, &y[1],
135 &one);
136 for (i = 1; i <= *n; ++i) {
137 d__[i] = 1.;
138 }
139 err = stepy_(n, p, &a[a_offset], &d__[1], &y[1], &ada[ada_offset], info);
140 if (err) {
141 return err;
142 }
143 dcopy_(n, &c__[1], &one, &s[1], &one);
144 dgemv_("T", p, n, &c_b13, &a[a_offset], p, &y[1], &one, &c_b4, &s[1],
145 &one);
146 for (i = 1; i <= *n; ++i) {
147 if ((d1 = s[i], fabs(d1)) < *eps) {
148 d1 = s[i];
149 z__[i] = max(d1,0.) + *eps;
150 d1 = -s[i];
151 w[i] = max(d1,0.) + *eps;
152 } else {
153 d1 = s[i];
154 z__[i] = max(d1, 0.);
155 d1 = -s[i];
156 w[i] = max(d1, 0.);
157 }
158 s[i] = u[i] - x[i];
159 }
160
161 gap = ddot_(n, &z__[1], &one, &x[1], &one) +
162 ddot_(n, &w[1], &one, &s[1], &one);
163
164 looptop:
165
166 if (callback != NULL && (main_iters++ % ITERSTEP == 0)) {
167 callback();
168 }
169
170 if (gap > *eps && nit[1] < 50) {
171 ++nit[1];
172 for (i = 1; i <= *n; ++i) {
173 d__[i] = 1. / (z__[i] / x[i] + w[i] / s[i]);
174 ds[i] = z__[i] - w[i];
175 dz[i] = d__[i] * ds[i];
176 }
177
178 dcopy_(p, &b[1], &one, &dy[1], &one);
179 dgemv_("N", p, n, &c_b13, &a[a_offset], p, &x[1], &one, &c_b4, &dy[1],
180 &one);
181 dgemv_("N", p, n, &c_b4, &a[a_offset], p, &dz[1], &one, &c_b4, &dy[1],
182 &one);
183 dcopy_(p, &dy[1], &one, &rhs[1], &one);
184 err = stepy_(n, p, &a[a_offset], &d__[1], &dy[1], &ada[ada_offset], info);
185 if (err) {
186 return err;
187 }
188
189 dgemv_("T", p, n, &c_b4, &a[a_offset], p, &dy[1], &one, &c_b13,
190 &ds[1], &one);
191 deltap = 1e20;
192 deltad = 1e20;
193
194 for (i = 1; i <= *n; ++i) {
195 dx[i] = d__[i] * ds[i];
196 ds[i] = -dx[i];
197 dz[i] = -z__[i] * (dx[i] / x[i] + 1.);
198 dw[i] = -w[i] * (ds[i] / s[i] + 1.);
199 if (dx[i] < 0.) {
200 d1 = deltap, d2 = -x[i] / dx[i];
201 deltap = min(d1,d2);
202 }
203 if (ds[i] < 0.) {
204 d1 = deltap, d2 = -s[i] / ds[i];
205 deltap = min(d1,d2);
206 }
207 if (dz[i] < 0.) {
208 d1 = deltad, d2 = -z__[i] / dz[i];
209 deltad = min(d1,d2);
210 }
211 if (dw[i] < 0.) {
212 d1 = deltad, d2 = -w[i] / dw[i];
213 deltad = min(d1,d2);
214 }
215 }
216 d1 = *beta * deltap;
217 deltap = min(d1,1.);
218 d1 = *beta * deltad;
219 deltad = min(d1,1.);
220 if (min(deltap,deltad) < 1.) {
221 ++nit[2];
222 mu = ddot_(n, &x[1], &one, &z__[1], &one) +
223 ddot_(n, &s[1], &one, &w[1], &one);
224 g = mu + deltap * ddot_(n, &dx[1], &one, &z__[1], &one) +
225 deltad * ddot_(n, &dz[1], &one, &x[1], &one) +
226 deltap * deltad * ddot_(n, &dz[1], &one, &dx[1], &one) +
227 deltap * ddot_(n, &ds[1], &one, &w[1], &one) +
228 deltad * ddot_(n, &dw[1], &one, &s[1], &one) +
229 deltap * deltad * ddot_(n, &ds[1], &one, &dw[1], &one);
230 d1 = g / mu;
231 mu = mu * (d1 * (d1 * d1)) / (double) (*n << 1);
232 for (i = 1; i <= *n; ++i) {
233 dr[i] = d__[i] * (mu * (1 / s[i] - 1 / x[i]) + dx[i]
234 * dz[i] / x[i] - ds[i] * dw[i] / s[i]);
235 }
236 dswap_(p, &rhs[1], &one, &dy[1], &one);
237 dgemv_("N", p, n, &c_b4, &a[a_offset], p, &dr[1], &one, &c_b4,
238 &dy[1], &one);
239 dpotrs_("U", p, &one, &ada[ada_offset], p, &dy[1], p, info);
240 if (*info != 0) {
241 fprintf(stderr, "lpfnb: dpotrs_ gave info = %d\n", *info);
242 }
243 dgemv_("T", p, n, &c_b4, &a[a_offset], p, &dy[1], &one, &zero,
244 &u[1], &one);
245 deltap = 1e20;
246 deltad = 1e20;
247 for (i = 1; i <= *n; ++i) {
248 dxdz = dx[i] * dz[i];
249 dsdw = ds[i] * dw[i];
250 dx[i] = d__[i] * (u[i] - z__[i] + w[i]) - dr[i];
251 ds[i] = -dx[i];
252 dz[i] = -z__[i] + (mu - z__[i] * dx[i] - dxdz) / x[i];
253 dw[i] = -w[i] + (mu - w[i] * ds[i] - dsdw) / s[i];
254 if (dx[i] < 0.) {
255 d1 = deltap, d2 = -x[i] / dx[i];
256 deltap = min(d1, d2);
257 }
258 if (ds[i] < 0.) {
259 d1 = deltap, d2 = -s[i] / ds[i];
260 deltap = min(d1, d2);
261 }
262 if (dz[i] < 0.) {
263 d1 = deltad, d2 = -z__[i] / dz[i];
264 deltad = min(d1, d2);
265 }
266 if (dw[i] < 0.) {
267 d1 = deltad, d2 = -w[i] / dw[i];
268 deltad = min(d1, d2);
269 }
270 }
271 d1 = *beta * deltap;
272 deltap = min(d1,1.);
273 d1 = *beta * deltad;
274 deltad = min(d1,1.);
275 }
276
277 daxpy_(n, &deltap, &dx[1], &one, &x[1], &one);
278 daxpy_(n, &deltap, &ds[1], &one, &s[1], &one);
279 daxpy_(p, &deltad, &dy[1], &one, &y[1], &one);
280 daxpy_(n, &deltad, &dz[1], &one, &z__[1], &one);
281 daxpy_(n, &deltad, &dw[1], &one, &w[1], &one);
282 gap = ddot_(n, &z__[1], &one, &x[1], &one) +
283 ddot_(n, &w[1], &one, &s[1], &one);
284
285 goto looptop;
286 }
287
288 daxpy_(n, &c_b13, &w[1], &one, &z__[1], &one);
289 dswap_(n, &z__[1], &one, &x[1], &one);
290
291 return err;
292 } /* end of lpfnb_ */
293
rqfnb_(integer * n,integer * p,double * a,double * y,double * rhs,double * d,double * u,double * beta,double * eps,double * wn,double * wp,integer * nit,integer * info,void (* callback)(void))294 int rqfnb_ (integer *n, integer *p, double *a, double *y,
295 double *rhs, double *d, double *u, double *beta,
296 double *eps, double *wn, double *wp, integer *nit,
297 integer *info, void (*callback)(void))
298 {
299 integer a_dim = *p, wn_dim = *n, wp_dim = *p;
300 int err;
301
302 /* Parameter adjustments */
303 wn -= 1 + wn_dim;
304 wp -= 1 + wp_dim;
305 a -= 1 + a_dim;
306
307 err = lpfnb_(n, p, &a[a_dim + 1], y, rhs, d, u, beta, eps,
308 &wn[wn_dim + 1], &wn[(wn_dim << 1) + 1], &wp[wp_dim + 1],
309 &wn[wn_dim * 3 + 1], &wn[(wn_dim << 2) + 1], &wn[wn_dim * 5 + 1],
310 &wn[wn_dim * 6 + 1], &wp[(wp_dim << 1) + 1], &wn[wn_dim * 7 + 1],
311 &wn[(wn_dim << 3) + 1], &wn[wn_dim * 9 + 1], &wp[wp_dim * 3 + 1],
312 &wp[(wp_dim << 2) + 1], nit, info, callback);
313
314 return err;
315 }
316