From 441540a99670b9134b759f8b34a03e7eeb19262b Mon Sep 17 00:00:00 2001 From: Michelle Date: Tue, 28 Jul 2026 17:32:49 +0200 Subject: [PATCH] Started on Random Library I made seed and uniform functions, but they are not omp friendly. They are runnning omp themself, but is not safe for threading/parallizing. Uniform is aproximated and is NOT validaded up against a real uniform distribution. --- include/random/seed.hpp | 110 +++++++ include/random/uniform.hpp | 196 +++++++++++ main.cpp | 649 +++++++++++++++++++++++++++++++++++++ src/io/print_tensor.cpp | 33 +- src/math/add.cpp | 4 +- src/random/seed.cpp | 122 +++++++ src/random/uniform.cpp | 260 +++++++++++++++ 7 files changed, 1359 insertions(+), 15 deletions(-) create mode 100644 include/random/seed.hpp create mode 100644 include/random/uniform.hpp create mode 100644 src/random/seed.cpp create mode 100644 src/random/uniform.cpp diff --git a/include/random/seed.hpp b/include/random/seed.hpp new file mode 100644 index 0000000..6b77268 --- /dev/null +++ b/include/random/seed.hpp @@ -0,0 +1,110 @@ +/**++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + * + * PANIC + * Portable Algorithms and Numerics In C++ + * + * Scientific computing from scratch, with feeling. + * + * Copyright (c) 2026 Michelle Bausager + * + * This file is part of PANIC. + * + * PANIC is free software licensed under the GNU General Public License v3.0 or later. + * You may redistribute and/or modify it under the terms of the GPL. + * + * PANIC is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; + * without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + * See the LICENSE file for the full license text. + * + * SPDX-License-Identifier: GPL-3.0-or-later + * + *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + * + * Project Name: PANIC + * Module Name: random + * File Name: seed.hpp + * Revision: 0.1.0 + * Date: 28-07-2026 + * Author: Michelle Bausager + * + * Description: + * Defines seed for use in other functions in random/ + * + *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ +#pragma once +//--------------------------------------------------------------------------------------------------------------------------- +// INCLUDE DESCRIPTION +//--------------------------------------------------------------------------------------------------------------------------- +#include // panic::uint_t, panic::int_t, and panic::real_t + + +//-------------------------------------------------------------------------------------------------------------------------- +// Function Name : panic::random::seed +// +// Description: +// base for random libary +//-------------------------------------------------------------------------------------------------------------------------- +namespace panic{ + namespace random{ + +/** + * @brief struct for seed object used in random/ + * + * The struct is used for PANIC random library. + */ +struct seed_t{ + + // Variable to store the value the random value is based on. + panic::types::uint_t value; + + /** + * @brief Empthy constructor + * + */ + seed_t(); + + /** + * @brief Contructor with seed. + * + * @param seed The seed that is used as the base for the object. + */ + seed_t(panic::types::uint_t seed); + + /** + * @brief Function to set the seed.. + * + * @param seed The seed that is used as the base for the object. + */ + bool set(panic::types::uint_t seed); + + // Creates a deterministic state from the seed and an index. + // This is OMP-friendly because it does not modify shared memory. + /** + * @brief Returns a random number based on seed and index + * + * @param index Value the random number will be drawn from. + * + * @Note this can be used with omp in a loop if the index in the loop + * is used as the input to this function. + */ + panic::types::uint_t state_at(panic::types::uint_t index) const; + + /** + * @brief Returns the base seed value. + * + */ + panic::types::uint_t get(); +}; + + + + } // namespace random +} // namespace panic + +//--------------------------------------------------------------------------------------------------------------------------- +// VARIABLE DESCRIPTION +//--------------------------------------------------------------------------------------------------------------------------- + +//--------------------------------------------------------------------------------------------------------------------------- +// FUNCTION PROTOTYPE +//--------------------------------------------------------------------------------------------------------------------------- \ No newline at end of file diff --git a/include/random/uniform.hpp b/include/random/uniform.hpp new file mode 100644 index 0000000..4951c2d --- /dev/null +++ b/include/random/uniform.hpp @@ -0,0 +1,196 @@ +/**++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + * + * PANIC + * Portable Algorithms and Numerics In C++ + * + * Scientific computing from scratch, with feeling. + * + * Copyright (c) 2026 Michelle Bausager + * + * This file is part of PANIC. + * + * PANIC is free software licensed under the GNU General Public License v3.0 or later. + * You may redistribute and/or modify it under the terms of the GPL. + * + * PANIC is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; + * without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + * See the LICENSE file for the full license text. + * + * SPDX-License-Identifier: GPL-3.0-or-later + * + *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + * + * Project Name: PANIC + * Module Name: random + * File Name: uniform.hpp + * Revision: 0.1.0 + * Date: 29-06-2026 + * Author: Michelle Bausager + * + * Description: + * Defindes the functions that returns values based on a uniform distribution. + * + *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ +#pragma once +//--------------------------------------------------------------------------------------------------------------------------- +// INCLUDE DESCRIPTION +//--------------------------------------------------------------------------------------------------------------------------- +#include // panic::uint_t, panic::int_t, and panic::real_t +#include +#include + + +//--------------------------------------------------------------------------------------------------------------------------- +// DEFINE DESCRIPTION +//--------------------------------------------------------------------------------------------------------------------------- + +//--------------------------------------------------------------------------------------------------------------------------- +// TYPE DESCRIPTION +//--------------------------------------------------------------------------------------------------------------------------- + +//-------------------------------------------------------------------------------------------------------------------------- +// Function Name : panic::random::seed +// +// Description: +// base for random libary +//-------------------------------------------------------------------------------------------------------------------------- +namespace panic{ + namespace random{ + +/** + * @brief Returns a value from a uniform distribution bewteen 0 and 1. + * + * Computes: + * @code + * c = uniform(); + * @endcode + * + * @return a value from a uniform distribution between 0 and 1 + * + * @note This is not a omp-safe function + */ +panic::types::real_t uniform(); + +/** + * @brief Returns a value from a uniform distribution with limits. + * + * Computes: + * @code + * c = uniform(0.1f, 42.0f); + * @endcode + * + * @tparam T Numeric element type. + * @param min minimum limit for return value. + * @param max maximum limit for return value. + * + * @return a value from a uniform distribution with limits + * + * @note This is not a omp-safe function + */ +template +T uniform(const T min, const T max); + + +/** + * @brief Filleds a vector with values from a uniform distribution from 0 to 1. + * + * Computes: + * @code + * panic::tensor::vector a(5); + * uniform(a); + * @endcode + * + * @tparam T Numeric element type. + * @param a Output vector filled iwht new values + * + * @return True of vector is filled correctly + * + * @note This is not a omp-safe function + */ +bool uniform(panic::tensor::real_vector& a); + +/** + * @brief Filleds a vector with values from a uniform distribution from min to max. + * + * Computes: + * @code + * panic::tensor::real_vector a(5); + * uniform(a, 1.0f, 42.0f); + * @endcode + * + * @tparam T Numeric element type. + * @param a Output vector filled iwht new values + * @param min minimum limit for return value. + * @param max maximum limit for return value. + * + * @return True of vector is filled correctly + * + * @note This is not a omp-safe function + */ +template +bool uniform(panic::tensor::vector& a, const T min, const T max); + + + +/** + * @brief Filleds a matrix with values from a uniform distribution from 0 to 1. + * + * Computes: + * @code + * panic::tensor::matrix A(5); + * uniform(A); + * @endcode + * + * @tparam T Numeric element type. + * @param a Output matrix filled with new values + * + * @return True of matrix is filled correctly + * + * @note This is not a omp-safe function + */ +bool uniform(panic::tensor::real_matrix& A); + + +/** + * @brief Filleds a matrix with values from a uniform distribution from min to max. + * + * Computes: + * @code + * panic::tensor::real_matrix A(5,5); + * uniform(A, 1.0f, 42.0f); + * @endcode + * + * @tparam T Numeric element type. + * @param a Output matrix filled with new values + * @param min minimum limit for return value. + * @param max maximum limit for return value. + * + * @return True of matrix is filled correctly + * + * @note This is not a omp-safe function + */ +template +bool uniform(panic::tensor::matrix& A, const T min, const T max); +/* + +panic::tensor::real_vector uniform_vector(panic::tensor::real_vector& a); +panic::tensor::real_vector uniform_vector(const panic::types::real_t min, const panic::types::real_t max); + +// Create and return vector + + + +// Create and return matrix + +*/ + + } // namespace random +} // namespace panic + +//--------------------------------------------------------------------------------------------------------------------------- +// VARIABLE DESCRIPTION +//--------------------------------------------------------------------------------------------------------------------------- + +//--------------------------------------------------------------------------------------------------------------------------- +// FUNCTION PROTOTYPE +//--------------------------------------------------------------------------------------------------------------------------- \ No newline at end of file diff --git a/main.cpp b/main.cpp index e992c04..71c36a6 100644 --- a/main.cpp +++ b/main.cpp @@ -34,6 +34,7 @@ //--------------------------------------------------------------------------------------------------------------------------- // INCLUDE DESCRIPTION //--------------------------------------------------------------------------------------------------------------------------- +#include #include // std::cout, std::endl #include // include types to use #include // Math constants @@ -43,6 +44,21 @@ #include #include #include +#include + + + +// For omp tesing: +#include +#include +#include +#include +#include +#include + +#if PANIC_HAS_OPENMP + #include +#endif //--------------------------------------------------------------------------------------------------------------------------- // DEFINE DESCRIPTION @@ -51,6 +67,601 @@ // #define TEST_FALG 1 + + + + + + + + + +//--------------------------------------------------------------------------------------------------------------------------- +// Function Name : benchmark_omp_min_work +// +// Description: +// Measures an operation with one thread and with all available threads. +// +// It then tests several possible omp_min_work values and estimates which +// threshold best matches the measured results. +//--------------------------------------------------------------------------------------------------------------------------- + +void benchmark_omp_min_work(){ + +#if !PANIC_HAS_OPENMP + + std::cout << "OpenMP is not enabled.\n"; + +#else + + //----------------------------------------------------------------------------------------------------------------------- + // OPENMP SETTINGS + //----------------------------------------------------------------------------------------------------------------------- + + omp_set_dynamic(0); + + const int parallel_thread_count = omp_get_max_threads(); + + + //----------------------------------------------------------------------------------------------------------------------- + // EDIT 1: + // Sizes to test. + // + // This list is suitable for vector operations. + //----------------------------------------------------------------------------------------------------------------------- + + const panic::types::uint_t sizes[] = { + 100, + 250, + 500, + 1000, + 2000, + 4000, + 8000, + 12000, + 16000, + 24000, + 32000, + 50000, + 75000, + 100000, + 150000, + 250000, + 500000, + 1000000, + 2000000, + 5000000 + }; + + /* + // Example sizes for square matrix multiplication: + + const panic::types::uint_t sizes[] = { + 4, + 6, + 8, + 10, + 12, + 16, + 20, + 24, + 32, + 48, + 64, + 96, + 128, + 192, + 256, + 384, + 512 + }; + */ + + const std::size_t number_of_sizes = + sizeof(sizes) / sizeof(sizes[0]); + + + //----------------------------------------------------------------------------------------------------------------------- + // EDIT 2: + // Possible omp_min_work values to test. + // + // These are work values, not necessarily vector or matrix sizes. + //----------------------------------------------------------------------------------------------------------------------- + + const std::uint64_t omp_min_work_values[] = { + 0, + 50, + 100, + 150, + 200, + 250, + 300, + 350, + 400, + 450, + 500, + 600, + 750, + 1000, + 1500, + 2000, + 4000, + 8000, + 16000, + 32000 + }; + + /* + // Example omp_min_work values for matrix multiplication: + + const std::uint64_t omp_min_work_values[] = { + 0, + 1000, + 5000, + 10000, + 25000, + 50000, + 100000, + 250000, + 500000, + 1000000, + 2000000, + 5000000, + 10000000, + 25000000, + 50000000, + 100000000 + }; + */ + + const std::size_t number_of_omp_min_work_values = + sizeof(omp_min_work_values) / + sizeof(omp_min_work_values[0]); + + + //----------------------------------------------------------------------------------------------------------------------- + // BENCHMARK SETTINGS + //----------------------------------------------------------------------------------------------------------------------- + + // The benchmark tries to perform approximately this much work for each size. + const std::uint64_t target_total_work = 50000000; + + // Minimum and maximum number of repetitions for each size. + const std::uint64_t minimum_repetitions = 3; + const std::uint64_t maximum_repetitions = 1000; + + // Prevent the compiler from treating all calculated results as unused. + volatile panic::types::real_t benchmark_sink = 0.0f; + + + //----------------------------------------------------------------------------------------------------------------------- + // ARRAYS FOR THE MEASURED RESULTS + //----------------------------------------------------------------------------------------------------------------------- + + std::uint64_t measured_work[number_of_sizes]; + + double one_thread_results[number_of_sizes]; + double parallel_results[number_of_sizes]; + + + //----------------------------------------------------------------------------------------------------------------------- + // PRINT BENCHMARK INFORMATION + //----------------------------------------------------------------------------------------------------------------------- + + std::cout + << std::right + << std::setw(12) << "Size" + << std::setw(16) << "Work" + << std::setw(14) << "Repetitions" + << std::setw(18) << "1 thread (us)" + << std::setw(18) << "Parallel (us)" + << std::setw(12) << "Speedup" + << std::setw(14) << "Fastest" + << "\n"; + + std::cout + << std::string(104, '-') + << "\n"; + + + //----------------------------------------------------------------------------------------------------------------------- + // TEST EVERY SIZE + //----------------------------------------------------------------------------------------------------------------------- + + for (std::size_t size_index = 0; + size_index < number_of_sizes; + ++size_index){ + + const panic::types::uint_t size = + sizes[size_index]; + + + //------------------------------------------------------------------------------------------------------------------- + // EDIT 3: + // Initialize the vectors or matrices used by the operation. + //------------------------------------------------------------------------------------------------------------------- + + // Vector-add example: + + panic::tensor::real_vector a(size, 1.0f); + panic::tensor::real_vector b(size, 2.0f); + panic::tensor::real_vector c(size); + + + /* + // Square-matrix multiplication example: + + panic::tensor::real_matrix a( + size, + size, + 0.01f + ); + + panic::tensor::real_matrix b( + size, + size, + 0.02f + ); + + panic::tensor::real_matrix c( + size, + size + ); + */ + + + //------------------------------------------------------------------------------------------------------------------- + // EDIT 4: + // Calculate work in the same way as the function being tested. + //------------------------------------------------------------------------------------------------------------------- + + // Vector operation: + + const std::uint64_t work = + static_cast(size); + + + /* + // Square matrix add: + + const std::uint64_t work = + static_cast(size) * + static_cast(size); + */ + + + /* + // Square matrix multiplication: + + const std::uint64_t work = + static_cast(size) * + static_cast(size) * + static_cast(size); + */ + + + //------------------------------------------------------------------------------------------------------------------- + // CALCULATE NUMBER OF REPETITIONS + //------------------------------------------------------------------------------------------------------------------- + + std::uint64_t repetitions = + target_total_work / work; + + if (repetitions < minimum_repetitions){ + repetitions = minimum_repetitions; + } + + if (repetitions > maximum_repetitions){ + repetitions = maximum_repetitions; + } + + + //------------------------------------------------------------------------------------------------------------------- + // MEASURE WITH ONE THREAD AND ALL THREADS + // + // test == 0: one OpenMP thread + // test == 1: all available OpenMP threads + //------------------------------------------------------------------------------------------------------------------- + + double measured_time_us[2] = { + 0.0, + 0.0 + }; + + for (int test = 0; test < 2; ++test){ + + if (test == 0){ + omp_set_num_threads(1); + } + else{ + omp_set_num_threads( + parallel_thread_count + ); + } + + + //---------------------------------------------------------------------------------------------------------------- + // WARM-UP + //---------------------------------------------------------------------------------------------------------------- + + for (std::size_t warmup = 0; + warmup < 2; + ++warmup){ + + //------------------------------------------------------------------------------------------------------------ + // EDIT 5: + // Put the operation being tested here. + //------------------------------------------------------------------------------------------------------------ + + panic::math::add(a, b, c); + + // Matrix multiplication: + // panic::math::matmul(a, b, c); + } + + + //---------------------------------------------------------------------------------------------------------------- + // TIMED LOOP + //---------------------------------------------------------------------------------------------------------------- + + const std::chrono::steady_clock::time_point start = + std::chrono::steady_clock::now(); + + for (std::uint64_t repetition = 0; + repetition < repetitions; + ++repetition){ + + //------------------------------------------------------------------------------------------------------------ + // EDIT 6: + // Put the same operation here. + //------------------------------------------------------------------------------------------------------------ + + panic::math::add(a, b, c); + + // Matrix multiplication: + // panic::math::matmul(a, b, c); + } + + const std::chrono::steady_clock::time_point end = + std::chrono::steady_clock::now(); + + + //---------------------------------------------------------------------------------------------------------------- + // CALCULATE AVERAGE TIME PER OPERATION + //---------------------------------------------------------------------------------------------------------------- + + const double total_time_us = + std::chrono::duration( + end - start + ).count(); + + measured_time_us[test] = + total_time_us / + static_cast(repetitions); + + + //---------------------------------------------------------------------------------------------------------------- + // READ ONE RESULT + // + // Change this if the output cannot be accessed with c[size / 2]. + //---------------------------------------------------------------------------------------------------------------- + + benchmark_sink += c[size / 2]; + + /* + // For a matrix: + + benchmark_sink += c( + size / 2, + size / 2 + ); + */ + } + + + //------------------------------------------------------------------------------------------------------------------- + // STORE RESULTS + //------------------------------------------------------------------------------------------------------------------- + + const double one_thread_us = + measured_time_us[0]; + + const double parallel_us = + measured_time_us[1]; + + measured_work[size_index] = work; + one_thread_results[size_index] = one_thread_us; + parallel_results[size_index] = parallel_us; + + + //------------------------------------------------------------------------------------------------------------------- + // PRINT RESULTS FOR THIS SIZE + //------------------------------------------------------------------------------------------------------------------- + + const double speedup = + one_thread_us / parallel_us; + + const char* fastest; + + if (parallel_us < one_thread_us){ + fastest = "parallel"; + } + else{ + fastest = "one_thread"; + } + + std::cout + << std::right + << std::setw(12) << size + << std::setw(16) << work + << std::setw(14) << repetitions + << std::setw(18) << std::fixed << std::setprecision(3) + << one_thread_us + << std::setw(18) << parallel_us + << std::setw(12) << speedup + << std::setw(14) << fastest + << "\n"; + } + + + //----------------------------------------------------------------------------------------------------------------------- + // TEST THE POSSIBLE OMP_MIN_WORK VALUES + //----------------------------------------------------------------------------------------------------------------------- + + std::uint64_t best_omp_min_work = 0; + + double best_score = + std::numeric_limits::max(); + + std::cout + << "\n" + << "omp_min_work" + << ", average_slowdown" + << "\n"; + + for (std::size_t threshold_index = 0; + threshold_index < number_of_omp_min_work_values; + ++threshold_index){ + + const std::uint64_t omp_min_work = + omp_min_work_values[threshold_index]; + + double score = 0.0; + + + //------------------------------------------------------------------------------------------------------------------- + // SEE WHICH VERSION THIS THRESHOLD WOULD SELECT + //------------------------------------------------------------------------------------------------------------------- + + for (std::size_t size_index = 0; + size_index < number_of_sizes; + ++size_index){ + + double selected_time_us; + + + //---------------------------------------------------------------------------------------------------------------- + // This uses > because the PANIC functions currently use: + // + // work > omp_min_work + //---------------------------------------------------------------------------------------------------------------- + + if (measured_work[size_index] > omp_min_work){ + + // This threshold would select OpenMP. + selected_time_us = + parallel_results[size_index]; + } + else{ + + // This threshold would select serial execution. + selected_time_us = + one_thread_results[size_index]; + } + + + //---------------------------------------------------------------------------------------------------------------- + // FIND THE FASTEST MEASURED VERSION FOR THIS SIZE + //---------------------------------------------------------------------------------------------------------------- + + double fastest_time_us = + one_thread_results[size_index]; + + if ( + parallel_results[size_index] < + fastest_time_us + ){ + fastest_time_us = + parallel_results[size_index]; + } + + + //---------------------------------------------------------------------------------------------------------------- + // CALCULATE HOW MUCH SLOWER THE SELECTED VERSION IS + // + // 1.0 means the threshold selected the fastest version. + // 1.1 means it was 10% slower than the fastest version. + //---------------------------------------------------------------------------------------------------------------- + + score += + selected_time_us / + fastest_time_us; + } + + + //------------------------------------------------------------------------------------------------------------------- + // CALCULATE THE AVERAGE SCORE + //------------------------------------------------------------------------------------------------------------------- + + score /= + static_cast(number_of_sizes); + + + //------------------------------------------------------------------------------------------------------------------- + // PRINT THIS OMP_MIN_WORK RESULT + //------------------------------------------------------------------------------------------------------------------- + + std::cout + << omp_min_work + << ", " + << std::fixed + << std::setprecision(4) + << score + << "\n"; + + + //------------------------------------------------------------------------------------------------------------------- + // SAVE THE BEST OMP_MIN_WORK + //------------------------------------------------------------------------------------------------------------------- + + if (score < best_score){ + + best_score = score; + best_omp_min_work = omp_min_work; + } + } + + + //----------------------------------------------------------------------------------------------------------------------- + // PRINT FINAL ESTIMATE + //----------------------------------------------------------------------------------------------------------------------- + + std::cout + << "\n" + << "Best estimated omp_min_work: " + << best_omp_min_work + << "\n"; + + std::cout + << "Average slowdown score: " + << std::fixed + << std::setprecision(4) + << best_score + << "\n"; + + std::cout + << "\n" + << "A score of 1.0000 means the threshold selected the\n" + << "fastest measured version for every tested size.\n"; + + + //----------------------------------------------------------------------------------------------------------------------- + // PRINT CHECKSUM + //----------------------------------------------------------------------------------------------------------------------- + + std::cout + << "Benchmark checksum: " + << benchmark_sink + << "\n"; + +#endif +} + + + + + + //--------------------------------------------------------------------------------------------------------------------------- // VARIABLE DESCRIPTION //--------------------------------------------------------------------------------------------------------------------------- @@ -73,6 +684,10 @@ panic::tensor::int_matrix C(2,2, 3); int main(void) { + //benchmark_omp_min_work(); + + // Comment out benchmark_omp_min_work() when it is not needed. + a[2] = 1.2; b[2] = 3; @@ -107,6 +722,40 @@ int main(void) { panic::neural_network::layer_dense dense(100,10); std::cout << dense.forward(B1) << std::endl; + std::cout << "random" << std::endl; + + for (int i = 0; i < 5; ++i) + { + std::cout << panic::random::uniform() << std::endl; + } + + std::cout << "random(min, max)" << std::endl; + + panic::types::real_t a1 = -15; + panic::types::real_t a2 = 15; + + for (int i = 0; i < 5; ++i) + { + std::cout << panic::random::uniform(a1,a2) << std::endl; + } + + panic::tensor::real_vector a3(2); + panic::random::uniform(a3); + panic::io::print_vector(a3); + + + panic::tensor::uint_vector a4(2); + panic::random::uniform(a4, static_cast(1), static_cast(3)); + panic::io::print_vector(a4); + + panic::tensor::real_matrix D1(3,3); + panic::random::uniform(D1); + panic::io::print_matrix(D1); + + panic::tensor::real_matrix D3(3,3); + panic::random::uniform(D3, 100.f, 200.f); + panic::io::print_matrix(D3); + return 0; diff --git a/src/io/print_tensor.cpp b/src/io/print_tensor.cpp index 6e45b83..007471f 100644 --- a/src/io/print_tensor.cpp +++ b/src/io/print_tensor.cpp @@ -99,20 +99,27 @@ template bool print_vector(const panic::tensor::vector bool print_matrix(const panic::tensor::matrix& A){ - std::cout << "["; - for (panic::types::uint_t i = 0; i < A.rows(); ++i){ - std::cout << "["; - for (panic::types::uint_t j = 0; j < A.cols(); ++j){ - std::cout << A(i,j) << ", "; - } - if (i < A.rows()-1){ - std::cout << A(A.rows()-1,A.cols()-1) << "]" << std::endl; - }else{ - std::cout << A(A.rows()-1,A.cols()-1) << "]]" << std::endl; - } - } - return true; + std::cout << "["; + for (panic::types::uint_t i = 0; i < A.rows(); ++i){ + std::cout << "["; + + for (panic::types::uint_t j = 0; j < A.cols(); ++j){ + std::cout << A(i, j); + + if (j + 1 < A.cols()){ + std::cout << ", "; + } + } + + std::cout << "]"; + + if (i + 1 < A.rows()) + std::cout << "," << std::endl; + } + + std::cout << "]" << std::endl; + return true; } //-------------------------------------------------------------------------------------------------------------------------- diff --git a/src/math/add.cpp b/src/math/add.cpp index a72fd84..ec49721 100644 --- a/src/math/add.cpp +++ b/src/math/add.cpp @@ -51,8 +51,8 @@ * Small vectors and matrices are kept serial because the overhead of starting * worker threads can be larger than the work itself. */ -static const panic::types::uint_t add_omp_min_work = 10000; - +//static const panic::types::uint_t add_omp_min_work = 10000; +static const panic::types::uint_t add_omp_min_work = 0; //--------------------------------------------------------------------------------------------------------------------------- // INPLEMENTATION //--------------------------------------------------------------------------------------------------------------------------- diff --git a/src/random/seed.cpp b/src/random/seed.cpp new file mode 100644 index 0000000..73bed7e --- /dev/null +++ b/src/random/seed.cpp @@ -0,0 +1,122 @@ +/**++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + * + * PANIC + * Portable Algorithms and Numerics In C++ + * + * Scientific computing from scratch, with feeling. + * + * Copyright (c) 2026 Michelle Bausager + * + * This file is part of PANIC. + * + * PANIC is free software licensed under the GNU General Public License v3.0 or later. + * You may redistribute and/or modify it under the terms of the GPL. + * + * PANIC is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; + * without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + * See the LICENSE file for the full license text. + * + * SPDX-License-Identifier: GPL-3.0-or-later + * + *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + * + * Project Name: PANIC + * Module Name: random + * File Name: seed.cpp + * Revision: 0.1.0 + * Date: 28-06-2026 + * Author: Michelle Bausager + * + * Description: + * Defines seed for use in other functions in random/ + * + *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ + +//--------------------------------------------------------------------------------------------------------------------------- +// INCLUDE DESCRIPTION +//--------------------------------------------------------------------------------------------------------------------------- +#include +#include + +//--------------------------------------------------------------------------------------------------------------------------- +// FUNCTION PROTOTYPE +//--------------------------------------------------------------------------------------------------------------------------- +namespace panic{ + namespace random{ + + +//-------------------------------------------------------------------------------------------------------------------------- +// Constructor Name : panic::random::seed::seed_t +// +// Description: +// Creates an seed object. +//-------------------------------------------------------------------------------------------------------------------------- +seed_t::seed_t(){ +value = 1; +} + +//-------------------------------------------------------------------------------------------------------------------------- +// Constructor Name : panic::random::seed::seed_t +// +// Description: +// Creates an seed object with a specific seed. +//-------------------------------------------------------------------------------------------------------------------------- +seed_t::seed_t(panic::types::uint_t seed){ + set(seed); +} + +//-------------------------------------------------------------------------------------------------------------------------- +// Constructor Name : panic::random::seed::set_seed +// +// Description: +// Sets the seed to a none-zero value +//-------------------------------------------------------------------------------------------------------------------------- +bool seed_t::set(panic::types::uint_t seed){ + + if (seed == 0){ + value = 1; + }else{ + value = seed; + } + return 1; +} + +//-------------------------------------------------------------------------------------------------------------------------- +// Constructor Name : panic::random::seed::get +// +// Description: +// Returns the seed. +//-------------------------------------------------------------------------------------------------------------------------- +panic::types::uint_t seed_t::get(){ + return value; +} + + +//-------------------------------------------------------------------------------------------------------------------------- +// Constructor Name : panic::random::seed::state_at +// +// Description: +// Returns a pshodo-random number based on the input index and the seed value. +//-------------------------------------------------------------------------------------------------------------------------- +panic::types::uint_t seed_t::state_at(panic::types::uint_t index) const{ + + panic::types::uint_t state = value + index; + + // Mix it a few times so close indexes do not start too similarly. + //next_state(state); + //next_state(state); + //next_state(state); + + state ^= state >> 16; + state *= static_cast(0x7feb352d); + state ^= state >> 15; + state *= static_cast(0x846ca68b); + state ^= state >> 16; + + return state; + +} + + + } // namespace tensor +} // namespace panic \ No newline at end of file diff --git a/src/random/uniform.cpp b/src/random/uniform.cpp new file mode 100644 index 0000000..1a35325 --- /dev/null +++ b/src/random/uniform.cpp @@ -0,0 +1,260 @@ +/**++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + * + * PANIC + * Portable Algorithms and Numerics In C++ + * + * Scientific computing from scratch, with feeling. + * + * Copyright (c) 2026 Michelle Bausager + * + * This file is part of PANIC. + * + * PANIC is free software licensed under the GNU General Public License v3.0 or later. + * You may redistribute and/or modify it under the terms of the GPL. + * + * PANIC is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; + * without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + * See the LICENSE file for the full license text. + * + * SPDX-License-Identifier: GPL-3.0-or-later + * + *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + * + * Project Name: PANIC + * Module Name: random + * File Name: uniform.cpp + * Revision: 0.1.0 + * Date: 29-06-2026 + * Author: Michelle Bausager + * + * Description: + * Defines the dense layers used in neural network + * + *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ +//--------------------------------------------------------------------------------------------------------------------------- +// INCLUDE DESCRIPTION +//--------------------------------------------------------------------------------------------------------------------------- +#include +#include +#include +#include + +#include + + + + +//--------------------------------------------------------------------------------------------------------------------------- +// VARIABLE DESCRIPTION +//--------------------------------------------------------------------------------------------------------------------------- +static const panic::types::uint_t uniform_omp_min_work = 5000; +static const panic::types::uint_t max_uint_t = ~static_cast(0); +static const panic::types::real_t real_max_unit = static_cast(max_uint_t); + +static panic::types::uint_t uniform_state = static_cast(1); + + +//--------------------------------------------------------------------------------------------------------------------------- +// FUNCTION PROTOTYPE +//--------------------------------------------------------------------------------------------------------------------------- +namespace panic{ + namespace random{ + + +//-------------------------------------------------------------------------------------------------------------------------- +// Function Name : panic::random::uniform +// +// Description: +// returns a value from a uniform distribution bewteen 0 and 1 +//-------------------------------------------------------------------------------------------------------------------------- +panic::types::real_t uniform(){ + panic::random::seed_t seed; + uniform_state += static_cast(1); + seed.set(uniform_state); + + return (static_cast(seed.state_at(seed.get()) / real_max_unit)); + +} + + + +//-------------------------------------------------------------------------------------------------------------------------- +// Function Name : panic::random::uniform +// +// Description: +// returns a value from a uniform distribution bewteen min and max +//-------------------------------------------------------------------------------------------------------------------------- +template +T uniform(const T min, const T max){ + return (min + static_cast(static_cast((max - min)) * uniform())); +} +//-------------------------------------------------------------------------------------------------------------------------- +// EXPLICIT TEMPLATE INSTANTIATION +// +// The implementation is in this .cpp file. +// Build the overload for the official PANIC numeric types. +//-------------------------------------------------------------------------------------------------------------------------- +template panic::types::uint_t + uniform(const panic::types::uint_t min, + const panic::types::uint_t mix +); +template panic::types::int_t + uniform(const panic::types::int_t min, + const panic::types::int_t max +); +template panic::types::real_t + uniform(const panic::types::real_t min, + const panic::types::real_t max +); + + + +//-------------------------------------------------------------------------------------------------------------------------- +// Function Name : panic::random::uniform +// +// Description: +// Fills a vector with a uniform distribution +//-------------------------------------------------------------------------------------------------------------------------- +bool uniform(panic::tensor::real_vector& a){ + + panic::random::seed_t seed; + panic::types::uint_t work = a.size(); + uniform_state += work; + seed.set(uniform_state); + + PANIC_OMP_PARALLEL_FOR_IF(work > uniform_omp_min_work) + for (panic::types::uint_t i = 0; i < work; ++i){ + a[i] = (static_cast(seed.state_at(i)) / real_max_unit); + } + + return true; + +} + + +//-------------------------------------------------------------------------------------------------------------------------- +// Function Name : panic::random::uniform +// +// Description: +// Fills a vector with a uniform distribution with limits +//-------------------------------------------------------------------------------------------------------------------------- +template +bool uniform(panic::tensor::vector& a, const T min, const T max){ + + panic::random::seed_t seed; + panic::types::uint_t work = a.size(); + + uniform_state += work; + seed.set(uniform_state); + + panic::types::real_t temp = static_cast(max - min); + + PANIC_OMP_PARALLEL_FOR_IF(work > uniform_omp_min_work) + for (panic::types::uint_t i = 0; i < work; ++i){ + a[i] = min + static_cast((temp*(static_cast(seed.state_at(i)) / real_max_unit))); + } + + return true; +} +//-------------------------------------------------------------------------------------------------------------------------- +// EXPLICIT TEMPLATE INSTANTIATION +// +// The implementation is in this .cpp file. +// Build the overload for the official PANIC numeric types. +//-------------------------------------------------------------------------------------------------------------------------- +template bool uniform(panic::tensor::vector& a, + const panic::types::uint_t min, + const panic::types::uint_t mix +); +template bool uniform(panic::tensor::vector& a, + const panic::types::int_t min, + const panic::types::int_t max +); +template bool uniform(panic::tensor::vector& a, + const panic::types::real_t min, + const panic::types::real_t max +); + + +//-------------------------------------------------------------------------------------------------------------------------- +// Function Name : panic::random::uniform +// +// Description: +// Fills a matrix with a uniform distribution +//-------------------------------------------------------------------------------------------------------------------------- +bool uniform(panic::tensor::real_matrix& A){ + + panic::random::seed_t seed; + + panic::types::uint_t rows = A.rows(); + panic::types::uint_t cols = A.cols(); + + panic::types::uint_t work = rows*cols; + + uniform_state += work; + seed.set(uniform_state); + + PANIC_OMP_PARALLEL_FOR_IF(work > uniform_omp_min_work) + for (panic::types::uint_t i = 0; i < rows; ++i){ + for (panic::types::uint_t j = 0; j < cols; ++j){ + panic::types::uint_t index = i * cols + j; + A(i,j) = (static_cast(seed.state_at(index)) / real_max_unit); + } + } + + + return true; +} + +//-------------------------------------------------------------------------------------------------------------------------- +// Function Name : panic::random::uniform +// +// Description: +// Fills a matrix with a uniform distribution with limits +//-------------------------------------------------------------------------------------------------------------------------- +template +bool uniform(panic::tensor::matrix& A, const T min, const T max){ + + panic::random::seed_t seed; + + panic::types::uint_t rows = A.rows(); + panic::types::uint_t cols = A.cols(); + + panic::types::uint_t work = rows*cols; + + uniform_state += work; + seed.set(uniform_state); + + panic::types::real_t temp = static_cast(max - min); + + PANIC_OMP_PARALLEL_FOR_IF(work > uniform_omp_min_work) + for (panic::types::uint_t i = 0; i < rows; ++i){ + for (panic::types::uint_t j = 0; j < cols; ++j){ + panic::types::uint_t index = i * cols + j; + A(i,j) = min + static_cast((temp*(static_cast(seed.state_at(index)) / real_max_unit))); + } + } + return true; +} +//-------------------------------------------------------------------------------------------------------------------------- +// EXPLICIT TEMPLATE INSTANTIATION +// +// The implementation is in this .cpp file. +// Build the overload for the official PANIC numeric types. +//-------------------------------------------------------------------------------------------------------------------------- +template bool uniform(panic::tensor::matrix& a, + const panic::types::uint_t min, + const panic::types::uint_t mix +); +template bool uniform(panic::tensor::matrix& a, + const panic::types::int_t min, + const panic::types::int_t max +); +template bool uniform(panic::tensor::matrix& a, + const panic::types::real_t min, + const panic::types::real_t max +); + + + } // namespace random +} // namespace panic \ No newline at end of file