Internal Function Documentation

is_fermionic(ops::NTuple{N,OpID}, n::Int, op_cache_vec::OpCacheVec) where {N} -> Bool

Compute the fermionic parity of the part of ops acting on sites >= n.

The result is true when an odd number of the operators in ops are marked fermionic in op_cache_vec, and false otherwise.

are_equal(op1::NTuple{N,OpID}, op2::NTuple{N,OpID}, n::Int) where {N} -> Bool

Check whether two products of OpID are equal from site n onward. The tuples must be sorted by decreasing site.

The comparison proceeds from left to right. If both tuples have already moved to sites < n at the same position, the remaining entries are ignored and the tuples are considered equal. Otherwise, entries must match exactly.

terms_eq_from(n::Int) -> Function

Return a comparator function that checks whether two products of OpID are equal from site n onward.

The returned function expects the same tuple ordering as are_equal: once both tuples have moved to sites < n at the same position, the remaining entries are ignored.

sort_fermion_perm!(ops::AbstractVector{<:OpID}, op_cache_vec::OpCacheVec) -> Int

Sort ops in place by site index and return the fermionic sign generated by the permutation. The sort is stable.

Whenever two fermionic operators are swapped past one another, the returned sign is multiplied by -1. Non-fermionic swaps do not affect the sign. Operators on the same site keep their original relative order.

determine_val_type(os::OpIDSum)

Return the numeric element type needed to represent the coefficients and cached operator matrices in os.

This is ComplexF64 if either the stored coefficients or any cached operator matrix has nonzero imaginary part; otherwise it is Float64.

for_equal_sites(f::Function, ops::AbstractVector{<:OpID}) -> Nothing

Call f(i, j) for each maximal contiguous block of operators in ops that act on the same site, where i and j denote the start and end of the contiguous block.

This assumes ops are already sorted by site, as happens after sorting terms inside OpIDSum.

rewrite_in_operator_basis!(os::OpIDSum, basis_op_cache_vec::OpCacheVec)

Rewrite each same-site operator product in os into a single operator drawn from basis_op_cache_vec.

For each contiguous block of operators acting on the same site, the product matrix is matched against the supplied basis up to an overall scalar factor. The term coefficient is updated by that factor, the block is replaced by the matched basis operator, and the remaining slots are zeroed. If a product is the zero matrix, the whole term is set to zero. An error is thrown if a nonzero product cannot be represented by a single basis operator. After rewriting, each nonzero term is resorted by site and os.op_cache_vec is replaced by basis_op_cache_vec.

op_sum_to_opID_sum(os::OpSum{C}, sites::Vector{<:Index}) -> OpIDSum

Convert an ITensor OpSum into an OpIDSum.

Each distinct local operator encountered on a site is assigned an integer id and cached in an OpCacheVec, with id 1 reserved for the identity. Terms are then added to the resulting OpIDSum using those cached identifiers. The storage width is the maximum number of local operators in any input term, with a minimum width of two to support the reinterpreted fixed-width layout.

check_os_for_errors(os::OpIDSum) -> Nothing

Validate structural assumptions for os.

This checks that the cached local operators are linearly independent on every site and that, within every term, operators are sorted by site, at most one operator acts on each site, all terms carry the same total QN flux, and every term has even fermion parity.

prepare_opID_sum!(os::OpIDSum, basis_op_cache_vec::Union{Nothing,OpCacheVec}) -> Nothing

Apply optional preprocessing to os before MPO construction.

Currently this rewrites terms into basis_op_cache_vec when a basis cache is provided and otherwise leaves os unchanged.

get_onsite_op(ops::NTuple{N,OpID{Ti}}, n::Int) where {N, Ti} -> Ti

Return the operator id in ops acting on site n.

If no operator acts on site n, this returns 1, the reserved identity operator id.

BipartiteGraph{L,R,C}

A weighted bipartite graph with typed left and right vertex payloads.

Left and right vertex ids are the positions of those payloads in left_vertices and right_vertices. Edges are stored as per-left-vertex adjacency lists: right_vertex_ids_from_left[lv_id][edge_id] gives the right-vertex id for one stored edge entry, and edge_weights_from_left[lv_id][edge_id] gives its corresponding weight. Duplicate right-vertex ids may appear in a single adjacency list.

Type Parameters

• L: The type of the left vertices.
• R: The type of the right vertices.
• C: The scalar edge weight type.

Fields:

