From 4ac2837b80c3126c89a79e8ee61a99f4177316e0 Mon Sep 17 00:00:00 2001 From: Ying-Jer Kao Date: Mon, 6 Jul 2026 10:44:08 +0800 Subject: [PATCH 1/2] fix(linalg): make GPU scalar in-place ops correct on CUDA (#988) Follow-up to #906/#980. Fixes the two GPU gaps that could not be addressed on a CUDA-less machine, plus two related memory-safety bugs found while validating on a CUDA device. - Regression: iMul/iDiv threw for a non-contiguous GPU LHS. Route a broadcast scalar through the kernel (the rconst path scales every physical element, so layout mappers are irrelevant). cuMul/cuDiv -- unlike cuAdd/cuSub -- ignore the mappers, so a genuine non-contiguous tensor*=tensor / tensor/=tensor now throws loudly instead of silently pairing mismatched elements. - Efficiency: the shape-{1} scalar wrapper now stays CPU-resident. For a length-1 RHS the GPU kernels read the scalar with a host-side dereference and pass it by value, so there is no per-call H2D copy (and the previous GPU-resident scalar made that host dereference UB). cudaSetDevice now targets Lt.device(), and the device check permits a host-resident length-1 scalar against a GPU LHS. - Safety: real op= complex on GPU reinterpreted the real output buffer as complex (out-of-bounds write). Guarded device-independently so it throws like the CPU path. - Safety: a LHS narrower than the promoted dtype (e.g. float_gpu *= 2.0) had the promoted-width result written into the narrower buffer (OOB; float += 1.0 yielded 0). Now computed in the promoted dtype and truncated back, matching CPU element-wise semantics. Adds tests/gpu Tensor.Gpu* covering the non-contiguous *=//= regression, host-scalar values, storage sharing, dtype preservation, real op= complex throwing, and the non-contiguous tensor-tensor guard (no underscores per #857). Verified on RTX 4070 Ti (CUDA 13): 7/7 new GPU tests, 8/8 CPU scalar tests, and the existing GPU iSub/iMul/iDiv all-types suites pass. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/Tensor.cpp | 110 +++++++++++++++++++------------------- src/linalg/iAdd.cpp | 33 ++++++++++-- src/linalg/iDiv.cpp | 62 ++++++++++++++++++--- src/linalg/iMul.cpp | 62 ++++++++++++++++++--- src/linalg/iSub.cpp | 33 ++++++++++-- tests/gpu/Tensor_test.cpp | 102 +++++++++++++++++++++++++++++++++++ 6 files changed, 327 insertions(+), 75 deletions(-) diff --git a/src/Tensor.cpp b/src/Tensor.cpp index 6a5fa7b38..90623662d 100644 --- a/src/Tensor.cpp +++ b/src/Tensor.cpp @@ -547,20 +547,20 @@ namespace cytnx { } namespace { - // Wrap a scalar as a shape-{1} Tensor on `device` so scalar in-place - // arithmetic reuses the iAdd/iSub/iMul/iDiv kernels (which broadcast a - // length-1 RHS and mutate LHS storage in place). See #906. + // Wrap a scalar as a host-resident shape-{1} Tensor so scalar in-place arithmetic reuses the + // iAdd/iSub/iMul/iDiv kernels (which broadcast a length-1 RHS and mutate LHS storage in + // place). See #906. The wrapper deliberately stays on the CPU even when the LHS is on a GPU: + // for a length-1 RHS the GPU kernels read the scalar with a host-side dereference and pass it + // in by value, so keeping it on the host avoids a per-call H2D copy. See #988. template - Tensor _scalar_as_rank1_tensor(const T &rc, const int device) { + Tensor _scalar_as_rank1_tensor(const T &rc) { Tensor s({1}, Type.cy_typeid(rc), Device.cpu); s.storage().at(0) = rc; - if (device != Device.cpu) s = s.to(device); return s; } - Tensor _scalar_as_rank1_tensor(const Scalar &rc, const int device) { + Tensor _scalar_as_rank1_tensor(const Scalar &rc) { Tensor s({1}, rc.dtype(), Device.cpu); s.item() = rc; // Sproxy assignment - if (device != Device.cpu) s = s.to(device); return s; } } // namespace @@ -579,62 +579,62 @@ namespace cytnx { } template <> Tensor &Tensor::operator+=(const cytnx_complex128 &rc) { - cytnx::linalg::iAdd(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iAdd(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator+=(const cytnx_complex64 &rc) { - cytnx::linalg::iAdd(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iAdd(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator+=(const cytnx_double &rc) { - cytnx::linalg::iAdd(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iAdd(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator+=(const cytnx_float &rc) { - cytnx::linalg::iAdd(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iAdd(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator+=(const cytnx_int64 &rc) { - cytnx::linalg::iAdd(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iAdd(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator+=(const cytnx_uint64 &rc) { - cytnx::linalg::iAdd(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iAdd(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator+=(const cytnx_int32 &rc) { - cytnx::linalg::iAdd(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iAdd(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator+=(const cytnx_uint32 &rc) { - cytnx::linalg::iAdd(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iAdd(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator+=(const cytnx_int16 &rc) { - cytnx::linalg::iAdd(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iAdd(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator+=(const cytnx_uint16 &rc) { - cytnx::linalg::iAdd(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iAdd(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator+=(const cytnx_bool &rc) { - cytnx::linalg::iAdd(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iAdd(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator+=(const Scalar &rc) { - cytnx::linalg::iAdd(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iAdd(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> @@ -654,62 +654,62 @@ namespace cytnx { } template <> Tensor &Tensor::operator-=(const cytnx_complex128 &rc) { - cytnx::linalg::iSub(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iSub(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator-=(const cytnx_complex64 &rc) { - cytnx::linalg::iSub(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iSub(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator-=(const cytnx_double &rc) { - cytnx::linalg::iSub(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iSub(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator-=(const cytnx_float &rc) { - cytnx::linalg::iSub(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iSub(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator-=(const cytnx_int64 &rc) { - cytnx::linalg::iSub(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iSub(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator-=(const cytnx_uint64 &rc) { - cytnx::linalg::iSub(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iSub(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator-=(const cytnx_int32 &rc) { - cytnx::linalg::iSub(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iSub(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator-=(const cytnx_uint32 &rc) { - cytnx::linalg::iSub(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iSub(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator-=(const cytnx_int16 &rc) { - cytnx::linalg::iSub(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iSub(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator-=(const cytnx_uint16 &rc) { - cytnx::linalg::iSub(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iSub(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator-=(const cytnx_bool &rc) { - cytnx::linalg::iSub(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iSub(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator-=(const Scalar &rc) { - cytnx::linalg::iSub(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iSub(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> @@ -729,62 +729,62 @@ namespace cytnx { } template <> Tensor &Tensor::operator*=(const cytnx_complex128 &rc) { - cytnx::linalg::iMul(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iMul(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator*=(const cytnx_complex64 &rc) { - cytnx::linalg::iMul(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iMul(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator*=(const cytnx_double &rc) { - cytnx::linalg::iMul(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iMul(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator*=(const cytnx_float &rc) { - cytnx::linalg::iMul(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iMul(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator*=(const cytnx_int64 &rc) { - cytnx::linalg::iMul(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iMul(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator*=(const cytnx_uint64 &rc) { - cytnx::linalg::iMul(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iMul(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator*=(const cytnx_int32 &rc) { - cytnx::linalg::iMul(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iMul(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator*=(const cytnx_uint32 &rc) { - cytnx::linalg::iMul(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iMul(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator*=(const cytnx_int16 &rc) { - cytnx::linalg::iMul(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iMul(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator*=(const cytnx_uint16 &rc) { - cytnx::linalg::iMul(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iMul(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator*=(const cytnx_bool &rc) { - cytnx::linalg::iMul(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iMul(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator*=(const Scalar &rc) { - cytnx::linalg::iMul(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iMul(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> @@ -805,62 +805,62 @@ namespace cytnx { } template <> Tensor &Tensor::operator/=(const cytnx_complex128 &rc) { - cytnx::linalg::iDiv(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iDiv(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator/=(const cytnx_complex64 &rc) { - cytnx::linalg::iDiv(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iDiv(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator/=(const cytnx_double &rc) { - cytnx::linalg::iDiv(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iDiv(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator/=(const cytnx_float &rc) { - cytnx::linalg::iDiv(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iDiv(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator/=(const cytnx_int64 &rc) { - cytnx::linalg::iDiv(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iDiv(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator/=(const cytnx_uint64 &rc) { - cytnx::linalg::iDiv(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iDiv(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator/=(const cytnx_int32 &rc) { - cytnx::linalg::iDiv(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iDiv(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator/=(const cytnx_uint32 &rc) { - cytnx::linalg::iDiv(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iDiv(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator/=(const cytnx_int16 &rc) { - cytnx::linalg::iDiv(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iDiv(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator/=(const cytnx_uint16 &rc) { - cytnx::linalg::iDiv(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iDiv(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator/=(const cytnx_bool &rc) { - cytnx::linalg::iDiv(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iDiv(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> Tensor &Tensor::operator/=(const Scalar &rc) { - cytnx::linalg::iDiv(*this, _scalar_as_rank1_tensor(rc, this->device())); + cytnx::linalg::iDiv(*this, _scalar_as_rank1_tensor(rc)); return *this; } template <> diff --git a/src/linalg/iAdd.cpp b/src/linalg/iAdd.cpp index 183163404..cf8c6ba31 100644 --- a/src/linalg/iAdd.cpp +++ b/src/linalg/iAdd.cpp @@ -11,8 +11,21 @@ namespace cytnx { namespace linalg { void iAdd(Tensor &Lt, const Tensor &Rt) { - cytnx_error_msg(Lt.device() != Rt.device(), + // A length-1 RHS that stays on the host is treated as a broadcast scalar: the GPU kernels + // read it with a host-side dereference and pass it into the kernel by value, so it needs + // neither a device match nor a per-call H2D copy of the scalar. See #988. + const bool rhs_is_host_scalar = + (Rt.device() == Device.cpu && Rt._impl->storage()._impl->size() == 1); + cytnx_error_msg(Lt.device() != Rt.device() && !rhs_is_host_scalar, "[iAdd] error, the two tensors have to be on the same device.%s", "\n"); + // In-place ops write the result back into the LHS storage, so a complex result cannot be + // stored in a real LHS. Guard here (device-independent) so the GPU path throws like the CPU + // path instead of reinterpreting the real output buffer as complex and corrupting it. See + // #988. + cytnx_error_msg(!Type.is_complex(Lt.dtype()) && Type.is_complex(Rt.dtype()), + "[iAdd] Cannot perform real += complex in-place: a complex result cannot be " + "stored in a real tensor.%s", + "\n"); if (!(Rt.shape().size() == 1 && Rt.shape()[0] == 1)) { cytnx_error_msg(Lt.shape() != Rt.shape(), "[iAdd] error, the two tensors do not have the same shape. Lt rank: [%d] " @@ -27,6 +40,20 @@ namespace cytnx { R = Rt; } + // GPU broadcast scalar with a LHS *narrower* than the promoted dtype (e.g. a Float tensor + // plus a Double scalar, or an integer tensor plus a fractional scalar): the in-place GPU + // kernels write the promoted-width result straight into the narrower LHS buffer, an + // out-of-bounds write that corrupts memory. Compute in the promoted dtype, then truncate + // back to the LHS dtype -- this matches the CPU element-wise semantics (compute in the + // promoted type, store into the LHS type). See #988. The CPU path already does this + // in place, and real += complex is rejected above. + if (rhs_is_host_scalar && Lt.device() != Device.cpu && Lt.dtype() > Rt.dtype()) { + Tensor promoted = Lt.astype(Type.type_promote(Lt.dtype(), Rt.dtype())); + iAdd(promoted, R); + Lt = promoted.astype(Lt.dtype()); + return; + } + static const std::vector empty_mapper; // if contiguous, then no need to calculate the mappers @@ -36,7 +63,7 @@ namespace cytnx { detail::DispatchInplaceArithmeticCPU<0>(Lt, R, empty_mapper, empty_mapper, empty_mapper); } else { #ifdef UNI_GPU - checkCudaErrors(cudaSetDevice(Rt.device())); + checkCudaErrors(cudaSetDevice(Lt.device())); Tensor tmpo; if (Lt.dtype() <= Rt.dtype()) tmpo = Lt; @@ -61,7 +88,7 @@ namespace cytnx { Rt._impl->invmapper()); } else { #ifdef UNI_GPU - checkCudaErrors(cudaSetDevice(Rt.device())); + checkCudaErrors(cudaSetDevice(Lt.device())); Tensor tmpo; if (Lt.dtype() <= Rt.dtype()) tmpo = Lt; diff --git a/src/linalg/iDiv.cpp b/src/linalg/iDiv.cpp index 89abe4b29..9ff1ba845 100644 --- a/src/linalg/iDiv.cpp +++ b/src/linalg/iDiv.cpp @@ -11,8 +11,21 @@ namespace cytnx { namespace linalg { void iDiv(Tensor &Lt, const Tensor &Rt) { - cytnx_error_msg(Lt.device() != Rt.device(), + // A length-1 RHS that stays on the host is treated as a broadcast scalar: the GPU kernels + // read it with a host-side dereference and pass it into the kernel by value, so it needs + // neither a device match nor a per-call H2D copy of the scalar. See #988. + const bool rhs_is_host_scalar = + (Rt.device() == Device.cpu && Rt._impl->storage()._impl->size() == 1); + cytnx_error_msg(Lt.device() != Rt.device() && !rhs_is_host_scalar, "[iDiv] The two tensors cannot be on different devices.%s", "\n"); + // In-place ops write the result back into the LHS storage, so a complex result cannot be + // stored in a real LHS. Guard here (device-independent) so the GPU path throws like the CPU + // path instead of reinterpreting the real output buffer as complex and corrupting it. See + // #988. + cytnx_error_msg(!Type.is_complex(Lt.dtype()) && Type.is_complex(Rt.dtype()), + "[iDiv] Cannot perform real /= complex in-place: a complex result cannot be " + "stored in a real tensor.%s", + "\n"); if (!(Rt.shape().size() == 1 && Rt.shape()[0] == 1)) { cytnx_error_msg(Lt.shape() != Rt.shape(), @@ -28,6 +41,20 @@ namespace cytnx { R = Rt; } + // GPU broadcast scalar with a LHS *narrower* than the promoted dtype (e.g. a Float tensor + // over a Double scalar, or an integer tensor over a fractional scalar): the in-place GPU + // kernels write the promoted-width result straight into the narrower LHS buffer, an + // out-of-bounds write that corrupts memory. Compute in the promoted dtype, then truncate + // back to the LHS dtype -- this matches the CPU element-wise semantics (compute in the + // promoted type, store into the LHS type). See #988. The CPU path already does this + // in place, and real /= complex is rejected above. + if (rhs_is_host_scalar && Lt.device() != Device.cpu && Lt.dtype() > Rt.dtype()) { + Tensor promoted = Lt.astype(Type.type_promote(Lt.dtype(), Rt.dtype())); + iDiv(promoted, R); + Lt = promoted.astype(Lt.dtype()); + return; + } + static const std::vector empty_mapper; // if contiguous, then no need to calculate the mappers if ((Lt.is_contiguous() && Rt.is_contiguous())) { @@ -36,7 +63,7 @@ namespace cytnx { detail::DispatchInplaceArithmeticCPU<3>(Lt, R, empty_mapper, empty_mapper, empty_mapper); } else { #ifdef UNI_GPU - checkCudaErrors(cudaSetDevice(Rt.device())); + checkCudaErrors(cudaSetDevice(Lt.device())); Tensor tmpo; if (Lt.dtype() <= Rt.dtype()) tmpo = Lt; @@ -61,11 +88,32 @@ namespace cytnx { Rt._impl->invmapper()); } else { #ifdef UNI_GPU - cytnx_error_msg(true, - "[iDiv][on GPU/CUDA] error two tensors must be contiguous. Call " - "Contiguous_() or Contiguous() first%s", - "\n"); - + if (R._impl->storage()._impl->size() == 1) { + // Broadcast scalar RHS: the rconst kernel path divides every element of the + // (possibly non-contiguous) LHS storage by the scalar in place, so the layout + // mappers are irrelevant and the result is correct. This is the #988 regression fix + // that lets e.g. `gpu_tensor.permute(...); gpu_tensor /= 2.0;` succeed again. + checkCudaErrors(cudaSetDevice(Lt.device())); + Tensor tmpo; + if (Lt.dtype() <= Rt.dtype()) + tmpo = Lt; + else + tmpo = Lt.clone(); + linalg_internal::lii.cuAri_ii[Lt.dtype()][Rt.dtype()]( + tmpo._impl->storage()._impl, Lt._impl->storage()._impl, R._impl->storage()._impl, + Lt._impl->storage()._impl->size(), Lt._impl->shape(), Lt._impl->invmapper(), + Rt._impl->invmapper(), 3); + if (Lt.dtype() > Rt.dtype()) Lt = tmpo; + } else { + // Genuine non-contiguous tensor/=tensor: unlike cuAdd/cuSub, the cuDiv GPU kernels + // ignore the layout mappers, so routing here would silently pair mismatched + // elements. Fail loudly instead of corrupting data. Call Contiguous_()/Contiguous() + // on the operands first. + cytnx_error_msg(true, + "[iDiv][on GPU/CUDA] non-contiguous tensor/=tensor is not supported. " + "Call Contiguous_() or Contiguous() on the operands first%s", + "\n"); + } #else cytnx_error_msg(true, "[Div] fatal error, the tensor is on GPU without CUDA support.%s", "\n"); diff --git a/src/linalg/iMul.cpp b/src/linalg/iMul.cpp index b54993ccf..b4d58e4d7 100644 --- a/src/linalg/iMul.cpp +++ b/src/linalg/iMul.cpp @@ -10,8 +10,21 @@ namespace cytnx { namespace linalg { void iMul(Tensor &Lt, const Tensor &Rt) { - cytnx_error_msg(Lt.device() != Rt.device(), + // A length-1 RHS that stays on the host is treated as a broadcast scalar: the GPU kernels + // read it with a host-side dereference and pass it into the kernel by value, so it needs + // neither a device match nor a per-call H2D copy of the scalar. See #988. + const bool rhs_is_host_scalar = + (Rt.device() == Device.cpu && Rt._impl->storage()._impl->size() == 1); + cytnx_error_msg(Lt.device() != Rt.device() && !rhs_is_host_scalar, "[iMul] The two tensors cannot be on different devices.%s", "\n"); + // In-place ops write the result back into the LHS storage, so a complex result cannot be + // stored in a real LHS. Guard here (device-independent) so the GPU path throws like the CPU + // path instead of reinterpreting the real output buffer as complex and corrupting it. See + // #988. + cytnx_error_msg(!Type.is_complex(Lt.dtype()) && Type.is_complex(Rt.dtype()), + "[iMul] Cannot perform real *= complex in-place: a complex result cannot be " + "stored in a real tensor.%s", + "\n"); if (!(Rt.shape().size() == 1 && Rt.shape()[0] == 1)) { cytnx_error_msg(Lt.shape() != Rt.shape(), "[iMul] The two tensors do not have the same shape. Lt rank: [%d] " @@ -28,6 +41,20 @@ namespace cytnx { R = Rt; } + // GPU broadcast scalar with a LHS *narrower* than the promoted dtype (e.g. a Float tensor + // times a Double scalar, or an integer tensor times a fractional scalar): the in-place GPU + // kernels write the promoted-width result straight into the narrower LHS buffer, an + // out-of-bounds write that corrupts memory. Compute in the promoted dtype, then truncate + // back to the LHS dtype -- this matches the CPU element-wise semantics (compute in the + // promoted type, store into the LHS type). See #988. The CPU path already does this + // in place, and real *= complex is rejected above. + if (rhs_is_host_scalar && Lt.device() != Device.cpu && Lt.dtype() > Rt.dtype()) { + Tensor promoted = Lt.astype(Type.type_promote(Lt.dtype(), Rt.dtype())); + iMul(promoted, R); + Lt = promoted.astype(Lt.dtype()); + return; + } + static const std::vector empty_mapper; // if contiguous, then no need to calculate the mappers if ((Lt.is_contiguous() && Rt.is_contiguous())) { @@ -36,7 +63,7 @@ namespace cytnx { detail::DispatchInplaceArithmeticCPU<1>(Lt, R, empty_mapper, empty_mapper, empty_mapper); } else { #ifdef UNI_GPU - checkCudaErrors(cudaSetDevice(Rt.device())); + checkCudaErrors(cudaSetDevice(Lt.device())); Tensor tmpo; if (Lt.dtype() <= Rt.dtype()) tmpo = Lt; @@ -61,11 +88,32 @@ namespace cytnx { Rt._impl->invmapper()); } else { #ifdef UNI_GPU - cytnx_error_msg(true, - "[iMul][on GPU/CUDA] error two tensors must be contiguous. Call " - "Contiguous_() or Contiguous() first%s", - "\n"); - + if (R._impl->storage()._impl->size() == 1) { + // Broadcast scalar RHS: the rconst kernel path scales every element of the + // (possibly non-contiguous) LHS storage in place, so the layout mappers are + // irrelevant and the result is correct. This is the #988 regression fix that lets + // e.g. `gpu_tensor.permute(...); gpu_tensor *= 2.0;` succeed again. + checkCudaErrors(cudaSetDevice(Lt.device())); + Tensor tmpo; + if (Lt.dtype() <= Rt.dtype()) + tmpo = Lt; + else + tmpo = Lt.clone(); + linalg_internal::lii.cuAri_ii[Lt.dtype()][Rt.dtype()]( + tmpo._impl->storage()._impl, Lt._impl->storage()._impl, R._impl->storage()._impl, + Lt._impl->storage()._impl->size(), Lt._impl->shape(), Lt._impl->invmapper(), + Rt._impl->invmapper(), 1); + if (Lt.dtype() > Rt.dtype()) Lt = tmpo; + } else { + // Genuine non-contiguous tensor*=tensor: unlike cuAdd/cuSub, the cuMul GPU kernels + // ignore the layout mappers, so routing here would silently pair mismatched + // elements. Fail loudly instead of corrupting data. Call Contiguous_()/Contiguous() + // on the operands first. + cytnx_error_msg(true, + "[iMul][on GPU/CUDA] non-contiguous tensor*=tensor is not supported. " + "Call Contiguous_() or Contiguous() on the operands first%s", + "\n"); + } #else cytnx_error_msg(true, "[Mul] fatal error, the tensor is on GPU without CUDA support.%s", "\n"); diff --git a/src/linalg/iSub.cpp b/src/linalg/iSub.cpp index f98a5ad45..ea31e196f 100644 --- a/src/linalg/iSub.cpp +++ b/src/linalg/iSub.cpp @@ -11,8 +11,21 @@ namespace cytnx { namespace linalg { void iSub(Tensor &Lt, const Tensor &Rt) { - cytnx_error_msg(Lt.device() != Rt.device(), + // A length-1 RHS that stays on the host is treated as a broadcast scalar: the GPU kernels + // read it with a host-side dereference and pass it into the kernel by value, so it needs + // neither a device match nor a per-call H2D copy of the scalar. See #988. + const bool rhs_is_host_scalar = + (Rt.device() == Device.cpu && Rt._impl->storage()._impl->size() == 1); + cytnx_error_msg(Lt.device() != Rt.device() && !rhs_is_host_scalar, "[iSub] error, the two tensors have to be on the same device.%s", "\n"); + // In-place ops write the result back into the LHS storage, so a complex result cannot be + // stored in a real LHS. Guard here (device-independent) so the GPU path throws like the CPU + // path instead of reinterpreting the real output buffer as complex and corrupting it. See + // #988. + cytnx_error_msg(!Type.is_complex(Lt.dtype()) && Type.is_complex(Rt.dtype()), + "[iSub] Cannot perform real -= complex in-place: a complex result cannot be " + "stored in a real tensor.%s", + "\n"); if (!(Rt.shape().size() == 1 && Rt.shape()[0] == 1)) { cytnx_error_msg(Lt.shape() != Rt.shape(), @@ -28,6 +41,20 @@ namespace cytnx { R = Rt; } + // GPU broadcast scalar with a LHS *narrower* than the promoted dtype (e.g. a Float tensor + // minus a Double scalar, or an integer tensor minus a fractional scalar): the in-place GPU + // kernels write the promoted-width result straight into the narrower LHS buffer, an + // out-of-bounds write that corrupts memory. Compute in the promoted dtype, then truncate + // back to the LHS dtype -- this matches the CPU element-wise semantics (compute in the + // promoted type, store into the LHS type). See #988. The CPU path already does this + // in place, and real -= complex is rejected above. + if (rhs_is_host_scalar && Lt.device() != Device.cpu && Lt.dtype() > Rt.dtype()) { + Tensor promoted = Lt.astype(Type.type_promote(Lt.dtype(), Rt.dtype())); + iSub(promoted, R); + Lt = promoted.astype(Lt.dtype()); + return; + } + static const std::vector empty_mapper; // if contiguous, then no need to calculate the mappers if ((Lt.is_contiguous() && Rt.is_contiguous())) { @@ -36,7 +63,7 @@ namespace cytnx { detail::DispatchInplaceArithmeticCPU<2>(Lt, R, empty_mapper, empty_mapper, empty_mapper); } else { #ifdef UNI_GPU - checkCudaErrors(cudaSetDevice(Rt.device())); + checkCudaErrors(cudaSetDevice(Lt.device())); Tensor tmpo; if (Lt.dtype() <= Rt.dtype()) tmpo = Lt; @@ -61,7 +88,7 @@ namespace cytnx { Rt._impl->invmapper()); } else { #ifdef UNI_GPU - checkCudaErrors(cudaSetDevice(Rt.device())); + checkCudaErrors(cudaSetDevice(Lt.device())); Tensor tmpo; if (Lt.dtype() <= Rt.dtype()) tmpo = Lt; diff --git a/tests/gpu/Tensor_test.cpp b/tests/gpu/Tensor_test.cpp index 394ce3b96..5e48a7ed4 100644 --- a/tests/gpu/Tensor_test.cpp +++ b/tests/gpu/Tensor_test.cpp @@ -263,3 +263,105 @@ TEST_F(TensorTest, gpu_set) { // EXPECT_FALSE(tar3456.approx_eq(tarcomplex3456)); // EXPECT_TRUE(tone3456.approx_eq(tone3456.astype(Type.ComplexFloat), 1e-5)); // } + +// --------------------------------------------------------------------------- +// GPU scalar in-place arithmetic (issue #988). These mirror the CPU-side +// Tensor.ScalarInplace* tests on CUDA tensors and pin the non-contiguous +// *= / /= behaviour on the GPU path. Test names avoid underscores (#857). +// --------------------------------------------------------------------------- + +// Regression (#988): scalar *= on a *non-contiguous* GPU tensor used to throw +// "[iMul][on GPU/CUDA] ... must be contiguous". It must now scale in place, +// keeping the shared storage, with the layout mappers irrelevant to a +// broadcast scalar. +TEST(Tensor, GpuScalarInplaceNoncontigMul) { + Tensor a = arange(6).reshape({2, 3}).to(Device.cuda); // Double, contiguous + Tensor v = a.permute({1, 0}); // distinct impl, shared storage, non-contiguous + ASSERT_FALSE(v.is_contiguous()); + ASSERT_TRUE(is(a.storage(), v.storage())); + v *= 2.0; // must not throw + EXPECT_TRUE(is(a.storage(), v.storage())); // still in place / shared + Tensor a_cpu = a.to(Device.cpu); + for (cytnx_uint64 i = 0; i < 6; i++) + EXPECT_DOUBLE_EQ(a_cpu.storage().at(i), 2.0 * i); +} + +// Regression (#988): scalar /= on a non-contiguous GPU tensor. +TEST(Tensor, GpuScalarInplaceNoncontigDiv) { + Tensor a = arange(6).reshape({2, 3}).to(Device.cuda); + Tensor v = a.permute({1, 0}); + ASSERT_FALSE(v.is_contiguous()); + ASSERT_TRUE(is(a.storage(), v.storage())); + v /= 2.0; // must not throw + EXPECT_TRUE(is(a.storage(), v.storage())); + Tensor a_cpu = a.to(Device.cpu); + for (cytnx_uint64 i = 0; i < 6; i++) + EXPECT_DOUBLE_EQ(a_cpu.storage().at(i), i / 2.0); +} + +// Contiguous scalar ops on a GPU tensor still produce the right values. This +// also exercises the #988 efficiency path: the scalar wrapper stays on the +// host and is read by the GPU kernel with a host-side dereference. +TEST(Tensor, GpuScalarInplaceContiguousValues) { + Tensor a = arange(6).to(Device.cuda); + a *= 3.0; + a += 1.0; + EXPECT_EQ(a.device(), Device.cuda); + Tensor a_cpu = a.to(Device.cpu); + for (cytnx_uint64 i = 0; i < 6; i++) + EXPECT_DOUBLE_EQ(a_cpu.storage().at(i), 3.0 * i + 1.0); +} + +// Scalar in-place ops mutate the LHS storage in place (never detach), so a +// second handle onto the same GPU storage observes the change. Mirrors +// Tensor.ScalarInplaceSubMulDivKeepStorageSharing on CUDA. +TEST(Tensor, GpuScalarInplaceKeepsStorageSharing) { + Tensor a = zeros({4}, Type.Double, Device.cuda); + Tensor b = Tensor::from_storage(a.storage()); // distinct impl, shared storage + ASSERT_TRUE(is(a.storage(), b.storage())); + a += 1.0; + a -= 0.5; + a *= 4.0; + a /= 2.0; // ((0 + 1 - 0.5) * 4) / 2 == 1 + EXPECT_TRUE(is(a.storage(), b.storage())); + Tensor b_cpu = b.to(Device.cpu); + EXPECT_DOUBLE_EQ(b_cpu.storage().at(0), 1.0); +} + +// A double scalar must not promote a Float GPU tensor to Double. Mirrors +// Tensor.ScalarInplaceOpsPreserveDtype on CUDA. +TEST(Tensor, GpuScalarInplacePreserveDtype) { + Tensor a = ones({2}, Type.Float, Device.cuda); + a += 1.0; // double scalar + a -= 0.5; + a *= 2.0; + a /= 3.0; // ((1 + 1 - 0.5) * 2) / 3 == 1 + EXPECT_EQ(a.dtype(), Type.Float); + EXPECT_EQ(a.device(), Device.cuda); + Tensor a_cpu = a.to(Device.cpu); + EXPECT_FLOAT_EQ(a_cpu.storage().at(0), 1.0f); +} + +// A real GPU tensor op= a complex scalar cannot store a complex result, so it +// must throw (as on CPU) rather than silently reinterpreting the real buffer +// as complex. Mirrors Tensor.ScalarInplaceRealPlusComplexThrows on CUDA. +TEST(Tensor, GpuScalarInplaceRealOpComplexThrows) { + Tensor a = zeros({2}, Type.Double, Device.cuda); + EXPECT_THROW(a += cytnx_complex128(0, 1), std::logic_error); + EXPECT_THROW(a -= cytnx_complex128(0, 1), std::logic_error); + EXPECT_THROW(a *= cytnx_complex128(0, 1), std::logic_error); + EXPECT_THROW(a /= cytnx_complex128(0, 1), std::logic_error); +} + +// The cuMul/cuDiv GPU kernels (unlike cuAdd/cuSub) do not consume layout +// mappers, so a genuine non-contiguous tensor*=tensor must fail loudly instead +// of silently pairing mismatched elements. (The scalar broadcast case above is +// still supported because it ignores the mappers.) +TEST(Tensor, GpuNoncontigTensorTensorMulDivThrows) { + Tensor a = arange(6).reshape({2, 3}).to(Device.cuda); + Tensor b = arange(6).reshape({3, 2}).to(Device.cuda).permute({1, 0}); // {2,3}, non-contiguous + ASSERT_EQ(a.shape(), b.shape()); + ASSERT_FALSE(b.is_contiguous()); + EXPECT_THROW(a *= b, std::logic_error); + EXPECT_THROW(a /= b, std::logic_error); +} From 3fda11fc431120cf774896a6fdb2b5e2c8dcc945 Mon Sep 17 00:00:00 2001 From: Pochung Chen Date: Wed, 8 Jul 2026 21:04:17 +0800 Subject: [PATCH 2/2] ci: retrigger checks after retarget to master Empty commit: the pull_request workflows last ran against the pre-merge #980 base branch (since deleted); this fires a fresh full run (BuildAndTest, DownstreamFindPackage, wheels) with base=master. Co-Authored-By: Claude Fable 5