Bug for contractions with Z2 x U(1) symmetry

[manuschneider]

Contractions of several UniTensors with Z2 x U(1) symmetry fail:

bond_a = cytnx.Bond(
    cytnx.BD_IN, [cytnx.Qs(0,0)>>1],
    [cytnx.Symmetry.Zn(2), cytnx.Symmetry.U1()])
bond_b = bond_a.redirect()
uT = cytnx.UniTensor([bond_a, bond_b],dtype=cytnx.Type.Double)
[Sf,Uf,Vtf] = cytnx.linalg.Svd(uT,True)
Tfnew = cytnx.Contract([Uf,Sf,Vtf])

If I run this in python, the kernel crashes.

A sequential contraction as follows works fine:

Tnew = cytnx.Contract(Uf,Sf)
Tnew = cytnx.Contract([Tnew,Vtf])

[manuschneider]

The bug seems to be related to the optimal contraction order. The following code works:

Uf.set_name("U")
Sf.set_name("S")
Vtf.set_name("V")
Tnew = cytnx.Contract([Uf,Sf,Vtf], order = "(( U , S ), V)" , optimal = False)

But if I set optimal = True, the kernel crashes.

[yingjerkao]

What is the version of cytnx?

[manuschneider]

I tested it in the current master branch (happened before and after I merged the fermion branch), so commit fab7f48

[manuschneider]

This might be a more general issue of contractions (at least with more than one symmetry involved). Here is another example that fails:

symms = [cytnx.Symmetry.U1(), cytnx.Symmetry.U1()];
bond_a = cytnx.Bond(
    cytnx.BD_IN, [cytnx.Qs(0,1)>>2, cytnx.Qs(1,-1)>>2], symms)
bond_b = cytnx.Bond(
    cytnx.BD_OUT, [cytnx.Qs(0,1)>>4, cytnx.Qs(0,1)>>3], symms)
bond_c = bond_a.redirect()
bond_d = bond_b.redirect()

uT = cytnx.UniTensor([bond_a, bond_b, bond_c, bond_d],dtype=cytnx.Type.Double,labels=["a","b","c","d"])
uT.permute_([1, 0, 2, 3])
[S, U, Vt] = cytnx.linalg.Svd(uT)
Tnew=cytnx.Contract(U,S)

The last line gives an error:

# Cytnx error occur at void cytnx::linalg::iAdd(cytnx::Tensor&, const cytnx::Tensor&)
# error: [iAdd] error, the two tensor does not have the same shape. Lt rank: [3] Rt rank: [3] 

# file : [Cytnx/src/linalg/iAdd.cpp) (17)
Stack trace:
[Cytnx_lib/cytnx/cytnx.cpython-313-x86_64-linux-gnu.so)(+0xcf8d89) [0x7fa7f34f8d89]
[Cytnx_lib/cytnx/cytnx.cpython-313-x86_64-linux-gnu.so)(_ZN5cytnx6linalg4iAddERNS_6TensorERKS1_+0x137) [0x7fa7f3124017]
[Cytnx_lib/cytnx/cytnx.cpython-313-x86_64-linux-gnu.so)(_ZN5cytnx6TensorpLIS0_EERS0_RKT_+0x11) [0x7fa7f3112001]
[Cytnx_lib/cytnx/cytnx.cpython-313-x86_64-linux-gnu.so)(_ZN5cytnx14BlockUniTensor8contractERKN5boost13intrusive_ptrINS_14UniTensor_baseEEERKbS8_+0x1e08) [0x7fa7f30bc8b8]
[Cytnx_lib/cytnx/cytnx.cpython-313-x86_64-linux-gnu.so)(_ZN5cytnx8ContractERKNS_9UniTensorES2_RKbS4_+0xf2) [0x7fa7f3093692]

If I have only two legs (bond_a and bond_c), the kernel crashes instead. Without the permutation things work fine.

[manuschneider]

Also occurs with only one symmetry present:

