From ba4a6a7eee7229277c4ccf216e8d7750a84f4903 Mon Sep 17 00:00:00 2001 From: "Christian B. Mendl" Date: Sat, 13 Jul 2024 23:53:38 +0200 Subject: [PATCH] basic DMRG demo applied to the Fermi-Hubbard model --- CMakeLists.txt | 7 +- examples/dmrg/basic_dmrg_fermi_hubbard.c | 124 +++++++++++++++++++++++ src/algorithm/dmrg.c | 2 +- src/algorithm/dmrg.h | 2 +- src/mps/mps.c | 93 ++++++++++++++++- src/mps/mps.h | 2 + src/operator/hamiltonian.c | 21 ++-- src/operator/hamiltonian.h | 9 ++ src/tensor/block_sparse_tensor.c | 16 +-- src/tensor/dense_tensor.c | 27 ++--- src/util/rng.c | 10 ++ src/util/rng.h | 2 + src/util/util.h | 20 ++++ 13 files changed, 284 insertions(+), 51 deletions(-) create mode 100644 examples/dmrg/basic_dmrg_fermi_hubbard.c diff --git a/CMakeLists.txt b/CMakeLists.txt index 8f69c78..62c692c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -13,9 +13,12 @@ set(CHEMTENSOR_DIRS "src" "src/tensor" "src/mps" "src/operator" "src/algorithm" 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}) - +add_executable( chemtensor_test ${CHEMTENSOR_SOURCES} ${TEST_SOURCES}) target_include_directories(chemtensor_test PRIVATE ${CHEMTENSOR_DIRS} ${BLAS_INCLUDE_DIRS} ${LAPACKE_INCLUDE_DIRS} ${HDF5_INCLUDE_DIRS}) target_link_libraries( chemtensor_test PRIVATE ${BLAS_LIBRARIES} ${LAPACKE_LIBRARIES} lapacke ${HDF5_LIBRARIES}) +add_executable( basic_dmrg_fermi_hubbard ${CHEMTENSOR_SOURCES} "examples/dmrg/basic_dmrg_fermi_hubbard.c") +target_include_directories(basic_dmrg_fermi_hubbard PRIVATE ${CHEMTENSOR_DIRS} ${BLAS_INCLUDE_DIRS} ${LAPACKE_INCLUDE_DIRS} ${HDF5_INCLUDE_DIRS}) +target_link_libraries( basic_dmrg_fermi_hubbard PRIVATE ${BLAS_LIBRARIES} ${LAPACKE_LIBRARIES} lapacke ${HDF5_LIBRARIES}) + add_test(NAME chemtensor_test COMMAND chemtensor_test) diff --git a/examples/dmrg/basic_dmrg_fermi_hubbard.c b/examples/dmrg/basic_dmrg_fermi_hubbard.c new file mode 100644 index 0000000..70d63d7 --- /dev/null +++ b/examples/dmrg/basic_dmrg_fermi_hubbard.c @@ -0,0 +1,124 @@ +#include +#include +#include +#include "hamiltonian.h" +#include "dmrg.h" +#include "aligned_memory.h" +#include "rng.h" +#include "util.h" + + +int main() +{ + // number of spin-endowed lattice sites (local dimension is 4) + const int nsites = 6; + printf("number of spin-endowed lattice sites: %i -> Hilbert space dimension: %li\n", nsites, ipow(4, nsites)); + + // Hamiltonian parameters + const double t = 1.0; + const double u = 4.0; + const double mu = 1.5; + printf("Fermi-Hubbard Hamiltonian parameters: t = %g, u = %g, mu = %g\n", t, u, mu); + + // overall quantum number sector of quantum state (particle number and spin) + const qnumber q_pnum = 7; + const qnumber q_spin = 1; + printf("psi restricted to overall quantum number sector with particle number %i and spin %i (integer convention)\n", q_pnum, q_spin); + const qnumber qnum_sector = fermi_hubbard_encode_quantum_numbers(q_pnum, q_spin); + + // maximum allowed MPS bond dimension + const long max_vdim = 32; + + struct mpo hamiltonian; + { + struct mpo_assembly assembly; + construct_fermi_hubbard_1d_mpo_assembly(nsites, t, u, mu, &assembly); + mpo_from_assembly(&assembly, &hamiltonian); + delete_mpo_assembly(&assembly); + } + if (!mpo_is_consistent(&hamiltonian)) { + fprintf(stderr, "internal consistency check for Fermi-Hubbard Hamiltonian MPO failed"); + return -1; + } + + // initial state vector as MPS + struct mps psi; + { + struct rng_state rng_state; + seed_rng_state(42, &rng_state); + + construct_random_mps(hamiltonian.a[0].dtype, hamiltonian.nsites, hamiltonian.d, hamiltonian.qsite, qnum_sector, max_vdim, &rng_state, &psi); + if (!mps_is_consistent(&psi)) { + fprintf(stderr, "internal MPS consistency check failed"); + return -1; + } + } + + // run two-site DMRG + const int num_sweeps = 4; + const int maxiter_lanczos = 25; + double tol_split = 1e-8; + double* en_sweeps = aligned_alloc(MEM_DATA_ALIGN, num_sweeps * sizeof(double)); + double* entropy = aligned_alloc(MEM_DATA_ALIGN, (nsites - 1) * sizeof(double)); + printf("Running two-site DMRG (using num_sweeps = %i, maxiter_lanczos = %i, tol_split = %g)... ", num_sweeps, maxiter_lanczos, tol_split); + if (dmrg_twosite(&hamiltonian, num_sweeps, maxiter_lanczos, tol_split, max_vdim, &psi, en_sweeps, entropy) < 0) { + fprintf(stderr, "'dmrg_twosite' failed internally"); + return -2; + } + printf("Done.\n"); + + printf("||psi||: %g\n", mps_norm(&psi)); + + printf("psi bond dimensions: "); + for (int l = 0; l < nsites + 1; l++) { + printf("%li ", mps_bond_dim(&psi, l)); + } + printf("\n"); + printf("splitting entropies for each bond: "); + for (int l = 0; l < nsites - 1; l++) { + printf("%g ", entropy[l]); + } + printf("\n"); + + // matrix representation of Hamiltonian + struct dense_tensor hmat; + { + struct block_sparse_tensor hsparse; + mpo_to_matrix(&hamiltonian, &hsparse); + + // convert to dense tensor + block_sparse_to_dense_tensor(&hsparse, &hmat); + delete_block_sparse_tensor(&hsparse); + + // dimension and data type consistency checks + assert(hmat.ndim == 4); + // dummy virtual bond dimensions are retained + assert(hmat.dim[0] == 1 && hmat.dim[3] == 1); + assert(hmat.dim[1] == ipow(hamiltonian.d, hamiltonian.nsites)); + assert(hmat.dim[1] == hmat.dim[2]); + assert(hmat.dtype == DOUBLE_REAL); + } + + // reference eigenvalues (based on exact diagonalization of matrix representation) + double* w_ref = aligned_alloc(MEM_DATA_ALIGN, hmat.dim[1] * sizeof(double)); + int info = LAPACKE_dsyev(LAPACK_ROW_MAJOR, 'N', 'U', hmat.dim[1], hmat.data, hmat.dim[2], w_ref); + if (info != 0) { + fprintf(stderr, "LAPACK function 'dsyev' failed, return value: %i\n", info); + return info; + } + printf("reference ground state energy: %g\n", w_ref[0]); + + printf("energy after each DMRG sweep:\n"); + for (int i = 0; i < num_sweeps; i++) { + printf("%i: %g, diff to reference: %g\n", i + 1, en_sweeps[i], en_sweeps[i] - w_ref[0]); + } + + aligned_free(w_ref); + delete_dense_tensor(&hmat); + aligned_free(entropy); + aligned_free(en_sweeps); + delete_mps(&psi); + delete_mpo(&hamiltonian); + + return 0; +} diff --git a/src/algorithm/dmrg.c b/src/algorithm/dmrg.c index 3b38ccb..9c97507 100644 --- a/src/algorithm/dmrg.c +++ b/src/algorithm/dmrg.c @@ -246,7 +246,7 @@ int dmrg_singlesite(const struct mpo* hamiltonian, const int num_sweeps, const i /// The input 'psi' is used as starting state and is updated in-place during the optimization. /// int dmrg_twosite(const struct mpo* hamiltonian, const int num_sweeps, const int maxiter_lanczos, const double tol_split, const long max_vdim, - struct mps* psi, double* en_sweeps, double* restrict entropy) + struct mps* psi, double* restrict en_sweeps, double* restrict entropy) { // number of lattice sites const int nsites = hamiltonian->nsites; diff --git a/src/algorithm/dmrg.h b/src/algorithm/dmrg.h index 83c23ee..17378f8 100644 --- a/src/algorithm/dmrg.h +++ b/src/algorithm/dmrg.h @@ -10,4 +10,4 @@ int dmrg_singlesite(const struct mpo* hamiltonian, const int num_sweeps, const int maxiter_lanczos, struct mps* psi, double* en_sweeps); int dmrg_twosite(const struct mpo* hamiltonian, const int num_sweeps, const int maxiter_lanczos, const double tol_split, const long max_vdim, - struct mps* psi, double* en_sweeps, double* restrict entropy); + struct mps* psi, double* restrict en_sweeps, double* restrict entropy); diff --git a/src/mps/mps.c b/src/mps/mps.c index aef7ff9..133d0b0 100644 --- a/src/mps/mps.c +++ b/src/mps/mps.c @@ -55,6 +55,97 @@ void delete_mps(struct mps* mps) } +//________________________________________________________________________________________________________________________ +/// +/// \brief Construct a matrix product state with random normal tensor entries, given an overall quantum number sector and maximum virtual bond dimension. +/// +void construct_random_mps(const enum numeric_type dtype, const int nsites, const long d, const qnumber* qsite, const qnumber qnum_sector, const long max_vdim, struct rng_state* rng_state, struct mps* mps) +{ + assert(nsites >= 1); + assert(d >= 1); + + // virtual bond quantum numbers + long* dim_bonds = aligned_alloc(MEM_DATA_ALIGN, (nsites + 1) * sizeof(long)); + qnumber** qbonds = aligned_alloc(MEM_DATA_ALIGN, (nsites + 1) * sizeof(qnumber*)); + // dummy left virtual bond; set quantum number to zero + dim_bonds[0] = 1; + qbonds[0] = aligned_alloc(MEM_DATA_ALIGN, sizeof(qnumber)); + qbonds[0][0] = 0; + // dummy right virtual bond; set quantum number to overall quantum number sector + dim_bonds[nsites] = 1; + qbonds[nsites] = aligned_alloc(MEM_DATA_ALIGN, sizeof(qnumber)); + qbonds[nsites][0] = qnum_sector; + // virtual bond quantum numbers on left half + for (int l = 1; l < (nsites + 1) / 2; l++) + { + // enumerate all combinations of left bond quantum numbers and local physical quantum numbers + const long dim_full = dim_bonds[l - 1] * d; + qnumber* qnums_full = aligned_alloc(MEM_DATA_ALIGN, dim_full * sizeof(qnumber)); + for (long i = 0; i < dim_bonds[l - 1]; i++) { + for (long j = 0; j < d; j++) { + qnums_full[i*d + j] = qbonds[l - 1][i] + qsite[j]; + } + } + dim_bonds[l] = lmin(dim_full, max_vdim); + qbonds[l] = aligned_alloc(MEM_DATA_ALIGN, dim_bonds[l] * sizeof(qnumber)); + if (dim_full <= max_vdim) { + memcpy(qbonds[l], qnums_full, dim_bonds[l] * sizeof(qnumber)); + } + else { + // randomly select quantum numbers + for (long i = 0; i < dim_bonds[l]; i++) { + qbonds[l][i] = qnums_full[rand_interval(dim_full, rng_state)]; + } + } + aligned_free(qnums_full); + } + // virtual bond quantum numbers on right half + for (int l = nsites - 1; l >= (nsites + 1) / 2; l--) + { + // enumerate all combinations of right bond quantum numbers and local physical quantum numbers + const long dim_full = dim_bonds[l + 1] * d; + qnumber* qnums_full = aligned_alloc(MEM_DATA_ALIGN, dim_full * sizeof(qnumber)); + for (long i = 0; i < dim_bonds[l + 1]; i++) { + for (long j = 0; j < d; j++) { + qnums_full[i*d + j] = qbonds[l + 1][i] - qsite[j]; + } + } + dim_bonds[l] = lmin(dim_full, max_vdim); + qbonds[l] = aligned_alloc(MEM_DATA_ALIGN, dim_bonds[l] * sizeof(qnumber)); + if (dim_full <= max_vdim) { + memcpy(qbonds[l], qnums_full, dim_bonds[l] * sizeof(qnumber)); + } + else { + // randomly select quantum numbers + for (long i = 0; i < dim_bonds[l]; i++) { + qbonds[l][i] = qnums_full[rand_interval(dim_full, rng_state)]; + } + } + aligned_free(qnums_full); + } + + allocate_mps(dtype, nsites, d, qsite, dim_bonds, (const qnumber**)qbonds, mps); + + for (int l = 0; l < nsites + 1; l++) { + aligned_free(qbonds[l]); + } + aligned_free(qbonds); + aligned_free(dim_bonds); + + // fill MPS tensor entries with pseudo-random numbers, scaled by 1 / sqrt("number of entries") + for (int l = 0; l < nsites; l++) + { + // logical number of entries in MPS tensor + const long nelem = integer_product(mps->a[l].dim_logical, mps->a[l].ndim); + // ensure that 'alpha' is large enough to store any numeric type + dcomplex alpha; + assert(mps->a[l].dtype == dtype); + numeric_from_double(1.0 / sqrt(nelem), dtype, &alpha); + block_sparse_tensor_fill_random_normal(&alpha, numeric_zero(dtype), rng_state, &mps->a[l]); + } +} + + //________________________________________________________________________________________________________________________ /// /// \brief Internal consistency check of the MPS data structure. @@ -318,7 +409,7 @@ void mps_local_orthonormalize_qr(struct block_sparse_tensor* restrict a, struct //________________________________________________________________________________________________________________________ /// -/// \brief Right-orthonormalize a local MPS site tensor by RQ decomposition, and update tensor at next site. +/// \brief Right-orthonormalize a local MPS site tensor by RQ decomposition, and update tensor at previous site. /// void mps_local_orthonormalize_rq(struct block_sparse_tensor* restrict a, struct block_sparse_tensor* restrict a_prev) { diff --git a/src/mps/mps.h b/src/mps/mps.h index 90e8664..183896a 100644 --- a/src/mps/mps.h +++ b/src/mps/mps.h @@ -29,6 +29,8 @@ void allocate_mps(const enum numeric_type dtype, const int nsites, const long d, void delete_mps(struct mps* mps); +void construct_random_mps(const enum numeric_type dtype, const int nsites, const long d, const qnumber* qsite, const qnumber qnum_sector, const long max_vdim, struct rng_state* rng_state, struct mps* mps); + bool mps_is_consistent(const struct mps* mps); diff --git a/src/operator/hamiltonian.c b/src/operator/hamiltonian.c index b2ad328..8d2a77a 100644 --- a/src/operator/hamiltonian.c +++ b/src/operator/hamiltonian.c @@ -14,7 +14,7 @@ //________________________________________________________________________________________________________________________ /// -/// \brief Construct a Hamiltonian as MPO based on local operator chains, which are shifted along a 1D lattice. +/// \brief Construct a Hamiltonian as MPO operator graph based on local operator chains, which are shifted along a 1D lattice. /// static void local_opchains_to_mpo_graph(const int nsites, const struct op_chain* lopchains, const int nlopchains, struct mpo_graph* graph) { @@ -244,10 +244,10 @@ void construct_fermi_hubbard_1d_mpo_assembly(const int nsites, const double t, c const qnumber qs[4] = { 0, -1, 1, 0 }; assembly->d = 4; assembly->qsite = aligned_alloc(MEM_DATA_ALIGN, assembly->d * sizeof(qnumber)); - assembly->qsite[0] = (qn[0] << 16) + qs[0]; - assembly->qsite[1] = (qn[1] << 16) + qs[1]; - assembly->qsite[2] = (qn[2] << 16) + qs[2]; - assembly->qsite[3] = (qn[3] << 16) + qs[3]; + assembly->qsite[0] = fermi_hubbard_encode_quantum_numbers(qn[0], qs[0]); + assembly->qsite[1] = fermi_hubbard_encode_quantum_numbers(qn[1], qs[1]); + assembly->qsite[2] = fermi_hubbard_encode_quantum_numbers(qn[2], qs[2]); + assembly->qsite[3] = fermi_hubbard_encode_quantum_numbers(qn[3], qs[3]); assembly->dtype = DOUBLE_REAL; @@ -343,16 +343,13 @@ void construct_fermi_hubbard_1d_mpo_assembly(const int nsites, const double t, c assembly->coeffmap = aligned_alloc(MEM_DATA_ALIGN, sizeof(coeffmap)); memcpy(assembly->coeffmap, coeffmap, sizeof(coeffmap)); - // cast to unsigned integer to avoid compiler warning when bit-shifting - const unsigned int n1 = -1; - // local two-site and single-site terms // spin-up kinetic hopping - int oids_c0[] = { 3, 2 }; qnumber qnums_c0[] = { 0, ( 1 << 16) + 1, 0 }; - int oids_c1[] = { 4, 1 }; qnumber qnums_c1[] = { 0, (n1 << 16) - 1, 0 }; + int oids_c0[] = { 3, 2 }; qnumber qnums_c0[] = { 0, fermi_hubbard_encode_quantum_numbers( 1, 1), 0 }; + int oids_c1[] = { 4, 1 }; qnumber qnums_c1[] = { 0, fermi_hubbard_encode_quantum_numbers(-1, -1), 0 }; // spin-down kinetic hopping - int oids_c2[] = { 5, 8 }; qnumber qnums_c2[] = { 0, ( 1 << 16) - 1, 0 }; - int oids_c3[] = { 6, 7 }; qnumber qnums_c3[] = { 0, (n1 << 16) + 1, 0 }; + int oids_c2[] = { 5, 8 }; qnumber qnums_c2[] = { 0, fermi_hubbard_encode_quantum_numbers( 1, -1), 0 }; + int oids_c3[] = { 6, 7 }; qnumber qnums_c3[] = { 0, fermi_hubbard_encode_quantum_numbers(-1, 1), 0 }; // number operator - mu (n_up + n_dn) and interaction u (n_up-1/2) (n_dn-1/2) int oids_c4[] = { 9 }; qnumber qnums_c4[] = { 0, 0 }; int oids_c5[] = { 10 }; qnumber qnums_c5[] = { 0, 0 }; diff --git a/src/operator/hamiltonian.h b/src/operator/hamiltonian.h index 71e9cce..e3fbb17 100644 --- a/src/operator/hamiltonian.h +++ b/src/operator/hamiltonian.h @@ -12,6 +12,15 @@ void construct_heisenberg_xxz_1d_mpo_assembly(const int nsites, const double J, void construct_bose_hubbard_1d_mpo_assembly(const int nsites, const long d, const double t, const double u, const double mu, struct mpo_assembly* assembly); +//________________________________________________________________________________________________________________________ +/// +/// \brief Encode a particle number and spin quantum number for the Fermi-Hubbard Hamiltonian into a single quantum number. +/// +static inline qnumber fermi_hubbard_encode_quantum_numbers(const qnumber q_pnum, const qnumber q_spin) +{ + return (q_pnum << 16) + q_spin; +} + void construct_fermi_hubbard_1d_mpo_assembly(const int nsites, const double t, const double u, const double mu, struct mpo_assembly* assembly); void construct_molecular_hamiltonian_mpo_assembly(const struct dense_tensor* restrict tkin, const struct dense_tensor* restrict vint, struct mpo_assembly* assembly); diff --git a/src/tensor/block_sparse_tensor.c b/src/tensor/block_sparse_tensor.c index 44a9971..e39ee58 100644 --- a/src/tensor/block_sparse_tensor.c +++ b/src/tensor/block_sparse_tensor.c @@ -1744,16 +1744,6 @@ void block_sparse_tensor_dot(const struct block_sparse_tensor* restrict s, const } -//________________________________________________________________________________________________________________________ -/// -/// \brief Minimum of two integers. -/// -static inline long minl(const long a, const long b) -{ - return (a <= b) ? a : b; -} - - //________________________________________________________________________________________________________________________ /// /// \brief Compute the logical QR decomposition of a block-sparse matrix. @@ -1790,7 +1780,7 @@ int block_sparse_tensor_qr(const struct block_sparse_tensor* restrict a, struct assert(b->ndim == 2); assert(b->dtype == a->dtype); - const long k = minl(b->dim[0], b->dim[1]); + const long k = lmin(b->dim[0], b->dim[1]); // append a sequence of logical quantum numbers of length k for (long l = 0; l < k; l++) { @@ -1929,7 +1919,7 @@ int block_sparse_tensor_rq(const struct block_sparse_tensor* restrict a, struct assert(b->ndim == 2); assert(b->dtype == a->dtype); - const long k = minl(b->dim[0], b->dim[1]); + const long k = lmin(b->dim[0], b->dim[1]); // append a sequence of logical quantum numbers of length k for (long l = 0; l < k; l++) { @@ -2068,7 +2058,7 @@ int block_sparse_tensor_svd(const struct block_sparse_tensor* restrict a, struct assert(b->ndim == 2); assert(b->dtype == a->dtype); - const long k = minl(b->dim[0], b->dim[1]); + const long k = lmin(b->dim[0], b->dim[1]); // append a sequence of logical quantum numbers of length k for (long l = 0; l < k; l++) { diff --git a/src/tensor/dense_tensor.c b/src/tensor/dense_tensor.c index 192458b..660ebfa 100644 --- a/src/tensor/dense_tensor.c +++ b/src/tensor/dense_tensor.c @@ -1670,21 +1670,6 @@ void dense_tensor_kronecker_product(const struct dense_tensor* restrict s, const } -//________________________________________________________________________________________________________________________ -/// -/// \brief Minimum of two integers. -/// -static inline long minl(const long a, const long b) -{ - if (a <= b) { - return a; - } - else { - return b; - } -} - - //________________________________________________________________________________________________________________________ /// /// \brief Compute the QR decomposition of the matrix 'a', and store the result in 'q' and 'r' (will be allocated). @@ -1697,7 +1682,7 @@ int dense_tensor_qr(const struct dense_tensor* restrict a, struct dense_tensor* const long m = a->dim[0]; const long n = a->dim[1]; - const long k = minl(m, n); + const long k = lmin(m, n); const long dim_q[2] = { m, k }; const long dim_r[2] = { k, n }; @@ -1725,7 +1710,7 @@ int dense_tensor_qr_fill(const struct dense_tensor* restrict a, struct dense_ten const long m = a->dim[0]; const long n = a->dim[1]; - const long k = minl(m, n); + const long k = lmin(m, n); assert(q->dim[0] == m); assert(q->dim[1] == k); @@ -2029,7 +2014,7 @@ int dense_tensor_rq(const struct dense_tensor* restrict a, struct dense_tensor* const long m = a->dim[0]; const long n = a->dim[1]; - const long k = minl(m, n); + const long k = lmin(m, n); const long dim_r[2] = { m, k }; const long dim_q[2] = { k, n }; @@ -2057,7 +2042,7 @@ int dense_tensor_rq_fill(const struct dense_tensor* restrict a, struct dense_ten const long m = a->dim[0]; const long n = a->dim[1]; - const long k = minl(m, n); + const long k = lmin(m, n); assert(r->dim[0] == m); assert(r->dim[1] == k); @@ -2361,7 +2346,7 @@ int dense_tensor_svd(const struct dense_tensor* restrict a, struct dense_tensor* const long m = a->dim[0]; const long n = a->dim[1]; - const long k = minl(m, n); + const long k = lmin(m, n); const long dim_u[2] = { m, k }; const long dim_s[1] = { k }; @@ -2392,7 +2377,7 @@ int dense_tensor_svd_fill(const struct dense_tensor* restrict a, struct dense_te const long m = a->dim[0]; const long n = a->dim[1]; - const long k = minl(m, n); + const long k = lmin(m, n); assert( u->dim[0] == m); assert( u->dim[1] == k); diff --git a/src/util/rng.c b/src/util/rng.c index 6c4f745..436fff4 100644 --- a/src/util/rng.c +++ b/src/util/rng.c @@ -45,6 +45,16 @@ uint64_t rand_uint64(struct rng_state* state) } +//________________________________________________________________________________________________________________________ +/// +/// \brief Generate a uniformly distributed integer random number from the interval [0, bound). +/// +uint64_t rand_interval(const uint64_t bound, struct rng_state* state) +{ + return pcg32x2_boundedrand_r(bound, &state->pcgstate); +} + + //________________________________________________________________________________________________________________________ /// /// \brief Generate a uniformly distributed pseudo-random number from the interval [0, 1) diff --git a/src/util/rng.h b/src/util/rng.h index 671b93d..f575736 100644 --- a/src/util/rng.h +++ b/src/util/rng.h @@ -29,6 +29,8 @@ uint32_t rand_uint32(struct rng_state* state); uint64_t rand_uint64(struct rng_state* state); +uint64_t rand_interval(const uint64_t bound, struct rng_state* state); + //________________________________________________________________________________________________________________________ // diff --git a/src/util/util.h b/src/util/util.h index 5f55f60..a0d9ca5 100644 --- a/src/util/util.h +++ b/src/util/util.h @@ -29,6 +29,26 @@ static inline int imax(const int a, const int b) } +//________________________________________________________________________________________________________________________ +/// +/// \brief Minimum of two long integers. +/// +static inline long lmin(const long a, const long b) +{ + return (a <= b) ? a : b; +} + + +//________________________________________________________________________________________________________________________ +/// +/// \brief Maximum of two long integers. +/// +static inline long lmax(const long a, const long b) +{ + return (a >= b) ? a : b; +} + + //________________________________________________________________________________________________________________________ /// /// \brief Square function x -> x^2.