342 execution issue with custom preconditioning function in conjugate gradient non blocking backend (#344)

This PCG algorithm allows for user-defined preconditioners to be supplied by the user. This functionality, in turn, was then used by the solver and KML_Solver transition paths to offer ready-made libraries for which users could supply their own preconditioners via function pointers. Those preconditioners would operate on standard raw C arrays, thus making use of the ALP native interface -- however, the PCG algorithm in ALP/GraphBLAS implicitly assumed the preconditioner action was implemented as an ALP function; i.e., not using the native interface.

This mismatch causes issues in combination with the nonblocking backend, which guarantees valid nonblocking execution for pure ALP programs, but in combination with native code requires calls to `grb::wait`.

This MR introduces a template argument to indicate whether the preconditioner may make use of native calls. This is done by extending the already-existing Boolean argument that previously indicated whether to apply preconditioning-- this is now become an `int` rather than `bool`, with `0` indicating no preconditioning, `1` indicating non-native (ALP-based) user-defined preconditioning, and `2` indicating possibly-native user-defined preconditioning. Only in the latter case does this MR furthermore activate the new calls to `grb::wait` that are necessary for native preconditioning functions.

An alternative design has that the authors of user-defined preconditioners should call `grb::wait` before taking and operating on native views of ALP/GraphBLAS vectors. In this design, the reported issue #342 would in fact indicate a bug in the related transition libraries. Comparing the implemented solution with this alternative design has that more care should be taken by those supplying user-defined preconditioners, and we instead prefer to make the life of any ALP user easier; hence the present solution was taken instead.

The bug that this MR fixes affected the PCG transition paths only; i.e., users of `solver.h` and/or `kml_iss.h`, only when using custom preconditioners, and only non-diagonal preconditioner actions.

Author: anyzelman

include/graphblas/algorithms/conjugate_gradient.hpp

@@ -45,12 +45,19 @@ namespace grb {
 		 * positive definite.
 		 *
 		 * @tparam descr          The user descriptor
-		 * @tparam preconditioned Whether to apply any given preconditioners.
+		 * @tparam preconditioned Nonzero to apply any given preconditioners, zero to
+		 *                        not apply preconditioning and ignore related
+		 *                        arguments, two if the preconditioner action may be
+		 *                        native.
 		 *
-		 * \note The default value for \a preconditioned is <tt>true</tt> and it is
+		 * \note The default value for \a preconditioned is <tt>1</tt> and it is
 		 *       normally not necessary to override this value: if wishing to call a
 		 *       non-preconditioned CG, please see #conjugate_gradient.
 		 *
+		 * \note Native in the above means the preconditioner action may not consist
+		 *       of pure ALP code. In this case a nonblocking implementation should
+		 *       take care to not attempt to fuse across the preconditioner action.
+		 *
 		 * @tparam IOType         The input/output vector nonzero type
 		 * @tparam ResidualType   The type of the residual
 		 * @tparam NonzeroType    The matrix nonzero type
@@ -98,8 +105,8 @@ namespace grb {
 		 * achieved by constructing \a Minv so that the condition number of
 		 * \f$ M^{-1}A \f$ is much smaller than that of \f$ A \f$.
 		 *
-		 * @param[in]     Minv           The preconditioner action if
-		 *                               \a preconditioned equals <tt>true</tt>.
+		 * @param[in]     Minv           The preconditioner action, used only if
+		 *                               \a preconditioned is nonzero.
 		 *
 		 * If \a A is \f$ n \times n \f$, then \a x and \a b must have matching length
 		 * \f$ n \f$. The vector \a x furthermore must have a capacity of \f$ n \f$.
@@ -123,7 +130,7 @@ namespace grb {
 		 * @param[in,out] temp           A temporary vector of the same size as \a x.
 		 * @param[in,out] temp_precond   A temporary vector of the same size as \a x.
 		 *
-		 * \note If \a preconditioned is <tt>false</tt>, then both \a Minv and
+		 * \note If \a preconditioned is <tt>0</tt>, then both \a Minv and
 		 *       \a temp_precond are ignored. In this case, \a temp_precond need not
 		 *       have the same length as \a x, nor need it have full capacity.
 		 *
@@ -175,7 +182,7 @@ namespace grb {
 		 */
 		template<
 			Descriptor descr = descriptors::no_operation,
-			bool preconditioned = true,
+			int preconditioned = 1,
 			typename IOType,
 			typename ResidualType,
 			typename NonzeroType,
@@ -344,6 +351,9 @@ namespace grb {
 			// z = M^-1r
 			if( preconditioned ) {
 				ret = ret ? ret : grb::set( z, 0 ); // also ensures z is dense, henceforth
+				if( preconditioned == 2 ) {
+					ret = ret ? ret : grb::wait( z, r );
+				}
 				ret = ret ? ret : Minv( z, r );
 			} // else, z equals r (by reference)
 
@@ -443,6 +453,9 @@ namespace grb {
 				// beta = r' * z
 				if( preconditioned ) {
 					beta = zero;
+					if( preconditioned == 2 ) {
+						ret = ret ? ret : grb::wait( z, r );
+					}
 					ret = ret ? ret : Minv( z, r ); assert( ret == grb::SUCCESS );
 					ret = ret ? ret : grb::dot< descr_dense >(
 							beta,
@@ -572,7 +585,7 @@ namespace grb {
 			grb::Vector< IOType, backend > dummy_buffer( 0 );
 
 			// call PCG with preconditioning disabled
-			return preconditioned_conjugate_gradient< descr, false >(
+			return preconditioned_conjugate_gradient< descr, 0 >(
 					x, A, b,
 					dummy_preconditioner,
 					max_iterations, tol,

src/transition/solver.cpp

@@ -261,7 +261,7 @@ class CG_Data {
 					workspace[ 0 ], workspace[ 1 ], workspace[ 2 ]
 				);
 			} else {
-				return grb::algorithms::preconditioned_conjugate_gradient< descr >(
+				return grb::algorithms::preconditioned_conjugate_gradient< descr, 2 >(
 					x, matrix, b,
 					alpified_preconditioner(),
 					max_iter, tolerance,