symms = [cytnx.Symmetry.U1()]
bond_a = cytnx.Bond(cytnx.BD_IN, [cytnx.Qs(0)>>2, cytnx.Qs(1)>>2, cytnx.Qs(-1)>>2], symms)
bond_b = cytnx.Bond(cytnx.BD_IN, [cytnx.Qs(0)>>1, cytnx.Qs(1)>>4, cytnx.Qs(-1)>>3], symms)
bond_c = bond_a.redirect()
bond_d = bond_b.redirect()
Tf = cytnx.UniTensor([bond_a, bond_b, bond_c, bond_d],dtype=cytnx.Type.Double,
                        labels=["a","b","c","d"])
Tf.permute_([2, 3, 1, 0])
[Sf, Uf, Vtf] = cytnx.linalg.Svd(Tf)
Tnew=cytnx.Contract(Uf,Sf)

Gives the error

# Cytnx error occur at void cytnx::linalg::iAdd(cytnx::Tensor&, const cytnx::Tensor&)
# error: [iAdd] error, the two tensors do not have the same shape. Lt rank: [3] Rt rank: [3] 

[manuschneider]

The error persists in the current main branch.

[IvanaGyro]

Investigation: does #1027 (the #853 contraction-order fix) also fix this?

Short answer: no evidence that it does — none of the three repros posted in this thread currently reproduce a crash or error on current master in my test environment, with or without #1027's fix. Posting the investigation in case it's useful for whoever picks this back up.

Why I looked here

Contracts()/Contract([...], optimal=True) routes through Network::setOrder → RegularNetwork::getOptimalOrder → SearchTree::search_order → OptimalTreeSolver::solve — the exact function #1027 fixes (a heap-buffer-overflow when a connected component needs more than one sequential contraction). The first comment in this thread ("the bug seems to be related to the optimal contraction order... optimal=False with explicit order works, optimal=True crashes") pointed at exactly that code path, so it seemed worth checking directly.

Method

Built an ASan-instrumented Python extension at pre-#1027 master, and validated the setup first against a known case: a plain 4-tensor chain contraction through Contract([...]) reliably triggers AddressSanitizer's heap-buffer-overflow, with a full stack trace through OptimalTreeSolver::solve — confirming the harness would catch this bug if present.

Then ran all three repros posted in this issue against that same validated build:

1. The original repro (3-tensor Contract([Uf,Sf,Vtf]), Z2×U1, bond_a = cytnx.Bond(cytnx.BD_IN, [cytnx.Qs(0,0)>>1], [cytnx.Symmetry.Zn(2), cytnx.Symmetry.U1()])) — completes cleanly, no ASan report.
2. The set_name + explicit-order comparison from the first comment — both the optimal=False and optimal=True calls complete cleanly.
3. The U1×U1, 4-bond, permute_ + Contract(U, S) repro from a later comment — also completes cleanly (no iAdd shape-mismatch error reproduced).

Why the original repro can't hit #853's bug, structurally

#853/#1027's OOB requires an earlier merged node's index to land as the first operand of adjacencyMatrix[i] (an out-of-bounds vector index) in a later contraction round — which only happens once two merged nodes coexist simultaneously, requiring at least 4 tensors in one connected component. A 3-tensor U-S-V chain (repro 1) only ever produces one merged node, which is always appended last and is only ever read as the second operand (a safe in-range bitset bit-test, not a vector index) — it cannot structurally reach the OOB regardless of how it's triggered.

Repro 3 (Contract(U, S)) is only 2 tensors — the direct pairwise contract API, which doesn't go through SearchTree at all, so whatever it was hitting is unrelated to #853 by construction even if it reproduces elsewhere.

Open question

I don't have a build matching the original report's exact environment/cytnx version, and my tensors are zero-filled (vs. presumably non-degenerate content in the original use case) — either could matter. If someone can reproduce a crash on current master with one of these repros (or a fuller version of the original workflow), I'd be glad to check whether #1027 changes anything for that specific case.

---

Generated by Claude Code

[pcchen]

Posted by Claude Code on behalf of @pcchen

Cannot reproduce on current master — including with random-filled tensors

Following up on the earlier investigation in this thread: I ran every repro posted here against current master (aa62024), and additionally reran them with random-filled tensors plus a numeric reconstruction check, which addresses the open question about zero-filled tensors possibly masking the bug. All pass, and the contraction results are numerically correct, not just crash-free.

