diff --git a/.claude/skills/test-coverage-benchmark-workflow/SKILL.md b/.claude/skills/test-coverage-benchmark-workflow/SKILL.md
new file mode 100644
index 0000000..6183c59
--- /dev/null
+++ b/.claude/skills/test-coverage-benchmark-workflow/SKILL.md
@@ -0,0 +1,118 @@
+---
+name: test-coverage-benchmark-workflow
+description: Run AbstractOperators.jl's long-running test, coverage, and benchmark workflow — filtered TestItemRunner runs, coverage capture, and local/CI AirspeedVelocity benchmark comparisons. Use when running the full test suite, generating coverage, or comparing performance against master.
+---
+
+1. Start from the smallest relevant test scope; prefer a persistent Julia REPL for repeated
+   filtered `TestItemRunner.run_tests(...)` calls.
+2. Fix real implementation bugs in source instead of weakening tests; rerun the same filtered
+   slice until green, then expand to adjacent slices, then run the full suite.
+3. Capture all run logs under `.temp/`.
+4. For performance-sensitive changes, benchmark before and after; run focused ASV filters
+   first, then a single full ASV comparison for final validation. Treat
+   `speedup + uncertainty < 0.95` (master/dirty ratio) as a significant regression.
+5. Prefer representative large inputs for linear and nonlinear operators to reduce
+   microbenchmark noise, but wrap only fast operators in calculus operators to measure the
+   calculus overhead itself.
+6. Use AirspeedVelocity with an explicit script path when comparing against revisions that do
+   not yet contain the benchmark file.
+
+Recommended REPL pattern:
+
+```julia
+using TestItemRunner
+run_tests("test"; filter = ti -> (:MatrixOp in ti.tags) && (:linearoperator in ti.tags))
+```
+
+Main package coverage (also exercises subpackage and extension code — DSPOperators,
+FFTWOperators, NFFTOperators, and WaveletOperators have no standalone `test/` directory, and
+extensions are exercised through the parent package's tests):
+
+```sh
+julia --project=test --code-coverage=user test/runtests.jl
+```
+
+Process coverage after a local run:
+
+```sh
+julia -e 'using Coverage; Coverage.LCOV.writefile("lcov.info", Coverage.process_folder())'
+```
+
+Filtered test run:
+
+```julia
+using TestItemRunner
+TestItemRunner.run_tests(pwd(); filter = ti -> :MatrixOp in ti.tags) # by tag
+TestItemRunner.run_tests(pwd(); filter = ti -> ti.name == "DCT")     # by test name
+```
+
+### Local benchmark comparison with AirspeedVelocity
+
+```sh
+mkdir -p .temp/asv
+benchpkg \
+  --path . \
+  --rev master,dirty \
+  --script benchmark/benchmarks.jl \
+  --output-dir .temp/asv \
+  --exeflags="--threads=4"
+```
+
+Filtered comparison for a single benchmark family:
+
+```sh
+mkdir -p .temp/asv
+benchpkg \
+  --path . \
+  --rev master,dirty \
+  --script benchmark/benchmarks.jl \
+  --output-dir .temp/asv \
+  --exeflags="--threads=4" \
+  --add RecursiveArrayTools \
+  --filter MIMOFilt
+```
+
+Render a comparison table:
+
+```sh
+benchpkgtable \
+  --path . \
+  --rev master,dirty \
+  --input-dir .temp/asv \
+  --ratio \
+  --mode time,memory
+```
+
+### CI benchmark comparison (GitHub Actions)
+
+The GitHub Actions benchmark CI does **not** use the AirspeedVelocity action because the
+root-level Julia workspace (`[workspace]` in `Project.toml`) causes that action's
+revision-management to mis-resolve the monorepo subprojects. Instead, two workflows implement
+a fork-safe two-stage approach:
+
+- **`benchmark.yml`** — unprivileged `pull_request` job that checks out both the base and head
+  revisions, runs `benchmark/compare.jl` against explicit worktree paths, and uploads
+  `body.md`, `pr_number.txt`, and `julia_version.txt` as an artifact.
+- **`post_benchmark_comment.yml`** — privileged `workflow_run` job that downloads the artifact
+  and creates or updates the PR comment.
+
+The comparison table mirrors AirspeedVelocity output with separate Time and Memory sections,
+base/head columns, a ratio column, and emoji indicators:
+
+- 🚀 significant speedup: `ratio − ratio_err > 1.2` (time) or `ratio < 0.5` (memory)
+- 🐢 significant slowdown: `ratio + ratio_err < 0.8` (time) or `ratio > 1.5` (memory)
+
+To run the comparison locally with the same script used by CI:
+
+```sh
+git worktree add .temp/base master
+
+julia --project=benchmark benchmark/compare.jl \
+  --base-dir  .temp/base \
+  --head-dir  . \
+  --output-dir .temp/bench-compare \
+  --pr        0 \
+  --julia-version "$(julia -e 'print(VERSION)')"
+
+cat .temp/bench-compare/body.md
+```
diff --git a/.github/agents/julia.agent.md b/.github/agents/julia.agent.md
deleted file mode 100644
index b06708f..0000000
--- a/.github/agents/julia.agent.md
+++ /dev/null
@@ -1,116 +0,0 @@
----
-description: "Use when improving Julia code quality with very long test suites, slow CI, TestItemRunner tagging/filtering, iterative fix-and-rerun loops, or flaky tests. Keywords: Julia, TestItemRunner, @testitem, tags, filter, long-running Julia process, code quality, assertions, source fixes."
-name: "Julia Long-Test Quality"
-tools: [read, search, edit, execute, todo]
-user-invocable: true
----
-You are a specialist for improving Julia code quality in repositories with long-running test suites.
-
-## Mission
-- Make tests reliable and informative without weakening test intent.
-- Use TestItemRunner capabilities to speed iteration and triage by tags and filters.
-- Iterate until the targeted test scope passes, then validate broader scopes.
-
-## Hard Constraints
-- Never remove assertions to make tests pass.
-- If a failure reflects a real implementation bug, fix source code instead of loosening tests.
-- Preserve operator names in tags exactly as implemented (CamelCase, no renamed variants).
-- Keep changes minimal and localized; avoid unrelated refactors.
-
-## Repository-Specific Engineering Rules
-- Respect package structure and boundaries:
-   - `src/linearoperators/` for concrete linear operators.
-   - `src/nonlinearoperators/` for nonlinear operators.
-   - `src/calculus/` for operator calculus/composition.
-   - `src/batching/` for batch operators.
-- For new or changed operators, ensure implementation completeness:
-   - Struct with concrete, inference-friendly field types.
-   - Constructors for dimension tuple and/or data-driven construction.
-   - Forward path `mul!(y, op, x)` and adjoint path dispatch via `AdjointOperator`.
-   - Trait and property behavior remains consistent (`is_linear`, `is_diagonal`, rank/invertibility traits).
-   - Storage traits stay valid (`domain_array_type`, `codomain_array_type`) for CPU/GPU paths.
-- Prefer `copy_operator(op; array_type=nothing, threaded=nothing)` behavior when changing copy semantics:
-   - Deep-copy mutable working buffers only.
-   - Share immutable and read-only references.
-- Keep test files standalone-capable and aligned with TestItems setup modules.
-- Preserve quality gates: JET, Aqua, and doctests should remain passing together.
-- Use Runic formatting checks when editing Julia source or tests.
-
-## Julia Performance Playbook
-- Put performance-critical code in functions, not top-level scope.
-- Avoid untyped globals in hot paths; use function arguments and `const` globals where appropriate.
-- Prefer concrete field/container types; avoid abstract fields like `Function`, `AbstractArray`, or `Integer` in performance-sensitive structs.
-- Maintain type stability:
-   - Avoid variable type changes within loops.
-   - Use `zero(x)`, `oneunit(T)`, and stable return types.
-   - Use function barriers for setup-vs-kernel separation.
-- Measure, do not guess:
-   - Use `BenchmarkTools` for benchmarks.
-   - Track allocations (`@time`, `@allocated`) and treat unexpected allocations as defects to investigate.
-   - Use `@code_warntype` and JET to diagnose inference issues.
-- Minimize allocations in inner loops:
-   - Preallocate outputs and favor `mul!`/in-place APIs.
-   - Use broadcast fusion (`@.` / dotted ops) when beneficial.
-   - Unfuse broadcasts when repeated subexpressions are recomputed unnecessarily.
-   - Use `@views` for slicing when copy cost dominates.
-- Iterate arrays in memory-friendly order (column-major access patterns).
-- For threaded Julia code that also calls BLAS, avoid oversubscription (often `OPENBLAS_NUM_THREADS=1` is best with multithreaded Julia; validate on workload).
-- Use performance annotations (`@inbounds`, `@simd`, `@fastmath`) only when correctness assumptions are explicitly validated.
-
-## Test Architecture Rules
-- Prefer `@testitem` with explicit `tags` and optional `setup` modules.
-- Use tags that encode both operator and test type.
-- Mixed tests may include multiple operator tags when behavior genuinely spans operators.
-- Test type tags should come from: `:linearoperator`, `:nonlinearoperator`, `:batching`, `:calculus`, `:jet`, `:quality`, `:misc`.
-- Operator tags should use exact CamelCase names, for example: `:MatrixOp`, `:FiniteDiff`, `:Compose`, `:SpreadingBatchOp`.
-- Use `@run_package_tests filter=ti->...` to run focused slices.
-- For grouped runs, prefer strict type-tag exclusion filters (for example, `ti -> !(:jet in ti.tags)`).
-
-## JET.jl Requirements
-- Treat JET coverage as mandatory for all public API.
-- Ensure JET test coverage includes all three modes:
-   - `JET.test_package(...)` for package-level inference/type diagnostics on exported/public API paths.
-   - `@test_opt ...` for representative public operations and constructors.
-   - `@test_call ...` for key public call signatures and runtime-like call paths.
-- Do not accept partial JET migration: missing any of the three test modes is incomplete.
-- When adding or changing public API, update JET tests in the same change.
-
-## Fast Iteration Workflow
-1. Start one long-running Julia REPL in the package test environment.
-2. Load TestItemRunner once.
-3. Run filtered test slices repeatedly (by operator/type tags).
-4. Fix failures immediately; rerun the same filtered slice until green.
-5. Expand to adjacent slices, then run full suite.
-6. Capture outputs from each run into `.temp/` files for traceability.
-
-Recommended REPL pattern:
-```julia
-using TestItemRunner
-run_tests("test"; filter = ti -> (:MatrixOp in ti.tags) && (:linearoperator in ti.tags))
-```
-
-Recommended shell pattern for captured logs:
-```sh
-mkdir -p .temp
-julia --project=test test/runtests.jl > .temp/test_runtests.log 2>&1
-julia --project=test test/jet/test_package.jl > .temp/test_jet_package.log 2>&1
-```
-
-## Failure Triage
-1. Read the exact failing assertion and stacktrace first.
-2. Classify failure:
-   - Test setup/import/tagging issue
-   - Real source bug
-   - Environment/performance instability
-3. For real bugs, patch source and keep/assert expected behavior in tests.
-4. For flaky perf tests, stabilize methodology (workload, sampling, thresholds) without dropping coverage.
-5. Re-run the smallest relevant filtered subset before broad reruns.
-
-## Output Requirements
-- Report what was changed and why.
-- List files touched.
-- Provide exact filtered test commands used.
-- State pass/fail counts for the final run.
-- Call out remaining risks or follow-up items.
-- IMPORTANT! Store all temporary run outputs only under `.temp/` inside the repository (no temp scripts and logs elsewhere).
-- When performance work is included, report allocation deltas and the exact benchmark commands used.
diff --git a/.github/instructions/julia-operator-engineering.instructions.md b/.github/instructions/julia-operator-engineering.instructions.md
deleted file mode 100644
index 498297e..0000000
--- a/.github/instructions/julia-operator-engineering.instructions.md
+++ /dev/null
@@ -1,29 +0,0 @@
----
-description: "Use when editing Julia operator implementations in AbstractOperators.jl. Covers operator completeness, package boundaries, storage traits, and behavior-preserving refactors."
-name: "Julia Operator Engineering"
-applyTo: "src/**/*.jl"
----
-
-# Julia Operator Engineering
-
-- Respect package boundaries:
-  - `src/linearoperators/` for concrete linear operators.
-  - `src/nonlinearoperators/` for nonlinear operators.
-  - `src/calculus/` for composition and operator calculus.
-  - `src/batching/` for batch operators.
-- For new or changed operators, keep implementation complete:
-  - constructors,
-  - forward `mul!`,
-  - adjoint `mul!` where applicable,
-  - size/domain/codomain/storage traits,
-  - property traits such as linearity, diagonal structure, and rank-related predicates.
-- `check` utility function must be called in all effective `mul!` paths to ensure consistent argument validation and error messages.
-- Preserve `domain_array_type` and `codomain_array_type` semantics and dispatch compatibility.
-- Constructors should expose an `array_type` keyword where storage backend selection is meaningful.
-- `domain_array_type`/`codomain_array_type` must remain consistent with constructor-selected storage.
-- When storage checks become stricter, fix operator traits and tests instead of relaxing `check`.
-- Prefer behavior-preserving refactors: extract helpers, separate setup from kernels, reduce method size, but do not weaken checks.
-- If modifying copy semantics, preserve the package convention that immutable/read-only arrays are shared while mutable working buffers are copied deliberately.
-- Keep source formatted with Runic-compatible Julia style.
-- GPU extensions live under `ext/GpuExt/` (triggered by `GPUArrays`). Override `mul!` there for any operator whose base implementation uses scalar indexing loops (`@nloops`, `@nref`, `@inbounds y[i] = b[j]`); replace with broadcast-over-view (`y .= view(b, idx...)`).
-- When overriding a threaded operator (e.g. `Variation{..., true}`) for GPU, delegate to the non-threaded variant (`Variation{..., false}`) — the threading strategy is CPU-only.
diff --git a/.github/instructions/julia-performance.instructions.md b/.github/instructions/julia-performance.instructions.md
deleted file mode 100644
index 3f0beb4..0000000
--- a/.github/instructions/julia-performance.instructions.md
+++ /dev/null
@@ -1,28 +0,0 @@
----
-description: "Use when optimizing Julia code Covers type stability, allocations, preallocation, views, threading, BLAS oversubscription, and measurement discipline."
-name: "Julia Performance"
-applyTo: "src/**/*.jl,benchmark/**/*.jl"
----
-
-# Julia Performance
-
-## Core Principles
-
-- Keep hot paths inside functions, not top-level scope.
-- Avoid untyped globals in performance-sensitive paths; use arguments and `const` globals when needed.
-- Favor concrete field and container types; avoid abstractly typed hot fields.
-- Preserve type stability:
-  - avoid changing variable type in loops,
-  - prefer stable return types,
-  - use function barriers to separate setup from kernels.
-- Prefer in-place APIs and preallocation over repeated temporary allocations.
-- Use `@views` when slicing would otherwise allocate unnecessarily.
-- Respect Julia's column-major memory order when writing loops.
-- Use `@inbounds`, `@simd`, and `@fastmath` only when their correctness assumptions are justified.
-- For threaded Julia code that also uses BLAS, avoid oversubscription and benchmark with explicit thread settings.
-- Measure performance changes instead of guessing:
-  - benchmark representative workloads,
-  - inspect allocations,
-  - use JET and `@code_warntype` for inference issues.
-- For benchmark harnesses, derive element types robustly when operator type traits may return wrapped array types.
-- Keep benchmark setup deterministic (`Random.seed!(0)`) and validate key benchmark states with one smoke `mul!` path before full runs.
diff --git a/.github/instructions/julia-testing-and-jet.instructions.md b/.github/instructions/julia-testing-and-jet.instructions.md
deleted file mode 100644
index f90f31f..0000000
--- a/.github/instructions/julia-testing-and-jet.instructions.md
+++ /dev/null
@@ -1,31 +0,0 @@
----
-description: "Use when editing Julia tests, TestItemRunner suites, JET coverage, Aqua checks, or doctests in AbstractOperators.jl. Covers tags, logging, and quality gates."
-name: "Julia Testing And JET"
-applyTo: "test/**/*.jl,docs/**/*.md"
----
-
-# Julia Testing And JET
-
-- Prefer `@testitem` with explicit tags and optional setup modules.
-- Use type tags from: `:linearoperator`, `:nonlinearoperator`, `:batching`, `:calculus`, `:jet`, `:quality`, `:misc`.
-- Operator tags must use exact CamelCase type names, for example `:MatrixOp`, `:FiniteDiff`, `:Compose`, `:SpreadingBatchOp`.
-- Mixed tests may use multiple operator tags when the behavior genuinely spans operators.
-- Use strict TestItemRunner filters when slicing the suite.
-- Treat JET as mandatory for all public API:
-  - `JET.test_package(...)`
-  - `@test_opt`
-  - `@test_call`
-- Public API changes must update JET tests in the same change.
-- Keep Aqua and doctests passing alongside functional tests.
-- Never remove assertions to force green tests.
-- All temporary test and benchmark outputs must go under `.temp/` only.
-- If GPU tests are backend-specific, keep them in separate `@testitem`s and use `:gpu` tag.
-- When `VERB` is enabled, print each running testitem name at test-runner filter time.
-- For local coverage, mirror CI with `julia --project=test --code-coverage=user test/runtests.jl`, then process `*.cov` / `*.info` artifacts into `lcov.info` if needed.
-- Subpackages (DSPOperators, FFTWOperators, NFFTOperators, WaveletOperators) have no standalone `test/` directory; they are tested and their coverage is gathered exclusively through the parent package's `test/` project. Do not attempt a separate subpackage coverage run.
-- Extension coverage should be gathered through the parent-package tests that load the relevant trigger packages; do not assume a separate extension-only coverage run exists.
-- JET `@test_opt` flags `array_type::Type` (unparameterized keyword) as a source of runtime dispatch. Use `array_type::Type{<:AbstractArray}` and avoid kwarg-to-kwarg forwarding; use a typed positional-arg helper (e.g., `_make_eye(T, dims, S)`) so JET can resolve dispatch statically.
-- When Aqua reports "Unexpected Pass" on a `@test_broken`/`broken=true` check, the underlying issue is now fixed — remove the workaround and use `Aqua.test_all(pkg)` unconditionally.
-- Agent sub-tasks frequently generate `Eye(T, dims, array_type)` (3 positional args) instead of `Eye(T, dims; array_type=...)` (keyword). Always verify agent output for this pattern.
-- Stochastic test assertions `op * randn(n) ≈ other_op * (op * randn(n))` are wrong when the two `randn` calls produce different vectors; always capture into a variable first.
-- When testing GPU storage-type propagation, add `@test domain_array_type(op) <: CUDA.CuArray` / `<: AMDGPU.ROCArray` assertions directly in the per-operator CUDA/AMDGPU `@testitem`.
diff --git a/.github/skills/julia-gpu-implementation/SKILL.md b/.github/skills/julia-gpu-implementation/SKILL.md
deleted file mode 100644
index cbaff6b..0000000
--- a/.github/skills/julia-gpu-implementation/SKILL.md
+++ /dev/null
@@ -1,50 +0,0 @@
----
-name: julia-gpu-implementation
-description: 'Use for GPU operator implementations, GPU extension fixes, backend-specific testitems, and GPU benchmark validation in AbstractOperators.jl.'
-argument-hint: 'Describe the operator, GPU backend, or benchmark you want to implement or validate'
-user-invocable: true
----
-
-# Julia GPU Implementation
-
-## When To Use
-
-- Implementing or fixing GPU overrides under `ext/GpuExt/`.
-- Adding or updating CUDA/AMDGPU testitems.
-- Debugging backend-specific dispatch, storage traits, or array conversion issues.
-- Extending benchmark coverage for GPU behavior.
-- Checking whether a CPU operator should get a GPU path or stay CPU-only.
-
-## Implementation Rules
-
-- Julia package extensions can only `import` the parent package, trigger package(s), and stdlib; if extension code needs a parent dependency API, expose it from the parent module first.
-- For FFT plans, prefer `inv(plan)` (AbstractFFTs-generic) over backend-specific `FFTW.plan_inv(...)` to keep CUDA/AMDGPU compatibility.
-- With JLArrays/GPUArrays, avoid `copyto!(gpu, cpu_view)` where the source is a `SubArray`; materialize first, for example with `src[1:n]`, or copy from a plain array.
-- Preserve backend storage semantics and trait dispatch when adding GPU methods.
-- Keep CPU-only implementation details out of GPU overrides unless the backend truly supports them.
-- For GPU `GetIndex` overrides, keep boolean-mask and integer-vector fancy indexing in CPU paths unless the backend support is verified.
-- When overriding a threaded operator for GPU, delegate to the non-threaded variant; threading strategy is CPU-only.
-- Prefer direct `CuArray(arr)` / `CUDA.zeros(...)` / `AMDGPU.ROCArray(arr)` / `AMDGPU.zeros(...)` calls over intermediate conversion variables.
-- Benchmark setup code should normalize wrapped domain and codomain type traits to scalar element types before calling `randn` or `zeros`.
-
-## Testing Rules
-
-- For honest GPU coverage, keep JLArray checks separate from real device checks and add backend-specific tags such as `:cuda` and `:amdgpu` plus runtime skip guards.
-- In `test/runtests.jl`, filter backend-tagged testitems when the runtime is unavailable, but keep per-test safety checks too.
-- Add explicit tests for `domain_array_type` and `codomain_array_type`, and verify that `op * x` allocates on the active backend.
-- When adding CUDA/AMDGPU companion tests, prefer direct backend array construction instead of temporary conversion variables.
-- For GPU `GetIndex` tests, restrict indices to ranges, colons, and scalar integers; bool-mask and integer-vector `view` forms are not universally supported across GPU backends.
-- Migrate GPU-backend storage-type assertions from central quality files into each operator's own CUDA/AMDGPU `@testitem` so they run with the functional tests.
-- Use direct `import CUDA` / `import AMDGPU` plus `functional()` guards in testitems; avoid try/catch gating.
-
-## Benchmarking Rules
-
-- Benchmark scripts under `benchmark/` must prefer local workspace package paths over registry-installed copies, otherwise GPU fixes in sibling packages can be silently skipped.
-- Use representative large inputs for GPU crossover studies and keep the measurement setup deterministic.
-- Capture benchmark logs and generated reports under `.temp/`.
-
-## Tooling Reminders
-
-- Agent sub-tasks frequently generate `Eye(T, dims, array_type)` with three positional arguments instead of `Eye(T, dims; array_type=...)` with a keyword; verify this pattern.
-- JET `@test_opt` catches runtime dispatch from `array_type::Type` when it is unparameterized; use `array_type::Type{<:AbstractArray}` and avoid kwarg-to-kwarg forwarding by routing through an internal helper.
-- When fixing an "unexpected pass" Aqua error, remove the workaround and use `Aqua.test_all(pkg)` once the underlying issue is fixed.
diff --git a/.github/skills/julia-long-test-workflow/SKILL.md b/.github/skills/julia-long-test-workflow/SKILL.md
deleted file mode 100644
index 5a185f3..0000000
--- a/.github/skills/julia-long-test-workflow/SKILL.md
+++ /dev/null
@@ -1,141 +0,0 @@
----
-name: julia-long-test-workflow
-description: 'Use for long-running Julia test suites, TestItemRunner filtering, JET triage, benchmark-driven refactoring, and AirspeedVelocity branch-vs-master comparisons in AbstractOperators.jl.'
-argument-hint: 'Describe the operator, test group, or benchmark comparison you want to run'
-user-invocable: true
----
-
-# Julia Long Test Workflow
-
-## When To Use
-
-- Iterating on a failing Julia test suite that is too slow to rerun wholesale.
-- Narrowing failures with TestItemRunner tags or filenames.
-- Verifying JET coverage for public API.
-- Refactoring performance-sensitive operator code and checking for regressions.
-- Comparing the current branch against `master` using AirspeedVelocity.
-
-## Workflow
-
-1. Start from the smallest relevant test scope.
-2. Prefer a persistent Julia REPL for repeated filtered `TestItemRunner.run_tests(...)` calls.
-3. Fix real implementation bugs in source instead of weakening tests.
-4. Capture all run logs under `.temp/`.
-5. For performance-sensitive changes, benchmark before and after.
-6. Run focused ASV filters first, then a single full ASV comparison for final validation.
-7. Treat `speedup + uncertainty < 0.95` (master/dirty ratio) as a significant regression.
-8. Prefer representative large inputs for linear and nonlinear operators to reduce microbenchmark noise, but wrap only fast operators in calculus operators to measure the calculus overhead itself.
-9. Use AirspeedVelocity with an explicit script path when comparing against revisions that do not yet contain the benchmark file.
-
-## Common Commands
-
-Main package coverage:
-
-```sh
-julia --project=test --code-coverage=user test/runtests.jl
-```
-
-Subpackage coverage (DSPOperators, FFTWOperators, NFFTOperators, WaveletOperators have **no** standalone `test/` directory):
-
-> All subpackage code and their GPU extensions are exercised by the parent package's
-> `test/` project. Run the same coverage command above; the `.cov` files under each
-> subpackage's `src/` will be populated automatically.
-
-Process coverage after a local run:
-
-```sh
-julia -e 'using Coverage; Coverage.LCOV.writefile("lcov.info", Coverage.process_folder())'
-```
-
-Filtered test run:
-
-```julia
-using TestItemRunner
-TestItemRunner.run_tests(pwd(); filter = ti -> :MatrixOp in ti.tags) # example of filtering by tag
-TestItemRunner.run_tests(pwd(); filter = ti -> ti.name == "DCT") # example of filtering by test name instead of tags
-```
-
-### Local benchmark comparison with AirspeedVelocity
-
-AirspeedVelocity works well for local branch-vs-master comparisons and is the
-recommended tool for interactive performance investigation:
-
-```sh
-mkdir -p .temp/asv
-benchpkg \
-  --path . \
-  --rev master,dirty \
-  --script benchmark/benchmarks.jl \
-  --output-dir .temp/asv \
-  --exeflags="--threads=4"
-```
-
-Filtered AirSpeedVelocity comparison for a single benchmark family:
-
-```sh
-mkdir -p .temp/asv
-benchpkg \
-  --path . \
-  --rev master,dirty \
-  --script benchmark/benchmarks.jl \
-  --output-dir .temp/asv \
-  --exeflags="--threads=4" \
-  --add RecursiveArrayTools \
-  --filter MIMOFilt
-```
-
-Render a comparison table:
-
-```sh
-benchpkgtable \
-  --path . \
-  --rev master,dirty \
-  --input-dir .temp/asv \
-  --ratio \
-  --mode time,memory
-```
-
-> **Note:** Use AirspeedVelocity with an explicit `--script` path when comparing
-> against revisions that do not yet contain the benchmark file.
-
-### CI benchmark comparison (GitHub Actions)
-
-The GitHub Actions benchmark CI does **not** use the AirspeedVelocity action
-because the root-level Julia workspace (`[workspace]` in `Project.toml`) causes
-that action's revision-management to mis-resolve the monorepo subprojects.
-Instead, two workflows implement a fork-safe two-stage approach:
-
-- **`benchmark.yml`** – unprivileged `pull_request` job that checks out both
-  the base and head revisions, runs `benchmark/compare.jl` against explicit
-  worktree paths, and uploads `body.md`, `pr_number.txt`, and
-  `julia_version.txt` as an artifact.
-- **`post_benchmark_comment.yml`** – privileged `workflow_run` job that
-  downloads the artifact and creates or updates the PR comment.
-
-The comparison table mirrors AirspeedVelocity output with separate Time and
-Memory sections, base/head columns, a ratio column, and emoji indicators:
-- 🚀 significant speedup: `ratio − ratio_err > 1.2` (time) or `ratio < 0.5` (memory)
-- 🐢 significant slowdown: `ratio + ratio_err < 0.8` (time) or `ratio > 1.5` (memory)
-
-To run the comparison locally with the same script used by CI:
-
-```sh
-# Check out base separately, e.g. in a worktree:
-git worktree add .temp/base master
-
-julia --project=benchmark benchmark/compare.jl \
-  --base-dir  .temp/base \
-  --head-dir  . \
-  --output-dir .temp/bench-compare \
-  --pr        0 \
-  --julia-version "$(julia -e 'print(VERSION)')"
-
-cat .temp/bench-compare/body.md
-```
-
-## Done Criteria
-
-- Targeted tests pass.
-- JET coverage remains complete for public API touched.
-- Benchmark deltas are measured and reported.
-- Logs are saved under `.temp/`.
diff --git a/.gitignore b/.gitignore
index bcc7eb8..6b3911c 100644
--- a/.gitignore
+++ b/.gitignore
@@ -18,5 +18,8 @@ docs/Manifest.toml
 Manifest.toml
 Manifest-v*.toml
 .temp/
+
+# Julia coverage files
+*.cov
 test/gpu_env/
 benchmark/gpu_env/
diff --git a/AGENTS.md b/AGENTS.md
index 8b50d61..a28b3a8 100644
--- a/AGENTS.md
+++ b/AGENTS.md
@@ -1,25 +1,115 @@
 # AGENTS.md
 
-This repository uses layered guidance. Follow it in this order:
+Guidance for agents working in AbstractOperators.jl.
 
-1. Read this file first.
-2. Read any applicable files under `.github/instructions/` whose `applyTo` pattern matches the files you will edit.
-3. Read the matching skill under `.github/skills/` when the task clearly matches a skill's scope.
-4. Then inspect the target source files before editing.
+## Mission
 
-## How to choose guidance
+Make changes to operators and their tests reliable and informative without weakening test
+intent. Never remove assertions to force green tests. If a failure reflects a real
+implementation bug, fix the source instead of loosening the test. Keep changes minimal and
+localized; avoid unrelated refactors.
 
-- Use `julia-operator-engineering.instructions.md` for changes under `src/**/*.jl`.
-- Use `julia-performance.instructions.md` for code under `src/**/*.jl` and `benchmark/**/*.jl` when performance is relevant.
-- Use `julia-testing-and-jet.instructions.md` for `test/**/*.jl` and docs-backed test guidance.
-- Use `.github/skills/julia-long-test-workflow/SKILL.md` for long Julia test runs, filtered `TestItemRunner` work, JET triage, and AirspeedVelocity comparisons.
-- Use `.github/skills/julia-gpu-implementation/SKILL.md` for GPU operator implementations, GPU extensions, GPU-specific tests, and GPU benchmark validation.
+## Repository Layout & Operator Engineering
 
-## Working rules
+- For new or changed operators, keep the implementation complete:
+  - constructors,
+  - forward `mul!`,
+  - adjoint `mul!` where applicable,
+  - size/domain/codomain/storage traits,
+  - property traits such as linearity, diagonal structure, and rank-related predicates.
+- `check` utility function must be called in all effective `mul!` paths to ensure consistent
+  argument validation and error messages.
+- Preserve `domain_array_type`/`codomain_array_type` semantics and dispatch compatibility;
+  keep them consistent with constructor-selected storage.
+- Constructors should expose an `array_type` keyword where storage backend selection is
+  meaningful.
+- When storage checks become stricter, fix operator traits and tests instead of relaxing
+  `check`.
+- Prefer behavior-preserving refactors: extract helpers, separate setup from kernels, reduce
+  method size, but do not weaken checks.
+- If modifying copy semantics, preserve the convention that immutable/read-only arrays are
+  shared while mutable working buffers are copied deliberately (see
+  `copy_operator(op; array_type=nothing, threaded=nothing)`).
+- Keep source formatted with Runic-compatible Julia style.
 
-- Prefer the smallest skill and instruction set that fully covers the task.
-- Do not ignore a matching instruction file because a skill also exists; use both when they apply.
-- If multiple instruction files match, combine them rather than choosing only one.
-- If a task touches both implementation and tests, read both the source and test instruction files before editing.
-- Keep temporary artifacts under `.temp/`.
-- When in doubt, inspect the relevant files before making changes.
+GPU extension conventions live in `ext/GpuExt/CLAUDE.md`.
+
+## Performance
+
+- Measure, don't guess: use `BenchmarkTools`, track allocations (`@time`, `@allocated`) and
+  treat unexpected allocations as defects, use `@code_warntype` and JET to diagnose inference
+  issues.
+- Minimize allocations in inner loops: preallocate outputs, favor `mul!`/in-place APIs, use
+  broadcast fusion (`@.`) when beneficial, unfuse broadcasts when repeated subexpressions are
+  recomputed unnecessarily, use `@views` for slicing when copy cost dominates.
+- For threaded Julia code that also calls BLAS, avoid oversubscription (often
+  `OPENBLAS_NUM_THREADS=1` is best with multithreaded Julia; validate on workload).
+- Use `@inbounds`/`@simd`/`@fastmath` only when correctness assumptions are explicitly
+  validated.
+- Benchmark setup code should normalize wrapped domain and codomain type traits to scalar
+  element types before calling `randn`/`zeros`, use representative large inputs for GPU
+  crossover studies, and keep the measurement setup deterministic (`Random.seed!(0)`).
+
+## Testing & JET
+
+- Prefer `@testitem` with explicit tags and optional setup modules; keep test files
+  standalone-capable and aligned with TestItems setup modules.
+- Use type tags from: `:linearoperator`, `:nonlinearoperator`, `:batching`, `:calculus`,
+  `:jet`, `:quality`, `:misc`.
+- Operator tags must use exact CamelCase type names, e.g. `:MatrixOp`, `:FiniteDiff`,
+  `:Compose`, `:SpreadingBatchOp`. Mixed tests may use multiple operator tags when the
+  behavior genuinely spans operators.
+- Use `@run_package_tests filter=ti->...` / `TestItemRunner.run_tests(...)` for focused
+  slices; use strict tag-exclusion filters for grouped runs (e.g.
+  `ti -> !(:jet in ti.tags)`).
+- Treat JET as mandatory for all public API, across all three modes in the same change:
+  - `JET.test_package(...)` for package-level inference/type diagnostics,
+  - `@test_opt` for representative public operations and constructors,
+  - `@test_call` for key public call signatures and runtime-like call paths.
+  Missing any of the three is an incomplete migration. Public API changes must update JET
+  tests in the same change.
+- JET `@test_opt` flags `array_type::Type` (unparameterized keyword) as a source of runtime
+  dispatch. Use `array_type::Type{<:AbstractArray}` and avoid kwarg-to-kwarg forwarding; route
+  through a typed positional-arg helper (e.g. `_make_eye(T, dims, S)`) so JET can resolve
+  dispatch statically.
+- Keep Aqua and doctests passing alongside functional tests. When Aqua reports "Unexpected
+  Pass" on a `@test_broken`/`broken=true` check, the underlying issue is fixed — remove the
+  workaround and use `Aqua.test_all(pkg)` unconditionally.
+- If GPU tests are backend-specific, keep them in separate `@testitem`s with `:cuda`/`:amdgpu`
+  tags. In non-FFTW/non-DSP operator tests, prefer JLArray backend checks over CUDA/AMDGPU
+  device checks. Use direct `import CUDA`/`import AMDGPU` + `functional()` guards in
+  testitems; avoid try/catch gating. Restrict GPU `GetIndex` test indices to ranges, colons,
+  and scalar integers — bool-mask and integer-vector `view` forms are not universally
+  supported across GPU backends. Add `domain_array_type`/`codomain_array_type` tests and
+  verify `op * x` allocates on the active backend. Migrate GPU-backend storage-type assertions
+  into each operator's own CUDA/AMDGPU `@testitem` (e.g.
+  `@test domain_array_type(op) <: CUDA.CuArray`) so they run with the functional tests.
+- Stochastic test assertions like `op * randn(n) ≈ other_op * (op * randn(n))` are wrong when
+  the two `randn` calls produce different vectors — always capture into a variable first.
+- Agent sub-tasks frequently generate `Eye(T, dims, array_type)` (3 positional args) instead
+  of `Eye(T, dims; array_type=...)` (keyword). Always verify agent output for this pattern.
+- All temporary test and benchmark outputs must go under `.temp/` only.
+- When `VERB` is enabled, print each running testitem name at test-runner filter time.
+
+The long-running test/coverage/benchmark workflow (filtered TestItemRunner runs, coverage
+capture, local and CI AirspeedVelocity comparisons) is documented in the
+`test-coverage-benchmark-workflow` skill.
+
+## Failure Triage
+
+1. Read the exact failing assertion and stacktrace first.
+2. Classify the failure: test setup/import/tagging issue, real source bug, or
+   environment/performance instability.
+3. For real bugs, patch source and keep/assert expected behavior in tests.
+4. For flaky perf tests, stabilize methodology (workload, sampling, thresholds) without
+   dropping coverage.
+5. Re-run the smallest relevant filtered subset before broad reruns.
+
+## Output Requirements
+
+- Report what was changed and why; list files touched.
+- Provide the exact filtered test commands used and state pass/fail counts for the final run.
+- Call out remaining risks or follow-up items.
+- Store all temporary run outputs only under `.temp/` inside the repository.
+- When performance work is included, report allocation deltas and the exact benchmark commands
+  used.
diff --git a/CLAUDE.md b/CLAUDE.md
new file mode 120000
index 0000000..47dc3e3
--- /dev/null
+++ b/CLAUDE.md
@@ -0,0 +1 @@
+AGENTS.md
\ No newline at end of file
diff --git a/DSPOperators/ext/GpuExt/GpuExt.jl b/DSPOperators/ext/GpuExt/GpuExt.jl
index f6daeee..6dc5207 100644
--- a/DSPOperators/ext/GpuExt/GpuExt.jl
+++ b/DSPOperators/ext/GpuExt/GpuExt.jl
@@ -40,7 +40,7 @@ end
 
 function _make_gpu_filt(ref::AbstractGPUArray{T}, cpu_op::AbstractFilt{T, N}) where {T <: Real, N}
     fftlen = nextpow(2, cpu_op.dim_in[1] + length(cpu_op.b) - 1)
