399 demo complex cg using std complex

include/graphblas/algorithms/conjugate_gradient.hpp

@@ -23,8 +23,11 @@
 #define _H_GRB_ALGORITHMS_CONJUGATE_GRADIENT
 
 #include <cstdio>
+#include <complex>
 
 #include <graphblas.hpp>
+#include <graphblas/utils/iscomplex.hpp>
+
 
 namespace grb {
 
@@ -138,7 +141,8 @@ namespace grb {
 			class Minus = operators::subtract< IOType >,
 			class Divide = operators::divide< IOType >
 		>
-		grb::RC conjugate_gradient( grb::Vector< IOType > &x,
+		grb::RC conjugate_gradient(
+			grb::Vector< IOType > &x,
 			const grb::Matrix< NonzeroType > &A,
 			const grb::Vector< InputType > &b,
 			const size_t max_iterations,
@@ -190,7 +194,8 @@ namespace grb {
 			);
 
 			constexpr const Descriptor descr_dense = descr | descriptors::dense;
-			const ResidualType zero = ring.template getZero< ResidualType >();
+			const ResidualType zero_residual = ring.template getZero< ResidualType >();
+			const IOType zero = ring.template getZero< IOType >();
 			const size_t n = grb::ncols( A );
 
 			// dynamic checks
@@ -223,7 +228,7 @@ namespace grb {
 				}
 
 				// others
-				if( tol <= zero ) {
+				if( tol <= zero_residual ) {
 					std::cerr << "Error: tolerance input to CG must be strictly positive\n";
 					return ILLEGAL;
 				}
@@ -244,7 +249,7 @@ namespace grb {
 				assert( nnz( x ) == n );
 			}
 
-			ResidualType alpha, sigma, bnorm;
+			IOType sigma, bnorm, alpha, beta;
 
 			// temp = 0
 			grb::RC ret = grb::set( temp, 0 );
@@ -269,16 +274,33 @@ namespace grb {
 
 			// sigma = r' * r;
 			sigma = zero;
-			ret = ret ? ret : grb::dot< descr_dense >( sigma, r, r, ring );
+			if( grb::utils::is_complex< IOType >::value ) {
+				ret = ret ? ret : grb::eWiseLambda( [&temp,&r]( const size_t i ) {
+						temp[ i ] = grb::utils::is_complex< IOType >::conjugate( r[ i ] );
+					}, temp
+				);
+				ret = ret ? ret : grb::dot< descr_dense >( sigma, temp, r, ring );
+			} else {
+				ret = ret ? ret : grb::dot< descr_dense >( sigma, r, r, ring );
+			}
+
 			assert( ret == SUCCESS );
 
 			// bnorm = b' * b;
 			bnorm = zero;
-			ret = ret ? ret : grb::dot< descr_dense >( bnorm, b, b, ring );
+			if( grb::utils::is_complex< IOType >::value ) {
+				ret = ret ? ret : grb::eWiseLambda( [&temp,&b]( const size_t i ) {
+						temp[ i ] = grb::utils::is_complex< IOType >::conjugate( b[ i ] );
+					}, temp
+				);
+				ret = ret ? ret : grb::dot< descr_dense >( bnorm, temp, b, ring );
+			} else {
+				ret = ret ? ret : grb::dot< descr_dense >( bnorm, b, b, ring );
+			}
 			assert( ret == SUCCESS );
 
 			if( ret == SUCCESS ) {
-				tol *= sqrt( bnorm );
+				tol *= sqrt( grb::utils::is_complex< IOType >::modulus( bnorm ) );
 			}
 
 			size_t iter = 0;
@@ -292,13 +314,25 @@ namespace grb {
 				ret = ret ? ret : grb::mxv< descr_dense >( temp, A, u, ring );
 				assert( ret == SUCCESS );
 
-				// residual = u' * temp
-				residual = zero;
-				ret = ret ? ret : grb::dot< descr_dense >( residual, temp, u, ring );
+				// beta = u' * temp
+				beta = zero;
+				if( grb::utils::is_complex< IOType >::value ) {
+					ret = ret ? ret : grb::eWiseLambda( [&u]( const size_t i ) {
+							u[ i ] = grb::utils::is_complex< IOType >::conjugate( u[ i ] );
+						}, u
+					);
+				}
+				ret = ret ? ret : grb::dot< descr_dense >( beta, temp, u, ring );
+				if( grb::utils::is_complex< IOType >::value ) {
+					ret = ret ? ret : grb::eWiseLambda( [&u]( const size_t i ) {
+							u[ i ] = grb::utils::is_complex< IOType >::conjugate( u[ i ] );
+						}, u
+					);
+				}
 				assert( ret == SUCCESS );
 
-				// alpha = sigma / residual;
-				ret = ret ? ret : grb::apply( alpha, sigma, residual, divide );
+				// alpha = sigma / beta;
+				ret = ret ? ret : grb::apply( alpha, sigma, beta, divide );
 				assert( ret == SUCCESS );
 
 				// x = x + alpha * u;
@@ -314,9 +348,18 @@ namespace grb {
 				ret = ret ? ret : grb::foldl< descr_dense >( r, temp, minus );
 				assert( ret == SUCCESS );
 
-				// residual = r' * r;
-				residual = zero;
-				ret = ret ? ret : grb::dot< descr_dense >( residual, r, r, ring );
+				// beta = r' * r;
+				beta = zero;
+				if( grb::utils::is_complex< IOType >::value ) {
+					ret = ret ? ret : grb::eWiseLambda( [&temp,&r]( const size_t i ) {
+							temp[ i ] = grb::utils::is_complex< IOType >::conjugate( r[ i ] );
+						}, temp
+					);
+					ret = ret ? ret : grb::dot< descr_dense >( beta, temp, r, ring );
+				} else {
+					ret = ret ? ret : grb::dot< descr_dense >( beta, r, r, ring );
+				}
+				residual = grb::utils::is_complex< IOType >::modulus( beta );
 				assert( ret == SUCCESS );
 
 				if( ret == SUCCESS ) {
@@ -325,8 +368,8 @@ namespace grb {
 					}
 				}
 
-				// alpha = residual / sigma;
-				ret = ret ? ret : grb::apply( alpha, residual, sigma, divide );
+				// alpha = beta / sigma;
+				ret = ret ? ret : grb::apply( alpha, beta, sigma, divide );
 				assert( ret == SUCCESS );
 
 				// temp = r + alpha * u;
@@ -339,8 +382,7 @@ namespace grb {
 				// u = temp
 				std::swap( u, temp );
 
-				// sigma = residual;
-				sigma = residual;
+				sigma = beta;
 
 			} while( iter++ < max_iterations && ret == SUCCESS );
 

include/graphblas/utils.hpp

@@ -30,6 +30,7 @@
 #include <type_traits>
 
 #include <graphblas/descriptors.hpp>
+#include <graphblas/utils/iscomplex.hpp>
 
 
 namespace grb {
@@ -59,7 +60,12 @@ namespace grb {
 		 * @returns Whether a == b.
 		 */
 		template< typename T >
-		static bool equals( const T & a, const T & b, typename std::enable_if< ! std::is_floating_point< T >::value >::type * = NULL ) {
+		static bool equals(
+			const T &a, const T &b,
+			typename std::enable_if<
+				!std::is_floating_point< T >::value
+			>::type * = nullptr
+		) {
 			// simply do standard compare
 			return a == b;
 		}
@@ -83,18 +89,16 @@ namespace grb {
 		 */
 		template< typename T, typename U >
 		static bool equals( const T &a, const T &b, const U epsilons,
-			typename std::enable_if< std::is_floating_point< T >::value >::type * = NULL
+			typename std::enable_if<
+				std::is_floating_point< T >::value
+			>::type * = nullptr
 		) {
 			assert( epsilons >= 1 );
 
 			// if they are bit-wise equal, it's easy
 			if( a == b ) {
 #ifdef _DEBUG
- #ifndef _GRB_NO_STDIO
 				std::cout << "\t Bit-wise equal\n";
- #else
-				printf( "\t Bit-wise equal\n" );
- #endif
 #endif
 				return true;
 			}
@@ -106,22 +110,21 @@ namespace grb {
 			const T absPlus = absA + absB;
 
 			// find the effective epsilon
-			const T eps = static_cast< T >( epsilons ) * std::numeric_limits< T >::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.
+			// 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 ) {
 #ifdef _DEBUG
- #ifndef _GRB_NO_STDIO
 				std::cout << "\t Zero or close to zero difference\n";
- #else
-				printf( "\t Zero or close to zero difference\n" );
- #endif
 #endif
 				return absDiff < eps * min;
 			}
@@ -131,33 +134,21 @@ namespace grb {
 			if( absA > absB ) {
 				if( absB > max - absA ) {
 #ifdef _DEBUG
- #ifndef _GRB_NO_STDIO
 					std::cout << "\t Normalising absolute difference by max (I)\n";
- #else
-					printf( "\t Normalising absolute difference by max (I)\n" );
- #endif
 #endif
 					return absDiff / max < eps;
 				}
 			} else {
 				if( absA > max - absB ) {
 #ifdef _DEBUG
- #ifndef _GRB_NO_STDIO
 					std::cout << "\t Normalising absolute difference by max (II)\n";
- #else
-					printf( "\t Normalising absolute difference by max (II)\n" );
- #endif
 #endif
 					return absDiff / max < eps;
 				}
 			}
 			// use of relative error should be safe
 #ifdef _DEBUG
- #ifndef _GRB_NO_STDIO
 			std::cout << "\t Using relative error\n";
- #else
-			printf( "\t Using relative error\n" );
- #endif
 #endif
 			return absDiff / absPlus < eps;
 		}
@@ -176,7 +167,7 @@ namespace grb {
 		 * return value is a constexpr. (This was fixed in C++14.)
 		 */
 		template< typename T >
-		constexpr const T & static_min( const T & a, const T & b ) {
+		constexpr const T & static_min( const T &a, const T &b ) {
 			return a < b ? a : b;
 		}
 
@@ -186,19 +177,19 @@ namespace grb {
 		 */
 		template< typename T >
 		class SizeOf {
-		public:
-			/**
-			 * If \a T is <tt>void</tt>, this value equals 0 and
-			 * equal to <tt>sizeof(T)</tt> otherwise.
-			 */
-			static constexpr const size_t value = sizeof( T );
+			public:
+				/**
+				 * If \a T is <tt>void</tt>, this value equals 0 and
+				 * equal to <tt>sizeof(T)</tt> otherwise.
+				 */
+				static constexpr const size_t value = sizeof( T );
 		};
 
-		// void-specialisation of the above class
+		/** \internal void-specialisation of the above class */
 		template<>
 		class SizeOf< void > {
-		public:
-			static constexpr const size_t value = 0;
+			public:
+				static constexpr const size_t value = 0;
 		};
 
 		/**
@@ -237,7 +228,11 @@ namespace grb {
 		 * compile.
 		 */
 		template< Descriptor descriptor, typename T >
-		static bool interpretMask( const bool & assigned, const T * const val, const size_t offset ) {
+		static bool interpretMask(
+			const bool &assigned,
+			const T * const val,
+			const size_t offset
+		) {
 			// set default mask to false
 			bool ret = false;
 			// if we request a structural mask, decide only on passed assigned variable
@@ -252,20 +247,52 @@ namespace grb {
 			}
 			// check whether we should return the inverted value
 			if( descriptor & descriptors::invert_mask ) {
-				return ! ret;
+				return !ret;
+			} else {
+				return ret;
+			}
+		}
+
+		/** Specialisation for complex-valued masks */
+		template< Descriptor descriptor, typename T >
+		static bool interpretMask(
+				const bool &assigned,
+				const std::complex<T> * const val,
+				const size_t offset
+		) {
+			// set default mask to false
+			bool ret = false;
+			// if we request a structural mask, decide only on passed assigned variable
+			if( descriptor & descriptors::structural ) {
+				ret = assigned;
+			} else {
+				// if based on value, if there is a value, cast it to bool
+				if( assigned ) {
+					ret = static_cast< bool >( real( val [ offset ] ) ) ||
+					       static_cast< bool >( imag( val [ offset ] ) );
+				}
+				// otherwise there is no value and false is assumed
+			}
+			// check whether we should return the inverted value
+			if( descriptor & descriptors::invert_mask ) {
+				return !ret;
 			} else {
 				return ret;
 			}
 		}
 
 		/** Specialisation for void-valued masks */
 		template< Descriptor descriptor >
-		static bool interpretMask( const bool & assigned, const void * const, const size_t ) {
+		static bool interpretMask(
+			const bool &assigned,
+			const void * const,
+			const size_t
+		) {
 			// set default mask to false
 			bool ret = assigned;
 			// check whether we should return the inverted value
 			if( descriptor & descriptors::invert_mask ) {
-				return ! ret;
+				return !ret;
 			} else {
 				return ret;
 			}
@@ -276,3 +303,4 @@ namespace grb {
 } // namespace grb
 
 #endif
+