Environment: macOS 26.3.1 (Apple Silicon), Apple clang 21, openblas-apple preset with USE_HPTT=OFF, Python 3.14.

Part 1 — exact repros from this thread

Each repro runs in a separate subprocess so a hard kernel crash in one cannot mask the others.

Repro	Originally reported	Current master
Z2×U1, Contract([Uf,Sf,Vtf]) (issue body)	kernel crash	passes
Same with optimal=True vs optimal=False	crash with optimal=True	both pass
U1×U1, 4 bonds, permute_ + Svd + Contract(U,S)	iAdd shape-mismatch error	passes
Single U1 variant	iAdd shape-mismatch error	passes

repro613.py

"""Repros from issue #613, run in separate subprocesses so a hard crash
in one does not mask the others."""

import subprocess
import sys
import textwrap

REPROS = {
    "repro1-original-Z2xU1-3tensor-Contract-list": """
        import cytnx
        bond_a = cytnx.Bond(
            cytnx.BD_IN, [cytnx.Qs(0,0)>>1],
            [cytnx.Symmetry.Zn(2), cytnx.Symmetry.U1()])
        bond_b = bond_a.redirect()
        uT = cytnx.UniTensor([bond_a, bond_b], dtype=cytnx.Type.Double)
        [Sf, Uf, Vtf] = cytnx.linalg.Svd(uT, True)
        Tfnew = cytnx.Contract([Uf, Sf, Vtf])
        print("OK: shape", Tfnew.shape())
    """,
    "repro1b-optimal-True-vs-False": """
        import cytnx
        bond_a = cytnx.Bond(
            cytnx.BD_IN, [cytnx.Qs(0,0)>>1],
            [cytnx.Symmetry.Zn(2), cytnx.Symmetry.U1()])
        bond_b = bond_a.redirect()
        uT = cytnx.UniTensor([bond_a, bond_b], dtype=cytnx.Type.Double)
        [Sf, Uf, Vtf] = cytnx.linalg.Svd(uT, True)
        Uf.set_name("U"); Sf.set_name("S"); Vtf.set_name("V")
        t1 = cytnx.Contract([Uf, Sf, Vtf], order="(( U , S ), V)", optimal=False)
        print("OK optimal=False: shape", t1.shape())
        t2 = cytnx.Contract([Uf, Sf, Vtf], order="", optimal=True)
        print("OK optimal=True: shape", t2.shape())
    """,
    "repro2-U1xU1-4bond-permute-Contract-U-S": """
        import cytnx
        symms = [cytnx.Symmetry.U1(), cytnx.Symmetry.U1()]
        bond_a = cytnx.Bond(
            cytnx.BD_IN, [cytnx.Qs(0,1)>>2, cytnx.Qs(1,-1)>>2], symms)
        bond_b = cytnx.Bond(
            cytnx.BD_OUT, [cytnx.Qs(0,1)>>4, cytnx.Qs(0,1)>>3], symms)
        bond_c = bond_a.redirect()
        bond_d = bond_b.redirect()
        uT = cytnx.UniTensor([bond_a, bond_b, bond_c, bond_d],
                             dtype=cytnx.Type.Double, labels=["a","b","c","d"])
        uT.permute_([1, 0, 2, 3])
        [S, U, Vt] = cytnx.linalg.Svd(uT)
        Tnew = cytnx.Contract(U, S)
        print("OK: shape", Tnew.shape())
    """,
    "repro3-single-U1-permute-Contract-U-S": """
        import cytnx
        symms = [cytnx.Symmetry.U1()]
        bond_a = cytnx.Bond(cytnx.BD_IN,
            [cytnx.Qs(0)>>2, cytnx.Qs(1)>>2, cytnx.Qs(-1)>>2], symms)
        bond_b = cytnx.Bond(cytnx.BD_IN,
            [cytnx.Qs(0)>>1, cytnx.Qs(1)>>4, cytnx.Qs(-1)>>3], symms)
        bond_c = bond_a.redirect()
        bond_d = bond_b.redirect()
        Tf = cytnx.UniTensor([bond_a, bond_b, bond_c, bond_d],
                             dtype=cytnx.Type.Double, labels=["a","b","c","d"])
        Tf.permute_([2, 3, 1, 0])
        [Sf, Uf, Vtf] = cytnx.linalg.Svd(Tf)
        Tnew = cytnx.Contract(Uf, Sf)
        print("OK: shape", Tnew.shape())
    """,
}


