Merged updates from netlib up to:

20241108
  solvers2.tgz:  fix some bugs in parallel evaluations; new
possibilities for doing Hessian evaluation and/or sparsity
determination in parallel, controlled by new values in psinfo.h:

	int hesevalth;	/* number of threads for sphes_ew_ASL */
	int hessetupth;	/* number of threads for sphes_setup_ew_ASL */
	int hesvecth;	/* number of threads for Hessian-vector products */
	int sph_opts;	/* options affecting threaded sphes and sphes_setup */

Partially separable structure: Hessian is the sum of products
		matrix-transpose * internal Hessian * matrix.

Current sph_opts bits:
	1 = assume internal Hessians are dense
	2 = print dimensions internal Hessians
	4 = suppress parallel setup for Hessian evaluations
	8 = suppress parallel Hessian evaluations
	16 = 0x10 = print number of threads for setting up Hessian evaluations
	32 = 0x20 = print number of threads for doing Hessian evaluations
	64 = 0x40 = print number of threads used for Hessian-vector products
	128 = 0x80 = suppress parallel Hessian-vector products
	256 = 0x100 = allow parallel evaluation of internal Hessians
			used in Hessian-vector products.  The default is not
			to do this in parallel because it slows evaluations
			on a large example.

Author: mapgccv

src/examples/lin1.c

@@ -26,7 +26,7 @@ THIS SOFTWARE.
 
 #include "asl.h"
 
