Non-square matrix CCS assignment bugfix

Author: byjtew

include/graphblas/reference/blas3.hpp

@@ -1062,11 +1062,19 @@ namespace grb {
 
 			const size_t nzc = capacity( C );
 
+#ifdef _H_GRB_REFERENCE_OMP_BLAS3
+#pragma omp parallel for simd
+#endif
+			for( size_t i = 0; i <= n_A; i++ ) {
+				C_ccs_raw.col_start[ i ] = 0;
+			}
+
 			C_crs_raw.col_start[ 0 ] = 0;
-			C_ccs_raw.col_start[ 0 ] = 0;
 			// Prefix sum computation into L.CRS.col_start
 #ifdef _H_GRB_REFERENCE_OMP_BLAS3
-#pragma omp parallel for default( none ) shared( B_raw, A_raw, C_crs_raw, std::cout ) firstprivate( m_A )
+#pragma omp parallel for default( none ) \
+			shared( B_raw, A_raw, C_crs_raw, C_ccs_raw, std::cout ) \
+			firstprivate( m_A )
 #endif
 			for( size_t i = 0; i < m_A; i++ ) {
 				auto B_k = B_raw.col_start[ i ];
@@ -1086,48 +1094,65 @@ namespace grb {
 					}
 					if( B_raw.row_index[ B_k ] == j ) {
 						cumul += 1;
+						C_ccs_raw.col_start[ j + 1 ] += 1;
 					}
 				}
 				C_crs_raw.col_start[ i + 1 ] = cumul;
 			}
 
 #ifdef _DEBUG
 			// Print the CRS prefix sum
-			std::cout << "CRS prefix sum: ";
+			std::cout << "before nCRS prefix sum: ";
 			for( size_t i = 0; i <= m_A; i++ ) {
 				std::cout << C_crs_raw.col_start[ i ] << " ";
 			}
 			std::cout << "\n";
+			// Print the CCS prefix sum
+			std::cout << "before nCCS prefix sum: ";
+			for( size_t i = 0; i <= n_A; i++ ) {
+				std::cout << C_ccs_raw.col_start[ i ] << " ";
+			}
+			std::cout << "\n";
 #endif
 
 			// Apply the prefix sum
 			for( size_t i = 1; i <= m_A; i++ ) {
 				C_crs_raw.col_start[ i ] += C_crs_raw.col_start[ i - 1 ];
-				C_ccs_raw.col_start[ i ] = C_crs_raw.col_start[ i ];
 			}
+			for ( size_t i = 1; i <= n_A; i++ ) {
+				C_ccs_raw.col_start[ i ] += C_ccs_raw.col_start[ i - 1 ];
+			}
+
+#ifdef _DEBUG
+			// Print the CRS prefix sum
+			std::cout << "after nCRS prefix sum: ";
+			for( size_t i = 0; i <= m_A; i++ ) {
+				std::cout << C_crs_raw.col_start[ i ] << " ";
+			}
+			std::cout << "\n";
+			// Print the CCS prefix sum
+			std::cout << "after nCCS prefix sum: ";
+			for( size_t i = 0; i <= n_A; i++ ) {
+				std::cout << C_ccs_raw.col_start[ i ] << " ";
+			}
+			std::cout << "\n";
+#endif
 
 			// Check if the number of nonzeros is greater than the capacity