• left_vertices: metadata stored for each left vertex.
• right_vertices: metadata stored for each right vertex.
• right_vertex_ids_from_left: adjacency list from each left vertex to the connected right-vertex ids.
• edge_weights_from_left: edge weights stored entry-by-entry in parallel with right_vertex_ids_from_left.

BipartiteGraphConnectedComponents

Connected-component data for a BipartiteGraph.

Components that do not contain at least one edge are not retained.

Fields:

• lvs_of_component: for each retained component, the global left-vertex ids contained in that component.
• position_of_rvs_in_component: for each global right-vertex id, its local right-vertex position within its retained component. Right vertices that are not in a retained component keep the sentinel value typemax(Int).
• rv_size_of_component: number of right vertices in each retained component.

left_size(g::BipartiteGraph) -> Int

Return the number of left vertices in g.

left_size(ccs::BipartiteGraphConnectedComponents, cc::Int)

Return the number of left vertices in the retained connected component cc.

right_size(g::BipartiteGraph) -> Int

Return the number of right vertices in g.

left_vertex(g::BipartiteGraph, lv_id::Integer)

Return the payload associated with the left vertex id lv_id.

right_vertex(g::BipartiteGraph, rv_id::Integer)

Return the payload associated with the right vertex id rv_id.

num_edges(g::BipartiteGraph) -> Int

Return the number of stored edge entries in g.

weighted_edge_iterator(g::BipartiteGraph, lv_id::Integer)

Return a lazy iterator over the edges from the left vertex id lv_id.

Each item is a (rv_id, weight) pair, where rv_id is a right vertex id and weight is the corresponding edge weight. The iterator views the graph's parallel adjacency and weight vectors directly rather than copying them.

clear_edges_from_left!(g::BipartiteGraph, lv_id::Integer) -> Nothing

Remove all stored edge entries incident on the left vertex id lv_id.

This clears both the adjacent right-vertex ids and their corresponding edge weights, and releases retained capacity in those per-vertex adjacency lists.

combine_duplicate_adjacent_right_vertices!(g::BipartiteGraph, eq::Function) -> Vector{Int}

Combine adjacent right vertices in g that compare equal under eq.

The right vertices are assumed to already be grouped so that duplicate payloads appear contiguously. The predicate eq is called on adjacent right-vertex payloads. The first vertex of each equal run is kept, later vertices in the run are removed, and every right-vertex id stored in the left adjacency lists is remapped to the surviving right-vertex id. Edge weights are left unchanged, so duplicate edge entries may remain after the remapping. The returned vector maps each original right-vertex id to its new position after compaction.

compute_connected_components(g::BipartiteGraph) -> BipartiteGraphConnectedComponents
compute_connected_components(g::BipartiteGraph, workspace::Vector{Int}) -> BipartiteGraphConnectedComponents

Compute the connected components of g, retaining only components that contain at least one edge.

The returned object records, for each retained component, which global left-vertex ids it contains and how global right-vertex ids map into local right-vertex positions within that component.

The two-argument method uses workspace as scratch storage and as the backing storage for position_of_rvs_in_component in the returned object. It must have length at least right_size(g), and it should not be mutated while the returned component data is still in use.

num_connected_components(ccs::BipartiteGraphConnectedComponents)

Return the number of retained connected components in ccs.

get_cc_matrix(g::BipartiteGraph, ccs::BipartiteGraphConnectedComponents, cc::Int; clear_edges=false)
    -> Tuple{SparseMatrixCSC,Vector{Int},Vector{Int}}

Extract the retained connected component cc as a sparse weighted adjacency matrix.

The returned tuple contains:

• the component adjacency matrix, whose rows are local left-vertex ids and whose columns are local right-vertex ids,
• left_map, mapping local row indices back to global left-vertex ids,
• right_map, mapping local column indices back to global right-vertex ids.

If clear_edges=true, the consumed adjacency lists are emptied from g as the matrix is assembled. Duplicate stored edge entries with the same local row and column are summed, matching the semantics of sparse(row, col, val).

_minimum_vertex_cover_from_matching(
  global_right_vertex_ids_from_left,
  left_map,
  position_of_rvs_in_component,
  matched_right_of_left,
  matched_left_of_right,
) -> Tuple{Vector{Int},Vector{Int}}

Construct a minimum vertex cover from a maximum matching for one connected component.

