diff --git a/benchmark/run.jl b/benchmark/run.jl
index 574bb4ea7..cb748ccde 100644
--- a/benchmark/run.jl
+++ b/benchmark/run.jl
@@ -34,7 +34,7 @@ function ensure_project_ready()
         Pkg.PackageSpec(; path=joinpath(@__DIR__, "..", "lib", "CNPreferences")),
         Pkg.PackageSpec(; path=joinpath(@__DIR__, "..")),
     ])
-    Pkg.instantiate()
+    return Pkg.instantiate()
 end
 
 # default env name mirrors install_cupynumeric.sh: cupynumeric-bench-<major>.<minor>
@@ -46,7 +46,7 @@ function cupynumeric_env_name()
         info.version === nothing && continue
         return "cupynumeric-bench-$(info.version.major).$(info.version.minor)"
     end
-    error("could not resolve cupynumeric_jll version; set CUPYNUMERIC_ENV explicitly")
+    return error("could not resolve cupynumeric_jll version; set CUPYNUMERIC_ENV explicitly")
 end
 
 function dispatch(; gpus, cpus, name, T, N, M, n_iter, n_warmup, n_trial,
diff --git a/benchmark/src/benchmarks/gemm.jl b/benchmark/src/benchmarks/gemm.jl
index a43567926..4f85d3e24 100644
--- a/benchmark/src/benchmarks/gemm.jl
+++ b/benchmark/src/benchmarks/gemm.jl
@@ -8,7 +8,7 @@ dims(g::GEMM) = (g.N, g.M)
 data(g::GEMM{T}) where {T} = "GEMM with T=$(T), N=$(g.N), M=$(g.M)"
 
 function allowed_types(::Type{GEMM})
-    Union{cuNumeric.SUPPORTED_FLOAT_TYPES,cuNumeric.SUPPORTED_INT_TYPES}
+    return Union{cuNumeric.SUPPORTED_FLOAT_TYPES,cuNumeric.SUPPORTED_INT_TYPES}
 end
 
 total_flops(s::GEMM) = s.N * s.N * ((2*s.M) - 1)
diff --git a/benchmark/src/benchmarks/grayscott.jl b/benchmark/src/benchmarks/grayscott.jl
index 29f4151c6..3ba6e6398 100644
--- a/benchmark/src/benchmarks/grayscott.jl
+++ b/benchmark/src/benchmarks/grayscott.jl
@@ -150,7 +150,8 @@ let body = quote
         v_new[end, :] = v[2, :]
     end
     @eval _gs_step!(b::GrayScottBaseline, u, v, u_new, v_new, args::GSParams) = $body
-    @eval _gs_step!(b::GrayScottLifetimes, u, v, u_new, v_new, args::GSParams) = @analyze_lifetimes $body
+    @eval _gs_step!(b::GrayScottLifetimes, u, v, u_new, v_new, args::GSParams) =
+        @analyze_lifetimes $body
 end
 
 function run!(b::AbstractGrayScott, st::GrayScottState)
diff --git a/benchmark/src/benchmarks/montecarlo.jl b/benchmark/src/benchmarks/montecarlo.jl
index 1df91c97f..978e5906f 100644
--- a/benchmark/src/benchmarks/montecarlo.jl
+++ b/benchmark/src/benchmarks/montecarlo.jl
@@ -5,7 +5,7 @@ end
 name(::MonteCarloIntegration) = "montecarlo"
 dims(mci::MonteCarloIntegration) = (mci.n_samples, 1)
 function data(mci::MonteCarloIntegration{T}) where {T}
-    "Monte Carlo Integration with T=$(T), n_samples=$(mci.n_samples)"
+    return "Monte Carlo Integration with T=$(T), n_samples=$(mci.n_samples)"
 end
 
 allowed_types(::Type{MonteCarloIntegration}) = cuNumeric.SUPPORTED_FLOAT_TYPES
@@ -25,7 +25,7 @@ run!(mci::MonteCarloIntegration, x) = _domain_volume(mci) * sum(exp.(-x .^ 2))
 
 # n_samples comes in as N; M is unused.
 function build_benchmark(::Type{MonteCarloIntegration}, ::Type{T}, N, M) where {T}
-    MonteCarloIntegration{T}(; n_samples=N)
+    return MonteCarloIntegration{T}(; n_samples=N)
 end
 
 register_benchmark("montecarlo", MonteCarloIntegration)
diff --git a/benchmark/src/core.jl b/benchmark/src/core.jl
index 5d9b3abbb..526ce4eb1 100644
--- a/benchmark/src/core.jl
+++ b/benchmark/src/core.jl
@@ -43,11 +43,11 @@ function run! end
 # registers itself via `register_benchmark`.
 const BENCHMARKS = Dict{String,Type}()
 function register_benchmark(key::AbstractString, ::Type{B}) where {B<:AbstractBenchmark}
-    BENCHMARKS[key] = B
+    return BENCHMARKS[key] = B
 end
 
 function build_benchmark(::Type{B}, ::Type{T}, N, M) where {B<:AbstractBenchmark,T}
-    B{T}(; N=N, M=M)
+    return B{T}(; N=N, M=M)
 end
 
 #########################################
diff --git a/benchmark/src/parse_benchmarks.jl b/benchmark/src/parse_benchmarks.jl
index cff2eba45..28cad96d4 100644
--- a/benchmark/src/parse_benchmarks.jl
+++ b/benchmark/src/parse_benchmarks.jl
@@ -23,7 +23,7 @@ function parse_fusion(x)
     s = lowercase(string(x))
     s in ("on", "true") && return true
     s in ("off", "false") && return false
-    error("fusion must be on/off (or true/false); got $(repr(x))")
+    return error("fusion must be on/off (or true/false); got $(repr(x))")
 end
 
 # Value of a zipped field for sweep position `i`. length==1 field broadcasts.
diff --git a/benchmark/src/single.jl b/benchmark/src/single.jl
index e37f1d399..56ecd171b 100644
--- a/benchmark/src/single.jl
+++ b/benchmark/src/single.jl
@@ -57,7 +57,7 @@ function run_single(
     @printf("[%s] Mean Run Time: %.5f ± %.5f ms\n", label, mean(br.times_ms), _std(br.times_ms))
     @printf("[%s] FLOPS: %.5f ± %.5f GFLOPS\n", label, mean(br.gflops), _std(br.gflops))
     println("[$(label)] Correctness: $(br.correctness)")
-    save_result(br, gpus; mod=save_as)
+    return save_result(br, gpus; mod=save_as)
 end
 
 gpus = parse(Int, ARGS[1])
diff --git a/deps/build.jl b/deps/build.jl
index 9bcdd2ad0..0e0553adb 100644
--- a/deps/build.jl
+++ b/deps/build.jl
@@ -38,7 +38,7 @@ function build_cpp_wrapper(
     @info "libcunumeric_jl_wrapper: Building C++ Wrapper Library"
     isdir(install_root) && (rm(install_root; recursive=true); mkdir(install_root))
     bld_command = `$(joinpath(repo_root, "scripts/build_cpp_wrapper.sh")) $repo_root $cupynumeric_loc $legate_loc $blas_loc $install_root $(Threads.nthreads())`
-    BuildTools.run_build_wrapper_script(
+    return BuildTools.run_build_wrapper_script(
         repo_root, bld_command; cuda_root, cuda_enabled, log_dir=@__DIR__
     )
 end
@@ -64,7 +64,7 @@ function build_deps(pkg_root, cupynumeric_root, blas_root; cuda_root=nothing, cu
         install_lib;
         cuda_root, cuda_enabled,
     )
-    BuildTools.set_jll_artifact_override(:cunumeric_jl_wrapper_jll, install_lib)
+    return BuildTools.set_jll_artifact_override(:cunumeric_jl_wrapper_jll, install_lib)
 end
 
 function build(::CNPreferences.JLL)
@@ -90,7 +90,7 @@ function build(::CNPreferences.Conda)
     #!TODO SET LocalPreferences.toml to use local CUDA libraries
 
     is_cupynumeric_installed(cupynumeric_root; throw_errors=true)
-    build_deps(pkg_root, cupynumeric_root, cupynumeric_root)
+    return build_deps(pkg_root, cupynumeric_root, cupynumeric_root)
 end
 
 function build(::CNPreferences.Developer)
@@ -110,7 +110,7 @@ function build(::CNPreferences.Developer)
     end
 
     blas_lib = something(blas_lib, BuildTools.find_jll_artifact_dir(:OpenBLAS32_jll))
-    build_deps(pkg_root, cupynumeric_root, up_dir(blas_lib); cuda_root, cuda_enabled)
+    return build_deps(pkg_root, cupynumeric_root, up_dir(blas_lib); cuda_root, cuda_enabled)
 end
 
 const mode_str = load_preference(CNPreferences, "cunumeric_mode", CNPreferences.MODE_JLL)
diff --git a/docs/src/linalg.md b/docs/src/linalg.md
index 12742fad4..103c2c40b 100644
--- a/docs/src/linalg.md
+++ b/docs/src/linalg.md
@@ -1,10 +1,16 @@
 # Linear Algebra
 
-cuNumeric.jl supports a small set of linear algebra operations on `NDArray`. This page covers matrix multiply, batched solve, and related helpers. Related autodocs also appear under [NDArray Reference](./api.md).
+cuNumeric.jl provides matrix multiplication, batched solves, SVD, QR, and related
+helpers for `NDArray`.
+
+`solve`, `svd`, and `qr` accept `Float32`, `Float64`, `ComplexF32`, and
+`ComplexF64`. Integer and `Bool` inputs require `@allowpromotion` or
+`allowpromotion` and produce `Float64` outputs.
 
 ## Matrix multiply
 
-For two 2D arrays, `*` is matrix multiplication (GEMM), not elementwise multiply. Use `.*` when you want an elementwise product of matrices.
+For two 2D arrays, `*` performs matrix multiplication; use `.*` for an
+elementwise product.
 
 ```julia
 using LinearAlgebra
@@ -25,43 +31,50 @@ Filter = t -> t isa Function && nameof(t) === :mul!
 
 ## Solve (batched)
 
-`cuNumeric.solve(A, b)` solves linear systems. It is not Julia's `\`.
-
-```@docs
-cuNumeric.solve
-```
-
-Shapes follow the batched signature:
-
-- `A` is `(..., m, m)` (last two dims square)
-- `b` is `(..., m)` or `(..., m, n)`
-- result is `(..., m)` or `(..., m, n)`
-
-A 1D right-hand side is reshaped internally to a single column, then reshaped back.
+`cuNumeric.solve(A, b)` solves linear systems and returns an array with the same
+shape as `b`. `A` has shape `(..., m, m)` and `b` has shape `(..., m)` or
+`(..., m, n)`.
 
 ```julia
-using cuNumeric
-
-# Single system: (m, m) and (m,)
 A = cuNumeric.rand(Float32, 64, 64)
 b = cuNumeric.rand(Float32, 64)
 x = cuNumeric.solve(A, b)
 
-# Several right-hand sides: (m, m) and (m, n)
 B = cuNumeric.rand(Float32, 64, 4)
 X = cuNumeric.solve(A, B)
 
-# Batched systems: (batch, m, m) and (batch, m, n)
 As = cuNumeric.rand(Float32, 8, 32, 32)
 Bs = cuNumeric.rand(Float32, 8, 32, 2)
 Xs = cuNumeric.solve(As, Bs)
 ```
 
-Notes:
+## Singular value decomposition
+
+`cuNumeric.svd(A, full_matrices=true)` returns `(U, S, Vh)` for a 2D `m × n`
+array. With `k = min(m, n)`, the output shapes are:
+
+- Full: `U` is `m × m`, `S` has length `k`, and `Vh` is `n × n`.
+- Thin: `U` is `m × k`, `S` has length `k`, and `Vh` is `k × n`.
+
+```julia
+A = cuNumeric.rand(Float32, 128, 64)
+U, S, Vh = cuNumeric.svd(A, false)
+```
+
+`S` is real-valued for both real and complex inputs.
+
+## QR decomposition
+
+`cuNumeric.qr(A)` returns the economy-size factors `(Q, R)` for a 2D `m × n`
+array. With `k = min(m, n)`, `Q` is `m × k` and `R` is `k × n`.
+
+```julia
+A = cuNumeric.rand(Float32, 128, 64)
+Q, R = cuNumeric.qr(A)
+```
 
-- Accepted types: `Float32`, `Float64`, `ComplexF32`, `ComplexF64`. Integer or `Bool` inputs promote to `Float64` only when promotion is allowed (`@allowpromotion` / `allowpromotion`).
-- The implementation always goes through a batched Legate `SOLVE` task, including the 2D case.
-- Batch dimensions are supported in the API. Coverage for higher-rank batches in the test suite is still thin, so start with 2D and small batches when validating new code.
+SVD and QR currently accept only 2D arrays; batched decompositions are not
+supported.
 
 ## Helpers
 
@@ -74,4 +87,5 @@ These helpers live on `NDArray` and are also listed in the Public API:
 
 ## Not available yet
 
-There is no public `svd`, `qr`, `cholesky`, `eig`, `lu`, matrix `inv`, or `ldiv!` in cuNumeric.jl yet. Elementwise `inv` / `^-1` exist as unary ops; those are not matrix inverse.
+There is no public `cholesky`, `eig`, `lu`, matrix `inv`, or `ldiv!` yet.
+Elementwise `inv` / `^-1` are unary operations, not matrix inverse.
diff --git a/lib/cunumeric_jl_wrapper/src/types.cpp b/lib/cunumeric_jl_wrapper/src/types.cpp
index f51369d6d..82c0d5016 100644
--- a/lib/cunumeric_jl_wrapper/src/types.cpp
+++ b/lib/cunumeric_jl_wrapper/src/types.cpp
@@ -168,4 +168,10 @@ void wrap_linalg_ops(jlcxx::Module& mod) {
                 legate::LocalTaskID{CuPyNumericOpCode::CUPYNUMERIC_SOLVE});
   mod.set_const("MP_SOLVE",
                 legate::LocalTaskID{CuPyNumericOpCode::CUPYNUMERIC_MP_SOLVE});
+  mod.set_const("SVD", legate::LocalTaskID{CuPyNumericOpCode::CUPYNUMERIC_SVD});
+  mod.set_const("CQR", legate::LocalTaskID{CuPyNumericOpCode::CUPYNUMERIC_QR});
+  mod.set_const("SYEV",
+                legate::LocalTaskID{CuPyNumericOpCode::CUPYNUMERIC_SYEV});
+  mod.set_const("GEEV",
+                legate::LocalTaskID{CuPyNumericOpCode::CUPYNUMERIC_GEEV});
 }
diff --git a/src/cuNumeric.jl b/src/cuNumeric.jl
index 498948382..eb2d0f6b3 100644
--- a/src/cuNumeric.jl
+++ b/src/cuNumeric.jl
@@ -70,7 +70,8 @@ const SUPPORTED_NUMERIC_TYPES = Union{
 
 # solve has no integer backend kernel
 const SUPPORTED_SOLVE_TYPES = Union{SUPPORTED_FLOAT_TYPES,SUPPORTED_COMPLEX_TYPES}
-
+const SUPPORTED_SVD_TYPES = Union{SUPPORTED_FLOAT_TYPES,SUPPORTED_COMPLEX_TYPES}
+const SUPPORTED_QR_TYPES = Union{SUPPORTED_FLOAT_TYPES,SUPPORTED_COMPLEX_TYPES}
 const SUPPORTED_ARRAY_TYPES = Union{Bool,SUPPORTED_NUMERIC_TYPES}
 const SUPPORTED_TYPES = Union{SUPPORTED_ARRAY_TYPES,String}
 
@@ -150,6 +151,7 @@ const TASK_SCOPE_NAMES = CNPreferences.TASK_SCOPE_NAMES
 
 # NDArray internal
 include("ndarray/detail/ndarray.jl")
+include("ndarray/detail/linalg.jl")
 
 # Utilities
 include("cuda/strided_device_array.jl")
diff --git a/src/cuda/cuda_ptx_task.jl b/src/cuda/cuda_ptx_task.jl
index 8e34fcb45..40a08c9ca 100644
--- a/src/cuda/cuda_ptx_task.jl
+++ b/src/cuda/cuda_ptx_task.jl
@@ -69,13 +69,6 @@ function check_sz(arr, maxshape)
     end
 end
 
-# Unused by Launch (which uses `_add_task_array!` + eager finalize). If revived,
-# callers must finalize the returned LogicalArray's handle after add_input/output.
-function nda_to_logical_array(arr::NDArray{T,N}) where {T,N}
-    st_handle = cuNumeric.get_store(arr)
-    return Legate.LogicalArray{T,N}(st_handle, size(arr))
-end
-
 # `get_store` returns a Julia-owned `LogicalArrayImplAllocated` that shares the
 # underlying Legate store with the NDArray. `add_input`/`add_output` copy that
 # array into the task; if we leave the temporary alive until GC, store refcounts
diff --git a/src/cuda/cuda_util.jl b/src/cuda/cuda_util.jl
index 68b141c77..fc2d5ea03 100644
--- a/src/cuda/cuda_util.jl
+++ b/src/cuda/cuda_util.jl
@@ -14,7 +14,7 @@ end
 # Dense @cuda_task / RunPTXTask — MUST match CUDA.jl CuDeviceArray layout.
 # Other memory types: https://github.com/JuliaGPU/CUDA.jl/blob/345c1600ebd561135148bb04ee2657f521a40e25/CUDACore/src/device/pointer.jl#L7
 function ndarray_cuda_type(::Type{<:NDArray{T,N}}) where {T,N}
-    CUDACore.CuDeviceArray{T,N,CUDACore.AS.Global}
+    return CUDACore.CuDeviceArray{T,N,CUDACore.AS.Global}
 end
 
 function ndarray_cuda_type(::Type{T}) where {T}
diff --git a/src/cuda/strided_device_array.jl b/src/cuda/strided_device_array.jl
index 58efd836e..576281b7a 100644
--- a/src/cuda/strided_device_array.jl
+++ b/src/cuda/strided_device_array.jl
@@ -36,12 +36,12 @@ Base.length(a::CuStridedDeviceArray) = a.len
 Base.IndexStyle(::Type{<:CuStridedDeviceArray}) = IndexLinear()
 
 function Base.pointer(a::CuStridedDeviceArray{T,<:Any,A}) where {T,A}
-    Base.unsafe_convert(CUDACore.LLVMPtr{T,A}, a)
+    return Base.unsafe_convert(CUDACore.LLVMPtr{T,A}, a)
 end
 function Base.unsafe_convert(
     ::Type{CUDACore.LLVMPtr{T,A}}, a::CuStridedDeviceArray{T,<:Any,A}
 ) where {T,A}
-    a.ptr
+    return a.ptr
 end
 
 # 0-based element offset from a 1-based linear index in Julia column-major order
@@ -90,9 +90,10 @@ CUDACore.@device_function @inline function _strided_arrayset(
     return A
 end
 
-Base.@propagate_inbounds Base.getindex(A::CuStridedDeviceArray{T}, i::Integer) where {T} = _strided_arrayref(
-    A, i
-)
+Base.@propagate_inbounds Base.getindex(A::CuStridedDeviceArray{T}, i::Integer) where {T} =
+    _strided_arrayref(
+        A, i
+    )
 Base.@propagate_inbounds function Base.setindex!(
     A::CuStridedDeviceArray{T}, x, i::Integer
 ) where {T}
diff --git a/src/memory.jl b/src/memory.jl
index 4a3442940..b71866288 100644
--- a/src/memory.jl
+++ b/src/memory.jl
@@ -22,12 +22,6 @@ const post_gc_host_bytes = Atomic{Int64}(0)
 # how much new memory must accumulate before GC fires again
 const gc_hysteresis_frac = Ref{Float64}(0.05)
 
-# memory measured right after the last GC
-const post_gc_device_bytes = Atomic{Int64}(0)
-const post_gc_host_bytes = Atomic{Int64}(0)
-# how much new memory must accumulate before GC fires again
-const gc_hysteresis_frac = Ref{Float64}(0.05)
-
 @doc"""
     init_gc!()
 
@@ -38,7 +32,7 @@ function init_gc!()
     total_device_bytes[] = query_total_device_memory()
     total_host_bytes[] = query_total_host_memory()
     # @info "[cuNumeric GC] $(total_device_bytes[]) framebuffer available"
-    AUTO_GC_ENABLE[] = true
+    return AUTO_GC_ENABLE[] = true
 end
 
 @doc"""
diff --git a/src/ndarray/detail/linalg.jl b/src/ndarray/detail/linalg.jl
new file mode 100644
index 000000000..0362a84f5
--- /dev/null
+++ b/src/ndarray/detail/linalg.jl
@@ -0,0 +1,222 @@
+function choose_nd_color_shape(shape::NTuple{N,Int}) where {N}
+    color_shape = Base.ones(Int, N)
+    if N > 2
+        color_shape[1] = Legate.num_procs()
+        done = false
+        while !done && color_shape[1] % 2 == 0
+            weight_per_dim = [shape[i] / color_shape[i] for i in 1:(N - 2)]
+            max_weight, idx = findmax(weight_per_dim)
+            if weight_per_dim[idx] > 2 * weight_per_dim[1]
+                color_shape[1] ÷= 2
+                color_shape[idx] *= 2
+            else
+                done = true
+            end
+        end
+    end
+    return Tuple(color_shape)
+end
+
+function prepare_manual_task_for_batched_matrices(full_shape::NTuple{N,Int}) where {N}
+    initial_color_shape = choose_nd_color_shape(full_shape)
+    tilesize = Tuple(
+        (full_shape[i] + initial_color_shape[i] - 1) ÷ initial_color_shape[i] for i in 1:N
+    )
+    color_shape = Tuple((full_shape[i] + tilesize[i] - 1) ÷ tilesize[i] for i in 1:N)
+    return tilesize, color_shape
+end
+
+function solve_batched(a::NDArray{T,N}, b::NDArray, x::NDArray) where {T,N}
+    nrhs = size(b)[end]
+    full_shape = size(a)
+    tilesize_a, color_shape = prepare_manual_task_for_batched_matrices(full_shape)
+    tilesize_b = (tilesize_a[1:(end - 1)]..., nrhs)
+
+    store_a = nda_to_logical_store(a)
+    store_b = nda_to_logical_store(b)
+    store_x = nda_to_logical_store(x)
+
+    tiled_a = Legate.partition_by_tiling(store_a, collect(tilesize_a))
+    tiled_b = Legate.partition_by_tiling(store_b, collect(tilesize_b))
+    tiled_x = Legate.partition_by_tiling(store_x, collect(tilesize_b))
+
+    @task_scope "solve" begin
+        rt = Legate.get_runtime()
+        domain = Legate.domain_from_shape(Legate.Shape(Legate.to_cxx_vector(color_shape)))
+        lib = cuNumeric.get_lib()
+        task = Legate.create_manual_task(rt, lib, cuNumeric.SOLVE, domain)
+
+        Legate.add_input(task, tiled_a)
+        Legate.add_input(task, tiled_b)
+        Legate.add_output(task, tiled_x)
+
+        Legate.submit_manual_task(rt, task)
+    end
+end
+
+# solve runs in floating point:
+# int/bool inputs promote to Float64 (matching cupynumeric)
+const _SOLVE_PROMOTABLE = Union{SUPPORTED_INT_TYPES,Bool}
+const _SOLVE_ACCEPTED = Union{SUPPORTED_SOLVE_TYPES,_SOLVE_PROMOTABLE}
+_solve_eltype(::Type{T}) where {T<:_SOLVE_PROMOTABLE} = Float64
+_solve_eltype(::Type{T}) where {T<:SUPPORTED_SOLVE_TYPES} = T
+
+# `a` must be at least 2D, `b` at least 1D.
+function _solve_check_a_dims(a::NDArray{<:Any,0}, b::NDArray)
+    throw(ArgumentError("0-dimensional array given. Array must be at least two-dimensional"))
+end
+function _solve_check_a_dims(a::NDArray{<:Any,1}, b::NDArray)
+    throw(ArgumentError("1-dimensional array given. Array must be at least two-dimensional"))
+end
+_solve_check_a_dims(a::NDArray, b::NDArray) = _solve_check_b_dims(a, b)
+
+function _solve_check_b_dims(a::NDArray, b::NDArray{<:Any,0})
+    throw(ArgumentError("0-dimensional array given. Array must be at least one-dimensional"))
+end
+_solve_check_b_dims(a::NDArray, b::NDArray) = _solve(a, b)
+
+# 2D case: (m,m),(m)->(m).
+# Backend needs rhs "b" to be 2D. We reshape b from (n,) to (n,1)
+function _solve(a::NDArray{T,2}, b::NDArray{S,1}) where {T,S}
+    m = size(b)[1]
+    return reshape(_solve(a, reshape(b, (m, 1))), (m,))
+end
+
+# 2D (m,m),(m,n)->(m,n) and batched (...,m,m),(...,m,n)->(...,m,n)
+function _solve(a::NDArray{T,N}, b::NDArray{S,N}) where {T,S,N}
+    size(a)[end - 1] != size(a)[end] &&
+        throw(ArgumentError("Last 2 dimensions of the array must be square"))
+    size(a)[end] != size(b)[end - 1] &&
+        throw(
+            ArgumentError(
+                "Input operand 1 has a mismatch in its dimension " *
+                "$(N-2), with signature (...,m,m),(...,m,n)->(...,m,n)" *
+                " (size $(size(b)[end-1]) is different from $(size(a)[end]))",
+            ),
+        )
+    prod(size(a)) == 0 || prod(size(b)) == 0 && return zeros(T, size(b)...)
+    x = zeros(T, size(b)...)
+    solve_batched(a, b, x)
+    return x
+end
+
+# Mismatched batch dimensions
+function _solve(a::NDArray{T,N}, b::NDArray{S,M}) where {T,N,S,M}
+    throw(ArgumentError("Batched matrices require signature (...,m,m),(...,m,n)->(...,m,n)"))
+end
+
+function svd_single(a::NDArray{T,N}, u::NDArray, s::NDArray, vh::NDArray) where {T,N}
+    rt = Legate.get_runtime()
+    lib = cuNumeric.get_lib()
+    task = Legate.create_auto_task(rt, lib, cuNumeric.SVD)
+
+    l_a = nda_to_logical_array(a)
+    l_u = nda_to_logical_array(u)
+    l_s = nda_to_logical_array(s)
+    l_vh = nda_to_logical_array(vh)
+
+    Legate.add_input(task, l_a)
+    Legate.add_output(task, l_u)
+    Legate.add_output(task, l_s)
+    Legate.add_output(task, l_vh)
+
+    Legate.add_broadcast(task, l_a)
+    Legate.add_broadcast(task, l_u)
+    Legate.add_broadcast(task, l_s)
+    Legate.add_broadcast(task, l_vh)
+
+    return Legate.submit_auto_task(rt, task)
+end
+
+function _svd(a::NDArray{T,2}, full_matrices::Bool) where {T}
+    m, n = size(a)
+    k = min(m, n)
+    S = real(T)
+    # cuSolver requires full square buffers regardless of full_matrices
+    u_buf = zeros(T, m, m)
+    s = zeros(S, k)
+    vh_buf = zeros(T, n, n)
+    svd_single(a, u_buf, s, vh_buf)
+    # Backend factors are logically ordered; only thin strided views need materialization.
+    u = full_matrices ? u_buf : copy(u_buf[:, 1:k])
+    vh = full_matrices ? vh_buf : copy(vh_buf[1:k, :])
+    return u, s, vh
+end
+
+# svd runs on float/complex only — no integer backend
+const _SVD_PROMOTABLE = Union{SUPPORTED_INT_TYPES,Bool}
+const _SVD_ACCEPTED = Union{SUPPORTED_SVD_TYPES,_SVD_PROMOTABLE}
+_svd_eltype(::Type{T}) where {T<:_SVD_PROMOTABLE} = Float64
+_svd_eltype(::Type{T}) where {T<:SUPPORTED_SVD_TYPES} = T
+
+function _svd_check_dims(a::NDArray{<:Any,0}, full_matrices::Bool)
+    throw(ArgumentError("0-dimensional array given. Array must be at least two-dimensional"))
+end
+
+function _svd_check_dims(a::NDArray{<:Any,1}, full_matrices::Bool)
+    throw(ArgumentError("1-dimensional array given. Array must be at least two-dimensional"))
+end
+
+function _svd_check_dims(a::NDArray{<:Any,2}, full_matrices::Bool)
+    return _svd(a, full_matrices)
+end
+
+function _svd_check_dims(a::NDArray, full_matrices::Bool)
+    throw(ArgumentError("cuNumeric does not yet support stacked 2d arrays"))
+end
+
+# qr
+
+function qr_single(a::NDArray{T,N}, q::NDArray, r::NDArray) where {T,N}
+    rt = Legate.get_runtime()
+    lib = cuNumeric.get_lib()
+    task = Legate.create_auto_task(rt, lib, cuNumeric.CQR)
+
+    l_a = nda_to_logical_array(a)
+    l_q = nda_to_logical_array(q)
+    l_r = nda_to_logical_array(r)
+
+    Legate.add_input(task, l_a)
+    Legate.add_output(task, l_q)
+    Legate.add_output(task, l_r)
+
+    Legate.add_broadcast(task, l_a)
+    Legate.add_broadcast(task, l_q)
+    Legate.add_broadcast(task, l_r)
+
+    return Legate.submit_auto_task(rt, task)
+end
+
+function _qr(a::NDArray{T,2}) where {T}
+    m, n = size(a)
+    k = min(m, n)
+    # cuSolver requires full square buffers regardless of output shape
+    q_buf = zeros(T, m, m)
+    r_buf = zeros(T, n, n)
+    qr_single(a, q_buf, r_buf)
+    # Host conversion assumes contiguous storage, so materialize the economy slices.
+    q = copy(q_buf[:, 1:k])
+    r = copy(r_buf[1:k, :])
+    return q, r
+end
+
+const _QR_PROMOTABLE = Union{SUPPORTED_INT_TYPES,Bool}
+const _QR_ACCEPTED = Union{SUPPORTED_QR_TYPES,_QR_PROMOTABLE}
+_qr_eltype(::Type{T}) where {T<:_QR_PROMOTABLE} = Float64
+_qr_eltype(::Type{T}) where {T<:SUPPORTED_QR_TYPES} = T
+
+function _qr_check_dims(a::NDArray{<:Any,0})
+    throw(ArgumentError("0-dimensional array given. Array must be at least two-dimensional"))
+end
+
+function _qr_check_dims(a::NDArray{<:Any,1})
+    throw(ArgumentError("1-dimensional array given. Array must be at least two-dimensional"))
+end
+
+function _qr_check_dims(a::NDArray{<:Any,2})
+    return _qr(a)
+end
+
+function _qr_check_dims(a::NDArray)
+    throw(ArgumentError("cuNumeric does not yet support stacked 2d arrays"))
+end
diff --git a/src/ndarray/detail/ndarray.jl b/src/ndarray/detail/ndarray.jl
index dfae001fc..c8232e63a 100644
--- a/src/ndarray/detail/ndarray.jl
+++ b/src/ndarray/detail/ndarray.jl
@@ -599,3 +599,8 @@ function nda_to_logical_store(arr::NDArray{T,N}) where {T,N}
     finalize(la_handle)
     return Legate.LogicalStore{T,N}(st_handle, size(arr))
 end
+
+function nda_to_logical_array(arr::NDArray{T,N}) where {T,N}
+    st_handle = cuNumeric.get_store(arr)
+    return Legate.LogicalArray{T,N}(st_handle, size(arr))
+end
diff --git a/src/ndarray/linalg.jl b/src/ndarray/linalg.jl
index 6893b0ca5..075820f85 100644
--- a/src/ndarray/linalg.jl
+++ b/src/ndarray/linalg.jl
@@ -1,75 +1,3 @@
-function choose_nd_color_shape(shape::NTuple{N,Int}) where {N}
-    color_shape = Base.ones(Int, N)
-    if N > 2
-        color_shape[1] = Legate.num_procs()
-        done = false
-        while !done && color_shape[1] % 2 == 0
-            weight_per_dim = [shape[i] / color_shape[i] for i in 1:(N - 2)]
-            max_weight, idx = findmax(weight_per_dim)
-            if weight_per_dim[idx] > 2 * weight_per_dim[1]
-                color_shape[1] ÷= 2
-                color_shape[idx] *= 2
-            else
-                done = true
-            end
-        end
-    end
-    return Tuple(color_shape)
-end
-
-function prepare_manual_task_for_batched_matrices(full_shape::NTuple{N,Int}) where {N}
-    initial_color_shape = choose_nd_color_shape(full_shape)
-    tilesize = Tuple(
-        (full_shape[i] + initial_color_shape[i] - 1) ÷ initial_color_shape[i] for i in 1:N
-    )
-    color_shape = Tuple((full_shape[i] + tilesize[i] - 1) ÷ tilesize[i] for i in 1:N)
-    return tilesize, color_shape
-end
-
-function solve_batched(a::NDArray{T,N}, b::NDArray, x::NDArray) where {T,N}
-    nrhs = size(b)[end]
-    full_shape = size(a)
-    tilesize_a, color_shape = prepare_manual_task_for_batched_matrices(full_shape)
-    tilesize_b = (tilesize_a[1:(end - 1)]..., nrhs)
-
-    store_a = nda_to_logical_store(a)
-    store_b = nda_to_logical_store(b)
-    store_x = nda_to_logical_store(x)
-
-    tiled_a = Legate.partition_by_tiling(store_a, collect(tilesize_a))
-    tiled_b = Legate.partition_by_tiling(store_b, collect(tilesize_b))
-    tiled_x = Legate.partition_by_tiling(store_x, collect(tilesize_b))
-    # Same Legate Julia-wrapper pin class as Launch `_add_task_array!` / get_store
-    # temps: one could finalize store_/tiled_ handles here after partition/add_*
-    # copies ownership into the task. Not enabled yet — weak linalg test coverage.
-    # finalize(store_a.handle)
-    # finalize(store_b.handle)
-    # finalize(store_x.handle)
-
-    @task_scope "solve" begin
-        rt = Legate.get_runtime()
-        domain = Legate.domain_from_shape(Legate.Shape(Legate.to_cxx_vector(color_shape)))
-        lib = cuNumeric.get_lib()
-        task = Legate.create_manual_task(rt, lib, cuNumeric.SOLVE, domain)
-
-        Legate.add_input(task, tiled_a)
-        # finalize(tiled_a.handle)
-        Legate.add_input(task, tiled_b)
-        # finalize(tiled_b.handle)
-        Legate.add_output(task, tiled_x)
-        # finalize(tiled_x.handle)
-
-        Legate.submit_manual_task(rt, task)
-    end
-end
-
-# solve runs in floating point:
-# int/bool inputs promote to Float64 (matching cupynumeric)
-const _SOLVE_PROMOTABLE = Union{SUPPORTED_INT_TYPES,Bool}
-const _SOLVE_ACCEPTED = Union{SUPPORTED_SOLVE_TYPES,_SOLVE_PROMOTABLE}
-_solve_eltype(::Type{T}) where {T<:_SOLVE_PROMOTABLE} = Float64
-_solve_eltype(::Type{T}) where {T<:SUPPORTED_SOLVE_TYPES} = T
-
 # Type/dim guards dispatch on one argument at a time, then forward to `_solve`.
 """
     cuNumeric.solve(A, b)
@@ -97,46 +25,24 @@ function solve(a::NDArray, b::NDArray)
     throw(ArgumentError("array type $bad is unsupported in solve"))
 end
 
-# `a` must be at least 2D, `b` at least 1D.
-function _solve_check_a_dims(a::NDArray{<:Any,0}, b::NDArray)
-    throw(ArgumentError("0-dimensional array given. Array must be at least two-dimensional"))
-end
-function _solve_check_a_dims(a::NDArray{<:Any,1}, b::NDArray)
-    throw(ArgumentError("1-dimensional array given. Array must be at least two-dimensional"))
-end
-_solve_check_a_dims(a::NDArray, b::NDArray) = _solve_check_b_dims(a, b)
-
-function _solve_check_b_dims(a::NDArray, b::NDArray{<:Any,0})
-    throw(ArgumentError("0-dimensional array given. Array must be at least one-dimensional"))
+function svd(a::NDArray{<:_SVD_ACCEPTED}, full_matrices::Bool=true)
+    A = eltype(a)
+    O = _svd_eltype(A)
+    A <: _SVD_PROMOTABLE && assertpromotion(svd, A, O)
+    return _svd_check_dims(unchecked_promote_arr(a, O), full_matrices)
 end
-_solve_check_b_dims(a::NDArray, b::NDArray) = _solve(a, b)
 
-# 2D case: (m,m),(m)->(m).
-# Backend needs rhs "b" to be 2D. We reshape b from (n,) to (n,1)
-function _solve(a::NDArray{T,2}, b::NDArray{S,1}) where {T,S}
-    m = size(b)[1]
-    return reshape(_solve(a, reshape(b, (m, 1))), (m,))
+function svd(a::NDArray, full_matrices::Bool=true)
+    throw(ArgumentError("array type $(eltype(a)) is unsupported in svd"))
 end
 
-# 2D (m,m),(m,n)->(m,n) and batched (...,m,m),(...,m,n)->(...,m,n)
-function _solve(a::NDArray{T,N}, b::NDArray{S,N}) where {T,S,N}
-    size(a)[end - 1] != size(a)[end] &&
-        throw(ArgumentError("Last 2 dimensions of the array must be square"))
-    size(a)[end] != size(b)[end - 1] &&
-        throw(
-            ArgumentError(
-                "Input operand 1 has a mismatch in its dimension " *
-                "$(N-2), with signature (...,m,m),(...,m,n)->(...,m,n)" *
-                " (size $(size(b)[end-1]) is different from $(size(a)[end]))",
-            ),
-        )
-    prod(size(a)) == 0 || prod(size(b)) == 0 && return zeros(T, size(b)...)
-    x = zeros(T, size(b)...)
-    solve_batched(a, b, x)
-    return x
+function qr(a::NDArray{<:_QR_ACCEPTED})
+    A = eltype(a)
+    O = _qr_eltype(A)
+    A <: _QR_PROMOTABLE && assertpromotion(qr, A, O)
+    return _qr_check_dims(unchecked_promote_arr(a, O))
 end
 
-# Mismatched batch dimensions
-function _solve(a::NDArray{T,N}, b::NDArray{S,M}) where {T,N,S,M}
-    throw(ArgumentError("Batched matrices require signature (...,m,m),(...,m,n)->(...,m,n)"))
+function qr(a::NDArray)
+    throw(ArgumentError("array type $(eltype(a)) is unsupported in qr"))
 end
diff --git a/src/ndarray/ndarray.jl b/src/ndarray/ndarray.jl
index d06a65bef..e058a014d 100644
--- a/src/ndarray/ndarray.jl
+++ b/src/ndarray/ndarray.jl
@@ -146,7 +146,7 @@ end
 function (::Type{<:Array{A}})(arr::NDArray{B,0}) where {A,B}
     out = Array{A}(undef)
     allowscalar() do
-        out[] = convert(A, arr[])
+        return out[] = convert(A, arr[])
     end
     return out
 end
@@ -279,33 +279,33 @@ Base.IndexStyle(::NDArray) = IndexCartesian()
 
 function Base.show(io::IO, arr::NDArray{T,0}) where {T}
     allowscalar() do
-        print(io, "NDArray{$(T),0}(", repr(arr[]), ")")
+        return print(io, "NDArray{$(T),0}(", repr(arr[]), ")")
     end
 end
 
 function Base.show(io::IO, ::MIME"text/plain", arr::NDArray{T,0}) where {T}
     println(io, "0-dimensional NDArray{$(T),0}")
     allowscalar() do
-        print(io, arr[])
+        return print(io, arr[])
     end
 end
 
 function Base.show(io::IO, arr::NDArray{T,N}) where {T,N}
-    print(io, "NDArray{$(T),$(N)} with size ", size(arr))
+    return print(io, "NDArray{$(T),$(N)} with size ", size(arr))
 end
 
 function Base.show(io::IO, ::MIME"text/plain", arr::NDArray{T,N}) where {T,N}
     println(io, "NDArray{$(T),$(N)} with size ", size(arr))
-    Base.print_array(io, Array(arr))
+    return Base.print_array(io, Array(arr))
 end
 
 function Base.print(arr::NDArray{T}) where {T}
-    Base.show(stdout, arr)
+    return Base.show(stdout, arr)
 end
 
 function Base.println(arr::NDArray{T}) where {T}
     Base.show(stdout, arr)
-    print("\n")
+    return print("\n")
 end
 #### ARRAY INDEXING AND SLICES ####
 
@@ -370,39 +370,39 @@ end
 #! TODO SUPPORT CONVERSION OF VALUES
 function Base.setindex!(arr::NDArray{T,N}, value::T, idxs::Vararg{Int,N}) where {T,N}
     assertscalar("setindex!")
-    _setindex!(Val{N}(), arr, value, idxs...)
+    return _setindex!(Val{N}(), arr, value, idxs...)
 end
 
 function Base.setindex!(arr::NDArray{Complex{T},N}, value::T, idxs::Vararg{Int,N}) where {T,N}
     assertscalar("setindex!")
-    _setindex!(Val{N}(), arr, Complex{T}(value), idxs...)
+    return _setindex!(Val{N}(), arr, Complex{T}(value), idxs...)
 end
 
 function Base.setindex!(arr::NDArray{T,N}, value, idxs::Vararg{Int,N}) where {T,N}
     assertscalar("setindex!")
-    _setindex!(Val{N}(), arr, convert(T, value), idxs...)
+    return _setindex!(Val{N}(), arr, convert(T, value), idxs...)
 end
 
 function _setindex!(::Val{0}, arr::NDArray{T,0}, value::T) where {T<:SUPPORTED_NUMERIC_TYPES}
     acc = NDArrayAccessor{T,1}()
-    write(acc, arr.ptr, StdVector(UInt64[0]), value)
+    return write(acc, arr.ptr, StdVector(UInt64[0]), value)
 end
 
 function _setindex!(::Val{0}, arr::NDArray{Bool,0}, value::Bool)
     acc = NDArrayAccessor{CxxWrap.CxxBool,1}()
-    write(acc, arr.ptr, StdVector(UInt64[0]), value)
+    return write(acc, arr.ptr, StdVector(UInt64[0]), value)
 end
 
 function _setindex!(
     ::Val{N}, arr::NDArray{T,N}, value::T, idxs::Vararg{Int,N}
 ) where {T<:SUPPORTED_NUMERIC_TYPES,N}
     acc = NDArrayAccessor{T,N}()
-    write(acc, arr.ptr, to_cpp_index(idxs), value)
+    return write(acc, arr.ptr, to_cpp_index(idxs), value)
 end
 
 function _setindex!(::Val{N}, arr::NDArray{Bool,N}, value::Bool, idxs::Vararg{Int,N}) where {N}
     acc = NDArrayAccessor{CxxWrap.CxxBool,N}()
-    write(acc, arr.ptr, to_cpp_index(idxs), value)
+    return write(acc, arr.ptr, to_cpp_index(idxs), value)
 end
 
 #### START OF SLICING ####
@@ -417,35 +417,35 @@ function _setindex_slice!(lhs::NDArray, rhs::NDArray, slices)
 end
 
 function Base.setindex!(lhs::NDArray, rhs::NDArray, i::Colon, j::Int64)
-    _setindex_slice!(lhs, rhs, slice_array((0, Base.size(lhs, 1)), (j-1, j)))
+    return _setindex_slice!(lhs, rhs, slice_array((0, Base.size(lhs, 1)), (j-1, j)))
 end
 
 function Base.setindex!(lhs::NDArray, rhs::NDArray, i::Int64, j::Colon)
-    _setindex_slice!(lhs, rhs, slice_array((i-1, i)))
+    return _setindex_slice!(lhs, rhs, slice_array((i-1, i)))
 end
 
 function Base.setindex!(lhs::NDArray, rhs::NDArray, i::UnitRange, j::Colon)
-    _setindex_slice!(
+    return _setindex_slice!(
         lhs, rhs, slice_array((first(i) - 1, last(i)), (0, Base.size(lhs, 2)))
     )
 end
 
 function Base.setindex!(lhs::NDArray, rhs::NDArray, i::Colon, j::UnitRange)
-    _setindex_slice!(
+    return _setindex_slice!(
         lhs, rhs, slice_array((0, Base.size(lhs, 1)), (first(j) - 1, last(j)))
     )
 end
 
 function Base.setindex!(lhs::NDArray, rhs::NDArray, i::UnitRange, j::Int64)
-    _setindex_slice!(lhs, rhs, slice_array((first(i) - 1, last(i)), (j-1, j)))
+    return _setindex_slice!(lhs, rhs, slice_array((first(i) - 1, last(i)), (j-1, j)))
 end
 
 function Base.setindex!(lhs::NDArray, rhs::NDArray, i::Int64, j::UnitRange)
-    _setindex_slice!(lhs, rhs, slice_array((i-1, i), (first(j) - 1, last(j))))
+    return _setindex_slice!(lhs, rhs, slice_array((i-1, i), (first(j) - 1, last(j))))
 end
 
 function Base.setindex!(lhs::NDArray, rhs::NDArray, i::UnitRange, j::UnitRange)
-    _setindex_slice!(
+    return _setindex_slice!(
         lhs, rhs, slice_array((first(i) - 1, last(i)), (first(j) - 1, last(j)))
     )
 end
@@ -492,19 +492,19 @@ end
 
 Base.getindex(arr::NDArray{T}, c::Vararg{Colon,N}) where {T,N} = Base.copy(arr)
 function Base.setindex!(arr::NDArray{T}, rhs::NDArray{T}, c::Vararg{Colon,N}) where {T,N}
-    Base.copyto!(arr, rhs)
+    return Base.copyto!(arr, rhs)
 end
 
 function Base.setindex!(arr::NDArray{T,2}, val::T, i::Colon, j::Int64) where {T}
     s = nda_get_slice(arr, slice_array((0, Base.size(arr, 1)), (j-1, j)))
     nda_fill_array(s, val)
-    destroy!(s)
+    return destroy!(s)
 end
 
 function Base.setindex!(arr::NDArray{T,2}, val::T, i::Int64, j::Colon) where {T}
     s = nda_get_slice(arr, slice_array((i-1, i)))
     nda_fill_array(s, val)
-    destroy!(s)
+    return destroy!(s)
 end
 
 Base.fill!(arr::NDArray{T}, val::T) where {T} = nda_fill_array(arr, val)
@@ -653,7 +653,7 @@ Fill `arr` in-place with uniform random `Float64` values.
 """
 Random.rand!(arr::NDArray{Float64}) = cuNumeric.nda_random(arr, 0)
 function Random.rand!(arr::NDArray{T}) where {T}
-    error("rand! only supports NDArray{Float64} for now. Cast with cuNumeric.as_type.")
+    return error("rand! only supports NDArray{Float64} for now. Cast with cuNumeric.as_type.")
 end
 
 # Backend only generates Float64. Same-type path needs no cast; other floats
diff --git a/src/ndarray/promotion.jl b/src/ndarray/promotion.jl
index b9dcbac0e..b13437fbb 100644
--- a/src/ndarray/promotion.jl
+++ b/src/ndarray/promotion.jl
@@ -34,12 +34,12 @@ end
 function __checked_promote_op(
     f::typeof(Base.literal_pow), a::Type{Tuple{_,ARR_TYPE,Val{-1}}}
 ) where {_,ARR_TYPE}
-    __recip_type(ARR_TYPE)
+    return __recip_type(ARR_TYPE)
 end
 function __checked_promote_op(
     f::typeof(Base.literal_pow), a::Type{Tuple{_,ARR_TYPE,Val{2}}}
 ) where {_,ARR_TYPE}
-    ARR_TYPE
+    return ARR_TYPE
 end
 
 __checked_promote_op(::typeof(Base.inv), ::Type{Tuple{A}}) where {A} = __recip_type(A)
@@ -74,7 +74,7 @@ end
 # Base.promote_op(^, Float32, Int64) == Float32
 # Base.promote_op(^, Int32, Int64) == Int32
 function __my_promote_type(::Type{typeof(^)}, ::Type{A}, ::Type{Val{V}}) where {A,V}
-    __checked_promote_op(Base.:(^), A, typeof(V))
+    return __checked_promote_op(Base.:(^), A, typeof(V))
 end
 
 #! Not exaclty Julia behavior, but it it makes life easier...
diff --git a/src/util.jl b/src/util.jl
index 5ca776267..17a39b7af 100644
--- a/src/util.jl
+++ b/src/util.jl
@@ -33,7 +33,7 @@ Insert a Legate execution fence. `block=true` waits until prior ops finish;
 issue_execution_fence(; block::Bool=false) = Legate.issue_execution_fence(block)
 
 function Experimental(setting::Bool)
-    task_local_storage(:Experimental, setting)
+    return task_local_storage(:Experimental, setting)
 end
 
 function assert_experimental()
diff --git a/test/tests/linalg.jl b/test/tests/linalg.jl
index 8d5ededb5..10703f7c2 100644
--- a/test/tests/linalg.jl
+++ b/test/tests/linalg.jl
@@ -185,3 +185,176 @@ end
         end
     end
 end
+
+
+function check_svd_reconstruction(ref_A::AbstractMatrix, u, s, vh, tol_a, tol_r)
+    U  = Array(u)
+    S  = Array(s)
+    Vh = Array(vh)
+    A_rec = U * Diagonal(S) * Vh
+    return isapprox(ref_A, A_rec; atol=tol_a, rtol=tol_r)
+end
+
+function check_svd_orthonormality(u, vh, tol_a, tol_r)
+    U  = Array(u)
+    Vh = Array(vh)
+    ku = size(U, 2)
+    kv = size(Vh, 1)
+    ok_u  = isapprox(U'  * U,  Matrix{eltype(U)}(I, ku, ku);  atol=tol_a, rtol=tol_r)
+    ok_vh = isapprox(Vh * Vh', Matrix{eltype(Vh)}(I, kv, kv); atol=tol_a, rtol=tol_r)
+    return ok_u && ok_vh
+end
+
+@testset "svd square matrix" begin
+    @testset verbose=true for T in Base.uniontypes(cuNumeric.SUPPORTED_SVD_TYPES)
+        A_ref = my_rand(T, 5, 5)
+        nda   = cuNumeric.NDArray(A_ref)
+        u, s, vh = cuNumeric.svd(nda)
+        allowscalar() do
+            @test check_svd_reconstruction(A_ref, u, s, vh, atol(T), rtol(T))
+            @test check_svd_orthonormality(u, vh, atol(T), rtol(T))
+        end
+    end
+end
+
+@testset "svd tall matrix (m > n)" begin
+    @testset verbose=true for T in Base.uniontypes(cuNumeric.SUPPORTED_SVD_TYPES)
+        A_ref = my_rand(T, 6, 4)
+        nda   = cuNumeric.NDArray(A_ref)
+        u, s, vh = cuNumeric.svd(nda, false)  # thin SVD for reconstruction test
+        allowscalar() do
+            @test check_svd_reconstruction(A_ref, u, s, vh, atol(T), rtol(T))
+            @test check_svd_orthonormality(u, vh, atol(T), rtol(T))
+        end
+    end
+end
+ 
+@testset "svd thin output shapes (full_matrices=false)" begin
+    @testset verbose=true for T in Base.uniontypes(cuNumeric.SUPPORTED_SVD_TYPES)
+        m, n = 6, 4
+        k     = min(m, n)
+        A_ref = my_rand(T, m, n)
+        nda   = cuNumeric.NDArray(A_ref)
+        u, s, vh = cuNumeric.svd(nda, false)
+        allowscalar() do
+            @test size(Array(u))  == (m, k)
+            @test size(Array(s))  == (k,)
+            @test size(Array(vh)) == (k, n)
+            @test check_svd_reconstruction(A_ref, u, s, vh, atol(T), rtol(T))
+        end
+    end
+end
+ 
+@testset "svd full output shapes (full_matrices=true)" begin
+    @testset verbose=true for T in Base.uniontypes(cuNumeric.SUPPORTED_SVD_TYPES)
+        m, n = 6, 4
+        A_ref = my_rand(T, m, n)
+        nda   = cuNumeric.NDArray(A_ref)
+        u, s, vh = cuNumeric.svd(nda, true)
+        allowscalar() do
+            @test size(Array(u))  == (m, m)
+            @test size(Array(s))  == (min(m, n),)
+            @test size(Array(vh)) == (n, n)
+        end
+    end
+end
+ 
+@testset "svd singular values non-negative and sorted" begin
+    @testset verbose=true for T in Base.uniontypes(cuNumeric.SUPPORTED_SVD_TYPES)
+        A_ref = my_rand(T, 5, 5)
+        nda   = cuNumeric.NDArray(A_ref)
+        _, s, _ = cuNumeric.svd(nda)
+        allowscalar() do
+            sv = Array(s)
+            @test all(sv .>= 0)
+            @test issorted(sv; rev=true)
+        end
+    end
+end
+ 
+@testset "svd identity matrix" begin
+    @testset verbose=true for T in Base.uniontypes(cuNumeric.SUPPORTED_SVD_TYPES)
+        n     = 4
+        A_ref = Matrix{T}(I, n, n)
+        nda   = cuNumeric.NDArray(A_ref)
+        _, s, _ = cuNumeric.svd(nda)
+        allowscalar() do
+            @test cuNumeric.compare(ones(T, n), s, atol(T), rtol(T))
+        end
+    end
+end
+ 
+@testset "svd rank-1 matrix" begin
+    @testset verbose=true for T in Base.uniontypes(cuNumeric.SUPPORTED_SVD_TYPES)
+        # outer product of two vectors: exactly one nonzero singular value
+        # 5x4 satisfies the M >= N constraint
+        v1    = T.(collect(1:5))
+        v2    = T.(collect(1:4))
+        A_ref = v1 * v2'
+        nda   = cuNumeric.NDArray(A_ref)
+        _, s, _ = cuNumeric.svd(nda)
+        allowscalar() do
+            sv = Array(s)
+            @test sv[1] > atol(T)
+            @test all(sv[2:end] .< sqrt(atol(T)) * 100)
+        end
+    end
+end
+
+@testset "svd promotion" begin
+    @testset verbose=true for T in (Int32, Int64, Bool)
+        vals = T == Bool ? T[1 0; 0 1] : reshape(T.(collect(1:4)), 2, 2)
+        A = cuNumeric.NDArray(vals)
+        @test_throws "Implicit promotion" cuNumeric.svd(A)
+        allowpromotion() do
+            u, s, vh = cuNumeric.svd(A)
+            @test eltype(Array(u)) == Float64
+        end
+    end
+end
+
+@testset "qr reconstruction" begin
+    @testset verbose=true for T in Base.uniontypes(cuNumeric.SUPPORTED_QR_TYPES)
+        A_ref = my_rand(T, 6, 4)
+        nda   = cuNumeric.NDArray(A_ref)
+        q, r  = cuNumeric.qr(nda)
+        allowscalar() do
+            Q = Array(q)
+            R = Array(r)
+            @test size(Q) == (6, 4)
+            @test size(R) == (4, 4)
+            @test isapprox(A_ref, Q * R; atol=atol(T), rtol=rtol(T))
+            @test isapprox(Q' * Q, Matrix{eltype(Q)}(I, 4, 4); atol=atol(T), rtol=rtol(T))
+        end
+    end
+end
+
+@testset "qr square matrix" begin
+    @testset verbose=true for T in Base.uniontypes(cuNumeric.SUPPORTED_QR_TYPES)
+        A_ref = my_rand(T, 5, 5)
+        nda   = cuNumeric.NDArray(A_ref)
+        q, r  = cuNumeric.qr(nda)
+        allowscalar() do
+            Q = Array(q)
+            R = Array(r)
+            @test size(Q) == (5, 5)
+            @test size(R) == (5, 5)
+            @test isapprox(A_ref, Q * R; atol=atol(T), rtol=rtol(T))
+        end
+    end
+end
+
+@testset "qr promotion" begin
+    @testset verbose=true for T in (Int32, Int64, Bool)
+        vals = T == Bool ? T[1 0; 0 1] : reshape(T.(collect(1:4)), 2, 2)
+        A = cuNumeric.NDArray(vals)
+        @test_throws "Implicit promotion" cuNumeric.qr(A)
+        allowpromotion() do
+            q, r = cuNumeric.qr(A)
+            allowscalar() do
+                @test eltype(Array(q)) == Float64
+                @test isapprox(Float64.(vals), Array(q) * Array(r); atol=atol(Float64), rtol=rtol(Float64))
+            end
+        end
+    end
+end
\ No newline at end of file
diff --git a/test/tests/stability.jl b/test/tests/stability.jl
index 8eaa15a2d..b000bd522 100644
--- a/test/tests/stability.jl
+++ b/test/tests/stability.jl
@@ -119,6 +119,35 @@ end
     end
 end
 
+@testset verbose = true "svd" begin
+    @testset "$(T)" for T in Base.uniontypes(cuNumeric.SUPPORTED_SVD_TYPES)
+        A = cuNumeric.NDArray(T[1 0; 0 1])
+        @test @inferred(cuNumeric.svd(A)) !== nothing
+        @test @inferred(cuNumeric.svd(A, false)) !== nothing
+    end
+
+    @testset "promote $(T)" for T in (Int32, Int64, Bool)
+        A = cuNumeric.NDArray(T[1 0; 0 1])
+        allowpromotion() do
+            @test @inferred(cuNumeric.svd(A)) !== nothing
+        end
+    end
+end
+
+@testset verbose = true "qr" begin
+    @testset "$(T)" for T in Base.uniontypes(cuNumeric.SUPPORTED_QR_TYPES)
+        A = cuNumeric.NDArray(T[1 0; 0 1])
+        @test @inferred(cuNumeric.qr(A)) !== nothing
+    end
+
+    @testset "promote $(T)" for T in (Int32, Int64, Bool)
+        A = cuNumeric.NDArray(T[1 0; 0 1])
+        allowpromotion() do
+            @test @inferred(cuNumeric.qr(A)) !== nothing
+        end
+    end
+end
+
 @testset verbose = true "linalg ops" begin
     @testset "$(T)" for T in Base.uniontypes(cuNumeric.SUPPORTED_NUMERIC_TYPES)
         M = cuNumeric.zeros(T, 4, 3)
diff --git a/test/tests/util.jl b/test/tests/util.jl
index a477d0576..17b52e2e6 100644
--- a/test/tests/util.jl
+++ b/test/tests/util.jl
@@ -44,7 +44,7 @@ is_same(arr1::Array, arr2::NDArray) = @allowscalar (arr1 == arr2)[1]
 is_same(arr1::Array, arr2::Array) = (arr1 == arr2)
 
 function my_rand(::Type{F}, dims...; L=F(-1000), R=F(1000)) where {F<:AbstractFloat}
-    L .+ (R-L) .* rand(F, dims...)
+    return L .+ (R-L) .* rand(F, dims...)
 end
 function my_rand(::Type{I}, dims...; L=nothing, R=nothing) where {I<:Integer}
     L_default = I <: Unsigned ? 0 : max(-255, Int64(typemin(I)))
@@ -55,7 +55,7 @@ function my_rand(::Type{I}, dims...; L=nothing, R=nothing) where {I<:Integer}
     return floor.(I, res)
 end
 function my_rand(::Type{CT}, dims...; L=T(-100), R=T(100)) where {T,CT<:Complex{T}}
-    Complex.(my_rand(T, dims...; L=L, R=R), my_rand(T, dims...; L=L, R=R))
+    return Complex.(my_rand(T, dims...; L=L, R=R), my_rand(T, dims...; L=L, R=R))
 end
 my_rand(::Type{Bool}, dims...) = rand(Bool, dims...)