Skip to content

Commit

Permalink
pseudo-random number generation
Browse files Browse the repository at this point in the history
  • Loading branch information
cmendl committed Jul 11, 2024
1 parent 5f54ec5 commit 11b9527
Show file tree
Hide file tree
Showing 10 changed files with 494 additions and 61 deletions.
2 changes: 1 addition & 1 deletion CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@ find_package(LAPACK REQUIRED)
find_package(HDF5 REQUIRED COMPONENTS C)

set(CHEMTENSOR_DIRS "src" "src/tensor" "src/mps" "src/operator" "src/algorithm" "src/util")
set(CHEMTENSOR_SOURCES "src/tensor/dense_tensor.c" "src/tensor/block_sparse_tensor.c" "src/tensor/qnumber.c" "src/tensor/clebsch_gordan.c" "src/tensor/su2_recoupling.c" "src/tensor/su2_tree.c" "src/tensor/su2_tensor.c" "src/mps/mps.c" "src/operator/op_chain.c" "src/operator/local_op.c" "src/operator/mpo_graph.c" "src/operator/mpo.c" "src/operator/ttno_graph.c" "src/operator/ttno.c" "src/operator/hamiltonian.c" "src/algorithm/bond_ops.c" "src/algorithm/operation.c" "src/algorithm/dmrg.c" "src/algorithm/gradient.c" "src/util/util.c" "src/util/queue.c" "src/util/linked_list.c" "src/util/hash_table.c" "src/util/abstract_graph.c" "src/util/bipartite_graph.c" "src/util/integer_linear_algebra.c" "src/util/krylov.c")
set(CHEMTENSOR_SOURCES "src/tensor/dense_tensor.c" "src/tensor/block_sparse_tensor.c" "src/tensor/qnumber.c" "src/tensor/clebsch_gordan.c" "src/tensor/su2_recoupling.c" "src/tensor/su2_tree.c" "src/tensor/su2_tensor.c" "src/mps/mps.c" "src/operator/op_chain.c" "src/operator/local_op.c" "src/operator/mpo_graph.c" "src/operator/mpo.c" "src/operator/ttno_graph.c" "src/operator/ttno.c" "src/operator/hamiltonian.c" "src/algorithm/bond_ops.c" "src/algorithm/operation.c" "src/algorithm/dmrg.c" "src/algorithm/gradient.c" "src/util/util.c" "src/util/queue.c" "src/util/linked_list.c" "src/util/hash_table.c" "src/util/abstract_graph.c" "src/util/bipartite_graph.c" "src/util/integer_linear_algebra.c" "src/util/krylov.c" "src/util/pcg_basic.c" "src/util/rng.c")
set(TEST_SOURCES "test/tensor/test_dense_tensor.c" "test/tensor/test_block_sparse_tensor.c" "test/tensor/test_clebsch_gordan.c" "test/tensor/test_su2_tree.c" "test/tensor/test_su2_tensor.c" "test/mps/test_mps.c" "test/operator/test_mpo_graph.c" "test/operator/test_mpo.c" "test/operator/test_ttno_graph.c" "test/operator/test_ttno.c" "test/operator/test_hamiltonian.c" "test/algorithm/test_bond_ops.c" "test/algorithm/test_operation.c" "test/algorithm/test_dmrg.c" "test/algorithm/numerical_gradient.c" "test/algorithm/test_gradient.c" "test/util/test_queue.c" "test/util/test_linked_list.c" "test/util/test_hash_table.c" "test/util/test_bipartite_graph.c" "test/util/test_integer_linear_algebra.c" "test/util/test_krylov.c" "test/run_tests.c")

add_executable(chemtensor_test ${CHEMTENSOR_SOURCES} ${TEST_SOURCES})
Expand Down
17 changes: 17 additions & 0 deletions src/tensor/block_sparse_tensor.c
Original file line number Diff line number Diff line change
Expand Up @@ -445,6 +445,23 @@ void conjugate_block_sparse_tensor(struct block_sparse_tensor* t)
}


//________________________________________________________________________________________________________________________
///
/// \brief Fill the non-zero blocks of a block-sparse tensor with random normal entries.
///
void block_sparse_tensor_fill_random_normal(const void* alpha, const void* shift, struct rng_state* rng_state, struct block_sparse_tensor* t)
{
const long nblocks = integer_product(t->dim_blocks, t->ndim);
for (long k = 0; k < nblocks; k++)
{
struct dense_tensor* b = t->blocks[k];
if (b != NULL) {
dense_tensor_fill_random_normal(alpha, shift, rng_state, b);
}
}
}


