[tkralphs]

Problem

In solvers2, hvcompd and hvcomps (hvpcompd_ew_ASL, hvpcomps_ew_ASL) can return the Hessian-vector product for an earlier point. If no gradient was computed before, they return 0.

The cause is the test that decides whether the gradient at the current x is still needed. It uses ncxval[co] for constraints and noxval[no] for objectives, but those count function evaluations. Gradient evaluations are counted in ncxgval and noxgval, which xpsg_check_ASL uses for the same decision.

The routines compute the function value at x first when it isn't known. After that the test always says the gradient is current, so it isn't computed. The product is then built from whatever the last gradient evaluation left.

A smaller slip appears in both solvers and solvers2. The test for whether the objective value is known uses ncxval[no], a constraint's counter, instead of noxval[no]. In solvers the gradient test uses f->nxval, so the result is right there.

Fix

The fix changes 8 lines in pshvprod.c:

• solvers2: use ncxgval[co] and noxgval[no] for the gradient tests, and noxval[no] for the objective value;
• solvers: use noxval[no] for the objective value.

Reproducer

This uses sqrt0.nl from #31: min sqrt(x^2+y^2) s.t. sqrt(x^2+y^2) >= 0. The program computes values and gradients at (1,1), then only the values at (3,4), then calls hvcompd and hvcomps with p = (1,2).

stale.c

/* stale: hvcompd and hvcomps for sqrt0.nl at x = (3,4), after the values
 * and gradients were computed at (1,1) and only the values at (3,4).
 * The right product with p = (1,2) is (-0.064, 0.048); at (1,1) it would
 * be (-0.3536, 0.3536).  Compile with -DASL2 for solvers2.
 * usage: stale stub obj|con */
#include "asl_pfgh.h"
#include <string.h>

#ifdef ASL2
#define A asl->i.Ew0
#else
#define A (ASL*)asl
#endif

int main(int argc, char **argv)
{
	ASL *asl;
	FILE *nl;
	real x1[2] = {1., 1.}, x2[2] = {3., 4.}, c[1], G[2], H[2], p[2] = {1., 2.};
	fint ne = 0;
	int co;
	varno_t z[2];

	asl = ASL_alloc(ASL_read_pfgh);
	nl = jac0dim(argv[1], (fint)strlen(argv[1]));
	pfgh_read(nl, 0);
	co = strcmp(argv[2], "obj") ? 0 : -1;
	conval(x1, c, &ne);
	objval(0, x1, &ne);
	jacval(x1, G, &ne);
	objgrd(0, x1, G, &ne);
	conval(x2, c, &ne);
	objval(0, x2, &ne);
	asl->p.Hvcompd(A, H, p, co);
	printf("hvcompd %s: %.17g %.17g\n", argv[2], H[0], H[1]);
	asl->p.Hvcomps(A, H, p, co, 2, z);
	printf("hvcomps %s: %.17g %.17g\n", argv[2], H[0], H[1]);
	return 0;
	}

The correct product at (3,4) is (−0.064, 0.048). At (1,1) it is (−0.3536, 0.3536).

The same values came out for both routines, for the objective (obj) and for the constraint (con):

Library	Before	After
solvers2 (20260107)	−0.35355339059327373, 0.35355339059327373	−0.064000000000000001, 0.048000000000000043
solvers (20260109)	−0.064000000000000001, 0.048000000000000043	unchanged

Other checks

I also ran a larger set of test cases at x = (3,4), read with pfgh_read:

• every Hessian routine: sphes, fullhes, duthes, hvinit/hvcomp, hvcompd and hvcomps, plus the nerror variants;
• for the objective and for the constraint;
• after five different sets of earlier calls: none, the values, the constraint values only, the values and gradients, and the values and hvinit.

Results:

• solvers2 before: hvcompd/hvcomps returned 0 in 14 of the 20 cases. They were right only after the gradients had been computed, or, for the objective, by accident after only the constraint values.
• solvers2 after: all 20 cases are right.
• Everything else: unchanged in both libraries. The full set is 132 cases in solvers and 129 in solvers2, and includes every result of the other Hessian routines, the jacval/congrd/objgrd tests from Restore the caller's err_jmp on return from the evaluation routines #31, and the error tests at x = (0,0).

This applies to master by itself and doesn't depend on #31. I found it while testing the Hessian routines for #31.

The build and test setup is the same as in #31. That means Linux x86_64 with gcc 11.4, compiling the library sources directly with gcc -O2 (solvers2 from the a = list in makefile.u) rather than through CMake.

I didn't add an entry to solvers/changes, so that this PR doesn't conflict with the one in #31. Suggested text:

20261029
  pshvprod.c (solvers2):  fix bugs with hvcompd and hvcomps at a point
where the function values but not the gradients had been computed:
they returned the Hessian-vector product at the point of the last
gradient evaluation (or 0).  pshvprod.c (solvers and solvers2):  test
noxval rather than ncxval for the objective value.

Related, not changed here

I noticed these in solvers while testing:

• fullhese macro: in solvers/asl.h it calls p.Fulhes, the variant without nerror, with seven arguments. As a result, fullhese(...) doesn't compile ("too many arguments to function 'asl->p.Fulhes'"). It should call p.Fulhese, as in solvers2/asl.h.
• Missing macro: solvers/asl.h has no sphese macro. solvers2/asl.h has one.
• Hvcompe never installed: pfg_read.c doesn't set p.Hvcompe, although hvpcompe_ASL exists in pshvprod.c and Hvinite is set. With pfgh_read, asl->p.Hvcompe stops with "hvcompe called before pfgh_read or fgh_read." solvers2's pfghread.c does set it.
• Wrong name in the error message: when one of the Hessian routines finds a derivative error while evaluating values or gradients itself, the solvers message reads Error evaluating "var =" definition -1: can't evaluate sqrt'(0). This happens for the objective and for a constraint. solvers2 names the objective gradient or the constraint gradient.

🤖 Generated with Claude Code