diff --git a/include/graphblas/base/blas3.hpp b/include/graphblas/base/blas3.hpp index c22bfba71..32c8d25ac 100644 --- a/include/graphblas/base/blas3.hpp +++ b/include/graphblas/base/blas3.hpp @@ -442,6 +442,96 @@ namespace grb { return ret == SUCCESS ? UNSUPPORTED : ret; } + /** + * Masked outer product of two vectors. Assuming vectors \a u and \a v are + * oriented column-wise, the result matrix \a A will contain \f$ uv^T \f$, + * masked to non-zero values from \a mask. This is an out-of-place + * function and will be updated soon to be in-place instead. + * + * \internal Implemented via mxm as a multiplication of a column vector with + * a row vector, while following the given mask. + * + * @tparam descr The descriptor to be used. Optional; the default is + * #grb::descriptors::no_operation. + * @tparam Operator The operator to use. + * @tparam InputType1 The value type of the left-hand vector. + * @tparam InputType2 The value type of the right-hand vector. + * @tparam MaskType The value type of the mask matrix. + * @tparam OutputType The value type of the ouput matrix + * + * @param[out] A The output matrix. + * @param[out] mask The mask matrix. + * @param[in] u The left-hand input vector. + * @param[in] v The right-hand input vector. + * @param[in] op The operator. + * @param[in] phase The #grb::Phase the call should execute. Optional; the + * default parameter is #grb::EXECUTE. + * + * @return #grb::SUCCESS On successful completion of this call. + * @return #grb::MISMATCH Whenever the dimensions of the inputs and/or outputs + * do not match. All input data containers are left + * untouched if this exit code is returned; it will be + * be as though this call was never made. + * @return #grb::FAILED If \a phase is #grb::EXECUTE, indicates that the + * capacity of \a A was insufficient. The output + * matrix \a A is cleared, and the call to this function + * has no further effects. + * @return #grb::OUTOFMEM If \a phase is #grb::RESIZE, indicates an + * out-of-memory exception. The call to this function + * shall have no other effects beyond returning this + * error code; the previous state of \a A is retained. + * @return #grb::PANIC A general unmitigable error has been encountered. If + * returned, ALP enters an undefined state and the user + * program is encouraged to exit as quickly as possible. + * + * \par Performance semantics + * + * Each backend must define performance semantics for this primitive. + * + * @see perfSemantics + */ + template< + Descriptor descr = descriptors::no_operation, + class Operator, + typename InputType1, typename InputType2, + typename OutputType, typename MaskType, + typename Coords, + typename RIT, typename CIT, typename NIT, + Backend backend + > + RC outer( + Matrix< OutputType, backend, RIT, CIT, NIT > &A, + const Matrix< MaskType, backend, RIT, CIT, NIT > &mask, + const Vector< InputType1, backend, Coords > &u, + const Vector< InputType2, backend, Coords > &v, + const Operator &op = Operator(), + const Phase &phase = EXECUTE, + const typename std::enable_if< + grb::is_operator< Operator >::value && + !grb::is_object< InputType1 >::value && + !grb::is_object< InputType2 >::value && + !grb::is_object< MaskType >::value && + !grb::is_object< OutputType >::value, + void >::type * const = nullptr + ) { + (void) A; + (void) mask; + (void) u; + (void) v; + (void) op; + (void) phase; +#ifdef _DEBUG + std::cerr << "Selected backend does not implement grb::outer\n"; +#endif +#ifndef NDEBUG + const bool selected_backend_does_not_support_masked_outer = false; + assert( selected_backend_does_not_support_masked_outer ); +#endif + const RC ret = grb::clear( A ); + return ret == SUCCESS ? UNSUPPORTED : ret; + } + + /** * @} */ diff --git a/include/graphblas/bsp1d/blas3.hpp b/include/graphblas/bsp1d/blas3.hpp index 386beb164..6531dbeea 100644 --- a/include/graphblas/bsp1d/blas3.hpp +++ b/include/graphblas/bsp1d/blas3.hpp @@ -205,6 +205,50 @@ namespace grb { return internal::checkGlobalErrorStateOrClear( C, ret ); } + /** \internal Simply delegates to process-local backend */ + template< + Descriptor descr = descriptors::no_operation, + class Operator, + typename InputType1, typename InputType2, + typename MaskType, typename OutputType, + typename Coords, + typename RIT, typename CIT, typename NIT + > + RC outer( + Matrix< OutputType, BSP1D, RIT, CIT, NIT > &A, + const Matrix< MaskType, BSP1D, RIT, CIT, NIT > &mask, + const Vector< InputType1, BSP1D, Coords > &u, + const Vector< InputType2, BSP1D, Coords > &v, + const Operator &mul = Operator(), + const Phase &phase = EXECUTE, + const typename std::enable_if< + grb::is_operator< Operator >::value && + !grb::is_object< InputType1 >::value && + !grb::is_object< InputType2 >::value && + !grb::is_object< MaskType >::value && + !grb::is_object< OutputType >::value, + void >::type * const = nullptr + ) { + assert( phase != TRY ); + RC ret = outer< descr, Operator >( + internal::getLocal( A ), + internal::getLocal( mask ), + internal::getLocal( u ), + internal::getLocal( v ), + mul, + phase + ); + if( phase == RESIZE ) { + if( collectives<>::allreduce( ret, operators::any_or< RC >() ) != SUCCESS ) { + return PANIC; + } else { + return SUCCESS; + } + } + assert( phase == EXECUTE ); + return internal::checkGlobalErrorStateOrClear( A, ret ); + } + } // namespace grb #endif diff --git a/include/graphblas/hyperdags/blas3.hpp b/include/graphblas/hyperdags/blas3.hpp index ee0c10f36..16a52b9b9 100644 --- a/include/graphblas/hyperdags/blas3.hpp +++ b/include/graphblas/hyperdags/blas3.hpp @@ -257,6 +257,59 @@ namespace grb { return ret; } + template< + Descriptor descr = descriptors::no_operation, + class Operator, + typename InputType1, typename InputType2, + typename MaskType, typename OutputType, + typename Coords, + typename RIT, typename CIT, typename NIT + > + RC outer( + Matrix< OutputType, hyperdags, RIT, CIT, NIT > &A, + const Matrix< MaskType, hyperdags, RIT, CIT, NIT > &mask, + const Vector< InputType1, hyperdags, Coords > &u, + const Vector< InputType2, hyperdags, Coords > &v, + const Operator &mul = Operator(), + const Phase &phase = EXECUTE, + const typename std::enable_if< + grb::is_operator< Operator >::value && + !grb::is_object< InputType1 >::value && + !grb::is_object< InputType2 >::value && + !grb::is_object< MaskType >::value && + !grb::is_object< OutputType >::value, + void >::type * const = nullptr + ) { + const RC ret = outer< descr >( + internal::getMatrix( A ), + internal::getMatrix( mask ), + internal::getVector( u ), + internal::getVector( v ), + mul, + phase + ); + if( ret != SUCCESS ) { return ret; } + if( phase != EXECUTE ) { return ret; } + if( nrows( A ) == 0 || ncols( A ) == 0 ) { return ret; } + std::array< const void *, 0 > sourcesP{}; + std::array< uintptr_t, 4 > sourcesC{ + getID( internal::getVector(u) ), + getID( internal::getVector(v) ), + getID( internal::getMatrix(mask) ), + getID( internal::getMatrix(A) ) + }; + std::array< uintptr_t, 1 > destinations{ + getID( internal::getMatrix(A) ) + }; + internal::hyperdags::generator.addOperation( + internal::hyperdags::MASKED_OUTER, + sourcesP.begin(), sourcesP.end(), + sourcesC.begin(), sourcesC.end(), + destinations.begin(), destinations.end() + ); + return ret; + } + template< Descriptor descr = descriptors::no_operation, typename OutputType, typename InputType1, typename InputType2, diff --git a/include/graphblas/hyperdags/hyperdags.hpp b/include/graphblas/hyperdags/hyperdags.hpp index 4ef0e0059..01bc272c8 100644 --- a/include/graphblas/hyperdags/hyperdags.hpp +++ b/include/graphblas/hyperdags/hyperdags.hpp @@ -380,12 +380,16 @@ namespace grb { SET_MATRIX_MATRIX_INPUT2, + SET_MATRIX_MATRIX_MASKED, + MXM_MATRIX_MATRIX_MATRIX_SEMIRING, MXM_MATRIX_MATRIX_MATRIX_MONOID, OUTER, + MASKED_OUTER, + UNZIP_VECTOR_VECTOR_VECTOR, ZIP_MATRIX_VECTOR_VECTOR_VECTOR, @@ -493,7 +497,7 @@ namespace grb { }; /** \internal How many operation vertex types exist. */ - const constexpr size_t numOperationVertexTypes = 106; + const constexpr size_t numOperationVertexTypes = 107; /** \internal An array of all operation vertex types. */ const constexpr enum OperationVertexType @@ -553,6 +557,7 @@ namespace grb { MXM_MATRIX_MATRIX_MATRIX_SEMIRING, MXM_MATRIX_MATRIX_MATRIX_MONOID, OUTER, + MASKED_OUTER, UNZIP_VECTOR_VECTOR_VECTOR, ZIP_MATRIX_VECTOR_VECTOR_VECTOR, ZIP_MATRIX_VECTOR_VECTOR, diff --git a/include/graphblas/hyperdags/io.hpp b/include/graphblas/hyperdags/io.hpp index 901d20c48..034d41ae4 100644 --- a/include/graphblas/hyperdags/io.hpp +++ b/include/graphblas/hyperdags/io.hpp @@ -403,6 +403,41 @@ namespace grb { return ret; } + template< + Descriptor descr = descriptors::no_operation, + typename OutputType, typename MaskType, typename InputType, + typename RIT1, typename CIT1, typename NIT1, + typename RIT2, typename CIT2, typename NIT2 + > + RC set( + Matrix< OutputType, hyperdags, RIT1, CIT1, NIT1 > &C, + const Matrix< MaskType, hyperdags, RIT2, CIT2, NIT2 > &M, + const Matrix< InputType, hyperdags, RIT2, CIT2, NIT2 > &A, + const Phase &phase = EXECUTE + ) { + const RC ret = set< descr >( + internal::getMatrix( C ), internal::getMatrix( M ), + internal::getMatrix( A ), phase + ); + if( ret != SUCCESS ) { return ret; } + if( phase != EXECUTE ) { return ret; } + if( nrows( A ) == 0 || ncols( A ) == 0 ) { return ret; } + std::array< const void *, 0 > sourcesP{}; + std::array< uintptr_t, 3 > sourcesC{ + getID( internal::getMatrix(A) ), + getID( internal::getMatrix(M) ), + getID( internal::getMatrix(C) ) + }; + std::array< uintptr_t, 1 > destinations{ getID( internal::getMatrix(C) ) }; + internal::hyperdags::generator.addOperation( + internal::hyperdags::SET_MATRIX_MATRIX_MASKED, + sourcesP.begin(), sourcesP.end(), + sourcesC.begin(), sourcesC.end(), + destinations.begin(), destinations.end() + ); + return ret; + } + template< typename DataType, typename Coords > RC clear( Vector< DataType, hyperdags, Coords > &x ) { const RC ret = clear( internal::getVector( x ) ); diff --git a/include/graphblas/nonblocking/blas3.hpp b/include/graphblas/nonblocking/blas3.hpp index c52c9aaff..9cb925efd 100644 --- a/include/graphblas/nonblocking/blas3.hpp +++ b/include/graphblas/nonblocking/blas3.hpp @@ -416,6 +416,52 @@ namespace grb { ); } + template< + Descriptor descr = descriptors::no_operation, + class Operator, + typename InputType1, typename InputType2, + typename OutputType, typename MaskType, + typename Coords, + typename RIT, typename CIT, typename NIT + > + RC outer( + Matrix< OutputType, nonblocking, RIT, CIT, NIT > &A, + const Matrix< MaskType, nonblocking, RIT, CIT, NIT > &mask, + const Vector< InputType1, nonblocking, Coords > &u, + const Vector< InputType2, nonblocking, Coords > &v, + const Operator &mul = Operator(), + const Phase &phase = EXECUTE, + const typename std::enable_if< + grb::is_operator< Operator >::value && + !grb::is_object< InputType1 >::value && + !grb::is_object< InputType2 >::value && + !grb::is_object< MaskType >::value && + !grb::is_object< OutputType >::value, + void >::type * const = nullptr + ) { + if( internal::NONBLOCKING::warn_if_not_native && + config::PIPELINE::warn_if_not_native + ) { + std::cerr << "Warning: outer (nonblocking) currently delegates " + << "to a blocking implementation.\n" + << " Further similar such warnings will be suppressed.\n"; + internal::NONBLOCKING::warn_if_not_native = false; + } + + // nonblocking execution is not supported + // first, execute any computation that is not completed + internal::le.execution(); + + // second, delegate to the reference backend + return outer< descr, Operator >( + internal::getRefMatrix( A ), + internal::getRefMatrix( mask ), + internal::getRefVector( u ), + internal::getRefVector( v ), + mul, phase + ); + } + namespace internal { template< diff --git a/include/graphblas/nonblocking/io.hpp b/include/graphblas/nonblocking/io.hpp index 64e53172d..a429d8dd2 100644 --- a/include/graphblas/nonblocking/io.hpp +++ b/include/graphblas/nonblocking/io.hpp @@ -1217,6 +1217,41 @@ namespace grb { } } + + template< + Descriptor descr = descriptors::no_operation, + typename OutputType, typename MaskType, typename InputType, + typename RIT1, typename CIT1, typename NIT1, + typename RIT2, typename CIT2, typename NIT2 + > + RC set( + Matrix< OutputType, nonblocking, RIT1, CIT1, NIT1 > &C, + const Matrix< MaskType, nonblocking, RIT2, CIT2, NIT2 > &M, + const Matrix< InputType, nonblocking, RIT2, CIT2, NIT2 > &A, + const Phase &phase = EXECUTE + ) noexcept { + if( internal::NONBLOCKING::warn_if_not_native && + config::PIPELINE::warn_if_not_native + ) { + std::cerr << "Warning: masked set (nonblocking) currently delegates " + << "to a blocking implementation.\n" + << " Further similar such warnings will be suppressed.\n"; + internal::NONBLOCKING::warn_if_not_native = false; + } + + // nonblocking execution is not supported + // first, execute any computation that is not completed + internal::le.execution(); + + // second, delegate to the reference backend + return set< descr >( + internal::getRefMatrix( C ), + internal::getRefMatrix( M ), + internal::getRefMatrix( A ), + phase + ); + } + template< Descriptor descr = descriptors::no_operation, typename InputType, diff --git a/include/graphblas/reference/blas3.hpp b/include/graphblas/reference/blas3.hpp index 868547183..cb044fc4f 100644 --- a/include/graphblas/reference/blas3.hpp +++ b/include/graphblas/reference/blas3.hpp @@ -799,6 +799,75 @@ namespace grb { return internal::matrix_zip_generic< descr, true >( A, x, y, x, phase ); } + /** + * Selecting a submatrix of matrix \a B based on the given vector of \a rows and \a cols + * Depends on the masked outer product and masked set. Currently, only the structural version is supported. + */ + + template< + Descriptor descr = descriptors::no_operation, + class Operator, + typename InputType, typename MaskType, typename OutputType, + typename Coords, + typename RIT, typename CIT, typename NIT + > + RC selectSubmatrix( + Matrix< OutputType, reference, RIT, CIT, NIT > &B, + const Matrix< InputType, reference, RIT, CIT, NIT > &A, + const Vector< MaskType, reference, Coords > &rows, + const Vector< MaskType, reference, Coords > &cols, + const typename std::enable_if< + grb::is_operator< Operator >::value && + !grb::is_object< InputType >::value && + !grb::is_object< MaskType >::value && + !grb::is_object< OutputType >::value, + void >::type * const = nullptr + ) { + //static asserts + static_assert( + ! ( descr & descriptors::transpose_matrix ), + "grb::selectSubmatrix can not be called with descriptors::transpose_matrix" + ); + static_assert( + ( descr & descriptors::structural ), + "Only the structural version is supported for grb::selectSubmatrix" + ); +#ifdef _DEBUG + std::cout << "In grb::selectSubmatrix (reference)\n"; +#endif + + const size_t nrows = size( rows ); + const size_t ncols = size( cols ); + if( nrows != grb::nrows( A ) || nrows != grb::nrows( B ) ) { + return MISMATCH; + } + + if( ncols != grb::ncols( A ) || ncols != grb::ncols( B ) ) { + return MISMATCH; + } + + //mask contains only those values that need to be selected from A + Matrix< MaskType, reference, RIT, CIT, NIT > mask( nrows, ncols ); + + RC ret = outer( mask, A, rows, cols, grb::operators::zip< MaskType, MaskType, reference >(), Phase::RESIZE ); + if( ret != SUCCESS ) { + return ret; + } + ret = outer( mask, A, rows, cols, grb::operators::zip< MaskType, MaskType, reference >() ); + if( ret != SUCCESS ) { + return ret; + } + + ret = set( B, mask, A, Phase::RESIZE ); + if( ret != SUCCESS ) { + return ret; + } + ret = set( B, mask, A ); + return ret; + } + + + /** * Outer product of two vectors. Assuming vectors \a u and \a v are oriented * column-wise, the result matrix \a A will contain \f$ uv^T \f$. This is an @@ -916,6 +985,243 @@ namespace grb { return ret; } + + /** + * A masked outer product of two vectors. Assuming vectors \a u and \a v are oriented + * column-wise, the result matrix \a A will contain \f$ uv^T \f$, masked to non-zero values from \a mask. + */ + template< + Descriptor descr = descriptors::no_operation, + class Operator, + typename InputType1, typename InputType2, + typename MaskType, typename OutputType, + typename Coords, + typename RIT, typename CIT, typename NIT + > + RC outer( + Matrix< OutputType, reference, RIT, CIT, NIT > &A, + const Matrix< MaskType, reference, RIT, CIT, NIT > &mask, + const Vector< InputType1, reference, Coords > &u, + const Vector< InputType2, reference, Coords > &v, + const Operator &mul = Operator(), + const Phase &phase = EXECUTE, + const typename std::enable_if< + grb::is_operator< Operator >::value && + !grb::is_object< InputType1 >::value && + !grb::is_object< InputType2 >::value && + !grb::is_object< MaskType >::value && + !grb::is_object< OutputType >::value, + void >::type * const = nullptr + ) { + // static checks + NO_CAST_ASSERT( ( !(descr & descriptors::no_casting) || + std::is_same< typename Operator::D1, InputType1 >::value + ), "grb::outer", + "called with a prefactor vector that does not match the first domain " + "of the given multiplication operator" ); + NO_CAST_ASSERT( ( !(descr & descriptors::no_casting) || + std::is_same< typename Operator::D2, InputType2 >::value + ), "grb::outer", + "called with a postfactor vector that does not match the first domain " + "of the given multiplication operator" ); + NO_CAST_ASSERT( ( !(descr & descriptors::no_casting) || + std::is_same< typename Operator::D3, OutputType >::value + ), "grb::outer", + "called with an output matrix that does not match the output domain of " + "the given multiplication operator" ); + static_assert( + ! ( descr & descriptors::structural && descr & descriptors::invert_mask ), + "grb::outer can not be called with both descriptors::structural " + "and descriptors::invert_mask in the masked variant" + ); +#ifdef _DEBUG + std::cout << "In grb::outer (reference)\n"; +#endif + + const size_t nrows = size( u ); + const size_t ncols = size( v ); + + const size_t m = grb::nrows( mask ); + const size_t n = grb::ncols( mask ); + + if( m == 0 || n == 0 ) { + // If the mask has a null size, it will be ignored + return outer< descr >( A, u, v, mul, phase ); + } + + constexpr bool crs_only = descr & descriptors::force_row_major; + + assert( phase != TRY ); + if( nrows != grb::nrows( A ) || nrows != m ) { + return MISMATCH; + } + + if( ncols != grb::ncols( A ) || ncols != n ) { + return MISMATCH; + } + + if( nnz( u ) == 0 || nnz( v ) == 0 ) { + clear( A ); + return SUCCESS; + } + + + const auto &mask_raw = internal::getCRS( mask ); + + char * mask_arr = nullptr; + char * mask_buf = nullptr; + MaskType * mask_valbuf = nullptr; + internal::getMatrixBuffers( mask_arr, mask_buf, mask_valbuf, 1, mask ); + + internal::Coordinates< reference > mask_coors; + mask_coors.set( mask_arr, false, mask_buf, ncols ); + + size_t nzc = 0; + +#ifdef _H_GRB_REFERENCE_OMP_BLAS3 + #pragma omp parallel for reduction(+:nzc) +#endif + for( size_t i = 0; i < nrows; ++i ) { + if( internal::getCoordinates( u ).assigned( i ) ) { + for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { + const auto k_col = mask_raw.row_index[ k ]; + if( + internal::getCoordinates( v ).assigned( k_col ) && + utils::interpretMatrixMask< descr, MaskType >( true, mask_raw.getValues(), k ) + ) { + + nzc++; + } + } + } + } + + if( phase == RESIZE ) { + return resize( A, nzc ); + } + + assert( phase == EXECUTE ); + if( capacity( A ) < nzc ) { +#ifdef _DEBUG + std::cout << "\t insufficient capacity to complete " + "requested masked outer-product computation\n"; +#endif + const RC clear_rc = clear( A ); + if( clear_rc != SUCCESS ) { + return PANIC; + } else { + return FAILED; + } + } + + RC ret = SUCCESS; + if( phase == EXECUTE ) { + ret = grb::clear( A ); + } + assert( nnz( A ) == 0 ); + + auto &CRS_raw = internal::getCRS( A ); + auto &CCS_raw = internal::getCCS( A ); + + const InputType1 * __restrict__ const x = internal::getRaw( u ); + const InputType2 * __restrict__ const y = internal::getRaw( v ); + config::NonzeroIndexType * A_col_index = internal::template + getReferenceBuffer< typename config::NonzeroIndexType >( ncols + 1 ); + + CRS_raw.col_start[ 0 ] = 0; + + if( !crs_only ) { + +#ifdef _H_GRB_REFERENCE_OMP_BLAS3 + #pragma omp parallel for simd +#endif + for( size_t j = 0; j <= ncols; ++j ) { + CCS_raw.col_start[ j ] = 0; + } + } + + + nzc = 0; + + for( size_t i = 0; i < nrows; ++i ) { + if( internal::getCoordinates( u ).assigned( i ) ) { + for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { + const auto k_col = mask_raw.row_index[ k ]; + if( + internal::getCoordinates( v ).assigned( k_col ) && + utils::interpretMatrixMask< descr, MaskType >( true, mask_raw.getValues(), k ) + ) { + nzc++; + if( !crs_only ) { + CCS_raw.col_start[ k_col + 1 ]++; + } + } + } + } + CRS_raw.col_start[ i + 1 ] = nzc; + } + + if( !crs_only ) { + for( size_t j = 1; j < ncols; ++j ) { + CCS_raw.col_start[ j + 1 ] += CCS_raw.col_start[ j ]; + } + + +#ifdef _H_GRB_REFERENCE_OMP_BLAS3 + #pragma omp parallel for simd +#endif + for( size_t j = 0; j < ncols; ++j ) { + A_col_index[ j ] = 0; + } + } + + // use previously computed CCS offset array to update CCS during the + // computational phase + nzc = 0; + for( size_t i = 0; i < nrows; ++i ) { + if( internal::getCoordinates( u ).assigned( i ) ) { + for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { + const auto k_col = mask_raw.row_index[ k ]; + if( + internal::getCoordinates( v ).assigned( k_col ) && + utils::interpretMatrixMask< descr, MaskType >( true, mask_raw.getValues(), k ) + ) { + OutputType val; + grb::apply( val, + x[ i ], + y[ k_col ], + mul ); + CRS_raw.row_index[ nzc ] = k_col; + CRS_raw.setValue( nzc, val ); + // update CCS + if( !crs_only ) { + const auto CCS_index = A_col_index[ k_col ] + CCS_raw.col_start[ k_col ]; + A_col_index[ k_col ]++; + CCS_raw.row_index[ CCS_index ] = i; + CCS_raw.setValue( CCS_index, val ); + } + // update count + nzc++; + } + } + } + CRS_raw.col_start[ i + 1 ] = nzc; + } + +#ifndef NDEBUG + if( !crs_only ) { + for( size_t j = 0; j < ncols; ++j ) { + assert( CCS_raw.col_start[ j + 1 ] - CCS_raw.col_start[ j ] == + A_col_index[ j ] ); + } + } +#endif + + internal::setCurrentNonzeroes( A, nzc ); + + return ret; + } + namespace internal { /** diff --git a/include/graphblas/reference/compressed_storage.hpp b/include/graphblas/reference/compressed_storage.hpp index fa1aac97c..5d7939b87 100644 --- a/include/graphblas/reference/compressed_storage.hpp +++ b/include/graphblas/reference/compressed_storage.hpp @@ -425,6 +425,15 @@ namespace grb { return values; } + /** + * @returns The value array. + * + * \warning Does not check for NULL pointers. + */ + D * getValues() const noexcept { + return values; + } + /** * @returns The index array. * @@ -1105,6 +1114,13 @@ namespace grb { return nullptr; } + /** + * @returns A null pointer (since this is a pattern matrix). + */ + char * getValues() const noexcept { + return nullptr; + } + /** * @returns The index array. * diff --git a/include/graphblas/reference/io.hpp b/include/graphblas/reference/io.hpp index 891b13488..127394e95 100644 --- a/include/graphblas/reference/io.hpp +++ b/include/graphblas/reference/io.hpp @@ -1156,6 +1156,208 @@ namespace grb { } } + template< + Descriptor descr = descriptors::no_operation, + typename OutputType, typename MaskType, typename InputType, + typename RIT1, typename CIT1, typename NIT1, + typename RIT2, typename CIT2, typename NIT2 + > + RC set( + Matrix< OutputType, reference, RIT1, CIT1, NIT1 > &C, + const Matrix< MaskType, reference, RIT2, CIT2, NIT2 > &M, + const Matrix< InputType, reference, RIT2, CIT2, NIT2 > &A, + const Phase &phase = EXECUTE + ) noexcept { + static_assert( !std::is_void< OutputType >::value, + "grb::set (masked set to matrix): cannot have a pattern " + "matrix as output" ); + static_assert( std::is_convertible< InputType, OutputType >::value, + "grb::set (masked set to matrix): input type cannot be " + "converted to output type" + ); + static_assert( + ! ( descr & descriptors::structural && descr & descriptors::invert_mask ), + "grb::set can not be called with both descriptors::structural " + "and descriptors::invert_mask in the masked variant" + ); +#ifdef _DEBUG + std::cout << "Called grb::set (matrix-to-matrix-masked, reference)\n"; +#endif + // static checks + NO_CAST_ASSERT( ( !(descr & descriptors::no_casting) || + std::is_same< InputType, OutputType >::value + ), "grb::set", + "called with non-matching value types" + ); + + // dynamic checks + assert( phase != TRY ); + + const size_t nrows = grb::nrows( C ); + const size_t ncols = grb::ncols( C ); + + const size_t m = grb::nrows( M ); + const size_t n = grb::ncols( M ); + + /*grb::Monoid< + grb::operators::left_assign< OutputType >, + grb::identities::zero + > dummyMonoid;*/ + + if( m == 0 || n == 0 ) { + // If the mask has a null size, it will be ignored + return set< descr >( C, A, phase ); + } + + if( nrows != grb::nrows( A ) || nrows != m ) { + return MISMATCH; + } + + if( ncols != grb::ncols( A ) || ncols != n ) { + return MISMATCH; + } + + const auto &mask_raw = internal::getCRS( M ); + + const auto &A_raw = internal::getCRS( A ); + + size_t nzc = 0; + + char * mask_arr = nullptr; + char * mask_buf = nullptr; + MaskType * mask_valbuf = nullptr; + internal::getMatrixBuffers( mask_arr, mask_buf, mask_valbuf, 1, M ); + + internal::Coordinates< reference > mask_coors; + mask_coors.set( mask_arr, false, mask_buf, ncols ); + + for( size_t i = 0; i < nrows; ++i ) { + mask_coors.clear(); + for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { + const auto k_col = mask_raw.row_index[ k ]; + if( utils::interpretMatrixMask< descr, MaskType >( true, mask_raw.getValues(), k ) ) { + mask_coors.assign( k_col ); + } + } +#ifdef _H_GRB_REFERENCE_OMP_BLAS3 + #pragma omp parallel for reduction(+:nzc) +#endif + for( auto k = A_raw.col_start[ i ]; k < A_raw.col_start[ i + 1 ]; ++k ) { + const auto k_col = A_raw.row_index[ k ]; + if( mask_coors.assigned( k_col ) ) { + nzc++; + } + } + } + + if( phase == RESIZE ) { + return resize( C, nzc ); + } + + assert( phase == EXECUTE ); + + if( capacity( C ) < nzc ) { +#ifdef _DEBUG + std::cout << "\t insufficient capacity to complete " + "requested masked set matrix to matrix computation\n"; +#endif + const RC clear_rc = clear( C ); + if( clear_rc != SUCCESS ) { + return PANIC; + } else { + return FAILED; + } + } + + auto &CRS_raw = internal::getCRS( C ); + auto &CCS_raw = internal::getCCS( C ); + + config::NonzeroIndexType * C_col_index = internal::template + getReferenceBuffer< typename config::NonzeroIndexType >( ncols + 1 ); + + CRS_raw.col_start[ 0 ] = 0; + +#ifdef _H_GRB_REFERENCE_OMP_BLAS3 + #pragma omp parallel for simd +#endif + for( size_t j = 0; j <= ncols; ++j ) { + CCS_raw.col_start[ j ] = 0; + } + + + nzc = 0; + for( size_t i = 0; i < nrows; ++i ) { + mask_coors.clear(); + for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { + const auto k_col = mask_raw.row_index[ k ]; + if( utils::interpretMask< descr, MaskType >( true, mask_raw.getValues(), k ) ) { + mask_coors.assign( k_col ); + } + } + for( auto k = A_raw.col_start[ i ]; k < A_raw.col_start[ i + 1 ]; ++k ) { + const auto k_col = A_raw.row_index[ k ]; + if( mask_coors.assigned( k_col ) ) { + nzc++; + CCS_raw.col_start[ k_col + 1 ]++; + } + } + CRS_raw.col_start[ i + 1 ] = nzc; + } + + for( size_t j = 1; j < ncols; ++j ) { + CCS_raw.col_start[ j + 1 ] += CCS_raw.col_start[ j ]; + } + +#ifdef _H_GRB_REFERENCE_OMP_BLAS3 + #pragma omp parallel for simd +#endif + for( size_t j = 0; j < ncols; ++j ) { + C_col_index[ j ] = 0; + } + + + // use previously computed CCS offset array to update CCS during the + // computational phase + nzc = 0; + for( size_t i = 0; i < nrows; ++i ) { + mask_coors.clear(); + for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { + const auto k_col = mask_raw.row_index[ k ]; + if( utils::interpretMatrixMask< descr, MaskType >( true, mask_raw.getValues(), k ) ) { + mask_coors.assign( k_col ); + } + } + for( auto k = A_raw.col_start[ i ]; k < A_raw.col_start[ i + 1 ]; ++k ) { + const auto k_col = A_raw.row_index[ k ]; + if( mask_coors.assigned( k_col ) ) { + OutputType val = A_raw.getValue( k, *( (OutputType*) nullptr ) ); + CRS_raw.row_index[ nzc ] = k_col; + CRS_raw.setValue( nzc, val ); + const size_t CCS_index = C_col_index[ k_col ] + CCS_raw.col_start[ k_col ]; + C_col_index[ k_col ]++; + CCS_raw.row_index[ CCS_index ] = i; + CCS_raw.setValue( CCS_index, val ); + nzc++; + } + } + } +#ifndef NDEBUG + for( size_t j = 0; j < ncols; ++j ) { + assert( CCS_raw.col_start[ j + 1 ] - CCS_raw.col_start[ j ] == + C_col_index[ j ] ); + } +#endif + + internal::setCurrentNonzeroes( C, nzc ); + + return SUCCESS; + } + + + + + + /** * Ingests raw data into a GraphBLAS vector. * diff --git a/include/graphblas/utils.hpp b/include/graphblas/utils.hpp index c5239afdc..e62a4fb8b 100644 --- a/include/graphblas/utils.hpp +++ b/include/graphblas/utils.hpp @@ -323,7 +323,10 @@ namespace grb { ret = assigned; } else { // if based on value, if there is a value, cast it to bool - if( assigned ) { + if( !assigned ) { + return false; + } + else { ret = static_cast< bool >( val[ offset ] ); } // otherwise there is no value and false is assumed @@ -350,7 +353,10 @@ namespace grb { ret = assigned; } else { // if based on value, if there is a value, cast it to bool - if( assigned ) { + if( !assigned ) { + return false; + } + else { ret = static_cast< bool >( real( val [ offset ] ) ) || static_cast< bool >( imag( val [ offset ] ) ); } @@ -374,7 +380,40 @@ namespace grb { // set default mask to false bool ret = assigned; // check whether we should return the inverted value - if( descriptor & descriptors::invert_mask ) { + if( ( descriptor & descriptors::structural_complement ) == descriptors::structural_complement ) { + return !ret; + } else { + return ret; + } + } + + /** Specialisation for void-valued matrice's masks */ + template< Descriptor descriptor, typename MatrixDataType, typename ValuesType > + static bool interpretMatrixMask( + const bool &assigned, + const ValuesType * const values, + const size_t k, + typename std::enable_if< + !std::is_void< MatrixDataType >::value + >::type * = nullptr + ) { + return interpretMask< descriptor, ValuesType >( assigned, values, k ); + } + + /** Specialisation for void-valued matrice's masks */ + template< Descriptor descriptor, typename MatrixDataType, typename ValuesType > + static bool interpretMatrixMask( + const bool &assigned, + const ValuesType * const, + const size_t, + typename std::enable_if< + std::is_void< MatrixDataType >::value + >::type * = nullptr + ) { + // set default mask to false + bool ret = assigned; + // check whether we should return the inverted value + if( ( descriptor & descriptors::structural_complement ) == descriptors::structural_complement ) { return !ret; } else { return ret; diff --git a/src/graphblas/hyperdags/hyperdags.cpp b/src/graphblas/hyperdags/hyperdags.cpp index 6000f3af7..94bd0a98e 100644 --- a/src/graphblas/hyperdags/hyperdags.cpp +++ b/src/graphblas/hyperdags/hyperdags.cpp @@ -215,11 +215,17 @@ std::string grb::internal::hyperdags::toString( case SET_MATRIX_MATRIX_INPUT2: return "set( matrix, matrix, scalar )"; + case SET_MATRIX_MATRIX_MASKED: + return "set( matrix, matrix, matrix )"; + case MXM_MATRIX_MATRIX_MATRIX_MONOID: return "mxm( matrix, matrix, matrix, monoid, scalar, scalar )"; case OUTER: - return "outer( matrix, vector, vector, scalar, scalar )"; + return "outer( matrix, vector, vector, operation )"; + + case MASKED_OUTER: + return "outer( matrix, matrix, vector, vector, operation )"; case MXV_VECTOR_VECTOR_MATRIX_VECTOR_VECTOR_R: return "mxv( vector, vector, matrix, vector, vector, ring )"; diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt index ab0dac498..0d61e2521 100644 --- a/tests/unit/CMakeLists.txt +++ b/tests/unit/CMakeLists.txt @@ -259,6 +259,10 @@ add_grb_executables( outer outer.cpp BACKENDS reference reference_omp hyperdags nonblocking ) +add_grb_executables( maskedOuter maskedOuter.cpp + BACKENDS reference reference_omp bsp1d hybrid hyperdags nonblocking +) + # must generate the golden output for other tests force_add_grb_executable( mxv mxv.cpp BACKEND reference diff --git a/tests/unit/maskedOuter.cpp b/tests/unit/maskedOuter.cpp new file mode 100644 index 000000000..4256a9a3d --- /dev/null +++ b/tests/unit/maskedOuter.cpp @@ -0,0 +1,513 @@ + +/* + * Copyright 2021 Huawei Technologies Co., Ltd. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include + +#include "graphblas.hpp" + + +using namespace grb; + +// sample data + +static const double m1[ 4 ] = { 0.5, 3.4, 5, 0 }; + +static const size_t I1[ 4 ] = { 0, 1, 2, 0 }; +static const size_t J1[ 4 ] = { 0, 1, 2, 2 }; + +static const double m2[ 8 ] = { 1, 1, 0, 0, 0, 0, 0, 0 }; + +static const size_t I2[ 8 ] = { 0, 2, 0, 0, 1, 1, 2, 2 }; +static const size_t J2[ 8 ] = { 0, 2, 1, 2, 0, 2, 0, 1 }; + +static const double mask_test1_expect[ 3 ] = { 4, 10, 18 }; +static const double mask_test2_expect[ 3 ] = { 11, 20, 27 }; + +void grb_program( const void *, const size_t in_size, int &error ) { + error = 0; + + if( in_size != 0 ) { + (void)fprintf( stderr, "Unit tests called with unexpected input\n" ); + error = 1; + return; + } + + // allocate + grb::Vector< double > u = { 1, 2, 3 }; + grb::Vector< double > v = { 4, 5, 6 }; + grb::Matrix< double > Result1( 3, 3 ); + grb::Matrix< double > Result2( 3, 3 ); + grb::Matrix< double > Mask1( 3, 3 ); + grb::Matrix< double > Mask2( 3, 3 ); + grb::Vector< double > test1 = { 1, 1, 1 }; + grb::Vector< double > out1( 3 ); + grb::Vector< double > test2 = { 1, 1, 1 }; + grb::Vector< double > out2( 3 ); + + // semiring + grb::Semiring< + grb::operators::add< double >, grb::operators::mul< double >, + grb::identities::zero, grb::identities::one + > ring; + + grb::RC rc = grb::buildMatrixUnique( Mask1, &( I1[ 0 ] ), &( J1[ 0 ] ), m1, 4, SEQUENTIAL ); + + if( rc != SUCCESS ) { + std::cerr << "\t first mask buildMatrix FAILED\n"; + error = 5; + } + + if( !error ) { + rc = grb::buildMatrixUnique( Mask2, &( I2[ 0 ] ), &( J2[ 0 ] ), m2, 8, SEQUENTIAL ); + if( rc != SUCCESS ) { + std::cerr << "\t second mask buildMatrix FAILED\n"; + error = 10; + } + } + + + + if( !error ) { + rc = grb::outer( Result1, Mask1, u, v, ring.getMultiplicativeOperator(), RESIZE ); + rc = rc ? rc : grb::outer( Result1, Mask1, u, v, ring.getMultiplicativeOperator() ); + if( rc != grb::SUCCESS ) { + std::cerr << "Unexpected return code from grb::outer: " + << toString( rc ) << ".\n"; + error = 15; + } + } + + + if( !error && grb::nnz( Result1 ) != 3 ) { + std::cerr << "\t Unexpected number of nonzeroes in matrix: " + << grb::nnz( Result1 ) << ", expected 3.\n"; + error = 20; + } + + if( !error ) { + rc = grb::outer< descriptors::force_row_major | descriptors::invert_mask >( Result2, Mask2, u, v, ring.getMultiplicativeOperator(), RESIZE ); + rc = rc ? rc : grb::outer< descriptors::force_row_major | descriptors::invert_mask >( Result2, Mask2, u, v, ring.getMultiplicativeOperator() ); + if( rc != grb::SUCCESS ) { + std::cerr << "Unexpected return code from grb::outer: " + << toString( rc ) << ".\n"; + error = 25; + } + } + + + if( !error && grb::nnz( Result2 ) != 6 ) { + std::cerr << "\t Unexpected number of nonzeroes in matrix: " + << grb::nnz( Result2 ) << ", expected 6.\n"; + error = 30; + } + + if( !error ) { + rc = grb::vxm( out1, test1, Result1, ring ); + if( rc != grb::SUCCESS ) { + std::cerr << "Unexpected return code from premultiplying Result1 by a vector (vxm): " + << toString( rc ) << ".\n"; + error = 35; + } + } + + if( !error ) { + if( grb::nnz( out1 ) != 3 ) { + std::cerr << "\t Unexpected number of nonzeroes (premultiply): " + << grb::nnz( out1 ) << ", expected 3\n"; + error = 40; + } + for( const auto &pair : out1 ) { + size_t i = pair.first; + if( pair.second != mask_test1_expect[ i ] ) { + std::cerr << "Premultiplying Result1 by a vector of all ones, " + << "unexpected value " << pair.second << " " + << "at coordinate " << i << ", expected " + << mask_test1_expect[ i ] << ".\n"; + error = 45; + break; + } + } + } + + if( !error ) { + rc = grb::vxm< grb::descriptors::transpose_matrix >( out2, test2, Result2, ring ); + if( rc != grb::SUCCESS ) { + std::cerr << "Unexpected return code from postmultiplying Result2 by a vector (vxm): " + << toString( rc ) << ".\n"; + error = 60; + } + } + + if( !error ) { + if( grb::nnz( out2 ) != 3 ) { + std::cerr << "\t Unexpected number of nonzeroes (postmultiply): " + << grb::nnz( out1 ) << ", expected 3\n"; + error = 65; + } + for( const auto &pair : out2 ) { + size_t i = pair.first; + if( pair.second != mask_test2_expect[ i ] ) { + std::cerr << "Postmultiplying Result2 by a vector of all ones, " + << "unexpected value " << pair.second << " " + << "at coordinate " << i << ", " + << "expected " << mask_test2_expect[ i ] << ".\n"; + error = 70; + break; + } + } + } +} + +void grb_program_custom_size( const size_t &n, int &error ) { + grb::Semiring< + grb::operators::add< double >, grb::operators::mul< double >, + grb::identities::zero, grb::identities::one + > ring; + + // initialize test + grb::Matrix< double > Result( n, n ); + grb::Matrix< double > Mask( n, n ); + size_t I[ 2 * n - 1 ], J[ 2 * n - 1 ]; + double mask_in [ 2 * n - 1 ]; + double u_in[ n ]; + double v_in[ n ]; + double test_in[ n ]; + double expect[ n ]; + grb::Vector< double > u( n ); + grb::Vector< double > v( n ); + grb::Vector< double > test( n ); + grb::Vector< double > out( n ); + for( size_t k = 0; k < n; ++k ) { + I[ k ] = J[ k ] = k; + mask_in[ k ] = 1; + u_in[ k ] = k + 1; + v_in[ k ] = k + 1; + test_in[ k ] = 1; + if( k < n - 1 ) { + I[ n + k ] = k; + J[ n + k ] = k + 1; + mask_in [ n + k ] = 1; + } + if( k == 0 ) { + expect [ k ] = 1; + } + else { + expect [ k ] = ( k + 1 ) * ( 2 * k + 1 ); + } + } + + const double * vec_iter = &(u_in[ 0 ]); + grb::RC rc = grb::buildVector( u, vec_iter, vec_iter + n, SEQUENTIAL ); + if( rc != SUCCESS ) { + std::cerr << "\t buildVector of u vector FAILED\n"; + error = 5; + } + + if( !error ) { + vec_iter = &(v_in[ 0 ]); + rc = grb::buildVector( v, vec_iter, vec_iter + n, SEQUENTIAL ); + if( rc != SUCCESS ) { + std::cerr << "\t buildVector of v vector FAILED\n"; + error = 10; + } + } + + if( !error ) { + vec_iter = &(test_in[ 0 ]); + rc = grb::buildVector( test, vec_iter, vec_iter + n, SEQUENTIAL ); + if( rc != SUCCESS ) { + std::cerr << "\t buildVector of test input vector FAILED\n"; + error = 15; + } + } + + if( !error ) { + rc = grb::resize( Mask, 2 * n - 1 ); + if( rc != SUCCESS ) { + std::cerr << "\t mask matrix resize FAILED\n"; + error = 20; + } + } + + if( !error ) { + rc = grb::buildMatrixUnique( Mask, I, J, mask_in, 2 * n - 1, SEQUENTIAL ); + if( rc != SUCCESS ) { + std::cerr << "\t buildMatrixUnique of mask matrix FAILED\n"; + error = 25; + } + } + + if( !error ) { + rc = grb::outer< descriptors::structural >( Result, Mask, u, v, ring.getMultiplicativeOperator(), RESIZE ); + rc = rc ? rc : grb::outer< descriptors::structural >( Result, Mask, u, v, ring.getMultiplicativeOperator() ); + if( rc != grb::SUCCESS ) { + std::cerr << "Unexpected return code from grb::outer: " + << toString( rc ) << ".\n"; + error = 30; + } + } + + + if( !error && grb::nnz( Result ) != 2 * n -1 ) { + std::cerr << "\t Unexpected number of nonzeroes in matrix: " + << grb::nnz( Result ) << ", expected " + << 2 * n - 1 <<".\n"; + error = 35; + } + + if( !error ) { + rc = grb::vxm( out, test, Result, ring ); + if( rc != grb::SUCCESS ) { + std::cerr << "Unexpected return code from premultiplying Result by a vector (vxm): " + << toString( rc ) << ".\n"; + error = 40; + } + } + + if( !error ) { + if( grb::nnz( out ) != n ) { + std::cerr << "\t Unexpected number of nonzeroes (premultiply): " + << grb::nnz( out ) << ", expected " + << n << ".\n"; + error = 45; + } + for( const auto &pair : out ) { + size_t i = pair.first; + if( pair.second != expect[ i ] ) { + std::cerr << "Premultiplying Result by a vector of all ones, " + << "unexpected value " << pair.second << " " + << "at coordinate " << i << ", expected " + << expect[ i ] << ".\n"; + error = 50; + break; + } + } + } + +} + +static const double mask_void_test1_expect[ 3 ] = { 4, 10, 24 }; +static const double mask_void_test2_expect[ 3 ] = { 15, 20, 45 }; + +void grb_program_void( const void *, const size_t in_size, int &error ) { + error = 0; + + if( in_size != 0 ) { + (void)fprintf( stderr, "Unit tests called with unexpected input\n" ); + error = 1; + return; + } + + // allocate + grb::Vector< double > u = { 1, 2, 3 }; + grb::Vector< double > v = { 4, 5, 6 }; + grb::Matrix< double > Result1( 3, 3 ); + grb::Matrix< double > Result2( 3, 3 ); + grb::Matrix< void > Mask1( 3, 3 ); + grb::Matrix< void > Mask2( 3, 3 ); + grb::Vector< double > test1 = { 1, 1, 1 }; + grb::Vector< double > out1( 3 ); + grb::Vector< double > test2 = { 1, 1, 1 }; + grb::Vector< double > out2( 3 ); + + // semiring + grb::Semiring< + grb::operators::add< double >, grb::operators::mul< double >, + grb::identities::zero, grb::identities::one + > ring; + + grb::RC rc = grb::buildMatrixUnique( Mask1, &( I1[ 0 ] ), &( J1[ 0 ] ), 4, SEQUENTIAL ); + + if( rc != SUCCESS ) { + std::cerr << "\t first mask buildMatrix FAILED\n"; + error = 5; + } + + if( !error ) { + rc = grb::buildMatrixUnique( Mask2, &( I2[ 0 ] ), &( J2[ 0 ] ), 8, SEQUENTIAL ); + if( rc != SUCCESS ) { + std::cerr << "\t second mask buildMatrix FAILED\n"; + error = 10; + } + } + + + + if( !error ) { + rc = grb::outer( Result1, Mask1, u, v, ring.getMultiplicativeOperator(), RESIZE ); + rc = rc ? rc : grb::outer( Result1, Mask1, u, v, ring.getMultiplicativeOperator() ); + if( rc != grb::SUCCESS ) { + std::cerr << "Unexpected return code from grb::outer: " + << toString( rc ) << ".\n"; + error = 15; + } + } + + + if( !error && grb::nnz( Result1 ) != 4 ) { + std::cerr << "\t Unexpected number of nonzeroes in matrix: " + << grb::nnz( Result1 ) << ", expected 4.\n"; + error = 20; + } + + if( !error ) { + rc = grb::outer< descriptors::force_row_major | descriptors::invert_mask >( Result2, Mask2, u, v, ring.getMultiplicativeOperator(), RESIZE ); + rc = rc ? rc : grb::outer< descriptors::force_row_major | descriptors::invert_mask >( Result2, Mask2, u, v, ring.getMultiplicativeOperator() ); + if( rc != grb::SUCCESS ) { + std::cerr << "Unexpected return code from grb::outer: " + << toString( rc ) << ".\n"; + error = 25; + } + } + + + if( !error && grb::nnz( Result2 ) != 8 ) { + std::cerr << "\t Unexpected number of nonzeroes in matrix: " + << grb::nnz( Result2 ) << ", expected 8.\n"; + error = 30; + } + + if( !error ) { + rc = grb::vxm( out1, test1, Result1, ring ); + if( rc != grb::SUCCESS ) { + std::cerr << "Unexpected return code from premultiplying Result1 by a vector (vxm): " + << toString( rc ) << ".\n"; + error = 35; + } + } + + if( !error ) { + if( grb::nnz( out1 ) != 3 ) { + std::cerr << "\t Unexpected number of nonzeroes (premultiply): " + << grb::nnz( out1 ) << ", expected 3\n"; + error = 40; + } + for( const auto &pair : out1 ) { + size_t i = pair.first; + if( pair.second != mask_void_test1_expect[ i ] ) { + std::cerr << "Premultiplying Result1 by a vector of all ones, " + << "unexpected value " << pair.second << " " + << "at coordinate " << i << ", expected " + << mask_void_test1_expect[ i ] << ".\n"; + error = 45; + break; + } + } + } + + if( !error ) { + rc = grb::vxm< grb::descriptors::transpose_matrix >( out2, test2, Result2, ring ); + if( rc != grb::SUCCESS ) { + std::cerr << "Unexpected return code from postmultiplying Result2 by a vector (vxm): " + << toString( rc ) << ".\n"; + error = 60; + } + } + + if( !error ) { + if( grb::nnz( out2 ) != 3 ) { + std::cerr << "\t Unexpected number of nonzeroes (postmultiply): " + << grb::nnz( out1 ) << ", expected 3\n"; + error = 65; + } + for( const auto &pair : out2 ) { + size_t i = pair.first; + if( pair.second != mask_void_test2_expect[ i ] ) { + std::cerr << "Postmultiplying Result2 by a vector of all ones, " + << "unexpected value " << pair.second << " " + << "at coordinate " << i << ", " + << "expected " << mask_void_test2_expect[ i ] << ".\n"; + error = 70; + break; + } + } + } +} + +int main( int argc, char ** argv ) { + (void)argc; + std::cout << "Functional test executable: " << argv[ 0 ] << "\n"; + + bool printUsage = false; + size_t n = 100; + + // error checking + if( argc > 2 ) { + printUsage = true; + } + if( argc == 2 ) { + size_t read; + std::istringstream ss( argv[ 1 ] ); + if( ! ( ss >> read ) ) { + std::cerr << "Error parsing first argument\n"; + printUsage = true; + } else if( ! ss.eof() ) { + std::cerr << "Error parsing first argument\n"; + printUsage = true; + } else { + // all OK + n = read; + } + } + if( printUsage ) { + std::cerr << "Usage: " << argv[ 0 ] << " [n]\n"; + std::cerr << " -n (optional, default is 100): an integer, the " + "test size.\n"; + return 1; + } + + + + int error = 0; + grb::Launcher< AUTOMATIC > launcher; + + if( launcher.exec( &grb_program, nullptr, 0, error ) != SUCCESS ) { + std::cerr << "Test 1 failed to launch\n"; + error = 255; + } + if( error != 0 ) { + std::cerr << std::flush; + std::cout << "Test 1 FAILED\n" << std::endl; + return error; + } + + if( launcher.exec( &grb_program_custom_size, n, error ) != SUCCESS ) { + std::cerr << "Launching test 2 FAILED\n"; + error = 255; + } + if( error != 0 ) { + std::cerr << std::flush; + std::cout << "Test 2 FAILED\n" << std::endl; + return error; + } + + if( launcher.exec( &grb_program_void, nullptr, 0, error ) != SUCCESS ) { + std::cerr << "Test 3 failed to launch\n"; + error = 255; + } + if( error != 0 ) { + std::cerr << std::flush; + std::cout << "Test 3 FAILED\n" << std::endl; + return error; + } + + // done + std::cout << "Test OK\n" << std::endl; + return 0; +} + diff --git a/tests/unit/matrixSet.cpp b/tests/unit/matrixSet.cpp index dfa9c0c8c..15eac3bd8 100644 --- a/tests/unit/matrixSet.cpp +++ b/tests/unit/matrixSet.cpp @@ -43,6 +43,12 @@ void grb_program( const size_t &n, grb::RC &rc ) { grb::Matrix< void > C( n, n ); grb::Matrix< void > D( n, n ); grb::Matrix< unsigned int > E( n, n ); + grb::Matrix< int > mask( n, n ); + grb::Matrix< int > output( n, n ); + grb::Matrix< int > input( n, n ); + + + rc = grb::resize( A, 15 ); if( rc == SUCCESS ) { rc = grb::buildMatrixUnique( A, I, J, data1, 15, SEQUENTIAL ); @@ -61,6 +67,35 @@ void grb_program( const size_t &n, grb::RC &rc ) { } } } + + size_t I_mask[ 2 * n - 1 ], J_mask[ 2 * n - 1 ]; + int mask_vals [ 2 * n - 1 ]; + int input_vals [ 2 * n - 1 ]; + + for( size_t k = 0; k < n; ++k ) { + I_mask[ k ] = J_mask[ k ] = k; + mask_vals[ k ] = 1; + input_vals[ k ] = k; + if( k < n - 1 ) { + I_mask[ n + k ] = k; + J_mask[ n + k ] = k + 1; + mask_vals [ n + k ] = 0; + input_vals[ n + k ] = k; + } + } + + rc = grb::buildMatrixUnique( mask, I_mask, J_mask, mask_vals, 2 * n - 1, SEQUENTIAL ); + if( rc != SUCCESS ) { + std::cerr << "\t buildMatrixUnique of mask matrix FAILED\n"; + return; + } + + rc = grb::buildMatrixUnique( input, I_mask, J_mask, input_vals, 2 * n - 1, SEQUENTIAL ); + if( rc != SUCCESS ) { + std::cerr << "\t buildMatrixUnique of input matrix FAILED\n"; + return; + } + if( rc == SUCCESS ) { rc = grb::resize( B, 15 ); } @@ -73,6 +108,9 @@ void grb_program( const size_t &n, grb::RC &rc ) { if( rc == SUCCESS ) { rc = grb::resize( E, 15 ); } + if( rc == SUCCESS ) { + rc = grb::resize(output, 2 * n - 1 ); + } if( rc != SUCCESS || grb::nnz( A ) != 15 ) { std::cerr << "\tinitialisation FAILED\n"; return; @@ -201,6 +239,90 @@ void grb_program( const size_t &n, grb::RC &rc ) { if( rc != SUCCESS ) { return; } + + //check masked matrix set + + rc = grb::set< descriptors::structural >( output, mask, input ); + if( rc != SUCCESS ) { + std::cerr << "\t grb::set structural (matrix to matrix masked) FAILED\n"; + return; + } + + if( grb::nnz( output ) != 2 * n - 1 ) { + std::cerr << "\t unexpected number of output elements ( " << grb::nnz( output ) << " ), expected " << 2 * n - 1 <<".\n"; + rc = FAILED; + } + + for( const auto & triplet : output ) { + if( triplet.first.first != triplet.first.second && triplet.first.first != triplet.first.second - 1 ) { + std::cerr << "\tunexpected entry at ( " << triplet.first.first << ", " << triplet.first.second << " ), value " << triplet.second << ".\n"; + rc = FAILED; + } if( (int) triplet.first.first != triplet.second ) { + std::cerr << "\tunexpected entry at ( " << triplet.first.first << ", " << triplet.first.second << " ) with value " << triplet.second; + std::cerr << ", expected value "<< triplet.first.first <<".\n"; + rc = FAILED; + } + } + + if( rc != SUCCESS ) { + return; + } + + + + rc = grb::set( output, mask, input ); + if( rc != SUCCESS ) { + std::cerr << "\t grb::set (matrix to matrix masked) FAILED\n"; + return; + } + + if( grb::nnz( output ) != n ) { + std::cerr << "\t unexpected number of output elements ( " << grb::nnz( output ) << " ), expected " << n <<".\n"; + rc = FAILED; + } + + for( const auto & triplet : output ) { + if( triplet.first.first != triplet.first.second ) { + std::cerr << "\tunexpected entry at ( " << triplet.first.first << ", " << triplet.first.second << " ), value " << triplet.second << ".\n"; + rc = FAILED; + } if( (int) triplet.first.first != triplet.second ) { + std::cerr << "\tunexpected entry at ( " << triplet.first.first << ", " << triplet.first.second << " ) with value " << triplet.second; + std::cerr << ", expected value "<< triplet.first.first <<".\n"; + rc = FAILED; + } + } + + if( rc != SUCCESS ) { + return; + } + + + + rc = grb::set< descriptors::invert_mask >( output, mask, input ); + if( rc != SUCCESS ) { + std::cerr << "\t grb::set invert mask (matrix to matrix masked) FAILED\n"; + return; + } + + if( grb::nnz( output ) != n - 1 ) { + std::cerr << "\t unexpected number of output elements ( " << grb::nnz( output ) << " ), expected " << n - 1 <<".\n"; + rc = FAILED; + } + + for( const auto & triplet : output ) { + if( triplet.first.first != triplet.first.second - 1 ) { + std::cerr << "\tunexpected entry at ( " << triplet.first.first << ", " << triplet.first.second << " ), value " << triplet.second << ".\n"; + rc = FAILED; + } if( (int) triplet.first.first != triplet.second ) { + std::cerr << "\tunexpected entry at ( " << triplet.first.first << ", " << triplet.first.second << " ) with value " << triplet.second; + std::cerr << ", expected value "<< triplet.first.first <<".\n"; + rc = FAILED; + } + } + + if( rc != SUCCESS ) { + return; + } } int main( int argc, char ** argv ) { diff --git a/tests/unit/unittests.sh b/tests/unit/unittests.sh index 5b7dfc16a..f3db533a0 100755 --- a/tests/unit/unittests.sh +++ b/tests/unit/unittests.sh @@ -643,6 +643,12 @@ for MODE in ${MODES}; do grep 'Test OK' ${TEST_OUT_DIR}/outer_${MODE}_${BACKEND}_${P}_${T}.log || echo "Test FAILED" echo " " + echo ">>> [x] [ ] Testing grb::maskedOuter on a small matrix" + $runner ${TEST_BIN_DIR}/maskedOuter_${MODE}_${BACKEND} &> ${TEST_OUT_DIR}/maskedOuter_${MODE}_${BACKEND}_${P}_${T}.log + head -1 ${TEST_OUT_DIR}/maskedOuter_${MODE}_${BACKEND}_${P}_${T}.log + grep 'Test OK' ${TEST_OUT_DIR}/maskedOuter_${MODE}_${BACKEND}_${P}_${T}.log || echo "Test FAILED" + echo " " + echo ">>> [x] [ ] Testing vector times matrix using the normal (+,*)" echo " semiring over integers on a diagonal matrix" echo " "