//________________________________________________________________________________________________________________________
///
/// \brief Convert a block-sparse to an equivalent dense tensor.
Expand Down
2 changes: 2 additions & 0 deletions src/tensor/block_sparse_tensor.h
Original file line number Diff line number Diff line change
Expand Up @@ -74,6 +74,8 @@ void rscale_block_sparse_tensor(const void* alpha, struct block_sparse_tensor* t

void conjugate_block_sparse_tensor(struct block_sparse_tensor* t);

void block_sparse_tensor_fill_random_normal(const void* alpha, const void* shift, struct rng_state* rng_state, struct block_sparse_tensor* t);


//________________________________________________________________________________________________________________________
//
Expand Down
59 changes: 59 additions & 0 deletions src/tensor/dense_tensor.c
Original file line number Diff line number Diff line change
Expand Up @@ -517,6 +517,65 @@ void dense_tensor_set_identity(struct dense_tensor* t)
}


//________________________________________________________________________________________________________________________
///
/// \brief Fill the entries of a dense tensor with random normal entries.
///
void dense_tensor_fill_random_normal(const void* alpha, const void* shift, struct rng_state* rng_state, struct dense_tensor* t)
{
const long nelem = dense_tensor_num_elements(t);

switch (t->dtype)
{
case SINGLE_REAL:
{
float* data = t->data;
const float a = *((float*)alpha);
const float s = *((float*)shift);
for (long j = 0; j < nelem; j++) {
data[j] = a * randnf(rng_state) + s;
}
break;
}
case DOUBLE_REAL:
{
double* data = t->data;
const double a = *((double*)alpha);
const double s = *((double*)shift);
for (long j = 0; j < nelem; j++) {
data[j] = a * randn(rng_state) + s;
}
break;
}
case SINGLE_COMPLEX:
{
scomplex* data = t->data;
const scomplex a = *((scomplex*)alpha);
const scomplex s = *((scomplex*)shift);
for (long j = 0; j < nelem; j++) {
data[j] = a * crandnf(rng_state) + s;
}
break;
}
case DOUBLE_COMPLEX:
{
dcomplex* data = t->data;
const dcomplex a = *((dcomplex*)alpha);
const dcomplex s = *((dcomplex*)shift);
for (long j = 0; j < nelem; j++) {
data[j] = a * crandn(rng_state) + s;
}
break;
}
default:
{
// unknown data type
assert(false);
}
}
}


//________________________________________________________________________________________________________________________
///
/// \brief Temporary data structure storing tensor dimensions and a permutation of them.
Expand Down
3 changes: 3 additions & 0 deletions src/tensor/dense_tensor.h
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@
#include <assert.h>
#include "numeric.h"
#include "util.h"
#include "rng.h"


//________________________________________________________________________________________________________________________
Expand Down Expand Up @@ -117,6 +118,8 @@ void conjugate_dense_tensor(struct dense_tensor* t);

void dense_tensor_set_identity(struct dense_tensor* t);

void dense_tensor_fill_random_normal(const void* alpha, const void* shift, struct rng_state* rng_state, struct dense_tensor* t);


//________________________________________________________________________________________________________________________
//
Expand Down
151 changes: 151 additions & 0 deletions src/util/pcg_basic.c
Original file line number Diff line number Diff line change
@@ -0,0 +1,151 @@
/// \file pcg_basic.c
/// \brief PCG random number generation.

// PCG Random Number Generation for C.
//
// Copyright 2014 Melissa O'Neill <[email protected]>
//
// 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.
//
// For additional information about the PCG random number generation scheme,
// including its license and other licensing options, visit
//
// http://www.pcg-random.org


#include "pcg_basic.h"


//________________________________________________________________________________________________________________________
///
/// \brief Seed the 32-bit RNG. Specified in two parts, state initializer and a sequence selection constant (a.k.a. stream id).
///
void pcg32_srandom_r(const uint64_t initstate, const uint64_t initseq, struct pcg32_random* rng)
{
rng->state = 0U;
rng->inc = (initseq << 1u) | 1u;
pcg32_random_r(rng);
rng->state += initstate;
pcg32_random_r(rng);
}


//________________________________________________________________________________________________________________________
///
/// \brief Generate a uniformly distributed 32-bit random number.
///
uint32_t pcg32_random_r(struct pcg32_random* rng)
{
uint64_t oldstate = rng->state;
rng->state = oldstate * 6364136223846793005ULL + rng->inc;
uint32_t xorshifted = ((oldstate >> 18u) ^ oldstate) >> 27u;
uint32_t rot = oldstate >> 59u;
return (xorshifted >> rot) | (xorshifted << ((-rot) & 31));
}


