random_array_generator.cpp

Go to the documentation of this file.

/*
 *  SPDX-FileCopyrightText: Copyright 2021, Siavash Ameli <sameli@berkeley.edu>
 *  SPDX-License-Identifier: BSD-3-Clause
 *  SPDX-FileType: SOURCE
 *
 *  This program is free software: you can redistribute it and/or modify it
 *  under the terms of the license found in the LICENSE.txt file in the root
 *  directory of this source tree.
 */
 
 
// =======
// Headers
// =======
 
#include "./random_array_generator.h"
#include <omp.h>  // omp_set_num_threads
#include <stdint.h>  // uint64_t
#include "../_c_trace_estimator/c_orthogonalization.h"  // cOrthogonalization
#include "../_c_basic_algebra/c_vector_operations.h"  // cVectorOperations
 
 
// =====================
// generate random array
// =====================
 
 
template <typename DataType>
void RandomArrayGenerator<DataType>::generate_random_array(
        RandomNumberGenerator& random_number_generator,
        DataType* array,
        const LongIndexType array_size,
        const IndexType num_threads)
{
    // Set the number of threads only if num_threads is non-zero. If
    // num_threads is zero, it indicates to not create new threads in this
    // function, rather, this function was called from another caller function
    // that was already parallelized and this function should be executed in
    // one of those threads.
    if (num_threads > 0)
    {
        // num_threads parallel threads will be created in this function.
        omp_set_num_threads(num_threads);
    }
 
    // Finding the thread id. It depends where we call omp_get_thread_num. If
    // we call it here (below), it gives the thread id of the parent function.
    // But if we call it inside a #pragma omp parallel, it gives the newer
    // thread id that is "created" by the pragma.
    int thread_id = 0;
    if (num_threads == 0)
    {
        // If num_threads is zero (which means we will not create a new thread
        // in this function), we get the thread id of the parent function.
        thread_id = omp_get_thread_num();
    }
 
    // Number of bits to generate in each call of the random generator. This is
    // the number of bits in a uint64_t integer.
    const int bits_per_byte = 8;
    const int num_bits = sizeof(uint64_t) * bits_per_byte;
 
    // Number of chunks of divided array each in num_bits length
    const LongIndexType num_chunks = \
            static_cast<LongIndexType>(array_size/num_bits);
 
    // Shared-memory parallelism over individual row vectors. The parallel
    // section only works if num_threads is non-zero. Otherwise it runs in
    // serial order.
    #pragma omp parallel if (num_threads > 0)
    {
        // If num_threads is zero, the following thread_id is the thread id
        // that the parent (caller) function created outside of this function.
        // But, if num_thread is non-zero, the following thread id is the
        // thread id that is created inside this parallel loop.
        if (num_threads > 0)
        {
            thread_id = omp_get_thread_num();
        }
 
        #pragma omp for schedule(static)
        for (LongIndexType i=0; i < num_chunks; ++i)
        {
            // Generate 64 bits (one integer)
            uint64_t bits = random_number_generator.next(thread_id);
 
            // Fill 64 elements of array with +1 or -1 depending on the bits
            for (IndexType j=0; j < num_bits; ++j)
            {
                // Check if the j-th bit (from right to left) is 1 or 0
                if (bits & ( uint64_t(1) << j))
                {
                    // Bit is 1. Write +1.0 in array
                    array[i*num_bits + j] = 1.0;
                }
                else
                {
                    // Bit is 0. Write -1.0 in array
                    array[i*num_bits + j] = -1.0;
                }
            }
        }
    }
 
    // The previous for loop (the above) does not fill all elements of array
    // since it only iterates on the multiples of 64 (num_bits). We fill the
    // rest of the array. There are less than 64 elements remained to fill. So,
    // it suffice to generate only 64 bits (one integer) for the rest of the
    // elements of the array
    uint64_t bits = random_number_generator.next(thread_id);
    IndexType j = 0;
 
    // This loop should have less than 64 iterations.
    for (LongIndexType i = num_chunks * num_bits; i < array_size; ++i)
    {
        // Check if the j-th bit (from right to left) is 1 or 0
        if (bits & ( uint64_t(1) << j))
        {
            // Bit is 1. Write +1.0 in array
            array[i] = 1.0;
        }
        else
        {
            // Bit is 0. Write -1.0 in array
            array[i] = -1.0;
        }
 
        // Increment bit counter
        ++j;
    }
}
 
 
// ===============================
// Explicit template instantiation
// ===============================
 
template class RandomArrayGenerator<float>;
template class RandomArrayGenerator<double>;
template class RandomArrayGenerator<long double>;