The matching vectors use component-local ids, with 0 denoting an unmatched vertex. This applies Konig's theorem: starting from unmatched left vertices, it traverses alternating paths through unmatched left-to-right edges and matched right-to-left edges. The cover is the unreached local left vertices together with the reached local right vertices. Returned ids are component-local and sorted in ascending order.

_hopcroft_karp_maximum_matching(
  global_right_vertex_ids_from_left,
  left_map,
  position_of_rvs_in_component,
  num_right,
) -> Tuple{Vector{Int},Vector{Int}}

Compute a maximum matching for one connected component using Hopcroft-Karp.

left_map maps component-local left ids to global left-vertex ids, and position_of_rvs_in_component maps global right-vertex ids to component-local right ids. The returned matched_right_of_left has one entry per local left vertex and stores the matched local right id, or 0 if unmatched. matched_left_of_right is the analogous local-right-to-local-left map. Duplicate stored edge entries are harmless and have the same effect as a single edge for matching purposes.

_hopcroft_karp_layered_bfs!(
  global_right_vertex_ids_from_left,
  left_map,
  position_of_rvs_in_component,
  matched_right_of_left,
  matched_left_of_right,
  dist,
  queue,
) -> Int

Build the Hopcroft-Karp BFS layers from all unmatched local left vertices.

dist is filled with the layer distance for reachable local left vertices and typemax(Int) for unreachable ones. queue is reused as scratch storage. The return value is the shortest augmenting-path distance measured by the next left-to-right step, or typemax(Int) when no augmenting path remains.

_hopcroft_karp_augment_from!(
  global_right_vertex_ids_from_left,
  left_map,
  position_of_rvs_in_component,
  start_lv_id,
  unmatched_distance,
  matched_right_of_left,
  matched_left_of_right,
  dist,
  next_edge_id,
  path_left,
  path_right,
) -> Bool

Search one layered augmenting path from the unmatched local left vertex start_lv_id.

The search is iterative to avoid recursion on large components. On success it mutates matched_right_of_left and matched_left_of_right in place to apply the augmentation and returns true; otherwise it returns false. next_edge_id, path_left, and path_right are scratch buffers reused across searches within one Hopcroft-Karp phase.

minimum_vertex_cover(g::BipartiteGraph, ccs::BipartiteGraphConnectedComponents, cc::Int)
    -> Tuple{Vector{Int},Vector{Int}}

Compute a minimum vertex cover for retained connected component cc of g.

The result is (local_left_ids, local_right_ids), where both vectors contain component-local 1-based ids. Use ccs.lvs_of_component[cc][local_left_id] to map a returned left id back to a global left-vertex id. To map a returned right id back to a global right-vertex id, invert or scan ccs.position_of_rvs_in_component.

The cover is computed by finding a maximum matching with Hopcroft-Karp and then applying Konig's theorem. Duplicate stored edge entries do not change the cover, and this function does not mutate g or ccs.

BlockSparseMatrix{C}

Dictionary-backed block-sparse matrix representation used for intermediate MPO tensor storage.

Keys are (left_link, right_link) pairs and values are dense local operator matrices for that block. The right link ids are local to the component being processed; at_site! later returns offsets that place component-local right-link ids into the full outgoing MPO bond.

_qn_position_map(ind::ITensors.QNIndex) -> Tuple{Vector{Int},Vector{Int}}

Map dense index positions of a QN index to block coordinates.

For each 1-based dense position in ind, the first returned vector stores the QN block number containing that position and the second returned vector stores the offset within that block. These lookup tables are used when translating dense link and site coordinates into BlockSparseTensor block coordinates.

_blockid(left_block, right_block, site_prime_block, site_block,
         num_right_blocks, num_site_prime_blocks, num_site_blocks) -> Int

Encode a 4-index block coordinate into a dense 1-based block id.

The id orders blocks lexicographically with left_block slowest and site_block fastest. It is only a temporary lookup key; _block_from_id performs the inverse conversion.

_block_from_id(block_id, num_right_blocks, num_site_prime_blocks, num_site_blocks)
    -> Block{4}

Decode the dense block id produced by _blockid into an ITensors Block{4}.

_fill_splitblocks!(blocks, block_values, offsets, block_sparse_matrices) -> Nothing

Fill split-block destination blocks and values from component block matrices.

This helper is used after all QN indices have been split into one-dimensional blocks. For every nonzero matrix entry, it writes the corresponding Block((left_link, right_link + offset, site_prime, site)) destination and the entry value into the preallocated blocks and block_values arrays. Zero entries are skipped.

