fix(linalg): replace lower-enum-index dtype selection with Type.type_promote

[yingjerkao]

Summary

Follow-up to #982: that PR fixed Type_class::type_promote (mixed complex/real pairs now promote by precision, e.g. ComplexFloat x Double -> ComplexDouble) and converted Vectordot. This PR converts the remaining nine linalg entry points that still hand-rolled "pick the lower enum index" dtype selection and therefore diverged from Add/Sub/Mul/Div and Vectordot:

Dot (Mat-Vec path), Matmul, Matmul_dg, Tensordot (cuTensor path), Axpy/Axpy_, Gemm/Gemm_, Gemm_Batch, Ger, and Mod.

All sites now follow the Vectordot reference pattern: promote via Type.type_promote, astype both operands (the promoted dtype can differ from both inputs; astype is a no-op when the dtype already matches), allocate the output at the promoted dtype, and dispatch the kernel by it.

Latent bugs fixed along the way

• Axpy/Ger floored the dtype with an inverted comparison (fin_dtype < Type.Float instead of >): complex inputs were clamped onto the Double kernel (all-ComplexDouble Axpy threw "Cannot cast complex128 to real"), and integer inputs indexed past the size-5 kernel tables — a reproducible segfault. Integer/bool now floors to Double after promotion.
• Gemm (out-of-place) was missing the integer→Double floor that Gemm_ has; integer inputs segfaulted the same way.
• Ger built its implicit alpha = 1 scalar before the floor could change fin_dtype; it is now created at the final dtype.
• Mod's buffer and kernel disagreed: ModInternalImpl already computes in Type_class::type_promote_t<TL,TR>, but the output buffer was allocated with the old min-index rule, so e.g. Uint64 % Int32 wrote Int64 data into a Uint64-labeled buffer. All buffer-dtype sites (Tensor-Tensor, Scalar, and the typed-scalar specializations) now use Type.type_promote.

Tests