-    array_type = _wrapper_type(ref){T}
+    storage_type = _wrapper_type(ref){T}
     buf = similar(ref, T, fftlen)
     buf_fft = _gpu_complex_buffer(ref, fftlen)
     buf_out = similar(ref, T, fftlen)
@@ -57,7 +57,7 @@ function _make_gpu_filt(ref::AbstractGPUArray{T}, cpu_op::AbstractFilt{T, N}) wh
     h_fft_conj = similar(h_fft)
     h_fft_conj .= conj.(h_fft)
 
-    return GpuFilt{T, N, array_type, typeof(buf), typeof(buf_fft), typeof(plan_fwd), typeof(plan_inv)}(
+    return GpuFilt{T, N, storage_type, typeof(buf), typeof(buf_fft), typeof(plan_fwd), typeof(plan_inv)}(
         cpu_op.dim_in,
         collect(cpu_op.b),
         collect(cpu_op.a),
diff --git a/DSPOperators/src/Xcorr.jl b/DSPOperators/src/Xcorr.jl
index c1e1714..a825d4b 100644
--- a/DSPOperators/src/Xcorr.jl
+++ b/DSPOperators/src/Xcorr.jl
@@ -1,22 +1,5 @@
 export Xcorr
 
-# Adjoint FFT state for non-CPU array backends.
-# CPU uses a tiled FIR loop instead; adj_fft is Nothing for CPU arrays.
-# Hfft and Hc are separate type params because plan_rfft on a CPU-backed GPU
-# mock array (e.g. JLArrays) may return a plain Vector from `R * buf`, while
-# similar(h, Complex{T}, ...) returns the array's own complex type.
-struct XcorrAdjFFT{
-        Hfft <: AbstractVector, H <: AbstractVector, Hc <: AbstractVector,
-        P3 <: AbstractFFTs.Plan, P4 <: AbstractFFTs.Plan,
-    }
-    fftlen::Int
-    h_fft::Hfft  # rfft(h padded); type may differ from buf_c
-    buf::H       # scratch buffer size fftlen
-    buf_c::Hc    # complex scratch buffer
-    R::P3        # rfft/fft plan
-    I::P4        # irfft/ifft plan
-end
-
 """
 	Xcorr([domain_type=Float64::Type,] dim_in::Tuple, h::AbstractVector)
 	Xcorr(x::AbstractVector, h::AbstractVector)
@@ -33,7 +16,8 @@ julia> Xcorr(Float64, (10,), [1.0, 0.5, 0.2])
 """
 struct Xcorr{
         T, H <: AbstractVector{T}, Hc <: AbstractVector,
-        P1 <: AbstractFFTs.Plan, P2 <: AbstractFFTs.Plan, Adj,
+        P1 <: AbstractFFTs.Plan, P2 <: AbstractFFTs.Plan,
+        P3 <: AbstractFFTs.Plan, P4 <: AbstractFFTs.Plan,
     } <: LinearOperator
     dim_in::Tuple{Int}
     h::H
@@ -45,8 +29,13 @@ struct Xcorr{
     buf_fwd_c::Hc    # complex scratch buffer
     R_fwd::P1        # rfft plan, fftlen_fwd
     I_fwd::P2        # irfft plan, fftlen_fwd
-    # Adjoint: XcorrAdjFFT{...} for GPU backends, Nothing for CPU
-    adj_fft::Adj
+    # Adjoint pass (conv(b, h) and slice)
+    fftlen_adj::Int
+    h_fft_adj::Hc    # rfft(h padded to fftlen_adj)
+    buf_adj::H       # scratch buffer size fftlen_adj
+    buf_adj_c::Hc    # complex scratch buffer
+    R_adj::P3        # rfft plan, fftlen_adj
+    I_adj::P4        # irfft plan, fftlen_adj
 end
 
 # FFT planning flags: FFTW.MEASURE only for CPU Arrays; no flags for GPU backends.
@@ -62,6 +51,7 @@ function Xcorr(domain_type::Type, DomainDim::NTuple{N, Int}, h::H) where {H <: A
     m = length(h)
     padlen = max(n, m)
     outlen = 2 * padlen - 1
+
     plan_kw = _xcorr_plan_kwargs(H)
 
     # Forward pass plans
@@ -82,38 +72,31 @@ function Xcorr(domain_type::Type, DomainDim::NTuple{N, Int}, h::H) where {H <: A
     h_fft_conj = conj.(R_fwd * buf_fwd)
     fill!(buf_fwd, zero(domain_type))
 
-    # Adjoint: CPU uses tiled FIR — no FFT state needed.
-    # GPU backends allocate FFT plans; same fftlen as forward pass is correct
-    # (wrap-around from h only affects positions < m ≤ padlen, outside the
-    # extracted range padlen..padlen+n-1).
-    if H <: Array
-        adj_fft = nothing
+    # Adjoint pass: CPU uses tiled FIR — no FFT state needed.
+    # GPU backends allocate FFT plans; same fftlen as forward pass is correct.
+    fftlen_adj = fftlen_fwd
+    buf_adj = similar(h, fftlen_adj)
+    if domain_type <: Real
+        R_adj = plan_rfft(buf_adj; plan_kw...)
+        buf_adj_c = similar(h, Complex{domain_type}, fftlen_adj ÷ 2 + 1)
+        I_adj = plan_irfft(buf_adj_c, fftlen_adj; plan_kw...)
     else
-        fftlen_adj = fftlen_fwd
-        buf_adj = similar(h, fftlen_adj)
-        if domain_type <: Real
-            R_adj = plan_rfft(buf_adj; plan_kw...)
-            buf_adj_c = similar(h, Complex{domain_type}, fftlen_adj ÷ 2 + 1)
-            I_adj = plan_irfft(buf_adj_c, fftlen_adj; plan_kw...)
-        else
-            R_adj = plan_fft(buf_adj; plan_kw...)
-            buf_adj_c = similar(buf_adj)
-            I_adj = inv(R_adj)
-        end
-        fill!(buf_adj, zero(domain_type))
-        copyto!(view(buf_adj, 1:m), h)
-        h_fft_adj = R_adj * buf_adj
-        fill!(buf_adj, zero(domain_type))
-        adj_fft = XcorrAdjFFT(fftlen_adj, h_fft_adj, buf_adj, buf_adj_c, R_adj, I_adj)
+        R_adj = plan_fft(buf_adj; plan_kw...)
+        buf_adj_c = similar(buf_adj)
+        I_adj = inv(R_adj)
     end
+    fill!(buf_adj, zero(domain_type))
+    copyto!(view(buf_adj, 1:m), h)
+    h_fft_adj = R_adj * buf_adj
+    fill!(buf_adj, zero(domain_type))
 
     return Xcorr{
         domain_type, typeof(h), typeof(buf_fwd_c),
-        typeof(R_fwd), typeof(I_fwd), typeof(adj_fft),
+        typeof(R_fwd), typeof(I_fwd), typeof(R_adj), typeof(I_adj),
     }(
         DomainDim, h,
         fftlen_fwd, padlen, h_fft_conj, buf_fwd, buf_fwd_c, R_fwd, I_fwd,
-        adj_fft,
+        fftlen_adj, h_fft_adj, buf_adj, buf_adj_c, R_adj, I_adj,
     )
 end
 
@@ -149,19 +132,19 @@ function mul!(y, L::AdjointOperator{<:Xcorr{T, <:Array{T}}}, b) where {T}
     return y
 end
 
-# GPU adjoint: FFT-based conv via XcorrAdjFFT
-function mul!(y, L::AdjointOperator{<:Xcorr{T, <:Any, <:Any, <:Any, <:Any, <:XcorrAdjFFT}}, b) where {T}
+# GPU adjoint: FFT-based conv
+function mul!(y, L::AdjointOperator{<:Xcorr{T}}, b) where {T}
     check(y, L, b)
     A = L.A
-    adj = A.adj_fft
     n = length(y)
     outlen = length(b)
-    fill!(adj.buf, zero(T))
-    copyto!(view(adj.buf, 1:outlen), b)
-    mul!(adj.buf_c, adj.R, adj.buf)
-    adj.buf_c .*= adj.h_fft
-    mul!(adj.buf, adj.I, adj.buf_c)
-    y .= @view(adj.buf[A.padlen:(A.padlen + n - 1)])
+    fill!(A.buf_adj, zero(T))
+    copyto!(view(A.buf_adj, 1:outlen), b)
+    mul!(A.buf_adj_c, A.R_adj, A.buf_adj)
+    A.buf_adj_c .*= A.h_fft_adj
+    mul!(A.buf_adj, A.I_adj, A.buf_adj_c)
+    padlen = A.padlen
+    y .= @view(A.buf_adj[padlen:(padlen + n - 1)])
     return y
 end
 
diff --git a/FFTWOperators/src/Shift.jl b/FFTWOperators/src/Shift.jl
index 4463de1..adb523d 100644
--- a/FFTWOperators/src/Shift.jl
+++ b/FFTWOperators/src/Shift.jl
@@ -439,13 +439,14 @@ function _is_dft_op(op, side)
         return true
     elseif op isa Compose
         subops = AbstractOperators.get_operators(op)
-        if all(o -> is_diagonal(o) || _is_dft_op(o, side), subops)
-            # it is an elementwise modification of a DFT/IDFT
-            return true
+        if side == :domain
+            # Domain shift: innermost (first) op must be DFT-like; all outer ops must be
+            # diagonal so they commute with SignAlternation.
+            return _is_dft_op(first(subops), side) && all(is_diagonal, subops[2:end])
         else
-            # alternatively, it is enough to check the first/last operator
-            op = size == :domain ? first(subops) : last(subops)
-            return _is_dft_op(op, side)
+            # Codomain shift: outermost (last) op must be DFT-like; all inner ops must be
+            # diagonal so they commute with SignAlternation.
+            return _is_dft_op(last(subops), side) && all(is_diagonal, subops[1:(end - 1)])
         end
     else
         return false
diff --git a/docs/src/gpu.md b/docs/src/gpu.md
index b86f5e0..cfb2a12 100644
--- a/docs/src/gpu.md
+++ b/docs/src/gpu.md
@@ -144,10 +144,6 @@ dct_op = CpuOperatorWrapper(DCT(Float32, (64,)); array_type = CuArray{Float32})
 y_gpu = dct_op * x_gpu  # GPU in → CPU DCT → GPU out
 ```
 
-### WaveletOperators CPU-only status
-
-WaveletOperators.jl currently relies on CPU execution. Its operators do not yet support GPU arrays, so wavelet transforms should remain on CPU or be wrapped explicitly as CPU operators when building mixed CPU/GPU pipelines.
-
 ## CpuOperatorWrapper
 
 For operators that do not natively support GPU arrays (e.g., FFTWOperators DCT, custom CPU-only operators), use `CpuOperatorWrapper`. This wrapper preallocates CPU buffers for the operator's domain and codomain, allowing GPU arrays to be passed in and out while the computation happens on CPU.
@@ -166,7 +162,7 @@ using AbstractOperators, CUDA
 op = FiniteDiff(Float32, (64,))  # or any FFTWOperators, etc.
 
 # Wrap it — preallocates CPU buffers for domain and codomain
-wrapper = CpuOperatorWrapper(op; array_type = CuArray{Float32})  # specify GPU array type for buffers
+wrapper = CpuOperatorWrapper(op)
 
 x_gpu = CUDA.randn(Float32, 64)
 y_gpu = similar(x_gpu, 63)
@@ -174,3 +170,9 @@ y_gpu = similar(x_gpu, 63)
 mul!(y_gpu, wrapper, x_gpu)   # GPU in → CPU compute → GPU out
 mul!(x_gpu, wrapper', y_gpu)  # GPU in → CPU adjoint → GPU out
 ```
+
+The wrapper preallocates CPU buffers (`dom_buf`, `cod_buf`) to avoid allocations during `mul!`. For parallel use, create independent copies per thread:
+
+```julia
+wrappers = [CpuOperatorWrapper(op) for _ in 1:Threads.nthreads()]
+```
diff --git a/ext/GpuExt/AGENTS.md b/ext/GpuExt/AGENTS.md
new file mode 100644
index 0000000..87e5fe7
--- /dev/null
+++ b/ext/GpuExt/AGENTS.md
@@ -0,0 +1,20 @@
+# GPU Implementation
+
+- Julia package extensions can only `import` the parent package, trigger package(s), and
+  stdlib; if extension code needs a parent dependency API, expose it from the parent module
+  first.
+- Override `mul!` in `ext/GpuExt/` for any operator whose base implementation uses scalar
+  indexing loops (`@nloops`, `@nref`, `@inbounds y[i] = b[j]`); replace with broadcast-over-view
+  (`y .= view(b, idx...)`).
+- When overriding a threaded operator (e.g. `Variation{..., true}`) for GPU, delegate to the
+  non-threaded variant (`Variation{..., false}`) — threading strategy is CPU-only.
+- For FFT plans, prefer `inv(plan)` (AbstractFFTs-generic) over backend-specific
+  `FFTW.plan_inv(...)` to keep CUDA/AMDGPU compatibility.
+- With JLArrays/GPUArrays, avoid `copyto!(gpu, cpu_view)` where the source is a `SubArray`;
+  materialize first (e.g. `src[1:n]`), or copy from a plain array.
+- Keep CPU-only implementation details out of GPU overrides unless the backend truly supports
+  them.
+- For GPU `GetIndex` overrides, keep boolean-mask and integer-vector fancy indexing in CPU
+  paths unless the backend support is verified.
+- Prefer direct `CuArray(arr)` / `CUDA.zeros(...)` / `AMDGPU.ROCArray(arr)` / `AMDGPU.zeros(...)`
+  calls over intermediate conversion variables.
diff --git a/ext/GpuExt/CLAUDE.md b/ext/GpuExt/CLAUDE.md
new file mode 120000
index 0000000..47dc3e3
--- /dev/null
+++ b/ext/GpuExt/CLAUDE.md
@@ -0,0 +1 @@
+AGENTS.md
\ No newline at end of file
diff --git a/ext/GpuExt/linearoperators/getindex.jl b/ext/GpuExt/linearoperators/getindex.jl
index 805688d..fbcfb0c 100644
--- a/ext/GpuExt/linearoperators/getindex.jl
+++ b/ext/GpuExt/linearoperators/getindex.jl
@@ -3,10 +3,6 @@ function _to_gpu_indices(ref_array::AbstractGPUArray, cpu_idx::AbstractVector{<:
     return ArrayT(Vector{Int}(cpu_idx))
 end
 
-function _to_gpu_indices(ref_array::AbstractGPUArray, gpu_idx::AbstractGPUArray{<:Integer})
-    return gpu_idx
-end
-
 function _to_gpu_indices(array_type::Type, cpu_idx::AbstractVector{<:Integer})
     ArrayT = Base.typename(array_type).wrapper
     return ArrayT(Vector{Int}(cpu_idx))
diff --git a/src/AbstractOperators.jl b/src/AbstractOperators.jl
index 000608d..bbb7a29 100644
--- a/src/AbstractOperators.jl
+++ b/src/AbstractOperators.jl
@@ -60,6 +60,8 @@ include("linearoperators/LBFGS.jl")
 
 # Calculus rules
 
+# Calculus rules
+
 # Batch operators
 include("batching/BatchOp.jl")
 include("batching/SimpleBatchOp.jl")
diff --git a/src/batching/SpreadingBatchOp.jl b/src/batching/SpreadingBatchOp.jl
index 3a7ee32..fa91a30 100644
--- a/src/batching/SpreadingBatchOp.jl
+++ b/src/batching/SpreadingBatchOp.jl
@@ -109,7 +109,17 @@ function BatchOp(
         threaded::Bool = nthreads() > 1,
         threading_strategy::Symbol = ThreadingStrategy.AUTO,
     )
-    return BatchOp(operators, (); threaded, threading_strategy)
+    op_domain_dims = ndims(operators[1], 2)
+    op_codomain_dims = ndims(operators[1], 1)
+    spreading_dims = ndims(operators)
+    batch_domain_dim_mask = get_batch_dim_mask(op_domain_dims, spreading_dims, ())
+    batch_codomain_dim_mask = get_batch_dim_mask(op_codomain_dims, spreading_dims, ())
+    return BatchOp(
+        operators,
+        batch_domain_dim_mask => batch_codomain_dim_mask;
+        threaded,
+        threading_strategy,
+    )
 end
 
 function BatchOp(
diff --git a/src/calculus/BroadCast.jl b/src/calculus/BroadCast.jl
index 17c75f9..e8a420f 100644
--- a/src/calculus/BroadCast.jl
+++ b/src/calculus/BroadCast.jl
@@ -244,19 +244,19 @@ function permute(R::OperatorBroadCast{T, N, M, true}, p::AbstractVector{Int}) wh
 end
 
 function _copy_operator_impl(
-        op::NoOperatorBroadCast{T, N, M, Th, S}; array_type = nothing, threaded = nothing
+        op::NoOperatorBroadCast{T, N, M, Th, S}; storage_type = nothing, threaded = nothing
     ) where {T, N, M, Th, S}
     new_threaded = threaded === nothing ? Th : threaded
-    new_S = array_type === nothing ? S : array_type
+    new_S = storage_type === nothing ? S : storage_type
     return NoOperatorBroadCast(T, new_S, op.dim_in, op.reshaped_dim_in, op.dim_out; threaded = new_threaded)
 end
 
 function _copy_operator_impl(
-        op::OperatorBroadCast{T, N, M, Th}; array_type = nothing, threaded = nothing
+        op::OperatorBroadCast{T, N, M, Th}; storage_type = nothing, threaded = nothing
     ) where {T, N, M, Th}
     new_threaded = threaded === nothing ? Th : threaded
     inner_op = Th ? op.A[1] : op.A
-    new_op = copy_operator(inner_op; array_type, threaded)
+    new_op = copy_operator(inner_op; storage_type, threaded)
     return BroadCast(new_op, op.dim_out; threaded = new_threaded)
 end
 
diff --git a/src/calculus/Compose.jl b/src/calculus/Compose.jl
index ebe7d8e..39671f2 100644
--- a/src/calculus/Compose.jl
+++ b/src/calculus/Compose.jl
@@ -303,8 +303,8 @@ remove_displacement(C::Compose) = Compose(remove_displacement.(C.A), C.buf)
 
 get_operators(C::Compose) = C.A
 
-function _copy_operator_impl(op::Compose; array_type = nothing, threaded = nothing)
-    new_bufs = tuple([_convert_buffer(b, array_type) for b in op.buf]...)
-    new_ops = tuple([copy_operator(a; array_type, threaded) for a in op.A]...)
+function _copy_operator_impl(op::Compose; storage_type = nothing, threaded = nothing)
+    new_bufs = tuple([_convert_buffer(b, storage_type) for b in op.buf]...)
+    new_ops = tuple([copy_operator(a; storage_type, threaded) for a in op.A]...)
     return Compose(new_ops, new_bufs)
 end
diff --git a/src/calculus/DCAT.jl b/src/calculus/DCAT.jl
index 061b451..2accd14 100644
--- a/src/calculus/DCAT.jl
+++ b/src/calculus/DCAT.jl
@@ -72,14 +72,6 @@ end
     return first(t) == v ? i : _dcat_find_in_tuple(Base.tail(t), v, i + 1)
 end
 
-# Apply the inverse permutation encoded in idxs to natural, using only tuple
-# operations so that no Vector is allocated on the hot path.
-function _dcat_apply_invperm(natural::Tuple, idxs::Tuple)
-    N = length(natural)
-    p = _dcat_flatten_idxs(idxs)
-    return ntuple(j -> natural[_dcat_find_in_tuple(p, j)], Val(N))
-end
-
 # Constructors
 DCAT(A::AbstractOperator) = A
 
@@ -203,6 +195,14 @@ function get_normal_op(H::DCAT)
     return DCAT(tuple([get_normal_op(H.A[i]) for i in eachindex(H.A)]...), idxs, idxs)
 end
 
+# Apply inverse permutation encoded in idxs to natural, using only tuple
+# operations so that no Vector is allocated on the hot path.
+function _dcat_apply_invperm(natural::Tuple, idxs::Tuple)
+    N = length(natural)
+    p = _dcat_flatten_idxs(idxs)
+    return ntuple(j -> natural[_dcat_find_in_tuple(p, j)], Val(N))
+end
+
 # Properties
 Base.:(==)(H1::DCAT{N, L1, P1, P2}, H2::DCAT{N, L2, P1, P2}) where {N, L1, L2, P1, P2} = H1.A == H2.A && H1.idxD == H2.idxD && H1.idxC == H2.idxC
 
diff --git a/src/calculus/HCAT.jl b/src/calculus/HCAT.jl
index ba1ce10..d5e0e86 100644
--- a/src/calculus/HCAT.jl
+++ b/src/calculus/HCAT.jl
@@ -327,12 +327,7 @@ is_sliced(L::HCAT) = any(is_sliced.(L.A))
 function get_slicing_expr(L::HCAT)
     exprs = ()
     for i in eachindex(L.A)
-        expr = get_slicing_expr(L[i])
-        if expr isa Tuple && all(e -> e isa Tuple, expr)
-            exprs = (exprs..., expr...)
-        else
-            exprs = (exprs..., expr)
-        end
+        exprs = (exprs..., get_slicing_expr(L.A[i]))
     end
     if length(exprs) == 1
         return exprs[1]
@@ -364,8 +359,8 @@ end
 
 remove_displacement(H::HCAT) = HCAT(remove_displacement.(H.A), H.idxs, H.buf)
 
-function _copy_operator_impl(op::HCAT; array_type = nothing, threaded = nothing)
-    new_buf = _convert_buffer(op.buf, array_type)
-    new_ops = tuple([copy_operator(a; array_type, threaded) for a in op.A]...)
+function _copy_operator_impl(op::HCAT; storage_type = nothing, threaded = nothing)
+    new_buf = _convert_buffer(op.buf, storage_type)
+    new_ops = tuple([copy_operator(a; storage_type, threaded) for a in op.A]...)
     return HCAT(new_ops, op.idxs, new_buf)
 end
diff --git a/src/calculus/HadamardProd.jl b/src/calculus/HadamardProd.jl
index 712b929..8238cfb 100644
--- a/src/calculus/HadamardProd.jl
+++ b/src/calculus/HadamardProd.jl
@@ -112,21 +112,21 @@ function remove_displacement(P::HadamardProd)
     )
 end
 
-function _copy_operator_impl(op::HadamardProd; array_type = nothing, threaded = nothing)
-    new_bufA = _convert_buffer(op.bufA, array_type)
-    new_bufB = _convert_buffer(op.bufB, array_type)
-    new_bufD = _convert_buffer(op.bufD, array_type)
-    new_A = copy_operator(op.A; array_type, threaded)
-    new_B = copy_operator(op.B; array_type, threaded)
+function _copy_operator_impl(op::HadamardProd; storage_type = nothing, threaded = nothing)
+    new_bufA = _convert_buffer(op.bufA, storage_type)
+    new_bufB = _convert_buffer(op.bufB, storage_type)
+    new_bufD = _convert_buffer(op.bufD, storage_type)
+    new_A = copy_operator(op.A; storage_type, threaded)
+    new_B = copy_operator(op.B; storage_type, threaded)
     return HadamardProd(new_A, new_B, new_bufA, new_bufB, new_bufD)
 end
 
-function _copy_operator_impl(op::HadamardProdJac; array_type = nothing, threaded = nothing)
-    new_bufA = _convert_buffer(op.bufA, array_type)
-    new_bufB = _convert_buffer(op.bufB, array_type)
-    new_bufD = _convert_buffer(op.bufD, array_type)
-    new_A = copy_operator(op.A; array_type, threaded)
-    new_B = copy_operator(op.B; array_type, threaded)
+function _copy_operator_impl(op::HadamardProdJac; storage_type = nothing, threaded = nothing)
+    new_bufA = _convert_buffer(op.bufA, storage_type)
+    new_bufB = _convert_buffer(op.bufB, storage_type)
+    new_bufD = _convert_buffer(op.bufD, storage_type)
+    new_A = copy_operator(op.A; storage_type, threaded)
+    new_B = copy_operator(op.B; storage_type, threaded)
     return HadamardProdJac{typeof(new_A), typeof(new_B), typeof(new_bufA), typeof(new_bufD)}(
         new_A, new_B, new_bufA, new_bufB, new_bufD
     )
diff --git a/src/calculus/OperatorWrapper.jl b/src/calculus/OperatorWrapper.jl
index e7a1e26..600f829 100644
--- a/src/calculus/OperatorWrapper.jl
+++ b/src/calculus/OperatorWrapper.jl
@@ -61,8 +61,6 @@ function OperatorWrapper(op::AbstractOperator; array_type::Type = Array)
     S = _array_wrapper_type(array_type)
     T_dom = domain_type(op)
     T_cod = codomain_type(op)
-    N_dom = ndims(dom_buf)
-    N_cod = ndims(cod_buf)
     DS = S{T_dom}
     CS = S{T_cod}
     return OperatorWrapper{typeof(op), typeof(dom_buf), typeof(cod_buf), DS, CS}(op, dom_buf, cod_buf)
@@ -127,9 +125,9 @@ displacement(A::OperatorWrapper) = displacement(A.op)
 remove_displacement(A::OperatorWrapper) = OperatorWrapper(remove_displacement(A.op))
 
 function _copy_operator_impl(
-        A::OperatorWrapper{Op, DB, CB, DS, CS}; array_type = nothing, threaded = nothing
+        A::OperatorWrapper{Op, DB, CB, DS, CS}; storage_type = nothing, threaded = nothing
     ) where {Op, DB, CB, DS, CS}
-    new_op = copy_operator(A.op; array_type = nothing, threaded)
+    new_op = copy_operator(A.op; storage_type = nothing, threaded)
     return OperatorWrapper{typeof(new_op), DB, CB, DS, CS}(
         new_op, similar(A.dom_buf), similar(A.cod_buf)
     )
diff --git a/src/calculus/Sum.jl b/src/calculus/Sum.jl
index e51f373..02eec7d 100644
--- a/src/calculus/Sum.jl
+++ b/src/calculus/Sum.jl
@@ -186,10 +186,10 @@ end
 
 remove_displacement(S::Sum) = Sum(remove_displacement.(S.A), S.bufC, S.bufD)
 
-function _copy_operator_impl(op::Sum; array_type = nothing, threaded = nothing)
-    new_bufC = _convert_buffer(op.bufC, array_type)
-    new_bufD = _convert_buffer(op.bufD, array_type)
-    new_ops = tuple([copy_operator(a; array_type, threaded) for a in op.A]...)
+function _copy_operator_impl(op::Sum; storage_type = nothing, threaded = nothing)
+    new_bufC = _convert_buffer(op.bufC, storage_type)
+    new_bufD = _convert_buffer(op.bufD, storage_type)
+    new_ops = tuple([copy_operator(a; storage_type, threaded) for a in op.A]...)
     K = length(new_ops)
     L = typeof(new_ops)
     return Sum{K, typeof(new_bufC), typeof(new_bufD), L}(new_ops, new_bufC, new_bufD)
diff --git a/src/calculus/VCAT.jl b/src/calculus/VCAT.jl
index 6579b7f..5895cdb 100644
--- a/src/calculus/VCAT.jl
+++ b/src/calculus/VCAT.jl
@@ -35,12 +35,7 @@ struct VCAT{
         CS <: AbstractArray,  # codomain storage type (fixed at construction)
     } <: AbstractOperator
     A::L     # tuple of AbstractOperators
-    idxs::P  # indices
-    # H = VCAT(Eye(n),VCAT(Eye(n),Eye(n))) has H.idxs = (1,2,3)
-    # `AbstractOperators` are flatten
-    # H = VCAT(Eye(n),Compose(MatrixOp(randn(n,n)),VCAT(Eye(n),Eye(n))))
-    # has H.idxs = (1,(2,3))
-    # `AbstractOperators` are stack
+    idxs::P  # indices; always NTuple{N, Int} since inner VCATs are flattened at construction
     buf::C   # buffer memory
     function VCAT(
             A::L, idxs::P, buf::C
@@ -86,27 +81,16 @@ function VCAT(A::Vararg{AbstractOperator})
     return VCAT(AA, buf)
 end
 
-# compile-time codomain ndoms for VCAT's sub-operators
-_ndoms_from_type(::Type{<:VCAT{N}}, dim::Int) where {N} = dim == 1 ? N : 1
-
 @generated function VCAT(AA::NTuple{N, AbstractOperator}, buf) where {N}
-    N == 1 && return :(AA[1])
-    # Build idxs at compile time using operator element types
-    K = 0
-    idx_exprs = []
-    for i in 1:N
-        nd = _ndoms_from_type(fieldtype(AA, i), 1)
-        if nd == 1
-            K += 1
-            push!(idx_exprs, K)
-        else
-            K0 = K
-            push!(idx_exprs, ntuple(j -> K0 + j, nd))
-            K += nd
-        end
+    if N isa Int
+        N == 1 && return :(AA[1])
+        # Build idxs at compile time: inner VCATs are always flattened, so all elements have nd=1
+        idxs_literal = Expr(:tuple, (1:N)...)
+        return :(VCAT(AA, $idxs_literal, buf))
+    else
+        # N is not statically known (e.g. built up in a loop); fall back to runtime length
+        return :(VCAT(AA, ntuple(identity, length(AA)), buf))
     end
-    idxs_literal = Expr(:tuple, idx_exprs...)
-    return :(VCAT(AA, $idxs_literal, buf))
 end
 
 VCAT(A::AbstractOperator) = A
@@ -116,17 +100,8 @@ VCAT(A::AbstractOperator) = A
 @generated function mul!(y::ArrayPartition, H::VCAT{N, L, P}, b::AbstractArray) where {N, L, P}
     ex = :(check(y, H, b))
     for i in 1:N
-        if fieldtype(P, i) <: Int
-            # flatten operator
-            # build mul!(y.x[H.idxs[i]], H.A[i], b)
-            yy = :(y.x[H.idxs[$i]])
-        else
-            # stacked operator
-            # build mul!(ArrayPartition( y[.xH.idxs[i][1]], y.x[H.idxs[i][2]] ...  ), H.A[i], b)
-            yy = [:(y.x[H.idxs[$i][$ii]]) for ii in eachindex(fieldnames(fieldtype(P, i)))]
-            yy = :(ArrayPartition($(yy...)))
-        end
-        ex = :($ex; mul!($yy, H.A[$i], b))
+        # P always has Int elements (inner VCATs are flattened at construction)
+        ex = :($ex; mul!(y.x[H.idxs[$i]], H.A[$i], b))
     end
     ex = :($ex; return y)
     return ex
@@ -137,30 +112,11 @@ end
     ) where {N, L, P}
     ex = :(check(y, A, b); H = A.A)
 
-    if fieldtype(P, 1) <: Int
-        # flatten operator
-        # build mul!(y, H.A[1]', b.x[H.idxs[1]])
-        bb = :(b.x[H.idxs[1]])
-    else
-        # stacked operator
-        # build mul!(y, H.A[1]',ArrayPartition( b.x[H.idxs[1][1]], b.x[H.idxs[1][2]] ...  ))
-        bb = [:(b.x[H.idxs[1][$ii]]) for ii in eachindex(fieldnames(fieldtype(P, 1)))]
-        bb = :(ArrayPartition($(bb...)))
-    end
-    ex = :($ex; mul!(y, H.A[1]', $bb)) # write on y
+    # P always has Int elements (inner VCATs are flattened at construction)
+    ex = :($ex; mul!(y, H.A[1]', b.x[H.idxs[1]])) # write on y
 
     for i in 2:N
-        if fieldtype(P, i) <: Int
-            # flatten operator
-            # build mul!(H.buf, H.A[i]', b.x[H.idxs[i]])
-            bb = :(b.x[H.idxs[$i]])
-        else
-            # stacked operator
-            # build mul!(H.buf, H.A[i]',( b.x[H.idxs[i][1]], b.x[H.idxs[i][2]] ...  ))
-            bb = [:(b.x[H.idxs[$i][$ii]]) for ii in eachindex(fieldnames(fieldtype(P, i)))]
-            bb = :(ArrayPartition($(bb...)))
-        end
-        ex = :($ex; mul!(H.buf, H.A[$i]', $bb)) # write on H.buf
+        ex = :($ex; mul!(H.buf, H.A[$i]', b.x[H.idxs[$i]])) # write on H.buf
         # sum H.buf with y
         ex = :($ex; y .+= H.buf)
     end
@@ -175,17 +131,8 @@ function Base.:(==)(H1::VCAT{N, L1, P1, C}, H2::VCAT{N, L2, P2, C}) where {N, L1
 end
 
 @generated function size(H::VCAT{N, L, P}) where {N, L, P}
-    exprs = []
-    for i in 1:N
-        Pi = fieldtype(P, i)
-        if Pi <: Integer
-            push!(exprs, :(size(H.A[$i], 1)))
-        else
-            for ii in eachindex(fieldnames(Pi))
-                push!(exprs, :(size(H.A[$i], 1)[$ii]))
-            end
-        end
-    end
+    # P always has Int elements (inner VCATs are flattened at construction)
+    exprs = [:(size(H.A[$i], 1)) for i in 1:N]
     natural_expr = Expr(:tuple, exprs...)
     return :((_vcat_apply_invperm($natural_expr, H.idxs), size(H.A[1], 2)))
 end
@@ -203,17 +150,8 @@ end
 
 domain_type(L::VCAT) = domain_type.(Ref(L.A[1]))
 @generated function codomain_type(H::VCAT{N, L, P}) where {N, L, P}
-    exprs = []
-    for i in 1:N
-        Pi = fieldtype(P, i)
-        if Pi <: Integer
-            push!(exprs, :(codomain_type(H.A[$i])))
-        else
-            for ii in eachindex(fieldnames(Pi))
-                push!(exprs, :(codomain_type(H.A[$i])[$ii]))
-            end
-        end
-    end
+    # P always has Int elements (inner VCATs are flattened at construction)
+    exprs = [:(codomain_type(H.A[$i])) for i in 1:N]
     natural_expr = Expr(:tuple, exprs...)
     return :(_vcat_apply_invperm($natural_expr, H.idxs))
 end
@@ -230,7 +168,7 @@ function get_slicing_expr(L::VCAT)
     return get_slicing_expr.(Tuple(L.A[i] for i in eachindex(L.A)))
 end
 function remove_slicing(L::VCAT)
-    new_ops = remove_slicing.(L[i] for i in eachindex(L.A))
+    new_ops = collect(map(remove_slicing, L.A))
     if !any(a -> a isa HCAT, new_ops) && all(i -> i isa Int, L.idxs)
         return DCAT(new_ops[collect(L.idxs)]...)
     elseif all(a -> a isa HCAT, L.A) && any(a -> any(is_null, a.A), L.A) && any(op -> size(op, 2) != size(new_ops[1], 2), new_ops)
@@ -284,8 +222,8 @@ end
 
 remove_displacement(V::VCAT) = VCAT(remove_displacement.(V.A), V.idxs, V.buf)
 
-function _copy_operator_impl(op::VCAT; array_type = nothing, threaded = nothing)
-    new_buf = _convert_buffer(op.buf, array_type)
-    new_ops = tuple([copy_operator(a; array_type, threaded) for a in op.A]...)
+function _copy_operator_impl(op::VCAT; storage_type = nothing, threaded = nothing)
+    new_buf = _convert_buffer(op.buf, storage_type)
+    new_ops = tuple([copy_operator(a; storage_type, threaded) for a in op.A]...)
     return VCAT(new_ops, op.idxs, new_buf)
 end
diff --git a/src/combination_rules.jl b/src/combination_rules.jl
index b05aad2..9a0dd19 100644
--- a/src/combination_rules.jl
+++ b/src/combination_rules.jl
@@ -284,8 +284,6 @@ end
 combine_matrix(L::AbstractMatrix, R::AbstractMatrix) = L * R
 combine_matrix(L::AbstractMatrix, R::AbstractVector) = L * Diagonal(R)
 combine_matrix(L::AbstractVector, R::AbstractMatrix) = Diagonal(L) * R
-combine_matrix(L::Number, R::AbstractMatrix) = L * R
-combine_matrix(L::AbstractMatrix, R::Number) = R * L
 function combine(T1::DiagOp, T2::MatrixOp)
     return MatrixOp(domain_type(T2), size(T2, 2), combine_matrix(T1.d, T2.A))
 end
diff --git a/src/linearoperators/DiagOp.jl b/src/linearoperators/DiagOp.jl
index f24fc80..651c6e2 100644
--- a/src/linearoperators/DiagOp.jl
+++ b/src/linearoperators/DiagOp.jl
@@ -111,10 +111,10 @@ domain_type(::DiagOp{<:Any, D}) where {D} = D
 codomain_type(::DiagOp{<:Any, <:Any, C}) where {C} = C
 is_thread_safe(::DiagOp) = true
 
-function _copy_operator_impl(op::DiagOp{B}; array_type = nothing, threaded = nothing) where {B}
+function _copy_operator_impl(op::DiagOp{B}; storage_type = nothing, threaded = nothing) where {B}
     new_threaded = threaded === nothing ? (B == FastBroadcast.True()) : threaded
-    new_d = array_type === nothing ? op.d : _convert_buffer(op.d, array_type)
-    new_at = array_type === nothing ? _array_wrapper(op.d) : array_type
+    new_d = storage_type === nothing ? op.d : _convert_buffer(op.d, storage_type)
+    new_at = storage_type === nothing ? _array_wrapper(op.d) : storage_type
     return DiagOp(domain_type(op), op.dim_in, new_d; threaded = new_threaded, array_type = new_at)
 end
 
diff --git a/src/linearoperators/Variation.jl b/src/linearoperators/Variation.jl
index a067e8a..6ea2658 100644
--- a/src/linearoperators/Variation.jl
+++ b/src/linearoperators/Variation.jl
@@ -214,10 +214,10 @@ codomain_array_type(::Variation{T, N, Th, S}) where {T, N, Th, S} = S
 is_thread_safe(::Variation) = true
 
 function _copy_operator_impl(
-        op::Variation{T, N, Th, S}; array_type = nothing, threaded = nothing
+        op::Variation{T, N, Th, S}; storage_type = nothing, threaded = nothing
     ) where {T, N, Th, S}
     new_threaded = threaded === nothing ? Th : threaded
-    new_at = array_type === nothing ? _array_wrapper_type(S) : array_type
+    new_at = storage_type === nothing ? _array_wrapper_type(S) : storage_type
     return Variation(T, op.dim_in; threaded = new_threaded, array_type = new_at)
 end
 
diff --git a/src/linearoperators/Zeros.jl b/src/linearoperators/Zeros.jl
index 8628191..5c6b0de 100644
--- a/src/linearoperators/Zeros.jl
+++ b/src/linearoperators/Zeros.jl
@@ -60,7 +60,7 @@ function Zeros(
         codomain_type::NTuple{NN, Type},
         dim_out::NTuple{NN, Tuple},
     ) where {NN}
-    return VCAT([Zeros(domain_type, dim_in, codomain_type[i], dim_out[i]) for i in 1:NN]...)
+    return VCAT(ntuple(i -> Zeros(domain_type, dim_in, codomain_type[i], dim_out[i]), Val(NN))...)
 end
 
 Zeros(A::AbstractOperator) = Zeros(domain_type(A), size(A, 2), codomain_type(A), size(A, 1))
diff --git a/src/properties.jl b/src/properties.jl
index 0d878a0..b2d7bbc 100644
--- a/src/properties.jl
+++ b/src/properties.jl
@@ -74,7 +74,7 @@ RecursiveArrayTools.ArrayPartition{ComplexF64, Tuple{Array{ComplexF64}, Array{Co
 ```
 """
 function domain_array_type(L::AbstractOperator)
-    return _array_type_for_elem(domain_type(L))
+    return _storage_type_for_elem(domain_type(L))
 end
 
 """
@@ -91,11 +91,11 @@ RecursiveArrayTools.ArrayPartition{ComplexF64, Tuple{Array{ComplexF64}, Array{Co
 ```
 """
 function codomain_array_type(L::AbstractOperator)
-    return _array_type_for_elem(codomain_type(L))
+    return _storage_type_for_elem(codomain_type(L))
 end
 
-_array_type_for_elem(T::Type) = Array{T}
-function _array_type_for_elem(dt::Tuple)
+_storage_type_for_elem(T::Type) = Array{T}
+function _storage_type_for_elem(dt::Tuple)
     arrayTypes = Tuple{[Array{t} for t in dt]...}
     return ArrayPartition{promote_type(dt...), arrayTypes}
 end
@@ -108,6 +108,7 @@ end
 function _normalize_array_type(array_type::Type{A}, elem_type::Type{T}) where {A <: AbstractArray, T}
     return _array_wrapper_type(A){T}
 end
+
 _storage_eltype(::Type{<:AbstractArray{T}}) where {T} = T
 
 function allocate_in_domain(L::AbstractOperator, dims... = size(L, 2)...)
@@ -342,8 +343,6 @@ julia> AbstractOperators.combine(Eye(10), DiagOp(rand(10)))
 function combine(L, R)
     if is_eye(L)
         return R
-    elseif is_eye(R)
-        return L
     elseif is_null(L)
         if size(R, 1) == size(R, 2) && domain_type(R) == codomain_type(R)
             return L
@@ -486,28 +485,28 @@ function string_dom(dm::Tuple, sz::Tuple)
 end
 
 """
-    copy_operator(op::AbstractOperator; array_type=nothing, threaded=nothing)
+    copy_operator(op::AbstractOperator; storage_type=nothing, threaded=nothing)
 
 Create a copy of `op` suitable for parallel use.
 
 - Immutable fields (operator arrays, type params) are **shared** (no copy).
 - Mutable buffer fields are **deep-copied**.
-- `array_type`: if provided (e.g., `CuArray`), convert buffer arrays to that storage.
+- `storage_type`: if provided (e.g., `CuArray`), convert buffer arrays to that storage.
 - `threaded`: if provided (`true`/`false`), toggle threading for operators that support it.
 
-When `array_type` is `nothing` and `threaded` is `nothing`, equivalent to the old `copy_op`
+When `storage_type` is `nothing` and `threaded` is `nothing`, equivalent to the old `copy_op`
 but more efficient (shares immutable data).
 """
-function copy_operator(op::AbstractOperator; array_type = nothing, threaded = nothing)
-    if is_thread_safe(op) && threaded === nothing && array_type === nothing
+function copy_operator(op::AbstractOperator; storage_type = nothing, threaded = nothing)
+    if is_thread_safe(op) && threaded === nothing && storage_type === nothing
         return op  # safe to share
     end
-    return _copy_operator_impl(op; array_type, threaded)
+    return _copy_operator_impl(op; storage_type, threaded)
 end
 
 # Default implementation: just deepcopy (fallback)
 function _copy_operator_impl(
-        op::T; array_type = nothing, threaded = nothing
+        op::T; storage_type = nothing, threaded = nothing
     ) where {T <: AbstractOperator}
     return deepcopy(op)
 end
@@ -516,8 +515,8 @@ end
 function _convert_buffer(buf::AbstractArray, ::Nothing)
     return similar(buf)  # same type, new allocation
 end
-function _convert_buffer(buf::AbstractArray{T}, array_type::Type) where {T}
-    return similar(array_type{T}, size(buf))
+function _convert_buffer(buf::AbstractArray{T}, storage_type::Type) where {T}
+    return similar(storage_type{T}, size(buf))
 end
 
 _should_thread(::Number) = false
diff --git a/src/syntax.jl b/src/syntax.jl
index ad21e64..a9a0d2e 100644
--- a/src/syntax.jl
+++ b/src/syntax.jl
@@ -56,8 +56,6 @@ function Base.getindex(A::Compose, idx...)
     if ndoms(A, 2) == 1
         Gout = GetIndex(codomain_type(A), size(A, 1), idx)
         return Gout * A
-    elseif all(is_diagonal, A.A[2:end])
-        return Compose((getindex(A.A[1], idx...), A.A[2:end]...), A.buf)
     else
         error("cannot split operator of type $(typeof(A))")
     end
@@ -103,11 +101,7 @@ function Base.getindex(H::VCAT, idx::Union{AbstractArray, Int})
         for i in idx
             for ii in eachindex(H.idxs)
                 if i in H.idxs[ii]
-                    if typeof(H.idxs[ii]) <: Int
-                        new_H = (new_H..., H.A[ii])
-                    else
-                        error("cannot split operator: $H")
-                    end
+                    new_H = (new_H..., H.A[ii])
                 end
             end
         end
diff --git a/test/batching/test_SimpleBatchOp.jl b/test/batching/test_SimpleBatchOp.jl
index d50d69d..2600144 100644
--- a/test/batching/test_SimpleBatchOp.jl
+++ b/test/batching/test_SimpleBatchOp.jl
@@ -1,5 +1,5 @@
 @testmodule SimpleBatchOpHelpers begin
-    using Random, BenchmarkTools, LinearAlgebra, AbstractOperators, Test
+    using Random, BenchmarkTools, LinearAlgebra, AbstractOperators, JLArrays, Test
 
     function test_simple_batchop(op, batch_op, x, y, z, threaded)
         if threaded && Threads.nthreads() > 1
@@ -171,6 +171,52 @@ end
     end
 end
 
+@testitem "SimpleBatchOpMultiThreaded properties" tags = [:batching, :SimpleBatchOp] setup = [TestUtils] begin
+    using Random, LinearAlgebra, AbstractOperators
+    Random.seed!(0)
+    # Directly construct SimpleBatchOpMultiThreaded to test property/diag methods
+    # without requiring nthreads() > 1 at test time.
+    op = DiagOp([1.0, 2.0])
+    st = BatchOp(op, (2,); threaded = false)  # creates SimpleBatchOpSingleThreaded
+    @assert st isa AbstractOperators.SimpleBatchOpSingleThreaded
+    # Build MultiThreaded variant with same shape, 2 operator copies
+    mt = let T = typeof(st)
+        dT = T.parameters[1]
+        cT = T.parameters[2]
+        dM = T.parameters[3]
+        cM = T.parameters[4]
+        opT = typeof(op)
+        N = length(st.domain_size)
+        M = length(st.codomain_size)
+        C = 2
+        ops = (op, copy_operator(op))
+        AbstractOperators.SimpleBatchOpMultiThreaded{dT, cT, dM, cM, opT, N, M, C}(
+            ops, st.domain_size, st.codomain_size, CartesianIndices(st.batch_size)
+        )
+    end
+    @test diag_AAc(mt) == diag_AAc(st)
+    @test diag_AcA(mt) == diag_AcA(st)
+    @test diag(mt) == diag(st)
+    @test AbstractOperators.has_optimized_normalop(mt) == AbstractOperators.has_optimized_normalop(st)
+    @test opnorm(mt) == opnorm(st)
+    @test estimate_opnorm(mt) == estimate_opnorm(st)
+    # Eye operator: scalar diag paths
+    eye_op = Eye(Float64, (2,))
+    eye_st = BatchOp(eye_op, (2,); threaded = false)
+    eye_mt = let T = typeof(eye_st)
+        dT, cT, dM, cM = T.parameters[1], T.parameters[2], T.parameters[3], T.parameters[4]
+        opT = typeof(eye_op)
+        N, M, C = length(eye_st.domain_size), length(eye_st.codomain_size), 2
+        ops = (eye_op, copy_operator(eye_op))
+        AbstractOperators.SimpleBatchOpMultiThreaded{dT, cT, dM, cM, opT, N, M, C}(
+            ops, eye_st.domain_size, eye_st.codomain_size, CartesianIndices(eye_st.batch_size)
+        )
+    end
+    @test diag(eye_mt) == 1.0
+    @test diag_AcA(eye_mt) == 1.0
+    @test diag_AAc(eye_mt) == 1.0
+end
+
 @testitem "SimpleBatchOp benchmark" tags = [:batching, :SimpleBatchOp] setup = [TestUtils, SimpleBatchOpHelpers] begin
     using Random
     Random.seed!(0)
diff --git a/test/batching/test_SpreadingBatchOp.jl b/test/batching/test_SpreadingBatchOp.jl
index 02bfc7e..4fece5d 100644
--- a/test/batching/test_SpreadingBatchOp.jl
+++ b/test/batching/test_SpreadingBatchOp.jl
@@ -1,5 +1,5 @@
 @testmodule SpreadingBatchOpHelpers begin
-    using Random, BenchmarkTools, LinearAlgebra, AbstractOperators, Test
+    using Random, BenchmarkTools, LinearAlgebra, AbstractOperators, JLArrays, Test
 
     function test_spreading_batchop(operators, batch_op, x, y, z, threaded)
         if !threaded
@@ -319,3 +319,32 @@ end
         @test all(Array(y_gpu)[2, :, :] .≈ 2.0)
     end
 end
+
+@testitem "BatchOp without explicit sizes (lines 107, 112)" tags = [:batching, :SpreadingBatchOp] setup = [TestUtils, SpreadingBatchOpHelpers] begin
+    using Random, AbstractOperators
+    Random.seed!(0)
+    # BatchOp(operators) with no size args → calls BatchOp(operators, ()) → line 112
+    n = 5
+    ops = [Eye(Float64, (n,)) for _ in 1:3]
+    bop = BatchOp(ops; threaded = false)
+    @test bop isa AbstractOperators.SpreadingBatchOp
+    x = randn(n, 3)
+    y = bop * x
+    @test size(y) == (n, 3)
+    @test y ≈ x
+end
+
+@testitem "BatchOp unsupported threading strategy for non-thread-safe ops (line 404)" tags = [:batching, :SpreadingBatchOp] setup = [TestUtils, SpreadingBatchOpHelpers] begin
+    using Random, AbstractOperators
+    Random.seed!(0)
+    if Threads.nthreads() > 1
+        n = 5
+        # LBFGS is not thread-safe; an unknown strategy reaches the else branch at line 404
+        ops = [LBFGS(zeros(n), 3) for _ in 1:3]
+        @test_throws ArgumentError BatchOp(
+            ops, (4,), (:b, :s, :_);
+            threaded = true,
+            threading_strategy = :UNKNOWN_STRATEGY,
+        )
+    end
+end
diff --git a/test/calculus/test_Ax_mul_Bx.jl b/test/calculus/test_Ax_mul_Bx.jl
index 1b0e987..80a2153 100644
--- a/test/calculus/test_Ax_mul_Bx.jl
+++ b/test/calculus/test_Ax_mul_Bx.jl
@@ -101,7 +101,7 @@
     @test Jacobian(Ax_mul_Bx(A, B), x) == Jacobian(Ax_mul_Bx(A, B), x)
 end
 
-@testitem "Ax_mul_Bx (GPU)" tags = [:gpu, :calculus, :Ax_mul_Bx] setup = [TestUtils] begin
+@testitem "Ax_mul_Bx (GPU)" tags = [:gpu, :calculus, :Ax_mul_Bx] setup = [TestUtils, GPUNLTestUtils] begin
     using Random, AbstractOperators, GPUEnv
 
     for backend in gpu_backends()
@@ -117,7 +117,8 @@ end
         test_NLop_gpu(P, x, r, false)
 
         n2 = 3
-        P2 = Ax_mul_Bx(Sin(gpu_zeros(backend, Float64, n2, n2)), Cos(gpu_zeros(backend, Float64, n2, n2)))
+        AT = gpu_wrapper(backend, Float64, n2, n2)
+        P2 = Ax_mul_Bx(Sin(Float64, (n2, n2); array_type = AT), Cos(Float64, (n2, n2); array_type = AT))
         x2 = gpu_randn(backend, n2, n2)
         r2 = gpu_randn(backend, n2, n2)
         test_NLop_gpu(P2, x2, r2, false)
diff --git a/test/calculus/test_Ax_mul_Bxt.jl b/test/calculus/test_Ax_mul_Bxt.jl
index 45fdfc9..05adf76 100644
--- a/test/calculus/test_Ax_mul_Bxt.jl
+++ b/test/calculus/test_Ax_mul_Bxt.jl
@@ -1,6 +1,7 @@
 @testitem "Ax_mul_Bxt: basic mul" tags = [:calculus, :Ax_mul_Bxt] setup = [TestUtils] begin
     using Random, AbstractOperators
     Random.seed!(0)
+    verb && println(" --- Testing Ax_mul_Bxt: basic mul --- ")
 
     n = 10
     A, B = Eye(n), Sin(n)
@@ -38,6 +39,7 @@ end
 @testitem "Ax_mul_Bxt: HCAT and permute" tags = [:calculus, :Ax_mul_Bxt] setup = [TestUtils] begin
     using Random, AbstractOperators
     Random.seed!(0)
+    verb && println(" --- Testing Ax_mul_Bxt: HCAT and permute --- ")
 
     # testing with HCAT
     m, n = 3, 5
@@ -72,6 +74,7 @@ end
 @testitem "Ax_mul_Bxt: error paths and equality" tags = [:calculus, :Ax_mul_Bxt] setup = [TestUtils] begin
     using Random, AbstractOperators
     Random.seed!(0)
+    verb && println(" --- Testing Ax_mul_Bxt: error paths and equality --- ")
 
     # ndims==2 branch with mismatched second codomain dimension
     struct AxDummy2D <: AbstractOperator
@@ -109,7 +112,7 @@ end
     @test Jacobian(Ax_mul_Bxt(A, B), x) == Jacobian(Ax_mul_Bxt(A, B), x)
 end
 
-@testitem "Ax_mul_Bxt (GPU)" tags = [:gpu, :calculus, :Ax_mul_Bxt] setup = [TestUtils] begin
+@testitem "Ax_mul_Bxt (GPU)" tags = [:gpu, :calculus, :Ax_mul_Bxt] setup = [TestUtils, GPUNLTestUtils] begin
     using Random, AbstractOperators, GPUEnv
 
     for backend in gpu_backends()
@@ -118,10 +121,18 @@ end
         n = 10
         P = Ax_mul_Bxt(
             Eye(Float64, (n,); array_type = gpu_wrapper(backend, Float64, n)),
-            Sin(gpu_zeros(backend, Float64, n)),
+            Sin(Float64, (n,); array_type = gpu_wrapper(backend, Float64, n)),
         )
         x = gpu_randn(backend, n)
         r = gpu_randn(backend, n, n)
         test_NLop_gpu(P, x, r, false)
     end
 end
+
+@testitem "Ax_mul_Bxt: size mismatch error (1D)" tags = [:calculus, :Ax_mul_Bxt] setup = [TestUtils] begin
+    using AbstractOperators
+    # Two 1D operators with different sizes -> triggers size(A) != size(B) in inner constructor
+    A = MatrixOp(randn(3, 4))   # size = ((3,), (4,))
+    B = MatrixOp(randn(5, 4))   # size = ((5,), (4,)) -- different codomain
+    @test_throws DimensionMismatch Ax_mul_Bxt(A, B)
+end
diff --git a/test/calculus/test_Axt_mul_Bx.jl b/test/calculus/test_Axt_mul_Bx.jl
index 1546321..7192e87 100644
--- a/test/calculus/test_Axt_mul_Bx.jl
+++ b/test/calculus/test_Axt_mul_Bx.jl
@@ -1,6 +1,7 @@
 @testitem "Axt_mul_Bx: basic mul" tags = [:calculus, :Axt_mul_Bx] setup = [TestUtils] begin
     using Random, AbstractOperators
     Random.seed!(0)
+    verb && println(" --- Testing Axt_mul_Bx: basic mul --- ")
 
     n = 10
     A, B = Eye(n), Sin(n)
@@ -40,6 +41,7 @@ end
 @testitem "Axt_mul_Bx: HCAT and permute" tags = [:calculus, :Axt_mul_Bx] setup = [TestUtils] begin
     using Random, AbstractOperators
     Random.seed!(0)
+    verb && println(" --- Testing Axt_mul_Bx: HCAT and permute --- ")
 
     # testing with HCAT
     m, n = 3, 5
@@ -74,6 +76,7 @@ end
 @testitem "Axt_mul_Bx: error paths and equality" tags = [:calculus, :Axt_mul_Bx] setup = [TestUtils] begin
     using Random, AbstractOperators
     Random.seed!(0)
+    verb && println(" --- Testing Axt_mul_Bx: error paths and equality --- ")
 
     # ndims==2 branch with mismatched first codomain dimension
     struct AxtDummy2D <: AbstractOperator
@@ -139,3 +142,28 @@ end
         @test Array(y2) ≈ Ax' * Bx
     end
 end
+
+@testitem "Axt_mul_Bx 2D operator DimensionMismatch" tags = [:calculus, :Axt_mul_Bx] setup = [TestUtils] begin
+    using Random, AbstractOperators
+    Random.seed!(0)
+    # 2D operators with same domain but different first codomain dimension (line 39)
+    A = MatrixOp(randn(3, 5), 4)
+    B = MatrixOp(randn(6, 5), 4)
+    @test_throws DimensionMismatch Axt_mul_Bx(A, B)
+end
+
+@testitem "Axt_mul_Bx: size mismatch error (1D)" tags = [:calculus, :Axt_mul_Bx] setup = [TestUtils] begin
+    using AbstractOperators
+    # Two 1D operators with different sizes -> triggers size(A) != size(B) in inner constructor
+    A = MatrixOp(randn(3, 4))   # size = ((3,), (4,))
+    B = MatrixOp(randn(5, 4))   # size = ((5,), (4,)) -- different codomain
+    @test_throws DimensionMismatch Axt_mul_Bx(A, B)
+end
+
+@testitem "Axt_mul_Bx: size mismatch error (1D)" tags = [:calculus, :Axt_mul_Bx] setup = [TestUtils] begin
+    using AbstractOperators
+    # Two 1D operators with different sizes -> triggers size(A) != size(B) in inner constructor
+    A = MatrixOp(randn(3, 4))   # size = ((3,), (4,))
+    B = MatrixOp(randn(5, 4))   # size = ((5,), (4,)) -- different codomain
+    @test_throws DimensionMismatch Axt_mul_Bx(A, B)
+end
diff --git a/test/calculus/test_Jacobian.jl b/test/calculus/test_Jacobian.jl
index a357302..03f91e6 100644
--- a/test/calculus/test_Jacobian.jl
+++ b/test/calculus/test_Jacobian.jl
@@ -1,5 +1,6 @@
 @testitem "Jacobian: basic HCAT" tags = [:calculus, :Jacobian] setup = [TestUtils] begin
     using AbstractOperators
+    verb && println(" --- Testing Jacobian: basic HCAT --- ")
 
     m, n = 3, 5
     x = ArrayPartition(randn(m), randn(n))
@@ -19,6 +20,7 @@ end
 
 @testitem "Jacobian: LinearOperator and Scale paths" tags = [:calculus, :Jacobian] setup = [TestUtils] begin
     using AbstractOperators
+    verb && println(" --- Testing Jacobian: LinearOperator and Scale paths --- ")
 
     # 1. LinearOperator path (Jacobian of a LinearOperator returns itself)
     n_lin = 6
@@ -40,6 +42,7 @@ end
 
 @testitem "Jacobian: AffineAdd and Transpose paths" tags = [:calculus, :Jacobian] setup = [TestUtils] begin
     using AbstractOperators
+    verb && println(" --- Testing Jacobian: AffineAdd and Transpose paths --- ")
 
     # 3. AffineAdd path (Jacobian should drop displacement)
     n_aff = 4
@@ -58,6 +61,7 @@ end
 
 @testitem "Jacobian: Compose Sum VCAT paths" tags = [:calculus, :Jacobian] setup = [TestUtils] begin
     using AbstractOperators
+    verb && println(" --- Testing Jacobian: Compose Sum VCAT paths --- ")
 
     # 5. Compose path (single op) with tuple input to trigger second Compose method
     n_cp = 3
@@ -93,6 +97,7 @@ end
 
 @testitem "Jacobian: HCAT and DCAT paths" tags = [:calculus, :Jacobian] setup = [TestUtils] begin
     using AbstractOperators
+    verb && println(" --- Testing Jacobian: HCAT and DCAT paths --- ")
 
     # 8. HCAT path with multi-index block (length(idx) > 1) to cover else branch
     m_h1, m_h2 = 3, 2
@@ -121,6 +126,7 @@ end
 
 @testitem "Jacobian: Reshape and BroadCast paths" tags = [:calculus, :Jacobian] setup = [TestUtils] begin
     using AbstractOperators
+    verb && println(" --- Testing Jacobian: Reshape and BroadCast paths --- ")
 
     # 10. Reshape path
     n_rs = 6
@@ -142,6 +148,7 @@ end
 
 @testitem "Jacobian: equality and properties" tags = [:calculus, :Jacobian] setup = [TestUtils] begin
     using AbstractOperators
+    verb && println(" --- Testing Jacobian: equality and properties --- ")
 
     # 12. Equality and show output
     n_eq = 4
diff --git a/test/calculus/test_adjointoperator.jl b/test/calculus/test_adjointoperator.jl
index 1e6f64c..d38d9dd 100644
--- a/test/calculus/test_adjointoperator.jl
+++ b/test/calculus/test_adjointoperator.jl
@@ -96,7 +96,7 @@ end
         Random.seed!(0)
 
         n = 5
-        op = FiniteDiff(gpu_zeros(backend, Float64, n))
+        op = FiniteDiff(Float64, (n,); array_type = gpu_wrapper(backend, Float64, n))
         opT = AdjointOperator(op)
         test_op(opT, gpu_randn(backend, n - 1), gpu_randn(backend, n), false)
 
diff --git a/test/calculus/test_affineadd.jl b/test/calculus/test_affineadd.jl
index 8ff4276..b01442e 100644
--- a/test/calculus/test_affineadd.jl
+++ b/test/calculus/test_affineadd.jl
@@ -147,3 +147,21 @@ end
         @test collect(r_adj) ≈ collect(r_adj2)
     end
 end
+
+@testitem "AffineAdd: array type mismatch error" tags = [:calculus, :AffineAdd] setup = [TestUtils] begin
+    using AbstractOperators
+    n = 5
+    op = Eye(Float64, (n,))
+    # eltype(d) != codomain_type(op): ComplexF64 vs Float64 (line 39)
+    @test_throws ErrorException AffineAdd(op, randn(ComplexF64, n))
+    # Float32 vs Float64
+    @test_throws ErrorException AffineAdd(op, Float32.(randn(n)))
+end
+
+@testitem "AffineAdd: element type mismatch error" tags = [:calculus, :AffineAdd] setup = [TestUtils] begin
+    using AbstractOperators
+    n = 4
+    op = MatrixOp(randn(n, n))    # codomain_type = Float64
+    d = randn(Float32, n)          # eltype = Float32 != Float64
+    @test_throws ErrorException AffineAdd(op, d)
+end
diff --git a/test/calculus/test_broadcast.jl b/test/calculus/test_broadcast.jl
index 70c0ed1..ceb23d7 100644
--- a/test/calculus/test_broadcast.jl
+++ b/test/calculus/test_broadcast.jl
@@ -1,6 +1,7 @@
 @testitem "BroadCast: basic mul" tags = [:calculus, :BroadCast] setup = [TestUtils] begin
     using Random, AbstractOperators
     Random.seed!(0)
+    verb && println(" --- Testing BroadCast --- ")
 
     m, n = 8, 4
     dim_out = (m, 10)
@@ -180,3 +181,65 @@ end
         test_op(opR2, gpu_randn(backend, m2, n2), gpu_randn(backend, dim_out2...), false)
     end
 end
+
+@testitem "BroadCast same-size returns operator unchanged" tags = [:calculus, :BroadCast] setup = [TestUtils] begin
+    using Random, AbstractOperators
+    Random.seed!(0)
+    m, n = 5, 3
+    A = MatrixOp(randn(m, n))
+    # BroadCast with dim_out == size(A, 1) should return A unchanged (line 80)
+    result = BroadCast(A, size(A, 1))
+    @test result === A
+    # Same for 2D Eye
+    B = Eye(m, n)
+    result2 = BroadCast(B, size(B, 1))
+    @test result2 === B
+end
+
+@testitem "BroadCast non-compact adjoint reshape" tags = [:calculus, :BroadCast] setup = [TestUtils] begin
+    using Random, AbstractOperators, LinearAlgebra
+    Random.seed!(0)
+    m, n = 3, 2
+    # Reshape codomain to (1, m) so slices of dim_out=(4,m,5) won't match (line 144)
+    A = reshape(MatrixOp(randn(m, n)), 1, m)
+    dim_out = (4, m, 5)
+    B_noncompact = BroadCast(A, dim_out; threaded = false)
+    x = randn(n)
+    y = B_noncompact * x
+    @test size(y) == dim_out
+    y_test = randn(dim_out)
+    x_back = B_noncompact' * y_test
+    @test size(x_back) == (n,)
+end
+
+@testitem "BroadCast: copy_operator" tags = [:calculus, :BroadCast] setup = [TestUtils] begin
+    using Random, AbstractOperators, LinearAlgebra
+    Random.seed!(7)
+
+    m, n = 8, 4
+    dim_out = (m, 10)
+
+    # NoOperatorBroadCast branch (identity input)
+    opEye = Eye(m)
+    opNo = BroadCast(opEye, dim_out)
+    opNo2 = copy_operator(opNo; threaded = true)
+    @test opNo2 isa AbstractOperators.NoOperatorBroadCast
+    x = randn(m)
+    y1 = zeros(dim_out)
+    y2 = zeros(dim_out)
+    mul!(y1, opNo, x)
+    mul!(y2, opNo2, x)
+    @test y1 ≈ y2
+
+    # OperatorBroadCast branch (wraps another operator)
+    opA = MatrixOp(randn(m, n))
+    opWrapped = BroadCast(opA, dim_out)
+    opWrapped2 = copy_operator(opWrapped; threaded = true)
+    @test opWrapped2 isa AbstractOperators.OperatorBroadCast
+    x2 = randn(n)
+    y3 = zeros(dim_out)
+    y4 = zeros(dim_out)
+    mul!(y3, opWrapped, x2)
+    mul!(y4, opWrapped2, x2)
+    @test y3 ≈ y4
+end
diff --git a/test/calculus/test_combinations.jl b/test/calculus/test_combinations.jl
index c40cd7e..d8bf354 100644
--- a/test/calculus/test_combinations.jl
+++ b/test/calculus/test_combinations.jl
@@ -1,6 +1,7 @@
 @testitem "Combinations: HCAT and Compose" tags = [:calculus, :Combinations] setup = [TestUtils] begin
     using Random, AbstractOperators
     Random.seed!(42)
+    verb && println(" --- Testing Combinations: HCAT and Compose --- ")
 
     m1, m2, m3, m4 = 4, 7, 3, 2
     A1 = randn(m3, m1)
@@ -44,6 +45,7 @@ end
 @testitem "Combinations: VCAT and HCAT mixtures" tags = [:calculus, :Combinations] setup = [TestUtils] begin
     using Random, AbstractOperators
     Random.seed!(43)
+    verb && println(" --- Testing Combinations: VCAT/HCAT --- ")
 
     # VCAT of HCATs
     m1, m2, n1 = 4, 7, 3
@@ -98,6 +100,7 @@ end
 @testitem "Combinations: Sum structures" tags = [:calculus, :Combinations] setup = [TestUtils] begin
     using Random, AbstractOperators
     Random.seed!(44)
+    verb && println(" --- Testing Combinations: Sum --- ")
 
     # Sum of HCATs
     m, n1, n2, n3 = 4, 7, 5, 3
@@ -142,6 +145,7 @@ end
 @testitem "Combinations: Scale structures" tags = [:calculus, :Combinations] setup = [TestUtils] begin
     using Random, AbstractOperators
     Random.seed!(45)
+    verb && println(" --- Testing Combinations: Scale --- ")
 
     # Scale of DCAT
     m1, n1 = 4, 7
@@ -229,6 +233,7 @@ end
 @testitem "Combinations: Nonlinear" tags = [:calculus, :Combinations] setup = [TestUtils] begin
     using Random, AbstractOperators
     Random.seed!(46)
+    verb && println(" --- Testing Combinations: Nonlinear --- ")
 
     # Nonlinear HCAT of VCAT
     n, m1, m2, m3 = 4, 3, 2, 7
@@ -289,6 +294,208 @@ end
     @test norm(y - (sin.(exp.(x + d1) - d2) .+ d3)) < 1.0e-8
 end
 
+@testitem "Combinations: AffineAdd merging and Zeros" tags = [:calculus, :Combinations] setup = [TestUtils] begin
+    using Random, LinearAlgebra, AbstractOperators
+    Random.seed!(47)
+    verb && println(" --- Testing Combinations: AffineAdd merging and Zeros --- ")
+
+    n = 8
+
+    # AffineAdd(linear) * AffineAdd(linear) — S1==S2==true
+    d1 = randn(n)
+    d2 = randn(n)
+    T1 = AffineAdd(Eye(n), d1)        # x + d1, S=true
+    T2 = AffineAdd(Eye(n), d2)        # x + d2, S=true
+    Tc = T1 * T2                       # should combine, not stay as Compose
+    @test !(Tc isa AbstractOperators.Compose)
+    x = randn(n)
+    @test norm(Tc * x - (x + d2 + d1)) < 1.0e-12
+
+    # AffineAdd(linear) * AffineAdd(linear) — S1=true, S2=false
+    T3 = AffineAdd(Eye(n), d1, true)  # x + d1
+    T4 = AffineAdd(Eye(n), d2, false) # x - d2
+    Tc2 = T3 * T4
+    @test !(Tc2 isa AbstractOperators.Compose)
+    @test norm(Tc2 * x - (x - d2 + d1)) < 1.0e-12
+
+    # AffineAdd(linear) * AffineAdd(linear) — S1=false, S2=true
+    T5 = AffineAdd(Eye(n), d1, false) # x - d1
+    T6 = AffineAdd(Eye(n), d2, true)  # x + d2
+    Tc3 = T5 * T6
+    @test !(Tc3 isa AbstractOperators.Compose)
+    @test norm(Tc3 * x - (x + d2 - d1)) < 1.0e-12
+
+    # combine(linear_op, AffineAdd) with scalar displacement
+    m = 6
+    A = randn(m, n)
+    scalar_d = 2.5
+    opA = MatrixOp(A)
+    opAA_scalar = AffineAdd(Eye(n), scalar_d)  # scalar displacement: x + 2.5
+    Tc4 = opA * opAA_scalar                     # A * (x + 2.5) = A*x + A*fill(2.5,n)
+    @test !(Tc4 isa AbstractOperators.Compose)
+    @test norm(Tc4 * x - A * (x .+ scalar_d)) < 1.0e-12
+
+    # combine(L, Sum) with non-square L
+    n2, m2 = 5, 7
+    A_outer = MatrixOp(randn(n2, m2))  # n2 × m2 (non-square)
+    A1 = MatrixOp(randn(m2, n2))       # m2 × n2
+    A2 = MatrixOp(randn(m2, n2))       # m2 × n2
+    Sop = Sum(A1, A2)
+    combined_sum = A_outer * Sop       # should combine; non-square → Sum(ops...) branch
+    @test !(combined_sum isa AbstractOperators.Compose)
+    x2 = randn(n2)
+    y_ref = A_outer.A * (A1.A * x2 + A2.A * x2)
+    @test norm(combined_sum * x2 - y_ref) < 1.0e-12
+
+    # combine(Zeros, R) — is_null(L): square R with same types (returns L)
+    Z_sq = Zeros(Float64, (n,), Float64, (n,))
+    E_sq = Eye(n)
+    ZE = Z_sq * E_sq
+    @test is_null(ZE)
+    @test size(ZE) == size(Z_sq)
+
+    # combine(Zeros, R) — is_null(L): non-square R (creates new Zeros)
+    # Z_rect: domain=(p,), codomain=(q,); A_rect must have codomain=(p,) to compose with Z_rect
+    p, q = 4, 6
+    Z_rect = Zeros(Float64, (p,), Float64, (q,))  # domain=(p,), codomain=(q,)
+    A_rect = MatrixOp(randn(p, n))                 # domain=(n,), codomain=(p,)
+    ZA = Z_rect * A_rect
+    @test is_null(ZA)
+    @test size(ZA, 1) == (q,)   # codomain of Z_rect
+    @test size(ZA, 2) == (n,)   # domain of A_rect
+
+    # combine(L, Zeros) — is_null(R): square L (returns R)
+    EZ = E_sq * Z_sq
+    @test is_null(EZ)
+    @test size(EZ) == size(Z_sq)
+
+    # combine(L, Zeros) — is_null(R): non-square L (creates new Zeros)
+    A_ns = MatrixOp(randn(m2, n2))  # m2×n2 (non-square)
+    Z_ns = Zeros(Float64, (n2,), Float64, (n2,))  # n2×n2
+    AZ = A_ns * Z_ns
+    @test is_null(AZ)
+    @test size(AZ, 1) == (m2,)
+    @test size(AZ, 2) == (n2,)
+end
+
+@testitem "Combinations: Scale+Compose forwarding branches" tags = [:calculus, :Combinations] setup = [TestUtils] begin
+    using Random, LinearAlgebra, AbstractOperators
+    Random.seed!(48)
+
+    n, m = 4, 5
+    d = randn(n)
+    A = randn(n, m)
+    d2 = randn(n)
+
+    # combine(Scale, Compose): can_be_combined(L.A, R.A[end]) path (line 82)
+    # Compose(DiagOp(d2), MatrixOp(A)) stores as A=(MatrixOp,DiagOp), so A[end] = DiagOp(d2)
+    # can_be_combined(DiagOp(d), DiagOp(d2)) = true → forwarding branch
+    inner_comp = DiagOp(d2) * MatrixOp(A)
+    s_diag = Scale(2.0, DiagOp(d))
+    combined_sc = s_diag * inner_comp
+    x = randn(m)
+    @test combined_sc * x ≈ s_diag * (inner_comp * x)
+
+    # combine(Scale, MatrixOp): can_be_combined(T1.A, T2) = true path (line 199)
+    # Scale(DiagOp) * MatrixOp — can_be_combined(DiagOp, MatrixOp) = true
+    sm = Scale(3.0, DiagOp(d)) * MatrixOp(A)
+    @test sm * x ≈ 3.0 * (d .* (A * x))
+
+    # combine(Scale, AdjointMatrixOp): else branch (line 208)
+    # Scale(FiniteDiff) * MatrixOp(n×n)' — can_be_combined(FiniteDiff, AdjointMatrixOp) = false
+    A2 = randn(n, n)
+    sf = Scale(2.0, FiniteDiff((n,))) * MatrixOp(A2)'
+    xf = randn(n)
+    @test sf * xf ≈ 2.0 * (FiniteDiff((n,)) * (MatrixOp(A2)' * xf))
+
+    # combine(AdjointScale, DiagOp): can_be_combined forwarding branch (line 234)
+    sd = Scale(2.0, DiagOp(d))
+    adj_sd_diag = sd' * DiagOp(d2)
+    xd = randn(n)
+    @test adj_sd_diag * xd ≈ sd' * (DiagOp(d2) * xd)
+
+    # combine(DiagOp, Scale) forwarding branch (line 250)
+    ds = DiagOp(d) * Scale(2.0, DiagOp(d2))
+    @test ds * xd ≈ DiagOp(d) * (2.0 * (DiagOp(d2) * xd))
+
+    # combine(AdjointDiagOp, Scale) forwarding branch (line 258)
+    ads = DiagOp(d)' * Scale(2.0, DiagOp(d2))
+    @test ads * xd ≈ DiagOp(d)' * (2.0 * (DiagOp(d2) * xd))
+end
+
+@testitem "Combinations: AdjointScale and Compose*Scale forwarding branches" tags = [:calculus, :Combinations] setup = [TestUtils] begin
+    using Random, LinearAlgebra, AbstractOperators
+    Random.seed!(49)
+
+    n = 4
+    d, d2 = randn(n), randn(n)
+    A = randn(n, n)
+
+    # combine(AdjointScale, Compose): can_be_combined(Scale.A', Compose.A[end]) = true path (line 90)
+    # Scale(DiagOp(d))' * Compose(DiagOp(d2), MatrixOp(A))
+    # Compose stores as (MatrixOp, DiagOp), A[end] = DiagOp(d2)
+    # can_be_combined(DiagOp(d)', DiagOp(d2)) = true → forwarding
+    sc = Scale(2.0, DiagOp(d))
+    inner_comp = DiagOp(d2) * MatrixOp(A)
+    combined_adj = sc' * inner_comp
+    x = randn(n)
+    @test combined_adj * x ≈ sc' * (inner_comp * x)
+
+    # combine(Compose, Scale): can_be_combined(Compose.A[1], Scale.A) = true path (lines 98-99)
+    # Compose(MatrixOp(A), DiagOp(d)) stores as (DiagOp, MatrixOp), A[1] = DiagOp(d)
+    # can_be_combined(DiagOp(d), DiagOp(d2)) = true → forwarding
+    comp_sc = MatrixOp(A) * DiagOp(d)  # Compose(MatrixOp, DiagOp), stored (DiagOp, MatrixOp)
+    combined_cs = comp_sc * Scale(2.0, DiagOp(d2))
+    @test combined_cs * x ≈ comp_sc * (2.0 * (DiagOp(d2) * x))
+
+    # combine(Compose, AdjointScale): can_be_combined(Compose.A[1], Scale.A') = true path (lines 109-110)
+    # Same Compose(MatrixOp, DiagOp), Scale(DiagOp)' → can_be_combined(DiagOp, AdjointDiagOp) = true
+    combined_cas = comp_sc * Scale(2.0, DiagOp(d2))'
+    @test combined_cas * x ≈ comp_sc * (Scale(2.0, DiagOp(d2))' * x)
+end
+
+
+@testitem "Combinations: Scale+Compose else branches (lines 82, 90, 98-99, 109-110)" tags = [:calculus, :Combinations] setup = [TestUtils] begin
+    using Random, LinearAlgebra, AbstractOperators
+    using AbstractOperators: combine, can_be_combined
+    Random.seed!(0)
+    n = 5
+    # Build a 2-op Compose that doesn't get simplified: FD(n) * MatrixOp(n-1,n-1)
+    # A tuple order is (inner, outer) so comp.A = (FD(n), MatrixOp(n-1,n-1))
+    comp = MatrixOp(randn(n-1, n-1)) * FiniteDiff((n,))   # domain (n,) → codomain (n-1,)
+
+    # Line 82 else: combine(Scale, Compose) when can_be_combined(L.A, comp.A[end]) = false
+    # but can_be_combined(Scale, comp) = true via the all-linear+MatrixOp condition
+    L_82 = Scale(2.0, FiniteDiff((n-1,)))  # domain (n-1,) → codomain (n-2,)
+    @test can_be_combined(L_82, comp)
+    result_82 = combine(L_82, comp)
+    x1 = randn(n)
+    @test result_82 * x1 ≈ L_82 * (comp * x1)
+
+    # Line 90 else: combine(AdjointScale, Compose) when can_be_combined(FD(n)', comp.A[end]) = false
+    # Use MatrixOp(n-1,n-1) * FD(n,) so inner operator (FD') doesn't combine with outer (MatrixOp)
+    comp2 = MatrixOp(randn(n-1, n-1)) * FiniteDiff((n,))  # domain (n,) → codomain (n-1,)
+    adj_L = Scale(2.0, FiniteDiff((n,)))'   # domain (n-1,) → codomain (n,)
+    @test can_be_combined(adj_L, comp2)
+    result_90 = combine(adj_L, comp2)
+    x2 = randn(n)
+    @test result_90 * x2 ≈ adj_L * (comp2 * x2)
+
+    # Lines 98-99 else: combine(Compose, Scale) when can_be_combined(comp.A[1]=FD(n), FD(n+1)) = false
+    scale_inner = Scale(2.0, FiniteDiff((n+1,)))  # domain (n+1,) → codomain (n,)
+    @test can_be_combined(comp, scale_inner)
+    result_98 = combine(comp, scale_inner)
+    x3 = randn(n+1)
+    @test result_98 * x3 ≈ comp * (scale_inner * x3)
+
+    # Lines 109-110 else: combine(Compose, AdjointScale) when can_be_combined(comp.A[1]=FD(n), FD(n)') = false
+    adj_scale_inner = Scale(2.0, FiniteDiff((n,)))'  # domain (n-1,) → codomain (n,)
+    @test can_be_combined(comp, adj_scale_inner)
+    result_109 = combine(comp, adj_scale_inner)
+    x4 = randn(n-1)
+    @test result_109 * x4 ≈ comp * (adj_scale_inner * x4)
+end
+
 @testitem "Combinations (GPU)" tags = [:gpu, :calculus, :Combinations] setup = [TestUtils] begin
     using Random, AbstractOperators, GPUEnv
 
diff --git a/test/calculus/test_compose.jl b/test/calculus/test_compose.jl
index 92b3398..5f27d4f 100644
--- a/test/calculus/test_compose.jl
+++ b/test/calculus/test_compose.jl
@@ -110,6 +110,49 @@ end
     @test_throws DimensionMismatch AbstractOperators.Compose((B, A), ())
 end
 
+@testitem "Compose: type mismatch error" tags = [:calculus, :Compose] setup = [TestUtils] begin
+    using AbstractOperators
+    # Compose L1 ∘ L2 where domain_type(L1) ≠ codomain_type(L2)
+    L1 = MatrixOp(randn(Float64, 3, 4))
+    L2 = MatrixOp(randn(ComplexF64, 4, 5))
+    @test_throws DomainError Compose(L1, L2)
+end
+
+@testitem "Compose: L'*L shortcut to get_normal_op" tags = [:calculus, :Compose] setup = [TestUtils] begin
+    using Random, AbstractOperators
+    Random.seed!(0)
+    # DiagOp has has_optimized_normalop = true
+    # L' * L should short-circuit through line 118-119 in Compose(L1, L2)
+    d = randn(4)
+    op = DiagOp(d)
+    normal = op' * op
+    x = randn(4)
+    @test normal * x ≈ op' * (op * x)
+end
+
+@testitem "get_normal_op(Compose): optimized if branch" tags = [:calculus, :Compose] setup = [TestUtils] begin
+    using Random, AbstractOperators
+    Random.seed!(0)
+    # MatrixOp as outer/left operator has has_optimized_normalop = true
+    # L.A[end] = MatrixOp → hits the if branch in get_normal_op(L::Compose)
+    A = randn(4, 5)
+    L = Compose(MatrixOp(A), FiniteDiff((6,)))
+    N = AbstractOperators.get_normal_op(L)
+    x = randn(6)
+    @test N isa AbstractOperator
+    @test N * x ≈ L' * (L * x)
+end
+
+@testitem "Scale(1, Compose) returns Compose unchanged" tags = [:calculus, :Compose] setup = [TestUtils] begin
+    using Random, AbstractOperators
+    Random.seed!(0)
+    A = randn(4, 5)
+    d = randn(4)
+    comp = Compose(DiagOp(d), MatrixOp(A))
+    s1 = Scale(1.0, comp)
+    @test s1 === comp
+end
+
 @testitem "Adjacent adjoint optimized normal (GetIndex*GetIndex')" tags = [:calculus, :Compose] setup = [TestUtils] begin
     using Random, AbstractOperators
     Random.seed!(0)
@@ -340,6 +383,61 @@ end
     end
 end
 
+@testitem "get_normal_op(Compose) if branch" tags = [:calculus, :Compose] setup = [TestUtils] begin
+    using Random, LinearAlgebra, AbstractOperators
+    Random.seed!(0)
+    # get_normal_op lines 227-229: has_optimized_normalop(L.A[end]) == true
+    # Compose(MatrixOp, FiniteDiff) stores A = (FiniteDiff, MatrixOp); A[end] = MatrixOp
+    L = Compose(MatrixOp(randn(5, 4)), FiniteDiff((5,)))
+    @test AbstractOperators.has_optimized_normalop(L) == true
+    N = AbstractOperators.get_normal_op(L)
+    @test N isa AbstractOperator
+    x = randn(5)
+    @test N * x ≈ L' * (L * x)
+end
+
+@testitem "copy_operator on Compose" tags = [:calculus, :Compose] setup = [TestUtils] begin
+    using Random, AbstractOperators
+    Random.seed!(0)
+    # _copy_operator_impl lines 307-309: use FiniteDiff+MatrixOp (no combine rule → stays Compose)
+    A = randn(4, 4)
+    L = Compose(FiniteDiff((4,)), MatrixOp(A))
+    @test L isa Compose
+    L2 = copy_operator(L)
+    @test L2 isa Compose
+    @test length(L2.A) == length(L.A)
+    x = randn(4)
+    @test L2 * x ≈ L * x
+    y = randn(3)
+    @test L2' * y ≈ L' * y
+end
+
+@testitem "remove_slicing Compose: is_eye path" tags = [:calculus, :Compose] setup = [TestUtils] begin
+    using Random, AbstractOperators
+    Random.seed!(0)
+    # is_sliced(L.A[1]) && is_eye(remove_slicing(L.A[1])) branch (lines 267-273)
+    # AffineAdd(GetIndex): is_sliced=true, not isa GetIndex, remove_slicing returns Eye
+    G = GetIndex((5,), 2:4)
+    b = randn(3)
+    af = AffineAdd(G, b)
+    @test AbstractOperators.is_sliced(af)
+    @test !(af isa GetIndex)
+    @test AbstractOperators.is_eye(AbstractOperators.remove_slicing(af))
+    # length == 2 case: Compose((af, A2), buf) → remove_slicing returns A2
+    A2 = MatrixOp(randn(3, 3))
+    L2 = AbstractOperators.Compose((af, A2), (zeros(3),))
+    out = AbstractOperators.remove_slicing(L2)
+    @test out === A2
+    # length > 2 case: Compose((af, fd, M), bufs) → remove_slicing returns Compose(fd, M)
+    # use FiniteDiff + MatrixOp (no combine rule → stays Compose after remove_slicing)
+    fd = FiniteDiff((3,))
+    A3 = MatrixOp(randn(4, 3))
+    L3 = AbstractOperators.Compose((af, fd, A3), (zeros(3), zeros(3)))
+    out3 = AbstractOperators.remove_slicing(L3)
+    @test out3 isa Compose
+    @test length(out3.A) == 2
+end
+
 @testitem "Compose (GPU)" tags = [:gpu, :calculus, :Compose] setup = [TestUtils] begin
     using Random, AbstractOperators, GPUEnv
 
@@ -356,3 +454,58 @@ end
         test_op(opC2, gpu_randn(backend, n), gpu_randn(backend, n), false)
     end
 end
+
+@testitem "Compose: combine-at-i2 inlines nested Compose (line 73)" tags = [:calculus, :Compose] setup = [TestUtils] begin
+    using Random, AbstractOperators, LinearAlgebra
+    Random.seed!(0)
+    n = 6
+    fd = FiniteDiff((n,))                   # domain (n,) → codomain (n-1,)
+    M = MatrixOp(randn(n - 1, n - 1))       # domain (n-1,) → codomain (n-1,)
+    s = Scale(2.0, FiniteDiff((n - 1,)))    # domain (n-1,) → codomain (n-2,)
+    # At i=2: can_be_combined(s, M) → combine returns a Compose → triggers i -= 1 (line 73)
+    buf1 = zeros(n - 1)
+    buf2 = zeros(n - 1)
+    L = AbstractOperators.Compose((fd, M, s), (buf1, buf2))
+    @test L isa AbstractOperators.Compose
+    x = randn(n)
+    @test L * x ≈ s * (M * (fd * x))
+end
+
+@testitem "Compose: 4-op chain mid-pair combination triggers i-decrement (line 88)" tags = [:calculus, :Compose] setup = [TestUtils] begin
+    using Random, AbstractOperators, LinearAlgebra
+    Random.seed!(0)
+    n = 6
+    d1, d2 = randn(n - 1), randn(n - 1)
+    fd_n = FiniteDiff((n,))        # domain (n,) → codomain (n-1,)
+    diag1 = DiagOp(d1)             # domain/codomain (n-1,)
+    diag2 = DiagOp(d2)             # domain/codomain (n-1,)
+    fd_nm1 = FiniteDiff((n - 1,))  # domain (n-1,) → codomain (n-2,)
+    # At i=2: can_be_combined(diag2, diag1)=true → non-Compose result → elseif branch,
+    # i > 1 and buffers not equal → i -= 1 (line 88)
+    buf1, buf2, buf3 = zeros(n - 1), zeros(n - 1), zeros(n - 1)
+    L = AbstractOperators.Compose((fd_n, diag1, diag2, fd_nm1), (buf1, buf2, buf3))
+    @test L isa AbstractOperators.Compose
+    x = randn(n)
+    @test length(L * x) == n - 2
+    @test L * x ≈ fd_nm1 * (diag2 * (diag1 * (fd_n * x)))
+end
+
+@testitem "Compose: buffer adjacency check after combination (lines 81-82)" tags = [:calculus, :Compose] setup = [TestUtils] begin
+    using Random, AbstractOperators, LinearAlgebra
+    Random.seed!(0)
+    n = 6
+    d1, d2 = randn(n - 1), randn(n - 1)
+    fd_n = FiniteDiff((n,))
+    diag1 = DiagOp(d1)
+    diag2 = DiagOp(d2)
+    fd_nm1 = FiniteDiff((n - 1,))
+    # buf[1] === buf[3]: after DiagOp pair combines at i=2 and buffers are rebuilt,
+    # buf[i-1] === buf[i] is detected → reallocation triggered (lines 81-82)
+    shared_buf = zeros(n - 1)
+    mid_buf = zeros(n - 1)
+    L = AbstractOperators.Compose((fd_n, diag1, diag2, fd_nm1), (shared_buf, mid_buf, shared_buf))
+    @test L isa AbstractOperators.Compose
+    x = randn(n)
+    @test length(L * x) == n - 2
+    @test L * x ≈ fd_nm1 * (diag2 * (diag1 * (fd_n * x)))
+end
diff --git a/test/calculus/test_dcat.jl b/test/calculus/test_dcat.jl
index d6da91d..2b9804f 100644
--- a/test/calculus/test_dcat.jl
+++ b/test/calculus/test_dcat.jl
@@ -1,6 +1,31 @@
+@testitem "DCAT: nested DCAT" tags = [:calculus, :DCAT] setup = [TestUtils] begin
+    using Random, LinearAlgebra, AbstractOperators
+    Random.seed!(0)
+
+    # DCAT does NOT flatten inner DCATs: when an inner DCAT{2} is an element,
+    # _ndoms_from_type returns 2, creating tuple indices in idxD/idxC.
+    # This exercises the `else` (tuple-index) branches in @generated mul!.
+    inner = DCAT(Eye(2), Eye(3))
+    outer = DCAT(inner, Eye(4))
+
+    # The domain storage type expands the inner DCAT's ArrayPartition:
+    # DS = ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}}}
+    # so input is a flat 3-element ArrayPartition.
+    x = ArrayPartition(randn(2), randn(3), randn(4))
+    y_ref = ArrayPartition(randn(2), randn(3), randn(4))
+
+    # test_op verifies forward mul!, adjoint mul!, and adjoint consistency
+    y = test_op(outer, x, y_ref, verb)
+
+    # Since all sub-operators are identity, output == input
+    @test norm(collect(y) - collect(x)) <= 1.0e-12
+    @test norm(collect(outer' * y) - collect(x)) <= 1.0e-12
+end
+
 @testitem "DCAT: basic mul" tags = [:calculus, :DCAT] setup = [TestUtils] begin
     using Random, AbstractOperators
     Random.seed!(0)
+    verb && println(" --- Testing DCAT --- ")
 
     m1, n1, m2, n2, m3, n3 = 4, 7, 5, 2, 5, 5
     A1 = randn(m1, n1)
@@ -91,7 +116,7 @@ end
         Random.seed!(0)
 
         n1, n2 = 3, 4
-        opD = DCAT(DiagOp(gpu_ones(backend, Float64, n1)), DiagOp(2 .* gpu_ones(backend, Float64, n2)))
+        opD = DCAT(DiagOp(gpu_ones(backend, Float64, n1)), DiagOp(to_gpu(backend, 2 .* ones(n2))))
         test_op(
             opD,
             ArrayPartition(gpu_randn(backend, n1), gpu_randn(backend, n2)),
diff --git a/test/calculus/test_hadamardprod.jl b/test/calculus/test_hadamardprod.jl
index 7303971..86b69b2 100644
--- a/test/calculus/test_hadamardprod.jl
+++ b/test/calculus/test_hadamardprod.jl
@@ -1,5 +1,6 @@
 @testitem "HadamardProd: basic mul" tags = [:calculus, :HadamardProd] setup = [TestUtils] begin
     using AbstractOperators
+    verb && println(" --- Testing HadamardProd: basic mul --- ")
 
     # Basic square identity factors (Eye.*Eye)
     n = 3
@@ -37,6 +38,7 @@ end
 
 @testitem "HadamardProd: properties" tags = [:calculus, :HadamardProd] setup = [TestUtils] begin
     using AbstractOperators
+    verb && println(" --- Testing HadamardProd: properties --- ")
 
     # Re-create the HCAT-based P for remove_displacement and permute tests
     m, n = 3, 5
@@ -83,6 +85,7 @@ end
 
 @testitem "HadamardProd: equality and permute" tags = [:calculus, :HadamardProd] setup = [TestUtils] begin
     using AbstractOperators
+    verb && println(" --- Testing HadamardProd: equality and permute --- ")
 
     # Equality / inequality
     n = 3
@@ -132,7 +135,7 @@ end
     @test remove_displacement(Prd) == Prd
 end
 
-@testitem "HadamardProd (GPU)" tags = [:gpu, :calculus, :HadamardProd] setup = [TestUtils] begin
+@testitem "HadamardProd (GPU)" tags = [:gpu, :calculus, :HadamardProd] setup = [TestUtils, GPUNLTestUtils] begin
     using Random, AbstractOperators, GPUEnv
 
     for backend in gpu_backends()
@@ -148,7 +151,8 @@ end
         test_NLop_gpu(P, x, r, false)
 
         n2, l = 3, 2
-        P2 = HadamardProd(Sin(gpu_zeros(backend, Float64, n2, l)), Cos(gpu_zeros(backend, Float64, n2, l)))
+        AT = gpu_wrapper(backend, Float64, n2, l)
+        P2 = HadamardProd(Sin(Float64, (n2, l); array_type = AT), Cos(Float64, (n2, l); array_type = AT))
         x2 = gpu_randn(backend, n2, l)
         r2 = gpu_randn(backend, n2, l)
         test_NLop_gpu(P2, x2, r2, false)
@@ -174,4 +178,11 @@ end
     y3 = zeros(n)
     mul!(y3, P2, x2)
     @test y3 ≈ P * x2
+
+    # Explicit threaded kwarg: exercise _copy_operator_impl unambiguously
+    P3 = copy_operator(P; threaded = true, storage_type = nothing)
+    @test P3 isa HadamardProd
+    y4 = zeros(n)
+    mul!(y4, P3, x)
+    @test y4 ≈ y1
 end
diff --git a/test/calculus/test_hcat.jl b/test/calculus/test_hcat.jl
index c23e6e5..ed503af 100644
--- a/test/calculus/test_hcat.jl
+++ b/test/calculus/test_hcat.jl
@@ -226,7 +226,7 @@ end
         Random.seed!(0)
 
         n = 4
-        opH = HCAT(DiagOp(gpu_ones(backend, Float64, n)), DiagOp(2 .* gpu_ones(backend, Float64, n)))
+        opH = HCAT(DiagOp(gpu_ones(backend, Float64, n)), DiagOp(to_gpu(backend, 2 .* ones(n))))
         test_op(opH, ArrayPartition(gpu_randn(backend, n), gpu_randn(backend, n)), gpu_randn(backend, n), false)
 
         m, n1, n2 = 4, 7, 5
@@ -236,3 +236,54 @@ end
         test_op(opH2, ArrayPartition(gpu_randn(backend, n1), gpu_randn(backend, n2)), gpu_randn(backend, m), false)
     end
 end
+
+@testitem "HCAT fun_name reversed idxs (line 294)" tags = [:calculus, :HCAT] setup = [TestUtils] begin
+    using Random, AbstractOperators
+    Random.seed!(0)
+    n = 8
+    A1 = MatrixOp(randn(4, n))
+    A2 = MatrixOp(randn(4, n))
+    H = HCAT(A1, A2)
+    # permute swaps domain slot ordering → idxs[1] == 2 triggers reversed branch (line 294)
+    Hp = AbstractOperators.permute(H, [2, 1])
+    @test Hp isa HCAT
+    name = AbstractOperators.fun_name(Hp)
+    @test occursin(",", name)
+end
+
+@testitem "HCAT get_slicing_expr: single-expr return (line 337)" tags = [:calculus, :HCAT] setup = [TestUtils] begin
+    using Random, AbstractOperators
+    Random.seed!(0)
+    n = 10
+    # Single-element HCAT → length(exprs) == 1 → return exprs[1]
+    op_gi = GetIndex(Float64, (n,), (1:4,))
+    H_single = HCAT((op_gi,), (1,), zeros(4))
+    @test AbstractOperators.is_sliced(H_single)
+    expr_single = AbstractOperators.get_slicing_expr(H_single)
+    @test expr_single == (1:4,)
+end
+
+@testitem "HCAT get_slicing_expr: multi-element loop (line 330)" tags = [:calculus, :HCAT] setup = [TestUtils] begin
+    using AbstractOperators
+    n = 12
+    op1 = GetIndex(Float64, (n,), (1:4,))
+    op2 = GetIndex(Float64, (n,), (5:8,))
+    op3 = GetIndex(Float64, (n,), (9:12,))
+    H = HCAT(op1, op2, op3)
+    @test AbstractOperators.is_sliced(H)
+    exprs = AbstractOperators.get_slicing_expr(H)
+    @test exprs == ((1:4,), (5:8,), (9:12,))
+end
+
+@testitem "HCAT getindex: tuple-idxs error (line 87)" tags = [:calculus, :HCAT] setup = [TestUtils] begin
+    using Random, AbstractOperators
+    Random.seed!(0)
+    n = 4
+    # HCAT with a Compose(MatrixOp, HCAT) sub-operator gives tuple idxs
+    H_inner = HCAT(MatrixOp(randn(n, n)), MatrixOp(randn(n, n)))
+    C_sub = Compose(MatrixOp(randn(n, n)), H_inner)
+    H_outer = HCAT(MatrixOp(randn(n, n)), C_sub)
+    @test H_outer.idxs == (1, (2, 3))
+    # Selecting partial index into the tuple-idxs sub-op should error
+    @test_throws ErrorException H_outer[2]
+end
diff --git a/test/calculus/test_operatorwrapper.jl b/test/calculus/test_operatorwrapper.jl
index b197ece..d8199ec 100644
--- a/test/calculus/test_operatorwrapper.jl
+++ b/test/calculus/test_operatorwrapper.jl
@@ -59,6 +59,24 @@ end
     test_op(wrapper, randn(n), randn(n), verb)
 end
 
+@testitem "OperatorWrapper: copy_operator" tags = [:calculus, :OperatorWrapper] setup = [TestUtils] begin
+    using Random, LinearAlgebra, AbstractOperators
+    Random.seed!(43)
+
+    n = 8
+    cpu_op = FiniteDiff(Float64, (n,), 1)
+    wrapper = OperatorWrapper(cpu_op)
+    wrapper2 = copy_operator(wrapper; threaded = true)
+    @test wrapper2 isa OperatorWrapper
+
+    x = randn(n)
+    y1 = zeros(n - 1)
+    y2 = zeros(n - 1)
+    mul!(y1, wrapper, x)
+    mul!(y2, wrapper2, x)
+    @test y1 ≈ y2
+end
+
 @testitem "OperatorWrapper (GPU)" tags = [:gpu, :calculus, :OperatorWrapper] setup = [TestUtils] begin
     using Random, AbstractOperators, GPUEnv
 
@@ -66,21 +84,21 @@ end
         Random.seed!(42)
         n = 32
         op = FiniteDiff(Float32, (n,), 1)
-        array_type = gpu_wrapper(backend, Float32, n)
-        wrapper = OperatorWrapper(op; array_type = array_type)
+        storage_type = gpu_wrapper(backend, Float32, n)
+        wrapper = OperatorWrapper(op; array_type = storage_type)
         @test domain_array_type(wrapper) <: backend.array_type
         @test codomain_array_type(wrapper) <: backend.array_type
 
         x = gpu_randn(backend, Float32, n)
         y = gpu_zeros(backend, Float32, n - 1)
         mul!(y, wrapper, x)
-        @test y isa array_type
+        @test y isa storage_type
         @test collect(y) ≈ op * collect(x)
 
         r = gpu_randn(backend, Float32, n - 1)
         z = gpu_zeros(backend, Float32, n)
         mul!(z, wrapper', r)
-        @test z isa array_type
+        @test z isa storage_type
         ref = zeros(Float32, n)
         mul!(ref, op', collect(r))
         @test collect(z) ≈ ref
diff --git a/test/calculus/test_reshape.jl b/test/calculus/test_reshape.jl
index 942458a..9ee5f3c 100644
--- a/test/calculus/test_reshape.jl
+++ b/test/calculus/test_reshape.jl
@@ -1,6 +1,7 @@
 @testitem "Reshape: basic 1D->2D" tags = [:calculus, :Reshape] setup = [TestUtils] begin
     using Random, AbstractOperators
     Random.seed!(0)
+    verb && println(" --- Testing Reshape: basic 1D->2D --- ")
 
     m, n = 8, 4
     dim_out = (2, 2, 2)
@@ -29,6 +30,7 @@ end
 @testitem "Reshape: displacement and storage" tags = [:calculus, :Reshape] setup = [TestUtils] begin
     using Random, AbstractOperators
     Random.seed!(0)
+    verb && println(" --- Testing Reshape: displacement and storage --- ")
 
     # testing displacement
     m, n = 8, 4
@@ -62,6 +64,7 @@ end
 @testitem "Reshape: Scale mul" tags = [:calculus, :Reshape] setup = [TestUtils] begin
     using Random, AbstractOperators
     Random.seed!(0)
+    verb && println(" --- Testing Reshape: Scale mul --- ")
 
     m, n = 8, 4
     coeff = pi
@@ -90,6 +93,7 @@ end
 @testitem "Reshape: Scale properties" tags = [:calculus, :Reshape] setup = [TestUtils] begin
     using Random, AbstractOperators
     Random.seed!(0)
+    verb && println(" --- Testing Reshape: Scale properties --- ")
 
     m, n = 8, 4
     coeff = pi
@@ -147,6 +151,7 @@ end
 @testitem "Reshape: equality and adjoint" tags = [:calculus, :Reshape] setup = [TestUtils] begin
     using Random, AbstractOperators
     Random.seed!(0)
+    verb && println(" --- Testing Reshape: equality and adjoint --- ")
 
     # Equality / inequality
     m, n = 8, 4
@@ -191,6 +196,7 @@ end
 @testitem "Reshape: permute and nonlinear" tags = [:calculus, :Reshape] setup = [TestUtils] begin
     using Random, AbstractOperators
     Random.seed!(0)
+    verb && println(" --- Testing Reshape: permute and nonlinear --- ")
 
     # permute domain ordering (wrap HCAT to get multi-domain) and ensure same behavior when inputs permuted
     mH = 6
diff --git a/test/calculus/test_scale.jl b/test/calculus/test_scale.jl
index 0841910..fe74a39 100644
--- a/test/calculus/test_scale.jl
+++ b/test/calculus/test_scale.jl
@@ -332,3 +332,25 @@ end
         test_op(op, gpu_randn(backend, n), gpu_randn(backend, n), false)
     end
 end
+
+@testitem "Scale(coeff, AdjointMatrixOp) constructor paths" tags = [:calculus, :Scale] setup = [TestUtils] begin
+    using Random, AbstractOperators, LinearAlgebra
+    Random.seed!(0)
+    n = 5
+    A = MatrixOp(randn(n, n))
+    adj_A = A'
+    # coeff == 1 → returns the AdjointOperator unchanged (line 79)
+    result1 = Scale(1.0, adj_A)
+    @test result1 === adj_A
+    x = randn(n)
+    @test result1 * x ≈ adj_A * x
+    # coeff != 1 → returns AdjointOperator(Scale(conj(coeff), A)) (line 84)
+    result2 = Scale(2.0, adj_A)
+    @test result2 * x ≈ 2.0 * (adj_A * x)
+    # complex coefficient on complex adjoint MatrixOp
+    B = MatrixOp(randn(ComplexF64, n, n))
+    adj_B = B'
+    result3 = Scale(1.5 + 0.5im, adj_B)
+    xc = randn(ComplexF64, n)
+    @test result3 * xc ≈ (1.5 + 0.5im) * (adj_B * xc)
+end
diff --git a/test/calculus/test_sum.jl b/test/calculus/test_sum.jl
index 6253d76..9287f24 100644
--- a/test/calculus/test_sum.jl
+++ b/test/calculus/test_sum.jl
@@ -87,8 +87,7 @@ end
         Random.seed!(0)
 
         n = 5
-        x = gpu_ones(backend, Float64, n)
-        opS = Sum(Eye(x), Scale(2.0, Eye(x)))
+        opS = Sum(DiagOp(gpu_ones(backend, Float64, n)), DiagOp(to_gpu(backend, 2 .* ones(n))))
         test_op(opS, gpu_randn(backend, n), gpu_randn(backend, n), false)
 
         m, n2 = 5, 7
@@ -98,3 +97,33 @@ end
         test_op(opS2, gpu_randn(backend, n2), gpu_randn(backend, m), false)
     end
 end
+
+@testitem "Sum: all-Zeros degenerate case (line 79)" tags = [:calculus, :Sum] setup = [TestUtils] begin
+    using AbstractOperators
+    n, m = 5, 4
+    # Sum of two Zeros: @generated code returns A[1] when n_flat == 0 (line 79)
+    z1 = Zeros(Float64, (n,), Float64, (m,))
+    z2 = Zeros(Float64, (n,), Float64, (m,))
+    result = Sum(z1, z2)
+    @test result === z1
+    x = randn(n)
+    @test all(result * x .== 0)
+end
+
+@testitem "Sum: copy_operator" tags = [:calculus, :Sum] setup = [TestUtils] begin
+    using Random, AbstractOperators
+    Random.seed!(4)
+
+    m, n = 5, 7
+    A1 = randn(m, n)
+    A2 = randn(m, n)
+    opS = Sum(MatrixOp(A1), MatrixOp(A2))
+    opS2 = copy_operator(opS; threaded = true)
+    @test opS2 isa Sum
+    x = randn(n)
+    y1 = zeros(m)
+    y2 = zeros(m)
+    mul!(y1, opS, x)
+    mul!(y2, opS2, x)
+    @test y1 ≈ y2
+end
diff --git a/test/calculus/test_vcat.jl b/test/calculus/test_vcat.jl
index b3cac06..c459077 100644
--- a/test/calculus/test_vcat.jl
+++ b/test/calculus/test_vcat.jl
@@ -111,6 +111,28 @@ end
     @test exprs[1] == (1:5,) && exprs[2] == (6:10,)
     @test !is_sliced(AbstractOperators.remove_slicing(Vs))
 
+    # remove_slicing second branch: VCAT of HCATs with GetIndex of different output sizes
+    # Forces the elseif branch where new_ops have different domain sizes after remove_slicing
+    g3 = GetIndex(Float64, (5,), (1:3,))  # 5-dim → 3-dim
+    g4 = GetIndex(Float64, (5,), (1:4,))  # 5-dim → 4-dim
+    z3 = Zeros(Float64, (5,), Float64, (3,))  # domain (5,), codomain (3,)
+    z4 = Zeros(Float64, (5,), Float64, (4,))  # domain (5,), codomain (4,)
+    H1 = HCAT(g3, z3)  # 3×(5+5), both ops have codomain (3,)
+    H2 = HCAT(z4, g4)  # 4×(5+5), both ops have codomain (4,)
+    Vs2 = VCAT(H1, H2)  # 7×10
+    @test is_sliced(Vs2)
+    rs2 = AbstractOperators.remove_slicing(Vs2)
+    @test rs2 isa VCAT
+    @test !is_sliced(rs2)
+    # Result should act like DCAT(Eye(3), Eye(4)): maps (x1, x2) → (x1, x2) with no-op zeros
+    b3 = randn(3)
+    b4 = randn(4)
+    b_domain = ArrayPartition(b3, b4)
+    y_codomain = ArrayPartition(randn(3), randn(4))
+    mul!(y_codomain, rs2, b_domain)
+    @test y_codomain.x[1] ≈ b3
+    @test y_codomain.x[2] ≈ b4
+
     # fun_name and equality
     A1 = Eye(3)
     A2 = Eye(3)
@@ -134,7 +156,7 @@ end
         Random.seed!(0)
 
         n = 4
-        opV = VCAT(DiagOp(gpu_ones(backend, Float64, n)), DiagOp(2 .* gpu_ones(backend, Float64, n)))
+        opV = VCAT(DiagOp(gpu_ones(backend, Float64, n)), DiagOp(to_gpu(backend, 2 .* ones(n))))
         test_op(opV, gpu_randn(backend, n), ArrayPartition(gpu_randn(backend, n), gpu_randn(backend, n)), false)
 
         m1, m2, n = 4, 7, 5
@@ -145,6 +167,31 @@ end
     end
 end
 
+@testitem "VCAT: constructor errors" tags = [:calculus, :VCAT] setup = [TestUtils] begin
+    using Random, AbstractOperators
+    Random.seed!(0)
+
+    # Line 49: domain dimension mismatch → DimensionMismatch
+    @test_throws DimensionMismatch VCAT(MatrixOp(randn(3, 4)), MatrixOp(randn(3, 5)))
+
+    # Line 52: domain type mismatch → generic error (throw(error(...)))
+    @test_throws Exception VCAT(MatrixOp(randn(3, 4)), MatrixOp(ones(ComplexF64, 2, 4)))
+end
+
+@testitem "VCAT remove_slicing: unsupported VCAT error (line 200)" tags = [:calculus, :VCAT] setup = [TestUtils] begin
+    using AbstractOperators
+    n = 8
+    # VCAT of HCATs with no null operators: both removal branches fail → error at line 200
+    g1 = GetIndex(Float64, (n,), (1:4,))
+    g2 = GetIndex(Float64, (n,), (3:6,))
+    g3 = GetIndex(Float64, (n,), (5:8,))
+    H1 = HCAT(g1, g2)
+    H2 = HCAT(g2, g3)
+    Vs = VCAT(H1, H2)
+    @test AbstractOperators.is_sliced(Vs)
+    @test_throws ErrorException AbstractOperators.remove_slicing(Vs)
+end
+
 @testitem "VCAT: copy_operator" tags = [:calculus, :VCAT] setup = [TestUtils] begin
     using Random, AbstractOperators
     Random.seed!(4)
diff --git a/test/dsp/test_dsp_operators.jl b/test/dsp/test_dsp_operators.jl
index 65013a3..9430811 100644
--- a/test/dsp/test_dsp_operators.jl
+++ b/test/dsp/test_dsp_operators.jl
@@ -31,6 +31,31 @@
     @test is_full_column_rank(op) == true
 end
 
+@testitem "Conv complex domain" tags = [:dsp, :Conv] setup = [TestUtils] begin
+    using DSPOperators, DSP, LinearAlgebra, Random
+    Random.seed!(0)
+    n, m = 5, 6
+    h = randn(ComplexF64, m)
+    op = Conv(ComplexF64, (n,), h)
+    x1 = randn(ComplexF64, n)
+    y1 = test_op(op, x1, randn(ComplexF64, n + m - 1), verb)
+    y2 = conv(x1, h)
+    @test all(norm.(y1 .- y2) .<= 1.0e-10)
+end
+
+@testitem "Filt: a[1] != 1 normalization" tags = [:dsp, :Filt] setup = [TestUtils] begin
+    using DSPOperators, DSP, LinearAlgebra, Random
+    Random.seed!(0)
+    n = 10
+    b = [2.0; 0.0; 2.0; 0.0; 0.0]
+    a = [2.0; 2.0; 2.0]    # a[1] != 1, triggers normalization
+    op = Filt(Float64, (n,), b, a)
+    # after normalization b/2 and a/2 => same IIR
+    op_ref = Filt(Float64, (n,), b ./ 2, a ./ 2)
+    x1 = randn(n)
+    @test op * x1 ≈ op_ref * x1
+end
+
 @testitem "Filt: IIR and FIR mappings" tags = [:dsp, :Filt] setup = [TestUtils] begin
     using DSPOperators, DSP, LinearAlgebra, Random
 
@@ -210,7 +235,7 @@ end
 
 @testitem "Conv (GPU)" tags = [:gpu, :dsp, :Conv] setup = [TestUtils] begin
     using DSPOperators, DSP, GPUEnv, LinearAlgebra, Random
-    for backend in gpu_backends(supports_fftw = true)
+    for backend in gpu_backends(; include_jlarrays = false, supports_fftw = true)
         Random.seed!(0)
         n, m = 20, 6
         h_cpu = randn(m)
@@ -239,7 +264,7 @@ end
 
 @testitem "Xcorr (GPU)" tags = [:gpu, :dsp, :Xcorr] setup = [TestUtils] begin
     using DSPOperators, DSP, GPUEnv, LinearAlgebra, Random
-    for backend in gpu_backends(supports_fftw = true)
+    for backend in gpu_backends(; include_jlarrays = false, supports_fftw = true)
         Random.seed!(0)
         n, m = 15, 5
         h_cpu = randn(m)
@@ -269,7 +294,7 @@ end
 
 @testitem "Filt (GPU, FIR)" tags = [:gpu, :dsp, :Filt] setup = [TestUtils] begin
     using DSPOperators, DSP, GPUEnv, LinearAlgebra, Random
-    for backend in gpu_backends(supports_fftw = true)
+    for backend in gpu_backends(; include_jlarrays = false, supports_fftw = true)
         Random.seed!(42)
         n = 20
         b = randn(5)
@@ -297,7 +322,7 @@ end
 
 @testitem "MIMOFilt (GPU, FIR)" tags = [:gpu, :dsp, :MIMOFilt] setup = [TestUtils] begin
     using DSPOperators, DSP, GPUEnv, LinearAlgebra, Random
-    for backend in gpu_backends(supports_fftw = true)
+    for backend in gpu_backends(; include_jlarrays = false, supports_fftw = true)
         Random.seed!(7)
         m, n = 10, 3
         b = [randn(5), randn(3), randn(4), randn(5), randn(3), randn(4)]
@@ -323,3 +348,17 @@ end
         @test collect(z) ≈ z_cpu atol = 1.0e-10
     end
 end
+
+@testitem "Xcorr complex domain" tags = [:dsp, :Xcorr] setup = [TestUtils] begin
+    using DSPOperators, LinearAlgebra, Random
+    Random.seed!(0)
+    n, m = 5, 6
+    h = randn(ComplexF64, m)
+    op = Xcorr(ComplexF64, (n,), h)
+    x1 = randn(ComplexF64, n)
+    y1 = op * x1
+    @test length(y1) == 2 * max(n, m) - 1
+    z1 = op' * y1
+    @test length(z1) == n
+    @test z1 ≈ op' * (op * x1)
+end
diff --git a/test/fftw/test_fftw_operators.jl b/test/fftw/test_fftw_operators.jl
index dec66fe..e95059b 100644
--- a/test/fftw/test_fftw_operators.jl
+++ b/test/fftw/test_fftw_operators.jl
@@ -177,6 +177,20 @@ end
     @test norm(op * (op' * y1) - diag_AAc(op) * y1) <= 1.0e-12
 end
 
+@testitem "DFT ORTHO normalization" tags = [:fftw, :DFT] setup = [TestUtils] begin
+    using AbstractOperators, FFTW, LinearAlgebra, Random, FFTWOperators
+    Random.seed!(0)
+    n = 8
+    op = DFT(Float64, (n,); normalization = FFTWOperators.ORTHO)
+    x1 = randn(n)
+    y1 = op * x1
+    y2 = fft(x1) ./ sqrt(n)
+    @test norm(y1 .- y2) <= 1.0e-12
+    # Adjoint applies ORTHO scaling too
+    z1 = op' * y1
+    @test norm(z1 .- x1) <= 1.0e-12
+end
+
 @testitem "IDFT" tags = [:fftw, :IDFT] setup = [TestUtils] begin
     using AbstractOperators
     using FFTW, LinearAlgebra, Random, FFTWOperators
@@ -380,7 +394,7 @@ end
 @testitem "DFT/RDFT/IRDFT (GPU)" tags = [:gpu, :fftw, :DFT, :RDFT, :IRDFT] setup = [TestUtils] begin
     using FFTW, FFTWOperators, GPUEnv, LinearAlgebra, Random, AbstractOperators
 
-    for backend in gpu_backends(supports_fftw = true)
+    for backend in gpu_backends(; include_jlarrays = false, supports_fftw = true)
         Random.seed!(0)
         n, m = 8, 6
 
diff --git a/test/fftw/test_shift_operators.jl b/test/fftw/test_shift_operators.jl
index 4903555..fad8e3d 100644
--- a/test/fftw/test_shift_operators.jl
+++ b/test/fftw/test_shift_operators.jl
@@ -90,10 +90,34 @@ end
 
     y = similar(z)
     @test_throws ArgumentError alternate_sign!(y, z)
+
+    # too many dirs (N < M)
+    m2d = ones(2, 2)
+    @test_throws ArgumentError alternate_sign!(m2d, (1, 2, 3))
+    y2d = similar(m2d)
+    @test_throws ArgumentError alternate_sign!(y2d, m2d, (1, 2, 3))
+
+    # unsorted dirs
+    @test_throws ArgumentError alternate_sign!(m2d, (2, 1))
+    @test_throws ArgumentError alternate_sign!(y2d, m2d, (2, 1))
+
+    # out-of-range dirs
+    @test_throws ArgumentError alternate_sign!(z, 0)
+    @test_throws ArgumentError alternate_sign!(z, 2)  # 1D array, dir=2 > ndims
+
+    # non-threaded paths (threaded=false)
+    v2 = collect(1.0:4.0)
+    alternate_sign!(v2, 1; threaded = false)
+    @test v2 == [1.0, -2.0, 3.0, -4.0]
+
+    x2 = collect(reshape(1.0:4.0, 2, 2))
+    y2 = similar(x2)
+    alternate_sign!(y2, x2, 1, 2; threaded = false)
+    @test y2 == [1.0 -3.0; -2.0 4.0]
 end
 
 @testitem "fftshift/ifftshift wrappers" tags = [:fftw, :FFTShift] setup = [TestUtils] begin
-    using FFTW, LinearAlgebra, Random, FFTWOperators
+    using FFTW, LinearAlgebra, Random, FFTWOperators, AbstractOperators
     # Even length
     n = 4
     A = DFT(n)
@@ -132,6 +156,27 @@ end
 
     T6 = ifftshift_op(A; domain_shifts = (1,))
     @test (T6 * x) ≈ (A * FFTW.ifftshift(x, (1,)))
+
+    # Compose operators: DiagOp * DFT (all-diagonal/DFT → _is_dft_op true from all() branch)
+    n = 8
+    Random.seed!(42)
+    dft_c = DFT(ComplexF64, n)
+    d = randn(ComplexF64, n)
+    diag_op = DiagOp(d)
+    composed1 = diag_op * dft_c  # Compose: DFT applied first, then DiagOp
+    xc = randn(ComplexF64, n)
+    T7 = fftshift_op(composed1; domain_shifts = (1,))
+    @test T7 * xc ≈ composed1 * FFTW.fftshift(xc, (1,))
+    T8 = fftshift_op(composed1; codomain_shifts = (1,))
+    @test T8 * xc ≈ FFTW.fftshift(composed1 * xc, (1,))
+
+    # MatrixOp * DFT (else branch: not all-diagonal/DFT, but first subop is DFT)
+    mat_op = MatrixOp(randn(ComplexF64, n, n))
+    composed2 = mat_op * dft_c  # Compose: DFT applied first, then MatrixOp
+    T9 = fftshift_op(composed2; domain_shifts = (1,))
+    @test T9 * xc ≈ composed2 * FFTW.fftshift(xc, (1,))
+    T10 = fftshift_op(composed2; codomain_shifts = (1,))
+    @test T10 * xc ≈ FFTW.fftshift(composed2 * xc, (1,))
 end
 
 @testitem "Combination rules: FFTShift/IFFTShift with DFT/IDFT" tags = [:fftw, :CombinationRules] setup = [TestUtils] begin
diff --git a/test/gpu_nl_test_utils.jl b/test/gpu_nl_test_utils.jl
new file mode 100644
index 0000000..4773105
--- /dev/null
+++ b/test/gpu_nl_test_utils.jl
@@ -0,0 +1,34 @@
+@testsnippet GPUNLTestUtils begin
+    using Test
+    using LinearAlgebra
+    using RecursiveArrayTools
+    using AbstractOperators
+
+    _to_cpu(x::AbstractArray) = collect(x)
+    _to_cpu(x::RecursiveArrayTools.ArrayPartition) =
+        RecursiveArrayTools.ArrayPartition(collect.(x.x)...)
+
+    function _assert_cpu_approx(x, y; atol = 1.0e-8)
+        @test norm(_to_cpu(x) .- _to_cpu(y)) <= atol
+    end
+
+    function test_NLop_gpu(A::AbstractOperator, x, y, verb::Bool = false)
+        verb && (println(), println(A))
+
+        Ax = A * x
+        Ax2 = similar(Ax)
+        mul!(Ax2, A, x)
+        _assert_cpu_approx(Ax, Ax2)
+
+        @test_throws ErrorException A'
+
+        J = Jacobian(A, x)
+        grad = J' * y
+        mul!(Ax2, A, x)
+        grad2 = similar(grad)
+        mul!(grad2, J', y)
+        _assert_cpu_approx(grad, grad2; atol = 1.0e-8)
+
+        return Ax, grad
+    end
+end
diff --git a/test/linearoperators/test_diagop.jl b/test/linearoperators/test_diagop.jl
index e50eb88..41a0c60 100644
--- a/test/linearoperators/test_diagop.jl
+++ b/test/linearoperators/test_diagop.jl
@@ -103,12 +103,36 @@ end
     @test size(op) == ((n,), (n,))
 end
 
+@testitem "DiagOp: copy_operator" tags = [:linearoperator, :DiagOp] setup = [TestUtils] begin
+    using Random, AbstractOperators
+    Random.seed!(5)
+
+    n = 6
+    d = randn(n)
+    op = DiagOp(d)
+    op2 = copy_operator(op; threaded = true)
+    @test op2 isa DiagOp
+    # DiagOp's diagonal `d` is treated as read-only data and shared (not deep-copied)
+    # when storage_type is not given, so op2 may be egal to op for these immutable structs.
+    x = randn(n)
+    y1 = zeros(n)
+    y2 = zeros(n)
+    mul!(y1, op, x)
+    mul!(y2, op2, x)
+    @test y1 ≈ y2
+
+    op3 = copy_operator(op; storage_type = Array)
+    @test op3 isa DiagOp
+    @test op3.d !== op.d
+end
+
 @testitem "DiagOp (GPU)" tags = [:gpu, :linearoperator, :DiagOp] setup = [TestUtils, DiagOpTestHelper] begin
     using Random, AbstractOperators, GPUEnv
 
     for backend in gpu_backends()
         Random.seed!(0)
-        test_diagop_mul(x -> to_gpu(backend, x), false, test_op, to_cpu, norm)
+        conv = x -> to_gpu(backend, x)
+        test_diagop_mul(conv, false, test_op, collect, norm)
         x = gpu_randn(backend, 4)
         op = DiagOp(x)
         @test domain_array_type(op) <: backend.array_type
diff --git a/test/linearoperators/test_eye.jl b/test/linearoperators/test_eye.jl
index d37a3bb..3030565 100644
--- a/test/linearoperators/test_eye.jl
+++ b/test/linearoperators/test_eye.jl
@@ -18,11 +18,13 @@
 end  # @testmodule EyeTestHelper
 
 @testitem "Eye" tags = [:linearoperator, :Eye] setup = [TestUtils, EyeTestHelper] begin
-    using Random, AbstractOperators
+    using Random, AbstractOperators, JLArrays
     Random.seed!(0)
 
     test_eye_mul(identity, verb, test_op, to_cpu, norm)
 
+    test_eye_mul(identity, verb, test_op, to_cpu, norm)
+
     n = 4
     op = Eye(Float64, (n,))
     x1 = randn(n)
@@ -94,6 +96,7 @@ end
 
     for backend in gpu_backends()
         Random.seed!(0)
-        test_eye_mul(x -> to_gpu(backend, x), false, test_op, to_cpu, norm)
+        conv = x -> to_gpu(backend, x)
+        test_eye_mul(conv, false, test_op, collect, norm)
     end
 end
diff --git a/test/linearoperators/test_finitediff.jl b/test/linearoperators/test_finitediff.jl
index cbdb97b..4d0f6bd 100644
--- a/test/linearoperators/test_finitediff.jl
+++ b/test/linearoperators/test_finitediff.jl
@@ -25,6 +25,8 @@ end  # @testmodule FiniteDiffTestHelper
 
     test_finitediff_mul(identity, verb, test_op)
 
+    test_finitediff_mul(identity, verb, test_op)
+
     n = 10
     op = FiniteDiff(Float64, (n,))
     x1 = randn(n)
@@ -86,11 +88,21 @@ end
     @test size(FiniteDiff(Float64, (3, 4, 5), 2)) == ((3, 3, 5), (3, 4, 5))
 end
 
-@testitem "FiniteDiff (GPU)" tags = [:gpu, :linearoperator, :FiniteDiff] setup = [TestUtils, FiniteDiffTestHelper] begin
+@testitem "FiniteDiff (GPU)" tags = [:gpu, :linearoperator, :FiniteDiff] setup = [TestUtils] begin
     using Random, AbstractOperators, GPUEnv
 
     for backend in gpu_backends()
         Random.seed!(0)
-        test_finitediff_mul(x -> to_gpu(backend, x), false, test_op)
+
+        n = 10
+        op = FiniteDiff(Float64, (n,); array_type = gpu_wrapper(backend, Float64, n))
+        test_op(op, gpu_randn(backend, n), gpu_randn(backend, n - 1), false)
+
+        n, m = 10, 5
+        op = FiniteDiff(Float64, (n, m); array_type = gpu_wrapper(backend, Float64, n, m))
+        test_op(op, gpu_randn(backend, n, m), gpu_randn(backend, n - 1, m), false)
+
+        op2 = FiniteDiff(Float64, (n, m), 2; array_type = gpu_wrapper(backend, Float64, n, m))
+        test_op(op2, gpu_randn(backend, n, m), gpu_randn(backend, n, m - 1), false)
     end
 end
diff --git a/test/linearoperators/test_getindex.jl b/test/linearoperators/test_getindex.jl
index 4645acf..89bc1ef 100644
--- a/test/linearoperators/test_getindex.jl
+++ b/test/linearoperators/test_getindex.jl
@@ -107,7 +107,7 @@ end
     @test occursin("↓", String(take!(io)))
 end
 
-@testitem "GetIndex (GPU)" tags = [:gpu, :linearoperator, :GetIndex] setup = [TestUtils] begin
+@testitem "GetIndex (GPU)" tags = [:linearoperator, :GetIndex, :gpu] setup = [TestUtils] begin
     using Random, AbstractOperators, GPUEnv
 
     for backend in gpu_backends()
diff --git a/test/linearoperators/test_lbfgs.jl b/test/linearoperators/test_lbfgs.jl
index 0594e73..b281be1 100644
--- a/test/linearoperators/test_lbfgs.jl
+++ b/test/linearoperators/test_lbfgs.jl
@@ -1,6 +1,7 @@
 @testitem "L-BFGS: construction and basic mul" tags = [:linearoperator, :LBFGS] setup = [TestUtils] begin
     using AbstractOperators
     using AbstractOperators: LBFGS, update!, mul!, reset!
+    verb && println(" --- Testing L-BFGS: construction and basic mul --- ")
 
     mem = 3
     x = zeros(10)
@@ -25,6 +26,7 @@ end
 @testitem "L-BFGS: update and two-loop recursion" tags = [:linearoperator, :LBFGS] setup = [TestUtils] begin
     using AbstractOperators
     using AbstractOperators: LBFGS, update!, mul!, reset!
+    verb && println(" --- Testing L-BFGS: update and two-loop recursion --- ")
 
     Q = [
         32.0 13.1 -4.9 -3.0 6.0 2.2 2.6 3.4 -1.9 -7.5
@@ -130,6 +132,7 @@ end
 @testitem "L-BFGS: memory limit and reset" tags = [:linearoperator, :LBFGS] setup = [TestUtils] begin
     using AbstractOperators
     using AbstractOperators: LBFGS, update!, mul!, reset!
+    verb && println(" --- Testing L-BFGS: memory limit and reset --- ")
 
     Q = [
         32.0 13.1 -4.9 -3.0 6.0 2.2 2.6 3.4 -1.9 -7.5
diff --git a/test/linearoperators/test_lmatrixop.jl b/test/linearoperators/test_lmatrixop.jl
index a3c7499..392da54 100644
--- a/test/linearoperators/test_lmatrixop.jl
+++ b/test/linearoperators/test_lmatrixop.jl
@@ -2,11 +2,16 @@
     using Random, AbstractOperators
 
     Random.seed!(0)
+    verb && println(" --- Testing LMatrixOp: basic mul --- ")
 
-    n, m = 5, 6
-    b = randn(m)
-    op = LMatrixOp(Float64, (n, m), b)
-    test_op(op, randn(n, m), randn(n), verb)
+    function test_lmatrixop_mul(conv, verb)
+        n, m = 5, 6
+        b = randn(m)
+        op = LMatrixOp(Float64, (n, m), conv(b))
+        test_op(op, conv(randn(n, m)), conv(randn(n)), verb)
+    end
+
+    test_lmatrixop_mul(identity, verb)
 
     n, m = 5, 6
     b = randn(m)
@@ -54,6 +59,7 @@ end
 @testitem "LMatrixOp: other constructors" tags = [:linearoperator, :LMatrixOp] setup = [TestUtils] begin
     using Random, AbstractOperators
     Random.seed!(0)
+    verb && println(" --- Testing LMatrixOp: other constructors --- ")
 
     n, m, l = 5, 6, 7
 
@@ -95,6 +101,7 @@ end
 @testitem "LMatrixOp: scale and properties" tags = [:linearoperator, :LMatrixOp] setup = [TestUtils] begin
     using Random, AbstractOperators
     Random.seed!(0)
+    verb && println(" --- Testing LMatrixOp: scale and properties --- ")
 
     n, m, l = 5, 6, 7
     bvec = randn(m)
@@ -135,7 +142,7 @@ end
     @test occursin("(⋅)b", s)
 end
 
-@testitem "LMatrixOp (GPU)" tags = [:gpu, :linearoperator, :LMatrixOp] setup = [TestUtils] begin
+@testitem "LMatrixOp (GPU)" tags = [:linearoperator, :LMatrixOp, :gpu] setup = [TestUtils] begin
     using Random, AbstractOperators, GPUEnv
     for backend in gpu_backends()
         Random.seed!(0)
diff --git a/test/linearoperators/test_matrixop.jl b/test/linearoperators/test_matrixop.jl
index 4a9e3ba..86ace6f 100644
--- a/test/linearoperators/test_matrixop.jl
+++ b/test/linearoperators/test_matrixop.jl
@@ -114,11 +114,11 @@ end
     for backend in gpu_backends()
         Random.seed!(0)
         n, m = 5, 4
-        A = gpu_randn(backend, n, m)
-        test_op(MatrixOp(A), gpu_randn(backend, m), gpu_randn(backend, n), false)
-        Ac = gpu_randn(backend, ComplexF64, n, m)
+        A = randn(n, m)
+        test_op(MatrixOp(to_gpu(backend, A)), gpu_randn(backend, m), gpu_randn(backend, n), false)
+        Ac = randn(n, m) + im * randn(n, m)
         test_op(
-            MatrixOp(Ac),
+            MatrixOp(to_gpu(backend, Ac)),
             gpu_randn(backend, ComplexF64, m),
             gpu_randn(backend, ComplexF64, n),
             false,
diff --git a/test/linearoperators/test_mylinop.jl b/test/linearoperators/test_mylinop.jl
index 9555932..bfb08fa 100644
--- a/test/linearoperators/test_mylinop.jl
+++ b/test/linearoperators/test_mylinop.jl
@@ -60,6 +60,31 @@ end
     @test op2 * x1 ≈ A * x1
 end
 
+@testitem "AbstractOperator fallback properties" tags = [:misc, :MyLinOp] setup = [TestUtils] begin
+    using AbstractOperators, LinearAlgebra, Random
+    Random.seed!(0)
+    n, m = 5, 4
+    A = randn(n, m)
+    op = MyLinOp(Float64, (m,), (n,), (y, x) -> mul!(y, A, x), (y, x) -> mul!(y, A', x))
+
+    # is_thread_safe falls back to false for custom operators
+    @test is_thread_safe(op) == false
+    # is_sliced falls back to false
+    @test is_sliced(op) == false
+    # get_slicing_expr falls back to Colon() for non-null operators
+    @test AbstractOperators.get_slicing_expr(op) == Colon()
+    # get_slicing_mask throws for operators without specialization
+    @test_throws ErrorException AbstractOperators.get_slicing_mask(op)
+    # has_optimized_normalop falls back to false
+    @test AbstractOperators.has_optimized_normalop(op) == false
+    # get_normal_op falls back to L' * L
+    normal = AbstractOperators.get_normal_op(op)
+    x = randn(m)
+    @test normal * x ≈ A' * (A * x)
+    # diag throws for non-diagonal operators
+    @test_throws ErrorException diag(op)
+end
+
 @testitem "MyLinOp (GPU)" tags = [:gpu, :linearoperator, :MyLinOp] setup = [TestUtils] begin
     using Random, AbstractOperators, GPUEnv
 
diff --git a/test/linearoperators/test_variation.jl b/test/linearoperators/test_variation.jl
index d436bcf..f52eb90 100644
--- a/test/linearoperators/test_variation.jl
+++ b/test/linearoperators/test_variation.jl
@@ -1,6 +1,15 @@
 @testitem "Variation: basic mul" tags = [:linearoperator, :Variation] setup = [TestUtils] begin
     using Random, SparseArrays, LinearAlgebra, AbstractOperators
     Random.seed!(0)
+    verb && println(" --- Testing Variation: basic mul --- ")
+
+    function test_variation_mul(conv, verb)
+        n, m = 10, 5
+        op = Variation(conv(zeros(Float64, n, m)); threaded = false)
+        test_op(op, conv(randn(n, m)), conv(randn(n * m, 2)), verb)
+    end
+
+    test_variation_mul(identity, verb)
 
     for threaded in (false, true)
         n, m = 10, 5
@@ -42,6 +51,7 @@ end
 @testitem "Variation: 3D mul and constructors" tags = [:linearoperator, :Variation] setup = [TestUtils] begin
     using Random, AbstractOperators
     Random.seed!(0)
+    verb && println(" --- Testing Variation: 3D mul and constructors --- ")
 
     for threaded in (false, true)
         n, m, l = 100, 50, 30
@@ -74,6 +84,7 @@ end
 @testitem "Variation: adjoint and properties" tags = [:linearoperator, :Variation] setup = [TestUtils] begin
     using Random, LinearAlgebra, AbstractOperators
     Random.seed!(0)
+    verb && println(" --- Testing Variation: adjoint and properties --- ")
 
     for threaded in (false, true)
         n, m, l = 100, 50, 30
@@ -123,6 +134,22 @@ end
     end
 end
 
+@testitem "Variation: copy_operator" tags = [:linearoperator, :Variation] setup = [TestUtils] begin
+    using Random, AbstractOperators
+    Random.seed!(6)
+
+    n, m = 10, 5
+    op = Variation(zeros(Float64, n, m); threaded = false)
+    op2 = copy_operator(op; threaded = true)
+    @test op2 isa Variation
+    x = randn(n, m)
+    y1 = zeros(n * m, 2)
+    y2 = zeros(n * m, 2)
+    mul!(y1, op, x)
+    mul!(y2, op2, x)
+    @test y1 ≈ y2
+end
+
 @testitem "Variation (GPU)" tags = [:gpu, :linearoperator, :Variation] setup = [TestUtils] begin
     using Random, AbstractOperators, GPUEnv
 
diff --git a/test/linearoperators/test_zeropad.jl b/test/linearoperators/test_zeropad.jl
index b6f9d4f..f343020 100644
--- a/test/linearoperators/test_zeropad.jl
+++ b/test/linearoperators/test_zeropad.jl
@@ -103,6 +103,7 @@ end
 
     for backend in gpu_backends()
         Random.seed!(0)
-        test_zeropad_mul(x -> to_gpu(backend, x), false, test_op, to_cpu, norm)
+        conv = x -> to_gpu(backend, x)
+        test_zeropad_mul(conv, false, test_op, collect, norm)
     end
 end
diff --git a/test/linearoperators/test_zerosop.jl b/test/linearoperators/test_zerosop.jl
index fd49ad1..2b64b38 100644
--- a/test/linearoperators/test_zerosop.jl
+++ b/test/linearoperators/test_zerosop.jl
@@ -8,6 +8,17 @@
     y1 = test_op(op, randn(n), randn(m), verb)
     @test y1 == zeros(Float64, m)
 
+    function test_zeros_mul(conv, verb)
+        n = (3, 4)
+        m = (5, 2)
+        ST = Base.typename(typeof(conv(zeros(Float64, 1)))).wrapper
+        op = Zeros(Float64, n, Float64, m; array_type = ST)
+        y1 = test_op(op, conv(randn(n)), conv(randn(m)), verb)
+        @test to_cpu(y1) == zeros(Float64, m)
+    end
+
+    test_zeros_mul(identity, verb)
+
     n = (3, 4)
     D = Float64
     m = (5, 2)
diff --git a/test/test_LinearMapsExt.jl b/test/test_LinearMapsExt.jl
index a957911..436de9e 100644
--- a/test/test_LinearMapsExt.jl
+++ b/test/test_LinearMapsExt.jl
@@ -25,6 +25,7 @@
     @test is_full_column_rank(LM)
     @test is_positive_definite(LM) == all(d .> 0)
     @test is_positive_semidefinite(LM) == all(d .>= 0)
+    @test isposdef(LM) == all(d .> 0)
 
     # Complex Diagonal
     d = rand(ComplexF64, 10, 10)
diff --git a/test/test_combination_rules.jl b/test/test_combination_rules.jl
index b882756..8af736e 100644
--- a/test/test_combination_rules.jl
+++ b/test/test_combination_rules.jl
@@ -727,3 +727,22 @@ end
     c7 = combine(scl_adj, mat_adj)
     @test c7 * x ≈ scl_adj * (mat_adj * x)
 end
+
+@testitem "CR: generic combine branches" tags = [:calculus, :CombinationRules] begin
+    using AbstractOperators
+    using AbstractOperators: combine
+
+    n = 4
+    mat = MatrixOp(randn(n, n))
+    eye = Eye(Float64, (n,))
+    z = Zeros(Float64, (n,), Float64, (n,))
+
+    # combine(L, R) where is_eye(R) → returns L (line 345)
+    result = combine(mat, eye)
+    @test result === mat
+
+    # combine(L, R) where is_null(R), L linear, zero displacement, square → returns R (line 354)
+    diag_op = DiagOp(ones(n))  # square, linear, zero displacement
+    result2 = combine(diag_op, z)
+    @test result2 === z
+end
diff --git a/test/test_gpu_quality.jl b/test/test_gpu_quality.jl
index cc4cf4d..0a5f964 100644
--- a/test/test_gpu_quality.jl
+++ b/test/test_gpu_quality.jl
@@ -1,4 +1,4 @@
-@testitem "GpuExt Quality" tags = [:gpu, :quality] begin
+@testitem "GpuExt Quality" tags = [:quality, :gpu] begin
     using AbstractOperators, JLArrays, LinearAlgebra
 
     # Loading JLArrays triggers GpuExt (GPUArrays is transitive dep)
@@ -28,7 +28,7 @@
     @test y_alloc isa JLArrays.JLArray
 end
 
-@testitem "GpuExt JET" tags = [:gpu, :jet] begin
+@testitem "GpuExt JET" tags = [:jet, :gpu] begin
     using AbstractOperators, JET, JLArrays
 
     # Verify key GPU-dispatched functions are type-stable (no dynamic dispatch)
@@ -55,7 +55,10 @@ end
     @test_call target_modules = (AbstractOperators,) mul!(y10, op_zp, x)
     @test_call target_modules = (AbstractOperators,) mul!(x, op_zp', y10)
 
-    # _should_thread dispatches to false for GPU arrays
+    # _should_thread dispatches to false for GPU arrays (Type and instance overloads)
     @test AbstractOperators._should_thread(typeof(d)) == false
+    @test AbstractOperators._should_thread(d) == false
     @test AbstractOperators._should_thread(Array{Float64}) == (Threads.nthreads() > 1)
+    # storage_type_display_string returns GPU marker for GPU array types
+    @test AbstractOperators.array_type_display_string(typeof(d)) == "ᵍᵖᵘ"
 end
diff --git a/test/test_nonlinear_operators.jl b/test/test_nonlinear_operators.jl
index d13c839..4ad9d15 100644
--- a/test/test_nonlinear_operators.jl
+++ b/test/test_nonlinear_operators.jl
@@ -138,126 +138,126 @@ end
 
 # ─── GPU test items for nonlinear operators ───────────────────────────────────
 
-@testitem "NonlinearOp: Sigmoid (GPU)" tags = [:gpu, :nonlinearoperator, :Sigmoid] setup = [TestUtils] begin
+@testitem "NonlinearOp: Sigmoid (GPU)" tags = [:gpu, :nonlinearoperator, :Sigmoid] setup = [TestUtils, GPUNLTestUtils] begin
     using GPUEnv, Random, AbstractOperators
 
     for backend in gpu_backends()
         Random.seed!(0)
         n = 4
         x = gpu_randn(backend, n)
-        op = Sigmoid(x; gamma = 2.0)  # construct from GPU array to get GPU storage type
+        op = Sigmoid(Float64, (n,), 2.0; array_type = gpu_wrapper(backend, Float64, n))
         test_NLop_gpu(op, x, gpu_randn(backend, n), false)
     end
 end
 
-@testitem "NonlinearOp: SoftMax (GPU)" tags = [:gpu, :nonlinearoperator, :SoftMax] setup = [TestUtils] begin
+@testitem "NonlinearOp: SoftMax (GPU)" tags = [:gpu, :nonlinearoperator, :SoftMax] setup = [TestUtils, GPUNLTestUtils] begin
     using GPUEnv, Random, AbstractOperators
 
     for backend in gpu_backends()
         Random.seed!(0)
         n = 10
         x = gpu_randn(backend, n)
-        op = SoftMax(x)  # construct from GPU array so buffer is GPU-typed
+        op = SoftMax(Float64, (n,); array_type = gpu_wrapper(backend, Float64, n))
         test_NLop_gpu(op, x, gpu_randn(backend, n), false)
     end
 end
 
-@testitem "NonlinearOp: SoftPlus (GPU)" tags = [:gpu, :nonlinearoperator, :SoftPlus] setup = [TestUtils] begin
+@testitem "NonlinearOp: SoftPlus (GPU)" tags = [:gpu, :nonlinearoperator, :SoftPlus] setup = [TestUtils, GPUNLTestUtils] begin
     using GPUEnv, Random, AbstractOperators
 
     for backend in gpu_backends()
         Random.seed!(0)
         n = 10
         x = gpu_randn(backend, n)
-        op = SoftPlus(x)
+        op = SoftPlus(Float64, (n,); array_type = gpu_wrapper(backend, Float64, n))
         test_NLop_gpu(op, x, gpu_randn(backend, n), false)
     end
 end
 
-@testitem "NonlinearOp: Exp (GPU)" tags = [:gpu, :nonlinearoperator, :Exp] setup = [TestUtils] begin
+@testitem "NonlinearOp: Exp (GPU)" tags = [:gpu, :nonlinearoperator, :Exp] setup = [TestUtils, GPUNLTestUtils] begin
     using GPUEnv, Random, AbstractOperators
 
     for backend in gpu_backends()
         Random.seed!(0)
         n, m = 4, 5
         x = gpu_randn(backend, n, m)
-        op = Exp(x)
+        op = Exp(Float64, (n, m); array_type = gpu_wrapper(backend, Float64, n, m))
         test_NLop_gpu(op, x, gpu_randn(backend, n, m), false)
     end
 end
 
-@testitem "NonlinearOp: Sin (GPU)" tags = [:gpu, :nonlinearoperator, :Sin] setup = [TestUtils] begin
+@testitem "NonlinearOp: Sin (GPU)" tags = [:gpu, :nonlinearoperator, :Sin] setup = [TestUtils, GPUNLTestUtils] begin
     using GPUEnv, Random, AbstractOperators
 
     for backend in gpu_backends()
         Random.seed!(0)
         n, m = 4, 5
         x = gpu_randn(backend, n, m)
-        op = Sin(x)
+        op = Sin(Float64, (n, m); array_type = gpu_wrapper(backend, Float64, n, m))
         test_NLop_gpu(op, x, gpu_randn(backend, n, m), false)
     end
 end
 
-@testitem "NonlinearOp: Cos (GPU)" tags = [:gpu, :nonlinearoperator, :Cos] setup = [TestUtils] begin
+@testitem "NonlinearOp: Cos (GPU)" tags = [:gpu, :nonlinearoperator, :Cos] setup = [TestUtils, GPUNLTestUtils] begin
     using GPUEnv, Random, AbstractOperators
 
     for backend in gpu_backends()
         Random.seed!(0)
         n, m = 4, 5
         x = gpu_randn(backend, n, m)
-        op = Cos(x)
+        op = Cos(Float64, (n, m); array_type = gpu_wrapper(backend, Float64, n, m))
         test_NLop_gpu(op, x, gpu_randn(backend, n, m), false)
     end
 end
 
-@testitem "NonlinearOp: Atan (GPU)" tags = [:gpu, :nonlinearoperator, :Atan] setup = [TestUtils] begin
+@testitem "NonlinearOp: Atan (GPU)" tags = [:gpu, :nonlinearoperator, :Atan] setup = [TestUtils, GPUNLTestUtils] begin
     using GPUEnv, Random, AbstractOperators
 
     for backend in gpu_backends()
         Random.seed!(0)
         n = 10
         x = gpu_randn(backend, n)
-        op = Atan(x)
+        op = Atan(Float64, (n,); array_type = gpu_wrapper(backend, Float64, n))
         test_NLop_gpu(op, x, gpu_randn(backend, n), false)
     end
 end
 
-@testitem "NonlinearOp: Tanh (GPU)" tags = [:gpu, :nonlinearoperator, :Tanh] setup = [TestUtils] begin
+@testitem "NonlinearOp: Tanh (GPU)" tags = [:gpu, :nonlinearoperator, :Tanh] setup = [TestUtils, GPUNLTestUtils] begin
     using GPUEnv, Random, AbstractOperators
 
     for backend in gpu_backends()
         Random.seed!(0)
         n = 10
         x = gpu_randn(backend, n)
-        op = Tanh(x)
+        op = Tanh(Float64, (n,); array_type = gpu_wrapper(backend, Float64, n))
         test_NLop_gpu(op, x, gpu_randn(backend, n), false)
     end
 end
 
-@testitem "NonlinearOp: Sech (GPU)" tags = [:gpu, :nonlinearoperator, :Sech] setup = [TestUtils] begin
+@testitem "NonlinearOp: Sech (GPU)" tags = [:gpu, :nonlinearoperator, :Sech] setup = [TestUtils, GPUNLTestUtils] begin
     using GPUEnv, Random, AbstractOperators
 
     for backend in gpu_backends()
         Random.seed!(0)
         n = 10
         x = gpu_randn(backend, n)
-        op = Sech(x)
+        op = Sech(Float64, (n,); array_type = gpu_wrapper(backend, Float64, n))
         test_NLop_gpu(op, x, gpu_randn(backend, n), false)
     end
 end
 
-@testitem "NonlinearOp: Pow (GPU)" tags = [:gpu, :nonlinearoperator, :Pow] setup = [TestUtils] begin
+@testitem "NonlinearOp: Pow (GPU)" tags = [:gpu, :nonlinearoperator, :Pow] setup = [TestUtils, GPUNLTestUtils] begin
     using GPUEnv, Random, AbstractOperators
 
     for backend in gpu_backends()
         Random.seed!(0)
         n = 10
         x = gpu_randn(backend, n)
-        op = Pow(x, 2)
+        op = Pow(Float64, (n,), 2; array_type = gpu_wrapper(backend, Float64, n))
         test_NLop_gpu(op, x, gpu_randn(backend, n), false)
 
         x2 = abs.(gpu_randn(backend, n))
-        op2 = Pow(x2, 0.5)
+        op2 = Pow(Float64, (n,), 0.5; array_type = gpu_wrapper(backend, Float64, n))
         test_NLop_gpu(op2, x2, abs.(gpu_randn(backend, n)), false)
     end
 end
diff --git a/test/test_quality.jl b/test/test_quality.jl
index 8cf3a88..f0e7d83 100644
--- a/test/test_quality.jl
+++ b/test/test_quality.jl
@@ -15,5 +15,6 @@ end
 @testitem "Aqua" tags = [:quality] begin
     using Aqua, AbstractOperators
 
-    Aqua.test_all(AbstractOperators, persistent_tasks = VERSION >= v"1.11")
+    # persistent_tasks is disabled: it is unreliable and not relevant for this package.
+    Aqua.test_all(AbstractOperators, persistent_tasks = false)
 end
diff --git a/test/test_syntax.jl b/test/test_syntax.jl
index 48da474..cc92e4b 100644
--- a/test/test_syntax.jl
+++ b/test/test_syntax.jl
@@ -427,6 +427,15 @@ end
     sliced = S[2:(n - 1)]
     x = randn(n)
     @test sliced * x ≈ (S * x)[2:(n - 1)]
+
+    # Line 71: Sum with multi-domain input (ndoms > 1 → iterates over A.A)
+    # getindex on multi-domain Sum selects domain sub-operators (not codomain rows)
+    n2, m1, m2 = 4, 2, 3
+    Hbase = HCAT(MatrixOp(randn(n2, m1)), MatrixOp(randn(n2, m2)))
+    S_md = Sum(Hbase, Hbase)
+    sliced_md = S_md[1]   # selects first sub-op (domain = m1) from each HCAT member
+    x1 = randn(m1)
+    @test sliced_md * x1 ≈ (Hbase[1] + Hbase[1]) * x1
 end
 
 @testitem "Syntax: Scale getindex (ndoms == 1 branch)" tags = [:misc, :Syntax] setup = [TestUtils] begin
@@ -437,4 +446,69 @@ end
     sliced = s[1:2]
     x = randn(n)
     @test sliced * x ≈ (s * x)[1:2]
+
+    # Line 132: Scale wrapping multi-domain operator
+    # getindex on multi-domain Scale selects domain sub-operators (not codomain rows)
+    n2, m1, m2 = 4, 2, 3
+    Hbase = HCAT(MatrixOp(randn(n2, m1)), MatrixOp(randn(n2, m2)))
+    sc_md = Scale(2.0, Hbase)
+    sc_sliced = sc_md[1]   # selects first sub-op (domain = m1)
+    x1 = randn(m1)
+    @test sc_sliced * x1 ≈ 2.0 .* (Hbase[1] * x1)
+end
+
+@testitem "Syntax: check domain/codomain ArrayPartition errors" tags = [:misc, :Syntax] setup = [TestUtils] begin
+    using AbstractOperators, RecursiveArrayTools
+    n, m1, m2 = 4, 2, 3
+    op_multi_in = HCAT(MatrixOp(randn(n, m1)), MatrixOp(randn(n, m2)))
+    y = zeros(n)
+    # utils.jl:104 — multi-domain op with non-ArrayPartition input
+    @test_throws ArgumentError AbstractOperators.check(y, op_multi_in, randn(m1))
+
+    # utils.jl:128 — multi-codomain op with non-ArrayPartition output
+    op_multi_out = DCAT(MatrixOp(randn(m1, m1)), MatrixOp(randn(m2, m2)))
+    @test_throws ArgumentError AbstractOperators.check(
+        randn(m1), op_multi_out, ArrayPartition(randn(m1), randn(m2))
+    )
+end
+
+@testitem "Syntax: Scale complex coeff on real AdjointMatrixOp errors" tags = [:misc, :Syntax] setup = [TestUtils] begin
+    using AbstractOperators
+    # MatrixOp.jl:84 — real-codomain adjoint MatrixOp scaled by complex scalar
+    n = 4
+    op = MatrixOp(randn(n, n))'
+    @test_throws ErrorException Scale(1.0im, op)
+end
+
+@testitem "Syntax: Compose getindex with multi-domain errors (line 62)" tags = [:misc, :Syntax] setup = [TestUtils] begin
+    using AbstractOperators, RecursiveArrayTools
+    # syntax.jl:62: Compose with ndoms>1 cannot be split (error branch)
+    # diagonal * HCAT always simplifies to HCAT via combination rules,
+    # so any Compose with ndoms>1 has a non-diagonal tail and hits the error.
+    n = 4
+    M = MatrixOp(randn(n, n))
+    H = HCAT(DiagOp(randn(n)), DiagOp(randn(n)))
+    C = M * H   # non-diagonal outer: Compose with ndoms>1 and non-diagonal tail
+    @test ndoms(C, 2) == 2
+    @test_throws ErrorException C[1:2]
+end
+
+@testitem "copy_operator: fast path, slow path, default fallback" tags = [:misc, :Syntax] setup = [TestUtils] begin
+    using Random, AbstractOperators
+    Random.seed!(8)
+
+    n = 5
+    A = randn(n, n)
+    op = MatrixOp(A)
+    @test is_thread_safe(op) == true
+
+    # Fast path: thread-safe operator, no kwargs -> same object shared
+    op_shared = copy_operator(op)
+    @test op_shared === op
+
+    # Slow path: explicit kwarg forces the default (deepcopy) _copy_operator_impl fallback
+    op_copy = copy_operator(op; threaded = true)
+    @test op_copy isa MatrixOp
+    x = randn(n)
+    @test op_copy * x ≈ op * x
 end
diff --git a/test/utils.jl b/test/utils.jl
index de05bad..6d4880b 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -17,15 +17,13 @@
         Pkg.develop(path = normpath(joinpath(@__DIR__, "..", "FFTWOperators"))) # FFTWOperators
         Pkg.develop(path = normpath(joinpath(@__DIR__, "..", "NFFTOperators"))) # NFFTOperators
         Pkg.develop(path = normpath(joinpath(@__DIR__, "..", "WaveletOperators"))) # WaveletOperators
+    elseif VERSION >= v"1.11" && Base.find_package("AcceleratedDCTs") === nothing
+        Pkg.add(name = "AcceleratedDCTs", version = "0.4")
     end
 
     using GPUEnv
     GPUEnv.activate(; persist = true)
 
-    if VERSION >= v"1.11" && Base.find_package("AcceleratedDCTs") === nothing
-        Pkg.add(name = "AcceleratedDCTs", version = "0.4")
-    end
-
     const verb = get(ENV, "ABSTRACTOPERATORS_TEST_VERBOSE", "false") == "true"
 
     to_cpu(x::AbstractArray) = collect(x)
diff --git a/test/wavelets/test_wavelet_operators.jl b/test/wavelets/test_wavelet_operators.jl
index 4235671..0e274c9 100644
--- a/test/wavelets/test_wavelet_operators.jl
+++ b/test/wavelets/test_wavelet_operators.jl
@@ -18,3 +18,18 @@
 
     @test all(norm.(y1 .- y2) .<= 1.0e-12)
 end
+
+@testitem "WaveletOp constructor errors" tags = [:wavelet, :WaveletOp] setup = [TestUtils] begin
+    using Wavelets, WaveletOperators
+    wt = wavelet(WT.db4)
+
+    # 1D: odd dimension
+    @test_throws ArgumentError WaveletOp(Float64, wt, 5)
+    # 1D: too many levels
+    @test_throws ArgumentError WaveletOp(Float64, wt, 8, 100)
+
+    # ND: odd dimension in tuple
+    @test_throws ArgumentError WaveletOp(Float64, wt, (5, 8))
+    # ND: too many levels
+    @test_throws ArgumentError WaveletOp(Float64, wt, (8, 8), 100)
+end