Add early exit in preconditioned CG when initial residual is below tolerance (#350)

Previously, the (preconditioned) Conjugate Gradient (CG) algorithm applies the preconditioner before checking whether the current residual already met the required tolerance. The changes from this MR allows the CG to return immediately when the initial residual is small enough, avoiding the unnecessary cost of a preconditioner call. This change can improve performance in cases where the initial guess is already close to the solution. Thanks to @ChristosMatzoros for proposing this enhancement!

Author: ChristosMatzoros

NOTICE

@@ -21,7 +21,7 @@ to Huawei Technologies Co., Ltd. or one of its subsidiaries:
 
  - Aristeidis Mastoras, Huawei Technologies Switzerland AG; 2020-current.
 
- - Alberto Scolari, Huawei Technologies Switzerland AG; 2021-current.
+ - Alberto Scolari, Huawei Technologies Switzerland AG; 2021-2025.
 
  - Verner Vlacic, Huawei Technologies Switzerland AG; 2021-current.
 
@@ -34,6 +34,8 @@ to Huawei Technologies Co., Ltd. or one of its subsidiaries:
  - Denis Jelovina, Huawei Technologies Switzerland AG; 2022-current.
 
  - Benjamin Lozes, Huawei Technologies Switzerland AG; 2023.
+ 
+ - Christos Konstantinos Matzoros, Huawei Technologies Switzerland AG; 2024-current.
 
 The experimental banshee backend has been developed in collaboration with
 Prof. Luca Benini at ETH Zuerich and his group. In particular this backend

include/graphblas/algorithms/conjugate_gradient.hpp

@@ -348,6 +348,37 @@ namespace grb {
 			ret = ret ? ret : grb::foldl< descr_dense >( r, temp, minus );
 			assert( ret == grb::SUCCESS );
 
+			// bnorm = b' * b;
+			bnorm = zero;
+			ret = ret ? ret : grb::dot< descr_dense >(
+					bnorm,
+					b, b,
+					ring.getAdditiveMonoid(),
+					grb::operators::conjugate_left_mul< IOType >()
+				);
+			assert( ret == grb::SUCCESS );
+
+			// get effective tolerance
+			if( ret == grb::SUCCESS ) {
+				tol *= std::sqrt( grb::utils::is_complex< IOType >::modulus( bnorm ) );
+			}
+
+			// get residual
+			alpha = zero;
+			ret = ret ? ret : grb::dot< descr_dense >(
+					alpha,
+					r, r,
+					ring.getAdditiveMonoid(),
+					grb::operators::conjugate_left_mul< IOType >()
+				);
+			assert( ret == grb::SUCCESS );
+			residual = grb::utils::is_complex< IOType >::modulus( alpha );
+
+			// check residual for early exit
+			if( ret == grb::SUCCESS ) {
+				if( sqrt( residual ) < tol ) { return ret; }
+			}
+
 			// z = M^-1r
 			if( preconditioned ) {
 				ret = ret ? ret : grb::set( z, 0 ); // also ensures z is dense, henceforth
@@ -374,19 +405,8 @@ namespace grb {
 
 			assert( ret == grb::SUCCESS );
 
-			// bnorm = b' * b;
-			bnorm = zero;
-			ret = ret ? ret : grb::dot< descr_dense >(
-				bnorm,
-				b, b,
-				ring.getAdditiveMonoid(),
-				grb::operators::conjugate_left_mul< IOType >() );
-			assert( ret == grb::SUCCESS );
-
-			// get effective tolerance and exit on any error during prelude
-			if( ret == grb::SUCCESS ) {
-				tol *= std::sqrt( grb::utils::is_complex< IOType >::modulus( bnorm ) );
-			} else {
+			// exit on any error during prelude
+			if( ret != grb::SUCCESS ) {
 				std::cerr << "Warning: preconditioned CG caught error during prelude ("
 					<< grb::toString( ret ) << ")\n";
 				return ret;