_to_ITensor_splitblocks(offsets, block_sparse_matrices, llink, rlink, site) -> ITensor

Construct a QN-conserving MPO tensor using split one-dimensional QN blocks.

All three input indices are first passed through ITensors.splitblocks. Each nonzero scalar entry in block_sparse_matrices then becomes one structural Block{4} entry in the output tensor. offsets[i] places component i into its slice of the outgoing right-link space, so local component right-link ids are shifted to global right-link ids by right_link + offsets[i].

to_ITensor(offsets, block_sparse_matrices, llink::QNIndex, rlink::QNIndex, site::QNIndex;
           splitblocks=false) -> ITensor

Assemble a block-sparse ITensor MPO tensor from component-local block matrices.

block_sparse_matrices[i] stores the local MPO matrix blocks for component i as (left_link, right_link) => local_operator_matrix. The corresponding offsets[i] shifts those component-local right-link ids into the global outgoing bond space. The returned tensor has indices (dag(llink), rlink, prime(site), dag(site)).

When splitblocks=true, construction delegates to _to_ITensor_splitblocks. Otherwise the method discovers the QN blocks touched by nonzero entries, allocates only those structural blocks, zero-initializes their storage, and writes nonzero entries into the appropriate block views.

to_ITensor(offsets, block_sparse_matrices, llink::Index, rlink::Index, site::Index;
           splitblocks=false) -> ITensor

Assemble a dense-index ITensor MPO tensor from component-local block matrices.

This fallback is used for non-QN indices. It allocates a dense array with axes (llink, rlink, prime(site), site) and writes each local operator matrix into the slice selected by (left_link, right_link + offsets[i], :, :). The splitblocks keyword is accepted for interface compatibility but has no effect for dense indices.

add_to_local_matrix!(a, weight, local_op, needs_JW_string) -> Nothing

Accumulate a weighted local operator matrix into a.

If needs_JW_string is true, the contribution is multiplied by the diagonal Jordan-Wigner sign pattern expected for 2-state or 4-state fermionic sites. a is mutated in place and local_op is read without copying.

CoSorterElement{T1,T2}

Pair-like element used by CoSorter so that sorting one array carries along a second array with the same permutation.

x is the sortable value and y is the carried value.

CoSorter{T1,T2,S,C}

Lightweight view that exposes two arrays as a single sortable collection, ordering by sortarray and applying the same swaps to coarray.

This is used when sorting OpIDSum term storage while keeping the corresponding scalar coefficients aligned with the same permutation.

sparse_qr(A::SparseMatrixCSC, tol::Real, absolute_tol::Bool)
    -> Tuple{Q,R,prow,pcol,rank}

Compute a sparse QR factorization of A using SuiteSparse SPQR.

If absolute_tol is false, tol is interpreted relative to SPQR's default tolerance scale. The return values are the orthogonal factor Q, upper-triangular factor R, row and column permutations, and the numerical rank. SPQR uses the current BLAS thread count for its internal thread setting.

for_non_zeros_batch(f::Function, A::SparseMatrixCSC, max_col::Int) -> Nothing

Iterate over the nonzero entries of the first max_col columns of A, calling f(values, rows, col) once per nonempty column.

values and rows are views into the CSC storage of A for that column, so the callback must not retain them beyond the call.

for_non_zeros_batch(f::Function, Q::SparseArrays.SPQR.QRSparseQ, max_col::Int) -> Nothing

Iterate over the first max_col columns of a sparse SPQR Q factor, forming each column explicitly and calling f(column, col). The dense column workspace is reused between calls, so the callback must copy it if it needs to keep it.

at_site!(ValType, g, n, sites, tol, absolute_tol, op_cache_vec, alg;
         combine_qn_sectors, output_level=0)
    -> Tuple{MPOGraph,Vector{Int},Vector{BlockSparseMatrix{ValType}},Index}

Process one site of the MPO construction algorithm.

Equivalent adjacent right vertices are compacted, connected components are computed, and each component is processed by either the sparse-QR path (alg == "QR") or the minimum-vertex-cover path (alg == "VC"). The returned tuple contains:

• the graph for site n + 1,
• offsets locating each connected component inside the outgoing bond space,
• the block-sparse local MPO tensors for each component,
• the outgoing link Index, optionally grouped into merged QN sectors.

tol and absolute_tol are used only by the QR path. combine_qn_sectors merges adjacent outgoing link sectors with the same QN after component ranks are known.