def main() -> int:
    failures = 0
    for name, code in REPROS.items():
        print(f"=== {name} ===", flush=True)
        result = subprocess.run(
            [sys.executable, "-c", textwrap.dedent(code)],
            capture_output=True, text=True, timeout=300,
        )
        sys.stdout.write(result.stdout)
        if result.returncode != 0:
            failures += 1
            print(f"FAILED (exit {result.returncode})")
            sys.stdout.write(result.stderr[-2000:])
        print(flush=True)
    print(f"summary: {len(REPROS) - failures}/{len(REPROS)} passed")
    return 1 if failures else 0


if __name__ == "__main__":
    sys.exit(main())

=== repro1-original-Z2xU1-3tensor-Contract-list ===
OK: shape [1, 1]

=== repro1b-optimal-True-vs-False ===
OK optimal=False: shape [1, 1]
OK optimal=True: shape [1, 1]

=== repro2-U1xU1-4bond-permute-Contract-U-S ===
OK: shape [7, 4, 28]

=== repro3-single-U1-permute-Contract-U-S ===
OK: shape [6, 8, 48]

summary: 4/4 passed

Part 2 — random-filled tensors + numeric verification

The earlier investigation noted its tensors were zero-filled, leaving open whether nontrivial content could still trigger the bug. This part fills the same three structures with cytnx.random.uniform_ and verifies that Contract(Contract(U, S), Vt) reconstructs the original tensor. The diff is taken at the UniTensor level, so blocks are matched by quantum number rather than storage order.

repro613_random.py

"""Issue #613 repros with random-filled tensors plus numeric verification
that Contract results reconstruct the original tensor."""

import sys

import cytnx


def fill_random(ut, seed):
    cytnx.random.uniform_(ut, low=-1.0, high=1.0, seed=seed)


def reconstruction_error(original, s, u, vt):
    rec = cytnx.Contract(cytnx.Contract(u, s), vt)
    diff = original - rec.permute(original.labels())
    return max(float(b.Abs().Max().item()) for b in diff.get_blocks_())


def check(name, fn):
    try:
        fn()
        print(f"{name}: OK")
        return True
    except Exception as exc:
        print(f"{name}: FAILED - {type(exc).__name__}: {exc}")
        return False


def repro1_random():
    bond_a = cytnx.Bond(
        cytnx.BD_IN, [cytnx.Qs(0, 0) >> 1],
        [cytnx.Symmetry.Zn(2), cytnx.Symmetry.U1()])
    bond_b = bond_a.redirect()
    ut = cytnx.UniTensor([bond_a, bond_b], dtype=cytnx.Type.Double)
    fill_random(ut, seed=7)
    s, u, vt = cytnx.linalg.Svd(ut, True)
    combined = cytnx.Contract([u, s, vt])
    err = reconstruction_error(ut, s, u, vt)
    assert err < 1e-12, f"reconstruction error {err}"
    print(f"  3-tensor Contract shape {combined.shape()}, recon err {err:.2e}")


def repro2_random():
    symms = [cytnx.Symmetry.U1(), cytnx.Symmetry.U1()]
    bond_a = cytnx.Bond(cytnx.BD_IN, [cytnx.Qs(0, 1) >> 2, cytnx.Qs(1, -1) >> 2], symms)
    bond_b = cytnx.Bond(cytnx.BD_OUT, [cytnx.Qs(0, 1) >> 4, cytnx.Qs(0, 1) >> 3], symms)
    ut = cytnx.UniTensor(
        [bond_a, bond_b, bond_a.redirect(), bond_b.redirect()],
        dtype=cytnx.Type.Double, labels=["a", "b", "c", "d"])
    fill_random(ut, seed=11)
    ut.permute_([1, 0, 2, 3])
    s, u, vt = cytnx.linalg.Svd(ut)
    us = cytnx.Contract(u, s)
    err = reconstruction_error(ut, s, u, vt)
    assert err < 1e-10, f"reconstruction error {err}"
    print(f"  Contract(U,S) shape {us.shape()}, recon err {err:.2e}")


