blob 05f5028b9ce25f3d9d9711c8eab0a77711c925ae
1 /* minimize - minimization for Bayes routines in XLISP-STAT and S */
2 /* XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney */
3 /* Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz */
4 /* You may give out copies of this software; for conditions see the */
5 /* file COPYING included with this distribution. */
7 /*
8 * Nonlinear optimization modules adapted from Dennis and Schnabel,
9 * "Numerical Methods for Unconstrained Optimization and Nonlinear
10 * Equations."
13 #include <stdio.h>
14 #include "xmath.h"
15 extern char buf[200];
16 #define PRINTSTR(s) printf(s)
18 typedef struct {
19 size_t n, k;
20 void (*ffun)(), (*gfun)();
21 double f, typf, new_f;
22 double crit, new_crit;
23 double *x, *new_x, *sx, *delf, *new_delf, *qnstep, *F;
24 double **hessf, **H, **L, **new_delg;
26 int itncount, itnlimit, maxtaken, consecmax, retcode, termcode;
27 int use_line_search, verbose, values_supplied, change_sign;
29 } Iteration;
32 extern double macheps();
35 extern double choldecomp();
40 /*
41 This is a little more advanced section on functions, but is very useful. Take this for example:
42 int applyeqn(int F(int), int max, int min) {
43 int itmp;
45 itmp = F(int) + min;
46 itmp = itmp - max;
47 return itmp;
50 What does this function do if we call it with applyeqn(square(x),
51 y, z);? What happens is that the int F(int) is a reference to the
52 function that is passed in as a parameter. Thus inside applyeqn where
53 there is a call to F, it actually is a call to square! This is very
54 useful if we have one set function, but wish to vary the input
55 according to a particular function. So if we had a different function
56 called cube we could change how we call applyeqn by calling the
57 function by applyeqn(cube(x), y, z);. */
60 /************************************************************************/
61 /** **/
62 /** Definitions and Globals **/
63 /** **/
64 /************************************************************************/
66 # define nil 0L
67 # define FALSE 0
68 # define TRUE 1
71 # define CONSEC_MAX_LIMIT 5
72 # define ALPHA .0001
73 # define MAX_STEP_FACTOR 1000
74 # define GRADTOL_POWER 1.0 / 3.0
75 # define STEPTOL_POWER 2.0 / 3.0
76 # define ITNLIMIT 100
77 # define VERBOSE 0
78 # define USE_SEARCH TRUE
80 typedef double **RMatrix, *RVector;
83 static char *termcodes[] = {"not yet terminated",
85 "step size is less than step tolerance",
86 "no satisfactory step found in backtracking",
87 "iteration limit exceeded",
88 "maximum size step taken 5 iterations in a row"};
90 /************************************************************************/
91 /** **/
92 /** Utility Functions **/
93 /** **/
94 /************************************************************************/
96 static double Max(a, b)
97 double a, b;
99 return(a > b ? a : b);
102 static double Min(a, b)
103 double a, b;
105 return(a > b ? b : a);
108 /************************************************************************/
109 /** **/
110 /** Cholesky Solving Functions **/
111 /** **/
112 /************************************************************************/
115 /* solve Ly = b for y */
116 static void
117 lsolve(size_t n, double *b, double **L, double *y)
119 size_t i, j;
121 for (i = 0; i < n; i++) {
122 y[i] = b[i];
123 for (j = 0; j < i; j++) y[i] -= L[i][j] * y[j];
124 if (L[i][i] != 0) y[i] /= L[i][i];
128 /* solve (L^T)x = y for x */
129 static void
130 ltsolve(size_t n, double *y, double **L, double *x)
132 size_t i, j;
134 for (i = n - 1; i >= 0; i--) {
135 x[i] = y[i];
136 for (j = i + 1; j < n; j++) x[i] -= L[j][i] * x[j];
137 if (L[i][i] != 0) x[i] /= L[i][i];
141 /* solve (L L^T) s = -g for s */
142 static void
143 cholsolve(size_t n,
144 double *g, double **L,
145 double *s)
147 size_t i;
149 /* solve Ly = g */
150 lsolve(n, g, L, s);
152 /* solve L^Ts = y */
153 ltsolve(n, s, L, s);
155 for (i = 0; i < n; i++) s[i] = -s[i];
158 static void
159 modelhess(size_t n, double *sx, double **H, double **L, double *diagadd)
161 size_t i, j;
162 double sqrteps, maxdiag, mindiag, maxoff, maxoffl, maxposdiag, mu,
163 maxadd, maxev, minev, offrow, sdd;
165 /* scale H on both sides with sx */
166 for (i = 0; i < n; i++)
167 for (j = 0; j < n; j++) H[i][j] /= sx[i] * sx[j];
170 * find mu to add to diagonal so no diagonal elements are negative
171 * and largest diagonal dominates largest off diagonal element in H
173 sqrteps = sqrt(macheps());
174 for (maxdiag = H[0][0], i = 1; i < n; i++)
175 maxdiag = Max(maxdiag, H[i][i]);
176 for (mindiag = H[0][0], i = 1; i < n; i++)
177 mindiag = Min(mindiag, H[i][i]);
178 maxposdiag = Max(0.0, maxdiag);
180 if (mindiag <= sqrteps * maxposdiag) {
181 mu = 2 * (maxposdiag - mindiag) * sqrteps - mindiag;
182 maxdiag += mu;
184 else mu = 0.0;
186 if (n > 1) {
187 for (maxoff = fabs(H[0][1]), i = 0; i < n; i++)
188 for (j = i + 1; j < n; j++)
189 maxoff = Max(maxoff, fabs(H[i][j]));
191 else maxoff = 0.0;
193 if (maxoff * (1 + 2 * sqrteps) > maxdiag) {
194 mu += (maxoff - maxdiag) + 2 * sqrteps * maxoff;
195 maxdiag = maxoff * (1 + 2 * sqrteps);
198 if (maxdiag == 0.0) {
199 mu = 1;
200 maxdiag = 1;
203 if (mu > 0) for (i = 0; i < n; i++) H[i][i] += mu;
205 maxoffl = sqrt(Max(maxdiag, maxoff / n));
208 * compute the perturbed Cholesky decomposition
210 for (i = 0; i < n; i++)
211 for (j = 0; j < n; j++) L[i][j] = H[i][j];
215 * add something to diagonal, if needed, to make positive definite
216 * and recompute factorization
218 if (maxadd > 0) {
219 maxev = H[0][0];
220 minev = H[0][0];
221 for (i = 0; i < n; i++) {
222 for (offrow = 0.0, j = 0; j < n; j++)
223 if (i != j) offrow += fabs(H[i][j]);
224 maxev = Max(maxev, H[i][i] + offrow);
225 minev = Min(minev, H[i][i] - offrow);
227 sdd = (maxev - minev) * sqrteps - minev;
228 sdd = Max(sdd, 0.0);
230 for (i = 0; i < n; i++) H[i][i] += mu;
231 for (i = 0; i < n; i++)
232 for (j = 0; j < n; j++) L[i][j] = H[i][j];
238 /* unscale H and L */
239 for (i = 0; i < n; i++)
240 for (j = 0; j < n; j++) {
241 H[i][j] *= sx[i] * sx[j];
242 L[i][j] *= sx[i];
246 /************************************************************************/
247 /** **/
248 /** Stopping Criteria **/
249 /** **/
250 /************************************************************************/
252 static double
255 size_t n, i;
256 double size, term, crit, typf;
257 double *x, *delf, *sx;
259 n = iter->n + iter->k;
260 crit = iter->crit;
261 typf = iter->typf;
262 x = iter->x;
263 sx = iter->sx;
264 delf = (new) ? iter->new_delf : iter->delf;
266 for (i = 0, size = 0.0; i < n; i++) {
267 term = fabs(delf[i]) * Max(x[i], 1.0 / sx[i]) / Max(fabs(crit), typf);
268 size = Max(size, term);
270 return(size);
273 static double
274 incrsize(Iteration *iter)
276 size_t n, i;
277 double size, term;
278 double *x, *new_x, *sx;
280 new_x = iter->new_x;
281 n = iter->n + iter->k;
282 x = iter->x;
283 sx = iter->sx;
285 for (i = 0, size = 0.0; i < n; i++) {
286 term = fabs(x[i] - new_x[i]) / Max(fabs(x[i]), 1.0 / sx[i]);
287 size = Max(size, term);
289 return(size);
292 static int /* AJR:TESTME guess! */
293 stoptest0(Iteration *iter)
295 iter->consecmax = 0;
298 iter->termcode = 1;
299 else iter->termcode = 0;
301 return(iter->termcode);
304 static int
305 stoptest(Iteration *iter)
307 int termcode, retcode, itncount, itnlimit, maxtaken, consecmax;
310 retcode = iter->retcode;
312 steptol = iter->steptol;
313 itncount = iter->itncount;
314 itnlimit = iter->itnlimit;
315 maxtaken = iter->maxtaken;
316 consecmax = iter->consecmax;
318 termcode = 0;
319 if (retcode == 1) termcode = 3;
321 else if (incrsize(iter) <= steptol) termcode = 2;
322 else if (itncount >= itnlimit) termcode = 4;
323 else if (maxtaken) {
324 consecmax++;
325 if (consecmax >= CONSEC_MAX_LIMIT) termcode = 5;
327 else consecmax = 0;
329 iter->consecmax = consecmax;
330 iter->termcode = termcode;
332 return(termcode);
335 /************************************************************************/
336 /** **/
337 /** Function and Derivative Evaluation **/
338 /** **/
339 /************************************************************************/
341 static void
342 eval_funval(Iteration *iter)
344 int i;
346 (*(iter->ffun))(iter->x, &iter->f, nil, nil);
347 if (iter->k == 0) iter->crit = iter->f;
348 else {
350 for (i = 0, iter->crit = 0.0; i < iter->n + iter->k; i++)
351 iter->crit += 0.5 * pow(fabs(iter->delf[i]), 2.0);
355 static void
356 eval_next_funval(Iteration *iter)
358 int i;
360 (*(iter->ffun))(iter->new_x, &iter->new_f, nil, nil);
361 if (iter->k == 0) iter->new_crit = iter->new_f;
362 else {
364 for (i = 0, iter->new_crit = 0.0; i < iter->n + iter->k; i++)
365 iter->new_crit += 0.5 * pow(fabs(iter->new_delf[i]), 2.0);
369 static void
372 size_t i, j, n, k;
374 n = iter->n;
375 k = iter->k;
377 (*(iter->ffun))(iter->x, nil, iter->delf, nil);
378 if (k > 0) {
379 (*(iter->gfun))(iter->x, iter->delf + n, nil, nil);
380 (*(iter->gfun))(iter->x, nil, iter->new_delg, nil);
381 for (i = 0; i < n; i++) {
382 for (j = 0; j < k; j++)
383 iter->delf[i] += iter->x[n + j] * iter->new_delg[j][i];
384 for (j = 0; j < k; j++) {
385 iter->hessf[n + j][i] = iter->new_delg[j][i];
386 iter->hessf[i][n + j] = iter->new_delg[j][i];
392 static void
395 size_t i, j, n, k;
397 n = iter->n;
398 k = iter->k;
399 (*(iter->ffun))(iter->new_x, nil, iter->new_delf, nil);
400 if (k > 0) {
401 (*(iter->gfun))(iter->new_x, iter->new_delf + n, nil, nil);
402 (*(iter->gfun))(iter->new_x, nil, iter->new_delg, nil);
403 for (i = 0; i < n; i++) {
404 for (j = 0; j < k; j++)
405 iter->new_delf[i] += iter->new_x[n + j] * iter->new_delg[j][i];
410 static void eval_hessian(Iteration *iter)
412 size_t i, j, n, k;
414 n = iter->n;
415 k = iter->k;
416 (*(iter->ffun))(iter->x, nil, nil, iter->hessf);
417 for (i = n; i < n + k; i++)
418 for (j = n; j < n + k; j++) iter->hessf[i][j] = 0.0;
421 /************************************************************************/
422 /** **/
423 /** Backtracking Line Search **/
424 /** **/
425 /************************************************************************/
427 static void
428 linesearch(Iteration *iter)
430 size_t i, n;
431 double newtlen, maxstep, initslope, rellength, lambda, minlambda,
432 lambdatemp, lambdaprev, a, b, disc, critprev, f1, f2, a11, a12, a21, a22,
433 del;
434 double *qnstep, *delf, *x, *new_x, *sx;
436 n = iter->n + iter->k;
437 if (! iter->use_line_search) {
438 iter->maxtaken = FALSE;
439 for (i = 0; i < n; i++)
440 iter->new_x[i] = iter->x[i] + iter->qnstep[i];
441 eval_next_funval(iter);
442 iter->retcode = 0;
444 else{
445 qnstep = iter->qnstep;
446 maxstep = iter->maxstep;
447 delf = (iter->k == 0) ? iter->delf : iter->F;
448 x = iter->x;
449 new_x = iter->new_x;
450 sx = iter->sx;
452 iter->maxtaken = FALSE;
453 iter->retcode = 2;
455 for (i = 0, newtlen = 0.0; i < n; i++)
456 newtlen += pow(sx[i] * qnstep[i], 2.0);
457 newtlen = sqrt(newtlen);
459 if (newtlen > maxstep) {
460 for (i = 0; i < n; i++) qnstep[i] *= (maxstep / newtlen);
461 newtlen = maxstep;
464 for (i = 0, initslope = 0.0; i < n; i++) initslope += delf[i] * qnstep[i];
466 for (i = 0, rellength = 0.0; i < n; i++)
467 rellength = Max(rellength, fabs(qnstep[i]) / Max(fabs(x[i]), 1.0 / sx[i]));
469 minlambda = (rellength == 0.0) ? 2.0 : iter->steptol / rellength;
471 lambda = 1.0;
472 lambdaprev = 1.0; /* to make compilers happy */
473 critprev = 1.0; /* to make compilers happy */
474 while (iter->retcode > 1) {
475 for (i = 0; i < n; i++) new_x[i] = x[i] + lambda * qnstep[i];
476 eval_next_funval(iter);
477 if (iter->new_crit <= iter->crit + ALPHA * lambda * initslope) {
478 iter->retcode = 0;
479 if (lambda == 1.0 && newtlen > 0.99 * maxstep) iter->maxtaken = TRUE;
481 else if (lambda < minlambda) {
482 iter->retcode = 1;
483 iter->new_crit = iter->crit;
484 for (i = 0; i < n; i++) new_x[i] = x[i];
486 else {
487 if (lambda == 1.0) { /* first backtrack, quadratic fit */
488 lambdatemp = - initslope
489 / (2 * (iter->new_crit - iter->crit - initslope));
491 else { /* all subsequent backtracks, cubic fit */
492 del = lambda - lambdaprev;
493 f1 = iter->new_crit - iter->crit - lambda * initslope;
494 f2 = critprev - iter->crit - lambdaprev * initslope;
495 a11 = 1.0 / (lambda * lambda);
496 a12 = -1.0 / (lambdaprev * lambdaprev);
497 a21 = -lambdaprev * a11;
498 a22 = -lambda * a12;
499 a = (a11 * f1 + a12 * f2) / del;
500 b = (a21 * f1 + a22 * f2) / del;
501 disc = b * b - 3 * a * initslope;
502 if (a == 0) { /* cubic is a quadratic */
503 lambdatemp = - initslope / (2 * b);
505 else { /* legitimate cubic */
506 lambdatemp = (-b + sqrt(disc)) / (3 * a);
508 lambdatemp = Min(lambdatemp, 0.5 * lambda);
510 lambdaprev = lambda;
511 critprev = iter->new_crit;
512 lambda = Max(0.1 * lambda, lambdatemp);
513 if (iter->verbose > 0) {
514 sprintf(buf, "Backtracking: lambda = %g\n", lambda);
515 PRINTSTR(buf);
522 /************************************************************************/
523 /** **/
524 /** Status Printing Functions **/
525 /** **/
526 /************************************************************************/
528 static void
531 if (iter->verbose > 0) {
532 sprintf(buf, "Iteration %d.\n", iter->itncount);
533 PRINTSTR(buf);
537 static void
538 print_status(Iteration *iter)
540 size_t i, j;
542 if (iter->verbose > 0) {
543 sprintf(buf, "Criterion value = %g\n",
544 (iter->change_sign) ? -iter->crit : iter->crit);
545 PRINTSTR(buf);
546 if (iter->verbose > 1) {
547 PRINTSTR("Location = <");
548 for (i = 0; i < iter->n + iter->k; i++) {
549 sprintf(buf, (i < iter->n + iter->k - 1) ? "%g " : "%g>\n", iter->x[i]);
550 PRINTSTR(buf);
553 if (iter->verbose > 2) {
555 for (i = 0; i < iter->n + iter->k; i++) {
556 sprintf(buf, (i < iter->n + iter->k - 1) ? "%g " : "%g>\n",
557 (iter->change_sign) ? -iter->delf[i] : iter->delf[i]);
558 PRINTSTR(buf);
561 if (iter->verbose > 3) {
562 PRINTSTR("Hessian:\n");
563 for (i = 0; i < iter->n + iter->k; i++) {
564 for (j = 0; j < iter->n + iter->k; j++) {
565 sprintf(buf, (j < iter->n + iter->k - 1) ? "%g " : "%g\n",
566 (iter->change_sign) ? -iter->hessf[i][j] : iter->hessf[i][j]);
567 PRINTSTR(buf);
572 if (iter->termcode != 0 && iter->verbose > 0) {
573 sprintf(buf, "Reason for termination: %s.\n", termcodes[iter->termcode]);
574 PRINTSTR(buf);
578 /************************************************************************/
579 /** **/
580 /** Iteration Driver **/
581 /** **/
582 /************************************************************************/
584 static void
585 findqnstep(Iteration *iter)
587 size_t i, j, N, l;
589 if (iter->k == 0) {
590 modelhess(iter->n, iter->sx, iter->hessf, iter->L, &iter->diagadd);
591 cholsolve(iter->n, iter->delf, iter->L, iter->qnstep);
593 else {
594 N = iter->n + iter->k;
595 for (i = 0; i < N; i++) {
596 for (l = 0, iter->F[i] = 0.0; l < N; l++)
597 iter->F[i] += iter->hessf[i][l] * iter->delf[l];
598 for (j = 0; j < N; j++)
599 for (l = 0, iter->H[i][j] = 0.0; l < N; l++)
600 iter->H[i][j] += iter->hessf[i][l] * iter->hessf[j][l];
602 modelhess(N, iter->sx, iter->H, iter->L, &iter->diagadd);
603 cholsolve(N, iter->F, iter->L, iter->qnstep);
607 static void iterupdate(Iteration *iter)
609 size_t i, j, n, k;
611 n = iter->n;
612 k = iter->k;
613 iter->f = iter->new_f;
614 iter->crit = iter->new_crit;
615 for (i = 0; i < n + k; i++) {
616 iter->x[i] = iter->new_x[i];
617 iter->delf[i] = iter->new_delf[i];
619 for (i = 0; i < k; i++) {
620 for (j = 0; j < n; j++) {
621 iter->hessf[n + i][j] = iter->new_delg[i][j];
622 iter->hessf[j][n + i] = iter->new_delg[i][j];
627 static void
628 mindriver(Iteration *iter)
630 iter->consecmax = 0;
631 iter->itncount = 0;
632 iter->termcode = 0;
633 if (! iter->values_supplied) {
634 eval_funval(iter);
635 if (iter->k == 0) eval_gradient(iter);
636 eval_hessian(iter);
639 stoptest0(iter);
641 print_status(iter);
642 while (iter->termcode == 0) {
643 iter->itncount++;
645 findqnstep(iter);
646 linesearch(iter);
647 if (iter->k == 0) eval_next_gradient(iter);
648 stoptest(iter);
649 iterupdate(iter);
650 eval_hessian(iter);
651 print_status(iter);
655 /************************************************************************/
656 /** **/
657 /** External Interface Routines **/
658 /** **/
659 /************************************************************************/
661 static Iteration myiter;
663 size_t
664 minworkspacesize(size_t n, size_t k)
666 size_t size;
668 /* x, new_x, sx, delf, new_delf, qnstep and F */
669 size = 7 * sizeof(double) * (n + k);
671 /* hessf, H and L */
672 size += 3 * (sizeof(double *) * (n + k)
673 + sizeof(double) * (n + k) * (n + k));
675 /* delg and new_delg */
676 size += 2 * (sizeof(double *) * k + sizeof(double) * n * k);
678 return(size);
681 char *
682 minresultstring(int code)
684 if (code <= 0) return("bad input data");
685 else if (code <= 5) return(termcodes[code]);
686 else return("unknown return code");
689 void
690 minsetup(size_t n, size_t k,
691 void (*ffun)(), void (*gfun)(),
692 double *x, double typf, double *typx, char *work)
693 /* int (*ffun)(), (*gfun)();*/
695 Iteration *iter = &myiter;
696 size_t i, j;
697 double nx0, ntypx;
699 n = (n > 0) ? n : 0;
700 k = (k > 0) ? k : 0;
702 iter->n = n;
703 iter->k = k;
704 iter->ffun = ffun; /* FIXME:AJR */
705 iter->gfun = gfun;
707 iter->x = (double *) work; work += sizeof(double) * (n + k);
708 iter->new_x = (double *) work; work += sizeof(double) * (n + k);
709 iter->sx = (double *) work; work += sizeof(double) * (n + k);
710 iter->delf = (double *) work; work += sizeof(double) * (n + k);
711 iter->new_delf = (double *) work; work += sizeof(double) * (n + k);
712 iter->qnstep = (double *) work; work += sizeof(double) * (n + k);
713 iter->F = (double *) work; work += sizeof(double) * (n + k);
714 for (i = 0; i < n; i++) {
715 iter->x[i] = x[i];
716 iter->sx[i] = (typx != nil && typx[i] > 0.0) ? 1.0 / typx[i] : 1.0;
718 for (i = 0; i < k; i++) {
719 iter->x[n + i] = x[n + i];
720 iter->sx[n + i] = 1.0;
723 iter->hessf = (double **) work; work += sizeof(double *) * (n + k);
724 for (i = 0; i < n + k; i++) {
725 iter->hessf[i] = (double *) work;
726 work += sizeof(double) * (n + k);
728 for (i = 0; i < n + k; i++)
729 for (j = 0; j < n + k; j++) iter->hessf[i][j] = 0.0;
730 iter->L = (double **) work; work += sizeof(double *) * (n + k);
731 for (i = 0; i < n + k; i++) {
732 iter->L[i] = (double *) work;
733 work += sizeof(double) * (n + k);
735 iter->H = (double **) work; work += sizeof(double *) * (n + k);
736 for (i = 0; i < n + k; i++) {
737 iter->H[i] = (double *) work;
738 work += sizeof(double) * (n + k);
741 iter->new_delg = (k > 0) ? (double **) work : nil;
742 work += sizeof(double *) * k;
743 for (i = 0; i < k; i++) {
744 iter->new_delg[i] = (double *) work;
745 work += sizeof(double) * n;
748 iter->typf = (typf > 0.0) ? typf : 1.0;
750 iter->verbose = VERBOSE;
751 iter->use_line_search = USE_SEARCH;
752 iter->itncount = 0;
753 iter->itnlimit = ITNLIMIT;
755 iter->steptol = pow(macheps(), STEPTOL_POWER);
756 for (i = 0, nx0 = 0.0, ntypx = 0.0; i < iter->n; i++) {
757 nx0 += fabs(iter->new_x[i] / iter->sx[i]);
758 ntypx += fabs(1.0 / iter->sx[i]);
760 iter->maxstep = MAX_STEP_FACTOR * Max(nx0, ntypx);
761 iter->values_supplied = FALSE;
764 void
765 minsetoptions(double gradtol, double steptol, double maxstep,
766 int itnlimit, int verbose, int use_search, int change_sign)
768 Iteration *iter = &myiter;
771 if (steptol > 0.0) iter->steptol = steptol;
772 if (maxstep > 0.0) iter->maxstep = maxstep;
773 if (itnlimit >= 0) iter->itnlimit = itnlimit;
774 if (verbose >= 0) iter->verbose = verbose;
775 iter->use_line_search = use_search;
776 iter->change_sign = change_sign;
779 void
780 minsupplyvalues(double f, double *delf, double **hessf,
781 double *g, double **delg)
783 Iteration *iter = &myiter;
784 size_t i, j, n, k;
786 n = iter->n;
787 k = iter->k;
789 iter->f = f;
790 for (i = 0; i < n; i++) {
791 iter->delf[i] = delf[i];
792 for (j = 0; j < k; j++)
793 iter->delf[i] += iter->x[n + j] * delg[j][i];
794 for (j = 0; j < n; j++) iter->hessf[i][j] = hessf[i][j];
796 for (i = 0; i < k; i++) {
797 iter->delf[n + i] = g[i];
798 for (j = 0; j < n; j++) {
799 iter->hessf[n + i][j] = delg[i][j];
800 iter->hessf[j][n + i] = delg[i][j];
803 if (k == 0) iter->crit = f;
804 else {
805 for (i = 0, iter->crit = 0.0; i < n + k; i++)
806 iter->crit += 0.5 * pow(fabs(iter->delf[i]), 2.0);
808 iter->values_supplied = TRUE;
811 void
812 minimize()
814 mindriver(&myiter);
817 void
818 minresults(double *x, double *pf, double *pcrit, double *delf, double **hessf,
819 double *g, double **delg, int *pcount, int *ptermcode, double *diagadd)
821 Iteration *iter = &myiter;
822 size_t i, j, n, k;
824 n = iter->n;
825 k = iter->k;
827 if (pf != nil) *pf = iter->f;
828 if (pcrit != nil) *pcrit = iter->crit;
829 for (i = 0; i < n; i++) {
830 if (x != nil) x[i] = iter->x[i];
831 if (delf != nil) delf[i] = iter->delf[i];
832 for (j = 0; j < n; j++) if (hessf != nil) hessf[i][j] = iter->hessf[i][j];
834 for (i = 0; i < k; i++) {
835 if (x != nil) x[n + i] = iter->x[n + i];
836 if (g != nil) g[i] = iter->delf[n + i];
837 for (j = 0; j < n; j++)
838 if (delg != nil) delg[i][j] = iter->hessf[n + i][j];
840 if (pcount != nil) *pcount = iter->itncount;
841 if (ptermcode != nil) *ptermcode = iter->termcode;