Resident DF AO-to-MO transform (#805, #798)¶
The native molecular density-fitting method and its MP2 caller now use two checked index programs over the existing resident per-P provider. Each real FP64 panel is column-major. The left and right coefficient views share the AO frame, but carry independently bound orbital identities. Neither equal column counts nor a shared read-only coefficient owner identifies the orbital spaces. The method owns the whitened factor and the correlated caller owns the orbital selection. No change to fitting rank, reference, frozen core or scaling.
The old native expression evaluates left first. In ordered mode we keep that
association and sum each contracted AO axis in ascending order; we do not
promise bit identity with vendor GEMM. In reassociated mode we first contract
the smaller orbital extent (left on ties). For square AO panels this minimizes
both the first contraction’s work and its intermediate. There is no online
search or coefficient-value cache. Both orders can be forced for comparison.
A block spans at most 32 first-stage orbitals; the second stage accumulates
directly into row-major column p * n_right + q. A single first-stage orbital
uses two unthreaded BuiltIn GEMVs. The general leaf uses sequential Eigen
GEMM with explicitly supplied aligned A/B packs (block dimensions at most
64); older Eigen headers use the scalar fallback. This follows the layout and
packing distinctions of Matthews, sections 4-7
and GETT, sections 3-5. It does not add a
general GETT optimizer, BLAS dependency or complex adjoint convention.
The finite native interface charges the complete returned array, retained DF object/cache and numeric matrices, factor container capacity and panels, live coefficient allocations, O(P) tile descriptors, coefficient layout copies, per-worker intermediate, pack alignment and executor controls. Distinct input owners are counted once; the private result cannot alias them. No global AO copy, all-P numeric intermediate or four-index object is created by the transform. The ledger measures admitted storage, not process RSS. A cap below the one-worker requirement refuses before tile visits. With additional space, one OpenMP team may split independent P rows: its size is bounded by the cap, P and the requested worker count. Nested teams are suppressed. The MP2 and UMP2 callers propagate a finite method limit and take the larger of the method-phase and transform peaks. Both include their live orbital-energy vectors; beta includes the already live alpha output.
Measured transform cost¶
Single arm64 CPU, Eigen 5.0.1, C++17 release optimization, one worker and one BLAS thread. Times below are microseconds. Each row uses 21 rotated samples. Old is the predecessor including its output copy. Right is a layout-aware right-first Eigen reference writing the final product directly to the result. Auto includes new provider, domains, leases, planning, allocation, tile visits and execution on resident inputs, excluding Gaussian factor construction. The panel column is timed separately inside execution. No output conversion or coefficient packing occurs on these dense inputs.
P, AO, left, right |
Eigen old / right / auto |
BLAS build old / right / auto |
BLAS build first plan / warm plan / panel |
Tracked peak bytes |
|---|---|---|---|---|
5, 7, 2, 3 |
0.54 / 0.46 / 2.79 |
0.46 / 0.42 / 2.04 |
1.96 / 0.50 / 0.63 |
8,848 |
24, 31, 5, 19 |
25.38 / 44.92 / 26.63 |
29.54 / 37.58 / 27.63 |
1.67 / 0.79 / 24.42 |
231,704 |
32, 48, 31, 7 |
191.58 / 93.75 / 98.46 |
67.83 / 47.29 / 96.25 |
1.38 / 0.83 / 92.25 |
706,288 |
19, 63, 1, 47 |
18.71 / 240.13 / 21.21 |
17.42 / 77.54 / 18.54 |
1.83 / 0.75 / 16.17 |
641,272 |
13, 72, 35, 41 |
246.96 / 266.67 / 295.17 |
85.17 / 69.00 / 286.58 |
1.17 / 1.00 / 285.04 |
821,976 |
16, 128, 95, 21 |
1673.79 / 726.25 / 729.71 |
322.54 / 144.42 / 718.63 |
2.83 / 2.38 / 706.79 |
2,564,336 |
Every active run visits exactly P panels. For the last row the first complete
BLAS-build call took 717.13 us; median execution was 711.25 us. Forced left
and right medians were 1933.33 and 738.46 us. Each worker uses 87,584 bytes,
including 65,536 bytes of backend packs. The skinny case uses 1,048 bytes and
no GEMM packs. Maximum factor difference across both profiles was 1.14e-13.
In a separate strided-left run, the largest case packed 97,280 coefficient
bytes in 10.58 us (median); peak rose to 2,661,616 bytes.
Layout checks and planner setup cost more than the tiny contraction. The
largest case improves over Eigen’s left-first expression but loses to
Accelerate even at matched order. Opaque BLAS workspace is not admitted by
this leaf. Further backend qualification and crossovers remain open.
Reproduce with benchmarks/tensor_mo_transform.cpp, using the CPU ADR build
flags and a single thread; its CSV includes all orders, first call, planning,
packing, execution, panel time, workers and visits.
Numerical and workflow checks¶
Asymmetric sentinel factors, explicit double sums over two AO axes, unequal and selected spaces, signed strides, shared/expired owners, writer epochs, source and pattern revisions, zero/unavailable tiles, provider failure, overflow and exact/below-cap tests exercise the executor. The packed GEMM also runs with Eigen allocation disabled across a block edge and padded output. Sanitizers and a BLAS-enabled OpenMP build pass the four native cases, including nested-team suppression. The same four cases pass on Eigen 3.4.0. Five Python storage controls pass against a focused binding of the changed production method, including noncontiguous and empty columns and wrong-row refusal. The standard package test environment refused its incomplete basis overlay before collection and is not counted as a pass. A separately compiled Gaussian water fixture (def2-SVP, def2-universal-JKFIT, 24 AO and 113 auxiliary functions) reproduced independent raw-integral Cholesky whitening and left-first products in both 4-by-7 and 7-by-4 spaces to 6.67e-16. The public native call averaged 50.4 and 47.7 us respectively over 21 uses, with a 672,816-byte tracked peak.
A complete single-worker vq control used neutral singlet H4 at (0,0,0), (0,0,1.43), (0.2,3.8,0.15), (0.4,3.85,1.70) bohr, STO-3G, def2-SVP-RIFIT, zero frozen orbitals, exact restricted HF and unscaled DF-MP2. The installed native RHF supplied converged orbitals (energy/gradient targets 1e-13/1e-10); a separate executable compiled from the changed native MP2 and DF sources performed the affected method under a 64 MiB request. The comparison reuses the same reference for the predecessor. It does not qualify a packaged wheel.
Quantity (Hartree) |
New native result |
ORCA 6.1.1 |
|---|---|---|
HF |
-2.209069888186157 |
-2.20906988818442 |
Opposite spin |
-0.029472383917311 |
-0.029472384 |
Same spin |
-0.000150986978124 |
-0.000150987 |
Correlation |
-0.029623370895435 |
-0.029623371 |
Total |
-2.238693259081592 |
-2.238693259091 |
ORCA used VeryTightSCF, the matching def2-SVP/C auxiliary and no frozen core; an additional unit-opposite/zero-same-spin projection isolates its OS term. The ORCA SS entry is the difference of printed correlation energies. The native predecessor components agree within 3e-17 Eh. Independent Cholesky factors and explicit MP2 sums give maximum factor error 5.56e-17 and component error 4.34e-19 Eh. MP2 reported 104,064 workspace bytes. Runtime fingerprints and vq receipts are on #805/#798. No UMP2, CC, periodic, Python NumPy transform, or other method acceptance is inferred. Packaged binding integration and broader integral, SCF and tensor migrations remain open.