def repro3_random():
    symms = [cytnx.Symmetry.U1()]
    bond_a = cytnx.Bond(
        cytnx.BD_IN, [cytnx.Qs(0) >> 2, cytnx.Qs(1) >> 2, cytnx.Qs(-1) >> 2], symms)
    bond_b = cytnx.Bond(
        cytnx.BD_IN, [cytnx.Qs(0) >> 1, cytnx.Qs(1) >> 4, cytnx.Qs(-1) >> 3], symms)
    ut = cytnx.UniTensor(
        [bond_a, bond_b, bond_a.redirect(), bond_b.redirect()],
        dtype=cytnx.Type.Double, labels=["a", "b", "c", "d"])
    fill_random(ut, seed=13)
    ut.permute_([2, 3, 1, 0])
    s, u, vt = cytnx.linalg.Svd(ut)
    us = cytnx.Contract(u, s)
    err = reconstruction_error(ut, s, u, vt)
    assert err < 1e-10, f"reconstruction error {err}"
    print(f"  Contract(U,S) shape {us.shape()}, recon err {err:.2e}")


def main() -> int:
    results = [
        check("repro1-Z2xU1-random", repro1_random),
        check("repro2-U1xU1-random", repro2_random),
        check("repro3-U1-random", repro3_random),
    ]
    print(f"summary: {sum(results)}/{len(results)} passed")
    return 0 if all(results) else 1


if __name__ == "__main__":
    sys.exit(main())

  3-tensor Contract shape [1, 1], recon err 2.23e-308
repro1-Z2xU1-random: OK
  Contract(U,S) shape [7, 4, 28], recon err 2.66e-15
repro2-U1xU1-random: OK
  Contract(U,S) shape [6, 8, 48], recon err 9.05e-15
repro3-U1-random: OK
summary: 3/3 passed

