From a98f1c0a573d94e9e4c1362ae6192f6fd4e4cbd7 Mon Sep 17 00:00:00 2001 From: Denis Jelovina Date: Tue, 15 Jul 2025 18:03:59 +0200 Subject: [PATCH 01/18] Add Ising machine SB implementation and associated test files --- .../graphblas/algorithms/ising_machine_sb.hpp | 218 ++++++++++++++++++ tests/smoke/CMakeLists.txt | 5 + tests/smoke/ising_machine_sb.cpp | 121 ++++++++++ 3 files changed, 344 insertions(+) create mode 100644 include/graphblas/algorithms/ising_machine_sb.hpp create mode 100644 tests/smoke/ising_machine_sb.cpp diff --git a/include/graphblas/algorithms/ising_machine_sb.hpp b/include/graphblas/algorithms/ising_machine_sb.hpp new file mode 100644 index 000000000..de9e39599 --- /dev/null +++ b/include/graphblas/algorithms/ising_machine_sb.hpp @@ -0,0 +1,218 @@ + +/* + * Copyright 2025 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. + */ + +/** + * @file + * + * Implements the Ising machine SB + * + * @author Denis Jelovina + */ + +#ifndef _H_GRB_ALGORITHMS_ISING_MACHINE_SB +#define _H_GRB_ALGORITHMS_ISING_MACHINE_SB + +#include // master ALP/GraphBLAS header +#include // std::min/max +#include // std::abs, std::sqrt +#include +#include + +#include + + +namespace grb { + namespace algorithms { + + using IOType = double; // arithmetic scalar + using Vec = grb::Vector< IOType >; + using Mat = grb::Matrix< IOType >; + + // // Custom unary operator for sign extraction + // struct signum { + // constexpr IOType operator()( const IOType x ) const noexcept { + // return ( x > 0 ) - ( x < 0 ); + // } + // }; + + // // Custom unary operator for hard clipping to [-1,1] + // struct clip11 { + // constexpr IOType operator()( const IOType x ) const noexcept { + // return std::min( 1.0, std::max( -1.0, x ) ); + // } + // }; + + // // Unary predicate to build a structural mask |x|>1 + // struct abs_gt1 { + // constexpr bool operator()( const IOType x ) const noexcept { + // return std::abs( x ) > 1.0; + // } + // }; + + /*-------------------------------------------------------------* + * bSB — core optimisation routine * + *-------------------------------------------------------------*/ + template< + typename IsingHType, + typename IOType, + typename RSI, typename NZI, Backend backend, + class Ring = Semiring< + grb::operators::add< IOType >, grb::operators::mul< IOType >, + grb::identities::zero, grb::identities::one + >, + class Minus = operators::subtract< IOType >, + class Divide = operators::divide< IOType > + > + grb::RC bSB( + std::vector< IOType > & energies, // output length num_iters + grb::Vector< IOType, backend > & x_comp, // in/out, size N + grb::Vector< IOType, backend > & y_comp, // in/out, size N + const grb::Matrix< IsingHType, backend, RSI, RSI, NZI > & J, // NxN, symmetric + const grb::Vector< IOType, backend > & h, // size N + const IOType p_init, + const IOType p_end, + const std::size_t num_iters, + const IOType dt, + // default semiring, minus, divide + const Ring &ring = Ring(), + const Minus &minus = Minus(), + const Divide ÷ = Divide(), + const std::function< IOType( IOType ) > &sqrtX = + std_sqrt< IOType, IOType > + ) { + + const std::size_t N = grb::nrows(J); + + assert( grb::ncols(J) == N ); + assert( grb::nrows(h) == N ); + assert( grb::ncols(h) == 1 ); + assert( grb::nrows(x_comp) == N ); + assert( grb::ncols(x_comp) == 1 ); + assert( grb::nrows(y_comp) == N ); + assert( grb::ncols(y_comp) == 1 ); + // TODO: check that J is symmetric once properly implemented + //assert( grb::is_symmetric(J) ); + + if ( num_iters < 1 ) { + return grb::SUCCESS; + } + + grb::RC rc = grb::SUCCESS; + /* ---- pre-compute ---- */ + IOType sumJ2 = ring.template getZero< IOType >(); + rc = rc ? rc : grb::eWiseLambda( [&J, &sumJ2]( const size_t i, const size_t j, IOType& v ) { + (void) i; + (void) j; + // rewrite this to use graphblas language + sumJ2 += v*v; + }, J ); + if( rc != grb::SUCCESS ) { + std::cerr << "Error in eWiseLambda for sumJ2: " << rc << '\n'; + return rc; + } + // for debugging purposes, print sumJ2 + std::cout << "sumJ2: " << sumJ2 << '\n'; + + // rewrite this to use graphblas language + rc = rc ? rc : grb::foldl( sumJ2, static_cast( N - 1 ), divide ); + IOType xi = 0.5; + sumJ2 = sqrtX( sumJ2 ); + rc = rc ? rc : grb::foldl( xi, sumJ2, divide ); + // for debugging purposes, print xi + std::cout << "xi: " << xi << '\n'; + + /* ---- iteration variables ---- */ + IOType ps = p_init; + const IOType dps = ( p_end - p_init ) / static_cast( num_iters - 1 ); + // assert len of energies == N + assert( energies.size() == num_iters ); + + // Vec Jx( N ), temp( N ), mask( N ); // workspace vectors + grb::Vector< IOType, backend > Jx( N ), temp( N ), mask( N ); + + for ( std::size_t iter = 0; iter < num_iters; ++iter ) { + + // /* y_comp += ((-1+ps)*x_comp + xi*(Jx + h)) * dt */ + + // // Jx ← J * x_comp + // GRB_TRY( mxv( Jx, J, x_comp ) ); + + // // temp ← Jx + h + // GRB_TRY( eWiseApply< descriptors::no_operation >( + // temp, Jx, h, operators::add< IOType >() + // ) ); + + // // temp ← xi * temp + // GRB_TRY( apply( temp, temp, [&]( IOType v ){ return xi * v; } ) ); + + // // temp ← temp + (-1+ps) * x_comp + // const IOType scale = -1.0 + ps; + // GRB_TRY( eWiseApply( + // temp, temp, x_comp, + // [&]( IOType a, IOType b ){ return a + scale * b; } + // ) ); + + // // y_comp += dt * temp + // GRB_TRY( apply( temp, temp, [&]( IOType v ){ return v * dt; } ) ); + // GRB_TRY( eWiseApply( + // y_comp, y_comp, temp, operators::add< IOType >() + // ) ); + + // /* x_comp += dt * y_comp */ + // GRB_TRY( apply( temp, y_comp, [&]( IOType v ){ return v * dt; } ) ); + // GRB_TRY( eWiseApply( + // x_comp, x_comp, temp, operators::add< IOType >() + // ) ); + + // /* y_comp[ |x|>1 ] = 0 */ + // GRB_TRY( apply( mask, x_comp, abs_gt1() ) ); // bool structural mask + // GRB_TRY( eWiseApply< descriptors::structural >( + // y_comp, mask, [&]( IOType, bool ){ return 0.0; } + // ) ); + + // /* x_comp = clip( x_comp ) */ + // GRB_TRY( apply( x_comp, x_comp, clip11() ) ); + + // /* Energy evaluation */ + // Vec sol( N ); + // GRB_TRY( apply( sol, x_comp, signum() ) ); + + // // temp ← J * sol + // GRB_TRY( mxv( temp, J, sol ) ); + + // // e = -0.5 * sol^T * temp – h^T * sol + // IOType dot1 = 0.0, dot2 = 0.0; + // GRB_TRY( foldl( dot1, sol, temp, operators::add< IOType >(), + // operators::mul< IOType >() ) ); + // GRB_TRY( foldl( dot2, h, sol, operators::add< IOType >(), + // operators::mul< IOType >() ) ); + + energies[ iter ] = iter; // just for testing purposes + // energies[ iter ] = -0.5 * dot1 - dot2; + + ps += dps; + } + + return SUCCESS; + } + + } // algorithms namespace + +} // grb namespace + +#endif // end _H_GRB_ALGORITHMS_ISING_MACHINE_SB + diff --git a/tests/smoke/CMakeLists.txt b/tests/smoke/CMakeLists.txt index 7e7f2af4a..9ce5553c3 100644 --- a/tests/smoke/CMakeLists.txt +++ b/tests/smoke/CMakeLists.txt @@ -144,6 +144,11 @@ add_grb_executables( conjugate_gradient_complex conjugate_gradient.cpp COMPILE_DEFINITIONS _CG_COMPLEX ) +add_grb_executables( ising_machine_sb ising_machine_sb.cpp + BACKENDS reference reference_omp nonblocking + ADDITIONAL_LINK_LIBRARIES test_utils_headers +) + add_grb_executables( gmres gmres.cpp BACKENDS reference reference_omp bsp1d hybrid hyperdags nonblocking ADDITIONAL_LINK_LIBRARIES test_utils_headers diff --git a/tests/smoke/ising_machine_sb.cpp b/tests/smoke/ising_machine_sb.cpp new file mode 100644 index 000000000..48128b49d --- /dev/null +++ b/tests/smoke/ising_machine_sb.cpp @@ -0,0 +1,121 @@ + +/* + * 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 + + +#include + +#include + +#include + +#include + +#include +#include +#include + +#include + +#include + + +using namespace grb; +using namespace algorithms; + + + +using IOType = double; + +int main( int argc, char ** argv ) { + constexpr std::size_t N = 100; + constexpr std::size_t Nz = 400; + + /* --- Initialise ALP/GraphBLAS --- */ + if( grb::init() != grb::SUCCESS ) { + std::cerr << "ALP init failed\n"; + return -1; + } + + /* --- Problem setup (toy) --- */ + grb::Matrix J( N, N, Nz ); + grb::Vector h( N ); + grb::Vector x0( N ), y0( N ); // initialy + // ... populate J with random values + // generate random nonzero locations and values for J + std::vector row_indices(Nz); + std::vector col_indices(Nz); + std::vector nonzero_values(Nz); + // TODO: make sure there are not duplicate entries in J and J is symmetric + + IOType sumJ2_test = 0; + for( std::size_t i = 0; i < Nz; ++i ) { + row_indices[i] = rand() % N; + col_indices[i] = rand() % N; + nonzero_values[i] = static_cast(rand()) / RAND_MAX; // random value between 0 and 1 + sumJ2_test+=nonzero_values[i]*nonzero_values[i]; + } + std::cout << "sumJ2_test: " << sumJ2_test << '\n'; + IOType xi_test = 0.5 / std::sqrt( sumJ2_test / static_cast( N - 1 ) ); + std::cout << "xi_test: " << xi_test << '\n'; + + grb::RC rc = grb::SUCCESS; + rc = rc ? rc : buildMatrixUnique(J, row_indices.data(), col_indices.data(), nonzero_values.data(), Nz, grb::SEQUENTIAL); + if(rc != grb::SUCCESS) { + std::cerr << "matrix build failed\n"; + return grb::RC::PANIC; + } + + // Fill h, x0, y0 with random values using buildVector + std::vector h_buffer(N), x0_buffer(N), y0_buffer(N); + for( std::size_t i = 0; i < N; ++i ) { + h_buffer[i] = static_cast(rand()) / RAND_MAX; + x0_buffer[i] = static_cast(rand()) / RAND_MAX / 0.2 - 0.1; + y0_buffer[i] = static_cast(rand()) / RAND_MAX / 0.2 - 0.1; + } + rc = rc ? rc : buildVector(h, h_buffer.begin(), h_buffer.end(), grb::SEQUENTIAL); + rc = rc ? rc : buildVector(x0, x0_buffer.begin(), x0_buffer.end(), grb::SEQUENTIAL); + rc = rc ? rc : buildVector(y0, y0_buffer.begin(), y0_buffer.end(), grb::SEQUENTIAL); + if(rc != grb::SUCCESS) { + std::cerr << "Vector build failed\n"; + return grb::RC::PANIC; + } + + const IOType p0 = 0.1; + const IOType p1 = 1.0; + const std::size_t num_iters = 1000; + const IOType dt = 0.01; + + // energies is array of length num_iters, initialized to 0 + std::vector< IOType > energies( num_iters, 0 ); + + rc = rc ? rc : bSB( + energies, x0, y0, J, h, p0, p1, num_iters, dt + ); + + if( rc != grb::SUCCESS ) { + std::cerr << "bSB returned error code " << rc << '\n'; + } else { + std::cout << "Final energy = " << energies[num_iters-1] << '\n'; + } + + grb::finalize(); + return 0; +} From 298c55532871a1c6e107fdeec84b93e3ebf71e95 Mon Sep 17 00:00:00 2001 From: Denis Jelovina Date: Tue, 15 Jul 2025 19:27:13 +0200 Subject: [PATCH 02/18] Refactor Ising machine SB implementation and enhance debugging output in tests --- .../graphblas/algorithms/ising_machine_sb.hpp | 126 +++++++++++------- tests/smoke/ising_machine_sb.cpp | 6 + 2 files changed, 82 insertions(+), 50 deletions(-) diff --git a/include/graphblas/algorithms/ising_machine_sb.hpp b/include/graphblas/algorithms/ising_machine_sb.hpp index de9e39599..ff881118b 100644 --- a/include/graphblas/algorithms/ising_machine_sb.hpp +++ b/include/graphblas/algorithms/ising_machine_sb.hpp @@ -38,9 +38,7 @@ namespace grb { namespace algorithms { - using IOType = double; // arithmetic scalar - using Vec = grb::Vector< IOType >; - using Mat = grb::Matrix< IOType >; + // // Custom unary operator for sign extraction // struct signum { @@ -49,12 +47,13 @@ namespace grb { // } // }; - // // Custom unary operator for hard clipping to [-1,1] - // struct clip11 { - // constexpr IOType operator()( const IOType x ) const noexcept { - // return std::min( 1.0, std::max( -1.0, x ) ); - // } - // }; + // Custom unary operator for hard clipping to [-1,1] + template< typename IOType > + struct clip11 { + constexpr IOType operator()( const IOType x ) const noexcept { + return std::min( 1.0, std::max( -1.0, x ) ); + } + }; // // Unary predicate to build a structural mask |x|>1 // struct abs_gt1 { @@ -67,6 +66,7 @@ namespace grb { * bSB — core optimisation routine * *-------------------------------------------------------------*/ template< + Descriptor descr = descriptors::no_operation, typename IsingHType, typename IOType, typename RSI, typename NZI, Backend backend, @@ -95,6 +95,8 @@ namespace grb { std_sqrt< IOType, IOType > ) { + constexpr const Descriptor descr_dense = descr | descriptors::dense; + const std::size_t N = grb::nrows(J); assert( grb::ncols(J) == N ); @@ -124,16 +126,20 @@ namespace grb { std::cerr << "Error in eWiseLambda for sumJ2: " << rc << '\n'; return rc; } +#ifdef DEBUG // for debugging purposes, print sumJ2 std::cout << "sumJ2: " << sumJ2 << '\n'; +#endif // rewrite this to use graphblas language rc = rc ? rc : grb::foldl( sumJ2, static_cast( N - 1 ), divide ); IOType xi = 0.5; sumJ2 = sqrtX( sumJ2 ); rc = rc ? rc : grb::foldl( xi, sumJ2, divide ); +#ifdef DEBUG // for debugging purposes, print xi std::cout << "xi: " << xi << '\n'; +#endif /* ---- iteration variables ---- */ IOType ps = p_init; @@ -141,51 +147,71 @@ namespace grb { // assert len of energies == N assert( energies.size() == num_iters ); + // TODO: move these to the aggument list of bSB // Vec Jx( N ), temp( N ), mask( N ); // workspace vectors - grb::Vector< IOType, backend > Jx( N ), temp( N ), mask( N ); + grb::Vector< IOType, backend > Jx( N ), temp( N ); + grb::Vector< bool, backend > mask( N ); for ( std::size_t iter = 0; iter < num_iters; ++iter ) { - // /* y_comp += ((-1+ps)*x_comp + xi*(Jx + h)) * dt */ - - // // Jx ← J * x_comp - // GRB_TRY( mxv( Jx, J, x_comp ) ); - - // // temp ← Jx + h - // GRB_TRY( eWiseApply< descriptors::no_operation >( - // temp, Jx, h, operators::add< IOType >() - // ) ); - - // // temp ← xi * temp - // GRB_TRY( apply( temp, temp, [&]( IOType v ){ return xi * v; } ) ); - - // // temp ← temp + (-1+ps) * x_comp - // const IOType scale = -1.0 + ps; - // GRB_TRY( eWiseApply( - // temp, temp, x_comp, - // [&]( IOType a, IOType b ){ return a + scale * b; } - // ) ); - - // // y_comp += dt * temp - // GRB_TRY( apply( temp, temp, [&]( IOType v ){ return v * dt; } ) ); - // GRB_TRY( eWiseApply( - // y_comp, y_comp, temp, operators::add< IOType >() - // ) ); - - // /* x_comp += dt * y_comp */ - // GRB_TRY( apply( temp, y_comp, [&]( IOType v ){ return v * dt; } ) ); - // GRB_TRY( eWiseApply( - // x_comp, x_comp, temp, operators::add< IOType >() - // ) ); - - // /* y_comp[ |x|>1 ] = 0 */ - // GRB_TRY( apply( mask, x_comp, abs_gt1() ) ); // bool structural mask - // GRB_TRY( eWiseApply< descriptors::structural >( - // y_comp, mask, [&]( IOType, bool ){ return 0.0; } - // ) ); - - // /* x_comp = clip( x_comp ) */ - // GRB_TRY( apply( x_comp, x_comp, clip11() ) ); + /* y_comp += ((-1+ps)*x_comp + xi*(Jx + h)) * dt */ + + // Jx ← J * x_comp + rc = rc ? rc : grb::mxv< descr_dense >( Jx, J, x_comp, ring ); + assert( rc == grb::SUCCESS ); + + // temp ← Jx + h + rc = rc ? rc : grb::eWiseApply< descr_dense >( + temp, Jx, h, ring.getAdditiveMonoid() + ); + assert( rc == grb::SUCCESS ); + + // temp ← xi * temp + rc = rc ? rc : grb::foldl< descr_dense >( + temp, xi, ring.getMultiplicativeMonoid() + ); + assert( rc == grb::SUCCESS ); + + // temp ← temp + (-1+ps) * x_comp + const IOType scale = -1.0 + ps; + rc = rc ? rc : grb::eWiseMul< descr_dense >( temp, scale, x_comp, ring ); + assert( rc == grb::SUCCESS ); + + // y_comp += dt * temp + rc = rc ? rc : grb::eWiseMul< descr_dense >( y_comp, dt, temp, ring ); + assert( rc == grb::SUCCESS ); + + /* x_comp += dt * y_comp */ + rc = rc ? rc : grb::eWiseMul< descr_dense >( x_comp, dt, y_comp, ring ); + assert( rc == grb::SUCCESS ); + + /* y_comp[ |x|>1 ] = 0 */ + // mask = np.abs(x_comp) > 1 + rc = rc ? rc : grb::eWiseLambda< descr_dense >( [&mask, &x_comp]( const size_t i ) { + (void) i; + // rewrite this to use graphblas language + mask[i] = std::abs(x_comp[i]) > 1; + }, mask + ); + assert( rc == grb::SUCCESS ); + + // y_comp[ mask ] = 0 + rc = rc ? rc : grb::eWiseLambda< descr_dense >( [&mask, &y_comp]( const size_t i ) { + (void) i; + // rewrite this to use graphblas language + if(mask[i]) { + y_comp[i] = 0; + } + }, y_comp + ); + assert( rc == grb::SUCCESS ); + + /* x_comp = clip( x_comp ) */ + rc = rc ? rc : grb::eWiseLambda< descr_dense >( [&x_comp]( const size_t i ) { + (void) i; + x_comp[i] = clip11()( x_comp[i] ); }, x_comp + ); + assert( rc == grb::SUCCESS ); // /* Energy evaluation */ // Vec sol( N ); diff --git a/tests/smoke/ising_machine_sb.cpp b/tests/smoke/ising_machine_sb.cpp index 48128b49d..89d72d296 100644 --- a/tests/smoke/ising_machine_sb.cpp +++ b/tests/smoke/ising_machine_sb.cpp @@ -72,9 +72,15 @@ int main( int argc, char ** argv ) { nonzero_values[i] = static_cast(rand()) / RAND_MAX; // random value between 0 and 1 sumJ2_test+=nonzero_values[i]*nonzero_values[i]; } +#ifdef DEBUG + // for debugging purposes, print sumJ2_test std::cout << "sumJ2_test: " << sumJ2_test << '\n'; +#endif IOType xi_test = 0.5 / std::sqrt( sumJ2_test / static_cast( N - 1 ) ); +#ifdef DEBUG + // for debugging purposes, print xi_test std::cout << "xi_test: " << xi_test << '\n'; +#endif grb::RC rc = grb::SUCCESS; rc = rc ? rc : buildMatrixUnique(J, row_indices.data(), col_indices.data(), nonzero_values.data(), Nz, grb::SEQUENTIAL); From 417af03a3cde9741f276b6e18880d547056d4ae8 Mon Sep 17 00:00:00 2001 From: Denis Jelovina Date: Wed, 16 Jul 2025 14:11:19 +0200 Subject: [PATCH 03/18] Correct first version Add vector_print function and update test data for Ising machine SB --- .../graphblas/algorithms/ising_machine_sb.hpp | 190 +++++++++++++----- tests/smoke/ising_machine_sb.cpp | 69 +++---- 2 files changed, 165 insertions(+), 94 deletions(-) diff --git a/include/graphblas/algorithms/ising_machine_sb.hpp b/include/graphblas/algorithms/ising_machine_sb.hpp index ff881118b..ef4666376 100644 --- a/include/graphblas/algorithms/ising_machine_sb.hpp +++ b/include/graphblas/algorithms/ising_machine_sb.hpp @@ -38,6 +38,17 @@ namespace grb { namespace algorithms { + template< typename IOType, Backend backend > + void vector_print( grb::Vector< IOType, backend > & x_comp, const std::string & vector_name ) { + grb::PinnedVector< IOType > pinnedVector; + pinnedVector = grb::PinnedVector< IOType >( x_comp, grb::SEQUENTIAL ); + std::cout << "First 10 nonzeroes of " << vector_name << " are: ( "; + for( size_t k = 0; k < pinnedVector.nonzeroes() && k < 10; ++k ) { + const IOType & nonzeroValue = pinnedVector.getNonzeroValue( k ); + std::cout << nonzeroValue << " "; + } + std::cout << ")" << std::endl; + } // // Custom unary operator for sign extraction @@ -46,6 +57,9 @@ namespace grb { // return ( x > 0 ) - ( x < 0 ); // } // }; + inline int sign(double x) { + return (x > 0) - (x < 0); + } // Custom unary operator for hard clipping to [-1,1] template< typename IOType > @@ -63,37 +77,38 @@ namespace grb { // }; /*-------------------------------------------------------------* - * bSB — core optimisation routine * - *-------------------------------------------------------------*/ - template< - Descriptor descr = descriptors::no_operation, - typename IsingHType, - typename IOType, - typename RSI, typename NZI, Backend backend, - class Ring = Semiring< - grb::operators::add< IOType >, grb::operators::mul< IOType >, - grb::identities::zero, grb::identities::one - >, - class Minus = operators::subtract< IOType >, - class Divide = operators::divide< IOType > - > - grb::RC bSB( - std::vector< IOType > & energies, // output length num_iters - grb::Vector< IOType, backend > & x_comp, // in/out, size N - grb::Vector< IOType, backend > & y_comp, // in/out, size N - const grb::Matrix< IsingHType, backend, RSI, RSI, NZI > & J, // NxN, symmetric - const grb::Vector< IOType, backend > & h, // size N - const IOType p_init, - const IOType p_end, - const std::size_t num_iters, - const IOType dt, - // default semiring, minus, divide - const Ring &ring = Ring(), - const Minus &minus = Minus(), - const Divide ÷ = Divide(), - const std::function< IOType( IOType ) > &sqrtX = - std_sqrt< IOType, IOType > - ) { + * bSB — core optimisation routine * + *-------------------------------------------------------------*/ + template< Descriptor descr = descriptors::no_operation, + typename IsingHType, + typename IOType, + typename RSI, + typename NZI, + Backend backend, + class Ring = Semiring< + grb::operators::add< IOType >, + grb::operators::mul< IOType >, + grb::identities::zero, + grb::identities::one + >, + class Minus = operators::subtract< IOType >, + class Divide = operators::divide< IOType > + > + grb::RC bSB( std::vector< IOType > & energies, // output length num_iters + grb::Vector< IOType, backend > & x_comp, // in/out, size N + grb::Vector< IOType, backend > & y_comp, // in/out, size N + const grb::Matrix< IsingHType, backend, RSI, RSI, NZI > & J, // NxN, symmetric + // TODO: make h const + grb::Vector< IOType, backend > & h, // size N + const IOType p_init, + const IOType p_end, + const std::size_t num_iters, + const IOType dt, + // default semiring, minus, divide + const Ring & ring = Ring(), + const Minus & minus = Minus(), + const Divide & divide = Divide(), + const std::function< IOType( IOType ) > & sqrtX = std_sqrt< IOType, IOType > ) { constexpr const Descriptor descr_dense = descr | descriptors::dense; @@ -109,6 +124,17 @@ namespace grb { // TODO: check that J is symmetric once properly implemented //assert( grb::is_symmetric(J) ); + // print pinned vector x_comp + // for debugging purposes, print x_comp +//#ifdef DEBUG + vector_print( x_comp, "x_comp" ); + vector_print( y_comp, "y_comp" ); + vector_print( h, "h" ); +//#endif + + // assert that energies is of length num_iters + assert( energies.size() == num_iters ); + if ( num_iters < 1 ) { return grb::SUCCESS; } @@ -126,20 +152,20 @@ namespace grb { std::cerr << "Error in eWiseLambda for sumJ2: " << rc << '\n'; return rc; } -#ifdef DEBUG +//#ifdef DEBUG // for debugging purposes, print sumJ2 std::cout << "sumJ2: " << sumJ2 << '\n'; -#endif +//#endif // rewrite this to use graphblas language rc = rc ? rc : grb::foldl( sumJ2, static_cast( N - 1 ), divide ); IOType xi = 0.5; sumJ2 = sqrtX( sumJ2 ); rc = rc ? rc : grb::foldl( xi, sumJ2, divide ); -#ifdef DEBUG +//#ifdef DEBUG // for debugging purposes, print xi std::cout << "xi: " << xi << '\n'; -#endif +//#endif /* ---- iteration variables ---- */ IOType ps = p_init; @@ -151,6 +177,12 @@ namespace grb { // Vec Jx( N ), temp( N ), mask( N ); // workspace vectors grb::Vector< IOType, backend > Jx( N ), temp( N ); grb::Vector< bool, backend > mask( N ); + grb::Vector< IOType, backend > sol( N ); + + grb::set( Jx, ring.template getZero< IOType >() ); + grb::set( temp, ring.template getZero< IOType >() ); + grb::set( mask, false ); + grb::set( sol, 0 ); for ( std::size_t iter = 0; iter < num_iters; ++iter ) { @@ -159,31 +191,51 @@ namespace grb { // Jx ← J * x_comp rc = rc ? rc : grb::mxv< descr_dense >( Jx, J, x_comp, ring ); assert( rc == grb::SUCCESS ); +//#ifdef DEBUG + vector_print( Jx, "Jx" ); +//#endif // temp ← Jx + h rc = rc ? rc : grb::eWiseApply< descr_dense >( temp, Jx, h, ring.getAdditiveMonoid() ); assert( rc == grb::SUCCESS ); +//#ifdef DEBUG + vector_print( temp, "temp = Jx + h" ); +//#endif // temp ← xi * temp rc = rc ? rc : grb::foldl< descr_dense >( temp, xi, ring.getMultiplicativeMonoid() ); assert( rc == grb::SUCCESS ); +//#ifdef DEBUG + vector_print( temp, "xi * temp" ); +//#endif + // temp ← temp + (-1+ps) * x_comp const IOType scale = -1.0 + ps; rc = rc ? rc : grb::eWiseMul< descr_dense >( temp, scale, x_comp, ring ); assert( rc == grb::SUCCESS ); +//#ifdef DEBUG + std::cout << "scale: " << scale << '\n'; + vector_print( temp, "temp + (-1+ps) * x_comp" ); +//#endif // y_comp += dt * temp rc = rc ? rc : grb::eWiseMul< descr_dense >( y_comp, dt, temp, ring ); assert( rc == grb::SUCCESS ); +//#ifdef DEBUG + vector_print( y_comp, "y_comp" ); +//#endif /* x_comp += dt * y_comp */ rc = rc ? rc : grb::eWiseMul< descr_dense >( x_comp, dt, y_comp, ring ); assert( rc == grb::SUCCESS ); +//#ifdef DEBUG + vector_print( x_comp, "x_comp" ); +//#endif /* y_comp[ |x|>1 ] = 0 */ // mask = np.abs(x_comp) > 1 @@ -194,6 +246,9 @@ namespace grb { }, mask ); assert( rc == grb::SUCCESS ); +//#ifdef DEBUG + vector_print( mask, "mask" ); +//#endif // y_comp[ mask ] = 0 rc = rc ? rc : grb::eWiseLambda< descr_dense >( [&mask, &y_comp]( const size_t i ) { @@ -207,29 +262,60 @@ namespace grb { assert( rc == grb::SUCCESS ); /* x_comp = clip( x_comp ) */ - rc = rc ? rc : grb::eWiseLambda< descr_dense >( [&x_comp]( const size_t i ) { - (void) i; - x_comp[i] = clip11()( x_comp[i] ); }, x_comp + rc = rc ? rc : grb::eWiseLambda< descr_dense >( + [&x_comp]( const size_t i ) { + (void) i; + // TODO: rewrite in terms of foldl and foldr + x_comp[i] = std::min( 1.0, std::max( -1.0, x_comp[i] ) ); + }, + x_comp ); assert( rc == grb::SUCCESS ); - // /* Energy evaluation */ - // Vec sol( N ); - // GRB_TRY( apply( sol, x_comp, signum() ) ); - - // // temp ← J * sol - // GRB_TRY( mxv( temp, J, sol ) ); + /* Energy evaluation */ + // sol[i] = sign(x_comp[i]); which in graphblas is: + + + rc = rc ? rc : grb::eWiseLambda< descr_dense >( + [&sol,&x_comp]( const size_t i ) { + (void) i; + // TODO: rewrite in terms of foldl and foldr + sol[i] = sign(x_comp[i]); + }, + sol + ); + assert( rc == grb::SUCCESS ); +//#ifdef DEBUG + vector_print( sol, "sol" ); +//#endif - // // e = -0.5 * sol^T * temp – h^T * sol - // IOType dot1 = 0.0, dot2 = 0.0; - // GRB_TRY( foldl( dot1, sol, temp, operators::add< IOType >(), - // operators::mul< IOType >() ) ); - // GRB_TRY( foldl( dot2, h, sol, operators::add< IOType >(), - // operators::mul< IOType >() ) ); - energies[ iter ] = iter; // just for testing purposes - // energies[ iter ] = -0.5 * dot1 - dot2; + // temp ← J * sol + rc = rc ? rc : grb::set( temp, ring.template getZero< IOType >() ); + rc = rc ? rc : grb::mxv< descr_dense >( temp, J, sol, ring ); + assert( rc == grb::SUCCESS ); +//#ifdef DEBUG + vector_print( temp, "temp = J * sol" ); +//#endif + // e = -0.5 * sol.dot(temp) – h.dot(sol) + IOType dot1 = 0.0, dot2 = 0.0; + rc = rc ? rc : grb::dot< descr_dense >( dot1, sol, temp, ring ); + assert( rc == grb::SUCCESS ); +//#ifdef DEBUG + std::cout << "dot1: " << dot1 << '\n'; +//#endif + rc = rc ? rc : grb::dot< descr_dense >( dot2, h, sol, ring ); + assert( rc == grb::SUCCESS ); +//#ifdef DEBUG + std::cout << "dot2: " << dot2 << '\n'; +//#endif + + IOType e = -0.5 * dot1 - dot2; +//#ifdef DEBUG + std::cout << "e: " << e << '\n'; +//#endif + energies[ iter ] = e; ps += dps; } diff --git a/tests/smoke/ising_machine_sb.cpp b/tests/smoke/ising_machine_sb.cpp index 89d72d296..3397ee2a9 100644 --- a/tests/smoke/ising_machine_sb.cpp +++ b/tests/smoke/ising_machine_sb.cpp @@ -41,12 +41,22 @@ using namespace grb; using namespace algorithms; +// test data from python implementation +constexpr std::size_t N = 10; +constexpr std::size_t Nz = 54; +static const size_t i_arr[ Nz ] ={0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 9, 9}; +static const size_t j_arr[ Nz ] ={0, 2, 3, 4, 5, 1, 4, 5, 6, 7, 9, 0, 2, 4, 6, 9, 0, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 0, 1, 3, 5, 6, 8, 1, 2, 3, 5, 6, 1, 3, 7, 8, 9, 3, 5, 7, 8, 1, 2, 3, 7, 9}; +static const double v_arr[ Nz ]={-1, 1, 1, -1, -1, -1, -1, 1, 1, 1, 1, 1, -1, -1, -1, -1, 1,1, 1, 1, -1, -1, 1, 1, -1, -1, -1, 1, 1, -1, 1, 1, -1, -1,-1, 1, -1, -1, -1, -1, 1, -1, -1, 1, -1, 1, -1, 1, 1, 1, -1,1, -1, 1}; + +static const double h_arr[ N ]={ 1, -1, 1, -1, 1, 1, -1, 1, 1, 1}; +static const double x_arr[ N ]={-0.0996, -0.0315, 0.0572, 0.0630, 0.0087, -0.0143, -0.0170, -0.0411, 0.0433, -0.0298}; +static const double y_arr[ N ]={ 0.0373, 0.0540, 0.0486, -0.0877, -0.0418, -0.0261, 0.0018, -0.0710, 0.0507, -0.0483}; using IOType = double; +using JType = double; + int main( int argc, char ** argv ) { - constexpr std::size_t N = 100; - constexpr std::size_t Nz = 400; /* --- Initialise ALP/GraphBLAS --- */ if( grb::init() != grb::SUCCESS ) { @@ -58,56 +68,28 @@ int main( int argc, char ** argv ) { grb::Matrix J( N, N, Nz ); grb::Vector h( N ); grb::Vector x0( N ), y0( N ); // initialy - // ... populate J with random values - // generate random nonzero locations and values for J - std::vector row_indices(Nz); - std::vector col_indices(Nz); - std::vector nonzero_values(Nz); - // TODO: make sure there are not duplicate entries in J and J is symmetric - - IOType sumJ2_test = 0; - for( std::size_t i = 0; i < Nz; ++i ) { - row_indices[i] = rand() % N; - col_indices[i] = rand() % N; - nonzero_values[i] = static_cast(rand()) / RAND_MAX; // random value between 0 and 1 - sumJ2_test+=nonzero_values[i]*nonzero_values[i]; - } -#ifdef DEBUG - // for debugging purposes, print sumJ2_test - std::cout << "sumJ2_test: " << sumJ2_test << '\n'; -#endif - IOType xi_test = 0.5 / std::sqrt( sumJ2_test / static_cast( N - 1 ) ); -#ifdef DEBUG - // for debugging purposes, print xi_test - std::cout << "xi_test: " << xi_test << '\n'; -#endif - + // ... populate J with test (random) values grb::RC rc = grb::SUCCESS; - rc = rc ? rc : buildMatrixUnique(J, row_indices.data(), col_indices.data(), nonzero_values.data(), Nz, grb::SEQUENTIAL); + rc = rc ? rc : buildMatrixUnique( J, &( i_arr[ 0 ] ), &( j_arr[ 0 ] ), &( v_arr[ 0 ] ), Nz, grb::SEQUENTIAL ); + if(rc != grb::SUCCESS) { std::cerr << "matrix build failed\n"; return grb::RC::PANIC; } // Fill h, x0, y0 with random values using buildVector - std::vector h_buffer(N), x0_buffer(N), y0_buffer(N); - for( std::size_t i = 0; i < N; ++i ) { - h_buffer[i] = static_cast(rand()) / RAND_MAX; - x0_buffer[i] = static_cast(rand()) / RAND_MAX / 0.2 - 0.1; - y0_buffer[i] = static_cast(rand()) / RAND_MAX / 0.2 - 0.1; - } - rc = rc ? rc : buildVector(h, h_buffer.begin(), h_buffer.end(), grb::SEQUENTIAL); - rc = rc ? rc : buildVector(x0, x0_buffer.begin(), x0_buffer.end(), grb::SEQUENTIAL); - rc = rc ? rc : buildVector(y0, y0_buffer.begin(), y0_buffer.end(), grb::SEQUENTIAL); + rc = rc ? rc : buildVector(h, h_arr, h_arr + N, grb::SEQUENTIAL); + rc = rc ? rc : buildVector(x0, x_arr, x_arr + N, grb::SEQUENTIAL); + rc = rc ? rc : buildVector(y0, y_arr, y_arr + N, grb::SEQUENTIAL); if(rc != grb::SUCCESS) { std::cerr << "Vector build failed\n"; return grb::RC::PANIC; } - - const IOType p0 = 0.1; - const IOType p1 = 1.0; - const std::size_t num_iters = 1000; - const IOType dt = 0.01; + + const IOType p0 = 0.; + const IOType p1 = 1.1; + const std::size_t num_iters = 10; + const IOType dt = 0.25; // energies is array of length num_iters, initialized to 0 std::vector< IOType > energies( num_iters, 0 ); @@ -119,7 +101,10 @@ int main( int argc, char ** argv ) { if( rc != grb::SUCCESS ) { std::cerr << "bSB returned error code " << rc << '\n'; } else { - std::cout << "Final energy = " << energies[num_iters-1] << '\n'; + // print all energies + for (std::size_t i = 0; i < num_iters; ++i) { + std::cout << "Energy at iteration " << i << " = " << energies[i] << '\n'; + } } grb::finalize(); From ff50194b5f97b8d4982dcc4c01800a1f2a42d8e6 Mon Sep 17 00:00:00 2001 From: Denis Jelovina Date: Wed, 16 Jul 2025 15:44:01 +0200 Subject: [PATCH 04/18] Enhance Ising machine SB implementation with workspace vector initialization and update test cases for energy validation --- .../graphblas/algorithms/ising_machine_sb.hpp | 127 ++++++++++-------- tests/smoke/ising_machine_sb.cpp | 52 +++++-- 2 files changed, 110 insertions(+), 69 deletions(-) diff --git a/include/graphblas/algorithms/ising_machine_sb.hpp b/include/graphblas/algorithms/ising_machine_sb.hpp index ef4666376..eb6ec968d 100644 --- a/include/graphblas/algorithms/ising_machine_sb.hpp +++ b/include/graphblas/algorithms/ising_machine_sb.hpp @@ -42,12 +42,13 @@ namespace grb { void vector_print( grb::Vector< IOType, backend > & x_comp, const std::string & vector_name ) { grb::PinnedVector< IOType > pinnedVector; pinnedVector = grb::PinnedVector< IOType >( x_comp, grb::SEQUENTIAL ); - std::cout << "First 10 nonzeroes of " << vector_name << " are: ( "; + std::cout << "First 10 nonzeroes of " << vector_name << " = [ "; for( size_t k = 0; k < pinnedVector.nonzeroes() && k < 10; ++k ) { const IOType & nonzeroValue = pinnedVector.getNonzeroValue( k ); - std::cout << nonzeroValue << " "; + std::cout << nonzeroValue; + if(k!=pinnedVector.nonzeroes()-1) std::cout << ", "; } - std::cout << ")" << std::endl; + std::cout << "]" << std::endl; } @@ -104,6 +105,11 @@ namespace grb { const IOType p_end, const std::size_t num_iters, const IOType dt, + // workspace vectors + grb::Vector< IOType, backend > & Jx, + grb::Vector< IOType, backend > & temp, + grb::Vector< bool, backend > & mask, + grb::Vector< IOType, backend > & sol, // default semiring, minus, divide const Ring & ring = Ring(), const Minus & minus = Minus(), @@ -124,13 +130,19 @@ namespace grb { // TODO: check that J is symmetric once properly implemented //assert( grb::is_symmetric(J) ); + // initialize workspace vectors + grb::set( Jx, ring.template getZero< IOType >() ); + grb::set( temp, ring.template getZero< IOType >() ); + grb::set( mask, false ); + grb::set( sol, ring.template getZero< IOType >() ); + // print pinned vector x_comp // for debugging purposes, print x_comp -//#ifdef DEBUG +#ifdef DEBUG vector_print( x_comp, "x_comp" ); vector_print( y_comp, "y_comp" ); vector_print( h, "h" ); -//#endif +#endif // assert that energies is of length num_iters assert( energies.size() == num_iters ); @@ -142,30 +154,31 @@ namespace grb { grb::RC rc = grb::SUCCESS; /* ---- pre-compute ---- */ IOType sumJ2 = ring.template getZero< IOType >(); - rc = rc ? rc : grb::eWiseLambda( [&J, &sumJ2]( const size_t i, const size_t j, IOType& v ) { + rc = rc ? rc : grb::eWiseLambda( [&J, &sumJ2, &ring]( const size_t i, const size_t j, IOType& v ) { (void) i; (void) j; - // rewrite this to use graphblas language - sumJ2 += v*v; + IOType v2; + apply( v2, v, v, ring.getMultiplicativeOperator() ); + foldl( sumJ2, v2, ring.getAdditiveOperator() ); }, J ); if( rc != grb::SUCCESS ) { std::cerr << "Error in eWiseLambda for sumJ2: " << rc << '\n'; return rc; } -//#ifdef DEBUG +#ifdef DEBUG // for debugging purposes, print sumJ2 std::cout << "sumJ2: " << sumJ2 << '\n'; -//#endif +#endif // rewrite this to use graphblas language rc = rc ? rc : grb::foldl( sumJ2, static_cast( N - 1 ), divide ); IOType xi = 0.5; sumJ2 = sqrtX( sumJ2 ); rc = rc ? rc : grb::foldl( xi, sumJ2, divide ); -//#ifdef DEBUG +#ifdef DEBUG // for debugging purposes, print xi std::cout << "xi: " << xi << '\n'; -//#endif +#endif /* ---- iteration variables ---- */ IOType ps = p_init; @@ -173,69 +186,60 @@ namespace grb { // assert len of energies == N assert( energies.size() == num_iters ); - // TODO: move these to the aggument list of bSB - // Vec Jx( N ), temp( N ), mask( N ); // workspace vectors - grb::Vector< IOType, backend > Jx( N ), temp( N ); - grb::Vector< bool, backend > mask( N ); - grb::Vector< IOType, backend > sol( N ); - - grb::set( Jx, ring.template getZero< IOType >() ); - grb::set( temp, ring.template getZero< IOType >() ); - grb::set( mask, false ); - grb::set( sol, 0 ); - for ( std::size_t iter = 0; iter < num_iters; ++iter ) { /* y_comp += ((-1+ps)*x_comp + xi*(Jx + h)) * dt */ // Jx ← J * x_comp + grb::set( Jx, ring.template getZero< IOType >() ); rc = rc ? rc : grb::mxv< descr_dense >( Jx, J, x_comp, ring ); assert( rc == grb::SUCCESS ); -//#ifdef DEBUG +#ifdef DEBUG vector_print( Jx, "Jx" ); -//#endif +#endif // temp ← Jx + h + grb::set( temp, ring.template getZero< IOType >() ); rc = rc ? rc : grb::eWiseApply< descr_dense >( temp, Jx, h, ring.getAdditiveMonoid() ); assert( rc == grb::SUCCESS ); -//#ifdef DEBUG +#ifdef DEBUG vector_print( temp, "temp = Jx + h" ); -//#endif +#endif // temp ← xi * temp rc = rc ? rc : grb::foldl< descr_dense >( temp, xi, ring.getMultiplicativeMonoid() ); assert( rc == grb::SUCCESS ); -//#ifdef DEBUG +#ifdef DEBUG vector_print( temp, "xi * temp" ); -//#endif +#endif // temp ← temp + (-1+ps) * x_comp const IOType scale = -1.0 + ps; rc = rc ? rc : grb::eWiseMul< descr_dense >( temp, scale, x_comp, ring ); assert( rc == grb::SUCCESS ); -//#ifdef DEBUG +#ifdef DEBUG std::cout << "scale: " << scale << '\n'; vector_print( temp, "temp + (-1+ps) * x_comp" ); -//#endif +#endif // y_comp += dt * temp rc = rc ? rc : grb::eWiseMul< descr_dense >( y_comp, dt, temp, ring ); assert( rc == grb::SUCCESS ); -//#ifdef DEBUG +#ifdef DEBUG vector_print( y_comp, "y_comp" ); -//#endif +#endif /* x_comp += dt * y_comp */ rc = rc ? rc : grb::eWiseMul< descr_dense >( x_comp, dt, y_comp, ring ); assert( rc == grb::SUCCESS ); -//#ifdef DEBUG +#ifdef DEBUG vector_print( x_comp, "x_comp" ); -//#endif +#endif /* y_comp[ |x|>1 ] = 0 */ // mask = np.abs(x_comp) > 1 @@ -246,9 +250,9 @@ namespace grb { }, mask ); assert( rc == grb::SUCCESS ); -//#ifdef DEBUG +#ifdef DEBUG vector_print( mask, "mask" ); -//#endif +#endif // y_comp[ mask ] = 0 rc = rc ? rc : grb::eWiseLambda< descr_dense >( [&mask, &y_comp]( const size_t i ) { @@ -262,59 +266,64 @@ namespace grb { assert( rc == grb::SUCCESS ); /* x_comp = clip( x_comp ) */ - rc = rc ? rc : grb::eWiseLambda< descr_dense >( - [&x_comp]( const size_t i ) { - (void) i; - // TODO: rewrite in terms of foldl and foldr - x_comp[i] = std::min( 1.0, std::max( -1.0, x_comp[i] ) ); - }, - x_comp - ); + foldl( x_comp, static_cast(-1), grb::operators::max < IOType >() ); + foldl( x_comp, static_cast(1), grb::operators::min < IOType >() ); + // alternatively, we could use eWiseLambda: + // rc = rc ? rc : grb::eWiseLambda< descr_dense >( + // [&x_comp]( const size_t i ) { + // (void) i; + // foldl( x_comp[i], static_cast(-1), grb::operators::max < IOType >() ); + // foldl( x_comp[i], static_cast(1), grb::operators::min < IOType >() ); + // //x_comp[i] = std::min( 1.0, std::max( -1.0, x_comp[i] ) ); + // }, x_comp + // ); assert( rc == grb::SUCCESS ); +//#ifdef DEBUG + std::cout << "i = " << iter << "\n "; + vector_print( x_comp, "x_comp_alp " ); + vector_print( y_comp, "y_comp_alp" ); +//#endif /* Energy evaluation */ - // sol[i] = sign(x_comp[i]); which in graphblas is: - - + // sol[i] = sign(x_comp[i]); which in graphblas is: rc = rc ? rc : grb::eWiseLambda< descr_dense >( [&sol,&x_comp]( const size_t i ) { (void) i; - // TODO: rewrite in terms of foldl and foldr sol[i] = sign(x_comp[i]); }, sol ); assert( rc == grb::SUCCESS ); -//#ifdef DEBUG +#ifdef DEBUG vector_print( sol, "sol" ); -//#endif +#endif // temp ← J * sol rc = rc ? rc : grb::set( temp, ring.template getZero< IOType >() ); rc = rc ? rc : grb::mxv< descr_dense >( temp, J, sol, ring ); assert( rc == grb::SUCCESS ); -//#ifdef DEBUG +#ifdef DEBUG vector_print( temp, "temp = J * sol" ); -//#endif +#endif // e = -0.5 * sol.dot(temp) – h.dot(sol) IOType dot1 = 0.0, dot2 = 0.0; rc = rc ? rc : grb::dot< descr_dense >( dot1, sol, temp, ring ); assert( rc == grb::SUCCESS ); -//#ifdef DEBUG +#ifdef DEBUG std::cout << "dot1: " << dot1 << '\n'; -//#endif +#endif rc = rc ? rc : grb::dot< descr_dense >( dot2, h, sol, ring ); assert( rc == grb::SUCCESS ); -//#ifdef DEBUG +#ifdef DEBUG std::cout << "dot2: " << dot2 << '\n'; -//#endif +#endif IOType e = -0.5 * dot1 - dot2; -//#ifdef DEBUG +#ifdef DEBUG std::cout << "e: " << e << '\n'; -//#endif +#endif energies[ iter ] = e; ps += dps; } diff --git a/tests/smoke/ising_machine_sb.cpp b/tests/smoke/ising_machine_sb.cpp index 3397ee2a9..620d7e315 100644 --- a/tests/smoke/ising_machine_sb.cpp +++ b/tests/smoke/ising_machine_sb.cpp @@ -44,13 +44,34 @@ using namespace algorithms; // test data from python implementation constexpr std::size_t N = 10; constexpr std::size_t Nz = 54; -static const size_t i_arr[ Nz ] ={0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 9, 9}; -static const size_t j_arr[ Nz ] ={0, 2, 3, 4, 5, 1, 4, 5, 6, 7, 9, 0, 2, 4, 6, 9, 0, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 0, 1, 3, 5, 6, 8, 1, 2, 3, 5, 6, 1, 3, 7, 8, 9, 3, 5, 7, 8, 1, 2, 3, 7, 9}; -static const double v_arr[ Nz ]={-1, 1, 1, -1, -1, -1, -1, 1, 1, 1, 1, 1, -1, -1, -1, -1, 1,1, 1, 1, -1, -1, 1, 1, -1, -1, -1, 1, 1, -1, 1, 1, -1, -1,-1, 1, -1, -1, -1, -1, 1, -1, -1, 1, -1, 1, -1, 1, 1, 1, -1,1, -1, 1}; - -static const double h_arr[ N ]={ 1, -1, 1, -1, 1, 1, -1, 1, 1, 1}; -static const double x_arr[ N ]={-0.0996, -0.0315, 0.0572, 0.0630, 0.0087, -0.0143, -0.0170, -0.0411, 0.0433, -0.0298}; -static const double y_arr[ N ]={ 0.0373, 0.0540, 0.0486, -0.0877, -0.0418, -0.0261, 0.0018, -0.0710, 0.0507, -0.0483}; +static const size_t i_arr[ Nz ] = { + 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, + 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 9, 9 +}; +static const size_t j_arr[ Nz ] = { + 0, 2, 3, 4, 5, 1, 4, 5, 6, 7, 9, 0, 2, 4, 6, 9, 0, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, + 3, 4, 0, 1, 3, 5, 6, 8, 1, 2, 3, 5, 6, 1, 3, 7, 8, 9, 3, 5, 7, 8, 1, 2, 3, 7, 9 +}; +static const double v_arr[ Nz ] = { + -1, 1, 1, -1, -1, -1, -1, 1, 1, 1, 1, 1, -1, -1, -1, -1, 1, 1, 1, 1, -1, -1, 1, 1, -1, -1, -1, + 1, 1, -1, 1, 1, -1, -1, -1, 1, -1, -1, -1, -1, 1, -1, -1, 1, -1, 1, -1, 1, 1, 1, -1, 1, -1, 1 +}; + +static const double h_arr[ N ] = { 1, -1, 1, -1, 1, 1, -1, 1, 1, 1 }; +static const double x_arr[ N ] = { -0.0996, -0.0315, 0.0572, 0.0630, 0.0087, -0.0143, -0.0170, -0.0411, 0.0433, -0.0298 }; +static const double y_arr[ N ] = { 0.0373, 0.0540, 0.0486, -0.0877, -0.0418, -0.0261, 0.0018, -0.0710, 0.0507, -0.0483 }; + +const std::size_t num_iters = 100; +static const double energies_ref[ num_iters ] = { + -3, -3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, + -1, -9, -11, -11, -11, -11, -11, -11, -11, -11, -9, -9, -11, + -9, -9, -9, -9, -9, -9, -9, -9, -9, -9, -9, -9, -9, + -9, -9, -11, -11, -11, -9, -9, -9, -9, -9, -9, -9, -9, + -9, -9, -9, -9, -9, -9, -9, -11, -11, -11, -11, -11, -11, + -11, -11, -11, -11, -11, -11, -11, -11, -11, -11, -11, -11, -11, + -11, -11, -11, -11, -11, -11, -11, -11, -11, -11, -11, -11, -11, + -11, -11, -11, -11, -11, -11, -11, -11, -11 +}; using IOType = double; using JType = double; @@ -85,25 +106,36 @@ int main( int argc, char ** argv ) { std::cerr << "Vector build failed\n"; return grb::RC::PANIC; } - + const IOType p0 = 0.; const IOType p1 = 1.1; - const std::size_t num_iters = 10; const IOType dt = 0.25; // energies is array of length num_iters, initialized to 0 std::vector< IOType > energies( num_iters, 0 ); + grb::Vector< IOType > Jx( N ); + grb::Vector< IOType > temp( N ); + grb::Vector< bool > mask( N ); + grb::Vector< IOType > sol( N ); + rc = rc ? rc : bSB( - energies, x0, y0, J, h, p0, p1, num_iters, dt + energies, x0, y0, J, h, p0, p1, num_iters, dt, + Jx, temp, mask, sol ); + if( rc != grb::SUCCESS ) { std::cerr << "bSB returned error code " << rc << '\n'; } else { // print all energies for (std::size_t i = 0; i < num_iters; ++i) { std::cout << "Energy at iteration " << i << " = " << energies[i] << '\n'; + if( energies[i] != energies_ref[i]) { + std::cerr << "Error: Energy at iteration " << i << " does not match reference value.\n"; + std::cerr << "Expected: " << energies_ref[i] << ", got: " << energies[i] << '\n'; + return -1; + } } } From 83caeb586442a971d9e93ba18e5f7102b9acd2c7 Mon Sep 17 00:00:00 2001 From: Denis Jelovina Date: Wed, 16 Jul 2025 16:06:24 +0200 Subject: [PATCH 05/18] Refactor Ising machine SB implementation to use template types for sign function and update test cases to use integer types for matrices and vectors --- .../graphblas/algorithms/ising_machine_sb.hpp | 36 +++++++++---------- tests/smoke/ising_machine_sb.cpp | 18 ++++++---- 2 files changed, 29 insertions(+), 25 deletions(-) diff --git a/include/graphblas/algorithms/ising_machine_sb.hpp b/include/graphblas/algorithms/ising_machine_sb.hpp index eb6ec968d..350d07352 100644 --- a/include/graphblas/algorithms/ising_machine_sb.hpp +++ b/include/graphblas/algorithms/ising_machine_sb.hpp @@ -58,7 +58,8 @@ namespace grb { // return ( x > 0 ) - ( x < 0 ); // } // }; - inline int sign(double x) { + template< typename IType, typename ReturnType > + inline ReturnType sign(IType x) { return (x > 0) - (x < 0); } @@ -100,7 +101,7 @@ namespace grb { grb::Vector< IOType, backend > & y_comp, // in/out, size N const grb::Matrix< IsingHType, backend, RSI, RSI, NZI > & J, // NxN, symmetric // TODO: make h const - grb::Vector< IOType, backend > & h, // size N + grb::Vector< IsingHType, backend > & h, // size N const IOType p_init, const IOType p_end, const std::size_t num_iters, @@ -108,8 +109,9 @@ namespace grb { // workspace vectors grb::Vector< IOType, backend > & Jx, grb::Vector< IOType, backend > & temp, + grb::Vector< IsingHType, backend > & temp_int, grb::Vector< bool, backend > & mask, - grb::Vector< IOType, backend > & sol, + grb::Vector< IsingHType, backend > & sol, // default semiring, minus, divide const Ring & ring = Ring(), const Minus & minus = Minus(), @@ -130,11 +132,8 @@ namespace grb { // TODO: check that J is symmetric once properly implemented //assert( grb::is_symmetric(J) ); - // initialize workspace vectors - grb::set( Jx, ring.template getZero< IOType >() ); - grb::set( temp, ring.template getZero< IOType >() ); - grb::set( mask, false ); - grb::set( sol, ring.template getZero< IOType >() ); + grb::set( sol, ring.template getZero< IsingHType >() ); + grb::set( mask, ring.template getZero< bool >() ); // print pinned vector x_comp // for debugging purposes, print x_comp @@ -154,10 +153,10 @@ namespace grb { grb::RC rc = grb::SUCCESS; /* ---- pre-compute ---- */ IOType sumJ2 = ring.template getZero< IOType >(); - rc = rc ? rc : grb::eWiseLambda( [&J, &sumJ2, &ring]( const size_t i, const size_t j, IOType& v ) { + rc = rc ? rc : grb::eWiseLambda( [&J, &sumJ2, &ring]( const size_t i, const size_t j, IsingHType& v ) { (void) i; (void) j; - IOType v2; + IsingHType v2; apply( v2, v, v, ring.getMultiplicativeOperator() ); foldl( sumJ2, v2, ring.getAdditiveOperator() ); }, J ); @@ -278,18 +277,18 @@ namespace grb { // }, x_comp // ); assert( rc == grb::SUCCESS ); -//#ifdef DEBUG +#ifdef DEBUG std::cout << "i = " << iter << "\n "; vector_print( x_comp, "x_comp_alp " ); vector_print( y_comp, "y_comp_alp" ); -//#endif +#endif /* Energy evaluation */ // sol[i] = sign(x_comp[i]); which in graphblas is: rc = rc ? rc : grb::eWiseLambda< descr_dense >( [&sol,&x_comp]( const size_t i ) { (void) i; - sol[i] = sign(x_comp[i]); + sol[i] = sign(x_comp[i]); }, sol ); @@ -300,15 +299,16 @@ namespace grb { // temp ← J * sol - rc = rc ? rc : grb::set( temp, ring.template getZero< IOType >() ); - rc = rc ? rc : grb::mxv< descr_dense >( temp, J, sol, ring ); + rc = rc ? rc : grb::set( temp_int, ring.template getZero< IsingHType >() ); + rc = rc ? rc : grb::mxv< descr_dense >( temp_int, J, sol, ring ); assert( rc == grb::SUCCESS ); #ifdef DEBUG - vector_print( temp, "temp = J * sol" ); + vector_print( temp_int, "temp = J * sol" ); #endif // e = -0.5 * sol.dot(temp) – h.dot(sol) - IOType dot1 = 0.0, dot2 = 0.0; - rc = rc ? rc : grb::dot< descr_dense >( dot1, sol, temp, ring ); + IsingHType dot1 = 0; + IOType dot2 = 0; + rc = rc ? rc : grb::dot< descr_dense >( dot1, sol, temp_int, ring ); assert( rc == grb::SUCCESS ); #ifdef DEBUG std::cout << "dot1: " << dot1 << '\n'; diff --git a/tests/smoke/ising_machine_sb.cpp b/tests/smoke/ising_machine_sb.cpp index 620d7e315..e32d29aff 100644 --- a/tests/smoke/ising_machine_sb.cpp +++ b/tests/smoke/ising_machine_sb.cpp @@ -52,12 +52,12 @@ static const size_t j_arr[ Nz ] = { 0, 2, 3, 4, 5, 1, 4, 5, 6, 7, 9, 0, 2, 4, 6, 9, 0, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 0, 1, 3, 5, 6, 8, 1, 2, 3, 5, 6, 1, 3, 7, 8, 9, 3, 5, 7, 8, 1, 2, 3, 7, 9 }; -static const double v_arr[ Nz ] = { +static const int v_arr[ Nz ] = { -1, 1, 1, -1, -1, -1, -1, 1, 1, 1, 1, 1, -1, -1, -1, -1, 1, 1, 1, 1, -1, -1, 1, 1, -1, -1, -1, 1, 1, -1, 1, 1, -1, -1, -1, 1, -1, -1, -1, -1, 1, -1, -1, 1, -1, 1, -1, 1, 1, 1, -1, 1, -1, 1 }; -static const double h_arr[ N ] = { 1, -1, 1, -1, 1, 1, -1, 1, 1, 1 }; +static const int h_arr[ N ] = { 1, -1, 1, -1, 1, 1, -1, 1, 1, 1 }; static const double x_arr[ N ] = { -0.0996, -0.0315, 0.0572, 0.0630, 0.0087, -0.0143, -0.0170, -0.0411, 0.0433, -0.0298 }; static const double y_arr[ N ] = { 0.0373, 0.0540, 0.0486, -0.0877, -0.0418, -0.0261, 0.0018, -0.0710, 0.0507, -0.0483 }; @@ -74,7 +74,7 @@ static const double energies_ref[ num_iters ] = { }; using IOType = double; -using JType = double; +using JType = int; int main( int argc, char ** argv ) { @@ -86,8 +86,8 @@ int main( int argc, char ** argv ) { } /* --- Problem setup (toy) --- */ - grb::Matrix J( N, N, Nz ); - grb::Vector h( N ); + grb::Matrix J( N, N, Nz ); + grb::Vector h( N ); grb::Vector x0( N ), y0( N ); // initialy // ... populate J with test (random) values grb::RC rc = grb::SUCCESS; @@ -116,12 +116,14 @@ int main( int argc, char ** argv ) { grb::Vector< IOType > Jx( N ); grb::Vector< IOType > temp( N ); + grb::Vector< JType > temp_int( N ); grb::Vector< bool > mask( N ); - grb::Vector< IOType > sol( N ); + // TODO: make sol int type + grb::Vector< JType > sol( N ); rc = rc ? rc : bSB( energies, x0, y0, J, h, p0, p1, num_iters, dt, - Jx, temp, mask, sol + Jx, temp, temp_int, mask, sol ); @@ -137,6 +139,8 @@ int main( int argc, char ** argv ) { return -1; } } + std::cout << "All energies match reference values.\n"; + std::cout << "TEST OK\n"; } grb::finalize(); From e140618cdf09b9b2cb552032ab82d8529bbef635 Mon Sep 17 00:00:00 2001 From: Denis Jelovina Date: Wed, 16 Jul 2025 16:11:53 +0200 Subject: [PATCH 06/18] Update debug output macros in Ising machine SB implementation to use DEBUG_IMSB --- .../graphblas/algorithms/ising_machine_sb.hpp | 32 +++++++++---------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/include/graphblas/algorithms/ising_machine_sb.hpp b/include/graphblas/algorithms/ising_machine_sb.hpp index 350d07352..16b5c1707 100644 --- a/include/graphblas/algorithms/ising_machine_sb.hpp +++ b/include/graphblas/algorithms/ising_machine_sb.hpp @@ -137,7 +137,7 @@ namespace grb { // print pinned vector x_comp // for debugging purposes, print x_comp -#ifdef DEBUG +#ifdef DEBUG_IMSB vector_print( x_comp, "x_comp" ); vector_print( y_comp, "y_comp" ); vector_print( h, "h" ); @@ -164,7 +164,7 @@ namespace grb { std::cerr << "Error in eWiseLambda for sumJ2: " << rc << '\n'; return rc; } -#ifdef DEBUG +#ifdef DEBUG_IMSB // for debugging purposes, print sumJ2 std::cout << "sumJ2: " << sumJ2 << '\n'; #endif @@ -174,7 +174,7 @@ namespace grb { IOType xi = 0.5; sumJ2 = sqrtX( sumJ2 ); rc = rc ? rc : grb::foldl( xi, sumJ2, divide ); -#ifdef DEBUG +#ifdef DEBUG_IMSB // for debugging purposes, print xi std::cout << "xi: " << xi << '\n'; #endif @@ -193,7 +193,7 @@ namespace grb { grb::set( Jx, ring.template getZero< IOType >() ); rc = rc ? rc : grb::mxv< descr_dense >( Jx, J, x_comp, ring ); assert( rc == grb::SUCCESS ); -#ifdef DEBUG +#ifdef DEBUG_IMSB vector_print( Jx, "Jx" ); #endif @@ -203,7 +203,7 @@ namespace grb { temp, Jx, h, ring.getAdditiveMonoid() ); assert( rc == grb::SUCCESS ); -#ifdef DEBUG +#ifdef DEBUG_IMSB vector_print( temp, "temp = Jx + h" ); #endif @@ -212,7 +212,7 @@ namespace grb { temp, xi, ring.getMultiplicativeMonoid() ); assert( rc == grb::SUCCESS ); -#ifdef DEBUG +#ifdef DEBUG_IMSB vector_print( temp, "xi * temp" ); #endif @@ -221,7 +221,7 @@ namespace grb { const IOType scale = -1.0 + ps; rc = rc ? rc : grb::eWiseMul< descr_dense >( temp, scale, x_comp, ring ); assert( rc == grb::SUCCESS ); -#ifdef DEBUG +#ifdef DEBUG_IMSB std::cout << "scale: " << scale << '\n'; vector_print( temp, "temp + (-1+ps) * x_comp" ); #endif @@ -229,14 +229,14 @@ namespace grb { // y_comp += dt * temp rc = rc ? rc : grb::eWiseMul< descr_dense >( y_comp, dt, temp, ring ); assert( rc == grb::SUCCESS ); -#ifdef DEBUG +#ifdef DEBUG_IMSB vector_print( y_comp, "y_comp" ); #endif /* x_comp += dt * y_comp */ rc = rc ? rc : grb::eWiseMul< descr_dense >( x_comp, dt, y_comp, ring ); assert( rc == grb::SUCCESS ); -#ifdef DEBUG +#ifdef DEBUG_IMSB vector_print( x_comp, "x_comp" ); #endif @@ -249,7 +249,7 @@ namespace grb { }, mask ); assert( rc == grb::SUCCESS ); -#ifdef DEBUG +#ifdef DEBUG_IMSB vector_print( mask, "mask" ); #endif @@ -277,7 +277,7 @@ namespace grb { // }, x_comp // ); assert( rc == grb::SUCCESS ); -#ifdef DEBUG +#ifdef DEBUG_IMSB std::cout << "i = " << iter << "\n "; vector_print( x_comp, "x_comp_alp " ); vector_print( y_comp, "y_comp_alp" ); @@ -293,7 +293,7 @@ namespace grb { sol ); assert( rc == grb::SUCCESS ); -#ifdef DEBUG +#ifdef DEBUG_IMSB vector_print( sol, "sol" ); #endif @@ -302,7 +302,7 @@ namespace grb { rc = rc ? rc : grb::set( temp_int, ring.template getZero< IsingHType >() ); rc = rc ? rc : grb::mxv< descr_dense >( temp_int, J, sol, ring ); assert( rc == grb::SUCCESS ); -#ifdef DEBUG +#ifdef DEBUG_IMSB vector_print( temp_int, "temp = J * sol" ); #endif // e = -0.5 * sol.dot(temp) – h.dot(sol) @@ -310,18 +310,18 @@ namespace grb { IOType dot2 = 0; rc = rc ? rc : grb::dot< descr_dense >( dot1, sol, temp_int, ring ); assert( rc == grb::SUCCESS ); -#ifdef DEBUG +#ifdef DEBUG_IMSB std::cout << "dot1: " << dot1 << '\n'; #endif rc = rc ? rc : grb::dot< descr_dense >( dot2, h, sol, ring ); assert( rc == grb::SUCCESS ); -#ifdef DEBUG +#ifdef DEBUG_IMSB std::cout << "dot2: " << dot2 << '\n'; #endif IOType e = -0.5 * dot1 - dot2; -#ifdef DEBUG +#ifdef DEBUG_IMSB std::cout << "e: " << e << '\n'; #endif energies[ iter ] = e; From b56f0ee9725ca1ed3a8b30a04195903ea12be1f8 Mon Sep 17 00:00:00 2001 From: Denis Jelovina Date: Wed, 16 Jul 2025 16:14:50 +0200 Subject: [PATCH 07/18] Suppress unused parameter warnings in Ising machine SB implementation --- include/graphblas/algorithms/ising_machine_sb.hpp | 1 + tests/smoke/ising_machine_sb.cpp | 2 ++ 2 files changed, 3 insertions(+) diff --git a/include/graphblas/algorithms/ising_machine_sb.hpp b/include/graphblas/algorithms/ising_machine_sb.hpp index 16b5c1707..1858ba617 100644 --- a/include/graphblas/algorithms/ising_machine_sb.hpp +++ b/include/graphblas/algorithms/ising_machine_sb.hpp @@ -117,6 +117,7 @@ namespace grb { const Minus & minus = Minus(), const Divide & divide = Divide(), const std::function< IOType( IOType ) > & sqrtX = std_sqrt< IOType, IOType > ) { + (void)minus; // suppress unused parameter warning constexpr const Descriptor descr_dense = descr | descriptors::dense; diff --git a/tests/smoke/ising_machine_sb.cpp b/tests/smoke/ising_machine_sb.cpp index e32d29aff..579b09db6 100644 --- a/tests/smoke/ising_machine_sb.cpp +++ b/tests/smoke/ising_machine_sb.cpp @@ -78,6 +78,8 @@ using JType = int; int main( int argc, char ** argv ) { + (void)argc; + (void)argv; //suppress unused parameter warnings /* --- Initialise ALP/GraphBLAS --- */ if( grb::init() != grb::SUCCESS ) { From fb7196f64cd3965b63c76de98ef8cc553c66b8e0 Mon Sep 17 00:00:00 2001 From: Denis Jelovina Date: Wed, 16 Jul 2025 17:00:09 +0200 Subject: [PATCH 08/18] Fix bug, use size() for vectors --- include/graphblas/algorithms/ising_machine_sb.hpp | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/include/graphblas/algorithms/ising_machine_sb.hpp b/include/graphblas/algorithms/ising_machine_sb.hpp index 1858ba617..e23200952 100644 --- a/include/graphblas/algorithms/ising_machine_sb.hpp +++ b/include/graphblas/algorithms/ising_machine_sb.hpp @@ -124,12 +124,9 @@ namespace grb { const std::size_t N = grb::nrows(J); assert( grb::ncols(J) == N ); - assert( grb::nrows(h) == N ); - assert( grb::ncols(h) == 1 ); - assert( grb::nrows(x_comp) == N ); - assert( grb::ncols(x_comp) == 1 ); - assert( grb::nrows(y_comp) == N ); - assert( grb::ncols(y_comp) == 1 ); + assert( grb::size(h) == N ); + assert( grb::size(x_comp) == N ); + assert( grb::size(y_comp) == N ); // TODO: check that J is symmetric once properly implemented //assert( grb::is_symmetric(J) ); From e29da4a0d68b96f3afef7381ef466f512c76260a Mon Sep 17 00:00:00 2001 From: Denis Jelovina Date: Thu, 17 Jul 2025 16:36:53 +0200 Subject: [PATCH 09/18] fix bug in foldl and eWiseLambda calls trigered in nonlocking --- .../graphblas/algorithms/ising_machine_sb.hpp | 29 +++++++------------ 1 file changed, 11 insertions(+), 18 deletions(-) diff --git a/include/graphblas/algorithms/ising_machine_sb.hpp b/include/graphblas/algorithms/ising_machine_sb.hpp index e23200952..b793c284d 100644 --- a/include/graphblas/algorithms/ising_machine_sb.hpp +++ b/include/graphblas/algorithms/ising_machine_sb.hpp @@ -151,7 +151,7 @@ namespace grb { grb::RC rc = grb::SUCCESS; /* ---- pre-compute ---- */ IOType sumJ2 = ring.template getZero< IOType >(); - rc = rc ? rc : grb::eWiseLambda( [&J, &sumJ2, &ring]( const size_t i, const size_t j, IsingHType& v ) { + rc = rc ? rc : grb::eWiseLambda( [&sumJ2, &ring]( const size_t i, const size_t j, IsingHType& v ) { (void) i; (void) j; IsingHType v2; @@ -168,10 +168,10 @@ namespace grb { #endif // rewrite this to use graphblas language - rc = rc ? rc : grb::foldl( sumJ2, static_cast( N - 1 ), divide ); + rc = rc ? rc : grb::foldl< descr_dense >( sumJ2, static_cast( N - 1 ), divide ); IOType xi = 0.5; sumJ2 = sqrtX( sumJ2 ); - rc = rc ? rc : grb::foldl( xi, sumJ2, divide ); + rc = rc ? rc : grb::foldl< descr_dense >( xi, sumJ2, divide ); #ifdef DEBUG_IMSB // for debugging purposes, print xi std::cout << "xi: " << xi << '\n'; @@ -243,8 +243,8 @@ namespace grb { rc = rc ? rc : grb::eWiseLambda< descr_dense >( [&mask, &x_comp]( const size_t i ) { (void) i; // rewrite this to use graphblas language - mask[i] = std::abs(x_comp[i]) > 1; - }, mask + mask[i] = std::abs(x_comp[i]) <= 1; + }, mask, x_comp ); assert( rc == grb::SUCCESS ); #ifdef DEBUG_IMSB @@ -252,19 +252,14 @@ namespace grb { #endif // y_comp[ mask ] = 0 - rc = rc ? rc : grb::eWiseLambda< descr_dense >( [&mask, &y_comp]( const size_t i ) { - (void) i; - // rewrite this to use graphblas language - if(mask[i]) { - y_comp[i] = 0; - } - }, y_comp - ); + //rc = rc ? rc : grb::set< (descr_dense/* |descriptors::structural */) >( y_comp, mask, ring.template getZero< IOType >() ); + rc = rc ? rc : grb::foldl< descr_dense >( y_comp, mask, ring.getMultiplicativeOperator() ); + assert( rc == grb::SUCCESS ); /* x_comp = clip( x_comp ) */ - foldl( x_comp, static_cast(-1), grb::operators::max < IOType >() ); - foldl( x_comp, static_cast(1), grb::operators::min < IOType >() ); + foldl< descr_dense >( x_comp, static_cast(-1), grb::operators::max < IOType >() ); + foldl< descr_dense >( x_comp, static_cast(1), grb::operators::min < IOType >() ); // alternatively, we could use eWiseLambda: // rc = rc ? rc : grb::eWiseLambda< descr_dense >( // [&x_comp]( const size_t i ) { @@ -288,14 +283,12 @@ namespace grb { (void) i; sol[i] = sign(x_comp[i]); }, - sol + sol, x_comp ); assert( rc == grb::SUCCESS ); #ifdef DEBUG_IMSB vector_print( sol, "sol" ); #endif - - // temp ← J * sol rc = rc ? rc : grb::set( temp_int, ring.template getZero< IsingHType >() ); rc = rc ? rc : grb::mxv< descr_dense >( temp_int, J, sol, ring ); From e1a2d24cd096518cf681225fc9351a719b3e77ae Mon Sep 17 00:00:00 2001 From: Denis Jelovina Date: Fri, 18 Jul 2025 17:56:18 +0200 Subject: [PATCH 10/18] fix mask use --- .../graphblas/algorithms/ising_machine_sb.hpp | 20 ++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/include/graphblas/algorithms/ising_machine_sb.hpp b/include/graphblas/algorithms/ising_machine_sb.hpp index b793c284d..fc47a70bd 100644 --- a/include/graphblas/algorithms/ising_machine_sb.hpp +++ b/include/graphblas/algorithms/ising_machine_sb.hpp @@ -228,7 +228,7 @@ namespace grb { rc = rc ? rc : grb::eWiseMul< descr_dense >( y_comp, dt, temp, ring ); assert( rc == grb::SUCCESS ); #ifdef DEBUG_IMSB - vector_print( y_comp, "y_comp" ); + vector_print( y_comp, "y_comp (a)" ); #endif /* x_comp += dt * y_comp */ @@ -243,23 +243,25 @@ namespace grb { rc = rc ? rc : grb::eWiseLambda< descr_dense >( [&mask, &x_comp]( const size_t i ) { (void) i; // rewrite this to use graphblas language - mask[i] = std::abs(x_comp[i]) <= 1; + mask[i] = std::abs(x_comp[i]) > 1; }, mask, x_comp ); assert( rc == grb::SUCCESS ); #ifdef DEBUG_IMSB vector_print( mask, "mask" ); #endif - - // y_comp[ mask ] = 0 - //rc = rc ? rc : grb::set< (descr_dense/* |descriptors::structural */) >( y_comp, mask, ring.template getZero< IOType >() ); - rc = rc ? rc : grb::foldl< descr_dense >( y_comp, mask, ring.getMultiplicativeOperator() ); - + rc = rc ? rc : grb::foldl< descr_dense >( + y_comp, mask, ring.template getZero< IOType >(), + grb::operators::right_assign() + ); +#ifdef DEBUG_IMSB + vector_print( y_comp, "y_comp (b)" ); +#endif assert( rc == grb::SUCCESS ); /* x_comp = clip( x_comp ) */ - foldl< descr_dense >( x_comp, static_cast(-1), grb::operators::max < IOType >() ); - foldl< descr_dense >( x_comp, static_cast(1), grb::operators::min < IOType >() ); + rc = rc ? rc : foldl< descr_dense >( x_comp, static_cast(-1), grb::operators::max < IOType >() ); + rc = rc ? rc : foldl< descr_dense >( x_comp, static_cast(1), grb::operators::min < IOType >() ); // alternatively, we could use eWiseLambda: // rc = rc ? rc : grb::eWiseLambda< descr_dense >( // [&x_comp]( const size_t i ) { From b8d69ce755c8d185bc2c73cad7dc87ce9e7a498b Mon Sep 17 00:00:00 2001 From: Denis Jelovina Date: Mon, 28 Jul 2025 17:21:20 +0200 Subject: [PATCH 11/18] Fix right_assign type --- include/graphblas/algorithms/ising_machine_sb.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/graphblas/algorithms/ising_machine_sb.hpp b/include/graphblas/algorithms/ising_machine_sb.hpp index fc47a70bd..641b24e6c 100644 --- a/include/graphblas/algorithms/ising_machine_sb.hpp +++ b/include/graphblas/algorithms/ising_machine_sb.hpp @@ -252,7 +252,7 @@ namespace grb { #endif rc = rc ? rc : grb::foldl< descr_dense >( y_comp, mask, ring.template getZero< IOType >(), - grb::operators::right_assign() + grb::operators::right_assign() ); #ifdef DEBUG_IMSB vector_print( y_comp, "y_comp (b)" ); From 7d31d3d7f502c76da55835af9eb9f5de7e23492c Mon Sep 17 00:00:00 2001 From: Denis Jelovina Date: Tue, 19 Aug 2025 13:48:26 +0200 Subject: [PATCH 12/18] Use tem J2 matrix to calclate sqared marrix values --- examples/CMakeLists.txt | 11 +++ .../graphblas/algorithms/ising_machine_sb.hpp | 75 ++++++++++--------- tests/smoke/ising_machine_sb.cpp | 12 ++- 3 files changed, 61 insertions(+), 37 deletions(-) diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 9a6affa1a..4835b4423 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -31,3 +31,14 @@ if( WITH_OMP_BACKEND ) add_dependencies( examples sp_reference_omp ) endif() +if( WITH_REFERENCE_BACKEND ) + add_executable( bbs_reference bbs.cpp ) + target_link_libraries( bbs_reference backend_reference common_flags ) + add_dependencies( examples bbs_reference ) +endif() + +if( WITH_REFERENCE_BACKEND ) + add_executable( bsb_Nsize_reference bsb_Nsize.cpp ) + target_link_libraries( bsb_Nsize_reference backend_reference common_flags ) + add_dependencies( examples bsb_Nsize_reference ) +endif() \ No newline at end of file diff --git a/include/graphblas/algorithms/ising_machine_sb.hpp b/include/graphblas/algorithms/ising_machine_sb.hpp index 641b24e6c..80db39530 100644 --- a/include/graphblas/algorithms/ising_machine_sb.hpp +++ b/include/graphblas/algorithms/ising_machine_sb.hpp @@ -91,10 +91,16 @@ namespace grb { grb::operators::add< IOType >, grb::operators::mul< IOType >, grb::identities::zero, - grb::identities::one + grb::identities::one >, class Minus = operators::subtract< IOType >, - class Divide = operators::divide< IOType > + class Divide = operators::divide< IOType >, + class RingIType = Semiring< + grb::operators::add< IsingHType >, + grb::operators::mul< IsingHType >, + grb::identities::zero, + grb::identities::one + > > grb::RC bSB( std::vector< IOType > & energies, // output length num_iters grb::Vector< IOType, backend > & x_comp, // in/out, size N @@ -106,7 +112,8 @@ namespace grb { const IOType p_end, const std::size_t num_iters, const IOType dt, - // workspace vectors + // workspace + grb::Matrix< IsingHType, backend, RSI, RSI, NZI > & J2, grb::Vector< IOType, backend > & Jx, grb::Vector< IOType, backend > & temp, grb::Vector< IsingHType, backend > & temp_int, @@ -116,6 +123,9 @@ namespace grb { const Ring & ring = Ring(), const Minus & minus = Minus(), const Divide & divide = Divide(), + const IOType zero = 0, + const RingIType & ringIType = RingIType(), + const IsingHType zero_itype = 0, const std::function< IOType( IOType ) > & sqrtX = std_sqrt< IOType, IOType > ) { (void)minus; // suppress unused parameter warning @@ -127,10 +137,14 @@ namespace grb { assert( grb::size(h) == N ); assert( grb::size(x_comp) == N ); assert( grb::size(y_comp) == N ); + + assert( grb::capacity(J) == grb::capacity(J2) ); + assert( grb::ncols(J) == grb::ncols(J2) ); + assert( grb::nrows(J) == grb::nrows(J2) ); // TODO: check that J is symmetric once properly implemented //assert( grb::is_symmetric(J) ); - grb::set( sol, ring.template getZero< IsingHType >() ); + grb::set( sol, zero_itype ); grb::set( mask, ring.template getZero< bool >() ); // print pinned vector x_comp @@ -149,15 +163,19 @@ namespace grb { } grb::RC rc = grb::SUCCESS; + /* ---- pre-compute ---- */ - IOType sumJ2 = ring.template getZero< IOType >(); - rc = rc ? rc : grb::eWiseLambda( [&sumJ2, &ring]( const size_t i, const size_t j, IsingHType& v ) { + rc = rc ? rc : grb::set( J2, J ); + assert( rc == grb::SUCCESS ); + rc = rc ? rc : grb::eWiseLambda( [&ring]( const size_t i, const size_t j, IsingHType& v ) { (void) i; (void) j; - IsingHType v2; - apply( v2, v, v, ring.getMultiplicativeOperator() ); - foldl( sumJ2, v2, ring.getAdditiveOperator() ); - }, J ); + apply( v, v, v, ring.getMultiplicativeOperator() ); + }, J2 ); + assert( rc == grb::SUCCESS ); + + IsingHType sumJ2 = zero_itype; + rc = rc ? rc : grb::foldl( sumJ2, J2, ringIType.getAdditiveMonoid() ); if( rc != grb::SUCCESS ) { std::cerr << "Error in eWiseLambda for sumJ2: " << rc << '\n'; return rc; @@ -167,11 +185,11 @@ namespace grb { std::cout << "sumJ2: " << sumJ2 << '\n'; #endif - // rewrite this to use graphblas language rc = rc ? rc : grb::foldl< descr_dense >( sumJ2, static_cast( N - 1 ), divide ); IOType xi = 0.5; - sumJ2 = sqrtX( sumJ2 ); - rc = rc ? rc : grb::foldl< descr_dense >( xi, sumJ2, divide ); + IOType sqrt_sumJ2 = zero_itype; + sqrt_sumJ2 = sqrtX( static_cast( sumJ2 ) ); + rc = rc ? rc : grb::foldl< descr_dense >( xi, sqrt_sumJ2, divide ); #ifdef DEBUG_IMSB // for debugging purposes, print xi std::cout << "xi: " << xi << '\n'; @@ -187,17 +205,17 @@ namespace grb { /* y_comp += ((-1+ps)*x_comp + xi*(Jx + h)) * dt */ - // Jx ← J * x_comp - grb::set( Jx, ring.template getZero< IOType >() ); + // Jx <- J * x_comp + grb::set( Jx, zero ); rc = rc ? rc : grb::mxv< descr_dense >( Jx, J, x_comp, ring ); assert( rc == grb::SUCCESS ); #ifdef DEBUG_IMSB vector_print( Jx, "Jx" ); #endif - // temp ← Jx + h - grb::set( temp, ring.template getZero< IOType >() ); - rc = rc ? rc : grb::eWiseApply< descr_dense >( + // temp <- Jx + h + grb::set( temp, zero ); + rc = rc ? rc : grb::eWiseApply( temp, Jx, h, ring.getAdditiveMonoid() ); assert( rc == grb::SUCCESS ); @@ -205,7 +223,7 @@ namespace grb { vector_print( temp, "temp = Jx + h" ); #endif - // temp ← xi * temp + // temp <- xi * temp rc = rc ? rc : grb::foldl< descr_dense >( temp, xi, ring.getMultiplicativeMonoid() ); @@ -213,9 +231,7 @@ namespace grb { #ifdef DEBUG_IMSB vector_print( temp, "xi * temp" ); #endif - - - // temp ← temp + (-1+ps) * x_comp + // temp <- temp + (-1+ps) * x_comp const IOType scale = -1.0 + ps; rc = rc ? rc : grb::eWiseMul< descr_dense >( temp, scale, x_comp, ring ); assert( rc == grb::SUCCESS ); @@ -251,7 +267,7 @@ namespace grb { vector_print( mask, "mask" ); #endif rc = rc ? rc : grb::foldl< descr_dense >( - y_comp, mask, ring.template getZero< IOType >(), + y_comp, mask, zero, grb::operators::right_assign() ); #ifdef DEBUG_IMSB @@ -262,15 +278,6 @@ namespace grb { /* x_comp = clip( x_comp ) */ rc = rc ? rc : foldl< descr_dense >( x_comp, static_cast(-1), grb::operators::max < IOType >() ); rc = rc ? rc : foldl< descr_dense >( x_comp, static_cast(1), grb::operators::min < IOType >() ); - // alternatively, we could use eWiseLambda: - // rc = rc ? rc : grb::eWiseLambda< descr_dense >( - // [&x_comp]( const size_t i ) { - // (void) i; - // foldl( x_comp[i], static_cast(-1), grb::operators::max < IOType >() ); - // foldl( x_comp[i], static_cast(1), grb::operators::min < IOType >() ); - // //x_comp[i] = std::min( 1.0, std::max( -1.0, x_comp[i] ) ); - // }, x_comp - // ); assert( rc == grb::SUCCESS ); #ifdef DEBUG_IMSB std::cout << "i = " << iter << "\n "; @@ -291,8 +298,8 @@ namespace grb { #ifdef DEBUG_IMSB vector_print( sol, "sol" ); #endif - // temp ← J * sol - rc = rc ? rc : grb::set( temp_int, ring.template getZero< IsingHType >() ); + // temp <- J * sol + rc = rc ? rc : grb::set( temp_int, zero_itype ); rc = rc ? rc : grb::mxv< descr_dense >( temp_int, J, sol, ring ); assert( rc == grb::SUCCESS ); #ifdef DEBUG_IMSB diff --git a/tests/smoke/ising_machine_sb.cpp b/tests/smoke/ising_machine_sb.cpp index 579b09db6..c991e5ead 100644 --- a/tests/smoke/ising_machine_sb.cpp +++ b/tests/smoke/ising_machine_sb.cpp @@ -119,13 +119,19 @@ int main( int argc, char ** argv ) { grb::Vector< IOType > Jx( N ); grb::Vector< IOType > temp( N ); grb::Vector< JType > temp_int( N ); - grb::Vector< bool > mask( N ); + grb::Vector< bool > mask( N ); + grb::Matrix< JType > J2( N, N ); + rc = rc ? rc : grb::resize( J2, grb::nnz(J) ); + if(rc != grb::SUCCESS) { + std::cerr << "Matrix resize failed for J2\n"; + return grb::RC::PANIC; + } // TODO: make sol int type - grb::Vector< JType > sol( N ); + grb::Vector< JType > sol( N ); rc = rc ? rc : bSB( energies, x0, y0, J, h, p0, p1, num_iters, dt, - Jx, temp, temp_int, mask, sol + J2, Jx, temp, temp_int, mask, sol ); From e28991acffda9eda45ba1a7e5cbd26ae5e2d4369 Mon Sep 17 00:00:00 2001 From: Denis Jelovina Date: Wed, 20 Aug 2025 14:19:35 +0200 Subject: [PATCH 13/18] Refactor test data structure and improve input/output handling in ising_machine_sb.cpp --- tests/smoke/ising_machine_sb.cpp | 338 +++++++++++++++++++++++++------ 1 file changed, 279 insertions(+), 59 deletions(-) diff --git a/tests/smoke/ising_machine_sb.cpp b/tests/smoke/ising_machine_sb.cpp index c991e5ead..fff50ef1e 100644 --- a/tests/smoke/ising_machine_sb.cpp +++ b/tests/smoke/ising_machine_sb.cpp @@ -34,70 +34,105 @@ #include +#include #include - using namespace grb; using namespace algorithms; +using IOType = double; +using JType = int; -// test data from python implementation -constexpr std::size_t N = 10; -constexpr std::size_t Nz = 54; -static const size_t i_arr[ Nz ] = { - 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, - 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 9, 9 -}; -static const size_t j_arr[ Nz ] = { - 0, 2, 3, 4, 5, 1, 4, 5, 6, 7, 9, 0, 2, 4, 6, 9, 0, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, - 3, 4, 0, 1, 3, 5, 6, 8, 1, 2, 3, 5, 6, 1, 3, 7, 8, 9, 3, 5, 7, 8, 1, 2, 3, 7, 9 -}; -static const int v_arr[ Nz ] = { - -1, 1, 1, -1, -1, -1, -1, 1, 1, 1, 1, 1, -1, -1, -1, -1, 1, 1, 1, 1, -1, -1, 1, 1, -1, -1, -1, - 1, 1, -1, 1, 1, -1, -1, -1, 1, -1, -1, -1, -1, 1, -1, -1, 1, -1, 1, -1, 1, 1, 1, -1, 1, -1, 1 + +namespace test_data { + // test data from python implementation + constexpr const std::size_t max_iters = 100; + + constexpr std::size_t N = 10; + constexpr std::size_t Nz = 54; + static const size_t i_arr[ Nz ] = { + 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, + 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 9, 9 + }; + static const size_t j_arr[ Nz ] = { + 0, 2, 3, 4, 5, 1, 4, 5, 6, 7, 9, 0, 2, 4, 6, 9, 0, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, + 3, 4, 0, 1, 3, 5, 6, 8, 1, 2, 3, 5, 6, 1, 3, 7, 8, 9, 3, 5, 7, 8, 1, 2, 3, 7, 9 + }; + static const int v_arr[ Nz ] = { + -1, 1, 1, -1, -1, -1, -1, 1, 1, 1, 1, 1, -1, -1, -1, -1, 1, 1, 1, 1, -1, -1, 1, 1, -1, -1, -1, + 1, 1, -1, 1, 1, -1, -1, -1, 1, -1, -1, -1, -1, 1, -1, -1, 1, -1, 1, -1, 1, 1, 1, -1, 1, -1, 1 + }; + + static const int h_arr[ N ] = { 1, -1, 1, -1, 1, 1, -1, 1, 1, 1 }; + static const double x_arr[ N ] = { -0.0996, -0.0315, 0.0572, 0.0630, 0.0087, -0.0143, -0.0170, -0.0411, 0.0433, -0.0298 }; + static const double y_arr[ N ] = { 0.0373, 0.0540, 0.0486, -0.0877, -0.0418, -0.0261, 0.0018, -0.0710, 0.0507, -0.0483 }; + + const std::size_t num_iters = 100; + static const double energies_ref[ num_iters ] = { + -3, -3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, + -1, -9, -11, -11, -11, -11, -11, -11, -11, -11, -9, -9, -11, + -9, -9, -9, -9, -9, -9, -9, -9, -9, -9, -9, -9, -9, + -9, -9, -11, -11, -11, -9, -9, -9, -9, -9, -9, -9, -9, + -9, -9, -9, -9, -9, -9, -9, -11, -11, -11, -11, -11, -11, + -11, -11, -11, -11, -11, -11, -11, -11, -11, -11, -11, -11, -11, + -11, -11, -11, -11, -11, -11, -11, -11, -11, -11, -11, -11, -11, + -11, -11, -11, -11, -11, -11, -11, -11, -11 + }; + + static const int sol_ref_data[ N ] = { + 1, -1, 1, 1, 1, -1, -1, -1, 1, 1 + }; + + const IOType p0 = 0.; + const IOType p1 = 1.1; + const IOType dt = 0.25; + +} // namespace test_data + +struct input { + // contains the command line arguments + bool direct; + size_t rep; }; -static const int h_arr[ N ] = { 1, -1, 1, -1, 1, 1, -1, 1, 1, 1 }; -static const double x_arr[ N ] = { -0.0996, -0.0315, 0.0572, 0.0630, 0.0087, -0.0143, -0.0170, -0.0411, 0.0433, -0.0298 }; -static const double y_arr[ N ] = { 0.0373, 0.0540, 0.0486, -0.0877, -0.0418, -0.0261, 0.0018, -0.0710, 0.0507, -0.0483 }; - -const std::size_t num_iters = 100; -static const double energies_ref[ num_iters ] = { - -3, -3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, - -1, -9, -11, -11, -11, -11, -11, -11, -11, -11, -9, -9, -11, - -9, -9, -9, -9, -9, -9, -9, -9, -9, -9, -9, -9, -9, - -9, -9, -11, -11, -11, -9, -9, -9, -9, -9, -9, -9, -9, - -9, -9, -9, -9, -9, -9, -9, -11, -11, -11, -11, -11, -11, - -11, -11, -11, -11, -11, -11, -11, -11, -11, -11, -11, -11, -11, - -11, -11, -11, -11, -11, -11, -11, -11, -11, -11, -11, -11, -11, - -11, -11, -11, -11, -11, -11, -11, -11, -11 +struct output { + int error_code = 0; + size_t rep; + size_t iterations; + grb::utils::TimerResults times; + std::unique_ptr< PinnedVector< JType > > pinnedSolutionVector; + std::unique_ptr< PinnedVector< JType > > pinnedRefSolutionVector; }; -using IOType = double; -using JType = int; +void grbProgram( + const struct input &data_in, + struct output &out +) { + using namespace test_data; -int main( int argc, char ** argv ) { - (void)argc; - (void)argv; //suppress unused parameter warnings + // get user process ID + const size_t s = spmd<>::pid(); + assert( s < spmd<>::nprocs() ); - /* --- Initialise ALP/GraphBLAS --- */ - if( grb::init() != grb::SUCCESS ) { - std::cerr << "ALP init failed\n"; - return -1; - } + grb::utils::Timer timer; + timer.reset(); + // TODO: IO goes here + + // I/O done + out.times.io = timer.time(); + timer.reset(); - /* --- Problem setup (toy) --- */ + /* --- Problem setup --- */ grb::Matrix J( N, N, Nz ); grb::Vector h( N ); grb::Vector x0( N ), y0( N ); // initialy // ... populate J with test (random) values grb::RC rc = grb::SUCCESS; rc = rc ? rc : buildMatrixUnique( J, &( i_arr[ 0 ] ), &( j_arr[ 0 ] ), &( v_arr[ 0 ] ), Nz, grb::SEQUENTIAL ); - if(rc != grb::SUCCESS) { std::cerr << "matrix build failed\n"; - return grb::RC::PANIC; + return; } // Fill h, x0, y0 with random values using buildVector @@ -106,17 +141,13 @@ int main( int argc, char ** argv ) { rc = rc ? rc : buildVector(y0, y_arr, y_arr + N, grb::SEQUENTIAL); if(rc != grb::SUCCESS) { std::cerr << "Vector build failed\n"; - return grb::RC::PANIC; + return; } - const IOType p0 = 0.; - const IOType p1 = 1.1; - const IOType dt = 0.25; - // energies is array of length num_iters, initialized to 0 std::vector< IOType > energies( num_iters, 0 ); - grb::Vector< IOType > Jx( N ); + grb::Vector< IOType > Jx( N ); grb::Vector< IOType > temp( N ); grb::Vector< JType > temp_int( N ); grb::Vector< bool > mask( N ); @@ -124,15 +155,110 @@ int main( int argc, char ** argv ) { rc = rc ? rc : grb::resize( J2, grb::nnz(J) ); if(rc != grb::SUCCESS) { std::cerr << "Matrix resize failed for J2\n"; - return grb::RC::PANIC; + return; } - // TODO: make sol int type grb::Vector< JType > sol( N ); + grb::Vector< JType > sol_ref( N ); + rc = rc ? rc : buildVector(sol_ref, sol_ref_data, sol_ref_data + N, grb::SEQUENTIAL); + + out.times.preamble = timer.time(); + + rc = rc ? rc : wait(); + out.times.preamble = timer.time(); + + // by default, copy input requested repetitions to output repititions performed + out.rep = data_in.rep; + // time a single call + if( out.rep == 0 ) { + timer.reset(); + rc = bSB( + energies, x0, y0, J, h, p0, p1, num_iters, dt, + J2, Jx, temp, temp_int, mask, sol, out.iterations + ); + + rc = rc ? rc : wait(); + double single_time = timer.time(); + if( !(rc == SUCCESS || rc == FAILED) ) { + std::cerr << "Failure: call to Ising Machine SB did not succeed (" + << toString( rc ) << ")." << std::endl; + out.error_code = 20; + } + if( rc == FAILED ) { + std::cout << "Warning: call to Ising Machine SB did not converge\n"; + } + if( rc == SUCCESS ) { + rc = collectives<>::reduce( single_time, 0, operators::max< double >() ); + } + if( rc != SUCCESS ) { + out.error_code = 25; + } + out.times.useful = single_time; + out.rep = static_cast< size_t >( 1000.0 / single_time ) + 1; + if( rc == SUCCESS || rc == FAILED ) { + if( s == 0 ) { + if( rc == FAILED ) { + std::cout << "Info: cold Ising Machine SB did not converge within "; + } else { + std::cout << "Info: cold Ising Machine SB completed within "; + } + std::cout << out.iterations << " iterations. " + << "Time taken was " << single_time << " ms. " + << "Deduced inner repetitions parameter of " << out.rep << " " + << "to take 1 second or more per inner benchmark.\n"; + } + } + } else { + // do benchmark + timer.reset(); + for( size_t i = 0; i < out.rep && rc == SUCCESS; ++i ) { + if( rc == SUCCESS ) { + rc = bSB( + energies, x0, y0, J, h, p0, p1, num_iters, dt, + J2, Jx, temp, temp_int, mask, sol, out.iterations + ); + } + if( Properties<>::isNonblockingExecution ) { + rc = rc ? rc : wait(); + } + } + const double time_taken = timer.time(); + out.times.useful = time_taken / static_cast< double >( out.rep ); + // print timing at root process + if( grb::spmd<>::pid() == 0 ) { + std::cout << "Time taken for " << out.rep << " " + << "Ising Machine SB calls (hot start): " << out.times.useful << ". " + << "Error code is " << grb::toString( rc ) << std::endl; + std::cout << "\tnumber of IM-SB iterations: " << out.iterations << "\n"; + std::cout << "\tmilliseconds per iteration: " + << ( out.times.useful / static_cast< double >( out.iterations ) ) + << "\n"; + } + sleep( 1 ); + } + + // start postamble + timer.reset(); + + // set error code + if( rc == FAILED ) { + out.error_code = 30; + } else if( rc != SUCCESS ) { + std::cerr << "Benchmark run returned error: " << toString( rc ) << "\n"; + out.error_code = 35; + return; + } + + out.pinnedRefSolutionVector = std::unique_ptr< PinnedVector< JType > >( + new PinnedVector< JType >( sol_ref, SEQUENTIAL ) ); + + // output + out.pinnedSolutionVector = std::unique_ptr< PinnedVector< JType > >( + new PinnedVector< JType >( sol, SEQUENTIAL ) ); + + // finish timing + const double time_taken = timer.time(); + out.times.postamble = time_taken; - rc = rc ? rc : bSB( - energies, x0, y0, J, h, p0, p1, num_iters, dt, - J2, Jx, temp, temp_int, mask, sol - ); if( rc != grb::SUCCESS ) { @@ -140,17 +266,111 @@ int main( int argc, char ** argv ) { } else { // print all energies for (std::size_t i = 0; i < num_iters; ++i) { +#ifdef DEBUG_IMSB std::cout << "Energy at iteration " << i << " = " << energies[i] << '\n'; +#endif if( energies[i] != energies_ref[i]) { - std::cerr << "Error: Energy at iteration " << i << " does not match reference value.\n"; - std::cerr << "Expected: " << energies_ref[i] << ", got: " << energies[i] << '\n'; - return -1; +#ifdef DEBUG_IMSB + std::cerr << "Error: Energy at iteration " << i << " does not match reference value.\n"; + std::cerr << "Expected: " << energies_ref[i] << ", got: " << energies[i] << '\n'; +#endif + out.error_code = 40; + return ; } } std::cout << "All energies match reference values.\n"; std::cout << "TEST OK\n"; } - grb::finalize(); - return 0; + // set error code + out.error_code = rc; + +} + + +int main( int argc, char ** argv ) { + // TODO: add argument parsing for input file, direct/indirect addressing, etc. + // for now, just print the executable name + (void) argc; // unused + (void) argv; // unused + std::cout << "Test executable: " << argv[ 0 ] << std::endl; + + // the input struct + struct input in; + + // get inner number of iterations + in.rep = grb::config::BENCHMARKING::inner(); + + // get outer number of iterations + size_t outer = grb::config::BENCHMARKING::outer(); + + // check for verification of the output + // bool verification = false; + + std::cout << "Executable called with parameters " + << "inner repititions = " << in.rep << ", " + << "outer reptitions = " << outer << ", " + << std::endl; + + // set standard exit code + grb::RC rc = SUCCESS; + + // the output struct + struct output out; + + // launch estimator (if requested) + if( in.rep == 0 ) { + grb::Launcher< AUTOMATIC > launcher; + rc = launcher.exec( &grbProgram, in, out, true ); + if( rc == SUCCESS ) { + in.rep = out.rep; + } + if( rc != SUCCESS ) { + std::cerr << "launcher.exec returns with non-SUCCESS error code " + << (int)rc << std::endl; + return 80; + } + } + + // launch benchmark + if( rc == SUCCESS ) { + grb::Benchmarker< AUTOMATIC > benchmarker; + rc = benchmarker.exec( &grbProgram, in, out, 1, outer, true ); + } + if( rc != SUCCESS ) { + std::cerr << "benchmarker.exec returns with non-SUCCESS error code " + << grb::toString( rc ) << std::endl; + return 90; + } else if( out.error_code == 0 ) { + std::cout << "Benchmark completed successfully and took " << out.iterations + << ".\n"; + } + + std::cout << "Error code is " << out.error_code << ".\n"; + + if( !(out.pinnedSolutionVector) ) { + std::cerr << "no output vector to inspect" << std::endl; + } else { + const PinnedVector< JType > &solution = *(out.pinnedSolutionVector); + const PinnedVector< JType > &solution_ref = *(out.pinnedRefSolutionVector); + std::cout << "Size of x is " << solution.size() << std::endl; + if( solution.size() > 0 ) { + print_vector( solution, 30, "SOLUTION" ); + // expected solution from sol_ref_data + print_vector( solution_ref, 30, "EXPECTED SOLUTION" ); + } else { + std::cerr << "ERROR: solution contains no values" << std::endl; + } + } + + + if( out.error_code != 0 ) { + std::cerr << std::flush; + std::cout << "Test FAILED\n"; + } + + std::cout << "Test OK\n"; + + // done + return out.error_code; } From ba722f49773a6f3748335181a82cad6a0d4b237a Mon Sep 17 00:00:00 2001 From: Denis Jelovina Date: Wed, 20 Aug 2025 15:53:01 +0200 Subject: [PATCH 14/18] Add iteration tracking and I/O handling in Ising machine tests --- .../graphblas/algorithms/ising_machine_sb.hpp | 5 +- tests/smoke/ising_machine_sb.cpp | 160 +++++++++++++++--- 2 files changed, 144 insertions(+), 21 deletions(-) diff --git a/include/graphblas/algorithms/ising_machine_sb.hpp b/include/graphblas/algorithms/ising_machine_sb.hpp index 80db39530..3be5b6c43 100644 --- a/include/graphblas/algorithms/ising_machine_sb.hpp +++ b/include/graphblas/algorithms/ising_machine_sb.hpp @@ -119,6 +119,7 @@ namespace grb { grb::Vector< IsingHType, backend > & temp_int, grb::Vector< bool, backend > & mask, grb::Vector< IsingHType, backend > & sol, + size_t & iterations, // default semiring, minus, divide const Ring & ring = Ring(), const Minus & minus = Minus(), @@ -201,7 +202,7 @@ namespace grb { // assert len of energies == N assert( energies.size() == num_iters ); - for ( std::size_t iter = 0; iter < num_iters; ++iter ) { + for ( iterations = 0; iterations < num_iters; ++iterations ) { /* y_comp += ((-1+ps)*x_comp + xi*(Jx + h)) * dt */ @@ -324,7 +325,7 @@ namespace grb { #ifdef DEBUG_IMSB std::cout << "e: " << e << '\n'; #endif - energies[ iter ] = e; + energies[ iterations ] = e; ps += dps; } diff --git a/tests/smoke/ising_machine_sb.cpp b/tests/smoke/ising_machine_sb.cpp index fff50ef1e..496ced988 100644 --- a/tests/smoke/ising_machine_sb.cpp +++ b/tests/smoke/ising_machine_sb.cpp @@ -43,25 +43,59 @@ using namespace algorithms; using IOType = double; using JType = int; +/** Parser type */ +typedef grb::utils::MatrixFileReader< + JType, + std::conditional< + (sizeof(grb::config::RowIndexType) > sizeof(grb::config::ColIndexType)), + grb::config::RowIndexType, + grb::config::ColIndexType + >::type +> Parser; + +/** Nonzero type */ +typedef internal::NonzeroStorage< + grb::config::RowIndexType, + grb::config::ColIndexType, + JType +> NonzeroT; + +/** In-memory storage type using tuple */ +typedef grb::utils::Singleton< + std::tuple< + size_t, // n (rows/columns) + size_t, // nz (nonzeros) + std::vector, // matrix data + std::vector, // h vector + std::vector, // x vector + std::vector // y vector + // Add more types as needed + > +> Storage; + +// Access using std::get +// auto& n = std::get<0>(Storage::getData()); +// auto& nz = std::get<1>(Storage::getData()); +// auto& matrix_data = std::get<2>(Storage::getData()); namespace test_data { // test data from python implementation constexpr const std::size_t max_iters = 100; constexpr std::size_t N = 10; - constexpr std::size_t Nz = 54; - static const size_t i_arr[ Nz ] = { - 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, - 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 9, 9 - }; - static const size_t j_arr[ Nz ] = { - 0, 2, 3, 4, 5, 1, 4, 5, 6, 7, 9, 0, 2, 4, 6, 9, 0, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, - 3, 4, 0, 1, 3, 5, 6, 8, 1, 2, 3, 5, 6, 1, 3, 7, 8, 9, 3, 5, 7, 8, 1, 2, 3, 7, 9 - }; - static const int v_arr[ Nz ] = { - -1, 1, 1, -1, -1, -1, -1, 1, 1, 1, 1, 1, -1, -1, -1, -1, 1, 1, 1, 1, -1, -1, 1, 1, -1, -1, -1, - 1, 1, -1, 1, 1, -1, -1, -1, 1, -1, -1, -1, -1, 1, -1, -1, 1, -1, 1, -1, 1, 1, 1, -1, 1, -1, 1 - }; + // constexpr std::size_t Nz = 54; + // static const size_t i_arr[ Nz ] = { + // 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, + // 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 9, 9 + // }; + // static const size_t j_arr[ Nz ] = { + // 0, 2, 3, 4, 5, 1, 4, 5, 6, 7, 9, 0, 2, 4, 6, 9, 0, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, + // 3, 4, 0, 1, 3, 5, 6, 8, 1, 2, 3, 5, 6, 1, 3, 7, 8, 9, 3, 5, 7, 8, 1, 2, 3, 7, 9 + // }; + // static const int v_arr[ Nz ] = { + // -1, 1, 1, -1, -1, -1, -1, 1, 1, 1, 1, 1, -1, -1, -1, -1, 1, 1, 1, 1, -1, -1, 1, 1, -1, -1, -1, + // 1, 1, -1, 1, 1, -1, -1, -1, 1, -1, -1, -1, -1, 1, -1, -1, 1, -1, 1, -1, 1, 1, 1, -1, 1, -1, 1 + // }; static const int h_arr[ N ] = { 1, -1, 1, -1, 1, 1, -1, 1, 1, 1 }; static const double x_arr[ N ] = { -0.0996, -0.0315, 0.0572, 0.0630, 0.0087, -0.0143, -0.0170, -0.0411, 0.0433, -0.0298 }; @@ -91,6 +125,7 @@ namespace test_data { struct input { // contains the command line arguments + std::string filename; bool direct; size_t rep; }; @@ -105,6 +140,41 @@ struct output { }; +void ioProgram( const struct input &data_in, bool &success ) { + success = false; + // Parse and store matrix in singleton class + auto &data = std::get<2>(Storage::getData()); + try { + Parser parser( data_in.filename, data_in.direct ); + assert( parser.m() == parser.n() ); + std::get<0>(Storage::getData()) = parser.n(); + try { + std::get<1>(Storage::getData()) = parser.nz(); + } catch( ... ) { + std::get<1>(Storage::getData()) = parser.entries(); + } + /* Once internal issue #342 is resolved this can be re-enabled + for( + auto it = parser.begin( PARALLEL ); + it != parser.end( PARALLEL ); + ++it + ) { + data.push_back( *it ); + }*/ + for( + auto it = parser.begin( SEQUENTIAL ); + it != parser.end( SEQUENTIAL ); + ++it + ) { + data.push_back( NonzeroT( *it ) ); + } + } catch( std::exception &e ) { + std::cerr << "I/O program failed: " << e.what() << "\n"; + return; + } + success = true; +} + void grbProgram( const struct input &data_in, struct output &out @@ -124,16 +194,51 @@ void grbProgram( timer.reset(); /* --- Problem setup --- */ - grb::Matrix J( N, N, Nz ); grb::Vector h( N ); grb::Vector x0( N ), y0( N ); // initialy // ... populate J with test (random) values grb::RC rc = grb::SUCCESS; - rc = rc ? rc : buildMatrixUnique( J, &( i_arr[ 0 ] ), &( j_arr[ 0 ] ), &( v_arr[ 0 ] ), Nz, grb::SEQUENTIAL ); - if(rc != grb::SUCCESS) { - std::cerr << "matrix build failed\n"; - return; - } + // rc = rc ? rc : buildMatrixUnique( J, &( i_arr[ 0 ] ), &( j_arr[ 0 ] ), &( v_arr[ 0 ] ), Nz, grb::SEQUENTIAL ); + // if(rc != grb::SUCCESS) { + // std::cerr << "matrix build failed\n"; + // return; + // } + + + // load into GraphBLAS + const size_t n = std::get<0>(Storage::getData()); + grb::Matrix J( n, n ); + { + const auto &data = std::get<2>(Storage::getData()); + RC io_rc = buildMatrixUnique( + J, + utils::makeNonzeroIterator< + grb::config::RowIndexType, grb::config::ColIndexType, JType + >( data.cbegin() ), + utils::makeNonzeroIterator< + grb::config::RowIndexType, grb::config::ColIndexType, JType + >( data.cend() ), + SEQUENTIAL + ); + /* Once internal issue #342 is resolved this can be re-enabled + RC io_rc = buildMatrixUnique( + J, + utils::makeNonzeroIterator< + grb::config::RowIndexType, grb::config::ColIndexType, JType + >( data.cbegin() ), + utils::makeNonzeroIterator< + grb::config::RowIndexType, grb::config::ColIndexType, JType + >( data.cend() ), + PARALLEL + );*/ + io_rc = io_rc ? io_rc : wait(); + if( io_rc != SUCCESS ) { + std::cerr << "Failure: call to buildMatrixUnique did not succeed " + << "(" << toString( io_rc ) << ")." << std::endl; + out.error_code = 5; + return; + } + } // Fill h, x0, y0 with random values using buildVector rc = rc ? rc : buildVector(h, h_arr, h_arr + N, grb::SEQUENTIAL); @@ -297,6 +402,7 @@ int main( int argc, char ** argv ) { // the input struct struct input in; + in.filename = "/home/d/Scratch/SA/ising_machine.mtx"; // get inner number of iterations in.rep = grb::config::BENCHMARKING::inner(); @@ -315,6 +421,22 @@ int main( int argc, char ** argv ) { // set standard exit code grb::RC rc = SUCCESS; + // launch I/O + { + bool success; + grb::Launcher< AUTOMATIC > launcher; + rc = launcher.exec( &ioProgram, in, success, true ); + if( rc != SUCCESS ) { + std::cerr << "launcher.exec(I/O) returns with non-SUCCESS error code \"" + << grb::toString( rc ) << "\"\n"; + return 73; + } + if( !success ) { + std::cerr << "I/O program caught an exception\n"; + return 77; + } + } + // the output struct struct output out; From 99280d06f5b707688d6ec0449c9b675c2d66096c Mon Sep 17 00:00:00 2001 From: Denis Jelovina Date: Thu, 21 Aug 2025 16:46:29 +0200 Subject: [PATCH 15/18] Refactor test data structure and enhance I/O handling in ising_machine_sb.cpp --- tests/smoke/ising_machine_sb.cpp | 283 ++++++++++++++++++++++++------- 1 file changed, 222 insertions(+), 61 deletions(-) diff --git a/tests/smoke/ising_machine_sb.cpp b/tests/smoke/ising_machine_sb.cpp index 496ced988..f4fa2f105 100644 --- a/tests/smoke/ising_machine_sb.cpp +++ b/tests/smoke/ising_machine_sb.cpp @@ -72,7 +72,6 @@ typedef grb::utils::Singleton< // Add more types as needed > > Storage; - // Access using std::get // auto& n = std::get<0>(Storage::getData()); // auto& nz = std::get<1>(Storage::getData()); @@ -80,27 +79,49 @@ typedef grb::utils::Singleton< namespace test_data { // test data from python implementation - constexpr const std::size_t max_iters = 100; - - constexpr std::size_t N = 10; - // constexpr std::size_t Nz = 54; - // static const size_t i_arr[ Nz ] = { - // 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, - // 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 9, 9 - // }; - // static const size_t j_arr[ Nz ] = { - // 0, 2, 3, 4, 5, 1, 4, 5, 6, 7, 9, 0, 2, 4, 6, 9, 0, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, - // 3, 4, 0, 1, 3, 5, 6, 8, 1, 2, 3, 5, 6, 1, 3, 7, 8, 9, 3, 5, 7, 8, 1, 2, 3, 7, 9 - // }; - // static const int v_arr[ Nz ] = { - // -1, 1, 1, -1, -1, -1, -1, 1, 1, 1, 1, 1, -1, -1, -1, -1, 1, 1, 1, 1, -1, -1, 1, 1, -1, -1, -1, - // 1, 1, -1, 1, 1, -1, -1, -1, 1, -1, -1, -1, -1, 1, -1, -1, 1, -1, 1, -1, 1, 1, 1, -1, 1, -1, 1 - // }; - - static const int h_arr[ N ] = { 1, -1, 1, -1, 1, 1, -1, 1, 1, 1 }; - static const double x_arr[ N ] = { -0.0996, -0.0315, 0.0572, 0.0630, 0.0087, -0.0143, -0.0170, -0.0411, 0.0433, -0.0298 }; - static const double y_arr[ N ] = { 0.0373, 0.0540, 0.0486, -0.0877, -0.0418, -0.0261, 0.0018, -0.0710, 0.0507, -0.0483 }; + constexpr const std::size_t n = 10; + + const std::vector h_array_data = { 1, -1, 1, -1, 1, 1, -1, 1, 1, 1 }; + const std::vector x_array_data = { + -0.0996, -0.0315, 0.0572, 0.0630, 0.0087, + -0.0143, -0.0170, -0.0411, 0.0433, -0.0298 + }; + const std::vector y_array_data = { + 0.0373, 0.0540, 0.0486, -0.0877, -0.0418, + -0.0261, 0.0018, -0.0710, 0.0507, -0.0483 + }; + const std::vector sol_ref_data = { 1, -1, 1, 1, 1, -1, -1, -1, 1, 1 }; + const std::vector, JType > > j_matrix_data = { + {{1, 1}, -1}, {{2, 2}, -1}, {{3, 1}, 1}, {{3, 3}, -1}, + {{4, 1}, 1}, {{4, 4}, 1}, {{5, 1}, -1}, {{5, 2}, -1}, + {{5, 3}, -1}, {{5, 4}, 1}, {{5, 5}, 1}, + {{6, 1}, -1 }, {{6, 2}, 1}, {{6, 4}, 1}, + {{6, 6}, -1}, {{7, 2}, 1}, {{7, 3}, -1}, + {{7, 4}, -1}, {{7, 6}, -1}, + {{7, 7}, -1}, {{8, 2}, 1}, + {{8, 4}, -1}, {{8, 8}, -1}, + {{9, 4}, 1}, {{9, 6}, -1}, + {{9, 8}, 1}, {{9, 9}, 1}, + {{10, 2}, 1}, {{10, 3}, -1}, + {{10, 4}, 1}, {{10, 8}, -1}, + {{10, 10}, 1} + , + // since matrix is symmetric, we can add the symmetric entries + {{1, 3}, 1}, + {{1, 4}, 1}, {{1, 5}, -1}, {{2, 5}, -1}, + {{3, 5}, -1}, {{4, 5}, 1}, + {{1, 6}, -1 }, {{2, 6}, 1}, {{4, 6}, 1}, + {{2, 7}, 1}, {{3, 7}, -1}, + {{4, 7}, -1}, {{6, 7}, -1}, + {{2, 8}, 1}, + {{4, 8}, -1}, + {{4, 9}, 1}, {{6, 9}, -1}, + {{8, 9}, 1}, + {{2, 10}, 1}, {{3, 10}, -1}, + {{4, 10}, 1}, {{8, 10}, -1} + }; + constexpr const std::size_t max_iters = 100; const std::size_t num_iters = 100; static const double energies_ref[ num_iters ] = { -3, -3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, @@ -113,9 +134,6 @@ namespace test_data { -11, -11, -11, -11, -11, -11, -11, -11, -11 }; - static const int sol_ref_data[ N ] = { - 1, -1, 1, 1, 1, -1, -1, -1, 1, 1 - }; const IOType p0 = 0.; const IOType p1 = 1.1; @@ -124,12 +142,18 @@ namespace test_data { } // namespace test_data struct input { - // contains the command line arguments - std::string filename; + // contains the command line arguments related data + static bool use_default_data; + std::string filename_Jmatrix; + std::string filename_h; + std::string filename_x; + std::string filename_y; bool direct; size_t rep; }; +bool input::use_default_data = false; + struct output { int error_code = 0; size_t rep; @@ -139,13 +163,11 @@ struct output { std::unique_ptr< PinnedVector< JType > > pinnedRefSolutionVector; }; - -void ioProgram( const struct input &data_in, bool &success ) { - success = false; - // Parse and store matrix in singleton class - auto &data = std::get<2>(Storage::getData()); +template< typename Dtype > +void read_matrix_data(const std::string &filename, std::vector &data, bool direct) { + // Implementation for reading matrix data from file try { - Parser parser( data_in.filename, data_in.direct ); + Parser parser( filename, direct ); assert( parser.m() == parser.n() ); std::get<0>(Storage::getData()) = parser.n(); try { @@ -166,12 +188,114 @@ void ioProgram( const struct input &data_in, bool &success ) { it != parser.end( SEQUENTIAL ); ++it ) { - data.push_back( NonzeroT( *it ) ); + data.push_back( Dtype( *it ) ); + // print last data element + std::cout << "read_matrix_data_from_file: " << data.back().first.first << ", " + << data.back().first.second << ", " << data.back().second << "\n"; } } catch( std::exception &e ) { std::cerr << "I/O program failed: " << e.what() << "\n"; return; } +} + +template< typename NonzeroT, typename IType, typename VType > +void read_matrix_data_from_array( + const std::vector, VType > > &array, + std::vector &data +) { + // Implementation for reading matrix data from array + try { + for (const auto &entry : array) { + //std::cout << "read_matrix_data_from_array: " << entry.first.first << ", " << entry.first.second << ", " << entry.second << "\n"; + // since + // data.push_back( Dtype( *it ) ); + // print last data element + // std::cout << "read_matrix_data: " << it->first.first << ", " << it->first.second << ", " << it->second << "\n"; + // we want the same here + data.emplace_back( + NonzeroT( entry.first.first-1, entry.first.second-1, entry.second ) + ); + // print last data element from std::vector data + std::cout << "read_matrix_data_from_array: " << data.back().first.first << ", " + << data.back().first.second << ", " << data.back().second << "\n"; + } + std::get<0>(Storage::getData()) = test_data::n; + std::get<1>(Storage::getData()) = data.size(); + } catch (const std::exception &e) { + std::cerr << "Failed to read matrix data from array: " << e.what() << "\n"; + return; + } +} + +template< typename Dtype > +void read_vector_data(const std::string &filename, std::vector &data) { + // Implementation for reading vector data from file + try { + std::ifstream file( filename ); + if( !file.is_open() ) { + std::cerr << "Failed to open vector file: " << filename << "\n"; + return; + } + std::string line; + while( std::getline( file, line ) ) { + if( line.empty() ) continue; // skip empty lines + std::istringstream iss( line ); + Dtype v; + if( !(iss >> v) ) { + throw std::runtime_error( "Failed to parse line in vector file" ); + } + data.push_back( v ); + } + } catch( std::exception &e ) { + std::cerr << "I/O program failed: " << e.what() << "\n"; + return; + } +} + + +template< typename Dtype > +void read_vector_data_from_array( + const std::vector &array, std::vector &data +) { + // Implementation for reading vector data from array + try { + for (size_t i = 0; i < array.size(); ++i) { + data.push_back(array[i]); + } + } catch (const std::exception &e) { + std::cerr << "Failed to read vector data from array: " << e.what() << "\n"; + return; + } +} + +void ioProgram( const struct input &data_in, bool &success ) { + using namespace test_data; + success = false; + // Parse and store matrix in singleton class + auto &Jdata = std::get<2>(Storage::getData()); + // Read and store h vector in singleton class + auto &h = std::get<3>(Storage::getData()); + auto &x = std::get<4>(Storage::getData()); + auto &y = std::get<5>(Storage::getData()); + + if(data_in.use_default_data){ + // if no file provided, use default data from file_content + // TODO: for now use matrix data file + //read_matrix_data( data_in.filename_Jmatrix, Jdata, data_in.direct ); + read_matrix_data_from_array( test_data::j_matrix_data, Jdata ); + read_vector_data_from_array( test_data::h_array_data, h ); + read_vector_data_from_array( test_data::x_array_data, x ); + read_vector_data_from_array( test_data::y_array_data, y ); + + } else { + // read from files if provided + read_matrix_data( data_in.filename_Jmatrix, Jdata, data_in.direct ); + read_vector_data( data_in.filename_h, h ); + read_vector_data( data_in.filename_x, x ); + read_vector_data( data_in.filename_y, y ); + } + success = true; } @@ -194,19 +318,14 @@ void grbProgram( timer.reset(); /* --- Problem setup --- */ - grb::Vector h( N ); - grb::Vector x0( N ), y0( N ); // initialy + const size_t n = std::get<0>(Storage::getData()); + std::cout << "n = " << n << std::endl; + grb::Vector h( n ); + grb::Vector x0( n ), y0( n ); // initialy // ... populate J with test (random) values grb::RC rc = grb::SUCCESS; - // rc = rc ? rc : buildMatrixUnique( J, &( i_arr[ 0 ] ), &( j_arr[ 0 ] ), &( v_arr[ 0 ] ), Nz, grb::SEQUENTIAL ); - // if(rc != grb::SUCCESS) { - // std::cerr << "matrix build failed\n"; - // return; - // } - - // load into GraphBLAS - const size_t n = std::get<0>(Storage::getData()); + // load into GraphBLAS grb::Matrix J( n, n ); { const auto &data = std::get<2>(Storage::getData()); @@ -238,12 +357,47 @@ void grbProgram( out.error_code = 5; return; } + + // print matrix + if( s == 0 ) { + std::cout << "Matrix J:\n"; + print_matrix( J); + } } - // Fill h, x0, y0 with random values using buildVector - rc = rc ? rc : buildVector(h, h_arr, h_arr + N, grb::SEQUENTIAL); - rc = rc ? rc : buildVector(x0, x_arr, x_arr + N, grb::SEQUENTIAL); - rc = rc ? rc : buildVector(y0, y_arr, y_arr + N, grb::SEQUENTIAL); + // build vector h with data from singleton + { + const auto &h_data = std::get<3>(Storage::getData()); + rc = rc ? rc : buildVector( + h, + h_data.cbegin(), + h_data.cend(), + SEQUENTIAL + ); + } + + // build vector x with data from singleton + { + const auto &x_data = std::get<4>(Storage::getData()); + rc = rc ? rc : buildVector( + x0, + x_data.cbegin(), + x_data.cend(), + SEQUENTIAL + ); + } + + // build vector y with data from singleton + { + const auto &y_data = std::get<5>(Storage::getData()); + rc = rc ? rc : buildVector( + y0, + y_data.cbegin(), + y_data.cend(), + SEQUENTIAL + ); + } + if(rc != grb::SUCCESS) { std::cerr << "Vector build failed\n"; return; @@ -252,19 +406,20 @@ void grbProgram( // energies is array of length num_iters, initialized to 0 std::vector< IOType > energies( num_iters, 0 ); - grb::Vector< IOType > Jx( N ); - grb::Vector< IOType > temp( N ); - grb::Vector< JType > temp_int( N ); - grb::Vector< bool > mask( N ); - grb::Matrix< JType > J2( N, N ); + grb::Vector< IOType > Jx( n ); + grb::Vector< IOType > temp( n ); + grb::Vector< JType > temp_int( n ); + grb::Vector< bool > mask( n ); + grb::Matrix< JType > J2( n, n ); rc = rc ? rc : grb::resize( J2, grb::nnz(J) ); if(rc != grb::SUCCESS) { std::cerr << "Matrix resize failed for J2\n"; return; } - grb::Vector< JType > sol( N ); - grb::Vector< JType > sol_ref( N ); - rc = rc ? rc : buildVector(sol_ref, sol_ref_data, sol_ref_data + N, grb::SEQUENTIAL); + grb::Vector< JType > sol( n ); + // TODO: enable sol_ref + // grb::Vector< JType > sol_ref( N ); + // rc = rc ? rc : buildVector(sol_ref, sol_ref_data, sol_ref_data + N, grb::SEQUENTIAL); out.times.preamble = timer.time(); @@ -353,8 +508,9 @@ void grbProgram( return; } - out.pinnedRefSolutionVector = std::unique_ptr< PinnedVector< JType > >( - new PinnedVector< JType >( sol_ref, SEQUENTIAL ) ); + // TODO: enable sol_ref + // out.pinnedRefSolutionVector = std::unique_ptr< PinnedVector< JType > >( + // new PinnedVector< JType >( sol_ref, SEQUENTIAL ) ); // output out.pinnedSolutionVector = std::unique_ptr< PinnedVector< JType > >( @@ -371,9 +527,9 @@ void grbProgram( } else { // print all energies for (std::size_t i = 0; i < num_iters; ++i) { -#ifdef DEBUG_IMSB +//#ifdef DEBUG_IMSB std::cout << "Energy at iteration " << i << " = " << energies[i] << '\n'; -#endif +//#endif if( energies[i] != energies_ref[i]) { #ifdef DEBUG_IMSB std::cerr << "Error: Energy at iteration " << i << " does not match reference value.\n"; @@ -402,7 +558,11 @@ int main( int argc, char ** argv ) { // the input struct struct input in; - in.filename = "/home/d/Scratch/SA/ising_machine.mtx"; + in.filename_Jmatrix = "/home/d/Scratch/SA/ising_machine.mtx"; + in.filename_h = "/home/d/Scratch/SA/h_vector.dat"; + in.filename_x = "/home/d/Scratch/SA/x_vector.dat"; + in.filename_y = "/home/d/Scratch/SA/y_vector.dat"; + in.use_default_data = true; // use default data from files // get inner number of iterations in.rep = grb::config::BENCHMARKING::inner(); @@ -479,7 +639,8 @@ int main( int argc, char ** argv ) { if( solution.size() > 0 ) { print_vector( solution, 30, "SOLUTION" ); // expected solution from sol_ref_data - print_vector( solution_ref, 30, "EXPECTED SOLUTION" ); + // TODO: enable sol_ref + // print_vector( solution_ref, 30, "EXPECTED SOLUTION" ); } else { std::cerr << "ERROR: solution contains no values" << std::endl; } From e894de5ad040970541766b3c05a476c8b84e1247 Mon Sep 17 00:00:00 2001 From: Denis Jelovina Date: Fri, 22 Aug 2025 14:06:34 +0200 Subject: [PATCH 16/18] Enhance is_complex class to support arithmetic types and update Ising machine test to include verification against reference solutions --- include/graphblas/utils/iscomplex.hpp | 118 +++---- tests/smoke/CMakeLists.txt | 2 +- tests/smoke/ising_machine_sb.cpp | 474 ++++++++++++++++++++------ tests/smoke/smoketests.sh | 13 + 4 files changed, 439 insertions(+), 168 deletions(-) diff --git a/include/graphblas/utils/iscomplex.hpp b/include/graphblas/utils/iscomplex.hpp index 3a272d233..7d2e4bff3 100644 --- a/include/graphblas/utils/iscomplex.hpp +++ b/include/graphblas/utils/iscomplex.hpp @@ -49,65 +49,65 @@ namespace grb { template< typename C > class is_complex { - static_assert( - std::is_floating_point< C >::value, - "is_complex: C is not a floating point type" - ); - - public: - - /** - * If \a value is false, the type will be \a C. - * If \a value is true, the type will be C::value_type. - */ - typedef C type; - - /** Whether the type \a C is std::complex */ - static constexpr const bool value = false; - - /** - * @returns The conjugate of a given value if \a C is a complex type, or - * the given value if \a C is not complex. - */ - static C conjugate( const C &x ) noexcept { - return x; - } - - /** - * @returns The absolute value of a given value if \a C is a complex type, - * or the given value if \a C is not complex. - */ - static C modulus( const C &x ) noexcept { - return( x > 0 ? x : -x ); - } - - /** - * @returns The absolute value squared of a given value. - */ - static C norm( const C &x ) noexcept { - return x * x; - } - - /** - * @returns The polar coordinates of a given value. - * - * The first argument in the returned pair is the magnitude, while the - * second is the phase. - */ - static std::pair< C, C > polar( const C &x ) noexcept { - std::pair< C, C > ret{ std::abs( x ), x < 0 ? M_PI : 0 }; - return ret; - } - - /** - * @returns The multiplicative inverse of a given value. - */ - static C inverse( const C &x ) noexcept { - constexpr C one = 1; - return one / x; - } - - }; + static_assert( + std::is_arithmetic< C >::value, + "is_complex: C is not a numerical (arithmetic) type" + ); + + public: + + /** + * If \a value is false, the type will be \a C. + * If \a value is true, the type will be C::value_type. + */ + typedef C type; + + /** Whether the type \a C is std::complex */ + static constexpr const bool value = false; + + /** + * @returns The conjugate of a given value if \a C is a complex type, or + * the given value if \a C is not complex. + */ + static C conjugate( const C &x ) noexcept { + return x; + } + + /** + * @returns The absolute value of a given value if \a C is a complex type, + * or the given value if \a C is not complex. + */ + static C modulus( const C &x ) noexcept { + return( x > 0 ? x : -x ); + } + + /** + * @returns The absolute value squared of a given value. + */ + static C norm( const C &x ) noexcept { + return x * x; + } + + /** + * @returns The polar coordinates of a given value. + * + * The first argument in the returned pair is the magnitude, while the + * second is the phase. + */ + static std::pair< C, C > polar( const C &x ) noexcept { + std::pair< C, C > ret{ std::abs( x ), x < 0 ? M_PI : 0 }; + return ret; + } + + /** + * @returns The multiplicative inverse of a given value. + */ + static C inverse( const C &x ) noexcept { + constexpr C one = 1; + return one / x; + } + + }; /** \internal The specialisation for std::complex types. */ template< typename T > diff --git a/tests/smoke/CMakeLists.txt b/tests/smoke/CMakeLists.txt index 9ce5553c3..c1da83b8b 100644 --- a/tests/smoke/CMakeLists.txt +++ b/tests/smoke/CMakeLists.txt @@ -145,7 +145,7 @@ add_grb_executables( conjugate_gradient_complex conjugate_gradient.cpp ) add_grb_executables( ising_machine_sb ising_machine_sb.cpp - BACKENDS reference reference_omp nonblocking + BACKENDS reference reference_omp bsp1d hybrid hyperdags nonblocking ADDITIONAL_LINK_LIBRARIES test_utils_headers ) diff --git a/tests/smoke/ising_machine_sb.cpp b/tests/smoke/ising_machine_sb.cpp index f4fa2f105..6ba226b4a 100644 --- a/tests/smoke/ising_machine_sb.cpp +++ b/tests/smoke/ising_machine_sb.cpp @@ -18,7 +18,7 @@ #include #include #include - +#include #include @@ -68,7 +68,8 @@ typedef grb::utils::Singleton< std::vector, // matrix data std::vector, // h vector std::vector, // x vector - std::vector // y vector + std::vector, // y vector + std::vector // sol_ref vector // Add more types as needed > > Storage; @@ -91,33 +92,23 @@ namespace test_data { -0.0261, 0.0018, -0.0710, 0.0507, -0.0483 }; const std::vector sol_ref_data = { 1, -1, 1, 1, 1, -1, -1, -1, 1, 1 }; + // matrix in format of list of nested pairs ((i, j), value) const std::vector, JType > > j_matrix_data = { {{1, 1}, -1}, {{2, 2}, -1}, {{3, 1}, 1}, {{3, 3}, -1}, {{4, 1}, 1}, {{4, 4}, 1}, {{5, 1}, -1}, {{5, 2}, -1}, {{5, 3}, -1}, {{5, 4}, 1}, {{5, 5}, 1}, {{6, 1}, -1 }, {{6, 2}, 1}, {{6, 4}, 1}, {{6, 6}, -1}, {{7, 2}, 1}, {{7, 3}, -1}, - {{7, 4}, -1}, {{7, 6}, -1}, - {{7, 7}, -1}, {{8, 2}, 1}, - {{8, 4}, -1}, {{8, 8}, -1}, - {{9, 4}, 1}, {{9, 6}, -1}, - {{9, 8}, 1}, {{9, 9}, 1}, - {{10, 2}, 1}, {{10, 3}, -1}, - {{10, 4}, 1}, {{10, 8}, -1}, - {{10, 10}, 1} - , - // since matrix is symmetric, we can add the symmetric entries - {{1, 3}, 1}, - {{1, 4}, 1}, {{1, 5}, -1}, {{2, 5}, -1}, - {{3, 5}, -1}, {{4, 5}, 1}, - {{1, 6}, -1 }, {{2, 6}, 1}, {{4, 6}, 1}, - {{2, 7}, 1}, {{3, 7}, -1}, - {{4, 7}, -1}, {{6, 7}, -1}, - {{2, 8}, 1}, - {{4, 8}, -1}, - {{4, 9}, 1}, {{6, 9}, -1}, - {{8, 9}, 1}, - {{2, 10}, 1}, {{3, 10}, -1}, + {{7, 4}, -1}, {{7, 6}, -1}, {{7, 7}, -1}, {{8, 2}, 1}, + {{8, 4}, -1}, {{8, 8}, -1}, {{9, 4}, 1}, {{9, 6}, -1}, + {{9, 8}, 1}, {{9, 9}, 1}, {{10, 2}, 1}, {{10, 3}, -1}, + {{10, 4}, 1}, {{10, 8}, -1}, {{10, 10}, 1}, + // since matrix is symmetric, we add the symmetric entries + {{1, 3}, 1}, {{1, 4}, 1}, {{1, 5}, -1}, {{2, 5}, -1}, + {{3, 5}, -1}, {{4, 5}, 1}, {{1, 6}, -1 }, {{2, 6}, 1}, {{4, 6}, 1}, + {{2, 7}, 1}, {{3, 7}, -1}, {{4, 7}, -1}, {{6, 7}, -1}, + {{2, 8}, 1}, {{4, 8}, -1}, {{4, 9}, 1}, {{6, 9}, -1}, + {{8, 9}, 1}, {{2, 10}, 1}, {{3, 10}, -1}, {{4, 10}, 1}, {{8, 10}, -1} }; @@ -134,7 +125,6 @@ namespace test_data { -11, -11, -11, -11, -11, -11, -11, -11, -11 }; - const IOType p0 = 0.; const IOType p1 = 1.1; const IOType dt = 0.25; @@ -150,6 +140,17 @@ struct input { std::string filename_y; bool direct; size_t rep; + size_t outer; // number of outer repetitions for benchmarking + // number of iterations for Ising machine SB + size_t num_iters; + // p0, p1, dt parameters for Ising machine SB + IOType p0; + IOType p1; + IOType dt; + // whether to verify output against reference data + bool verify; + // filename of reference solution vector (not implemented) + std::string filename_ref_solution; }; bool input::use_default_data = false; @@ -189,9 +190,11 @@ void read_matrix_data(const std::string &filename, std::vector &data, boo ++it ) { data.push_back( Dtype( *it ) ); - // print last data element - std::cout << "read_matrix_data_from_file: " << data.back().first.first << ", " +#ifdef DEBUG_IMSB + // print last data element from std::vector data + std::cout << "read_matrix_data: " << data.back().first.first << ", " << data.back().first.second << ", " << data.back().second << "\n"; +#endif } } catch( std::exception &e ) { std::cerr << "I/O program failed: " << e.what() << "\n"; @@ -207,18 +210,14 @@ void read_matrix_data_from_array( // Implementation for reading matrix data from array try { for (const auto &entry : array) { - //std::cout << "read_matrix_data_from_array: " << entry.first.first << ", " << entry.first.second << ", " << entry.second << "\n"; - // since - // data.push_back( Dtype( *it ) ); - // print last data element - // std::cout << "read_matrix_data: " << it->first.first << ", " << it->first.second << ", " << it->second << "\n"; - // we want the same here data.emplace_back( NonzeroT( entry.first.first-1, entry.first.second-1, entry.second ) ); +#ifdef DEBUG_IMSB // print last data element from std::vector data std::cout << "read_matrix_data_from_array: " << data.back().first.first << ", " << data.back().first.second << ", " << data.back().second << "\n"; +#endif } std::get<0>(Storage::getData()) = test_data::n; std::get<1>(Storage::getData()) = data.size(); @@ -270,6 +269,7 @@ void read_vector_data_from_array( } void ioProgram( const struct input &data_in, bool &success ) { + using namespace test_data; success = false; // Parse and store matrix in singleton class @@ -278,22 +278,28 @@ void ioProgram( const struct input &data_in, bool &success ) { auto &h = std::get<3>(Storage::getData()); auto &x = std::get<4>(Storage::getData()); auto &y = std::get<5>(Storage::getData()); + auto &sol = std::get<6>(Storage::getData()); if(data_in.use_default_data){ // if no file provided, use default data from file_content - // TODO: for now use matrix data file - //read_matrix_data( data_in.filename_Jmatrix, Jdata, data_in.direct ); read_matrix_data_from_array( test_data::j_matrix_data, Jdata ); read_vector_data_from_array( test_data::h_array_data, h ); read_vector_data_from_array( test_data::x_array_data, x ); read_vector_data_from_array( test_data::y_array_data, y ); - + read_vector_data_from_array( test_data::sol_ref_data, sol ); } else { // read from files if provided read_matrix_data( data_in.filename_Jmatrix, Jdata, data_in.direct ); read_vector_data( data_in.filename_h, h ); read_vector_data( data_in.filename_x, x ); read_vector_data( data_in.filename_y, y ); + if(data_in.verify) { + if(data_in.filename_ref_solution.empty()) { + std::cerr << "Reference solution file not provided for verification\n"; + return; + } + } + read_vector_data( data_in.filename_ref_solution, sol ); } success = true; @@ -311,15 +317,10 @@ void grbProgram( grb::utils::Timer timer; timer.reset(); - // TODO: IO goes here - - // I/O done - out.times.io = timer.time(); - timer.reset(); /* --- Problem setup --- */ const size_t n = std::get<0>(Storage::getData()); - std::cout << "n = " << n << std::endl; + std::cout << "problem size n = " << n << "\n"; grb::Vector h( n ); grb::Vector x0( n ), y0( n ); // initialy // ... populate J with test (random) values @@ -358,11 +359,12 @@ void grbProgram( return; } - // print matrix - if( s == 0 ) { - std::cout << "Matrix J:\n"; - print_matrix( J); - } +#ifdef DEBUG_IMSB + if( s == 0 ) { + std::cout << "Matrix J:\n"; + print_matrix( J); + } +#endif } // build vector h with data from singleton @@ -417,11 +419,21 @@ void grbProgram( return; } grb::Vector< JType > sol( n ); - // TODO: enable sol_ref - // grb::Vector< JType > sol_ref( N ); - // rc = rc ? rc : buildVector(sol_ref, sol_ref_data, sol_ref_data + N, grb::SEQUENTIAL); - - out.times.preamble = timer.time(); + grb::Vector< JType > sol_ref( n ); + if(data_in.verify) { + // build vector sol_ref with data from singleton + const auto &sol_ref_data = std::get<6>(Storage::getData()); + rc = rc ? rc : buildVector( + sol_ref, + sol_ref_data.cbegin(), + sol_ref_data.cend(), + grb::SEQUENTIAL + ); + } + if(rc != grb::SUCCESS) { + std::cerr << "Vector build failed\n"; + return; + } rc = rc ? rc : wait(); out.times.preamble = timer.time(); @@ -477,7 +489,7 @@ void grbProgram( J2, Jx, temp, temp_int, mask, sol, out.iterations ); } - if( Properties<>::isNonblockingExecution ) { + if( grb::Properties<>::isNonblockingExecution ) { rc = rc ? rc : wait(); } } @@ -508,9 +520,8 @@ void grbProgram( return; } - // TODO: enable sol_ref - // out.pinnedRefSolutionVector = std::unique_ptr< PinnedVector< JType > >( - // new PinnedVector< JType >( sol_ref, SEQUENTIAL ) ); + out.pinnedRefSolutionVector = std::unique_ptr< PinnedVector< JType > >( + new PinnedVector< JType >( sol_ref, SEQUENTIAL ) ); // output out.pinnedSolutionVector = std::unique_ptr< PinnedVector< JType > >( @@ -525,23 +536,27 @@ void grbProgram( if( rc != grb::SUCCESS ) { std::cerr << "bSB returned error code " << rc << '\n'; } else { - // print all energies - for (std::size_t i = 0; i < num_iters; ++i) { -//#ifdef DEBUG_IMSB - std::cout << "Energy at iteration " << i << " = " << energies[i] << '\n'; -//#endif - if( energies[i] != energies_ref[i]) { -#ifdef DEBUG_IMSB - std::cerr << "Error: Energy at iteration " << i << " does not match reference value.\n"; - std::cerr << "Expected: " << energies_ref[i] << ", got: " << energies[i] << '\n'; -#endif - out.error_code = 40; - return ; - } - } - std::cout << "All energies match reference values.\n"; - std::cout << "TEST OK\n"; - } + if (data_in.verify) { + // print all energies + for (std::size_t i = 0; i < num_iters; ++i) { + #ifdef DEBUG_IMSB + std::cout << "Energy at iteration " << i << " = " << energies[i] << '\n'; + #endif + if( energies[i] != energies_ref[i]) { + #ifdef DEBUG_IMSB + std::cerr << "Error: Energy at iteration " << i << " does not match reference value.\n"; + std::cerr << "Expected: " << energies_ref[i] << ", got: " << energies[i] << '\n'; + #endif + out.error_code = 40; + return ; + } + } + std::cout << "All energies match reference values.\n"; + std::cout << "TEST OK\n"; + } else { + std::cout << "No verification performed (verification disabled).\n"; + } +} // set error code out.error_code = rc; @@ -549,40 +564,265 @@ void grbProgram( } -int main( int argc, char ** argv ) { - // TODO: add argument parsing for input file, direct/indirect addressing, etc. - // for now, just print the executable name - (void) argc; // unused - (void) argv; // unused - std::cout << "Test executable: " << argv[ 0 ] << std::endl; - - // the input struct - struct input in; - in.filename_Jmatrix = "/home/d/Scratch/SA/ising_machine.mtx"; - in.filename_h = "/home/d/Scratch/SA/h_vector.dat"; - in.filename_x = "/home/d/Scratch/SA/x_vector.dat"; - in.filename_y = "/home/d/Scratch/SA/y_vector.dat"; - in.use_default_data = true; // use default data from files +// supported command line arguments +void printhelp( char *progname ) { + std::cout << " Use: \n"; + std::cout << progname << " [--use-default-data] [--j-matrix-fname STR] [--h-fname STR] [--x-fname STR] [--y-fname STR] [--no-direct] [--num-iters INT] [--p0 FLOAT] [--p1 FLOAT] [--dt FLOAT] [--test-rep INT] [--test-outer-rep INT] [--verify] [--ref-solution-fname STR]\n"; + std::cout << "\n"; + std::cout << " --use-default-data (no argument): use hardcoded default data from the test_data namespace for internal tests\n"; + // input data parameters (mandatory if --use-default-data is not used) + std::cout << "\n"; + std::cout << "input data parameters (mandatory if --use-default-data is not used):\n"; + std::cout << " --j-matrix-fname STR: filename of J matrix in matrix market format\n"; + std::cout << " --h-fname STR: filename of h vector, where vector elements are stored line-by-line\n"; + std::cout << " --x-fname STR: filename of x vector, where vector elements are stored line-by-line\n"; + std::cout << " --y-fname STR: filename of y vector, where vector elements are stored line-by-line\n"; + std::cout << " --ref-solution-fname STR: mandatory if --verify is used and --use-default-data is not used. (not implemented) filename of reference solution vector, where vector elements are stored line-by-line\n"; + std::cout << " --no-direct (no argument): disable direct addressing\n"; + // either --use-default-data or input data parameters must be provided + std::cout << "\n"; + std::cout << "either --use-default-data or input data parameters must be provided\n"; + // solver parameters can be set via command-line arguments or will use defaults if not provided + std::cout << "\n"; + std::cout << "solver parameters (can be set via command-line arguments or will use defaults):\n"; + std::cout << " --p0 FLOAT: p0 parameter, default " << test_data::p0 << "\n"; + std::cout << " --p1 FLOAT: p1 parameter, default " << test_data::p1 << "\n"; + std::cout << " --dt FLOAT: dt parameter, default " << test_data::dt << "\n"; + std::cout << " --num-iters INT: number of iterations, default " << test_data::num_iters << " \n"; + // performance testing parameters (optional) + std::cout << "\n"; + std::cout << "performance testing parameters (optional):\n"; + std::cout << " --test-rep INT: consecutive test inner algorithm repetitions, default 1\n"; + std::cout << " --test-outer-rep INT: consecutive test outer (including IO) algorithm repetitions, default 1\n"; + // numerical verification parameters (optional) + std::cout << "\n"; + std::cout << "numerical verification parameters (optional):\n"; + std::cout << " --verify (no argument): verify output against reference solution\n"; + // other parameters + std::cout << "\n"; + std::cout << "other parameters:\n"; + std::cout << " --help, -h, -? (no argument): print this help message\n"; +} +bool parse_arguments( + input &in, + int argc, + char **argv +) { + in.filename_Jmatrix.clear(); + in.filename_h.clear(); + in.filename_x.clear(); + in.filename_y.clear(); + in.filename_ref_solution.clear(); + in.direct = true; // get inner number of iterations in.rep = grb::config::BENCHMARKING::inner(); - // get outer number of iterations - size_t outer = grb::config::BENCHMARKING::outer(); + in.outer = grb::config::BENCHMARKING::outer(); + in.p0 = 0.0; + in.p1 = 0.0; + in.dt = 0.0; + in.num_iters = 0; + in.verify = false; + + bool jmatrix_set = false, h_set = false, x_set = false, y_set = false, refsol_set = false; + + for( int i = 1; i < argc; ++i ) { + std::string arg( argv[ i ] ); + if( arg == "--use-default-data" ) { + in.use_default_data = true; + } else if( arg == "--j-matrix-fname" ) { + if( i + 1 < argc ) { + in.filename_Jmatrix = std::string( argv[ ++i ] ); + in.use_default_data = false; + jmatrix_set = true; + } else { + std::cerr << "--j-matrix-fname requires an argument\n"; + return false; + } + } else if( arg == "--h-fname" ) { + if( i + 1 < argc ) { + in.filename_h = std::string( argv[ ++i ] ); + in.use_default_data = false; + h_set = true; + } else { + std::cerr << "--h-fname requires an argument\n"; + return false; + } + } else if( arg == "--x-fname" ) { + if( i + 1 < argc ) { + in.filename_x = std::string( argv[ ++i ] ); + in.use_default_data = false; + x_set = true; + } else { + std::cerr << "--x-fname requires an argument\n"; + return false; + } + } else if( arg == "--y-fname" ) { + if( i + 1 < argc ) { + in.filename_y = std::string( argv[ ++i ] ); + in.use_default_data = false; + y_set = true; + } else { + std::cerr << "--y-fname requires an argument\n"; + return false; + } + } else if( arg == "--no-direct" ) { + in.direct = false; + } else if( arg == "--test-rep" ) { + if( i + 1 < argc ) { + int r = atoi( argv[ ++i ] ); + if( r < 0 ) { + std::cerr << "--test-rep requires a non-negative integer argument\n"; + return false; + } + in.rep = static_cast< size_t >( r ); + } else { + std::cerr << "--test-rep requires an argument\n"; + return false; + } + } else if( arg == "--test-outer-rep" ) { + if( i + 1 < argc ) { + int r = atoi( argv[ ++i ] ); + if( r < 1 ) { + std::cerr << "--test-outer-rep requires a positive integer argument\n"; + return false; + } + in.outer = static_cast< size_t >( r ); + } else { + std::cerr << "--test-outer-rep requires an argument\n"; + return false; + } + } else if( arg == "--p0" ) { + if( i + 1 < argc ) { + in.p0 = atof( argv[ ++i ] ); + } else { + std::cerr << "--p0 requires a FLOAT argument\n"; + return false; + } + } else if( arg == "--p1" ) { + if( i + 1 < argc ) { + in.p1 = atof( argv[ ++i ] ); + } else { + std::cerr << "--p1 requires a FLOAT argument\n"; + return false; + } + } else if( arg == "--dt" ) { + if( i + 1 < argc ) { + in.dt = atof( argv[ ++i ] ); + } else { + std::cerr << "--dt requires a FLOAT argument\n"; + return false; + } + } else if( arg == "--num-iters" ) { + if( i + 1 < argc ) { + int niter = atoi( argv[ ++i ] ); + if( niter < 1 ) { + std::cerr << "--num-iters requires a positive integer argument\n"; + return false; + } + in.num_iters = static_cast< size_t >( niter ); + } else { + std::cerr << "--num-iters requires an argument\n"; + return false; + } + } else if( arg == "--verify" ) { + in.verify = true; + } else if( arg == "--ref-solution-fname" ) { + if( i + 1 < argc ) { + in.filename_ref_solution = std::string( argv[ ++i ] ); + refsol_set = true; + } else { + std::cerr << "--ref-solution-fname requires an argument\n"; + return false; + } + } else if( arg == "--help" || arg == "-h" || arg == "-?" ) { + printhelp( argv[ 0 ] ); + return false; + } else { + std::cerr << "unknown command line argument \"" << arg << "\"\n"; + return false; + } + } + + // Check that either --use-default-data or all input data parameters are provided + if( !in.use_default_data ) { + if( !(jmatrix_set && h_set && x_set && y_set) ) { + std::cerr << "Error: Either --use-default-data or all input data parameters (--j-matrix-fname, --h-fname, --x-fname, --y-fname) must be provided.\n"; + return false; + } + } + + // Ensure reference solution filename is provided if verification is requested and not using default data + if( in.verify && !in.use_default_data ) { + if( !refsol_set || in.filename_ref_solution.empty() ) { + std::cerr << "Error: --ref-solution-fname STR is mandatory if --verify is used and --use-default-data is not used.\n"; + return false; + } + } + + // set defaults for solver parameters if not set + if( in.p0 == 0.0 && in.p1 == 0.0 && in.dt == 0.0 ) { + in.p0 = test_data::p0; + in.p1 = test_data::p1; + in.dt = test_data::dt; + } + if( in.num_iters == 0 ) { + in.num_iters = test_data::num_iters; + } + return true; +} + + +void print_cmd_paramaters( const input &in) { + // print paramerters + std::cout << "Parameters:\n"; + std::cout << " --use-default-data: " << (in.use_default_data ? "true" : "false") << "\n"; + if( !in.use_default_data ) { + std::cout << " --j-matrix-fname: " << in.filename_Jmatrix << "\n"; + std::cout << " --h-fname: " << in.filename_h << "\n"; + std::cout << " --x-fname: " << in.filename_x << "\n"; + std::cout << " --y-fname: " << in.filename_y << "\n"; + } + std::cout << " --no-direct: " << (in.direct ? "false" : "true") << "\n"; + std::cout << " --p0: " << in.p0 << "\n"; + std::cout << " --p1: " << in.p1 << "\n"; + std::cout << " --dt: " << in.dt << "\n"; + std::cout << " --num-iters: " << in.num_iters << "\n"; + std::cout << " --test-rep: " << in.rep << "\n"; + std::cout << " --test-outer-rep: " << in.outer << "\n"; + std::cout << " --verify: " << (in.verify ? "true" : "false") << "\n"; + if( in.verify ) { + std::cout << " --ref-solution-fname: " << in.filename_ref_solution << "\n"; + } + std::cout << "\n"; +} + +int main( int argc, char ** argv ) { + std::cout << "Test executable: " << argv[ 0 ] << std::endl; - // check for verification of the output - // bool verification = false; + input in; + output out; - std::cout << "Executable called with parameters " - << "inner repititions = " << in.rep << ", " - << "outer reptitions = " << outer << ", " - << std::endl; + if( !parse_arguments( in, argc, argv ) ) { + std::cerr << "error parsing command line arguments\n"; + printhelp( argv[0] ); + return 1; + } + +#ifdef DEBUG_IMSB + // print paramerters + print_cmd_paramaters( in ); +#endif // set standard exit code grb::RC rc = SUCCESS; // launch I/O { + grb::utils::Timer timer; + timer.reset(); + bool success; grb::Launcher< AUTOMATIC > launcher; rc = launcher.exec( &ioProgram, in, success, true ); @@ -595,11 +835,11 @@ int main( int argc, char ** argv ) { std::cerr << "I/O program caught an exception\n"; return 77; } + // I/O done + out.times.io = timer.time(); + timer.reset(); } - // the output struct - struct output out; - // launch estimator (if requested) if( in.rep == 0 ) { grb::Launcher< AUTOMATIC > launcher; @@ -617,7 +857,7 @@ int main( int argc, char ** argv ) { // launch benchmark if( rc == SUCCESS ) { grb::Benchmarker< AUTOMATIC > benchmarker; - rc = benchmarker.exec( &grbProgram, in, out, 1, outer, true ); + rc = benchmarker.exec( &grbProgram, in, out, 1, in.outer, true ); } if( rc != SUCCESS ) { std::cerr << "benchmarker.exec returns with non-SUCCESS error code " @@ -630,6 +870,7 @@ int main( int argc, char ** argv ) { std::cout << "Error code is " << out.error_code << ".\n"; + // inspect output vector if( !(out.pinnedSolutionVector) ) { std::cerr << "no output vector to inspect" << std::endl; } else { @@ -638,21 +879,38 @@ int main( int argc, char ** argv ) { std::cout << "Size of x is " << solution.size() << std::endl; if( solution.size() > 0 ) { print_vector( solution, 30, "SOLUTION" ); - // expected solution from sol_ref_data - // TODO: enable sol_ref - // print_vector( solution_ref, 30, "EXPECTED SOLUTION" ); } else { std::cerr << "ERROR: solution contains no values" << std::endl; } + if( in.verify && (solution_ref.size() > 0) ) { + print_vector( solution_ref, 30, "REFERENCE SOLUTION" ); + } } - - if( out.error_code != 0 ) { - std::cerr << std::flush; - std::cout << "Test FAILED\n"; - } - - std::cout << "Test OK\n"; + // verify output vector if requested + if( in.verify ) { + const PinnedVector< JType > &solution = *(out.pinnedSolutionVector); + const PinnedVector< JType > &solution_ref = *(out.pinnedRefSolutionVector); + assert(solution.size() == solution_ref.size()); + assert(solution.size() == std::get<0>(Storage::getData())); + JType norm2 = 0; + for( size_t i = 0; i < solution.size(); i++ ) { + JType diff = solution.getNonzeroValue(i) - solution_ref.getNonzeroValue(i); + norm2 += diff * diff; + } + out.error_code = ( norm2 == 0 ) ? 0 : 50; + if( out.error_code == 0 ) { + std::cout << "Output vector verificaton was successful!\n"; + std::cout << "Test OK\n"; + } else { + std::cerr << std::flush; + std::cout << "Verification FAILED\n"; + std::cout << "Test FAILED\n"; + } + } else { + std::cout << "Output vector verificaton was not requested!\n"; + std::cout << "Test OK\n"; + } // done return out.error_code; diff --git a/tests/smoke/smoketests.sh b/tests/smoke/smoketests.sh index ce6c3ce6b..77111a7ed 100755 --- a/tests/smoke/smoketests.sh +++ b/tests/smoke/smoketests.sh @@ -258,6 +258,19 @@ for BACKEND in ${BACKENDS[@]}; do fi echo " " + if [ -f "${TEST_BIN_DIR}/ising_machine_sb_${BACKEND}" ] + then + echo ">>> [x] [ ] Testing the Ising machine simulated bifurcation algorithm" + echo " for a small 10-spin problem with known solution. This test" + echo " verifies against a predefined solution vector. The test" + echo " employs the grb::Launcher in automatic mode. It uses" + echo " direct-mode file IO." + $runner ${TEST_BIN_DIR}/ising_machine_sb_${BACKEND} --use-default-data --verify &> ${TEST_OUT_DIR}/ising_machine_sb_${BACKEND}_${P}_${T}.log + head -1 ${TEST_OUT_DIR}/ising_machine_sb_${BACKEND}_${P}_${T}.log + grep 'Test OK' ${TEST_OUT_DIR}/ising_machine_sb_${BACKEND}_${P}_${T}.log || echo "Test FAILED" + echo " " + fi + NTEST=256 if [ -f "${TEST_BIN_DIR}/gmres_${BACKEND}" ] then From e6f3b247b7c32b53658e390e3b2b302aab506958 Mon Sep 17 00:00:00 2001 From: Denis Jelovina Date: Fri, 22 Aug 2025 14:19:48 +0200 Subject: [PATCH 17/18] Remove unused example from CMakeLists.txt --- examples/CMakeLists.txt | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 4835b4423..9a6affa1a 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -31,14 +31,3 @@ if( WITH_OMP_BACKEND ) add_dependencies( examples sp_reference_omp ) endif() -if( WITH_REFERENCE_BACKEND ) - add_executable( bbs_reference bbs.cpp ) - target_link_libraries( bbs_reference backend_reference common_flags ) - add_dependencies( examples bbs_reference ) -endif() - -if( WITH_REFERENCE_BACKEND ) - add_executable( bsb_Nsize_reference bsb_Nsize.cpp ) - target_link_libraries( bsb_Nsize_reference backend_reference common_flags ) - add_dependencies( examples bsb_Nsize_reference ) -endif() \ No newline at end of file From 2782a1dff426d990ae51a5b5d9062e557cd0d755 Mon Sep 17 00:00:00 2001 From: Denis Jelovina Date: Fri, 22 Aug 2025 15:03:32 +0200 Subject: [PATCH 18/18] Refactor is_complex class to improve code formatting and readability --- include/graphblas/utils/iscomplex.hpp | 118 +++++++++++++------------- 1 file changed, 59 insertions(+), 59 deletions(-) diff --git a/include/graphblas/utils/iscomplex.hpp b/include/graphblas/utils/iscomplex.hpp index 7d2e4bff3..1e5b38866 100644 --- a/include/graphblas/utils/iscomplex.hpp +++ b/include/graphblas/utils/iscomplex.hpp @@ -49,65 +49,65 @@ namespace grb { template< typename C > class is_complex { - static_assert( - std::is_arithmetic< C >::value, - "is_complex: C is not a numerical (arithmetic) type" - ); - - public: - - /** - * If \a value is false, the type will be \a C. - * If \a value is true, the type will be C::value_type. - */ - typedef C type; - - /** Whether the type \a C is std::complex */ - static constexpr const bool value = false; - - /** - * @returns The conjugate of a given value if \a C is a complex type, or - * the given value if \a C is not complex. - */ - static C conjugate( const C &x ) noexcept { - return x; - } - - /** - * @returns The absolute value of a given value if \a C is a complex type, - * or the given value if \a C is not complex. - */ - static C modulus( const C &x ) noexcept { - return( x > 0 ? x : -x ); - } - - /** - * @returns The absolute value squared of a given value. - */ - static C norm( const C &x ) noexcept { - return x * x; - } - - /** - * @returns The polar coordinates of a given value. - * - * The first argument in the returned pair is the magnitude, while the - * second is the phase. - */ - static std::pair< C, C > polar( const C &x ) noexcept { - std::pair< C, C > ret{ std::abs( x ), x < 0 ? M_PI : 0 }; - return ret; - } - - /** - * @returns The multiplicative inverse of a given value. - */ - static C inverse( const C &x ) noexcept { - constexpr C one = 1; - return one / x; - } - - }; + static_assert( + std::is_arithmetic< C >::value, + "is_complex: C is not a numerical (arithmetic) type" + ); + + public: + + /** + * If \a value is false, the type will be \a C. + * If \a value is true, the type will be C::value_type. + */ + typedef C type; + + /** Whether the type \a C is std::complex */ + static constexpr const bool value = false; + + /** + * @returns The conjugate of a given value if \a C is a complex type, or + * the given value if \a C is not complex. + */ + static C conjugate( const C &x ) noexcept { + return x; + } + + /** + * @returns The absolute value of a given value if \a C is a complex type, + * or the given value if \a C is not complex. + */ + static C modulus( const C &x ) noexcept { + return( x > 0 ? x : -x ); + } + + /** + * @returns The absolute value squared of a given value. + */ + static C norm( const C &x ) noexcept { + return x * x; + } + + /** + * @returns The polar coordinates of a given value. + * + * The first argument in the returned pair is the magnitude, while the + * second is the phase. + */ + static std::pair< C, C > polar( const C &x ) noexcept { + std::pair< C, C > ret{ std::abs( x ), x < 0 ? M_PI : 0 }; + return ret; + } + + /** + * @returns The multiplicative inverse of a given value. + */ + static C inverse( const C &x ) noexcept { + constexpr C one = 1; + return one / x; + } + + }; /** \internal The specialisation for std::complex types. */ template< typename T >