-			if( C_crs_raw.col_start[ m_A ] > nzc || C_ccs_raw.col_start[ m_A ] > nzc ) {
+			if( C_crs_raw.col_start[ m_A ] > nzc || C_ccs_raw.col_start[ n_A ] > nzc ) {
 #ifdef _DEBUG
 				std::cout << "Insufficient capacity detected for requested operation.\n"
-						  << "Requested " << C_crs_raw.col_start[ m_A ] << " nonzeros"
+						  << "Requested " << C_ccs_raw.col_start[ m_A ] << " nonzeros"
 						  << " but capacity is " << nzc << "\n";
 #endif
 				return MISMATCH;
 			}
 
-#ifdef _H_GRB_REFERENCE_OMP_BLAS3
-#pragma omp parallel for simd
-#endif
-			for( size_t i = 0; i < m_A; i++ ) {
-				C_crs_raw.row_index[ i ] = C_ccs_raw.row_index[ i ] = 0;
-			}
-
 			RC local_rc = rc;
 #ifdef _H_GRB_REFERENCE_OMP_BLAS3
 #pragma omp parallel default( none ) \
-			shared( C_ccs_raw, C_crs_raw, A_raw, B_raw, rc, std::cout ) \
-			firstprivate( local_rc, m_A, oper, A_identity, B_identity )
+				shared( C_ccs_raw, C_crs_raw, A_raw, B_raw, rc, std::cout ) \
+				firstprivate( local_rc, m_A, oper, A_identity, B_identity )
 #endif
 			{
 				size_t start_row = 0;
@@ -1144,6 +1169,7 @@ namespace grb {
 					const auto A_k_end = A_raw.col_start[ i + 1 ];
 					for( auto A_k = A_k_start; A_k < A_k_end; ++A_k ) {
 						const auto j = A_raw.row_index[ A_k ];
+
 						while(  B_k < B_raw.col_start[ i + 1 ]
 							    && B_raw.row_index[ B_k ] > j
 						) {
@@ -1165,16 +1191,17 @@ namespace grb {
 
 						C_crs_raw.row_index[ C_k ] = j;
 						C_crs_raw.setValue( C_k, c_val );
-						C_ccs_raw.row_index[ C_k ] = i;
-						C_ccs_raw.setValue( C_k, c_val );
+
+						C_ccs_raw.row_index[ C_ccs_raw.col_start[ j ] ] = i;
+						C_ccs_raw.setValue( C_ccs_raw.col_start[ j ], c_val );
 #ifdef _DEBUG
 						std::cout << "A( " + std::to_string( i ) + ";"
 									+ std::to_string( j ) + " ) = "
 									+ std::to_string( a_val ) + "\n";
 						std::cout << "B( " + std::to_string( i ) + ";"
 									+ std::to_string( j ) + " ) = "
 									+ std::to_string( b_val ) + "\n";
-						std::cout << "C.crs( " + std::to_string( i ) + ";"
+						std::cout << "C( " + std::to_string( i ) + ";"
 									+ std::to_string( j ) + " ) = "
 									+ std::to_string( c_val ) + "\n";
 #endif

tests/unit/eWiseApply_matrix.cpp

@@ -23,23 +23,24 @@
 
 using namespace grb;
 
-void grb_program( const int &, grb::RC &rc ) {
+void grb_program( const size_t &n, grb::RC &rc ) {
 
 	// large non-square mixed-domain matrix check
 	{
-		grb::Matrix< char > A( 10000000, 2000000 );
-		grb::Matrix< float > B( 10000000, 2000000 );
-		grb::Matrix< size_t > C( 10000000, 2000000 );
-		size_t * I = new size_t[ 2000000 ];
-		size_t * J = new size_t[ 2000000 ];
-		char * V = new char[ 2000000 ];
-		for( size_t k = 0; k < 2000000; ++k ) {
-			I[ k ] = J[ k ] = k;
+		grb::Matrix< char > A( n, 2*n );
+		grb::Matrix< float > B( n, 2*n );
+		grb::Matrix< size_t > C( n, 2*n );
+		size_t * I = new size_t[ n ];
+		size_t * J = new size_t[ n ];
+		char * V = new char[ n ];
+		for( size_t k = 0; k < n; ++k ) {
+			I[ k ] = k;
+			J[ k ] = k+n;
 			V[ k ] = 2;
 		}
-		rc = grb::buildMatrixUnique( A, I, J, V, 2000000, SEQUENTIAL );
-		rc = rc ? rc : grb::buildMatrixUnique( B, I, J, V, 2000000, SEQUENTIAL );
-		rc = rc ? rc : grb::buildMatrixUnique( C, I, J, V, 2000000, SEQUENTIAL );
+		rc = grb::buildMatrixUnique( A, I, J, V, n, SEQUENTIAL );
+		rc = rc ? rc : grb::buildMatrixUnique( B, I, J, V, n, SEQUENTIAL );
+		rc = rc ? rc : grb::buildMatrixUnique( C, I, J, V, n, SEQUENTIAL );
 		rc = rc ? rc : grb::eWiseApply( C, A, B,
 			grb::operators::add< float, size_t, char >(), RESIZE );
 		rc = rc ? rc : grb::eWiseApply( C, A, B,
@@ -49,13 +50,14 @@ void grb_program( const int &, grb::RC &rc ) {
 				<< "mixed-domain matrix check\n";
 			return;
 		}
+		
 		for( const auto &triple : C ) {
-			const size_t &i = triple.first.first;
-			const size_t &j = triple.first.second;
-			const size_t &v = triple.second;
-			if( i != j ) {
-				std::cout << "Unexpected entry at position ( " << i << ", " << j << " ) "
-					<< "-- only expected entries on the diagonal\n";
+			const auto &i = triple.first.first;
+			const auto &j = triple.first.second;
+			const auto &v = triple.second;
+			if( j != i+n ) {
+				std::cout << "Unexpected entry at position ( " << i << ", " << i+n << " ) "
+					<< "-- only expected entries on the n-th diagonal\n";
 				rc = FAILED;
 			}
 			if( v != 4 ) {
@@ -74,15 +76,14 @@ void grb_program( const int &, grb::RC &rc ) {
 
 int main( int argc, char ** argv ) {
 	// defaults
-	bool printUsage = false;
-	int input = 0; // unused
+	size_t input = 1000; // unused
 
 	// error checking
 	if( argc > 1 ) {
-		printUsage = true;
+		input = std::strtoul( argv[ 1 ], nullptr, 10 );
 	}
-	if( printUsage ) {
-		std::cerr << "Usage: " << argv[ 0 ] << "\n";
+	if( argc > 2 ) {
+		std::cerr << "Usage: " << argv[ 0 ] << "[n]\n";
 		return 1;
 	}