MPOGraph{N,C,Ti}

Type alias for the bipartite graph representation used during MPO construction.

Left vertices store LeftVertex metadata for the current site. Right vertices store fixed-width tuples of OpIDs describing the remaining operator content of a term, and edge weights carry scalar coefficients or decomposition weights.

pretty_print(g::MPOGraph, n::Int, op_cache_vec::OpCacheVec) -> Nothing

Print a human-readable summary of an MPO construction graph at site n.

The output includes left-vertex metadata, remaining right-vertex operator strings after site n, and weighted adjacency lists. This is intended for debugging and writes directly to standard output.

LeftVertex

Represents a left vertex of the bipartite MPO graph.

Each left vertex describes one incoming MPO link and the local operator that must be emitted at the current site while extending paths through the graph.

Fields:

• link: integer label for the incoming (left) MPO link.
• op_id: identifier specifying the operator acting on the current site.
• needs_JW_string: whether a Jordan-Wigner string must be inserted when connecting through this vertex.

add_to_next_graph!(next_graph, cur_graph, op_cache_vec, cur_site, cur_offset, next_edges) -> Nothing

Append the left vertices and adjacency lists described by next_edges to next_graph.

Each nonempty (bond_index, op_id) entry in next_edges creates one LeftVertex, reusing that entry's (right_vertex_ids, edge_weights) vectors as the outgoing adjacency list. The stored needs_JW_string flag is inferred from the fermionic parity of the connected right vertices. cur_offset shifts component-local bond indices into the full outgoing link.

find_first_eq_rv(g::MPOGraph, j::Int, n::Int) -> Int

Walk backward from right-vertex index j to the first earlier right vertex whose operator tuple is equivalent from site n onward.

Right vertices are expected to be ordered so equivalent suffixes are adjacent. This is used to merge equivalent right vertices after peeling off the operator acting on the current site.

merge_qn_sectors(qi_of_cc::Vector{Pair{QN,Int}})
    -> Tuple{Vector{Int},Vector{Pair{QN,Int}}}

Sort connected components by QN sector and merge adjacent entries with the same QN by summing their dimensions.

The first returned vector is the component order to use when laying out the outgoing bond space, and the second is the merged QN-sector description used to construct the outgoing link Index.

build_next_edges_specialization!(next_edges, g, cur_site, right_vertex_ids, edge_weights) -> Nothing

Fast path for building outgoing edges when a connected component has only a single left vertex.

For each current edge, this extracts the operator acting on cur_site + 1, finds the unique right vertex from cur_site + 2 onward, and stores the resulting id/weight entry in next_edges. next_edges must have one row and one column for each operator available on the next site.

process_single_left_vertex_cc!(
  matrix_of_cc, rank_of_cc, next_edges_of_cc, g, ccs, cc, n, sites, op_cache_vec
) -> Nothing

Handle the common connected-component case with exactly one left vertex.

This fills the local MPO tensor contribution for the component, records rank 1, and either applies the terminal scaling on the last site or builds the outgoing edges for the next site. For nonterminal sites, the consumed left adjacency list is cleared after the outgoing edges are built.

process_vertex_cover!(
  matrix_of_cc, rank_of_cc, next_edges_of_cc, g, ccs, n, sites, op_cache_vec
) -> Nothing

Process every connected component using the minimum-vertex-cover specialization.

For each component, the cover columns are laid out as [left_cover; right_cover] in the local bond dimension. Left-cover columns emit single local operators with unit graph weight; right-cover columns accumulate uncovered edges with their stored edge weights. On the final site, only the local tensor blocks are needed; otherwise this also builds next_edges_of_cc for the graph at site n + 1.

process_qr(
  matrix_of_cc, rank_of_cc, next_edges_of_cc, g, ccs, n, sites, tol, absolute_tol, op_cache_vec
) -> Nothing

Process every connected component using the sparse-QR path.

Single-left-vertex components use process_single_left_vertex_cc!; larger components are extracted as sparse weighted adjacency matrices, decomposed with sparse_qr, and translated into local tensor blocks plus outgoing graph edges. On the final site, the scalar R factor is applied directly to the local blocks.

@time_if output_level min_output_level msg ex

Conditionally time and print the execution of ex depending on the current verbosity level.

If output_level > min_output_level, this expands to a timed evaluation using @time, prefixing the printed message with min_output_level levels of two-space indentation followed by msg. Otherwise, ex is evaluated without timing output.