Subtype AbstractArray + LinearAlgebra.Diagonal Support

README.md

@@ -2,26 +2,41 @@
   <img src="docs/src/assets/logo.png" alt="cuNumeric.jl" width="50">
   <a href="https://julialegate.github.io/cuNumeric.jl/dev/">cuNumeric.jl</a>
 </h1>
+<h1>
+  <img src="docs/src/assets/logo.png" alt="cuNumeric.jl" width="50">
+  <a href="https://julialegate.github.io/cuNumeric.jl/dev/">cuNumeric.jl</a>
+</h1>
 
+[![Documentation dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://julialegate.github.io/cuNumeric.jl/dev/) [![codecov](https://codecov.io/github/julialegate/cuNumeric.jl/branch/main/graph/badge.svg)](https://app.codecov.io/github/JuliaLegate/cuNumeric.jl) [![License: MIT](https://img.shields.io/badge/License-MIT-green.svg)](https://opensource.org/licenses/MIT)
 [![Documentation dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://julialegate.github.io/cuNumeric.jl/dev/) [![codecov](https://codecov.io/github/julialegate/cuNumeric.jl/branch/main/graph/badge.svg)](https://app.codecov.io/github/JuliaLegate/cuNumeric.jl) [![License: MIT](https://img.shields.io/badge/License-MIT-green.svg)](https://opensource.org/licenses/MIT)
 
+cuNumeric.jl wraps and extends the [cuPyNumeric](https://github.com/nv-legate/cupynumeric) library from NVIDIA to bring distributed array computing on GPUs and CPUs to Julia. The central type is `NDArray`, which behaves like Julia's `Array` or the `CuArray` from [CUDA.jl](https://github.com/juliagpu/cuda.jl), but executes across multiple GPUs/CPUs. We implement array-level operations on `NDArray` which can be composed into larger programs without the need for explicit MPI calls or writing CUDA kernels.
 cuNumeric.jl wraps and extends the [cuPyNumeric](https://github.com/nv-legate/cupynumeric) library from NVIDIA to bring distributed array computing on GPUs and CPUs to Julia. The central type is `NDArray`, which behaves like Julia's `Array` or the `CuArray` from [CUDA.jl](https://github.com/juliagpu/cuda.jl), but executes across multiple GPUs/CPUs. We implement array-level operations on `NDArray` which can be composed into larger programs without the need for explicit MPI calls or writing CUDA kernels.
 
+cuNumeric.jl requires x86 Linux, an NVIDIA GPU, and Julia >= 1.10. If ARM support is of interest open an issue.
 cuNumeric.jl requires x86 Linux, an NVIDIA GPU, and Julia >= 1.10. If ARM support is of interest open an issue.
 
 ### Quick Start
 
+cuNumeric.jl can be installed with the Julia package manager. Activate your preferred environment and then from the Julia REPL run:
+
+
 cuNumeric.jl can be installed with the Julia package manager. Activate your preferred environment and then from the Julia REPL run:
 
 ```julia
 using Pkg
 Pkg.add(url = "https://github.com/JuliaLegate/cuNumeric.jl", rev = "main")
+using Pkg
+Pkg.add(url = "https://github.com/JuliaLegate/cuNumeric.jl", rev = "main")
 ```
 
 The first time might take awhile as it has to install multiple large dependencies such as the CUDA SDK (if you have an NVIDIA GPU). To use a local build of cupynumeric.so, see [Build Modes](./install.md).
 
+The first time might take awhile as it has to install multiple large dependencies such as the CUDA SDK (if you have an NVIDIA GPU). To use a local build of cupynumeric.so, see [Build Modes](./install.md).
+
 ```julia
 using cuNumeric
+using cuNumeric
 cuNumeric.versioninfo()
 ```
 
@@ -34,41 +49,57 @@ For more details, see [Hardware](./configuration/hardware.md).
 
 The semantics of `NDArray` closely mirror Julia's `Array`, and in most cases it is a drop-in replacement. You can use the same constructors (i.e., `zeros`, `ones`, `rand`), broadcasting, slicing, and linear algebra. Under the hood a few details differ from Base, and knowing them can help you write fast code.
 
-**Data may live across many devices.** An `NDArray` is a logical array whose physical buffers can be partitioned over GPUs and CPUs by the Legate runtime. You write ordinary array code and Legate decides where the data lives and how/when it is communicated between devices. As a result, elementwise indexing (i.e. `arr[1]`) is slow (and is prevented by default). Scalar indexing like this forces synchronization and blocks other tasks from executing.
+**Data may live across many devices.** An `NDArray` is a logical array whose physical buffers can be partitioned over GPUs and CPUs by the Legate runtime. You write ordinary array code and Legate decides where the data lives and how/when it is communicated between devices. As a result, elementwise indexing (i.e. `arr[1]`) is slow (and is prevented by default). Scalar indexing like this forces synchronization and blocks other tasks from executing. Functions like `println` result in data being copied to the host and can also be slow.
 
 **Slices are views.** Indexing an `NDArray` with ranges returns a view onto the same store, not a copy. That differs from Base Julia, where `A[1:n]` allocates a new `Array`. Mutations through an `NDArray` slice are visible through other aliases of the same data.
 
-**Reductions return arrays, not Julia scalars.** Reductions such as `sum(A)` produce a **0D or 1D** `NDArray` (axis reductions produce a lower-rank `NDArray`), rather than a bare `Float64` / `Float32`. That keeps the Legate task graph asynchronous instead of forcing synchronization to communite with the Julia runtime. When you need a plain Julia number, call `unwrap`:
+**Reductions return arrays, not Julia scalars.** Reductions such as `sum(A)` produce a **0D or 1D** `NDArray` (axis reductions produce a lower-rank `NDArray`), rather than a bare `Float64` / `Float32`. That keeps the Legate task graph asynchronous instead of forcing synchronization to communite with the Julia runtime. When you need a plain Julia number, call `unwrap` or `only`:
 
 ```julia
 s = sum(A)          # NDArray{T,0}
 x = unwrap(s)       # T, e.g. Float32
+x2 = only(s)
+```
+s = sum(A)          # NDArray{T,0}
+x = unwrap(s)       # T, e.g. Float32
 ```
 
+**The Legate runtime builds a DAG asynchronously.** Calling `cuNumeric.zeros` or `A .+ B` records work into Legate's task graph rather than blocking until every GPU kernel finishes. Results are materialized when you need them (for example `println`, `unwrap`, or converting with `Array(A)`). Hiding latency enables performant code.
 **The Legate runtime builds a DAG asynchronously.** Calling `cuNumeric.zeros` or `A .+ B` records work into Legate's task graph rather than blocking until every GPU kernel finishes. Results are materialized when you need them (for example `println`, `unwrap`, or converting with `Array(A)`). Hiding latency enables performant code.
 
+For API details see [Initialization](./api_initialization.md) and [NDArray Reference](./api.md). For anti-patterns that kill performance, see [Patterns to Avoid](./perf/patterns_to_avoid.md).
 For API details see [Initialization](./api_initialization.md) and [NDArray Reference](./api.md). For anti-patterns that kill performance, see [Patterns to Avoid](./perf/patterns_to_avoid.md).
 
+### Kernel Fusion
 ### Kernel Fusion
 
+Nested broadcast expressions fuse into a single kernel by default when on GPU. Prefer `@.` for multi-op elementwise code so every operator is dotted and the expression stays completely fused. Even just forgetting the `.` on unary negation (i.e., `y .= -a .+ b .* c`) will result in unfused code. Use the following pattern instead.
 Nested broadcast expressions fuse into a single kernel by default when on GPU. Prefer `@.` for multi-op elementwise code so every operator is dotted and the expression stays completely fused. Even just forgetting the `.` on unary negation (i.e., `y .= -a .+ b .* c`) will result in unfused code. Use the following pattern instead.
 
+```julia
+y .= @. -a + b * c
 ```julia
 y .= @. -a + b * c
 ```
 
 See [Kernel Fusion](./perf/kernel_fusion.md) and [Debugging](./debugging.md) for controls and pretty printers.
 
+See [Kernel Fusion](./perf/kernel_fusion.md) and [Debugging](./debugging.md) for controls and pretty printers.
+
 ### Helping the Garbage Collector
 
+Many calls such as array slicing and un-fused broadcasts allocate a new `NDArray`. The Legate runtime keeps track of all references to the underlying data and will not free the memory until Julia's GC frees the `NDArray` handles. Because Julia's GC runs on memory pressure and an `NDArray` only stores a pointer (i.e., Julia's GC does not know the true size), many dead buffers accumulate and can cause out-of-memory errors.
 Many calls such as array slicing and un-fused broadcasts allocate a new `NDArray`. The Legate runtime keeps track of all references to the underlying data and will not free the memory until Julia's GC frees the `NDArray` handles. Because Julia's GC runs on memory pressure and an `NDArray` only stores a pointer (i.e., Julia's GC does not know the true size), many dead buffers accumulate and can cause out-of-memory errors.
 
+`@analyze_lifetimes` performs a **static last-use analysis** at macro-expansion time and inserts eager calls to immediately free unused `NDArrays`. These buffers can then be reused by legate later for same-sized allocations.
 `@analyze_lifetimes` performs a **static last-use analysis** at macro-expansion time and inserts eager calls to immediately free unused `NDArrays`. These buffers can then be reused by legate later for same-sized allocations.
 
 ```julia
 @analyze_lifetimes begin
     result = @. A[1:end, :] + B[1:end, :]
     C .= @. result * 2.0f0
+    result = @. A[1:end, :] + B[1:end, :]
+    C .= @. result * 2.0f0
 end
 ```
 
@@ -100,3 +131,11 @@ More worked examples (initialization, Gray-Scott, …) are in the documentation
 ### Known Limitations
 
 - There is no support for `Float16` or `ComplexF16`
+- Arrays with 4 or more dimensions might have worse performance
+- Maximum array dimension is 6
+
+
+### Known Deviations from Base Julia
+- Reductions return 0D stores instead of scalars
+- Slices return views
+- `inv` does not throw `SingularException` for singular matrices

TODO.md

@@ -5,4 +5,58 @@
 - Support Ints on methods that takes floats
 - Programatic manipulation of Legate hardware config (not currently possible)
 - Add Aqua.jl to CI to ensure we didn't pirate any types
-- Fix CodeCov reports
+
+## Base
+
+Easy `Base` / `AbstractArray` gaps for `NDArray`. Module helpers
+(`cuNumeric.reshape` / `transpose` / `unique` / …) often exist, but the
+corresponding `Base` methods are missing, so calls fall through to
+`AbstractArray` and may scalar-index. Wire up `Base.*` when convenient.
+
+**P0**
+- `Base.reshape`, `Base.vec`
+- `fill!` (convertible eltypes)
+- `collect`
+
+**P0 done**
+- `iszero` / `isone` (on-device reductions → `Bool` via scalar sync; `isone` square 2D only)
+
+**P1**
+- `transpose` / `adjoint`
+- `unique`
+- `ones_like`
+- `dropdims`
+
+**P2**
+- `extrema`, `mean`
+- `count` / nonzero
+- `abs2`
+- numeric `all` / `any`
+- `sum(f, A)` specials
+
+**P3**
+- `copyto!(NDArray, AbstractArray)`
+- `floor` / `ceil` / `clamp`
+- 2D `permutedims`
+- `diff`
+
+## LinearAlgebra
+
+Starter list of easy/medium LA gaps. Prefer wiring `LinearAlgebra` entry
+points so they do not fall through to scalar-indexing Base paths.
+
+**BLAS-1 style (dense `NDArray`)**
+- `axpy!` / `axpby!` — common scale-and-add; examples use broadcast today
+- `scal!` and related in-place scale
+- `LinearAlgebra.dot` if not already Base-wired to `nda_dot`
+
+**Reductions / traces**
+- `LinearAlgebra.tr` for dense 2D `NDArray` — `cuNumeric.trace` exists; `tr` is
+  already wired for `Diagonal{<:NDArray}`
+
+**Diagonal vs fallthrough (context)**
+- Already on-device for `Diagonal`: `mul!` / `lmul!` / `rmul!`, `\` / `/`,
+  `tr`, `norm` / `opnorm`, many predicates — see `docs/src/linalg.md`
+- Still fallthrough / unsupported on `Diagonal` (e.g. `svd`, `pinv`,
+  `cholesky`, host `AbstractArray` RHS): leave alone unless fixing is cheap;
+  densify intentionally when needed

deps/build.jl

@@ -42,7 +42,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 8`
+    bld_command = `$(joinpath(repo_root, "scripts/build_cpp_wrapper.sh")) $repo_root $cupynumeric_loc $legate_loc $blas_loc $install_root $(Threads.nthreads())`
     return BuildTools.run_build_wrapper_script(
         repo_root, bld_command; cuda_root, cuda_enabled, log_dir=@__DIR__
     )

docs/src/api_initialization.md

@@ -2,49 +2,65 @@
 
 Constructors for new `NDArray`s. Default floating-point type is `Float32`.
 
-## zeros
+## Basic Initialization
+
+### zeros
 
 ```@docs
 cuNumeric.zeros
 ```
 
-## ones
+### ones
 
 ```@docs
 cuNumeric.ones
 ```
 
-## fill
+### fill
 
 ```@docs
 cuNumeric.fill
 ```
 
-## trues
+### trues
 
 ```@docs
 cuNumeric.trues
 ```
 
-## falses
+### falses
 
 ```@docs
 cuNumeric.falses
 ```
 
-## eye
+## Special Matrices
 
-```@docs
-cuNumeric.eye
+### Diagonal
+
+Construct a `Diagonal` matrix whose elements (diagonal only) are stored in an `NDArray`.
+`LinearAlgebra.I` can be used to construct dense identity matrices as well.
+
+```julia
+using LinearAlgebra
+using cuNumeric
+
+D = Diagonal(cuNumeric.ones(Float32, 5))  # preferred for diagonal work
+I32 = NDArray{Float32}(I, 5, 5)           # dense Float32 identity
+Ib = NDArray(I, 5, 5)                     # Bool identity
 ```
 
-## rand
+See [Linear Algebra](./linalg.md#diagonal-and-identity) for preferred patterns.
+
+## Random Numbers
+
+### rand
 
 ```@docs
 cuNumeric.rand
 ```
 
-## rand!
+### rand!
 
 ```@docs
 Random.rand!(::NDArray{<:cuNumeric.SUPPORTED_FLOAT_TYPES})

docs/src/examples/initialization.md

@@ -3,6 +3,7 @@
 Create `NDArray`s with the usual Julia-style constructors. The default element type is `Float32` unless you pass one.
 
 ```julia
+using LinearAlgebra
 using cuNumeric
 
 # Zeros / ones / fill
@@ -15,9 +16,10 @@ F = cuNumeric.fill(7.5f0, (2, 3))
 T = cuNumeric.trues(2, 3)
 Fbool = cuNumeric.falses(2, 3)
 
-# Identity
-I = cuNumeric.eye(5)
-I16 = cuNumeric.eye(Float32, 5)
+# Identity: prefer Diagonal / I; densify only when needed (no eye)
+D = Diagonal(cuNumeric.ones(Float32, 5))
+I32 = NDArray{Float32}(I, 5, 5)   # dense identity
+Ib = NDArray(I, 5, 5)             # Bool identity
 
 # Uniform / normal random values (native Float32 and Float64)
 R = cuNumeric.rand(4, 4)

docs/src/examples/special_mat.md

@@ -0,0 +1,23 @@
+# Special Matrices
+
+Currently cuNumeric only supports `LinearAlgebra.Diagonal`. Other special matrix types like `Tridiagonal` and `Symmetric` will follow.
+
+Diagonal matrices are common, require only storage of the diagonal elements and are often simple to compute operations on (i.e. `LinearAlgebra.inv`). `LinearAlgebra.Diagonal` matrices can be constructed from 1D or 2D `NDArray`s and have certain operations implemented (i.e., `eigen` and `inv`).
+
+
+```julia
+using cuNumeric
+using LinearAlgebra
+
+one_dim = cuNumeric.NDArray([1,2,3,4,5])
+two_dim = cuNumeric.rand(5,5)
+
+D1 = Diagonal(one_dim)
+D2 = Diagonal(two_dim)
+
+evals, evecs = eigen(D1)
+D1_inv = inv(D1)
+
+D1 ./= 2 # stays diagonal
+arr = D2 .+ two_dim # densifies because `two_dim` is not guranteed to be diagonal
+```