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. */
8 * Nonlinear optimization modules adapted from Dennis and Schnabel,
9 * "Numerical Methods for Unconstrained Optimization and Nonlinear
16 #define PRINTSTR(s) printf(s)
20 #define PRINTSTR(s) stdputstr(s)
23 extern double macheps();
25 /************************************************************************/
27 /** Definitions and Globals **/
29 /************************************************************************/
35 # define INIT_GRAD_FRAC .001
36 # define CONSEC_MAX_LIMIT 5
38 # define MAX_STEP_FACTOR 1000
39 # define GRADTOL_POWER 1.0 / 3.0
40 # define STEPTOL_POWER 2.0 / 3.0
43 # define USE_SEARCH TRUE
45 typedef double **RMatrix
, *RVector
;
49 int (*ffun
)(), (*gfun
)();
50 double f
, typf
, new_f
;
51 double crit
, new_crit
;
52 RVector x
, new_x
, sx
, delf
, new_delf
, qnstep
, F
;
53 RMatrix hessf
, H
, L
, new_delg
;
54 double gradtol
, steptol
, maxstep
;
55 int itncount
, itnlimit
, maxtaken
, consecmax
, retcode
, termcode
;
56 int use_line_search
, verbose
, values_supplied
, change_sign
;
60 static char *termcodes
[] = {"not yet terminated",
61 "gradient size is less than gradient tolerance",
62 "step size is less than step tolerance",
63 "no satisfactory step found in backtracking",
64 "iteration limit exceeded",
65 "maximum size step taken 5 iterations in a row"};
67 /************************************************************************/
69 /** Utility Functions **/
71 /************************************************************************/
73 static double Max(a
, b
)
76 return(a
> b
? a
: b
);
79 static double Min(a
, b
)
82 return(a
> b
? b
: a
);
85 /************************************************************************/
87 /** Cholesky Solving Functions **/
89 /************************************************************************/
91 /* solve (L L^T) s = -g for s */
92 static cholsolve(n
, g
, L
, s
)
105 for (i
= 0; i
< n
; i
++) s
[i
] = -s
[i
];
108 /* solve Ly = b for y */
109 static lsolve(n
, b
, L
, y
)
116 for (i
= 0; i
< n
; i
++) {
118 for (j
= 0; j
< i
; j
++) y
[i
] -= L
[i
][j
] * y
[j
];
119 if (L
[i
][i
] != 0) y
[i
] /= L
[i
][i
];
123 /* solve (L^T)x = y for x */
124 static ltsolve(n
, y
, L
, x
)
131 for (i
= n
- 1; i
>= 0; i
--) {
133 for (j
= i
+ 1; j
< n
; j
++) x
[i
] -= L
[j
][i
] * x
[j
];
134 if (L
[i
][i
] != 0) x
[i
] /= L
[i
][i
];
138 static modelhess(n
, sx
, H
, L
, diagadd
)
145 double sqrteps
, maxdiag
, mindiag
, maxoff
, maxoffl
, maxposdiag
, mu
,
146 maxadd
, maxev
, minev
, offrow
, sdd
;
148 /* scale H on both sides with sx */
149 for (i
= 0; i
< n
; i
++)
150 for (j
= 0; j
< n
; j
++) H
[i
][j
] /= sx
[i
] * sx
[j
];
153 * find mu to add to diagonal so no diagonal elements are negative
154 * and largest diagonal dominates largest off diagonal element in H
156 sqrteps
= sqrt(macheps());
157 for (maxdiag
= H
[0][0], i
= 1; i
< n
; i
++)
158 maxdiag
= Max(maxdiag
, H
[i
][i
]);
159 for (mindiag
= H
[0][0], i
= 1; i
< n
; i
++)
160 mindiag
= Min(mindiag
, H
[i
][i
]);
161 maxposdiag
= Max(0.0, maxdiag
);
163 if (mindiag
<= sqrteps
* maxposdiag
) {
164 mu
= 2 * (maxposdiag
- mindiag
) * sqrteps
- mindiag
;
170 for (maxoff
= fabs(H
[0][1]), i
= 0; i
< n
; i
++)
171 for (j
= i
+ 1; j
< n
; j
++)
172 maxoff
= Max(maxoff
, fabs(H
[i
][j
]));
176 if (maxoff
* (1 + 2 * sqrteps
) > maxdiag
) {
177 mu
+= (maxoff
- maxdiag
) + 2 * sqrteps
* maxoff
;
178 maxdiag
= maxoff
* (1 + 2 * sqrteps
);
181 if (maxdiag
== 0.0) {
186 if (mu
> 0) for (i
= 0; i
< n
; i
++) H
[i
][i
] += mu
;
188 maxoffl
= sqrt(Max(maxdiag
, maxoff
/ n
));
191 * compute the perturbed Cholesky decomposition
193 for (i
= 0; i
< n
; i
++)
194 for (j
= 0; j
< n
; j
++) L
[i
][j
] = H
[i
][j
];
195 choldecomp(L
, n
, maxoffl
, &maxadd
);
198 * add something to diagonal, if needed, to make positive definite
199 * and recompute factorization
204 for (i
= 0; i
< n
; i
++) {
205 for (offrow
= 0.0, j
= 0; j
< n
; j
++)
206 if (i
!= j
) offrow
+= fabs(H
[i
][j
]);
207 maxev
= Max(maxev
, H
[i
][i
] + offrow
);
208 minev
= Min(minev
, H
[i
][i
] - offrow
);
210 sdd
= (maxev
- minev
) * sqrteps
- minev
;
212 mu
= Min(maxadd
, sdd
);
213 for (i
= 0; i
< n
; i
++) H
[i
][i
] += mu
;
214 for (i
= 0; i
< n
; i
++)
215 for (j
= 0; j
< n
; j
++) L
[i
][j
] = H
[i
][j
];
216 choldecomp(L
, n
, maxoffl
, &maxadd
);
221 /* unscale H and L */
222 for (i
= 0; i
< n
; i
++)
223 for (j
= 0; j
< n
; j
++) {
224 H
[i
][j
] *= sx
[i
] * sx
[j
];
229 /************************************************************************/
231 /** Stopping Criteria **/
233 /************************************************************************/
235 static double gradsize(iter
, new)
240 double size
, term
, crit
, typf
;
243 n
= iter
->n
+ iter
->k
;
248 delf
= (new) ? iter
->new_delf
: iter
->delf
;
250 for (i
= 0, size
= 0.0; i
< n
; i
++) {
251 term
= fabs(delf
[i
]) * Max(x
[i
], 1.0 / sx
[i
]) / Max(fabs(crit
), typf
);
252 size
= Max(size
, term
);
257 static double incrsize(iter
)
262 RVector x
, new_x
, sx
;
265 n
= iter
->n
+ iter
->k
;
269 for (i
= 0, size
= 0.0; i
< n
; i
++) {
270 term
= fabs(x
[i
] - new_x
[i
]) / Max(fabs(x
[i
]), 1.0 / sx
[i
]);
271 size
= Max(size
, term
);
276 static stoptest0(iter
)
281 if (gradsize(iter
, FALSE
) <= INIT_GRAD_FRAC
* iter
->gradtol
)
283 else iter
->termcode
= 0;
285 return(iter
->termcode
);
288 static stoptest(iter
)
291 int termcode
, retcode
, itncount
, itnlimit
, maxtaken
, consecmax
;
292 double gradtol
, steptol
;
294 retcode
= iter
->retcode
;
295 gradtol
= iter
->gradtol
;
296 steptol
= iter
->steptol
;
297 itncount
= iter
->itncount
;
298 itnlimit
= iter
->itnlimit
;
299 maxtaken
= iter
->maxtaken
;
300 consecmax
= iter
->consecmax
;
303 if (retcode
== 1) termcode
= 3;
304 else if (gradsize(iter
, TRUE
) <= gradtol
) termcode
= 1;
305 else if (incrsize(iter
) <= steptol
) termcode
= 2;
306 else if (itncount
>= itnlimit
) termcode
= 4;
309 if (consecmax
>= CONSEC_MAX_LIMIT
) termcode
= 5;
313 iter
->consecmax
= consecmax
;
314 iter
->termcode
= termcode
;
319 /************************************************************************/
321 /** Function and Derivative Evaluation **/
323 /************************************************************************/
325 static eval_funval(iter
)
330 (*(iter
->ffun
))(iter
->x
, &iter
->f
, nil
, nil
);
331 if (iter
->k
== 0) iter
->crit
= iter
->f
;
334 for (i
= 0, iter
->crit
= 0.0; i
< iter
->n
+ iter
->k
; i
++)
335 iter
->crit
+= 0.5 * pow(fabs(iter
->delf
[i
]), 2.0);
339 static eval_next_funval(iter
)
344 (*(iter
->ffun
))(iter
->new_x
, &iter
->new_f
, nil
, nil
);
345 if (iter
->k
== 0) iter
->new_crit
= iter
->new_f
;
347 eval_next_gradient(iter
);
348 for (i
= 0, iter
->new_crit
= 0.0; i
< iter
->n
+ iter
->k
; i
++)
349 iter
->new_crit
+= 0.5 * pow(fabs(iter
->new_delf
[i
]), 2.0);
353 static eval_gradient(iter
)
361 (*(iter
->ffun
))(iter
->x
, nil
, iter
->delf
, nil
);
363 (*(iter
->gfun
))(iter
->x
, iter
->delf
+ n
, nil
, nil
);
364 (*(iter
->gfun
))(iter
->x
, nil
, iter
->new_delg
, nil
);
365 for (i
= 0; i
< n
; i
++) {
366 for (j
= 0; j
< k
; j
++)
367 iter
->delf
[i
] += iter
->x
[n
+ j
] * iter
->new_delg
[j
][i
];
368 for (j
= 0; j
< k
; j
++) {
369 iter
->hessf
[n
+ j
][i
] = iter
->new_delg
[j
][i
];
370 iter
->hessf
[i
][n
+ j
] = iter
->new_delg
[j
][i
];
376 static eval_next_gradient(iter
)
383 (*(iter
->ffun
))(iter
->new_x
, nil
, iter
->new_delf
, nil
);
385 (*(iter
->gfun
))(iter
->new_x
, iter
->new_delf
+ n
, nil
, nil
);
386 (*(iter
->gfun
))(iter
->new_x
, nil
, iter
->new_delg
, nil
);
387 for (i
= 0; i
< n
; i
++) {
388 for (j
= 0; j
< k
; j
++)
389 iter
->new_delf
[i
] += iter
->new_x
[n
+ j
] * iter
->new_delg
[j
][i
];
394 static eval_hessian(iter
)
401 (*(iter
->ffun
))(iter
->x
, nil
, nil
, iter
->hessf
);
402 for (i
= n
; i
< n
+ k
; i
++)
403 for (j
= n
; j
< n
+ k
; j
++) iter
->hessf
[i
][j
] = 0.0;
406 /************************************************************************/
408 /** Backtracking Line Search **/
410 /************************************************************************/
412 static linesearch(iter
)
416 double newtlen
, maxstep
, initslope
, rellength
, lambda
, minlambda
,
417 lambdatemp
, lambdaprev
, a
, b
, disc
, critprev
, f1
, f2
, a11
, a12
, a21
, a22
,
419 RVector qnstep
, delf
, x
, new_x
, sx
;
421 n
= iter
->n
+ iter
->k
;
422 if (! iter
->use_line_search
) {
423 iter
->maxtaken
= FALSE
;
424 for (i
= 0; i
< n
; i
++)
425 iter
->new_x
[i
] = iter
->x
[i
] + iter
->qnstep
[i
];
426 eval_next_funval(iter
);
430 qnstep
= iter
->qnstep
;
431 maxstep
= iter
->maxstep
;
432 delf
= (iter
->k
== 0) ? iter
->delf
: iter
->F
;
437 iter
->maxtaken
= FALSE
;
440 for (i
= 0, newtlen
= 0.0; i
< n
; i
++)
441 newtlen
+= pow(sx
[i
] * qnstep
[i
], 2.0);
442 newtlen
= sqrt(newtlen
);
444 if (newtlen
> maxstep
) {
445 for (i
= 0; i
< n
; i
++) qnstep
[i
] *= (maxstep
/ newtlen
);
449 for (i
= 0, initslope
= 0.0; i
< n
; i
++) initslope
+= delf
[i
] * qnstep
[i
];
451 for (i
= 0, rellength
= 0.0; i
< n
; i
++)
452 rellength
= Max(rellength
, fabs(qnstep
[i
]) / Max(fabs(x
[i
]), 1.0 / sx
[i
]));
454 minlambda
= (rellength
== 0.0) ? 2.0 : iter
->steptol
/ rellength
;
457 lambdaprev
= 1.0; /* to make compilers happy */
458 critprev
= 1.0; /* to make compilers happy */
459 while (iter
->retcode
> 1) {
460 for (i
= 0; i
< n
; i
++) new_x
[i
] = x
[i
] + lambda
* qnstep
[i
];
461 eval_next_funval(iter
);
462 if (iter
->new_crit
<= iter
->crit
+ ALPHA
* lambda
* initslope
) {
464 if (lambda
== 1.0 && newtlen
> 0.99 * maxstep
) iter
->maxtaken
= TRUE
;
466 else if (lambda
< minlambda
) {
468 iter
->new_crit
= iter
->crit
;
469 for (i
= 0; i
< n
; i
++) new_x
[i
] = x
[i
];
472 if (lambda
== 1.0) { /* first backtrack, quadratic fit */
473 lambdatemp
= - initslope
474 / (2 * (iter
->new_crit
- iter
->crit
- initslope
));
476 else { /* all subsequent backtracks, cubic fit */
477 del
= lambda
- lambdaprev
;
478 f1
= iter
->new_crit
- iter
->crit
- lambda
* initslope
;
479 f2
= critprev
- iter
->crit
- lambdaprev
* initslope
;
480 a11
= 1.0 / (lambda
* lambda
);
481 a12
= -1.0 / (lambdaprev
* lambdaprev
);
482 a21
= -lambdaprev
* a11
;
484 a
= (a11
* f1
+ a12
* f2
) / del
;
485 b
= (a21
* f1
+ a22
* f2
) / del
;
486 disc
= b
* b
- 3 * a
* initslope
;
487 if (a
== 0) { /* cubic is a quadratic */
488 lambdatemp
= - initslope
/ (2 * b
);
490 else { /* legitimate cubic */
491 lambdatemp
= (-b
+ sqrt(disc
)) / (3 * a
);
493 lambdatemp
= Min(lambdatemp
, 0.5 * lambda
);
496 critprev
= iter
->new_crit
;
497 lambda
= Max(0.1 * lambda
, lambdatemp
);
498 if (iter
->verbose
> 0) {
499 sprintf(buf
, "Backtracking: lambda = %g\n", lambda
);
507 /************************************************************************/
509 /** Status Printing Functions **/
511 /************************************************************************/
513 static print_header(iter
)
516 if (iter
->verbose
> 0) {
517 sprintf(buf
, "Iteration %d.\n", iter
->itncount
);
522 static print_status(iter
)
527 if (iter
->verbose
> 0) {
528 sprintf(buf
, "Criterion value = %g\n",
529 (iter
->change_sign
) ? -iter
->crit
: iter
->crit
);
531 if (iter
->verbose
> 1) {
532 PRINTSTR("Location = <");
533 for (i
= 0; i
< iter
->n
+ iter
->k
; i
++) {
534 sprintf(buf
, (i
< iter
->n
+ iter
->k
- 1) ? "%g " : "%g>\n", iter
->x
[i
]);
538 if (iter
->verbose
> 2) {
539 PRINTSTR("Gradient = <");
540 for (i
= 0; i
< iter
->n
+ iter
->k
; i
++) {
541 sprintf(buf
, (i
< iter
->n
+ iter
->k
- 1) ? "%g " : "%g>\n",
542 (iter
->change_sign
) ? -iter
->delf
[i
] : iter
->delf
[i
]);
546 if (iter
->verbose
> 3) {
547 PRINTSTR("Hessian:\n");
548 for (i
= 0; i
< iter
->n
+ iter
->k
; i
++) {
549 for (j
= 0; j
< iter
->n
+ iter
->k
; j
++) {
550 sprintf(buf
, (j
< iter
->n
+ iter
->k
- 1) ? "%g " : "%g\n",
551 (iter
->change_sign
) ? -iter
->hessf
[i
][j
] : iter
->hessf
[i
][j
]);
557 if (iter
->termcode
!= 0 && iter
->verbose
> 0) {
558 sprintf(buf
, "Reason for termination: %s.\n", termcodes
[iter
->termcode
]);
563 /************************************************************************/
565 /** Iteration Driver **/
567 /************************************************************************/
569 static findqnstep(iter
)
575 modelhess(iter
->n
, iter
->sx
, iter
->hessf
, iter
->L
, &iter
->diagadd
);
576 cholsolve(iter
->n
, iter
->delf
, iter
->L
, iter
->qnstep
);
579 N
= iter
->n
+ iter
->k
;
580 for (i
= 0; i
< N
; i
++) {
581 for (l
= 0, iter
->F
[i
] = 0.0; l
< N
; l
++)
582 iter
->F
[i
] += iter
->hessf
[i
][l
] * iter
->delf
[l
];
583 for (j
= 0; j
< N
; j
++)
584 for (l
= 0, iter
->H
[i
][j
] = 0.0; l
< N
; l
++)
585 iter
->H
[i
][j
] += iter
->hessf
[i
][l
] * iter
->hessf
[j
][l
];
587 modelhess(N
, iter
->sx
, iter
->H
, iter
->L
, &iter
->diagadd
);
588 cholsolve(N
, iter
->F
, iter
->L
, iter
->qnstep
);
592 static iterupdate(iter
)
599 iter
->f
= iter
->new_f
;
600 iter
->crit
= iter
->new_crit
;
601 for (i
= 0; i
< n
+ k
; i
++) {
602 iter
->x
[i
] = iter
->new_x
[i
];
603 iter
->delf
[i
] = iter
->new_delf
[i
];
605 for (i
= 0; i
< k
; i
++) {
606 for (j
= 0; j
< n
; j
++) {
607 iter
->hessf
[n
+ i
][j
] = iter
->new_delg
[i
][j
];
608 iter
->hessf
[j
][n
+ i
] = iter
->new_delg
[i
][j
];
613 static mindriver(iter
)
619 if (! iter
->values_supplied
) {
621 if (iter
->k
== 0) eval_gradient(iter
);
628 while (iter
->termcode
== 0) {
633 if (iter
->k
== 0) eval_next_gradient(iter
);
641 /************************************************************************/
643 /** External Interface Routines **/
645 /************************************************************************/
647 static Iteration myiter
;
649 minworkspacesize(n
, k
)
654 /* x, new_x, sx, delf, new_delf, qnstep and F */
655 size
= 7 * sizeof(double) * (n
+ k
);
658 size
+= 3 * (sizeof(double *) * (n
+ k
)
659 + sizeof(double) * (n
+ k
) * (n
+ k
));
661 /* delg and new_delg */
662 size
+= 2 * (sizeof(double *) * k
+ sizeof(double) * n
* k
);
667 char *minresultstring(code
)
670 if (code
<= 0) return("bad input data");
671 else if (code
<= 5) return(termcodes
[code
]);
672 else return("unknown return code");
675 minsetup(n
, k
, ffun
, gfun
, x
, typf
, typx
, work
)
676 int n
, k
, (*ffun
)(), (*gfun
)();
681 Iteration
*iter
= &myiter
;
693 iter
->x
= (RVector
) work
; work
+= sizeof(double) * (n
+ k
);
694 iter
->new_x
= (RVector
) work
; work
+= sizeof(double) * (n
+ k
);
695 iter
->sx
= (RVector
) work
; work
+= sizeof(double) * (n
+ k
);
696 iter
->delf
= (RVector
) work
; work
+= sizeof(double) * (n
+ k
);
697 iter
->new_delf
= (RVector
) work
; work
+= sizeof(double) * (n
+ k
);
698 iter
->qnstep
= (RVector
) work
; work
+= sizeof(double) * (n
+ k
);
699 iter
->F
= (RVector
) work
; work
+= sizeof(double) * (n
+ k
);
700 for (i
= 0; i
< n
; i
++) {
702 iter
->sx
[i
] = (typx
!= nil
&& typx
[i
] > 0.0) ? 1.0 / typx
[i
] : 1.0;
704 for (i
= 0; i
< k
; i
++) {
705 iter
->x
[n
+ i
] = x
[n
+ i
];
706 iter
->sx
[n
+ i
] = 1.0;
709 iter
->hessf
= (RMatrix
) work
; work
+= sizeof(double *) * (n
+ k
);
710 for (i
= 0; i
< n
+ k
; i
++) {
711 iter
->hessf
[i
] = (RVector
) work
;
712 work
+= sizeof(double) * (n
+ k
);
714 for (i
= 0; i
< n
+ k
; i
++)
715 for (j
= 0; j
< n
+ k
; j
++) iter
->hessf
[i
][j
] = 0.0;
716 iter
->L
= (RMatrix
) work
; work
+= sizeof(double *) * (n
+ k
);
717 for (i
= 0; i
< n
+ k
; i
++) {
718 iter
->L
[i
] = (RVector
) work
;
719 work
+= sizeof(double) * (n
+ k
);
721 iter
->H
= (RMatrix
) work
; work
+= sizeof(double *) * (n
+ k
);
722 for (i
= 0; i
< n
+ k
; i
++) {
723 iter
->H
[i
] = (RVector
) work
;
724 work
+= sizeof(double) * (n
+ k
);
727 iter
->new_delg
= (k
> 0) ? (RMatrix
) work
: nil
;
728 work
+= sizeof(double *) * k
;
729 for (i
= 0; i
< k
; i
++) {
730 iter
->new_delg
[i
] = (RVector
) work
;
731 work
+= sizeof(double) * n
;
734 iter
->typf
= (typf
> 0.0) ? typf
: 1.0;
736 iter
->verbose
= VERBOSE
;
737 iter
->use_line_search
= USE_SEARCH
;
739 iter
->itnlimit
= ITNLIMIT
;
740 iter
->gradtol
= pow(macheps(), GRADTOL_POWER
);
741 iter
->steptol
= pow(macheps(), STEPTOL_POWER
);
742 for (i
= 0, nx0
= 0.0, ntypx
= 0.0; i
< iter
->n
; i
++) {
743 nx0
+= fabs(iter
->new_x
[i
] / iter
->sx
[i
]);
744 ntypx
+= fabs(1.0 / iter
->sx
[i
]);
746 iter
->maxstep
= MAX_STEP_FACTOR
* Max(nx0
, ntypx
);
747 iter
->values_supplied
= FALSE
;
750 minsetoptions(gradtol
, steptol
, maxstep
, itnlimit
, verbose
, use_search
, change_sign
)
751 double gradtol
, steptol
, maxstep
;
752 int itnlimit
, verbose
, use_search
, change_sign
;
754 Iteration
*iter
= &myiter
;
756 if (gradtol
> 0.0) iter
->gradtol
= gradtol
;
757 if (steptol
> 0.0) iter
->steptol
= steptol
;
758 if (maxstep
> 0.0) iter
->maxstep
= maxstep
;
759 if (itnlimit
>= 0) iter
->itnlimit
= itnlimit
;
760 if (verbose
>= 0) iter
->verbose
= verbose
;
761 iter
->use_line_search
= use_search
;
762 iter
->change_sign
= change_sign
;
765 minsupplyvalues(f
, delf
, hessf
, g
, delg
)
770 Iteration
*iter
= &myiter
;
777 for (i
= 0; i
< n
; i
++) {
778 iter
->delf
[i
] = delf
[i
];
779 for (j
= 0; j
< k
; j
++)
780 iter
->delf
[i
] += iter
->x
[n
+ j
] * delg
[j
][i
];
781 for (j
= 0; j
< n
; j
++) iter
->hessf
[i
][j
] = hessf
[i
][j
];
783 for (i
= 0; i
< k
; i
++) {
784 iter
->delf
[n
+ i
] = g
[i
];
785 for (j
= 0; j
< n
; j
++) {
786 iter
->hessf
[n
+ i
][j
] = delg
[i
][j
];
787 iter
->hessf
[j
][n
+ i
] = delg
[i
][j
];
790 if (k
== 0) iter
->crit
= f
;
792 for (i
= 0, iter
->crit
= 0.0; i
< n
+ k
; i
++)
793 iter
->crit
+= 0.5 * pow(fabs(iter
->delf
[i
]), 2.0);
795 iter
->values_supplied
= TRUE
;
798 minimize() { mindriver(&myiter
); }
800 minresults(x
, pf
, pcrit
, delf
, hessf
, g
, delg
, pcount
, ptermcode
, diagadd
)
802 double *pf
, *pcrit
, *diagadd
;
804 int *pcount
, *ptermcode
;
806 Iteration
*iter
= &myiter
;
812 if (pf
!= nil
) *pf
= iter
->f
;
813 if (pcrit
!= nil
) *pcrit
= iter
->crit
;
814 for (i
= 0; i
< n
; i
++) {
815 if (x
!= nil
) x
[i
] = iter
->x
[i
];
816 if (delf
!= nil
) delf
[i
] = iter
->delf
[i
];
817 for (j
= 0; j
< n
; j
++) if (hessf
!= nil
) hessf
[i
][j
] = iter
->hessf
[i
][j
];
819 for (i
= 0; i
< k
; i
++) {
820 if (x
!= nil
) x
[n
+ i
] = iter
->x
[n
+ i
];
821 if (g
!= nil
) g
[i
] = iter
->delf
[n
+ i
];
822 for (j
= 0; j
< n
; j
++)
823 if (delg
!= nil
) delg
[i
][j
] = iter
->hessf
[n
+ i
][j
];
825 if (pcount
!= nil
) *pcount
= iter
->itncount
;
826 if (ptermcode
!= nil
) *ptermcode
= iter
->termcode
;
827 if (diagadd
!= nil
) *diagadd
= iter
->diagadd
;
831 double pdlogdet(n
, H
)
836 double logdet
, maxadd
;
838 choldecomp(H
, n
, 0.0, &maxadd
);
839 for (i
= 0, logdet
= 0.0; i
< n
; i
++) logdet
+= 2.0 * log(H
[i
][i
]);
844 return amount added to make pos definite
845 scaling
for constraints
846 alternate global strategies
847 callback function
for verbose mode