diff --git a/src/tensor/dense_tensor.c b/src/tensor/dense_tensor.c index 86023bb..1895098 100644 --- a/src/tensor/dense_tensor.c +++ b/src/tensor/dense_tensor.c @@ -2667,6 +2667,70 @@ bool dense_tensor_allclose(const struct dense_tensor* restrict s, const struct d } +//________________________________________________________________________________________________________________________ +/// +/// \brief Test whether a dense tensors is entrywise zero within tolerance 'tol'. +/// +/// To test whether the tensor is exactly zero, set 'tol = 0.0'. +/// +bool dense_tensor_is_zero(const struct dense_tensor* t, const double tol) +{ + const long nelem = dense_tensor_num_elements(t); + + switch (t->dtype) + { + case CT_SINGLE_REAL: + { + const float* data = t->data; + for (long i = 0; i < nelem; i++) { + if (fabsf(data[i]) > tol) { + return false; + } + } + break; + } + case CT_DOUBLE_REAL: + { + const double* data = t->data; + for (long i = 0; i < nelem; i++) { + if (fabs(data[i]) > tol) { + return false; + } + } + break; + } + case CT_SINGLE_COMPLEX: + { + const scomplex* data = t->data; + for (long i = 0; i < nelem; i++) { + if (cabsf(data[i]) > tol) { + return false; + } + } + break; + } + case CT_DOUBLE_COMPLEX: + { + const dcomplex* data = t->data; + for (long i = 0; i < nelem; i++) { + if (cabs(data[i]) > tol) { + return false; + } + } + break; + } + default: + { + // unknown data type + assert(false); + return false; + } + } + + return true; +} + + //________________________________________________________________________________________________________________________ /// /// \brief Test whether a dense tensors is close to the identity matrix within tolerance 'tol'. @@ -2746,6 +2810,7 @@ bool dense_tensor_is_identity(const struct dense_tensor* t, const double tol) return true; } + //________________________________________________________________________________________________________________________ /// /// \brief Test whether a dense tensors is an isometry within tolerance 'tol'. diff --git a/src/tensor/dense_tensor.h b/src/tensor/dense_tensor.h index 697f6f3..9fb9732 100644 --- a/src/tensor/dense_tensor.h +++ b/src/tensor/dense_tensor.h @@ -215,6 +215,8 @@ void dense_tensor_block(const struct dense_tensor* restrict t, const long* restr bool dense_tensor_allclose(const struct dense_tensor* restrict s, const struct dense_tensor* restrict t, const double tol); +bool dense_tensor_is_zero(const struct dense_tensor* t, const double tol); + bool dense_tensor_is_identity(const struct dense_tensor* t, const double tol); bool dense_tensor_is_isometry(const struct dense_tensor* t, const double tol, const bool transpose); diff --git a/src/tensor/su2_tensor.c b/src/tensor/su2_tensor.c index 7cf5066..a3c8a20 100644 --- a/src/tensor/su2_tensor.c +++ b/src/tensor/su2_tensor.c @@ -308,6 +308,372 @@ void su2_tensor_fmove(const struct su2_tensor* restrict t, const int i_ax, struc } +//________________________________________________________________________________________________________________________ +/// +/// \brief Contract two SU(2) tensors for the scenario that the to-be contracted axes form subtrees with matching topology and +/// that the resulting fusion-splitting tree is again simple. +/// +void su2_tensor_contract_simple(const struct su2_tensor* restrict s, const int* restrict i_ax_s, const struct su2_tensor* restrict t, const int* restrict i_ax_t, const int ndim_mult, struct su2_tensor* restrict r) +{ + assert(s->dtype == t->dtype); + + // dimension and quantum number compatibility checks + assert(ndim_mult >= 1); + for (int i = 0; i < ndim_mult; i++) + { + assert(0 <= i_ax_s[i] && i_ax_s[i] < s->ndim_logical + s->ndim_auxiliary); + assert(0 <= i_ax_t[i] && i_ax_t[i] < t->ndim_logical + t->ndim_auxiliary); + + // must both be logical or both be auxiliary axes + assert((i_ax_s[i] < s->ndim_logical && i_ax_t[i] < t->ndim_logical) || (i_ax_s[i] >= s->ndim_logical && i_ax_t[i] >= t->ndim_logical)); + + const struct su2_irreducible_list* irred_s = &s->outer_jlists[i_ax_s[i]]; + const struct su2_irreducible_list* irred_t = &t->outer_jlists[i_ax_t[i]]; + assert(su2_irreducible_list_equal(irred_s, irred_t)); + if (i_ax_s[i] < s->ndim_logical) { + for (int k = 0; k < irred_s->num; k++) { + // degeneracy dimensions for to-be contracted axes must match + assert(s->dim_degen[i_ax_s[i]][irred_s->jlist[k]] == + t->dim_degen[i_ax_t[i]][irred_t->jlist[k]]); + } + } + } + + const int ndim_outer_s = s->ndim_logical + s->ndim_auxiliary; + const int ndim_outer_t = t->ndim_logical + t->ndim_auxiliary; + + const struct su2_tree_node* subtree_s = su2_subtree_with_leaf_axes(s->tree.tree_split, i_ax_s, ndim_mult); + const struct su2_tree_node* subtree_t = su2_subtree_with_leaf_axes(t->tree.tree_fuse, i_ax_t, ndim_mult); + assert(subtree_s != NULL); + assert(subtree_t != NULL); + assert(su2_tree_equal_topology(subtree_s, subtree_t)); + // at least one of the subtrees must actually be the full tree, otherwise the resulting fusion-splitting tree is not simple + assert((subtree_s == s->tree.tree_split) || (subtree_t == t->tree.tree_fuse)); + // at least one of the subtree roots must be an internal axis (to ensure consistency of axis enumeration) + assert((subtree_s->i_ax >= ndim_outer_s) || (subtree_t->i_ax >= ndim_outer_t)); + + const int ndim_s = su2_tensor_ndim(s); + const int ndim_t = su2_tensor_ndim(t); + + int* i_ax_subtree_s = ct_malloc(ndim_s * sizeof(int)); + int* i_ax_subtree_t = ct_malloc(ndim_t * sizeof(int)); + const int ndim_subtree_s = su2_tree_axes_list(subtree_s, i_ax_subtree_s); + const int ndim_subtree_t = su2_tree_axes_list(subtree_t, i_ax_subtree_t); + assert(ndim_subtree_s == ndim_subtree_t); // also contains internal axes + // ensure that axes paired according to tree topology match the provided to-be contracted axes indices + for (int i = 0; i < ndim_mult; i++) { + bool found = false; + for (int j = 0; j < ndim_subtree_s; j++) { + if (i_ax_s[i] == i_ax_subtree_s[j]) { + assert(i_ax_t[i] == i_ax_subtree_t[j]); + found = true; + break; + } + } + assert(found); + } + + // map axes of 's' and 't' to axes of contracted tensor + int* axis_map_s = ct_calloc(ndim_s, sizeof(int)); + int* axis_map_t = ct_calloc(ndim_t, sizeof(int)); + // mark to-be contracted axes as invalid + for (int i = 0; i < ndim_subtree_s; i++) { + assert(0 <= i_ax_subtree_s[i] && i_ax_subtree_s[i] < ndim_s); + axis_map_s[i_ax_subtree_s[i]] = -1; + } + for (int i = 0; i < ndim_subtree_t; i++) { + assert(0 <= i_ax_subtree_t[i] && i_ax_subtree_t[i] < ndim_t); + axis_map_t[i_ax_subtree_t[i]] = -1; + } + // retain one of the subtree roots as axis + assert(axis_map_s[subtree_s->i_ax] == -1); + assert(axis_map_t[subtree_t->i_ax] == -1); + if (subtree_t->i_ax >= ndim_outer_t) { // subtree root in 't' is an internal axis, does not need to be retained + // retain subtree root in 's' + axis_map_s[subtree_s->i_ax] = 0; + } + else { + assert(subtree_s->i_ax >= ndim_outer_s); // subtree root in 's' is an internal axis, does not need to be retained + // retain subtree root in 't' + axis_map_t[subtree_t->i_ax] = 0; + } + // new logical axes + int ndim_logical_r = 0; + for (int i = 0; i < s->ndim_logical; i++) { + if (axis_map_s[i] != -1) { + axis_map_s[i] = ndim_logical_r++; + } + } + for (int i = 0; i < t->ndim_logical; i++) { + if (axis_map_t[i] != -1) { + axis_map_t[i] = ndim_logical_r++; + } + } + // new auxiliary axes + int ndim_auxiliary_r = 0; + for (int i = s->ndim_logical; i < s->ndim_logical + s->ndim_auxiliary; i++) { + if (axis_map_s[i] != -1) { + axis_map_s[i] = ndim_logical_r + ndim_auxiliary_r++; + } + } + for (int i = t->ndim_logical; i < t->ndim_logical + t->ndim_auxiliary; i++) { + if (axis_map_t[i] != -1) { + axis_map_t[i] = ndim_logical_r + ndim_auxiliary_r++; + } + } + // new internal axes + const int ndim_outer_r = ndim_logical_r + ndim_auxiliary_r; + int ndim_internal_r = 0; + for (int i = s->ndim_logical + s->ndim_auxiliary; i < ndim_s; i++) { + if (axis_map_s[i] != -1) { + axis_map_s[i] = ndim_outer_r + ndim_internal_r++; + } + } + for (int i = t->ndim_logical + t->ndim_auxiliary; i < ndim_t; i++) { + if (axis_map_t[i] != -1) { + axis_map_t[i] = ndim_outer_r + ndim_internal_r++; + } + } + assert(ndim_internal_r == ndim_outer_r - 3); + // identify roots of subtrees as same axis (after contracting subtrees to identity) + if (axis_map_t[subtree_t->i_ax] == -1) { + assert(axis_map_s[subtree_s->i_ax] >= 0); + axis_map_t[subtree_t->i_ax] = axis_map_s[subtree_s->i_ax]; + } + else { + assert(axis_map_s[subtree_s->i_ax] == -1); + axis_map_s[subtree_s->i_ax] = axis_map_t[subtree_t->i_ax]; + } + + struct su2_fuse_split_tree mapped_tree_s, mapped_tree_t; + copy_su2_fuse_split_tree(&s->tree, &mapped_tree_s); + copy_su2_fuse_split_tree(&t->tree, &mapped_tree_t); + su2_fuse_split_tree_update_axes_indices(&mapped_tree_s, axis_map_s); + su2_fuse_split_tree_update_axes_indices(&mapped_tree_t, axis_map_t); + + // construct the new fusion-splitting tree + struct su2_fuse_split_tree tree_r; + tree_r.ndim = ndim_logical_r + ndim_auxiliary_r + ndim_internal_r; + + if (subtree_s == s->tree.tree_split) + { + tree_r.tree_split = mapped_tree_t.tree_split; + + if (subtree_t == t->tree.tree_fuse) + { + tree_r.tree_fuse = mapped_tree_s.tree_fuse; + } + else + { + tree_r.tree_fuse = mapped_tree_t.tree_fuse; + // make a copy to avoid freeing same memory twice + struct su2_tree_node* tree_fuse_s = ct_malloc(sizeof(struct su2_tree_node)); + copy_su2_tree(mapped_tree_s.tree_fuse, tree_fuse_s); + struct su2_tree_node* old_subtree = su2_tree_replace_subtree(tree_r.tree_fuse, axis_map_t[subtree_t->i_ax], tree_fuse_s); + assert(old_subtree != NULL); + delete_su2_tree(old_subtree); + ct_free(old_subtree); // free actual node + } + } + else + { + assert(subtree_t == t->tree.tree_fuse); + tree_r.tree_fuse = mapped_tree_s.tree_fuse; + + tree_r.tree_split = mapped_tree_s.tree_split; + // make a copy to avoid freeing same memory twice + struct su2_tree_node* tree_split_t = ct_malloc(sizeof(struct su2_tree_node)); + copy_su2_tree(mapped_tree_t.tree_split, tree_split_t); + struct su2_tree_node* old_subtree = su2_tree_replace_subtree(tree_r.tree_split, axis_map_s[subtree_s->i_ax], tree_split_t); + assert(old_subtree != NULL); + delete_su2_tree(old_subtree); + ct_free(old_subtree); // free actual node + } + + struct su2_irreducible_list* outer_jlists_r = ct_calloc(ndim_outer_r, sizeof(struct su2_irreducible_list)); + for (int i = 0; i < ndim_outer_s; i++) { + if (axis_map_s[i] == -1) { + continue; + } + assert(axis_map_s[i] < ndim_outer_r); + copy_su2_irreducible_list(&s->outer_jlists[i], &outer_jlists_r[axis_map_s[i]]); + } + for (int i = 0; i < ndim_outer_t; i++) { + if (axis_map_t[i] == -1) { + continue; + } + assert(axis_map_t[i] < ndim_outer_r); + copy_su2_irreducible_list(&t->outer_jlists[i], &outer_jlists_r[axis_map_t[i]]); + } + + long** dim_degen_r = ct_calloc(ndim_logical_r, sizeof(long*)); + for (int is = 0; is < s->ndim_logical; is++) + { + if (axis_map_s[is] == -1) { + continue; + } + const int ir = axis_map_s[is]; + assert(0 <= ir && ir < ndim_logical_r); + assert(outer_jlists_r[ir].num > 0); + qnumber j_max = 0; + for (int k = 0; k < outer_jlists_r[ir].num; k++) { + j_max = imax(j_max, outer_jlists_r[ir].jlist[k]); + } + dim_degen_r[ir] = ct_malloc((j_max + 1) * sizeof(long)); + memcpy(dim_degen_r[ir], s->dim_degen[is], (j_max + 1) * sizeof(long)); + } + for (int it = 0; it < t->ndim_logical; it++) + { + if (axis_map_t[it] == -1) { + continue; + } + const int ir = axis_map_t[it]; + assert(0 <= ir && ir < ndim_logical_r); + assert(outer_jlists_r[ir].num > 0); + qnumber j_max = 0; + for (int k = 0; k < outer_jlists_r[ir].num; k++) { + j_max = imax(j_max, outer_jlists_r[ir].jlist[k]); + } + dim_degen_r[ir] = ct_malloc((j_max + 1) * sizeof(long)); + memcpy(dim_degen_r[ir], t->dim_degen[it], (j_max + 1) * sizeof(long)); + } + + allocate_su2_tensor(s->dtype, ndim_logical_r, ndim_auxiliary_r, &tree_r, outer_jlists_r, (const long**)dim_degen_r, r); + + // permutations of dense degeneracy tensors before contraction + int* perm_s = ct_malloc(s->ndim_logical * sizeof(int)); + int c = 0; + for (int i = 0; i < s->ndim_logical; i++) { + if (axis_map_s[i] != -1) { + perm_s[c++] = i; + } + } + for (int i = 0; i < ndim_mult; i++) { + // ignore auxiliary axes for permutation + if (i_ax_s[i] < s->ndim_logical) { + perm_s[c++] = i_ax_s[i]; + } + } + assert(c == s->ndim_logical); + int* perm_t = ct_malloc(t->ndim_logical * sizeof(int)); + c = 0; + for (int i = 0; i < ndim_mult; i++) { + // ignore auxiliary axes for permutation + if (i_ax_t[i] < t->ndim_logical) { + perm_t[c++] = i_ax_t[i]; + } + } + for (int i = 0; i < t->ndim_logical; i++) { + if (axis_map_t[i] != -1) { + perm_t[c++] = i; + } + } + assert(c == t->ndim_logical); + // number of to-be multiplied logical dimensions + int ndim_mult_logical = 0; + for (int i = 0; i < ndim_mult; i++) { + // ignore auxiliary axes + if (i_ax_s[i] < s->ndim_logical) { + ndim_mult_logical++; + } + } + + // contract degeneracy tensors + qnumber* jlist_r = ct_malloc(su2_tensor_ndim(r) * sizeof(qnumber)); + for (long cs = 0; cs < s->charge_sectors.nsec; cs++) + { + // current 'j' quantum numbers + const qnumber* jlist_s = &s->charge_sectors.jlists[cs * s->charge_sectors.ndim]; + + // corresponding "degeneracy" tensor + const struct dense_tensor* ds = s->degensors[cs]; + assert(ds->dtype == s->dtype); + assert(ds->ndim == s->ndim_logical); + struct dense_tensor ds_perm; + transpose_dense_tensor(perm_s, ds, &ds_perm); + + // fill 'j' quantum numbers for charge sector in to-be contracted tensor 'r' + for (int i = 0; i < ndim_s; i++) { + if (axis_map_s[i] == -1) { + continue; + } + jlist_r[axis_map_s[i]] = jlist_s[i]; + } + + for (long ct = 0; ct < t->charge_sectors.nsec; ct++) + { + // 'j' quantum numbers of current sector + const qnumber* jlist_t = &t->charge_sectors.jlists[ct * t->charge_sectors.ndim]; + + bool compatible_sector = true; + for (int i = 0; i < ndim_subtree_s; i++) { + if (jlist_s[i_ax_subtree_s[i]] != jlist_t[i_ax_subtree_t[i]]) { + compatible_sector = false; + break; + } + } + if (!compatible_sector) { + continue; + } + + // corresponding "degeneracy" tensor + const struct dense_tensor* dt = t->degensors[ct]; + assert(dt->dtype == t->dtype); + assert(dt->ndim == t->ndim_logical); + struct dense_tensor dt_perm; + transpose_dense_tensor(perm_t, dt, &dt_perm); + + // fill 'j' quantum numbers for charge sector in to-be contracted tensor 'r' + for (int i = 0; i < ndim_t; i++) { + if (axis_map_t[i] == -1) { + continue; + } + jlist_r[axis_map_t[i]] = jlist_t[i]; + } + + const long cr = charge_sector_index(&r->charge_sectors, jlist_r); + assert(cr != -1); + + // corresponding "degeneracy" tensor + struct dense_tensor* dr = r->degensors[cr]; + assert(dr->dtype == r->dtype); + assert(dr->ndim == r->ndim_logical); + + // actually multiply dense tensors and add result to 'dr' + dense_tensor_dot_update(numeric_one(s->dtype), &ds_perm, TENSOR_AXIS_RANGE_TRAILING, &dt_perm, TENSOR_AXIS_RANGE_LEADING, ndim_mult_logical, numeric_one(s->dtype), dr); + + delete_dense_tensor(&dt_perm); + } + + delete_dense_tensor(&ds_perm); + } + + ct_free(perm_t); + ct_free(perm_s); + + ct_free(jlist_r); + + for (int i = 0; i < ndim_logical_r; i++) { + ct_free(dim_degen_r[i]); + } + ct_free(dim_degen_r); + + for (int i = 0; i < ndim_outer_r; i++) { + delete_su2_irreducible_list(&outer_jlists_r[i]); + } + ct_free(outer_jlists_r); + + delete_su2_fuse_split_tree(&mapped_tree_t); + delete_su2_fuse_split_tree(&mapped_tree_s); + + ct_free(axis_map_t); + ct_free(axis_map_s); + + ct_free(i_ax_subtree_t); + ct_free(i_ax_subtree_s); +} + + //________________________________________________________________________________________________________________________ /// /// \brief Compute the lexicographically next quantum index. diff --git a/src/tensor/su2_tensor.h b/src/tensor/su2_tensor.h index a141526..1b80edd 100644 --- a/src/tensor/su2_tensor.h +++ b/src/tensor/su2_tensor.h @@ -69,6 +69,17 @@ enum tensor_axis_direction su2_tensor_logical_axis_direction(const struct su2_te bool su2_tensor_is_consistent(const struct su2_tensor* t); +//________________________________________________________________________________________________________________________ +// + +// flip interal fusin-splitting tree + +static inline void su2_tensor_flip_trees(struct su2_tensor* t) +{ + su2_fuse_split_tree_flip(&t->tree); +} + + //________________________________________________________________________________________________________________________ // @@ -77,6 +88,14 @@ bool su2_tensor_is_consistent(const struct su2_tensor* t); void su2_tensor_fmove(const struct su2_tensor* restrict t, const int i_ax, struct su2_tensor* restrict r); +//________________________________________________________________________________________________________________________ +// + +// contraction + +void su2_tensor_contract_simple(const struct su2_tensor* restrict s, const int* restrict i_ax_s, const struct su2_tensor* restrict t, const int* restrict i_ax_t, const int ndim_mult, struct su2_tensor* restrict r); + + //________________________________________________________________________________________________________________________ // diff --git a/src/tensor/su2_tree.c b/src/tensor/su2_tree.c index a2f418a..1034cbe 100644 --- a/src/tensor/su2_tree.c +++ b/src/tensor/su2_tree.c @@ -32,6 +32,51 @@ void delete_su2_irreducible_list(struct su2_irreducible_list* list) } +//________________________________________________________________________________________________________________________ +/// +/// \brief Whether two irreducible 'j' quantum number lists are logically equal. +/// +bool su2_irreducible_list_equal(const struct su2_irreducible_list* s, const struct su2_irreducible_list* t) +{ + if (s->num != t->num) { + return false; + } + + for (int k = 0; k < s->num; k++) { + if (s->jlist[k] != t->jlist[k]) { + return false; + } + } + + return true; +} + + +//________________________________________________________________________________________________________________________ +/// +/// \brief Comparison function used by 'qsort'. +/// +static int compare_su2_irreducible_lists(const void* a, const void* b) +{ + const struct su2_irreducible_list* x = a; + const struct su2_irreducible_list* y = b; + + assert(x->num == y->num); + for (int i = 0; i < x->num; i++) + { + if (x->jlist[i] < y->jlist[i]) { + return -1; + } + else if (x->jlist[i] > y->jlist[i]) { + return 1; + } + } + + // lists are equal + return 0; +} + + //________________________________________________________________________________________________________________________ /// /// \brief Allocate a charge sector array. @@ -55,6 +100,48 @@ void delete_charge_sectors(struct charge_sectors* sectors) } +//________________________________________________________________________________________________________________________ +/// +/// \brief Find the index of charge sector 'jlist', assuming that the list of charge sectors is sorted. +/// Returns -1 if charge sector cannot be found. +/// +long charge_sector_index(const struct charge_sectors* sectors, const qnumber* jlist) +{ + const struct su2_irreducible_list target = { + .jlist = (qnumber*)jlist, // cast to avoid compiler warning; we do not modify 'jlist' + .num = sectors->ndim, + }; + + // search interval: [lower, upper) + long lower = 0; + long upper = sectors->nsec; + while (true) + { + if (lower >= upper) { + return -1; + } + const long i = (lower + upper) / 2; + const struct su2_irreducible_list current = { + .jlist = §ors->jlists[i * sectors->ndim], + .num = sectors->ndim, + }; + const int c = compare_su2_irreducible_lists(&target, ¤t); + if (c < 0) { + // target < current + upper = i; + } + else if (c == 0) { + // target == current -> found it + return i; + } + else { + // target > current + lower = i + 1; + } + } +} + + //________________________________________________________________________________________________________________________ /// /// \brief Make a deep copy of an SU(2) symmetry tree. @@ -109,7 +196,7 @@ void delete_su2_tree(struct su2_tree_node* tree) /// /// \brief Whether two SU(2) symmetry trees are logically equal. /// -bool su2_tree_equal(const struct su2_tree_node* s, const struct su2_tree_node* t) +bool su2_tree_equal(const struct su2_tree_node* restrict s, const struct su2_tree_node* restrict t) { if (s->i_ax != t->i_ax) { return false; @@ -131,6 +218,28 @@ bool su2_tree_equal(const struct su2_tree_node* s, const struct su2_tree_node* t } +//________________________________________________________________________________________________________________________ +/// +/// \brief Whether two SU(2) symmetry trees have the same topology. +/// +bool su2_tree_equal_topology(const struct su2_tree_node* restrict s, const struct su2_tree_node* restrict t) +{ + // test whether both nodes are leaves + const bool leaf_s = su2_tree_node_is_leaf(s); + const bool leaf_t = su2_tree_node_is_leaf(t); + if (leaf_s != leaf_t) { + return false; + } + if (leaf_s) { + // both nodes are leaves + return true; + } + + // test whether subtrees have the same topology + return su2_tree_equal_topology(s->c[0], t->c[0]) && su2_tree_equal_topology(s->c[1], t->c[1]); +} + + //________________________________________________________________________________________________________________________ /// /// \brief Whether the SU(2) symmetry tree contains a leaf with axis index 'i_ax'. @@ -193,20 +302,176 @@ int su2_tree_num_nodes(const struct su2_tree_node* tree) //________________________________________________________________________________________________________________________ /// -/// \brief Fill the indicator for the axes indexed by the tree. +/// \brief Fill the list of axis indices contained in the tree, and return the number of axes. +/// +int su2_tree_axes_list(const struct su2_tree_node* tree, int* list) +{ + if (tree == NULL) { + return 0; + } + + // in-order traversal + int n0 = su2_tree_axes_list(tree->c[0], list); + list[n0] = tree->i_ax; + int n1 = su2_tree_axes_list(tree->c[1], list + n0 + 1); + + return n0 + 1 + n1; +} + + +//________________________________________________________________________________________________________________________ +/// +/// \brief Fill the indicator with the axes indexed by the tree. /// -void su2_tree_axes(const struct su2_tree_node* tree, bool* indicator) +void su2_tree_axes_indicator(const struct su2_tree_node* tree, bool* indicator) { if (tree == NULL) { return; } + assert(tree->i_ax >= 0); // must not be assigned yet assert(!indicator[tree->i_ax]); indicator[tree->i_ax] = true; - su2_tree_axes(tree->c[0], indicator); - su2_tree_axes(tree->c[1], indicator); + su2_tree_axes_indicator(tree->c[0], indicator); + su2_tree_axes_indicator(tree->c[1], indicator); +} + + +//________________________________________________________________________________________________________________________ +/// +/// \brief Update the axes indices using the provided map. +/// +void su2_tree_update_axes_indices(struct su2_tree_node* tree, const int* axis_map) +{ + if (tree == NULL) { + return; + } + + assert(tree->i_ax >= 0); + tree->i_ax = axis_map[tree->i_ax]; + + su2_tree_update_axes_indices(tree->c[0], axis_map); + su2_tree_update_axes_indices(tree->c[1], axis_map); +} + + +//________________________________________________________________________________________________________________________ +/// +/// \brief Temporary data structure for finding a subtree with leaf axis indices exactly matching a given list. +/// +struct subtree_with_leaf_axes_data +{ + const struct su2_tree_node* subtree; //!< pointer to found subtree, or NULL if not yet found + bool all_leaves_in_list; //!< whether axis indices of all leaves are contained in provided list + int num_leaves; //!< number of leaves in current subtree (if 'all_leaves_in_list' is true) +}; + + +//________________________________________________________________________________________________________________________ +/// +/// \brief Internal recursive function for finding a subtree with leaf axis indices exactly matching the provided 'i_ax' list. +/// +/// Entries in 'i_ax' must be pairwise different. +/// +static struct subtree_with_leaf_axes_data su2_tree_search_subtree_with_leaf_axes(const struct su2_tree_node* tree, const int* i_ax, const int num_axes) +{ + assert(tree != NULL); + + if (su2_tree_node_is_leaf(tree)) + { + struct subtree_with_leaf_axes_data data = { 0 }; + + for (int i = 0; i < num_axes; i++) { + if (tree->i_ax == i_ax[i]) { + data.all_leaves_in_list = true; + break; + } + } + + if (data.all_leaves_in_list) { + data.num_leaves = 1; + if (num_axes == 1) { + // searching for a single axis, which matches current leaf + data.subtree = tree; + } + } + + return data; + } + + struct subtree_with_leaf_axes_data data0 = su2_tree_search_subtree_with_leaf_axes(tree->c[0], i_ax, num_axes); + struct subtree_with_leaf_axes_data data1 = su2_tree_search_subtree_with_leaf_axes(tree->c[1], i_ax, num_axes); + + // fast return if subtree has already been found + if (data0.subtree != NULL) { + return data0; + } + if (data1.subtree != NULL) { + return data1; + } + + struct subtree_with_leaf_axes_data data = { 0 }; + + data.all_leaves_in_list = (data0.all_leaves_in_list && data1.all_leaves_in_list); + if (data.all_leaves_in_list) + { + data.num_leaves = data0.num_leaves + data1.num_leaves; + if (data.num_leaves == num_axes) { + // found subtree + data.subtree = tree; + } + } + + return data; +} + + +//________________________________________________________________________________________________________________________ +/// +/// \brief Find the subtree with the specified leaf axis indices. +/// +/// Ordering of axis indices is not relevant. Returns NULL if subtree does not exist. +/// +const struct su2_tree_node* su2_subtree_with_leaf_axes(const struct su2_tree_node* tree, const int* i_ax, const int num_axes) +{ + return su2_tree_search_subtree_with_leaf_axes(tree, i_ax, num_axes).subtree; +} + + +//________________________________________________________________________________________________________________________ +/// +/// \brief Replace the proper subtree with axis index 'i_ax' by 'subtree', and return a pointer to the old subtree, or NULL if 'i_ax' cannot be found. +/// +struct su2_tree_node* su2_tree_replace_subtree(struct su2_tree_node* tree, const int i_ax, struct su2_tree_node* subtree) +{ + if (tree == NULL) { + return NULL; + } + if (su2_tree_node_is_leaf(tree)) { + return NULL; + } + + if (tree->c[0]->i_ax == i_ax) { + struct su2_tree_node* old = tree->c[0]; + tree->c[0] = subtree; + return old; + } + + if (tree->c[1]->i_ax == i_ax) { + struct su2_tree_node* old = tree->c[1]; + tree->c[1] = subtree; + return old; + } + + // search in left subtree + struct su2_tree_node* old = su2_tree_replace_subtree(tree->c[0], i_ax, subtree); + if (old != NULL) { + return old; + } + // search in right subtree + return su2_tree_replace_subtree(tree->c[1], i_ax, subtree); } @@ -487,10 +752,15 @@ bool su2_fuse_split_tree_is_consistent(const struct su2_fuse_split_tree* tree) } const int i_ax_shared = tree->tree_fuse->i_ax; + // fusion and splitting trees cannot both be single leaves + if (su2_tree_node_is_leaf(tree->tree_fuse) && su2_tree_node_is_leaf(tree->tree_split)) { + return false; + } + bool* indicator_fuse = ct_calloc(tree->ndim, sizeof(bool)); bool* indicator_split = ct_calloc(tree->ndim, sizeof(bool)); - su2_tree_axes(tree->tree_fuse, indicator_fuse); - su2_tree_axes(tree->tree_split, indicator_split); + su2_tree_axes_indicator(tree->tree_fuse, indicator_fuse); + su2_tree_axes_indicator(tree->tree_split, indicator_split); for (int i = 0; i < tree->ndim; i++) { if (i == i_ax_shared) { @@ -531,6 +801,29 @@ bool su2_fuse_split_tree_equal(const struct su2_fuse_split_tree* s, const struct } +//________________________________________________________________________________________________________________________ +/// +/// \brief Flip the fuse and split trees. +/// +void su2_fuse_split_tree_flip(struct su2_fuse_split_tree* tree) +{ + struct su2_tree_node* tmp = tree->tree_fuse; + tree->tree_fuse = tree->tree_split; + tree->tree_split = tmp; +} + + +//________________________________________________________________________________________________________________________ +/// +/// \brief Update the axes indices of a fuse and split tree using the provided map. +/// +void su2_fuse_split_tree_update_axes_indices(struct su2_fuse_split_tree* tree, const int* axis_map) +{ + su2_tree_update_axes_indices(tree->tree_fuse, axis_map); + su2_tree_update_axes_indices(tree->tree_split, axis_map); +} + + //________________________________________________________________________________________________________________________ /// /// \brief Evaluate the fuse and split tree by contracting both subtrees, interpreting the nodes as Clebsch-Gordan coefficients. @@ -552,31 +845,6 @@ double su2_fuse_split_tree_eval_clebsch_gordan(const struct su2_fuse_split_tree* } -//________________________________________________________________________________________________________________________ -/// -/// \brief Comparison function used by 'qsort'. -/// -static int compare_charge_sectors(const void* a, const void* b) -{ - const struct su2_irreducible_list* x = a; - const struct su2_irreducible_list* y = b; - - assert(x->num == y->num); - for (int i = 0; i < x->num; i++) - { - if (x->jlist[i] < y->jlist[i]) { - return -1; - } - else if (x->jlist[i] > y->jlist[i]) { - return 1; - } - } - - // lists are equal - return 0; -} - - //________________________________________________________________________________________________________________________ /// /// \brief Enumerate all "charge sectors" ('j' quantum number configurations) for a given fuse and split tree layout. @@ -589,8 +857,8 @@ void su2_fuse_split_tree_enumerate_charge_sectors(const struct su2_fuse_split_tr bool* indicator_fuse = ct_calloc(tree->ndim, sizeof(bool)); bool* indicator_split = ct_calloc(tree->ndim, sizeof(bool)); - su2_tree_axes(tree->tree_fuse, indicator_fuse); - su2_tree_axes(tree->tree_split, indicator_split); + su2_tree_axes_indicator(tree->tree_fuse, indicator_fuse); + su2_tree_axes_indicator(tree->tree_split, indicator_split); // consistency check for (int i = 0; i < tree->ndim; i++) { if (i == i_ax_shared) { @@ -635,7 +903,7 @@ void su2_fuse_split_tree_enumerate_charge_sectors(const struct su2_fuse_split_tr } // sort lexicographically - qsort(merged_sectors, c, sizeof(struct su2_irreducible_list), compare_charge_sectors); + qsort(merged_sectors, c, sizeof(struct su2_irreducible_list), compare_su2_irreducible_lists); // copy data into output array allocate_charge_sectors(c, tree->ndim, sectors); diff --git a/src/tensor/su2_tree.h b/src/tensor/su2_tree.h index 21e86f7..1026244 100644 --- a/src/tensor/su2_tree.h +++ b/src/tensor/su2_tree.h @@ -23,6 +23,8 @@ void copy_su2_irreducible_list(const struct su2_irreducible_list* src, struct su void delete_su2_irreducible_list(struct su2_irreducible_list* list); +bool su2_irreducible_list_equal(const struct su2_irreducible_list* s, const struct su2_irreducible_list* t); + //________________________________________________________________________________________________________________________ /// @@ -39,6 +41,8 @@ void allocate_charge_sectors(const long nsec, const int ndim, struct charge_sect void delete_charge_sectors(struct charge_sectors* sectors); +long charge_sector_index(const struct charge_sectors* sectors, const qnumber* jlist); + //________________________________________________________________________________________________________________________ /// @@ -46,7 +50,7 @@ void delete_charge_sectors(struct charge_sectors* sectors); /// struct su2_tree_node { - int i_ax; //!< logical axis index (logical or auxiliary for leaf nodes, otherwise internal) + int i_ax; //!< tensor axis index (logical or auxiliary for leaf nodes, otherwise internal) struct su2_tree_node* c[2]; //!< pointer to left and right child nodes; NULL for a leaf node }; @@ -74,7 +78,9 @@ static inline bool su2_tree_node_is_leaf(const struct su2_tree_node* node) } } -bool su2_tree_equal(const struct su2_tree_node* s, const struct su2_tree_node* t); +bool su2_tree_equal(const struct su2_tree_node* restrict s, const struct su2_tree_node* restrict t); + +bool su2_tree_equal_topology(const struct su2_tree_node* restrict s, const struct su2_tree_node* restrict t); bool su2_tree_contains_leaf(const struct su2_tree_node* tree, const int i_ax); @@ -82,7 +88,15 @@ const struct su2_tree_node* su2_tree_find_parent_node(const struct su2_tree_node int su2_tree_num_nodes(const struct su2_tree_node* tree); -void su2_tree_axes(const struct su2_tree_node* tree, bool* indicator); +int su2_tree_axes_list(const struct su2_tree_node* tree, int* list); + +void su2_tree_axes_indicator(const struct su2_tree_node* tree, bool* indicator); + +void su2_tree_update_axes_indices(struct su2_tree_node* tree, const int* axis_map); + +const struct su2_tree_node* su2_subtree_with_leaf_axes(const struct su2_tree_node* tree, const int* i_ax, const int num_axes); + +struct su2_tree_node* su2_tree_replace_subtree(struct su2_tree_node* tree, const int i_ax, struct su2_tree_node* subtree); void su2_tree_fmove_left(struct su2_tree_node* tree); @@ -112,6 +126,10 @@ bool su2_fuse_split_tree_is_consistent(const struct su2_fuse_split_tree* tree); bool su2_fuse_split_tree_equal(const struct su2_fuse_split_tree* s, const struct su2_fuse_split_tree* t); +void su2_fuse_split_tree_flip(struct su2_fuse_split_tree* tree); + +void su2_fuse_split_tree_update_axes_indices(struct su2_fuse_split_tree* tree, const int* axis_map); + double su2_fuse_split_tree_eval_clebsch_gordan(const struct su2_fuse_split_tree* tree, const qnumber* restrict jlist, const int* restrict im_leaves); void su2_fuse_split_tree_enumerate_charge_sectors(const struct su2_fuse_split_tree* tree, const struct su2_irreducible_list* leaf_ranges, struct charge_sectors* sectors); diff --git a/src/util/rng.h b/src/util/rng.h index f575736..adf90a5 100644 --- a/src/util/rng.h +++ b/src/util/rng.h @@ -23,7 +23,7 @@ void seed_rng_state(const uint64_t seed, struct rng_state* state); //________________________________________________________________________________________________________________________ // -// uniformly distributed integer random numbers. +// uniformly distributed integer random numbers uint32_t rand_uint32(struct rng_state* state); diff --git a/test/run_tests.c b/test/run_tests.c index 642b565..27a78fc 100644 --- a/test/run_tests.c +++ b/test/run_tests.c @@ -44,6 +44,7 @@ char* test_clebsch_gordan_coefficients(); char* test_su2_tree_enumerate_charge_sectors(); char* test_su2_fuse_split_tree_enumerate_charge_sectors(); char* test_su2_tensor_fmove(); +char* test_su2_tensor_contract_simple(); char* test_su2_to_dense_tensor(); char* test_mps_vdot(); char* test_mps_orthonormalize_qr(); @@ -121,6 +122,7 @@ int main() TEST_FUNCTION_ENTRY(test_su2_tree_enumerate_charge_sectors), TEST_FUNCTION_ENTRY(test_su2_fuse_split_tree_enumerate_charge_sectors), TEST_FUNCTION_ENTRY(test_su2_tensor_fmove), + TEST_FUNCTION_ENTRY(test_su2_tensor_contract_simple), TEST_FUNCTION_ENTRY(test_su2_to_dense_tensor), TEST_FUNCTION_ENTRY(test_mps_vdot), TEST_FUNCTION_ENTRY(test_mps_orthonormalize_qr), diff --git a/test/tensor/test_su2_tensor.c b/test/tensor/test_su2_tensor.c index e1bd6a7..a579b8a 100644 --- a/test/tensor/test_su2_tensor.c +++ b/test/tensor/test_su2_tensor.c @@ -162,6 +162,352 @@ char* test_su2_tensor_fmove() } +char* test_su2_tensor_contract_simple() +{ + struct rng_state rng_state; + seed_rng_state(42, &rng_state); + + // construct the 's' tensor + struct su2_tensor s; + { + const int ndim = 9; + + // construct the fuse and split tree + // + // 5 0 2 + // \ \ / + // \ \/ fuse + // \ /6 + // \/ + // | + // |7 + // | + // /\ + // 8/ \ split + // /\ \ + // / \ \ + // 3 1 4 + // + struct su2_tree_node j0 = { .i_ax = 0, .c = { NULL, NULL } }; + struct su2_tree_node j1 = { .i_ax = 1, .c = { NULL, NULL } }; + struct su2_tree_node j2 = { .i_ax = 2, .c = { NULL, NULL } }; + struct su2_tree_node j3 = { .i_ax = 3, .c = { NULL, NULL } }; + struct su2_tree_node j4 = { .i_ax = 4, .c = { NULL, NULL } }; + struct su2_tree_node j5 = { .i_ax = 5, .c = { NULL, NULL } }; + struct su2_tree_node j6 = { .i_ax = 6, .c = { &j0, &j2 } }; + struct su2_tree_node j8 = { .i_ax = 8, .c = { &j3, &j1 } }; + struct su2_tree_node j7f = { .i_ax = 7, .c = { &j5, &j6 } }; + struct su2_tree_node j7s = { .i_ax = 7, .c = { &j8, &j4 } }; + + struct su2_fuse_split_tree tree = { .tree_fuse = &j7f, .tree_split = &j7s, .ndim = ndim }; + + if (!su2_fuse_split_tree_is_consistent(&tree)) { + return "internal consistency check for the fuse and split tree failed"; + } + + // outer (logical and auxiliary) 'j' quantum numbers + qnumber j0list[] = { 1, 3 }; + qnumber j1list[] = { 3 }; + qnumber j2list[] = { 0, 2 }; + qnumber j3list[] = { 2 }; + qnumber j4list[] = { 0 }; // auxiliary + qnumber j5list[] = { 0 }; // auxiliary + const struct su2_irreducible_list outer_jlists[6] = { + { .jlist = j0list, .num = ARRLEN(j0list) }, + { .jlist = j1list, .num = ARRLEN(j1list) }, + { .jlist = j2list, .num = ARRLEN(j2list) }, + { .jlist = j3list, .num = ARRLEN(j3list) }, + { .jlist = j4list, .num = ARRLEN(j4list) }, + { .jlist = j5list, .num = ARRLEN(j5list) }, + }; + + // degeneracy dimensions, indexed by 'j' quantum numbers + // j: 0 1 2 3 + const long dim_degen0[4] = { 0, 4, 0, 2 }; + const long dim_degen1[4] = { 0, 0, 0, 3 }; + const long dim_degen2[3] = { 5, 0, 2 }; + const long dim_degen3[3] = { 0, 0, 7 }; + const long* dim_degen[4] = { + dim_degen0, + dim_degen1, + dim_degen2, + dim_degen3, + }; + + const int ndim_logical = 4; + const int ndim_auxiliary = 2; + + allocate_su2_tensor(CT_DOUBLE_COMPLEX, ndim_logical, ndim_auxiliary, &tree, outer_jlists, dim_degen, &s); + + if (!su2_tensor_is_consistent(&s)) { + return "internal consistency check for SU(2) tensor failed"; + } + if (s.charge_sectors.nsec == 0) { + return "expecting at least one charge sector in SU(2) tensor"; + } + + // fill degeneracy tensors with random entries + for (long c = 0; c < s.charge_sectors.nsec; c++) + { + // corresponding "degeneracy" tensor + struct dense_tensor* d = s.degensors[c]; + assert(d != NULL); + assert(d->dtype == s.dtype); + assert(d->ndim == s.ndim_logical); + + dense_tensor_fill_random_normal(numeric_one(d->dtype), numeric_zero(d->dtype), &rng_state, d); + } + } + + // construct the 't' tensor + struct su2_tensor t; + { + const int ndim = 9; + + // construct the fuse and split tree + // + // 3 2 5 0 + // \ / / / + // \/ / / + // 6\ / / fuse + // \/ / + // 8\ / + // \/ + // | + // |7 + // | + // /\ + // / \ split + // 1 4 + // + struct su2_tree_node j0 = { .i_ax = 0, .c = { NULL, NULL } }; + struct su2_tree_node j1 = { .i_ax = 1, .c = { NULL, NULL } }; + struct su2_tree_node j2 = { .i_ax = 2, .c = { NULL, NULL } }; + struct su2_tree_node j3 = { .i_ax = 3, .c = { NULL, NULL } }; + struct su2_tree_node j4 = { .i_ax = 4, .c = { NULL, NULL } }; + struct su2_tree_node j5 = { .i_ax = 5, .c = { NULL, NULL } }; + struct su2_tree_node j6 = { .i_ax = 6, .c = { &j3, &j2 } }; + struct su2_tree_node j8 = { .i_ax = 8, .c = { &j6, &j5 } }; + struct su2_tree_node j7f = { .i_ax = 7, .c = { &j8, &j0 } }; + struct su2_tree_node j7s = { .i_ax = 7, .c = { &j1, &j4 } }; + + struct su2_fuse_split_tree tree = { .tree_fuse = &j7f, .tree_split = &j7s, .ndim = ndim }; + + if (!su2_fuse_split_tree_is_consistent(&tree)) { + return "internal consistency check for the fuse and split tree failed"; + } + + // outer (logical and auxiliary) 'j' quantum numbers + qnumber j0list[] = { 0, 2 }; + qnumber j1list[] = { 3, 5 }; + qnumber j2list[] = { 3 }; + qnumber j3list[] = { 2 }; + qnumber j4list[] = { 4 }; + qnumber j5list[] = { 0 }; // auxiliary + const struct su2_irreducible_list outer_jlists[6] = { + { .jlist = j0list, .num = ARRLEN(j0list) }, + { .jlist = j1list, .num = ARRLEN(j1list) }, + { .jlist = j2list, .num = ARRLEN(j2list) }, + { .jlist = j3list, .num = ARRLEN(j3list) }, + { .jlist = j4list, .num = ARRLEN(j4list) }, + { .jlist = j5list, .num = ARRLEN(j5list) }, + }; + + // degeneracy dimensions, indexed by 'j' quantum numbers + // j: 0 1 2 3 4 5 + const long dim_degen0[3] = { 5, 0, 2 }; + const long dim_degen1[6] = { 0, 0, 0, 2, 0, 4 }; + const long dim_degen2[4] = { 0, 0, 0, 3 }; + const long dim_degen3[3] = { 0, 0, 7 }; + const long dim_degen4[5] = { 0, 0, 0, 0, 3 }; + const long* dim_degen[5] = { + dim_degen0, + dim_degen1, + dim_degen2, + dim_degen3, + dim_degen4, + }; + + const int ndim_logical = 5; + const int ndim_auxiliary = 1; + + allocate_su2_tensor(CT_DOUBLE_COMPLEX, ndim_logical, ndim_auxiliary, &tree, outer_jlists, dim_degen, &t); + + if (!su2_tensor_is_consistent(&t)) { + return "internal consistency check for SU(2) tensor failed"; + } + if (t.charge_sectors.nsec == 0) { + return "expecting at least one charge sector in SU(2) tensor"; + } + + // fill degeneracy tensors with random entries + for (long c = 0; c < t.charge_sectors.nsec; c++) + { + // corresponding "degeneracy" tensor + struct dense_tensor* d = t.degensors[c]; + assert(d != NULL); + assert(d->dtype == t.dtype); + assert(d->ndim == t.ndim_logical); + + dense_tensor_fill_random_normal(numeric_one(d->dtype), numeric_zero(d->dtype), &rng_state, d); + } + } + + for (int flip = 0; flip < 2; flip++) + { + if (flip == 1) { + su2_tensor_flip_trees(&s); + su2_tensor_flip_trees(&t); + } + + // perform contraction + const int i_ax_s[3] = { 3, 1, 4 }; + const int i_ax_t[3] = { 3, 2, 5 }; + struct su2_tensor r; + if (flip == 0) { + su2_tensor_contract_simple(&s, i_ax_s, &t, i_ax_t, 3, &r); + } + else { + su2_tensor_contract_simple(&t, i_ax_t, &s, i_ax_s, 3, &r); + } + + if (!su2_tensor_is_consistent(&r)) { + return "internal consistency check for SU(2) tensor failed"; + } + if (r.ndim_logical != 5) { + return "contracted SU(2) tensor does not have expected number of logical dimensions"; + } + if (r.ndim_auxiliary != 1) { + return "contracted SU(2) tensor does not have expected number of auxiliary dimensions"; + } + + // reference SU(2) tree of contracted tensor + if (flip == 0) + { + // construct the fuse and split tree + // + // 5 0 1 2 + // \ \ / / + // \ \/ / + // \ /6 / fuse + // \/ / + // 7\ / + // \/ + // | + // |8 + // | + // /\ + // / \ split + // 3 4 + // + struct su2_tree_node j0 = { .i_ax = 0, .c = { NULL, NULL } }; + struct su2_tree_node j1 = { .i_ax = 1, .c = { NULL, NULL } }; + struct su2_tree_node j2 = { .i_ax = 2, .c = { NULL, NULL } }; + struct su2_tree_node j3 = { .i_ax = 3, .c = { NULL, NULL } }; + struct su2_tree_node j4 = { .i_ax = 4, .c = { NULL, NULL } }; + struct su2_tree_node j5 = { .i_ax = 5, .c = { NULL, NULL } }; + struct su2_tree_node j6 = { .i_ax = 6, .c = { &j0, &j1 } }; + struct su2_tree_node j7 = { .i_ax = 7, .c = { &j5, &j6 } }; + struct su2_tree_node j8f = { .i_ax = 8, .c = { &j7, &j2 } }; + struct su2_tree_node j8s = { .i_ax = 8, .c = { &j3, &j4 } }; + + struct su2_fuse_split_tree tree_ref = { .tree_fuse = &j8f, .tree_split = &j8s, .ndim = 9 }; + if (!su2_fuse_split_tree_is_consistent(&tree_ref)) { + return "internal consistency check for the fuse and split tree failed"; + } + + if (!su2_fuse_split_tree_equal(&r.tree, &tree_ref)) { + return "fuse and split tree of contracted SU(2) tensor does not match reference"; + } + } + else + { + // construct the fuse and split tree + // + // 1 2 + // \ / fuse + // \/ + // | + // |6 + // | + // /\ + // 7/ \ + // /\ \ + // / \8 \ split + // / /\ \ + // / / \ \ + // 5 3 4 0 + // + struct su2_tree_node j0 = { .i_ax = 0, .c = { NULL, NULL } }; + struct su2_tree_node j1 = { .i_ax = 1, .c = { NULL, NULL } }; + struct su2_tree_node j2 = { .i_ax = 2, .c = { NULL, NULL } }; + struct su2_tree_node j3 = { .i_ax = 3, .c = { NULL, NULL } }; + struct su2_tree_node j4 = { .i_ax = 4, .c = { NULL, NULL } }; + struct su2_tree_node j5 = { .i_ax = 5, .c = { NULL, NULL } }; + struct su2_tree_node j8 = { .i_ax = 8, .c = { &j3, &j4 } }; + struct su2_tree_node j7 = { .i_ax = 7, .c = { &j5, &j8 } }; + struct su2_tree_node j6f = { .i_ax = 6, .c = { &j1, &j2 } }; + struct su2_tree_node j6s = { .i_ax = 6, .c = { &j7, &j0 } }; + + struct su2_fuse_split_tree tree_ref = { .tree_fuse = &j6f, .tree_split = &j6s, .ndim = 9 }; + if (!su2_fuse_split_tree_is_consistent(&tree_ref)) { + return "internal consistency check for the fuse and split tree failed"; + } + + if (!su2_fuse_split_tree_equal(&r.tree, &tree_ref)) { + return "fuse and split tree of contracted SU(2) tensor does not match reference"; + } + } + + // convert to full dense tensors + struct dense_tensor s_dns; + struct dense_tensor t_dns; + struct dense_tensor r_dns; + su2_to_dense_tensor(&s, &s_dns); + su2_to_dense_tensor(&t, &t_dns); + su2_to_dense_tensor(&r, &r_dns); + + if (dense_tensor_is_zero(&s_dns, 0.) || dense_tensor_is_zero(&t_dns, 0.)) { + return "to-be contracted SU(2) tensors should not be zero"; + } + + struct dense_tensor s_dns_perm; + struct dense_tensor t_dns_perm; + const int perm_s[4] = { 0, 2, 3, 1 }; + const int perm_t[5] = { 3, 2, 0, 1, 4 }; + transpose_dense_tensor(perm_s, &s_dns, &s_dns_perm); + transpose_dense_tensor(perm_t, &t_dns, &t_dns_perm); + + // contract dense tensors, as reference + struct dense_tensor r_dns_ref; + if (flip == 0) { + dense_tensor_dot(&s_dns_perm, TENSOR_AXIS_RANGE_TRAILING, &t_dns_perm, TENSOR_AXIS_RANGE_LEADING, 2, &r_dns_ref); + } + else { + dense_tensor_dot(&t_dns_perm, TENSOR_AXIS_RANGE_LEADING, &s_dns_perm, TENSOR_AXIS_RANGE_TRAILING, 2, &r_dns_ref); + } + + // compare + if (!dense_tensor_allclose(&r_dns, &r_dns_ref, 1e-13)) { + return "tensor resulting from contraction of two SU(2) tensors does not match reference"; + } + + delete_dense_tensor(&r_dns_ref); + delete_dense_tensor(&t_dns_perm); + delete_dense_tensor(&s_dns_perm); + delete_dense_tensor(&r_dns); + delete_dense_tensor(&t_dns); + delete_dense_tensor(&s_dns); + delete_su2_tensor(&r); + } + + delete_su2_tensor(&t); + delete_su2_tensor(&s); + + return 0; +} + + char* test_su2_to_dense_tensor() { hid_t file = H5Fopen("../test/tensor/data/test_su2_to_dense_tensor.hdf5", H5F_ACC_RDONLY, H5P_DEFAULT); @@ -207,11 +553,11 @@ char* test_su2_to_dense_tensor() qnumber j3list[] = { 2, 4 }; qnumber j4list[] = { 0 }; // auxiliary const struct su2_irreducible_list outer_jlists[5] = { - { .jlist = j0list, .num = sizeof(j0list) / sizeof(qnumber) }, - { .jlist = j1list, .num = sizeof(j1list) / sizeof(qnumber) }, - { .jlist = j2list, .num = sizeof(j2list) / sizeof(qnumber) }, - { .jlist = j3list, .num = sizeof(j3list) / sizeof(qnumber) }, - { .jlist = j4list, .num = sizeof(j4list) / sizeof(qnumber) }, + { .jlist = j0list, .num = ARRLEN(j0list) }, + { .jlist = j1list, .num = ARRLEN(j1list) }, + { .jlist = j2list, .num = ARRLEN(j2list) }, + { .jlist = j3list, .num = ARRLEN(j3list) }, + { .jlist = j4list, .num = ARRLEN(j4list) }, }; // degeneracy dimensions, indexed by 'j' quantum numbers @@ -242,7 +588,7 @@ char* test_su2_to_dense_tensor() // fill degeneracy tensors with random entries struct rng_state rng_state; - seed_rng_state(42, &rng_state); + seed_rng_state(43, &rng_state); for (long c = 0; c < t.charge_sectors.nsec; c++) { // corresponding "degeneracy" tensor diff --git a/test/tensor/test_su2_tree.c b/test/tensor/test_su2_tree.c index 2b1bec7..0621102 100644 --- a/test/tensor/test_su2_tree.c +++ b/test/tensor/test_su2_tree.c @@ -3,6 +3,9 @@ #include "util.h" +#define ARRLEN(a) (sizeof(a) / sizeof(a[0])) + + char* test_su2_tree_enumerate_charge_sectors() { hid_t file = H5Fopen("../test/tensor/data/test_su2_tree_enumerate_charge_sectors.hdf5", H5F_ACC_RDONLY, H5P_DEFAULT); @@ -41,11 +44,11 @@ char* test_su2_tree_enumerate_charge_sectors() qnumber j2list[] = { 3, 7 }; qnumber j4list[] = { 0, 4, 10 }; const struct su2_irreducible_list leaf_ranges[5] = { - { .jlist = j0list, .num = sizeof(j0list) / sizeof(qnumber) }, - { .jlist = j1list, .num = sizeof(j1list) / sizeof(qnumber) }, - { .jlist = j2list, .num = sizeof(j2list) / sizeof(qnumber) }, - { .jlist = NULL, .num = 0 }, // not used - { .jlist = j4list, .num = sizeof(j4list) / sizeof(qnumber) }, + { .jlist = j0list, .num = ARRLEN(j0list) }, + { .jlist = j1list, .num = ARRLEN(j1list) }, + { .jlist = j2list, .num = ARRLEN(j2list) }, + { .jlist = NULL, .num = 0 }, // not used + { .jlist = j4list, .num = ARRLEN(j4list) }, }; struct charge_sectors sectors; @@ -140,13 +143,13 @@ char* test_su2_fuse_split_tree_enumerate_charge_sectors() qnumber j5list[] = { 0, 2, 6 }; qnumber j6list[] = { 1, 5 }; const struct su2_irreducible_list leaf_ranges[7] = { - { .jlist = j0list, .num = sizeof(j0list) / sizeof(qnumber) }, - { .jlist = j1list, .num = sizeof(j1list) / sizeof(qnumber) }, - { .jlist = j2list, .num = sizeof(j2list) / sizeof(qnumber) }, - { .jlist = j3list, .num = sizeof(j3list) / sizeof(qnumber) }, - { .jlist = j4list, .num = sizeof(j4list) / sizeof(qnumber) }, - { .jlist = j5list, .num = sizeof(j5list) / sizeof(qnumber) }, - { .jlist = j6list, .num = sizeof(j6list) / sizeof(qnumber) }, + { .jlist = j0list, .num = ARRLEN(j0list) }, + { .jlist = j1list, .num = ARRLEN(j1list) }, + { .jlist = j2list, .num = ARRLEN(j2list) }, + { .jlist = j3list, .num = ARRLEN(j3list) }, + { .jlist = j4list, .num = ARRLEN(j4list) }, + { .jlist = j5list, .num = ARRLEN(j5list) }, + { .jlist = j6list, .num = ARRLEN(j6list) }, }; struct charge_sectors sectors;