464 dense spmv test does not work for ecology2 mtx

include/graphblas/algorithms/bicgstab.hpp

@@ -314,7 +314,7 @@ namespace grb {
 #ifdef _DEBUG
 				std::cout << "\t\t rho = " << rho << "\n";
 #endif
-				if( ret == SUCCESS && utils::equals( rho, zero, 2*n-1 ) ) {
+				if( ret == SUCCESS && rho == zero ) {
 					std::cerr << "Error: BiCGstab detects r at iteration " << iterations <<
 						" is orthogonal to r-hat\n";
 					return FAILED;
@@ -346,7 +346,7 @@ namespace grb {
 				// alpha = rho / (rhat, v)
 				alpha = zero;
 				ret = ret ? ret : dot< dense_descr >( alpha, rhat, v, semiring );
-				if( utils::equals( alpha, zero, 2*n-1 ) ) {
+				if( alpha == zero ) {
 					std::cerr << "Error: BiCGstab detects rhat is orthogonal to v=Ap "
 						<< "at iteration " << iterations << ".\n";
 					return FAILED;
@@ -391,7 +391,7 @@ namespace grb {
 #ifdef _DEBUG
 				std::cout << "\t\t (t, s) = " << temp << "\n";
 #endif
-				if( ret == SUCCESS && utils::equals( rho, zero, 2*n-1 ) ) {
+				if( ret == SUCCESS && temp == zero ) {
 					std::cerr << "Error: BiCGstab detects As at iteration " << iterations <<
 						" is orthogonal to s\n";
 					return FAILED;

include/graphblas/utils.hpp

@@ -71,29 +71,83 @@ namespace grb {
 		}
 
 		/**
-		 * Checks whether two floating point values are equal.
+		 * Checks whether two floating point values, \a a and \a b, are equal.
 		 *
-		 * This function takes into account round-off errors due to machine precision.
+		 * The comparison is relative in the following sense.
 		 *
-		 * \warning This does not take into account accumulated numerical errors
-		 *          due to previous operations on the given values.
+		 * It is mandatory to supply an integer, \a epsilons, that indicates a
+		 * relative error. Here, \a epsilon represents the number of errors
+		 * accumulated during their computation, assuming that all arithmetic that
+		 * produced the two floating point arguments to this function are bounded by
+		 * the magnitudes of \a a and \a b.
 		 *
-		 * @\tparam T The numerical type.
+		 * Intuitively, one may take \a epsilons as the sum of the number of
+		 * operations that produced \a a and \a b, if the above assumption holds.
 		 *
-		 * @param a One of the two values to comare against.
-		 * @param b One of the two values to comare against.
-		 * @param epsilons How many floating point errors may have accumulated.
-		 *                 Should be chosen larger or equal to one.
+		 * The resulting bound can be tightened if the magnitudes encountered during
+		 * their computation are much smaller, and are (likely) too tight if those
+		 * magnitudes were much larger instead. In such cases, one should obtain or
+		 * compute a more appropriate error bound, and check whether the difference
+		 * is within such a more appropriate, absolute, error bound.
 		 *
-		 * @returns Whether a == b.
+		 * \warning When comparing for equality within a known absolute error bound,
+		 *          <b>do not this function</b>.
+		 *
+		 * \note If such a function is desired, please submit an issue on GitHub or
+		 *       Gitee. The absolute error bound can be provided again via a
+		 *       parameter \a epsilons, but then of a floating-point type instead of
+		 *       the current integer variant.
+		 *
+		 * @tparam T The numerical type.
+		 * @tparam U The integer type used for \a epsilons.
+		 *
+		 * @param[in] a        One of the two values to compare.
+		 * @param[in] b        One of the two values to compare.
+		 * @param[in] epsilons How many floating point errors may have accumulated;
+		 *                     must be chosen larger or equal to one.
+		 *
+		 * This function automatically adapts to the floating-point type used, and
+		 * takes into account the border cases where one or more of \a a and \a b may
+		 * be zero or subnormal. It also guards against overflow of the normalisation
+		 * strategy employed in its implementation.
+		 *
+		 * If one of \a a or \a b is zero, then an absolute acceptable error bound is
+		 * computed as \a epsilons times \f$ \epsilon \f$, where the latter is the
+		 * machine precision.
+		 *
+		 * \note This behaviour is consistent with the intuitive usage of this
+		 *       function as described above, assuming that the computations that led
+		 *       to zero involved numbers of zero orders of magnitude.
+		 *
+		 * \warning Typically it is hence \em not correct to rely on this function to
+		 *          compare some value \a a to zero by passing a zero \a b explicitly,
+		 *          and vice versa-- the magnitude of the values used while computing
+		 *          the nonzero matter.
+		 *
+		 * \note If a version of this function with an explicit relative error bound
+		 *       is desired, please submit an issue on GitHub or Gitee.
+		 *
+		 * @returns Whether \a a equals \a b, taking into account numerical drift
+		 *          within the effective range derived from \a epsilons and the
+		 *          magnitudes of \a a and \a b.
+		 *
+		 * This utility function is chiefly used by the ALP/GraphBLAS unit, smoke, and
+		 * performance tests.
 		 */
 		template< typename T, typename U >
-		static bool equals( const T &a, const T &b, const U epsilons,
+		static bool equals(
+			const T &a, const T &b,
+			const U epsilons,
 			typename std::enable_if<
-				std::is_floating_point< T >::value
+				std::is_floating_point< T >::value &&
+				std::is_integral< U >::value
 			>::type * = nullptr
 		) {
 			assert( epsilons >= 1 );
+#ifdef _DEBUG
+			std::cout << "Comparing " << a << " to " << b << " with relative tolerance "
+				<< "of " << epsilons << "\n";
+#endif
 
 			// if they are bit-wise equal, it's easy
 			if( a == b ) {
@@ -103,54 +157,68 @@ namespace grb {
 				return true;
 			}
 
+			// find the effective epsilon and absolute difference
+			const T eps = static_cast< T >( epsilons ) *
+				std::numeric_limits< T >::epsilon();
+			const T absDiff = fabs( a - b );
+
+			// If either a or b is zero, then we use the relative epsilons as an absolute
+			// error on the nonzero argument
+			if( a == 0 || b == 0 ) {
+				assert( a != 0 || b != 0 );
+#ifdef _DEBUG
+				std::cout << "\t One of the arguments is zero\n";
+#endif
+				return absDiff < eps;
+			}
+
+
 			// if not, we need to look at the absolute values
 			const T absA = fabs( a );
 			const T absB = fabs( b );
-			const T absDiff = fabs( a - b );
 			const T absPlus = absA + absB;
 
-			// find the effective epsilon
-			const T eps = static_cast< T >( epsilons ) *
-				std::numeric_limits< T >::epsilon();
-
 			// find the minimum and maximum *normal* values.
 			const T min = std::numeric_limits< T >::min();
-			const T max = std::numeric_limits< T >::max();
 
-			// if the difference is a subnormal number, it should be smaller than
-			// machine epsilon times min;
-			// if this is not the case, then we cannot safely conclude anything from this
-			// small a difference.
-			// The same is true if a or b are zero.
-			if( a == 0 || b == 0 || absPlus < min ) {
+			// If the combined magnitudes of a or b are subnormal, then scale by the
+			// smallest normal rather than the combined magnitude
+			if( absPlus < min ) {
 #ifdef _DEBUG
-				std::cout << "\t Zero or close to zero difference\n";
+				std::cout << "\t Subnormal comparison requested. Scaling the tolerance by "
+					<< "the smallest normal number rather than to the subnormal absolute "
+					<< "difference\n";
 #endif
-				return absDiff < eps * min;
+				return absDiff / min < eps;
 			}
 
-			// we wish to normalise absDiff by (absA + absB),
-			// However, absA + absB might overflow.
+			// we wish to normalise absDiff by (absA + absB), however, absA + absB might
+			// overflow. If it does, we normalise by the smallest of absA and absB
+			// instead of attempting to average.
+			const T max = std::numeric_limits< T >::max();
 			if( absA > absB ) {
 				if( absB > max - absA ) {
 #ifdef _DEBUG
-					std::cout << "\t Normalising absolute difference by max (I)\n";
+					std::cout << "\t Normalising absolute difference by smallest magnitude "
+						<< "(I)\n";
 #endif
-					return absDiff / max < eps;
+					return absDiff / absB < eps;
 				}
 			} else {
 				if( absA > max - absB ) {
 #ifdef _DEBUG
-					std::cout << "\t Normalising absolute difference by max (II)\n";
+					std::cout << "\t Normalising absolute difference by smallest magnitude "
+						<< "(II)\n";
 #endif
-					return absDiff / max < eps;
+					return absDiff / absA < eps;
 				}
 			}
-			// use of relative error should be safe
+
+			// at this point, the use of relative error vs. 0.5*absPlus should be safe
 #ifdef _DEBUG
 			std::cout << "\t Using relative error\n";
 #endif
-			return absDiff / absPlus < eps;
+			return (static_cast< T >(2) * (absDiff / absPlus)) < eps;
 		}
 
 		/**