(Repro 1's 2.23e-308 is a denormal — the 1×1 case reconstructs essentially exactly; the 4-bond cases reconstruct to normal double-precision roundoff.)

Conclusion

Combined with the earlier Linux/ASan check of the same repros in this thread, the evidence points to this having been fixed on master at some point after the "error persists" report of 2026-03-18 (plausibly by the block-contraction / contraction-order work, though I did not bisect). Caveat: this run is one platform (macOS arm64, OpenBLAS, no HPTT).

@manuschneider — could you retry on current master in your environment? If it no longer reproduces for you either, I'd suggest closing this as fixed — optionally after adding the Part 2 script as a pytests/ regression test so it stays fixed.

---

Generated by Claude Code

[IvanaGyro]

Investigate by Codex:

Verification result

You were right to ask for a real run. My previous guess (#1005) was too late and therefore wrong as the fixing PR for these repros.

After building and running the issue #613 examples across historical merge commits, I found:

• 2026-03-18 master / commit 92dcecb4 fails exactly as reported:

  • original Z2×U1 Contract([Uf, Sf, Vtf]) repro segfaults,

  • U1×U1 repro fails with the iAdd shape-mismatch,

  • single-U1 repro fails with the same iAdd shape-mismatch.

• PR Rsvd for symmetric tensors #744 / commit 928a4416 fixes only the original Z2×U1 crash, but the two later Contract(U, S) repros still fail.

• PR fixed reversed qnums infixed reversed qnums in linalg functions if combined incoming bond BdLeft is outgoing instead of incoming (#761 #763) linalg functions #764 / commit b9f3e12b is the first tested merge commit where all three examples pass.

So the corrected answer is:

PR #764 — “fixed reversed qnums in linalg functions if combined incoming bond BdLeft is outgoing instead of incoming” — is the first verified PR that fixes all examples from the issue comment.
Link: #764

There is a nuance:

• If you ask specifically about the first/original Z2×U1 Contract([Uf, Sf, Vtf]) crash, that was already fixed by PR Rsvd for symmetric tensors #744.

• If you ask about the full set of examples in the comment, including the two Contract(U, S) iAdd shape-mismatch repros, the first verified fixing PR is PR fixed reversed qnums infixed reversed qnums in linalg functions if combined incoming bond BdLeft is outgoing instead of incoming (#761 #763) linalg functions #764.

Why #764 fits the failure mode

The issue repros all go through linalg::Svd(...) before contraction. PR #764 changed how linalg code handles qnums when the combined left bond is outgoing. In current code, Svd.cpp builds combined left/right bonds, then chooses either normal qnums or reversed qnums depending on BdLeft.type(). That is the exact class of bug implied by the PR title and by the issue repros: bad auxiliary qnums from SVD cause later block contraction to accumulate incompatible block shapes.

The regression test added for this area explicitly covers a multi-sector SVD case where the combined left bond is BD_OUT, verifies SVD reconstruction, and checks that U, S, and vT agree on the auxiliary qnums and degeneracies.

Build/run matrix I verified

Commit / PR	Result
92dcecb — 2026-03-18 master	Fails 0/3: original repro segfaults; repro2/repro3 fail with iAdd shape mismatch
928a441 — PR #744	Partial pass 1/3: original repro passes; repro2/repro3 still fail
24e8cc6 — PR #760	Partial pass 1/3: original repro passes; repro2/repro3 still fail
b9f3e12 — PR #764	Passes 3/3
03b3901 — PR #830	Passes 3/3
cec15a0 — PR #859	Passes 3/3
a7e2175 — pre-#767 first-parent commit	Passes 3/3
dba9969 — PR #767	Passes 3/3
d08b5db — pre-#998 first-parent commit	Passes 3/3
246fa1a — PR #998	Passes 3/3

That confirms #998/#1005 are not the cause of the issue #613 examples becoming fixed; they are later unrelated fixes for shared-block metadata mutation.

Three Cases

Case 1: repro1-Z2xU1-list

import cytnx
bond_a = cytnx.Bond(cytnx.BD_IN, [cytnx.Qs(0,0)>>1], [cytnx.Symmetry.Zn(2), cytnx.Symmetry.U1()])
bond_b = bond_a.redirect()
uT = cytnx.UniTensor([bond_a, bond_b], dtype=cytnx.Type.Double)
[Sf, Uf, Vtf] = cytnx.linalg.Svd(uT, True)
Tfnew = cytnx.Contract([Uf, Sf, Vtf])
print("OK", Tfnew.shape())

Case 2: repro2-U1xU1-contract-US

import cytnx
symms = [cytnx.Symmetry.U1(), cytnx.Symmetry.U1()]
bond_a = cytnx.Bond(cytnx.BD_IN, [cytnx.Qs(0,1)>>2, cytnx.Qs(1,-1)>>2], symms)
bond_b = cytnx.Bond(cytnx.BD_OUT, [cytnx.Qs(0,1)>>4, cytnx.Qs(0,1)>>3], symms)
bond_c = bond_a.redirect(); bond_d = bond_b.redirect()
uT = cytnx.UniTensor([bond_a,bond_b,bond_c,bond_d], dtype=cytnx.Type.Double, labels=["a","b","c","d"])
uT.permute_([1,0,2,3])
[S,U,Vt] = cytnx.linalg.Svd(uT)
Tnew = cytnx.Contract(U,S)
print("OK", Tnew.shape())

Case 3: repro3-U1-contract-US

import cytnx
symms = [cytnx.Symmetry.U1()]
bond_a = cytnx.Bond(cytnx.BD_IN, [cytnx.Qs(0)>>2, cytnx.Qs(1)>>2, cytnx.Qs(-1)>>2], symms)
bond_b = cytnx.Bond(cytnx.BD_IN, [cytnx.Qs(0)>>1, cytnx.Qs(1)>>4, cytnx.Qs(-1)>>3], symms)
bond_c = bond_a.redirect(); bond_d = bond_b.redirect()
Tf = cytnx.UniTensor([bond_a,bond_b,bond_c,bond_d], dtype=cytnx.Type.Double, labels=["a","b","c","d"])
Tf.permute_([2,3,1,0])
[Sf,Uf,Vtf] = cytnx.linalg.Svd(Tf)
Tnew = cytnx.Contract(Uf,Sf)
print("OK", Tnew.shape())