• New tests/linalg_test/DtypePromotion_test.cpp (20 tests) covering the discriminating pairs (ComplexFloat x Double -> ComplexDouble, Uint64 x Int32 -> Int64) and the integer floors. All were observed failing before the fix: 16 assertion failures plus 3 segfaults. The Gemm_Batch case is #ifdef UNI_MKL-guarded like the existing suite; without MKL it still pins the promotion via the pre-dispatch cast of c plus the expected throw.
• Full suite on macOS/OpenBLAS: 1087 passed / 11 skipped / 4 failed — identical to baseline (the 4 are the known pre-existing BkUt_*Svd_truncate_return_err* failures addressed by fix(linalg): report only actually-dropped values in block SVD return_err #977). No existing test asserted the old lower-index dtype, so none needed updating.

Dependencies

Stacked on #982 (fix/type-promotion-cross-kind) — only the last commit (efa641fb) is new here. Merge #982 first; this PR retargets to master automatically when that branch is deleted on merge.

🤖 Generated with Claude Code

[gemini-code-assist]

Code Review

This pull request refactors several linear algebra operations (including Axpy, Dot, Gemm, Ger, Matmul, Mod, and Tensordot) to use proper dtype promotion via Type.type_promote instead of hand-rolled selection logic. It also ensures integer inputs to BLAS-only operations floor to Double and adds a comprehensive test suite for mixed-dtype promotion. The review feedback suggests optimizing Dot, Matmul_dg, and Mod by passing false to out.Init to skip redundant zero-initialization, since the output tensors are fully overwritten by their respective kernels.

Important

The consumer version of Gemini Code Assist on GitHub is being sunset. Starting June 18, 2026, new organization installations will be blocked, and all code review activity will officially cease on July 17, 2026.
For more details on the timeline and next steps, please review the Help Documentation.

[gemini-code-assist]

Since the Matvec_ii and cuMatvec_ii kernels fully overwrite the output tensor, we can skip the redundant zero-initialization by passing false as the fourth argument to out.Init. This matches the optimization pattern used in Matmul and Tensordot.

Suggested change

	out.Init(newshape, out_dtype, Tl.device());
	out.Init(newshape, out_dtype, Tl.device(), false);

[yingjerkao]

Applied in d711fbf. Verified the Matvec path fully overwrites out before adopting it: the integer kernel self-zeros each element (out[n] = 0; out[n] += ...) and the BLAS/cuBLAS gemv variants use beta = 0, so out is written without being read. DtypePromotion.Dot_matvec + the Dot suites pass under ASAN (which fills fresh allocations with 0xbe, so this isn't relying on coincidentally-zero memory).

[gemini-code-assist]

Since the Matmul_dg_ii and cuMatmul_dg_ii kernels fully overwrite the output tensor, we can skip the redundant zero-initialization by passing false as the fourth argument to out.Init. This matches the optimization pattern used in Matmul and Tensordot.

Suggested change

	out.Init({Tl.shape()[0], Tr.shape().back()}, out_dtype, Tl.device());
	out.Init({Tl.shape()[0], Tr.shape().back()}, out_dtype, Tl.device(), false);

[yingjerkao]

Not applying this one — it would corrupt results. Unlike Matmul (gemm, beta=0) and Tensordot, the CPU Matmul_dg diag_L path accumulates into out via axpy (_out[i*Nr] += _inl[i]*_inr[i*Nr], in Matmul_dg_internal.cpp), so it reads the existing out values and therefore requires the zero-initialization. (Only the diag_R driver and the GPU kernel overwrite directly.) Since one out.Init feeds both paths, dropping the zero-fill here would leave the diag_L BLAS result adding onto uninitialized memory. Leaving init_zero=true; confirmed DtypePromotion.Matmul_dg_complexfloat_double still passes.

[gemini-code-assist]

Since the element-wise modulo kernels fully overwrite the output tensor, we can skip the redundant zero-initialization by passing false as the fourth argument to out.Init. This matches the optimization pattern used in Matmul and Tensordot.

Suggested change

	if (Lt.shape().size() == 1 && Lt.shape()[0] == 1) {
	out.Init(Rt.shape(), Lt.dtype() < Rt.dtype() ? Lt.dtype() : Rt.dtype(), Lt.device());
	out.Init(Rt.shape(), Type.type_promote(Lt.dtype(), Rt.dtype()), Lt.device());
	icnst = true;
	} else if (Rt.shape().size() == 1 && Rt.shape()[0] == 1) {
	out.Init(Lt.shape(), Lt.dtype() < Rt.dtype() ? Lt.dtype() : Rt.dtype(), Lt.device());
	out.Init(Lt.shape(), Type.type_promote(Lt.dtype(), Rt.dtype()), Lt.device());
	icnst = true;
	} else {
	cytnx_error_msg(Lt.shape() != Rt.shape(),
	"[Mod] The two tensors do not have the same shape.%s", "\n");
	out.Init(Lt.shape(), Lt.dtype() < Rt.dtype() ? Lt.dtype() : Rt.dtype(), Lt.device());
	out.Init(Lt.shape(), Type.type_promote(Lt.dtype(), Rt.dtype()), Lt.device());
	}
	if (Lt.shape().size() == 1 && Lt.shape()[0] == 1) {
	out.Init(Rt.shape(), Type.type_promote(Lt.dtype(), Rt.dtype()), Lt.device(), false);
	icnst = true;
	} else if (Rt.shape().size() == 1 && Rt.shape()[0] == 1) {
	out.Init(Lt.shape(), Type.type_promote(Lt.dtype(), Rt.dtype()), Lt.device(), false);
	icnst = true;
	} else {
	cytnx_error_msg(Lt.shape() != Rt.shape(),
	"[Mod] The two tensors do not have the same shape.%s", "\n");
	out.Init(Lt.shape(), Type.type_promote(Lt.dtype(), Rt.dtype()), Lt.device(), false);
	}

[yingjerkao]

Applied in d711fbf (all three Init calls). ModInternalImpl (contiguous and non-contiguous) and cuMod assign every output element (_out[i] = ...), so out is fully overwritten. All DtypePromotion Mod cases pass under ASAN.

[ianmccul]

Ok, this is the follow-up that makes the type_promote change much safer for the old dtype-indexed linalg backends.

The key pattern here looks right to me: compute the promoted dtype, cast both operands to that dtype, allocate output at that dtype, then dispatch the old single-dtype kernel table. That avoids the dangerous case where only the output dtype changes while the backend still reinterprets input/output buffers according to the old enum-order rule.

I would still be careful about the claimed scope. This appears to cover the listed Tensor linalg entry points, but not all arithmetic paths. In particular, general Tensor in-place arithmetic (iAdd/iSub/iMul/iDiv) and UniTensor elementwise arithmetic still look like separate paths with old “choose a buffer, then mutate it” semantics. For example, linalg::Add(UniTensor, UniTensor) still appears to choose which operand to clone by enum ordering and then calls Add_ in-place.

I think the only safe way to do this change is to implement #941, followed by making dtype an opaque type (detect and fix every instance where dtype is computed by hand), and only then it might be safe to change type_promote.

Python code tends to do arithmetc on dtype as well, so the corresponding change to make dtype opaque in python (or switch to the numpy dtypes) should probably be done at the same time; see ianmccul/Cytnx-fixes#3

[ianmccul]

There are still float/complex promotion bugs at the tip of PR #984, even if we ignore integer arithmetic completely.

The clearest one is Outer.

Tensor out(new_shape, Type.type_promote(Tl.dtype(), Tr.dtype()), Tl.device());
lii.Outer_ii[Tl.dtype()][Tr.dtype()](out.storage, Tl.storage, Tr.storage, ...);

For ComplexFloat x Double, the front-end now allocates a ComplexDouble output, but the generated CPU kernel still dispatches to:

Outer_general<cytnx_complex64, cytnx_complex64, cytnx_double>(out, Lin, Rin, j1, j2);

So the output storage has ComplexDouble layout, but the kernel writes ComplexFloat values into it. That is a direct buffer-layout mismatch, not just a precision issue. The same dispatch pattern exists on GPU.

There are several other non-integer promotion paths still using ordinal dtype comparison rather than Type.type_promote:

• Lstsq.cpp still chooses A.dtype() < b.dtype() ? A.dtype() : b.dtype(), so ComplexFloat mixed with Double is cast to ComplexFloat, not ComplexDouble.
• Concatenate.cpp, Hstack.cpp, Vstack.cpp, and Directsum.cpp still use lower-enum dtype selection. Again, ComplexFloat + Double becomes ComplexFloat.
• Scalar binary arithmetic still uses enum order, for example this->dtype() < rhs.dtype(), so Scalar(ComplexFloat) + Scalar(Double) keeps/promotes to ComplexFloat.
• The UniTensor out-of-place arithmetic wrappers in Add.cpp, Sub.cpp, Mul.cpp, and Div.cpp still choose the output dtype by enum order before calling the in-place implementation.
• The existing Exp / Expf issue is still present: the output dtype is changed, then the kernel dispatch is selected from the output dtype while the input storage still has the original dtype. For example, Exp(Float) reads a float* buffer as double*, and Exp(ComplexFloat) reads complex<float>* as complex<double>*.

they were found with an AI search, there are probably others.

[codecov]

Codecov Report

❌ Patch coverage is 33.33333% with 112 lines in your changes missing coverage. Please review.
✅ Project coverage is 32.63%. Comparing base (d6dcd16) to head (d711fbf).
✅ All tests successful. No failed tests found.

Files with missing lines	Patch %	Lines
src/linalg/Mod.cpp	25.53%	2 Missing and 33 partials ⚠️
src/linalg/Sub.cpp	10.00%	3 Missing and 6 partials ⚠️
src/linalg/Div.cpp	20.00%	0 Missing and 8 partials ⚠️
src/linalg/Gemm.cpp	58.82%	0 Missing and 7 partials ⚠️
src/linalg/Add.cpp	37.50%	0 Missing and 5 partials ⚠️
src/linalg/Axpy.cpp	58.33%	0 Missing and 5 partials ⚠️
src/linalg/Mul.cpp	37.50%	0 Missing and 5 partials ⚠️
src/linalg/Outer.cpp	16.66%	0 Missing and 5 partials ⚠️
src/linalg/Directsum.cpp	0.00%	0 Missing and 4 partials ⚠️
src/linalg/Dot.cpp	20.00%	0 Missing and 4 partials ⚠️
... and 9 more

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #984      +/-   ##
==========================================
+ Coverage   31.82%   32.63%   +0.80%     
==========================================
  Files         230      230              
  Lines       33124    33021     -103     
  Branches    13852    13758      -94     
==========================================
+ Hits        10543    10775     +232     
+ Misses      15539    14972     -567     
- Partials     7042     7274     +232     

Flag	Coverage Δ
cpp	32.27% <33.33%> (+0.81%)	⬆️
python	61.84% <ø> (ø)

Flags with carried forward coverage won't be shown. Click here to find out more.

Components	Coverage Δ
C++ backend	33.48% <33.33%> (+0.93%)	⬆️
Python bindings	24.02% <ø> (ø)
Python package	61.84% <ø> (ø)

Files with missing lines	Coverage Δ
include/backend/Scalar.hpp	22.05% <100.00%> (-0.18%)	⬇️
src/linalg/Gemm_Batch.cpp	49.46% <100.00%> (+4.30%)	⬆️
src/linalg/Tensordot.cpp	29.09% <ø> (ø)
src/algo/Hstack.cpp	50.00% <33.33%> (+1.21%)	⬆️
src/algo/Vstack.cpp	43.58% <33.33%> (+1.08%)	⬆️
src/linalg/Exp.cpp	17.64% <0.00%> (-5.89%)	⬇️
src/linalg/Expf.cpp	0.00% <0.00%> (ø)
src/linalg/Lstsq.cpp	15.55% <0.00%> (+15.55%)	⬆️
src/algo/Concatenate.cpp	14.28% <0.00%> (-5.72%)	⬇️
src/linalg/Directsum.cpp	48.38% <0.00%> (+2.93%)	⬆️
... and 12 more

... and 6 files with indirect coverage changes

---

Continue to review full report in Codecov by Harness.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update d6dcd16...d711fbf. Read the comment docs.

🚀 New features to boost your workflow:

• 📦 JS Bundle Analysis: Save yourself from yourself by tracking and limiting bundle sizes in JS merges.

[yingjerkao]

Thanks @ianmccul — pushed 731294f addressing the concrete promotion bugs in the second comment.

Buffer-layout / reinterpret fixes

• Outer (CPU + GPU): now casts both operands to Type.type_promote(...) and dispatches by the promoted dtype, so it uses only the same-type diagonal kernels — no more writing complex64 into a complex128 allocation. (This also keeps Outer off the mixed-cuComplex kernel path.)
• Exp/Expf: the out-of-place versions fed Tin's original storage to the out.dtype()-keyed kernel (reading float* as double*). They now feed the already-cast out storage as the input; the kernels support in == out (as Exp_.cpp already relies on).

Lower-enum-index → Type.type_promote

• Lstsq, Concatenate, Hstack, Vstack, Directsum.
• Scalar binary arithmetic (radd/rmul/rsub/rdiv, both the Scalar and typed-scalar overloads).
• The UniTensor out-of-place Add/Sub/Mul/Div wrappers and their scalar variants (they chose the cloned operand / output dtype by enum order).

Tests: added DtypePromotion cases for every fixed entry point (the ComplexFloat × Double → ComplexDouble discriminator, plus Exp(Float) reading the correct buffer width), and updated the Hstack/Vstack/Directsum test helpers (CPU and GPU copies) that had hard-coded the old min-enum-index contract. Full CPU suite: 1101 passed / 11 skipped / 0 failed.

Tridiag also uses the min-index idiom but rejects complex inputs, so min ≡ type_promote there — left as-is.

On the broader point (first comment): you're right that this is whack-a-mole, and I hit concrete evidence for it — the GPU library doesn't currently compile on this branch: #982's type_promote change broke the cuKron device kernel (cuKron_internal.cuh:49, out[x] = Lin[x] * Rin[y] no longer type-checks because the promoted TO differs from the operands). I confirmed it reproduces on the branch tip with my changes stashed, so it's pre-existing and independent of this PR — but it means the "only the output dtype changed, the kernel still reinterprets by the old rule" hazard is already live on GPU, exactly as you described. The in-place iAdd/iSub/iMul/iDiv and UniTensor in-place paths you mention are also still separate. I don't think this PR should try to be the fix for all of that; the opaque-dtype / #941 direction is the real remedy.

I've filed the cuKron breakage as a separate must-fix: #999 (with root cause — Kron.cpp allocates the output via runtime Type.type_promote but derives the kernel TO from the compile-time type_promote_from_gpu_pointer_t trait, and #982 left the two inconsistent). It's also concrete evidence for doing #941 before further type_promote-dependent work.

[ianmccul]

I've filed the cuKron breakage as a separate must-fix: #999 (with root cause — Kron.cpp allocates the output via runtime Type.type_promote but derives the kernel TO from the compile-time type_promote_from_gpu_pointer_t trait, and #982 left the two inconsistent).

That cannot be the case, for sure type_promote_from_gpu_pointer_t uses the same type_promote mechanism, that was the whole point of the design: to localize the type promotion choice in one place. @pcchen already identified the true cause and a fix in #995 , the cuXXX have their own weird type promotion rules, cuFloatComplex * double does not promote properly, and there is no assignment operator for cuDoubleComplex = cuFloatComplex. GPU Kron itself is actually a good type-safe design, it is only the weird cuda types that caused it to break.

[pcchen]

Code Review

Replaces hand-rolled "pick the lower enum index" dtype selection with Type.type_promote across ~15 linalg/algo entry points (follow-up to the merged #982), fixing several real latent bugs along the way. Correct and well-tested.

Verified correct

I traced the headline fixes directly, and a full-breadth pass corroborated every other converted site.

• type_promote semantics — genuinely promotes across the real/complex boundary by precision (is_complex(L) != is_complex(R) → to_complex(promote(to_real, to_real)), e.g. ComplexFloat × Double → ComplexDouble), plus correct signed/unsigned handling. ✅
• Latent bugs fixed:
  • Axpy/Ger inverted floor — the old fin_dtype < Type.Float clamped complex inputs onto the Double kernel (the "Cannot cast complex128 to real" throw) and left integers to index past the size-5 kernel table (segfault). Now > Type.Float floors integers/bool — correct given the enum order (Float=4, ints ≥ 5). ✅
  • Gemm (out-of-place) — the previously-missing integer→Double floor is added, matching Gemm_. ✅
  • Ger alpha — now built after the floor, at the final dtype. ✅
  • Mod buffer/kernel disagreement — every buffer-dtype site (Tensor-Tensor, generic Scalar, all 10 typed-scalar specializations, and the reverse forms) now uses type_promote, matching what ModInternalImpl<TL,TR> computes (e.g. Uint64 % Int32 → Int64 buffer). ✅
• Every converted site casts both operands to the promoted dtype (no one-sided casts), allocates the output at it, and dispatches by it. Size-5 kernel tables (axpy/ger/Gemm) are only reached with floored/float dtypes; full tables (Matmul/Matvec/Outer/Mod) need no floor. Outer additionally fixes a 2-D-table-indexed-by-raw-dtype bug; Exp/Expf fix a separate reinterpret bug. ✅
• Base fix(Type): promote mixed complex/real dtypes by precision; add to_real/to_complex #982 is merged; CI is green on all platforms (BuildAndTest ×3 + DownstreamFindPackage), plus the 20 new DtypePromotion_test.cpp tests (all observed red pre-fix, incl. 3 segfaults).

Minor note (not a bug)

Gemm_Batch relies on per-operand dtype() > 4 error guards to reject integer inputs (they run before the dispatch, so the size-5 cuGemm_Batch_ii[promoted_dtype] is only reached with dtype ≤ 4 — no past-the-end access), rather than the floor-to-Double the other size-5 kernels use. That's a deliberate, safe asymmetry (batched GEMM doesn't support integer at all); a one-line comment noting it would help, but no change needed.

No findings requiring changes — a clean, high-value correctness fix that unifies dtype selection and squashes real segfaults/throws. LGTM.

Posted by Claude Code on behalf of @pcchen

[pcchen]

Approving. Verified across all ~15 converted sites: promotion correctly routed through Type.type_promote (real/complex by precision + signed/unsigned handling), both operands cast, outputs allocated at the promoted dtype, and the size-5 kernel tables only reached with floored/float dtypes. The four latent bugs are genuinely fixed — the Axpy/Ger inverted floor (complex mis-clamp + integer table-overrun segfault), the missing Gemm integer floor, Ger's alpha built after the floor, and Mod's buffer/kernel dtype disagreement. CI green on all platforms; 20 new promotion tests observed red pre-fix (incl. 3 segfaults).

Merge-order note: the Scalar.hpp hunks here are subsumed by #1011's variant rewrite — landing this first means #1011's rebase resolves that file trivially by taking its own version.

Posted by Claude Code on behalf of @pcchen