-main(int argc, char **argv)
+int main(int argc, char **argv)
 {
 	FILE *nl;
 	int i;

src/examples/lin2.c

@@ -26,7 +26,7 @@ THIS SOFTWARE.
 
 #include "asl.h"
 
-main(int argc, char **argv)
+int main(int argc, char **argv)
 {
 	FILE *nl;
 	int i, j, je;

src/examples/lin3.c

@@ -26,7 +26,7 @@ THIS SOFTWARE.
 
 #include "asl.h"
 
-main(int argc, char **argv)
+int main(int argc, char **argv)
 {
 	FILE *nl;
 	int i, j, je;

src/solvers/asldate.c

@@ -1 +1 @@
-long ASLdate_ASL = 20231111;
+long ASLdate_ASL = 20241111;

src/solvers/changes

@@ -3349,3 +3349,87 @@ source file jac2dim.c.
   solvers2.tgz: minor tweaks as in solvers.tgz, plus a bug fix to
 pfghread.c that affected Hessians of imported functions having
 several arguments involving variables.
+
+20240106
+  solvers2.tgz:  fix a bug with constant terms that seemingly involve
+variables, such as (v-v)^2; the terms were ignored instead of
+contributing to the objective or constraint in which they appeared,
+which mattered when the terms were nonzero, such as cos(v-v)^2.
+
+20240527
+  jacdim.c:  jac2dim_ should invoke jacpdim_ASL rather than (the
+unavailable) jac2dim_ASL.
+  mpec_adj.c:  fix a bug with complementarities involving double
+inequalities on a linear expression.
+
+20240528
+  pfghread.c in solvers2.tgz:  fix a bug (e.g., fault) seen in
+some Hessian computations or Hessian-vector products.
+
+0240601
+  mpec_adj.c in solvers2.tgz:  bug fix affecting ASL_cc_simplify.
+
+20240603
+  sphes.c in solvers2.tgz:  fix a glitch with -DALLOW_OPENMP.
+
+20240618
+  solvers2.tgz:  fix various bugs related to parallel evaluations.
+Now we assume pthreads if neither ALLOW_OPENMP nor NO_MBLK_LOCK is
+#defined.  Add the possibility of computing Hessian-vector products in
+parallel using at most asl->P.hesvecth processors.
+
+20240622
+  solvers2.tgz:  fix various more bugs related to parallel
+evaluations.  New int asl->P.thusedsetup records the number of threads
+used in the last call on sphes_setup (which is called automatically if
+needed and not yet called).  This joins asl->P.thused for the number
+of threads used in the most recent sphes call and asl->P.thusedhv for
+the number of threads used in the most recent Hessian-vector product.
+These numbers can be requested to be printed when they change by
+setting bits 16, 32, or 64 of asl->P.sphopts for setup, sphes, and
+(for 64) the various Hessian-vector product functions.
+
+20240926
+  Fix a bug in obj_adj.c that only matters when suf_sos is called
+and the model is of the form
+
+	...
+	var v;
+	s.t. vdef: v = ...;
+	minimize obj: v;
+
+(If instead one simply declared "minimize obj: ...;" then the bug did
+not bite.)  This fix affects both solvers.tgz and solvers2.tgz.
+
+20241018
+  Fix bugs in sphes.c and pshvprod.c of solvers2.tgz affecting
+parallel Hessian and Hessian-vector computations.
+
+20241108
+  solvers2.tgz:  fix some bugs in parallel evaluations; new
+possibilities for doing Hessian evaluation and/or sparsity
+determination in parallel, controlled by new values in psinfo.h:
+
+	int hesevalth;	/* number of threads for sphes_ew_ASL */
+	int hessetupth;	/* number of threads for sphes_setup_ew_ASL */
+	int hesvecth;	/* number of threads for Hessian-vector products */
+	int sph_opts;	/* options affecting threaded sphes and sphes_setup */
+
+Partially separable structure: Hessian is the sum of products
+		matrix-transpose * internal Hessian * matrix.
+
+Current sph_opts bits:
+	1 = assume internal Hessians are dense
+	2 = print dimensions internal Hessians
+	4 = suppress parallel setup for Hessian evaluations
+	8 = suppress parallel Hessian evaluations
+	16 = 0x10 = print number of threads for setting up Hessian evaluations
+	32 = 0x20 = print number of threads for doing Hessian evaluations
+	64 = 0x40 = print number of threads used for Hessian-vector products
+	128 = 0x80 = suppress parallel Hessian-vector products
+	256 = 0x100 = allow parallel evaluation of internal Hessians
+			used in Hessian-vector products.  The default is not
+			to do this in parallel because it slows evaluations
+			on a large example.
+
+Note that parallel Hessian evaluations are not available in solvers.tgz.

src/solvers/dtoa.c

@@ -28,7 +28,8 @@
  * does this with many compilers.  Whether this or another call is
  * appropriate depends on the compiler; for this to work, it may be
  * necessary to #include "float.h" or another system-dependent header
- * file.
+ * file.  When needed, this call avoids double rounding, which can
+ * cause one bit errors, e.g., with strtod on 8.3e26 or 6.3876e-16.
  */
 
 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
@@ -2716,13 +2717,14 @@ enum {	/* rounding values: same as FLT_ROUNDS */
 	};
 
  void
-gethex( const char **sp, U *rvp, int rounding, int sign MTd)
+gethex(const char **sp, U *rvp, int rounding, int sign MTd)
 {
 	Bigint *b;
+	char d;
 	const unsigned char *decpt, *s0, *s, *s1;
 	Long e, e1;
 	ULong L, lostbits, *x;
-	int big, denorm, esign, havedig, k, n, nbits, up, zret;
+	int big, denorm, esign, havedig, k, n, nb, nbits, nz, up, zret;
 #ifdef IBM
 	int j;
 #endif
@@ -2740,6 +2742,9 @@ gethex( const char **sp, U *rvp, int rounding, int sign MTd)
 #endif
 #endif /*}}*/
 		};
+#ifdef IEEE_Arith
+	int check_denorm = 0;
+#endif
 #ifdef USE_LOCALE
 	int i;
 #ifdef NO_LOCALE_CACHE
@@ -2891,7 +2896,7 @@ gethex( const char **sp, U *rvp, int rounding, int sign MTd)
 		k++;
 	b = Balloc(k MTa);
 	x = b->x;
-	n = 0;
+	havedig = n = nz = 0;
 	L = 0;
 #ifdef USE_LOCALE
 	for(i = 0; decimalpoint[i+1]; ++i);
@@ -2906,22 +2911,28 @@ gethex( const char **sp, U *rvp, int rounding, int sign MTd)
 		if (*--s1 == '.')
 			continue;
 #endif
+		if ((d = hexdig[*s1]))
+			havedig = 1;
+		else if (!havedig) {
+			e += 4;
+			continue;
+			}
 		if (n == ULbits) {
 			*x++ = L;
 			L = 0;
 			n = 0;
 			}
-		L |= (hexdig[*s1] & 0x0f) << n;
+		L |= (d & 0x0f) << n;
 		n += 4;
 		}
 	*x++ = L;
 	b->wds = n = x - b->x;
-	n = ULbits*n - hi0bits(L);
+	nb = ULbits*n - hi0bits(L);
 	nbits = Nbits;
 	lostbits = 0;
 	x = b->x;
-	if (n > nbits) {
-		n -= nbits;
+	if (nb > nbits) {
+		n = nb - nbits;
 		if (any_on(b,n)) {
 			lostbits = 1;
 			k = n - 1;
@@ -2934,8 +2945,8 @@ gethex( const char **sp, U *rvp, int rounding, int sign MTd)
 		rshift(b, n);
 		e += n;
 		}
-	else if (n < nbits) {
-		n = nbits - n;
+	else if (nb < nbits) {
+		n = nbits - nb;
 		b = lshift(b, n MTa);
 		e -= n;
 		x = b->x;
@@ -2990,12 +3001,49 @@ gethex( const char **sp, U *rvp, int rounding, int sign MTd)
 			return;
 			}
 		k = n - 1;
+#ifdef IEEE_Arith
+		if (!k) {
+			switch(rounding) {
+			  case Round_near:
+				if (((b->x[0] & 3) == 3) || (lostbits && (b->x[0] & 1))) {
+					multadd(b, 1, 1 MTa);
+ emin_check:
+					if (b->x[1] == (1 << (Exp_shift + 1))) {
+						rshift(b,1);
+						e = emin;
+						goto normal;
+						}
+					}
+				break;
+			  case Round_up:
+				if (!sign && (lostbits || (b->x[0] & 1))) {
+ incr_denorm:
+					multadd(b, 1, 2 MTa);
+					check_denorm = 1;
+					lostbits = 0;
+					goto emin_check;
+					}
+				break;
+			  case Round_down:
+				if (sign && (lostbits || (b->x[0] & 1)))
+					goto incr_denorm;
+				break;
+			  }
+			}
+#endif
 		if (lostbits)
 			lostbits = 1;
 		else if (k > 0)
 			lostbits = any_on(b,k);
+#ifdef IEEE_Arith
+		else if (check_denorm)
+			goto no_lostbits;
+#endif
 		if (x[k>>kshift] & 1 << (k & kmask))
 			lostbits |= 2;
+#ifdef IEEE_Arith
+ no_lostbits:
+#endif
 		nbits -= n;
 		rshift(b,n);
 		e = emin;
@@ -3020,16 +3068,9 @@ gethex( const char **sp, U *rvp, int rounding, int sign MTd)
 			k = b->wds;
 			b = increment(b MTa);
 			x = b->x;
-			if (denorm) {
-#if 0
-				if (nbits == Nbits - 1
-				 && x[nbits >> kshift] & 1 << (nbits & kmask))
-					denorm = 0; /* not currently used */
-#endif
-				}
-			else if (b->wds > k
+			if (!denorm && (b->wds > k
 			 || ((n = nbits & kmask) !=0
-			     && hi0bits(x[k-1]) < 32-n)) {
+			     && hi0bits(x[k-1]) < 32-n))) {
 				rshift(b,1);
 				if (++e > Emax)
 					goto ovfl;
@@ -3039,8 +3080,10 @@ gethex( const char **sp, U *rvp, int rounding, int sign MTd)
 #ifdef IEEE_Arith
 	if (denorm)
 		word0(rvp) = b->wds > 1 ? b->x[1] & ~0x100000 : 0;
-	else
+	else {
+ normal:
 		word0(rvp) = (b->x[1] & ~0x100000) | ((e + 0x3ff + 52) << 20);
+		}
 	word1(rvp) = b->x[0];
 #endif
 #ifdef IBM
@@ -3617,10 +3660,11 @@ strtod(const char *s00, char **se)
 				c = *++s;
 			if (c > '0' && c <= '9') {
 				L = c - '0';
-				s1 = s;
-				while((c = *++s) >= '0' && c <= '9')
-					L = 10*L + c - '0';
-				if (s - s1 > 8 || L > 19999)
+				while((c = *++s) >= '0' && c <= '9') {
+					if (L <= 19999)
+						L = 10*L + c - '0';
+					}
+				if (L > 19999)
 					/* Avoid confusion from exponents
 					 * so large that e might overflow.
 					 */
@@ -5416,7 +5460,8 @@ dtoa_r(double dd, int mode, int ndigits, int *decpt, int *sign, char **rve, char
 				&& (spec_case ? 4*res <= ulp : (2*res < ulp || dig & 1))) {
  ulp_reached:
 					if (ures < res
-					|| (ures == res && dig & 1))
+					|| (ures == res && dig & 1)
+					|| (dig == 9 && 2*ures <= ulp))
 						goto Roundup;
 					goto retc;
 					}
@@ -5662,7 +5707,7 @@ dtoa_r(double dd, int mode, int ndigits, int *decpt, int *sign, char **rve, char
 						}
 					}
 				if (!(zb & ures) && (ures-rb) << (1 - eulp) < ulp) {
-					if ((ures + rb) << (1 - eulp) < ulp)
+					if ((ures + rb) << (2 - eulp) < ulp)
 						goto Roundup;
 					goto Fast_failed1;
 					}

src/solvers/fpinit.c

@@ -49,6 +49,7 @@ int isatty_ASL; /* for use with "sw" under NT */
 #undef WIN32
 #define WIN32
 #else
+#define _GNU_SOURCE
 #include "fenv.h"
 #endif
 

src/solvers/getstub.c

@@ -616,6 +616,7 @@ badopt_ASL(Option_Info *oi)
 {
 	oi->n_badopts++;
 	oi->option_echo &= ~ASL_OI_echothis;
+	oi->need_nl = 0;
 	}
 
  char *

src/solvers/jac2dim.c

@@ -41,5 +41,5 @@ jac2dim_(const char *stub, fint *M, fint *N, fint *NO, fint *NZ,
 	if (cur_ASL)
 		return already_ASL("jacdim");
 	asl = ASL_alloc(ASL_read_pfgh);
-	return jac2dim_ASL(asl, stub, M, N, NO, NZ, MXROW, MXCOL, stub_len);
+	return jacpdim_ASL(asl, stub, M, N, NO, NZ, MXROW, MXCOL, stub_len);
 	}