A static class to generate random set of vectors. This class acts as a templated namespace, where all member methods are public and static. More...

#include <random_array_generator.h>

Static Public Member Functions
static void	generate_random_array (RandomNumberGenerator &random_number_generator, DataType *array, const LongIndexType array_size, const IndexType num_threads)
	Generates a pseudo-random array with Rademacher distribution where elements are either `+1` or `-1`.

Detailed Description

template<typename DataType>
class RandomArrayGenerator< DataType >

A static class to generate random set of vectors. This class acts as a templated namespace, where all member methods are public and static.

See also: Orthogonalization

Definition at line 37 of file random_array_generator.h.

Member Function Documentation

◆ 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
	)

static

Generates a pseudo-random array with Rademacher distribution where elements are either +1 or -1.

The Rademacher distribution is obtained from the Bernoulli distribution consisting of 0 and 1. To generate such distribution, a sequence of array_size/64 intergers, each with 64-bit, is generated using Xoshiro256** algorithm. The 64 bits of each integer are used to fill 64 elements of the array as follows. If the bit is 0, the array element is set to -1, and if the bit is 1, the array element is set to +1.

Thus, in this function, we use Xoshiro256** algorithm to generate 64 bits and use bits, not the integer itself. This approach is about ten times faster than convertng the random integer to double between [0,1] and then map them to +1 and -1.

Also, this function is more than a hundered times faster than using rand() function.

Parameters

[in]	random_number_generator	The random number generator object. This object should be initialized with `num_threads` by its constructor. On each parallel thread, an independent sequence of random numbers are generated.
[out]	array	1D array of the size `array_size`.
[out]	array_size	The size of the array.
[in]	num_threads	Number of OpenMP parallel threads. If `num_threads` is zero then no paralel thread is created inside this function, rather it is assumed that this functon is called inside a parallel region from the caller.

Definition at line 66 of file random_array_generator.cpp.

{
    // 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.
        #if defined(USE_OPENMP) && (USE_OPENMP == 1)
            omp_set_num_threads(num_threads);
        #endif
    }
 
    // 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.
        #if defined(USE_OPENMP) && (USE_OPENMP == 1)
            thread_id = omp_get_thread_num();
        #else
            thread_id = 0;
        #endif
    }
 
    // 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.
    #if defined(USE_OPENMP) && (USE_OPENMP == 1)
    #pragma omp parallel if (num_threads > 0)
    #endif
    {
        // 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)
        {
            #if defined(USE_OPENMP) && (USE_OPENMP == 1)
                thread_id = omp_get_thread_num();
            #else
                thread_id = 0;
            #endif
        }
 
        #if defined(USE_OPENMP) && (USE_OPENMP == 1)
        #pragma omp for schedule(static)
        #endif
        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;
    }
}

References RandomNumberGenerator::next(), omp_get_thread_num(), and omp_set_num_threads().

Referenced by cTraceEstimator< DataType >::_c_stochastic_lanczos_quadrature(), cuTraceEstimator< DataType >::_cu_stochastic_lanczos_quadrature(), cuOrthogonalization< DataType >::orthogonalize_vectors(), and cOrthogonalization< DataType >::orthogonalize_vectors().

Here is the call graph for this function:

Here is the caller graph for this function:

The documentation for this class was generated from the following files:

/home/runner/work/imate/imate/imate/_random_generator/random_array_generator.h
/home/runner/work/imate/imate/imate/_random_generator/random_array_generator.cpp

Static Public Member Functions

Detailed Description

Member Function Documentation

◆ generate_random_array()