//________________________________________________________________________________________________________________________
///
/// \brief Generate a uniformly distributed number, r, where 0 <= r < bound.
///
uint32_t pcg32_boundedrand_r(const uint32_t bound, struct pcg32_random* rng)
{
// To avoid bias, we need to make the range of the RNG a multiple of
// bound, which we do by dropping output less than a threshold.
// A naive scheme to calculate the threshold would be to do
//
// uint32_t threshold = 0x100000000ull % bound;
//
// but 64-bit div/mod is slower than 32-bit div/mod (especially on
// 32-bit platforms). In essence, we do
//
// uint32_t threshold = (0x100000000ull-bound) % bound;
//
// because this version will calculate the same modulus, but the LHS
// value is less than 2^32.

uint32_t threshold = -bound % bound;

// Uniformity guarantees that this loop will terminate. In practice, it
// should usually terminate quickly; on average (assuming all bounds are
// equally likely), 82.25% of the time, we can expect it to require just
// one iteration. In the worst case, someone passes a bound of 2^31 + 1
// (i.e., 2147483649), which invalidates almost 50% of the range. In
// practice, bounds are typically small and only a tiny amount of the range
// is eliminated.
for (;;) {
uint32_t r = pcg32_random_r(rng);
if (r >= threshold) {
return r % bound;
}
}
}


// This code shows how you can cope if you're on a 32-bit platform (or a
// 64-bit platform with a mediocre compiler) that doesn't support 128-bit math,
// or if you're using the basic version of the library which only provides
// 32-bit generation.
//
// Here we build a 64-bit generator by tying together two 32-bit generators.
// Note that we can do this because we set up the generators so that each
// 32-bit generator has a *totally different* different output sequence
// -- if you tied together two identical generators, that wouldn't be nearly
// as good.
//
// For simplicity, we keep the period fixed at 2^64. The state space is
// approximately 2^254 (actually 2^64 * 2^64 * 2^63 * (2^63 - 1)), which
// is huge.


//________________________________________________________________________________________________________________________
///
/// \brief Seed the 64-bit RNG.
///
void pcg32x2_srandom_r(uint64_t seed1, uint64_t seed2, uint64_t seq1, uint64_t seq2, struct pcg32x2_random* rng)
{
uint64_t mask = ~0ull >> 1;
// The stream for each of the two generators *must* be distinct
if ((seq1 & mask) == (seq2 & mask)) {
seq2 = ~seq2;
}
pcg32_srandom_r(seed1, seq1, &rng->gen[0]);
pcg32_srandom_r(seed2, seq2, &rng->gen[1]);
}


//________________________________________________________________________________________________________________________
///
/// \brief Generate a uniformly distributed 64-bit random number.
///
uint64_t pcg32x2_random_r(struct pcg32x2_random* rng)
{
return ((uint64_t)(pcg32_random_r(&rng->gen[0])) << 32)
| pcg32_random_r(&rng->gen[1]);
}


//________________________________________________________________________________________________________________________
///
/// \brief Generate a uniformly distributed 64-bit random number, r, where 0 <= r < bound.
///
uint64_t pcg32x2_boundedrand_r(const uint64_t bound, struct pcg32x2_random* rng)
{
uint64_t threshold = -bound % bound;
for (;;) {
uint64_t r = pcg32x2_random_r(rng);
if (r >= threshold) {
return r % bound;
}
}
}
63 changes: 63 additions & 0 deletions src/util/pcg_basic.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,63 @@
/// \file pcg_basic.h
/// \brief PCG random number generation.

// PCG Random Number Generation for C.
//
// Copyright 2014 Melissa O'Neill <[email protected]>
//
// 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.
//
// For additional information about the PCG random number generation scheme,
// including its license and other licensing options, visit
//
// http://www.pcg-random.org


#pragma once

#include <inttypes.h>


//________________________________________________________________________________________________________________________
///
/// \brief PCG 32-bit random number generator state.
///
struct pcg32_random
{
uint64_t state; //!< RNG state; all values are possible
uint64_t inc; //!< controls which RNG sequence (stream) is selected; must *always* be odd
};

void pcg32_srandom_r(const uint64_t initstate, const uint64_t initseq, struct pcg32_random* rng);


uint32_t pcg32_random_r(struct pcg32_random* rng);

uint32_t pcg32_boundedrand_r(const uint32_t bound, struct pcg32_random* rng);


//________________________________________________________________________________________________________________________
///
/// \brief PCG 64-bit random number generator state consisting of two 32-bit generators.
///
struct pcg32x2_random
{
struct pcg32_random gen[2]; //!< 32-bit RNG states
};

void pcg32x2_srandom_r(uint64_t seed1, uint64_t seed2, uint64_t seq1, uint64_t seq2, struct pcg32x2_random* rng);


uint64_t pcg32x2_random_r(struct pcg32x2_random* rng);

uint64_t pcg32x2_boundedrand_r(const uint64_t bound, struct pcg32x2_random* rng);
Loading

0 comments on commit 11b9527

Please sign in to comment.