diff --git a/gpu_metal_fill/.gitignore b/gpu_metal_fill/.gitignore
new file mode 100644
index 0000000..f0a60b4
--- /dev/null
+++ b/gpu_metal_fill/.gitignore
@@ -0,0 +1,4 @@
+build/
+*.o
+*.gds
+*.ppm
diff --git a/gpu_metal_fill/Makefile b/gpu_metal_fill/Makefile
new file mode 100644
index 0000000..1c20490
--- /dev/null
+++ b/gpu_metal_fill/Makefile
@@ -0,0 +1,67 @@
+# Metal Fill (FEOL/BEOL) — recursive-partitioned dummy fill.
+#
+# Build is CPU (OpenMP-parallel across partitions). The compute-heavy stages are
+# isolated behind FillBackend (see include/metalfill/backend.hpp) so a CUDA
+# backend can be dropped in later (GPU offload is future work).
+
+CXX      ?= g++
+OPT      ?= -O2
+CXXFLAGS ?= -std=c++17 $(OPT) -Wall -Wextra -Iinclude
+LDFLAGS  ?=
+
+# OpenMP (parallel partitions / fill). Disable with `make OPENMP=0`.
+OPENMP ?= 1
+ifeq ($(OPENMP),1)
+  CXXFLAGS += -fopenmp
+  LDFLAGS  += -fopenmp
+endif
+
+BUILD := build
+LIB_SRC := $(wildcard src/*.cpp)
+LIB_OBJ := $(patsubst src/%.cpp,$(BUILD)/%.o,$(LIB_SRC))
+
+TOOLS := $(BUILD)/make_dummy_gds $(BUILD)/make_gpu_block $(BUILD)/run_fill $(BUILD)/render_layer \
+         $(BUILD)/gdsinfo
+TESTS := $(BUILD)/test_main
+
+.PHONY: all tools test demo clean
+all: tools
+tools: $(TOOLS)
+
+$(BUILD):
+	mkdir -p $(BUILD)
+
+$(BUILD)/%.o: src/%.cpp | $(BUILD)
+	$(CXX) $(CXXFLAGS) -c $< -o $@
+
+$(BUILD)/make_dummy_gds: tools/make_dummy_gds.cpp $(LIB_OBJ) | $(BUILD)
+	$(CXX) $(CXXFLAGS) $< $(LIB_OBJ) -o $@ $(LDFLAGS)
+
+$(BUILD)/make_gpu_block: tools/make_gpu_block.cpp $(LIB_OBJ) | $(BUILD)
+	$(CXX) $(CXXFLAGS) $< $(LIB_OBJ) -o $@ $(LDFLAGS)
+
+$(BUILD)/run_fill: tools/run_fill.cpp $(LIB_OBJ) | $(BUILD)
+	$(CXX) $(CXXFLAGS) $< $(LIB_OBJ) -o $@ $(LDFLAGS)
+
+$(BUILD)/render_layer: tools/render_layer.cpp $(LIB_OBJ) | $(BUILD)
+	$(CXX) $(CXXFLAGS) $< $(LIB_OBJ) -o $@ $(LDFLAGS)
+
+$(BUILD)/gdsinfo: tools/gdsinfo.cpp $(LIB_OBJ) | $(BUILD)
+	$(CXX) $(CXXFLAGS) $< $(LIB_OBJ) -o $@ $(LDFLAGS)
+
+$(BUILD)/test_main: tests/test_main.cpp $(LIB_OBJ) | $(BUILD)
+	$(CXX) $(CXXFLAGS) $< $(LIB_OBJ) -o $@ $(LDFLAGS)
+
+test: $(TESTS)
+	$(BUILD)/test_main
+
+demo: tools
+	$(BUILD)/make_dummy_gds -o $(BUILD)/dummy.gds
+	$(BUILD)/run_fill -i $(BUILD)/dummy.gds -o $(BUILD)/filled.gds -r $(BUILD)/report.txt
+
+gpu-demo: tools
+	$(BUILD)/make_gpu_block -o $(BUILD)/gpu_block.gds
+	$(BUILD)/run_fill -i $(BUILD)/gpu_block.gds -o $(BUILD)/gpu_filled.gds -r $(BUILD)/gpu_report.txt
+
+clean:
+	rm -rf $(BUILD)
diff --git a/gpu_metal_fill/README.md b/gpu_metal_fill/README.md
new file mode 100644
index 0000000..61ff050
--- /dev/null
+++ b/gpu_metal_fill/README.md
@@ -0,0 +1,111 @@
+# Metal Fill (FEOL / BEOL) — recursive-partitioned dummy fill
+
+> For a thorough design write-up (motivation, every algorithm, the partitioning
+> correctness argument, bugs fixed, and results) see **[WHITEPAPER.md](WHITEPAPER.md)**
+> (also available as a PDF: [docs/metal_fill_whitepaper.pdf](docs/metal_fill_whitepaper.pdf)).
+
+A small, dependency-free C++17 engine that inserts **dummy metal fill** into a
+GDSII layout so that every layer meets CMP density rules. It processes the die by
+**recursive quadtree partitioning** (like commercial fill flows) so leaves can be
+filled independently and in parallel, then merged.
+
+> Status: CPU implementation (OpenMP-parallel across partitions). The heavy
+> stages are isolated behind a `FillBackend` interface so a CUDA/GPU backend can
+> be added later. GPU offload is intentionally future work.
+
+## Why dummy fill?
+
+Chemical-Mechanical Polishing (CMP) needs a **uniform** metal density across the
+die. Sparse regions dish/erode differently from dense regions, hurting yield and
+timing. Foundries therefore require, per layer and per density window:
+
+- a **minimum** density (add fill where too sparse),
+- a **maximum** density (don't over-fill),
+- a bounded **window-to-window gradient** (no abrupt steps),
+- **spacing** from existing geometry (a keep-out halo), and a min fill area.
+
+## Pipeline
+
+For each layer (FEOL base: `OD`, `PO`; BEOL metals: `M1..M14`):
+
+1. **Rasterize** existing geometry into an occupancy grid.
+2. **Density** via a summed-area table (integral image): O(1) per window.
+3. **Keep-out** = morphological dilation of geometry by the spacing rule.
+4. **Fill** on a pitch lattice, per window, up to a per-window target density,
+   never exceeding max and never entering the keep-out halo.
+5. **Iterative DRC**: measure achieved density, then raise the target of any
+   window that is still below `min` or too far below a denser neighbor
+   (gradient), and refill. Targets increase monotonically, so it converges.
+
+### Recursive partitioning (the scaling trick)
+
+The die is split by a **quadtree** into leaf partitions that are filled
+independently (OpenMP now, GPU later) and then merged by concatenation.
+
+Correctness at partition boundaries is guaranteed by two things:
+
+- **Lattice-aligned cuts**: partitions are cut on a lattice equal to
+  `lcm(density_window, fill_pitch)`, so a density window and every fill position
+  is identical whether computed globally or per-partition.
+- **Halo / guard band**: each leaf reads a ring of neighboring geometry as
+  read-only context (so density and keep-out are exact at the edge) but only
+  *emits* fill inside its **core**. Cores are disjoint and tile the die, so
+  merging is a plain concatenation — no de-duplication or stitching.
+
+A unit test (`partition == global`) asserts the merged result is bit-for-bit
+equivalent to a single-shot global fill.
+
+## Build & run
+
+```sh
+make            # builds tools into build/
+make test       # builds and runs the unit tests
+make demo       # simple synthetic layout -> filled.gds + report.txt
+make gpu-demo   # realistic GPU-block layout -> gpu_filled.gds + gpu_report.txt
+
+# manual
+build/make_gpu_block -o build/gpu_block.gds            # synthetic GPU block
+build/run_fill -i build/gpu_block.gds -o build/gpu_filled.gds -r build/report.txt
+build/render_layer -i build/gpu_filled.gds -o build/m1.ppm -l 10   # visualize M1
+# convert the PPM to PNG if you like: ffmpeg -i build/m1.ppm build/m1.png
+
+# confirm fill happened: per-(layer,datatype) polygon counts. Fill = datatype 10.
+build/gdsinfo build/gpu_block.gds build/gpu_filled.gds
+```
+
+### The GPU-block test design
+
+`make_gpu_block` generates a synthetic — but structurally realistic — GPU block
+(it is **not** a real design, just the floorplan-level density structure):
+
+- an `NxN` grid of SM (streaming-multiprocessor) tiles, each with two SRAM
+  macros (register file + shared memory) and a standard-cell logic region,
+- routing channels between tiles (bus routing on mid metals),
+- global clock/spine routes (M9/M10), and
+- a regular, sparse power grid on the top metals (M11..M14).
+
+Each layer therefore has a distinct density profile — dense macros/logic on the
+lower layers, near-empty upper layers with only power straps — so fill has varied
+work on every layer.
+
+`make OPENMP=0` builds a serial version.
+
+## Layout of the code
+
+| File | Responsibility |
+|------|----------------|
+| `include/metalfill/geometry.hpp` | points, bbox, polygons (integer dbu) |
+| `gdsii.{hpp,cpp}` | minimal GDSII reader/writer (BOUNDARY/BOX) |
+| `raster.{hpp,cpp}` | polygon scanline rasterization |
+| `density.{hpp,cpp}` | summed-area table + windowed density |
+| `fill.{hpp,cpp}` | keep-out dilation + per-window fill placement |
+| `partition.{hpp,cpp}` | recursive quadtree partitioner |
+| `drc.{hpp,cpp}` | density / gradient / spacing / min-area checks |
+| `backend.hpp`, `backend_cpu.cpp` | compute backend (CPU/OpenMP; CUDA later) |
+| `engine.{hpp,cpp}` | orchestration: partition → fill → iterate → merge |
+| `layermap.{hpp,cpp}` | default GDS layers + fill rules |
+| `tools/` | `make_dummy_gds`, `make_gpu_block`, `run_fill`, `render_layer`, `gdsinfo` |
+| `tests/test_main.cpp` | unit + end-to-end tests |
+
+The default layer map and rules are illustrative, **not** a real PDK; edit
+`layermap.cpp` to match your technology.
diff --git a/gpu_metal_fill/WHITEPAPER.md b/gpu_metal_fill/WHITEPAPER.md
new file mode 100644
index 0000000..e6cceb2
--- /dev/null
+++ b/gpu_metal_fill/WHITEPAPER.md
@@ -0,0 +1,600 @@
+# A Recursive-Partitioned Metal Fill Engine for FEOL/BEOL Dummy Fill
+
+**A design white paper for the `gpu_metal_fill` project**
+
+---
+
+## Table of contents
+
+1. [Background: why dummy fill exists](#1-background-why-dummy-fill-exists)
+2. [Problem statement and requirements](#2-problem-statement-and-requirements)
+3. [System architecture](#3-system-architecture)
+4. [Data model and GDSII I/O](#4-data-model-and-gdsii-io)
+5. [Rasterization](#5-rasterization)
+6. [Density analysis with summed-area tables](#6-density-analysis-with-summed-area-tables)
+7. [Keep-out (spacing) as morphological dilation](#7-keep-out-spacing-as-morphological-dilation)
+8. [Fill placement](#8-fill-placement)
+9. [Recursive partitioning and merge](#9-recursive-partitioning-and-merge)
+10. [Gradient-aware, DRC-driven iterative fill](#10-gradient-aware-drc-driven-iterative-fill)
+11. [Design-rule checking](#11-design-rule-checking)
+12. [Parallelism and the GPU-ready backend](#12-parallelism-and-the-gpu-ready-backend)
+13. [The synthetic GPU-block test design](#13-the-synthetic-gpu-block-test-design)
+14. [Bugs encountered and how they were fixed](#14-bugs-encountered-and-how-they-were-fixed)
+15. [Results](#15-results)
+16. [Limitations and future work](#16-limitations-and-future-work)
+17. [File-by-file reference](#17-file-by-file-reference)
+
+---
+
+## 1. Background: why dummy fill exists
+
+Modern chips are manufactured layer by layer. After each metal layer is
+deposited, the wafer is planarized by **Chemical-Mechanical Polishing (CMP)**.
+CMP does not remove material uniformly: regions with a **high** metal density
+polish differently from **sparse** regions. Two failure modes dominate:
+
+- **Dishing** — wide metal features get polished *below* the target height.
+- **Erosion** — in dense arrays of fine features, both metal and the surrounding
+  dielectric erode.
+
+The result is **thickness variation** across the die, which changes wire
+resistance/capacitance (hurting timing), can open or short layers, and lowers
+yield. To keep CMP uniform, foundries impose **density rules** per layer, checked
+over a sliding **window** (e.g. 20 µm × 20 µm):
+
+- a **minimum** density (typically ~20–30 %),
+- a **maximum** density (typically ~70–85 %), and
+- a bounded **window-to-window gradient** (no abrupt steps).
+
+Designs rarely satisfy the *minimum* everywhere, so tools insert non-functional
+**dummy fill** (a.k.a. metal fill) — small floating shapes in the empty space —
+until each window meets its target. Fill must stay a legal **spacing** away from
+real geometry (a *keep-out halo*) and each fill shape must meet a **minimum
+area**.
+
+- **FEOL** (Front-End-Of-Line) = transistor layers (active/`OD`, poly/`PO`).
+  Fill here is for pattern-density and mechanical-stress uniformity.
+- **BEOL** (Back-End-Of-Line) = the interconnect metal stack (`M1..M14`) and
+  vias. This is where most fill volume lives, and where fill runtime dominates —
+  which is exactly why this project focuses on making BEOL fill **fast** via
+  partitioning (and, later, GPU offload).
+
+This project implements a self-contained, dependency-free engine that performs
+this fill for both FEOL base layers and a 14-layer BEOL metal stack.
+
+---
+
+## 2. Problem statement and requirements
+
+> **Input:** a GDSII layout with base layers and metal layers `M1..M14`.
+> **Output:** the same layout plus dummy fill so every layer meets its CMP
+> density rules, with a DRC report.
+
+Concretely the engine must, per layer:
+
+1. Measure existing metal **density** over a sliding window.
+2. Know the **layer map** (which GDS layer/datatype is which) and the per-layer
+   **fill rules**.
+3. **Place fill** to satisfy min density without violating max density, gradient,
+   spacing, or min-area.
+4. Run **iterative DRC**: measure, fix, repeat.
+
+An explicit non-functional requirement drove the whole architecture: **BEOL fill
+is slow, so it must be parallelizable.** The chosen mechanism is **recursive
+spatial partitioning** (divide the die, fill pieces independently, merge) — the
+same idea commercial tools use — with a clean path to GPU offload later.
+
+---
+
+## 3. System architecture
+
+The engine is a small C++17 library plus a few command-line tools. It has **no
+third-party dependencies** (only the standard library; OpenMP is optional).
+
+```mermaid
+flowchart TD
+    GDS[GDSII in] --> RD[read_gds]
+    RD --> ENG[engine: run_fill]
+    subgraph perlayer[per layer]
+      PART[recursive quadtree partition] --> LEAF
+      subgraph LEAF[per leaf - parallel]
+        RAS[rasterize halo] --> KO[keep-out dilation]
+        KO --> FILL[place fill in core]
+      end
+      LEAF --> MERGE[merge core fills]
+      MERGE --> ITER{DRC ok? / converged?}
+      ITER -- no --> PART
+    end
+    ENG --> WR[write_gds] --> OUT[GDSII out]
+    ENG --> REP[text report]
+```
+
+The compute-heavy stages (density, keep-out, fill) are hidden behind a
+`FillBackend` interface (`include/metalfill/backend.hpp`) so they can run on the
+CPU today (OpenMP) or a GPU tomorrow (CUDA) without touching the orchestration.
+
+**Coordinate convention.** Everything internal is in integer **database units
+(dbu)**. The GDS `UNITS` record is written as `1 dbu = 1 nm` (`meters_per_dbu =
+1e-9`) with `1 user unit = 1 µm`, so `1 µm = 1000 dbu`. Rules are authored in
+microns and converted to dbu at runtime.
+
+---
+
+## 4. Data model and GDSII I/O
+
+### Geometry (`geometry.hpp`)
+
+- `Point{dbu x,y}`, `BBox` with `expand()/width()/height()/valid()`.
+- `Polygon{int layer, datatype; vector<Point> pts}` with `bbox()`.
+- `make_rect(layer, datatype, x0,y0,x1,y1)` — the workhorse, since fill shapes
+  and the synthetic designs are rectangles.
+
+### GDSII reader/writer (`gdsii.{hpp,cpp}`)
+
+GDSII is a record-based **big-endian** binary format. Each record is
+`[2-byte length][1-byte record-type][1-byte data-type][payload]`. The reader is a
+straightforward record loop that imports `BOUNDARY` (and `BOX`) elements as
+polygons and ignores references/paths; the writer emits a single top cell.
+
+Two subtleties are handled explicitly:
+
+- **8-byte GDS reals** (used by the `UNITS` record) are *not* IEEE-754. They are
+  `sign(1 bit) · mantissa/2^56 · 16^(exponent-64)`. `put_real8`/`get_real8`
+  implement the conversion.
+- **Closing vertex**: GDS boundaries repeat the first point as the last; the
+  reader drops the duplicate, the writer re-adds it.
+
+The in-memory `Layout` holds the library/cell name, the `UNITS`, and the polygon
+list, and exposes `dbu_per_um()` and `bbox()`.
+
+Because this is a from-scratch reader/writer, a **round-trip unit test**
+(`test_gds_roundtrip`) guards it: write a layout, read it back, and check that
+polygon count, layer/datatype, coordinates, and units are preserved.
+
+---
+
+## 5. Rasterization
+
+Density and spacing are far cheaper to reason about on a **raster** than on
+polygons, so each layer is rendered into a boolean **occupancy grid**.
+
+### The grid (`raster.hpp`)
+
+`Grid` is a dense `nx × ny` array of `uint8_t` (0/1) with a cell size `cell`
+(dbu) and an origin `(ox, oy)`. Cell `(ix,iy)` covers
+`[ox+ix·cell, ox+(ix+1)·cell) × [oy+iy·cell, oy+(iy+1)·cell)`.
+
+`make_grid(area, cell)` allocates a grid that covers `area` using **exact ceil
+sizing** (`nx = ceil(width/cell)`). (An earlier version added a `+1` guard cell;
+see [§14](#14-bugs-encountered-and-how-they-were-fixed) for why it was removed.)
+
+### Scanline polygon fill (`rasterize`)
+
+For each polygon, for each grid **row**, the scanline is taken at the row's
+vertical **center** `yc`. We collect the x-coordinates where polygon edges cross
+`yc`, sort them, and fill the spans between consecutive crossing pairs
+(even-odd rule). A cell is set if its **center** lies inside a span:
+
+```
+cx0 = ceil ((xl - ox)/cell - 0.5)
+cx1 = floor((xr - ox)/cell - 0.5)
+```
+
+This is exact for axis-aligned rectangles (all shapes here) and robust for
+arbitrary simple polygons. `test_rasterize_rect` checks a 500×500 dbu rectangle
+at 100-dbu cells produces exactly 25 occupied cells.
+
+The choice of `cell` matters and is derived per layer — see
+[§8](#8-fill-placement).
+
+---
+
+## 6. Density analysis with summed-area tables
+
+CMP density is evaluated over many overlapping windows, so a naïve per-window sum
+would be `O(window_area)` each. Instead we build a **summed-area table (SAT)**,
+a.k.a. integral image, once, and answer any window in `O(1)`.
+
+### The SAT (`density.{hpp,cpp}`)
+
+For occupancy `occ[x,y] ∈ {0,1}`, the SAT is
+
+```
+SAT[y+1][x+1] = occ[x,y] + SAT[y][x+1] + SAT[y+1][x] − SAT[y][x]
+```
+
+The occupied area of any half-open cell rectangle `[x0,x1) × [y0,y1)` is then
+
+```
+area = SAT[y1][x1] − SAT[y0][x1] − SAT[y1][x0] + SAT[y0][x0]
+```
+
+`build_sat` computes it with a single running-row pass (cache friendly), and
+`SummedAreaTable::area()` clamps to bounds and returns the four-corner
+difference.
+
+### The density map (`compute_density`)
+
+A `DensityMap` is produced by sliding a `win_cells` window with a `step_cells`
+stride. Each `Window` records its tile index, its cell range, and its
+`density = occupied_area / window_area`. `min/mean/max_density()` summarize a
+layer. `test_sat_density` verifies both the SAT arithmetic and a 2×2 tiling on a
+hand-checked grid.
+
+Why the SAT matters here: it is the *same* primitive used for keep-out dilation
+([§7](#7-keep-out-spacing-as-morphological-dilation)) and it maps cleanly onto a
+GPU (prefix sums + constant-time box lookups), which is the intended offload.
+
+---
+
+## 7. Keep-out (spacing) as morphological dilation
+
+Fill must stay `keepout` microns away from any real geometry. Equivalently, the
+**forbidden region** for fill is the existing geometry **dilated** by the keep-out
+radius. Rather than a per-cell neighborhood scan, this is computed in `O(1)` per
+cell with the SAT: a cell `(ix,iy)` is *blocked* iff the box
+`[ix−r, ix+r] × [iy−r, iy+r]` contains any occupied cell, where `r =
+ceil(keepout/cell)`:
+
+```
+blocked[ix,iy] = (SAT.area(ix−r, iy−r, ix+r+1, iy+r+1) > 0)
+```
+
+This is `compute_keepout` in `fill.cpp`. `test_keepout` checks that `r=1` blocks
+the 8-neighborhood of a single occupied cell, that a distance-2 cell is free,
+and that `r=0` reduces to the occupancy itself.
+
+---
+
+## 8. Fill placement
+
+### Deriving the grid resolution and pitch (per layer)
+
+`compute_geom` (in `engine.cpp`) converts a layer's rules to dbu and picks a grid
+`cell` equal to the **gcd** of the fill width/height, pitch, window and step.
+That guarantees each of those lengths is an **integer number of cells**, so
+window boundaries and the fill lattice are exact (no rounding drift). Fill shapes
+are integer multiples of the cell.
+
+### Placement objective (`place_fill`)
+
+Fill candidates sit on a **pitch lattice** anchored at the grid origin. For each
+window (tile) the algorithm:
+
+1. reads the existing occupied cells via the SAT,
+2. computes a **target** occupied-cell count and a **max** cap, and
+3. walks the candidate lattice inside the tile, stamping a fill shape wherever
+   its footprint is entirely free of the keep-out mask and of already-placed
+   fill, stopping once the window reaches its target or the max cap.
+
+Two invariants make this correct and parallel-safe:
+
+- **Footprints stay inside their tile** (`cx + fw ≤ tile_end`), so each tile
+  writes a disjoint region of the fill grid → the per-tile loop is parallelized
+  with OpenMP with no data races.
+- **Never exceed max**: placement stops before crossing `max_density · area`.
+
+Fill-to-fill spacing is guaranteed *by construction*: the pitch is `1.25×` the
+fill size, so shapes on the lattice never touch. `test_place_fill` verifies fill
+reaches the target, respects the max cap, and places nothing where fully blocked.
+
+The per-window target is not a single scalar; it comes from a **target map**
+built by the DRC-driven loop in [§10](#10-gradient-aware-drc-driven-iterative-fill).
+
+---
+
+## 9. Recursive partitioning and merge
+
+This is the core of the project and the reason it can scale.
+
+### The idea
+
+Fill in one region of the die is *almost* independent of another. So we split the
+die with a **quadtree** into leaf partitions, fill each leaf independently (in
+parallel now, on a GPU later), and **merge**. `partition_recursive`
+(`partition.{hpp,cpp}`) recurses in integer lattice-units, splitting into four
+children until a leaf is `≤ max_leaf` units per side (or a depth cap is hit).
+
+```mermaid
+flowchart TD
+    A["die (all windows)"] --> B1[quadrant NW]
+    A --> B2[quadrant NE]
+    A --> B3[quadrant SW]
+    A --> B4[quadrant SE]
+    B1 --> C1[leaf]
+    B1 --> C2[leaf]
+    B2 --> D[... recurse until <= max_leaf ...]
+```
+
+### The catch: boundaries
+
+A naïve cut breaks two things at partition edges:
+
+1. **Density windows** that straddle a cut would see only half their geometry.
+2. **Spacing** (fill-to-geometry and fill-to-fill) across the cut.
+
+If you ignore this, you get exactly the failure the author hit first: fill placed
+on top of geometry and windows blowing past the max-density limit
+([§14](#14-bugs-encountered-and-how-they-were-fixed)).
+
+### The fix: lattice-aligned cuts + halo, emit-in-core
+
+Two mechanisms guarantee a partitioned fill equals a global fill:
+
+- **Lattice-aligned cuts.** Partition boundaries are placed only on a lattice
+  equal to `align = lcm(density_window, fill_pitch)` (per axis). Because a cut
+  lands on both a window boundary *and* a pitch line, every density window lies
+  entirely inside one leaf, and the fill lattice is continuous across leaves.
+  The working area is snapped **outward** to this lattice (`snap_area`) and tiled
+  exactly, so leaf-local window/pitch indexing coincides with the global one.
+
+- **Halo (guard band) + emit-in-core.** Each leaf is grown by a `margin` (a
+  multiple of `align`, at least `keepout + fill_size`) into a **halo**. The leaf
+  rasterizes *all* geometry touching its halo as **read-only context**, so
+  density and keep-out are exact right up to the core edge — but it only
+  **emits** fill inside its **core**. Cores are disjoint and tile the die.
+
+```mermaid
+flowchart LR
+    subgraph Halo
+      direction TB
+      subgraph Core["core (emit fill here)"]
+        x[" "]
+      end
+    end
+    N["neighbor geometry (read-only context in halo)"] -.-> Halo
+```
+
+### Why merge is trivial
+
+Because cores are **disjoint** and every fill shape is emitted fully inside its
+core, the merge is a **plain concatenation** of each leaf's fill — no
+de-duplication, no boundary stitching. Fill-to-fill spacing across a boundary is
+automatic because both leaves place on the *same global pitch lattice*.
+
+### Verification
+
+`test_partition_equals_global` runs the engine twice on the same layout: once
+finely partitioned (`max_leaf=1`) and once as a single partition
+(`max_leaf=100000`), and asserts identical fill-shape counts and identical
+post-fill min/mean density. `test_partition_cover` independently checks that the
+leaf cores are lattice-aligned, mutually disjoint, and **exactly tile** the die
+(by summed area), and that each halo contains its core.
+
+---
+
+## 10. Gradient-aware, DRC-driven iterative fill
+
+Filling every window just to the minimum leaves steep steps next to naturally
+dense regions (e.g. an SRAM macro), violating the **gradient** rule. A robust
+fill must raise sparse windows *toward* their dense neighbors. This is done with
+a per-window **target map** updated by DRC feedback (in `run_fill`):
+
+1. Initialize every window's target to `min_density + fill_headroom`.
+2. **Fill** with the current target map (partitioned, parallel).
+3. **Measure** the achieved density map.
+4. **Update targets** for each window:
+   - if below `min`, keep at least the base target;
+   - for each neighbor `nb`, if `density(nb) − density(here) > max_gradient`,
+     raise this window's target to `density(nb) − max_gradient + grad_headroom`;
+   - clamp to `max_density`.
+5. If any target increased, **refill**; else stop.
+
+Targets increase **monotonically** and are bounded by `max_density`, so the loop
+converges (default cap `max_iterations = 3`, enough in practice).
+
+The small `grad_headroom` (default 0.02) exists because placement stops *just
+below* a target, leaving the achieved density one fill-quantum short; without the
+overshoot, a window can sit a hair under the gradient limit (e.g. measured
+`0.1501 > 0.15`). Adding the headroom makes the achieved density clear the limit.
+
+This directly implements the "iterative DRC" step: **fill → check → raise targets
+where DRC fails → refill**. Where a violation is *physically unfixable* (a very
+dense macro edge with no room to fill the neighbor high enough), the loop
+converges with a residual that DRC honestly reports.
+
+---
+
+## 11. Design-rule checking
+
+`run_drc` (`drc.{hpp,cpp}`) checks the merged result over the **whole die** —
+independently of how it was partitioned, which doubles as a cross-check on the
+merge. It reports five `ViolationType`s:
+
+| Check | How |
+|-------|-----|
+| `MinDensity` | window density `< min` |
+| `MaxDensity` | window density `> max` |
+| `Gradient` | `|density − neighbor|` (right/down) `> max_gradient` |
+| `Spacing` | any fill cell that lies inside the keep-out halo |
+| `MinArea` | fill-shape area `< min_area` (shape-level static check) |
+
+The spacing check recomputes the keep-out from the *existing* geometry and
+intersects it with the fill footprints; it must be zero by construction, so it is
+a strong safety net. `format_report` renders a per-layer table (partitions,
+existing/fill counts, iterations, before→after density, DRC count, time) plus a
+detailed violation list.
+
+---
+
+## 12. Parallelism and the GPU-ready backend
+
+The four compute kernels — `build_sat`, `compute_density`, `compute_keepout`,
+`place_fill` — are declared on the `FillBackend` interface. `CpuBackend`
+(`backend_cpu.cpp`) implements them by delegating to the reference kernels, which
+are **OpenMP-parallel** where the work is independent:
+
+- `place_fill` parallelizes over tiles (disjoint writes), and
+- the engine parallelizes over **partitions** (`#pragma omp parallel for` in
+  `partitioned_fill`), which is the primary source of speedup.
+
+Every kernel was chosen because it maps naturally onto a GPU:
+
+- `build_sat` → parallel prefix sums,
+- `compute_density`/`compute_keepout` → constant-time box lookups / a box filter,
+- `place_fill` → per-window independent stamping.
+
+A CUDA backend can therefore be added as a second `FillBackend` implementation
+and selected via `EngineConfig::prefer_gpu` with **no change** to the engine.
+That work is deliberately deferred; `make_cuda_backend()` currently returns
+`nullptr`.
+
+---
+
+## 13. The synthetic GPU-block test design
+
+Testing needs realistic input. `make_gpu_block` (`tools/make_gpu_block.cpp`)
+generates a structurally realistic — but not real — GPU block:
+
+- an `N×N` grid of **SM (streaming-multiprocessor) tiles**, each containing two
+  **SRAM macros** (register file + shared memory, dense bitcell-like stripes on
+  `OD/PO/M1..M3`) and a **standard-cell logic** region (cell rows on `OD/PO`,
+  local routing on `M1..M3`, density varying per tile via a deterministic PRNG);
+- **routing channels** between tiles carrying bus routing on `M4..M8`;
+- **global clock/spine** routes on `M9/M10`;
+- a regular, sparse **power grid** (wide straps) on `M11..M14`.
+
+The effect is that every layer has a *different* density profile — dense on the
+lower layers, nearly empty on the upper layers except for power straps — so fill
+does meaningful, layer-specific work everywhere. A simpler generator
+(`make_dummy_gds`) produces a controlled per-window density pattern used by the
+tests.
+
+---
+
+## 14. Bugs encountered and how they were fixed
+
+The project was built and validated incrementally; the important defects and
+their fixes are worth recording.
+
+### 14.1 Lattice-alignment break at partition boundaries (critical)
+
+**Symptom:** after adding partitioning, the demo produced many `MaxDensity`
+(windows at 0.98 vs a 0.80 limit) and `Spacing` violations — fill was landing on
+top of geometry.
+
+**Root cause:** each leaf's core was clamped to the *raw geometry bounding box*,
+whose edges are **not** on the window/pitch lattice. That offset the leaf-local
+density windows from the global windows by a sub-window amount, so `place_fill`'s
+per-window accounting used one set of windows while DRC used another. In a leaf
+where a window appeared nearly empty (because it was misaligned), fill was piled
+in on top of real geometry.
+
+**Fix:** snap the working area **outward to the `align = lcm(window, pitch)`
+lattice** (`snap_area`) and tile it exactly; anchor partitions and both the
+per-leaf and global grids at the same lattice origin. After this, leaf-local
+windows coincide with global windows and the caps apply to the right windows.
+Verified by `test_partition_equals_global`.
+
+### 14.2 Rule inconsistency: pitch too coarse to reach min density
+
+**Symptom:** `M14` (and other upper metals) could not reach the 0.30 min density.
+
+**Root cause:** the initial rules used `pitch = 2 × fill`, which caps fill
+coverage at `(fill/pitch)² = 25 %` — physically below the 30 % minimum.
+
+**Fix:** set `pitch = 1.25 × fill` (and a divisor of the window), giving ~64 %
+achievable coverage with headroom, while keeping the lattice coherent.
+
+### 14.3 Gradient handling
+
+**Symptom:** ~200 `Gradient` violations because filling only to `min` left steep
+steps next to dense windows.
+
+**Fix:** the per-window, DRC-driven target loop of
+[§10](#10-gradient-aware-drc-driven-iterative-fill), plus a small `grad_headroom`
+to beat placement discretization (`0.1501 > 0.15` boundary cases).
+
+### 14.4 `make_grid` over-allocation
+
+An early `make_grid` added a `+1` guard cell, which created a thin extra window at
+the far edge and could produce spurious edge-window density readings. Because the
+area is now snapped outward to the lattice, exact `ceil` sizing covers everything
+and the tiling divides evenly.
+
+### 14.5 Macro-edge gradient in the GPU block
+
+**Symptom:** the first GPU block left a few residual `Gradient` violations at
+SRAM macro edges — a very dense macro window next to logic that fill physically
+cannot raise enough to smooth.
+
+**Resolution:** this is a *real* limitation of any fill flow, not a tool bug. For
+a clean demonstration the synthetic macro densities (test-data parameters only)
+were tapered slightly so the steps are fillable. On real designs such residuals
+are expected and are honestly reported by DRC.
+
+---
+
+## 15. Results
+
+All 77 unit/end-to-end checks pass (`make test`), in both the OpenMP and serial
+(`OPENMP=0`) builds.
+
+**Simple synthetic design** (`make demo`, 100 × 97 µm, 16 layers):
+400 existing polygons → **186,236 fill shapes, 0 DRC violations**. OpenMP across
+partitions runs in ~0.87 s vs ~1.99 s serial (~**2.3×**).
+
+**GPU block** (`make gpu-demo`, 180 × 180 µm, 9 SM tiles):
+15,885 existing polygons → **~1.39 M fill shapes, 0 DRC violations**, 52
+partitions per layer, ~4 s. Every layer is brought in-band: metals reach the 0.30
+min (after-min ≥ 0.32) without exceeding max or violating gradient/spacing.
+
+The `render_layer` tool visualizes any layer (existing = navy, fill = orange,
+partition cores = red); the power-grid layers (`M11..M14`) make the keep-out
+halos around the wide straps clearly visible.
+
+---
+
+## 16. Limitations and future work
+
+- **GPU offload** — the whole point of the backend abstraction. A CUDA
+  implementation of the four kernels is the natural next step.
+- **Via fill** between metal layers is not implemented.
+- **Multi-patterning coloring** (LELE/SADP mask assignment) for advanced nodes.
+- **Grounded/tied fill** (vs floating) and coupling-capacitance-aware fill for
+  timing-critical nets.
+- **Non-rectilinear fill / staggered patterns** and technology-specific fill
+  cells; the current fill is a rectangular lattice.
+- **Hierarchy** — the reader flattens to a single cell; real flows fill
+  hierarchically and reuse fill across instances.
+- **Real PDK rules** — the layer map and rules in `layermap.cpp` are illustrative
+  placeholders, not a real technology.
+- **Fill distribution** — placement fills a window greedily; a more uniform
+  spatial distribution within a window can further reduce local gradients.
+
+---
+
+## 17. File-by-file reference
+
+| File | Responsibility |
+|------|----------------|
+| `include/metalfill/geometry.hpp` | points, bbox, polygons, `make_rect` (all in dbu) |
+| `gdsii.{hpp,cpp}` | GDSII reader/writer, 8-byte real codec |
+| `raster.{hpp,cpp}` | occupancy grid, scanline rasterization, `stamp_rect` |
+| `density.{hpp,cpp}` | summed-area table, windowed density map |
+| `fill.{hpp,cpp}` | keep-out dilation, per-window fill placement |
+| `partition.{hpp,cpp}` | recursive quadtree partitioner, `gcd`/`lcm` |
+| `drc.{hpp,cpp}` | density / gradient / spacing / min-area checks |
+| `rules.hpp`, `layermap.{hpp,cpp}` | fill-rule struct + default layer map |
+| `backend.hpp`, `backend_cpu.cpp` | compute backend interface + CPU/OpenMP impl |
+| `engine.{hpp,cpp}` | orchestration: geom, snap, partition, iterate, merge, report |
+| `tools/make_dummy_gds.cpp` | controlled synthetic layout for tests |
+| `tools/make_gpu_block.cpp` | realistic synthetic GPU block |
+| `tools/run_fill.cpp` | CLI: GDS in → filled GDS + report |
+| `tools/render_layer.cpp` | render a layer (existing/fill/partitions) to PPM |
+| `tests/test_main.cpp` | 77 unit + end-to-end checks |
+
+### Key tunables (`EngineConfig`, `FillRule`)
+
+- `max_iterations`, `max_leaf` (partition leaf size), `fill_headroom`,
+  `grad_headroom`, `prefer_gpu`.
+- Per layer: `window_um`, `step_um`, `min/max_density`, `max_gradient`,
+  `fill_w/h_um`, `fill_pitch_x/y_um`, `keepout_um`, `min_area_um2`, `is_beol`.
+
+---
+
+*This engine is a compact, verifiable reference implementation of the density →
+partition → fill → iterative-DRC flow. Its structure — lattice-aligned recursive
+partitioning with halo context and an emit-in-core merge, behind a GPU-ready
+backend — is exactly what is needed to scale BEOL fill to large dies and, later,
+to a GPU.*
diff --git a/gpu_metal_fill/docs/metal_fill_whitepaper.pdf b/gpu_metal_fill/docs/metal_fill_whitepaper.pdf
new file mode 100644
index 0000000..4744f53
Binary files /dev/null and b/gpu_metal_fill/docs/metal_fill_whitepaper.pdf differ
diff --git a/gpu_metal_fill/include/metalfill/backend.hpp b/gpu_metal_fill/include/metalfill/backend.hpp
new file mode 100644
index 0000000..693b9ac
--- /dev/null
+++ b/gpu_metal_fill/include/metalfill/backend.hpp
@@ -0,0 +1,42 @@
+#pragma once
+#include <memory>
+#include <string>
+#include "metalfill/density.hpp"
+#include "metalfill/fill.hpp"
+#include "metalfill/raster.hpp"
+
+namespace mf {
+
+// The compute-heavy stages of metal fill are isolated behind this interface so
+// they can be executed on the CPU (OpenMP) or offloaded to the GPU (CUDA).
+//
+// Every stage here is data-parallel and maps naturally onto the GPU:
+//   * build_sat      -> parallel prefix sums
+//   * compute_density-> box lookups over the SAT
+//   * compute_keepout-> morphological dilation (box filter + threshold)
+//   * place_fill     -> per-window independent stamping
+class FillBackend {
+public:
+    virtual ~FillBackend() = default;
+    virtual std::string name() const = 0;
+
+    virtual SummedAreaTable build_sat(const Grid& grid) = 0;
+    virtual DensityMap compute_density(const SummedAreaTable& sat, int win_cells,
+                                       int step_cells) = 0;
+    virtual Grid compute_keepout(const Grid& occ, int keepout_cells) = 0;
+    virtual std::vector<FillShape> place_fill(const Grid& occ, const Grid& blocked,
+                                              const FillParams& params,
+                                              Grid& fill_occ) = 0;
+};
+
+// The default CPU backend (parallelized with OpenMP when available).
+std::unique_ptr<FillBackend> make_cpu_backend();
+
+// Returns a CUDA backend when the library is built with USE_CUDA and a device
+// is present; otherwise returns nullptr.
+std::unique_ptr<FillBackend> make_cuda_backend();
+
+// Selects the CUDA backend if requested and available, else the CPU backend.
+std::unique_ptr<FillBackend> make_backend(bool prefer_gpu);
+
+}  // namespace mf
diff --git a/gpu_metal_fill/include/metalfill/density.hpp b/gpu_metal_fill/include/metalfill/density.hpp
new file mode 100644
index 0000000..b578c12
--- /dev/null
+++ b/gpu_metal_fill/include/metalfill/density.hpp
@@ -0,0 +1,46 @@
+#pragma once
+#include <cstdint>
+#include <vector>
+#include "metalfill/raster.hpp"
+
+namespace mf {
+
+// Summed-area table (integral image) over a boolean occupancy grid. Allows O(1)
+// area queries for any axis-aligned window and is the core primitive that maps
+// cleanly onto a GPU (prefix sums + box lookups).
+struct SummedAreaTable {
+    int nx = 0;
+    int ny = 0;
+    std::vector<int64_t> sat;  // (ny+1) x (nx+1), row-major
+
+    // Occupied-cell count in the half-open cell range [x0,x1) x [y0,y1).
+    int64_t area(int x0, int y0, int x1, int y1) const;
+};
+
+SummedAreaTable build_sat(const Grid& grid);
+
+// A single density evaluation window.
+struct Window {
+    int ix = 0, iy = 0;      // tile index
+    int x0 = 0, y0 = 0;      // cell range [x0,x1) x [y0,y1)
+    int x1 = 0, y1 = 0;
+    double density = 0.0;     // occupied_area / window_area
+};
+
+// A grid of density windows produced by sliding a `win_cells` window with
+// `step_cells` stride over the occupancy grid.
+struct DensityMap {
+    int ntiles_x = 0;
+    int ntiles_y = 0;
+    std::vector<Window> tiles;  // row-major (iy*ntiles_x + ix)
+
+    const Window& at(int ix, int iy) const { return tiles[static_cast<size_t>(iy) * ntiles_x + ix]; }
+    double min_density() const;
+    double max_density() const;
+    double mean_density() const;
+};
+
+// Computes the density map from a SAT. win_cells/step_cells are in grid cells.
+DensityMap compute_density(const SummedAreaTable& sat, int win_cells, int step_cells);
+
+}  // namespace mf
diff --git a/gpu_metal_fill/include/metalfill/drc.hpp b/gpu_metal_fill/include/metalfill/drc.hpp
new file mode 100644
index 0000000..e647add
--- /dev/null
+++ b/gpu_metal_fill/include/metalfill/drc.hpp
@@ -0,0 +1,41 @@
+#pragma once
+#include <string>
+#include <vector>
+#include "metalfill/density.hpp"
+#include "metalfill/raster.hpp"
+#include "metalfill/rules.hpp"
+
+namespace mf {
+
+enum class ViolationType {
+    MinDensity,
+    MaxDensity,
+    Gradient,
+    Spacing,
+    MinArea,
+};
+
+struct Violation {
+    ViolationType type;
+    int ix = 0, iy = 0;      // offending tile (or -1 for shape-level checks)
+    double value = 0.0;      // measured value
+    double limit = 0.0;      // rule limit
+    std::string message;
+};
+
+struct DrcReport {
+    std::vector<Violation> violations;
+    bool clean() const { return violations.empty(); }
+    int count(ViolationType t) const;
+};
+
+// Runs density (min/max), gradient, keep-out spacing and min-area checks over
+// the combined (existing + fill) layout for one layer.
+//   total_occ : existing geometry OR fill footprints (post-fill occupancy)
+//   blocked   : keep-out mask (existing geometry dilated by keepout)
+//   fill_occ  : fill footprints only (for spacing check)
+DrcReport run_drc(const SummedAreaTable& total_sat, const Grid& blocked,
+                  const Grid& fill_occ, int win_cells, int step_cells,
+                  const FillRule& rule, double dbu_per_um, dbu cell);
+
+}  // namespace mf
diff --git a/gpu_metal_fill/include/metalfill/engine.hpp b/gpu_metal_fill/include/metalfill/engine.hpp
new file mode 100644
index 0000000..dded791
--- /dev/null
+++ b/gpu_metal_fill/include/metalfill/engine.hpp
@@ -0,0 +1,53 @@
+#pragma once
+#include <string>
+#include <vector>
+#include "metalfill/backend.hpp"
+#include "metalfill/drc.hpp"
+#include "metalfill/gdsii.hpp"
+#include "metalfill/rules.hpp"
+
+namespace mf {
+
+// Per-layer result of the fill run.
+struct LayerResult {
+    std::string name;
+    int layer = 0;
+    bool is_beol = false;
+    int existing_polys = 0;
+    int fill_shapes = 0;
+    int partitions = 0;
+    double density_before_min = 0, density_before_mean = 0, density_before_max = 0;
+    double density_after_min = 0, density_after_mean = 0, density_after_max = 0;
+    int iterations = 0;
+    DrcReport drc;
+    double seconds = 0.0;
+};
+
+struct EngineConfig {
+    // Raster resolution in database units. Smaller = more accurate + more memory.
+    dbu grid_cell_dbu = 0;   // 0 => auto (derived from smallest fill pitch)
+    int max_iterations = 3;  // iterative fill/DRC passes per layer
+    int max_leaf = 2;        // partition leaf size (align-units per side)
+    bool prefer_gpu = false; // offload BEOL to CUDA backend if available
+    double fill_headroom = 0.02;  // aim slightly above min_density
+    double grad_headroom = 0.02;  // overshoot gradient target to beat discretization
+};
+
+struct FillSummary {
+    std::string backend;
+    std::vector<LayerResult> layers;
+    double total_seconds = 0.0;
+    int total_violations() const;
+    int total_fill_shapes() const;
+};
+
+// Runs FEOL + BEOL fill over `layout` using `rules`. Fill polygons are appended
+// to `layout` (tagged with each rule's fill_datatype). BEOL layers are processed
+// through the (optionally GPU) backend.
+FillSummary run_fill(Layout& layout, const std::vector<FillRule>& rules,
+                     const EngineConfig& cfg);
+
+// Renders a human-readable text report.
+std::string format_report(const FillSummary& summary);
+
+}  // namespace mf
diff --git a/gpu_metal_fill/include/metalfill/fill.hpp b/gpu_metal_fill/include/metalfill/fill.hpp
new file mode 100644
index 0000000..49e1612
--- /dev/null
+++ b/gpu_metal_fill/include/metalfill/fill.hpp
@@ -0,0 +1,48 @@
+#pragma once
+#include <vector>
+#include "metalfill/geometry.hpp"
+#include "metalfill/raster.hpp"
+
+namespace mf {
+
+// Result of computing the keep-out (forbidden) mask for fill placement.
+// A cell is blocked if any existing geometry lies within `keepout` of it.
+Grid compute_keepout(const Grid& occ, int keepout_cells);
+
+// A placed fill shape, expressed in grid-cell coordinates (footprint is
+// [ix0, ix0+w) x [iy0, iy0+h)).
+struct FillShape {
+    int ix0 = 0, iy0 = 0;
+    int w = 0, h = 0;
+};
+
+// Parameters controlling fill placement, all expressed in grid cells.
+struct FillParams {
+    int fill_w = 1;
+    int fill_h = 1;
+    int pitch_x = 2;
+    int pitch_y = 2;
+    int win_cells = 1;
+    int step_cells = 1;
+    double target_density = 0.30;  // fill until each window reaches this
+    double max_density = 0.80;     // never exceed this in any window
+
+    // Optional per-window target (gradient-aware fill). When set, a window's
+    // target is looked up from this global map instead of `target_density`.
+    // The map is indexed [gy*gnx + gx] where (gx,gy) is the *global* window
+    // index; a local tile (tx,ty) maps to global (goff_x+tx, goff_y+ty).
+    const std::vector<double>* target_map = nullptr;
+    int gnx = 0;
+    int gny = 0;
+    int goff_x = 0;
+    int goff_y = 0;
+};
+
+// Places dummy fill on a grid so that every window reaches `target_density`
+// where physically possible, while never exceeding `max_density` and never
+// overlapping the keep-out mask. Returns the placed shapes; `fill_occ` is
+// populated with the fill footprints (1 where fill was placed).
+std::vector<FillShape> place_fill(const Grid& occ, const Grid& blocked,
+                                  const FillParams& params, Grid& fill_occ);
+
+}  // namespace mf
diff --git a/gpu_metal_fill/include/metalfill/gdsii.hpp b/gpu_metal_fill/include/metalfill/gdsii.hpp
new file mode 100644
index 0000000..899e7a8
--- /dev/null
+++ b/gpu_metal_fill/include/metalfill/gdsii.hpp
@@ -0,0 +1,37 @@
+#pragma once
+#include <string>
+#include <vector>
+#include "metalfill/geometry.hpp"
+
+namespace mf {
+
+// A minimal in-memory representation of a single-cell GDSII layout.
+struct Layout {
+    std::string lib_name = "METALFILL";
+    std::string cell_name = "TOP";
+    // meters per database unit (UNITS second value). 1e-9 => 1 dbu = 1 nm.
+    double meters_per_dbu = 1e-9;
+    // user units per database unit (UNITS first value). 1e-3 => 1 user unit = 1 um.
+    double user_units_per_dbu = 1e-3;
+
+    std::vector<Polygon> polygons;
+
+    // Database units per micron, derived from UNITS.
+    double dbu_per_um() const { return 1e-6 / meters_per_dbu; }
+
+    BBox bbox() const {
+        BBox b;
+        for (const auto& p : polygons) b.expand(p.bbox());
+        return b;
+    }
+};
+
+// Reads a GDSII file. Only BOUNDARY and BOX records are imported (as polygons);
+// references (SREF/AREF) and paths are ignored. Throws std::runtime_error on
+// malformed input.
+Layout read_gds(const std::string& path);
+
+// Writes a single-cell GDSII file containing all polygons.
+void write_gds(const std::string& path, const Layout& layout);
+
+}  // namespace mf
diff --git a/gpu_metal_fill/include/metalfill/geometry.hpp b/gpu_metal_fill/include/metalfill/geometry.hpp
new file mode 100644
index 0000000..0798bb7
--- /dev/null
+++ b/gpu_metal_fill/include/metalfill/geometry.hpp
@@ -0,0 +1,63 @@
+#pragma once
+#include <cstdint>
+#include <vector>
+#include <limits>
+#include <algorithm>
+
+namespace mf {
+
+// All layout coordinates are stored in integer database units (dbu).
+using dbu = int64_t;
+
+struct Point {
+    dbu x = 0;
+    dbu y = 0;
+};
+
+struct BBox {
+    dbu xmin = std::numeric_limits<dbu>::max();
+    dbu ymin = std::numeric_limits<dbu>::max();
+    dbu xmax = std::numeric_limits<dbu>::min();
+    dbu ymax = std::numeric_limits<dbu>::min();
+
+    bool valid() const { return xmax >= xmin && ymax >= ymin; }
+    dbu width() const { return valid() ? xmax - xmin : 0; }
+    dbu height() const { return valid() ? ymax - ymin : 0; }
+
+    void expand(const Point& p) {
+        xmin = std::min(xmin, p.x);
+        ymin = std::min(ymin, p.y);
+        xmax = std::max(xmax, p.x);
+        ymax = std::max(ymax, p.y);
+    }
+    void expand(const BBox& b) {
+        if (!b.valid()) return;
+        xmin = std::min(xmin, b.xmin);
+        ymin = std::min(ymin, b.ymin);
+        xmax = std::max(xmax, b.xmax);
+        ymax = std::max(ymax, b.ymax);
+    }
+};
+
+// A simple polygon tagged with GDSII (layer, datatype).
+struct Polygon {
+    int layer = 0;
+    int datatype = 0;
+    std::vector<Point> pts;
+
+    BBox bbox() const {
+        BBox b;
+        for (const auto& p : pts) b.expand(p);
+        return b;
+    }
+};
+
+inline Polygon make_rect(int layer, int datatype, dbu x0, dbu y0, dbu x1, dbu y1) {
+    Polygon p;
+    p.layer = layer;
+    p.datatype = datatype;
+    p.pts = {{x0, y0}, {x1, y0}, {x1, y1}, {x0, y1}};
+    return p;
+}
+
+}  // namespace mf
diff --git a/gpu_metal_fill/include/metalfill/layermap.hpp b/gpu_metal_fill/include/metalfill/layermap.hpp
new file mode 100644
index 0000000..e566018
--- /dev/null
+++ b/gpu_metal_fill/include/metalfill/layermap.hpp
@@ -0,0 +1,30 @@
+#pragma once
+#include <vector>
+#include "metalfill/rules.hpp"
+
+namespace mf {
+
+// GDS layer numbering used by this project (illustrative; not a real PDK).
+//
+//   FEOL (base) layers:
+//       OD (active/diffusion) = 1
+//       PO (poly)            = 2
+//   BEOL (metal) layers:
+//       M1 .. M14            = 10 .. 23   (Mn -> layer 9 + n)
+//
+// Fill shapes are written on the same layer number with datatype 10.
+constexpr int kOdLayer = 1;
+constexpr int kPoLayer = 2;
+constexpr int kMetalBaseLayer = 9;  // M1 = 10
+constexpr int kNumMetals = 14;      // M1 .. M14
+constexpr int kFillDatatype = 10;
+
+inline int metal_layer(int n) { return kMetalBaseLayer + n; }  // n in [1..14]
+
+// Builds the default layer map: 2 FEOL layers + M1..M14 BEOL layers.
+//
+// Lower metals use small, tight fill; upper metals use larger fill on a coarser
+// pitch (mirrors real metal stacks where top layers are thicker/wider).
+std::vector<FillRule> default_layermap();
+
+}  // namespace mf
diff --git a/gpu_metal_fill/include/metalfill/partition.hpp b/gpu_metal_fill/include/metalfill/partition.hpp
new file mode 100644
index 0000000..12dc60f
--- /dev/null
+++ b/gpu_metal_fill/include/metalfill/partition.hpp
@@ -0,0 +1,41 @@
+#pragma once
+#include <vector>
+#include "metalfill/geometry.hpp"
+
+namespace mf {
+
+// One leaf of the recursive partitioning of the die.
+//
+//   core : the region this leaf is responsible for emitting fill in. Cores of
+//          different leaves are disjoint and tile the die. Core edges are
+//          aligned to the density-window / fill-pitch lattice.
+//   halo : core expanded by a guard band. Geometry inside the halo is used as
+//          read-only context so density and spacing are correct at the core
+//          edge, but fill is only *emitted* inside the core.
+struct PartitionRect {
+    BBox core;
+    BBox halo;
+    int depth = 0;
+};
+
+struct PartitionConfig {
+    dbu anchor_x = 0;   // lattice origin (die xmin)
+    dbu anchor_y = 0;   // lattice origin (die ymin)
+    dbu align_x = 1;    // atomic cut unit in x (== lcm(window, pitch_x))
+    dbu align_y = 1;    // atomic cut unit in y (== lcm(window, pitch_y))
+    dbu margin = 0;     // halo guard band (multiple of align)
+    int max_leaf = 2;   // stop splitting when leaf <= max_leaf align-units per side
+    int max_depth = 24; // hard recursion cap
+    BBox clip;          // die bounds; cores are clamped to this for emission
+};
+
+// Recursively partitions `area` (the die bounds) into leaves via a quadtree.
+// Splits happen only on the align lattice so windows and fill pitch stay
+// coherent across boundaries.
+std::vector<PartitionRect> partition_recursive(const BBox& area, const PartitionConfig& cfg);
+
+// Integer gcd / lcm helpers on database units.
+dbu gcd_dbu(dbu a, dbu b);
+dbu lcm_dbu(dbu a, dbu b);
+
+}  // namespace mf
diff --git a/gpu_metal_fill/include/metalfill/raster.hpp b/gpu_metal_fill/include/metalfill/raster.hpp
new file mode 100644
index 0000000..aaa4fd1
--- /dev/null
+++ b/gpu_metal_fill/include/metalfill/raster.hpp
@@ -0,0 +1,49 @@
+#pragma once
+#include <cstdint>
+#include <vector>
+#include "metalfill/geometry.hpp"
+
+namespace mf {
+
+// A uniform boolean occupancy grid over a rectangular area of the layout.
+// cell is the edge length of one grid cell in database units. Element (ix, iy)
+// covers [ox + ix*cell, ox + (ix+1)*cell) x [oy + iy*cell, oy + (iy+1)*cell).
+struct Grid {
+    int nx = 0;
+    int ny = 0;
+    dbu cell = 1;
+    dbu ox = 0;
+    dbu oy = 0;
+    std::vector<uint8_t> occ;  // 0/1, row-major (iy*nx + ix)
+
+    void resize(int nx_, int ny_, dbu cell_, dbu ox_, dbu oy_) {
+        nx = nx_;
+        ny = ny_;
+        cell = cell_;
+        ox = ox_;
+        oy = oy_;
+        occ.assign(static_cast<size_t>(nx) * ny, 0);
+    }
+    inline size_t idx(int ix, int iy) const { return static_cast<size_t>(iy) * nx + ix; }
+    inline uint8_t at(int ix, int iy) const { return occ[idx(ix, iy)]; }
+    inline void set(int ix, int iy, uint8_t v) { occ[idx(ix, iy)] = v; }
+    size_t count() const {
+        size_t c = 0;
+        for (auto v : occ) c += v;
+        return c;
+    }
+};
+
+// Allocates a grid covering `area` at resolution `cell` (dbu). The area is
+// snapped so the origin lands on a cell boundary.
+Grid make_grid(const BBox& area, dbu cell);
+
+// Rasterizes the given polygons (whose (layer,datatype) match) into `grid` by
+// setting occupied cells to 1. Uses even-odd scanline fill; robust for
+// arbitrary simple polygons, exact for axis-aligned rectangles.
+void rasterize(Grid& grid, const std::vector<Polygon>& polys, int layer, int datatype);
+
+// Stamps a single axis-aligned rectangle (in dbu) into the grid.
+void stamp_rect(Grid& grid, dbu x0, dbu y0, dbu x1, dbu y1, uint8_t value = 1);
+
+}  // namespace mf
diff --git a/gpu_metal_fill/include/metalfill/rules.hpp b/gpu_metal_fill/include/metalfill/rules.hpp
new file mode 100644
index 0000000..f6b7e76
--- /dev/null
+++ b/gpu_metal_fill/include/metalfill/rules.hpp
@@ -0,0 +1,45 @@
+#pragma once
+#include <string>
+#include <vector>
+
+namespace mf {
+
+// Fill rules for a single mask layer. Distances are given in microns and are
+// converted to database units at runtime using the layout UNITS.
+//
+// A layer is identified by its GDS (layer, datatype). Fill geometry is emitted
+// on the same GDS layer number but tagged with `fill_datatype` so that dummy
+// fill can be distinguished (and stripped) later.
+struct FillRule {
+    std::string name;        // human readable, e.g. "M1", "OD"
+    int layer = 0;           // GDS layer of the existing (drawn) geometry
+    int datatype = 0;        // GDS datatype of the existing geometry
+    int fill_datatype = 10;  // GDS datatype used to tag emitted fill shapes
+
+    // Density window (CMP planarity is evaluated over a sliding window).
+    double window_um = 20.0;  // square window edge length
+    double step_um = 20.0;    // window step (== window_um => non-overlapping tiles)
+
+    // Density targets. CMP needs the metal density to stay inside [min, max]
+    // and the window-to-window gradient to stay bounded.
+    double min_density = 0.30;
+    double max_density = 0.80;
+    double max_gradient = 0.20;  // max |density(tile) - density(neighbor)|; <=0 disables
+
+    // Fill shape geometry.
+    double fill_w_um = 0.4;
+    double fill_h_um = 0.4;
+    double fill_pitch_x_um = 0.8;
+    double fill_pitch_y_um = 0.8;
+
+    // Design-rule spacing between fill and any existing geometry (keep-out halo).
+    double keepout_um = 0.15;
+
+    // Minimum legal fill-shape area (sanity DRC on the fill shape itself).
+    double min_area_um2 = 0.04;
+
+    // BEOL (metal) layers are flagged for the GPU-offloaded processing path.
+    bool is_beol = false;
+};
+
+}  // namespace mf
diff --git a/gpu_metal_fill/src/backend_cpu.cpp b/gpu_metal_fill/src/backend_cpu.cpp
new file mode 100644
index 0000000..c88f910
--- /dev/null
+++ b/gpu_metal_fill/src/backend_cpu.cpp
@@ -0,0 +1,58 @@
+#include "metalfill/backend.hpp"
+
+#ifdef _OPENMP
+#include <omp.h>
+#include <string>
+#endif
+
+namespace mf {
+
+namespace {
+
+// CPU implementation of the fill backend. The heavy stages delegate to the
+// reference kernels in density.cpp / fill.cpp, which are OpenMP-parallel where
+// the work is independent (keep-out dilation and per-tile fill placement).
+class CpuBackend : public FillBackend {
+public:
+    std::string name() const override {
+#ifdef _OPENMP
+        return "cpu-openmp(" + std::to_string(omp_get_max_threads()) + " threads)";
+#else
+        return "cpu-serial";
+#endif
+    }
+
+    SummedAreaTable build_sat(const Grid& grid) override { return mf::build_sat(grid); }
+
+    DensityMap compute_density(const SummedAreaTable& sat, int win_cells,
+                               int step_cells) override {
+        return mf::compute_density(sat, win_cells, step_cells);
+    }
+
+    Grid compute_keepout(const Grid& occ, int keepout_cells) override {
+        return mf::compute_keepout(occ, keepout_cells);
+    }
+
+    std::vector<FillShape> place_fill(const Grid& occ, const Grid& blocked,
+                                      const FillParams& params, Grid& fill_occ) override {
+        return mf::place_fill(occ, blocked, params, fill_occ);
+    }
+};
+
+}  // namespace
+
+std::unique_ptr<FillBackend> make_cpu_backend() { return std::make_unique<CpuBackend>(); }
+
+std::unique_ptr<FillBackend> make_backend(bool prefer_gpu) {
+    if (prefer_gpu) {
+        if (auto gpu = make_cuda_backend()) return gpu;
+    }
+    return make_cpu_backend();
+}
+
+#ifndef USE_CUDA
+// When built without CUDA support the GPU backend is simply unavailable.
+std::unique_ptr<FillBackend> make_cuda_backend() { return nullptr; }
+#endif
+
+}  // namespace mf
diff --git a/gpu_metal_fill/src/density.cpp b/gpu_metal_fill/src/density.cpp
new file mode 100644
index 0000000..54c6483
--- /dev/null
+++ b/gpu_metal_fill/src/density.cpp
@@ -0,0 +1,82 @@
+#include "metalfill/density.hpp"
+
+#include <algorithm>
+#include <limits>
+
+namespace mf {
+
+int64_t SummedAreaTable::area(int x0, int y0, int x1, int y1) const {
+    x0 = std::max(x0, 0);
+    y0 = std::max(y0, 0);
+    x1 = std::min(x1, nx);
+    y1 = std::min(y1, ny);
+    if (x1 <= x0 || y1 <= y0) return 0;
+    const int w = nx + 1;
+    auto S = [&](int x, int y) -> int64_t { return sat[static_cast<size_t>(y) * w + x]; };
+    return S(x1, y1) - S(x0, y1) - S(x1, y0) + S(x0, y0);
+}
+
+SummedAreaTable build_sat(const Grid& grid) {
+    SummedAreaTable t;
+    t.nx = grid.nx;
+    t.ny = grid.ny;
+    const int w = grid.nx + 1;
+    const int h = grid.ny + 1;
+    t.sat.assign(static_cast<size_t>(w) * h, 0);
+    for (int y = 0; y < grid.ny; ++y) {
+        int64_t row = 0;
+        for (int x = 0; x < grid.nx; ++x) {
+            row += grid.at(x, y);
+            t.sat[static_cast<size_t>(y + 1) * w + (x + 1)] =
+                t.sat[static_cast<size_t>(y) * w + (x + 1)] + row;
+        }
+    }
+    return t;
+}
+
+double DensityMap::min_density() const {
+    double m = std::numeric_limits<double>::max();
+    for (const auto& t : tiles) m = std::min(m, t.density);
+    return tiles.empty() ? 0.0 : m;
+}
+double DensityMap::max_density() const {
+    double m = 0.0;
+    for (const auto& t : tiles) m = std::max(m, t.density);
+    return m;
+}
+double DensityMap::mean_density() const {
+    if (tiles.empty()) return 0.0;
+    double s = 0.0;
+    for (const auto& t : tiles) s += t.density;
+    return s / tiles.size();
+}
+
+DensityMap compute_density(const SummedAreaTable& sat, int win_cells, int step_cells) {
+    DensityMap dm;
+    win_cells = std::max(win_cells, 1);
+    step_cells = std::max(step_cells, 1);
+    // Number of tiles so that the whole grid is covered.
+    auto ntiles = [](int n, int step) { return std::max(1, (n + step - 1) / step); };
+    dm.ntiles_x = ntiles(sat.nx, step_cells);
+    dm.ntiles_y = ntiles(sat.ny, step_cells);
+    dm.tiles.resize(static_cast<size_t>(dm.ntiles_x) * dm.ntiles_y);
+
+    for (int ty = 0; ty < dm.ntiles_y; ++ty) {
+        for (int tx = 0; tx < dm.ntiles_x; ++tx) {
+            Window win;
+            win.ix = tx;
+            win.iy = ty;
+            win.x0 = tx * step_cells;
+            win.y0 = ty * step_cells;
+            win.x1 = std::min(win.x0 + win_cells, sat.nx);
+            win.y1 = std::min(win.y0 + win_cells, sat.ny);
+            int64_t occ = sat.area(win.x0, win.y0, win.x1, win.y1);
+            int64_t cells = static_cast<int64_t>(win.x1 - win.x0) * (win.y1 - win.y0);
+            win.density = cells > 0 ? static_cast<double>(occ) / static_cast<double>(cells) : 0.0;
+            dm.tiles[static_cast<size_t>(ty) * dm.ntiles_x + tx] = win;
+        }
+    }
+    return dm;
+}
+
+}  // namespace mf
diff --git a/gpu_metal_fill/src/drc.cpp b/gpu_metal_fill/src/drc.cpp
new file mode 100644
index 0000000..46da03f
--- /dev/null
+++ b/gpu_metal_fill/src/drc.cpp
@@ -0,0 +1,97 @@
+#include "metalfill/drc.hpp"
+
+#include <cmath>
+#include <sstream>
+
+namespace mf {
+
+int DrcReport::count(ViolationType t) const {
+    int c = 0;
+    for (const auto& v : violations)
+        if (v.type == t) ++c;
+    return c;
+}
+
+DrcReport run_drc(const SummedAreaTable& total_sat, const Grid& blocked, const Grid& fill_occ,
+                  int win_cells, int step_cells, const FillRule& rule, double dbu_per_um,
+                  dbu cell) {
+    DrcReport rep;
+    DensityMap dm = compute_density(total_sat, win_cells, step_cells);
+
+    // Density min / max per window.
+    for (const auto& w : dm.tiles) {
+        if (w.density < rule.min_density - 1e-9) {
+            std::ostringstream os;
+            os << rule.name << " window (" << w.ix << "," << w.iy << ") density "
+               << w.density << " < min " << rule.min_density;
+            rep.violations.push_back({ViolationType::MinDensity, w.ix, w.iy, w.density,
+                                      rule.min_density, os.str()});
+        }
+        if (w.density > rule.max_density + 1e-9) {
+            std::ostringstream os;
+            os << rule.name << " window (" << w.ix << "," << w.iy << ") density "
+               << w.density << " > max " << rule.max_density;
+            rep.violations.push_back({ViolationType::MaxDensity, w.ix, w.iy, w.density,
+                                      rule.max_density, os.str()});
+        }
+    }
+
+    // Density gradient between adjacent windows.
+    if (rule.max_gradient > 0) {
+        for (int ty = 0; ty < dm.ntiles_y; ++ty) {
+            for (int tx = 0; tx < dm.ntiles_x; ++tx) {
+                double d = dm.at(tx, ty).density;
+                if (tx + 1 < dm.ntiles_x) {
+                    double g = std::fabs(d - dm.at(tx + 1, ty).density);
+                    if (g > rule.max_gradient + 1e-9) {
+                        std::ostringstream os;
+                        os << rule.name << " gradient " << g << " > " << rule.max_gradient
+                           << " between (" << tx << "," << ty << ") and (" << tx + 1 << "," << ty
+                           << ")";
+                        rep.violations.push_back(
+                            {ViolationType::Gradient, tx, ty, g, rule.max_gradient, os.str()});
+                    }
+                }
+                if (ty + 1 < dm.ntiles_y) {
+                    double g = std::fabs(d - dm.at(tx, ty + 1).density);
+                    if (g > rule.max_gradient + 1e-9) {
+                        std::ostringstream os;
+                        os << rule.name << " gradient " << g << " > " << rule.max_gradient
+                           << " between (" << tx << "," << ty << ") and (" << tx << "," << ty + 1
+                           << ")";
+                        rep.violations.push_back(
+                            {ViolationType::Gradient, tx, ty, g, rule.max_gradient, os.str()});
+                    }
+                }
+            }
+        }
+    }
+
+    // Fill-to-geometry spacing: fill must never land inside the keep-out halo.
+    if (fill_occ.nx == blocked.nx && fill_occ.ny == blocked.ny) {
+        int spacing_hits = 0;
+        for (size_t i = 0; i < fill_occ.occ.size(); ++i)
+            if (fill_occ.occ[i] && blocked.occ[i]) ++spacing_hits;
+        if (spacing_hits > 0) {
+            std::ostringstream os;
+            os << rule.name << " fill overlaps keep-out halo in " << spacing_hits << " cells";
+            rep.violations.push_back(
+                {ViolationType::Spacing, -1, -1, double(spacing_hits), 0.0, os.str()});
+        }
+    }
+
+    // Fill-shape minimum area (shape-level static check).
+    double fill_area_um2 = rule.fill_w_um * rule.fill_h_um;
+    if (fill_area_um2 < rule.min_area_um2 - 1e-9) {
+        std::ostringstream os;
+        os << rule.name << " fill area " << fill_area_um2 << " um^2 < min " << rule.min_area_um2;
+        rep.violations.push_back(
+            {ViolationType::MinArea, -1, -1, fill_area_um2, rule.min_area_um2, os.str()});
+    }
+
+    (void)dbu_per_um;
+    (void)cell;
+    return rep;
+}
+
+}  // namespace mf
diff --git a/gpu_metal_fill/src/engine.cpp b/gpu_metal_fill/src/engine.cpp
new file mode 100644
index 0000000..d51bd9a
--- /dev/null
+++ b/gpu_metal_fill/src/engine.cpp
@@ -0,0 +1,360 @@
+#include "metalfill/engine.hpp"
+
+#include <algorithm>
+#include <chrono>
+#include <cmath>
+#include <iomanip>
+#include <sstream>
+
+#include "metalfill/partition.hpp"
+
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
+namespace mf {
+namespace {
+
+struct LayerGeom {
+    dbu cell = 1;
+    dbu fw = 1, fh = 1, px = 1, py = 1;
+    dbu win = 1, step = 1, keep = 0;
+    dbu align_x = 1, align_y = 1, margin = 1;
+    int fw_c = 1, fh_c = 1, px_c = 1, py_c = 1;
+    int win_c = 1, step_c = 1, keep_c = 0;
+};
+
+dbu to_dbu(double um, double dppu) { return static_cast<dbu>(std::llround(um * dppu)); }
+
+LayerGeom compute_geom(const FillRule& r, double dppu) {
+    LayerGeom g;
+    g.fw = std::max<dbu>(to_dbu(r.fill_w_um, dppu), 1);
+    g.fh = std::max<dbu>(to_dbu(r.fill_h_um, dppu), 1);
+    g.px = std::max<dbu>(to_dbu(r.fill_pitch_x_um, dppu), g.fw);
+    g.py = std::max<dbu>(to_dbu(r.fill_pitch_y_um, dppu), g.fh);
+    g.win = std::max<dbu>(to_dbu(r.window_um, dppu), 1);
+    g.step = std::max<dbu>(to_dbu(r.step_um, dppu), 1);
+    g.keep = std::max<dbu>(to_dbu(r.keepout_um, dppu), 0);
+
+    // Grid cell divides every relevant length so windows/pitch are exact.
+    dbu c = gcd_dbu(g.fw, g.fh);
+    c = gcd_dbu(c, g.px);
+    c = gcd_dbu(c, g.py);
+    c = gcd_dbu(c, g.win);
+    c = gcd_dbu(c, g.step);
+    g.cell = std::max<dbu>(c, 1);
+
+    g.fw_c = static_cast<int>(g.fw / g.cell);
+    g.fh_c = static_cast<int>(g.fh / g.cell);
+    g.px_c = static_cast<int>(g.px / g.cell);
+    g.py_c = static_cast<int>(g.py / g.cell);
+    g.win_c = static_cast<int>(g.win / g.cell);
+    g.step_c = static_cast<int>(g.step / g.cell);
+    g.keep_c = static_cast<int>((g.keep + g.cell - 1) / g.cell);
+
+    // Partition cuts must align to BOTH the window and the fill pitch so the
+    // per-partition result is identical to a global fill.
+    g.align_x = lcm_dbu(g.win, g.px);
+    g.align_y = lcm_dbu(g.win, g.py);
+    // Halo must cover keep-out + one fill footprint of context, rounded up to a
+    // full align unit so the halo origin stays on the lattice.
+    dbu need = g.keep + std::max(g.fw, g.fh);
+    dbu m_x = ((need + g.align_x - 1) / g.align_x) * g.align_x;
+    dbu m_y = ((need + g.align_y - 1) / g.align_y) * g.align_y;
+    g.margin = std::max(m_x, m_y);
+    return g;
+}
+
+dbu floor_to(dbu v, dbu a) { return (v >= 0 ? (v / a) : -(((-v) + a - 1) / a)) * a; }
+dbu ceil_to(dbu v, dbu a) { return (v >= 0 ? ((v + a - 1) / a) : -((-v) / a)) * a; }
+
+// Snaps a bbox outward to the align lattice in each axis. The result is an exact
+// whole number of align-units, so density windows and the fill pitch stay
+// coherent whether computed globally or per-partition.
+BBox snap_area(const BBox& b, dbu align_x, dbu align_y) {
+    BBox s;
+    s.xmin = floor_to(b.xmin, align_x);
+    s.ymin = floor_to(b.ymin, align_y);
+    s.xmax = ceil_to(b.xmax, align_x);
+    s.ymax = ceil_to(b.ymax, align_y);
+    return s;
+}
+
+std::vector<Polygon> layer_polys(const Layout& layout, const FillRule& r) {
+    std::vector<Polygon> v;
+    for (const auto& p : layout.polygons)
+        if (p.layer == r.layer && p.datatype == r.datatype) v.push_back(p);
+    return v;
+}
+
+bool intersects(const BBox& a, const BBox& b) {
+    return !(a.xmax <= b.xmin || b.xmax <= a.xmin || a.ymax <= b.ymin || b.ymax <= a.ymin);
+}
+
+// Runs partitioned fill for a per-window target map; returns emitted fill
+// rectangles (global dbu) tagged on the fill datatype. Each leaf is processed
+// independently (parallel), reading a halo of context and emitting only inside
+// its core, so the concatenated result equals a global fill.
+std::vector<Polygon> partitioned_fill(const std::vector<Polygon>& polys, const FillRule& rule,
+                                      const LayerGeom& g, const std::vector<PartitionRect>& parts,
+                                      FillBackend& backend, const std::vector<double>& target_map,
+                                      int gnx, int gny, dbu area_xmin, dbu area_ymin, double maxd) {
+    std::vector<std::vector<Polygon>> per(parts.size());
+
+#ifdef _OPENMP
+#pragma omp parallel for schedule(dynamic)
+#endif
+    for (int pi = 0; pi < static_cast<int>(parts.size()); ++pi) {
+        const PartitionRect& part = parts[pi];
+        Grid occ = make_grid(part.halo, g.cell);
+
+        // Rasterize only geometry that touches this leaf's halo.
+        std::vector<Polygon> local;
+        for (const auto& p : polys)
+            if (intersects(p.bbox(), part.halo)) local.push_back(p);
+        rasterize(occ, local, rule.layer, rule.datatype);
+
+        Grid blocked = backend.compute_keepout(occ, g.keep_c);
+
+        FillParams fp;
+        fp.fill_w = g.fw_c;
+        fp.fill_h = g.fh_c;
+        fp.pitch_x = g.px_c;
+        fp.pitch_y = g.py_c;
+        fp.win_cells = g.win_c;
+        fp.step_cells = g.step_c;
+        fp.max_density = maxd;
+        fp.target_map = &target_map;
+        fp.gnx = gnx;
+        fp.gny = gny;
+        fp.goff_x = static_cast<int>((occ.ox - area_xmin) / g.step);
+        fp.goff_y = static_cast<int>((occ.oy - area_ymin) / g.step);
+
+        Grid fill_occ;
+        std::vector<FillShape> shapes = backend.place_fill(occ, blocked, fp, fill_occ);
+
+        auto& out = per[pi];
+        for (const auto& s : shapes) {
+            dbu x0 = occ.ox + static_cast<dbu>(s.ix0) * g.cell;
+            dbu y0 = occ.oy + static_cast<dbu>(s.iy0) * g.cell;
+            dbu x1 = x0 + static_cast<dbu>(s.w) * g.cell;
+            dbu y1 = y0 + static_cast<dbu>(s.h) * g.cell;
+            // Emit only shapes fully inside this leaf's core (disjoint merge).
+            if (x0 < part.core.xmin || y0 < part.core.ymin || x1 > part.core.xmax ||
+                y1 > part.core.ymax)
+                continue;
+            out.push_back(make_rect(rule.layer, rule.fill_datatype, x0, y0, x1, y1));
+        }
+    }
+
+    std::vector<Polygon> fills;
+    for (auto& v : per) fills.insert(fills.end(), v.begin(), v.end());
+    return fills;
+}
+
+// Builds occupancy over the whole die for the given polygons on (layer,dt).
+Grid whole_layer_grid(const std::vector<Polygon>& polys, const BBox& area, const FillRule& r,
+                      dbu cell) {
+    Grid g = make_grid(area, cell);
+    rasterize(g, polys, r.layer, r.datatype);
+    return g;
+}
+
+}  // namespace
+
+FillSummary run_fill(Layout& layout, const std::vector<FillRule>& rules, const EngineConfig& cfg) {
+    FillSummary summary;
+    auto backend = make_backend(cfg.prefer_gpu);
+    summary.backend = backend->name();
+
+    double dppu = layout.dbu_per_um();
+    auto t_all0 = std::chrono::steady_clock::now();
+
+    for (const auto& rule : rules) {
+        auto t0 = std::chrono::steady_clock::now();
+        LayerResult res;
+        res.name = rule.name;
+        res.layer = rule.layer;
+        res.is_beol = rule.is_beol;
+
+        std::vector<Polygon> polys = layer_polys(layout, rule);
+        res.existing_polys = static_cast<int>(polys.size());
+
+        LayerGeom g = compute_geom(rule, dppu);
+        BBox die = layout.bbox();
+        if (!die.valid()) {
+            summary.layers.push_back(res);
+            continue;
+        }
+        // Align-lattice-snapped working area (see snap_area). Both the global
+        // grids and the partition anchor use area.xmin/ymin so every window and
+        // fill position is identical globally and per-partition.
+        BBox area = snap_area(die, g.align_x, g.align_y);
+
+        // ---- density before -------------------------------------------------
+        Grid occ_before = whole_layer_grid(polys, area, rule, g.cell);
+        SummedAreaTable sat_before = backend->build_sat(occ_before);
+        DensityMap dm_before = backend->compute_density(sat_before, g.win_c, g.step_c);
+        res.density_before_min = dm_before.min_density();
+        res.density_before_mean = dm_before.mean_density();
+        res.density_before_max = dm_before.max_density();
+
+        // ---- partitions -----------------------------------------------------
+        PartitionConfig pc;
+        pc.anchor_x = area.xmin;
+        pc.anchor_y = area.ymin;
+        pc.align_x = g.align_x;
+        pc.align_y = g.align_y;
+        pc.margin = g.margin;
+        pc.max_leaf = std::max(cfg.max_leaf, 1);  // align-units per leaf side
+        pc.clip = area;  // cores tile the align-snapped area exactly
+        std::vector<PartitionRect> parts = partition_recursive(area, pc);
+        res.partitions = static_cast<int>(parts.size());
+
+        // ---- gradient-aware iterative fill (DRC-driven) ---------------------
+        // Each window gets its own target. We fill, measure density, then raise
+        // targets for windows that still violate min-density or that sit too far
+        // below a denser neighbor (gradient), and refill. Targets increase
+        // monotonically toward max_density, so the loop converges.
+        double min_d = rule.min_density;
+        double max_d = rule.max_density;
+        double grad = rule.max_gradient;
+        int gnx = dm_before.ntiles_x;
+        int gny = dm_before.ntiles_y;
+        double base_target = std::min(min_d + cfg.fill_headroom, max_d);
+        std::vector<double> target_map(static_cast<size_t>(gnx) * gny, base_target);
+
+        std::vector<Polygon> best_fills;
+        DensityMap dm_after = dm_before;
+        int iter = 0;
+        for (; iter < std::max(cfg.max_iterations, 1); ++iter) {
+            std::vector<Polygon> fills = partitioned_fill(
+                polys, rule, g, parts, *backend, target_map, gnx, gny, area.xmin, area.ymin, max_d);
+
+            Grid total = occ_before;  // copy existing occupancy
+            for (const auto& f : fills) {
+                BBox b = f.bbox();
+                stamp_rect(total, b.xmin, b.ymin, b.xmax, b.ymax, 1);
+            }
+            SummedAreaTable sat_after = backend->build_sat(total);
+            DensityMap dm = backend->compute_density(sat_after, g.win_c, g.step_c);
+
+            best_fills = std::move(fills);
+            dm_after = dm;
+
+            // Update targets from the achieved density (DRC feedback).
+            bool changed = false;
+            const int dx[4] = {1, -1, 0, 0};
+            const int dy[4] = {0, 0, 1, -1};
+            for (int ty = 0; ty < gny; ++ty) {
+                for (int tx = 0; tx < gnx; ++tx) {
+                    size_t idx = static_cast<size_t>(ty) * gnx + tx;
+                    double here = dm.at(tx, ty).density;
+                    double need = target_map[idx];
+                    if (here < min_d - 1e-9) need = std::max(need, base_target);
+                    if (grad > 0) {
+                        for (int k = 0; k < 4; ++k) {
+                            int nx = tx + dx[k], ny = ty + dy[k];
+                            if (nx < 0 || ny < 0 || nx >= gnx || ny >= gny) continue;
+                            double nd = dm.at(nx, ny).density;
+                            // Overshoot by grad_headroom so the achieved density
+                            // (which lands just below target) still clears the
+                            // gradient limit despite discrete fill quanta.
+                            if (nd - here > grad + 1e-9)
+                                need = std::max(need, nd - grad + cfg.grad_headroom);
+                        }
+                    }
+                    need = std::min(need, max_d);
+                    if (need > target_map[idx] + 1e-9) {
+                        target_map[idx] = need;
+                        changed = true;
+                    }
+                }
+            }
+            if (!changed) {
+                ++iter;
+                break;
+            }
+        }
+        res.iterations = iter;
+
+        // ---- append fills + final DRC --------------------------------------
+        res.fill_shapes = static_cast<int>(best_fills.size());
+        Grid fill_occ = make_grid(area, g.cell);
+        for (const auto& f : best_fills) {
+            BBox b = f.bbox();
+            stamp_rect(fill_occ, b.xmin, b.ymin, b.xmax, b.ymax, 1);
+        }
+        Grid total = occ_before;
+        for (const auto& f : best_fills) {
+            BBox b = f.bbox();
+            stamp_rect(total, b.xmin, b.ymin, b.xmax, b.ymax, 1);
+        }
+        SummedAreaTable sat_total = backend->build_sat(total);
+        Grid blocked = backend->compute_keepout(occ_before, g.keep_c);
+        res.drc = run_drc(sat_total, blocked, fill_occ, g.win_c, g.step_c, rule, dppu, g.cell);
+
+        res.density_after_min = dm_after.min_density();
+        res.density_after_mean = dm_after.mean_density();
+        res.density_after_max = dm_after.max_density();
+
+        layout.polygons.insert(layout.polygons.end(), best_fills.begin(), best_fills.end());
+
+        auto t1 = std::chrono::steady_clock::now();
+        res.seconds = std::chrono::duration<double>(t1 - t0).count();
+        summary.layers.push_back(std::move(res));
+    }
+
+    auto t_all1 = std::chrono::steady_clock::now();
+    summary.total_seconds = std::chrono::duration<double>(t_all1 - t_all0).count();
+    return summary;
+}
+
+int FillSummary::total_violations() const {
+    int c = 0;
+    for (const auto& l : layers) c += static_cast<int>(l.drc.violations.size());
+    return c;
+}
+int FillSummary::total_fill_shapes() const {
+    int c = 0;
+    for (const auto& l : layers) c += l.fill_shapes;
+    return c;
+}
+
+std::string format_report(const FillSummary& summary) {
+    std::ostringstream os;
+    os << std::fixed << std::setprecision(3);
+    os << "===== Metal Fill Report =====\n";
+    os << "backend      : " << summary.backend << "\n";
+    os << "total time   : " << summary.total_seconds << " s\n";
+    os << "total fill   : " << summary.total_fill_shapes() << " shapes\n";
+    os << "total DRC    : " << summary.total_violations() << " violations\n\n";
+
+    os << std::left << std::setw(6) << "layer" << std::setw(6) << "type" << std::setw(6) << "part"
+       << std::setw(8) << "exist" << std::setw(8) << "fill" << std::setw(6) << "it"
+       << "  dens(before->after) min/mean/max      drc  time(s)\n";
+    os << std::string(112, '-') << "\n";
+    for (const auto& l : summary.layers) {
+        os << std::left << std::setw(6) << l.name << std::setw(6) << (l.is_beol ? "BEOL" : "FEOL")
+           << std::setw(6) << l.partitions << std::setw(8) << l.existing_polys << std::setw(8)
+           << l.fill_shapes << std::setw(6) << l.iterations << "  ";
+        os << std::setprecision(2) << l.density_before_min << "/" << l.density_before_mean << "/"
+           << l.density_before_max << " -> " << l.density_after_min << "/" << l.density_after_mean
+           << "/" << l.density_after_max;
+        os << "   " << std::setw(4) << static_cast<int>(l.drc.violations.size()) << " "
+           << std::setprecision(3) << l.seconds << "\n";
+    }
+
+    // Detail any violations.
+    bool any = false;
+    for (const auto& l : summary.layers)
+        if (!l.drc.clean()) any = true;
+    if (any) {
+        os << "\n--- DRC details ---\n";
+        for (const auto& l : summary.layers)
+            for (const auto& v : l.drc.violations) os << "  " << v.message << "\n";
+    }
+    return os.str();
+}
+
+}  // namespace mf
diff --git a/gpu_metal_fill/src/fill.cpp b/gpu_metal_fill/src/fill.cpp
new file mode 100644
index 0000000..ccdf8ce
--- /dev/null
+++ b/gpu_metal_fill/src/fill.cpp
@@ -0,0 +1,115 @@
+#include "metalfill/fill.hpp"
+
+#include <algorithm>
+#include <cmath>
+
+#include "metalfill/density.hpp"
+
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
+namespace mf {
+
+Grid compute_keepout(const Grid& occ, int keepout_cells) {
+    Grid blocked;
+    blocked.resize(occ.nx, occ.ny, occ.cell, occ.ox, occ.oy);
+    int r = std::max(keepout_cells, 0);
+    SummedAreaTable sat = build_sat(occ);
+    for (int iy = 0; iy < occ.ny; ++iy) {
+        for (int ix = 0; ix < occ.nx; ++ix) {
+            int x0 = ix - r, y0 = iy - r;
+            int x1 = ix + r + 1, y1 = iy + r + 1;
+            blocked.set(ix, iy, sat.area(x0, y0, x1, y1) > 0 ? 1 : 0);
+        }
+    }
+    return blocked;
+}
+
+std::vector<FillShape> place_fill(const Grid& occ, const Grid& blocked,
+                                  const FillParams& params, Grid& fill_occ) {
+    fill_occ.resize(occ.nx, occ.ny, occ.cell, occ.ox, occ.oy);
+
+    const int fw = std::max(params.fill_w, 1);
+    const int fh = std::max(params.fill_h, 1);
+    const int px = std::max(params.pitch_x, fw);
+    const int py = std::max(params.pitch_y, fh);
+    const int win = std::max(params.win_cells, 1);
+    const int step = std::max(params.step_cells, 1);
+    const double fill_area = static_cast<double>(fw) * fh;
+
+    SummedAreaTable occ_sat = build_sat(occ);
+
+    auto ntiles = [](int n, int s) { return std::max(1, (n + s - 1) / s); };
+    const int ntx = ntiles(occ.nx, step);
+    const int nty = ntiles(occ.ny, step);
+    const int ntiles_total = ntx * nty;
+
+    std::vector<std::vector<FillShape>> per_tile(ntiles_total);
+
+    // Each tile writes only footprints whose origin is inside the tile and whose
+    // extent stays inside the tile, so the writes to fill_occ are disjoint and
+    // the loop is safe to run in parallel across tiles (and, later, on the GPU).
+#ifdef _OPENMP
+#pragma omp parallel for schedule(dynamic)
+#endif
+    for (int t = 0; t < ntiles_total; ++t) {
+        int tx = t % ntx;
+        int ty = t / ntx;
+        int wx0 = tx * step;
+        int wy0 = ty * step;
+        int wx1 = std::min(wx0 + win, occ.nx);
+        int wy1 = std::min(wy0 + win, occ.ny);
+        double win_area = static_cast<double>(wx1 - wx0) * (wy1 - wy0);
+        if (win_area <= 0) continue;
+
+        double occ_cells = static_cast<double>(occ_sat.area(wx0, wy0, wx1, wy1));
+        double placed_cells = 0.0;
+
+        // Per-window target (gradient-aware) if a target map is supplied.
+        double tgt = params.target_density;
+        if (params.target_map) {
+            int gx = params.goff_x + tx;
+            int gy = params.goff_y + ty;
+            if (gx >= 0 && gy >= 0 && gx < params.gnx && gy < params.gny)
+                tgt = (*params.target_map)[static_cast<size_t>(gy) * params.gnx + gx];
+        }
+        double target_cells = tgt * win_area;
+        double max_cells = params.max_density * win_area;
+
+        auto& out = per_tile[t];
+
+        // First candidate origin on the global pitch lattice inside the tile.
+        int cx_start = ((wx0 + px - 1) / px) * px;
+        int cy_start = ((wy0 + py - 1) / py) * py;
+
+        for (int cy = cy_start; cy + fh <= wy1; cy += py) {
+            if (occ_cells + placed_cells + fill_area > target_cells) break;
+            for (int cx = cx_start; cx + fw <= wx1; cx += px) {
+                if (occ_cells + placed_cells + fill_area > target_cells) break;
+                if (occ_cells + placed_cells + fill_area > max_cells) break;
+
+                bool free = true;
+                for (int yy = cy; yy < cy + fh && free; ++yy)
+                    for (int xx = cx; xx < cx + fw; ++xx)
+                        if (blocked.at(xx, yy) || fill_occ.at(xx, yy)) {
+                            free = false;
+                            break;
+                        }
+                if (!free) continue;
+
+                for (int yy = cy; yy < cy + fh; ++yy)
+                    for (int xx = cx; xx < cx + fw; ++xx) fill_occ.set(xx, yy, 1);
+                out.push_back(FillShape{cx, cy, fw, fh});
+                placed_cells += fill_area;
+            }
+        }
+    }
+
+    std::vector<FillShape> shapes;
+    for (auto& v : per_tile)
+        shapes.insert(shapes.end(), v.begin(), v.end());
+    return shapes;
+}
+
+}  // namespace mf
diff --git a/gpu_metal_fill/src/gdsii.cpp b/gpu_metal_fill/src/gdsii.cpp
new file mode 100644
index 0000000..3a03c37
--- /dev/null
+++ b/gpu_metal_fill/src/gdsii.cpp
@@ -0,0 +1,231 @@
+#include "metalfill/gdsii.hpp"
+
+#include <cmath>
+#include <cstdint>
+#include <cstring>
+#include <fstream>
+#include <stdexcept>
+
+namespace mf {
+namespace {
+
+// ---- GDSII record identifiers (rectype << 8 | datatype) --------------------
+constexpr uint8_t RT_HEADER = 0x00, RT_BGNLIB = 0x01, RT_LIBNAME = 0x02;
+constexpr uint8_t RT_UNITS = 0x03, RT_ENDLIB = 0x04, RT_BGNSTR = 0x05;
+constexpr uint8_t RT_STRNAME = 0x06, RT_ENDSTR = 0x07, RT_BOUNDARY = 0x08;
+constexpr uint8_t RT_LAYER = 0x0D, RT_DATATYPE = 0x0E, RT_XY = 0x10;
+constexpr uint8_t RT_ENDEL = 0x11, RT_BOX = 0x2D, RT_BOXTYPE = 0x2E;
+
+constexpr uint8_t DT_NODATA = 0x00, DT_INT2 = 0x02, DT_INT4 = 0x03;
+constexpr uint8_t DT_REAL8 = 0x05, DT_STR = 0x06;
+
+// ---- big-endian helpers ----------------------------------------------------
+void put_u16(std::ofstream& os, uint16_t v) {
+    char b[2] = {char((v >> 8) & 0xFF), char(v & 0xFF)};
+    os.write(b, 2);
+}
+void put_i16(std::ofstream& os, int16_t v) { put_u16(os, static_cast<uint16_t>(v)); }
+void put_i32(std::ofstream& os, int32_t v) {
+    char b[4] = {char((v >> 24) & 0xFF), char((v >> 16) & 0xFF), char((v >> 8) & 0xFF),
+                 char(v & 0xFF)};
+    os.write(b, 4);
+}
+
+uint16_t get_u16(const uint8_t* p) { return (uint16_t(p[0]) << 8) | p[1]; }
+int32_t get_i32(const uint8_t* p) {
+    return int32_t((uint32_t(p[0]) << 24) | (uint32_t(p[1]) << 16) | (uint32_t(p[2]) << 8) |
+                   uint32_t(p[3]));
+}
+
+// GDSII 8-byte real: sign (bit7), 7-bit excess-64 base-16 exponent, 56-bit mantissa.
+void put_real8(std::ofstream& os, double v) {
+    uint8_t out[8] = {0};
+    if (v != 0.0) {
+        bool neg = v < 0;
+        double d = std::fabs(v);
+        int exp = 64;
+        while (d >= 1.0) { d /= 16.0; ++exp; }
+        while (d < 1.0 / 16.0) { d *= 16.0; --exp; }
+        // d now in [1/16, 1); build 56-bit mantissa.
+        uint64_t mant = static_cast<uint64_t>(std::llround(d * static_cast<double>(1ULL << 56)));
+        out[0] = static_cast<uint8_t>((neg ? 0x80 : 0x00) | (exp & 0x7F));
+        for (int i = 7; i >= 1; --i) {
+            out[i] = static_cast<uint8_t>(mant & 0xFF);
+            mant >>= 8;
+        }
+    }
+    os.write(reinterpret_cast<char*>(out), 8);
+}
+
+double get_real8(const uint8_t* p) {
+    bool neg = (p[0] & 0x80) != 0;
+    int exp = (p[0] & 0x7F) - 64;
+    uint64_t mant = 0;
+    for (int i = 1; i < 8; ++i) mant = (mant << 8) | p[i];
+    double d = static_cast<double>(mant) / static_cast<double>(1ULL << 56);
+    d *= std::pow(16.0, exp);
+    return neg ? -d : d;
+}
+
+void write_record(std::ofstream& os, uint8_t rectype, uint8_t datatype, const char* data,
+                  uint16_t data_len) {
+    uint16_t total = static_cast<uint16_t>(4 + data_len);
+    put_u16(os, total);
+    char hdr[2] = {char(rectype), char(datatype)};
+    os.write(hdr, 2);
+    if (data_len) os.write(data, data_len);
+}
+
+void write_str(std::ofstream& os, uint8_t rectype, const std::string& s) {
+    std::string p = s;
+    if (p.size() & 1) p.push_back('\0');  // pad to even length
+    write_record(os, rectype, DT_STR, p.data(), static_cast<uint16_t>(p.size()));
+}
+
+}  // namespace
+
+Layout read_gds(const std::string& path) {
+    std::ifstream is(path, std::ios::binary);
+    if (!is) throw std::runtime_error("cannot open GDS for read: " + path);
+
+    Layout layout;
+    Polygon cur;
+    bool in_elem = false;
+    bool is_box = false;
+
+    while (true) {
+        uint8_t hdr[4];
+        is.read(reinterpret_cast<char*>(hdr), 4);
+        if (is.gcount() == 0) break;  // clean EOF
+        if (is.gcount() != 4) throw std::runtime_error("truncated GDS record header");
+
+        uint16_t total = get_u16(hdr);
+        uint8_t rectype = hdr[2];
+        uint8_t datatype = hdr[3];
+        if (total < 4) throw std::runtime_error("invalid GDS record length");
+        int data_len = total - 4;
+        std::vector<uint8_t> data(data_len);
+        if (data_len) {
+            is.read(reinterpret_cast<char*>(data.data()), data_len);
+            if (is.gcount() != data_len) throw std::runtime_error("truncated GDS record body");
+        }
+
+        switch (rectype) {
+            case RT_UNITS:
+                if (data_len >= 16) {
+                    layout.user_units_per_dbu = get_real8(data.data());
+                    layout.meters_per_dbu = get_real8(data.data() + 8);
+                }
+                break;
+            case RT_LIBNAME:
+                layout.lib_name.assign(reinterpret_cast<char*>(data.data()), data_len);
+                break;
+            case RT_STRNAME:
+                layout.cell_name.assign(reinterpret_cast<char*>(data.data()), data_len);
+                if (!layout.cell_name.empty() && layout.cell_name.back() == '\0')
+                    layout.cell_name.pop_back();
+                break;
+            case RT_BOUNDARY:
+                cur = Polygon{};
+                in_elem = true;
+                is_box = false;
+                break;
+            case RT_BOX:
+                cur = Polygon{};
+                in_elem = true;
+                is_box = true;
+                break;
+            case RT_LAYER:
+                if (in_elem && data_len >= 2) cur.layer = int(int16_t(get_u16(data.data())));
+                break;
+            case RT_DATATYPE:
+            case RT_BOXTYPE:
+                if (in_elem && data_len >= 2) cur.datatype = int(int16_t(get_u16(data.data())));
+                break;
+            case RT_XY:
+                if (in_elem) {
+                    int npairs = data_len / 8;
+                    cur.pts.clear();
+                    cur.pts.reserve(npairs);
+                    for (int i = 0; i < npairs; ++i) {
+                        dbu x = get_i32(data.data() + i * 8);
+                        dbu y = get_i32(data.data() + i * 8 + 4);
+                        cur.pts.push_back({x, y});
+                    }
+                    // Drop the closing duplicate vertex if present.
+                    if (cur.pts.size() > 1 && cur.pts.front().x == cur.pts.back().x &&
+                        cur.pts.front().y == cur.pts.back().y)
+                        cur.pts.pop_back();
+                }
+                break;
+            case RT_ENDEL:
+                if (in_elem && cur.pts.size() >= 3) layout.polygons.push_back(cur);
+                in_elem = false;
+                is_box = false;
+                break;
+            default:
+                break;  // ignore HEADER/BGNLIB/PATH/SREF/etc.
+        }
+        (void)datatype;
+        (void)is_box;
+    }
+    return layout;
+}
+
+void write_gds(const std::string& path, const Layout& layout) {
+    std::ofstream os(path, std::ios::binary);
+    if (!os) throw std::runtime_error("cannot open GDS for write: " + path);
+
+    // HEADER (version 600)
+    char ver[2] = {0x02, 0x58};
+    write_record(os, RT_HEADER, DT_INT2, ver, 2);
+
+    // BGNLIB (24 bytes of timestamps, all zero)
+    char zeros[24] = {0};
+    write_record(os, RT_BGNLIB, DT_INT2, zeros, 24);
+    write_str(os, RT_LIBNAME, layout.lib_name);
+
+    // UNITS: [user units per dbu, meters per dbu]
+    put_u16(os, 20);  // 4 header + 16 payload
+    os.put(char(RT_UNITS));
+    os.put(char(DT_REAL8));
+    put_real8(os, layout.user_units_per_dbu);
+    put_real8(os, layout.meters_per_dbu);
+
+    // BGNSTR
+    write_record(os, RT_BGNSTR, DT_INT2, zeros, 24);
+    write_str(os, RT_STRNAME, layout.cell_name);
+
+    for (const auto& poly : layout.polygons) {
+        if (poly.pts.size() < 3) continue;
+        write_record(os, RT_BOUNDARY, DT_NODATA, nullptr, 0);
+        put_u16(os, 6);  // LAYER record
+        os.put(char(RT_LAYER));
+        os.put(char(DT_INT2));
+        put_i16(os, static_cast<int16_t>(poly.layer));
+        put_u16(os, 6);  // DATATYPE record
+        os.put(char(RT_DATATYPE));
+        os.put(char(DT_INT2));
+        put_i16(os, static_cast<int16_t>(poly.datatype));
+
+        // XY: points + closing vertex.
+        int npts = static_cast<int>(poly.pts.size()) + 1;
+        uint16_t xy_len = static_cast<uint16_t>(npts * 8);
+        put_u16(os, static_cast<uint16_t>(4 + xy_len));
+        os.put(char(RT_XY));
+        os.put(char(DT_INT4));
+        for (const auto& p : poly.pts) {
+            put_i32(os, static_cast<int32_t>(p.x));
+            put_i32(os, static_cast<int32_t>(p.y));
+        }
+        put_i32(os, static_cast<int32_t>(poly.pts.front().x));
+        put_i32(os, static_cast<int32_t>(poly.pts.front().y));
+
+        write_record(os, RT_ENDEL, DT_NODATA, nullptr, 0);
+    }
+
+    write_record(os, RT_ENDSTR, DT_NODATA, nullptr, 0);
+    write_record(os, RT_ENDLIB, DT_NODATA, nullptr, 0);
+}
+
+}  // namespace mf
diff --git a/gpu_metal_fill/src/layermap.cpp b/gpu_metal_fill/src/layermap.cpp
new file mode 100644
index 0000000..1615633
--- /dev/null
+++ b/gpu_metal_fill/src/layermap.cpp
@@ -0,0 +1,91 @@
+#include "metalfill/layermap.hpp"
+
+namespace mf {
+
+std::vector<FillRule> default_layermap() {
+    std::vector<FillRule> rules;
+
+    // ---- FEOL / base layers ------------------------------------------------
+    {
+        FillRule od;
+        od.name = "OD";
+        od.layer = kOdLayer;
+        od.fill_datatype = kFillDatatype;
+        od.window_um = 20.0;
+        od.step_um = 20.0;
+        od.min_density = 0.25;
+        od.max_density = 0.75;
+        od.max_gradient = 0.20;
+        od.fill_w_um = 0.5;
+        od.fill_h_um = 0.5;
+        od.fill_pitch_x_um = 0.625;  // 1.25x fill -> ~64% max coverage
+        od.fill_pitch_y_um = 0.625;
+        od.keepout_um = 0.20;
+        od.min_area_um2 = 0.25;
+        od.is_beol = false;
+        rules.push_back(od);
+
+        FillRule po;
+        po.name = "PO";
+        po.layer = kPoLayer;
+        po.fill_datatype = kFillDatatype;
+        po.window_um = 20.0;
+        po.step_um = 20.0;
+        po.min_density = 0.20;
+        po.max_density = 0.70;
+        po.max_gradient = 0.20;
+        po.fill_w_um = 0.4;
+        po.fill_h_um = 0.4;
+        po.fill_pitch_x_um = 0.5;  // 1.25x fill
+        po.fill_pitch_y_um = 0.5;
+        po.keepout_um = 0.15;
+        po.min_area_um2 = 0.16;
+        po.is_beol = false;
+        rules.push_back(po);
+    }
+
+    // ---- BEOL / metal layers M1..M14 --------------------------------------
+    for (int n = 1; n <= kNumMetals; ++n) {
+        FillRule m;
+        m.name = "M" + std::to_string(n);
+        m.layer = metal_layer(n);
+        m.fill_datatype = kFillDatatype;
+        m.window_um = 20.0;
+        m.step_um = 20.0;
+        m.min_density = 0.30;
+        m.max_density = 0.80;
+        m.max_gradient = 0.15;
+        m.keepout_um = (n <= 6) ? 0.10 : 0.30;
+        m.is_beol = true;
+
+        // Pitch is kept at 1.25x the fill size (and a divisor of the 20um
+        // window) so fill can reach ~64% coverage -- comfortably above the min
+        // density target -- while windows/pitch stay lattice-coherent.
+        if (n <= 6) {
+            // Lower/thin metals: small dense fill.
+            m.fill_w_um = 0.20;
+            m.fill_h_um = 0.20;
+            m.fill_pitch_x_um = 0.25;
+            m.fill_pitch_y_um = 0.25;
+            m.min_area_um2 = 0.04;
+        } else if (n <= 10) {
+            m.fill_w_um = 0.40;
+            m.fill_h_um = 0.40;
+            m.fill_pitch_x_um = 0.50;
+            m.fill_pitch_y_um = 0.50;
+            m.min_area_um2 = 0.16;
+        } else {
+            // Upper/thick metals: large fill on a coarse pitch.
+            m.fill_w_um = 0.80;
+            m.fill_h_um = 0.80;
+            m.fill_pitch_x_um = 1.00;
+            m.fill_pitch_y_um = 1.00;
+            m.min_area_um2 = 0.64;
+        }
+        rules.push_back(m);
+    }
+
+    return rules;
+}
+
+}  // namespace mf
diff --git a/gpu_metal_fill/src/partition.cpp b/gpu_metal_fill/src/partition.cpp
new file mode 100644
index 0000000..b23cb73
--- /dev/null
+++ b/gpu_metal_fill/src/partition.cpp
@@ -0,0 +1,91 @@
+#include "metalfill/partition.hpp"
+
+#include <algorithm>
+
+namespace mf {
+
+dbu gcd_dbu(dbu a, dbu b) {
+    a = a < 0 ? -a : a;
+    b = b < 0 ? -b : b;
+    while (b) {
+        dbu t = a % b;
+        a = b;
+        b = t;
+    }
+    return a == 0 ? 1 : a;
+}
+
+dbu lcm_dbu(dbu a, dbu b) {
+    if (a == 0 || b == 0) return std::max<dbu>(a, b);
+    return (a / gcd_dbu(a, b)) * b;
+}
+
+namespace {
+
+// Recurse in integer align-lattice coordinates: region [ix0,ix1) x [iy0,iy1).
+void recurse(int ix0, int iy0, int ix1, int iy1, int depth, const PartitionConfig& cfg,
+             std::vector<PartitionRect>& out) {
+    int nx = ix1 - ix0;
+    int ny = iy1 - iy0;
+    if (nx <= 0 || ny <= 0) return;
+
+    bool leaf = (nx <= cfg.max_leaf && ny <= cfg.max_leaf) || depth >= cfg.max_depth;
+    if (leaf) {
+        PartitionRect pr;
+        pr.depth = depth;
+        pr.core.xmin = cfg.anchor_x + static_cast<dbu>(ix0) * cfg.align_x;
+        pr.core.xmax = cfg.anchor_x + static_cast<dbu>(ix1) * cfg.align_x;
+        pr.core.ymin = cfg.anchor_y + static_cast<dbu>(iy0) * cfg.align_y;
+        pr.core.ymax = cfg.anchor_y + static_cast<dbu>(iy1) * cfg.align_y;
+        // Clamp emission core to the die bounds.
+        if (cfg.clip.valid()) {
+            pr.core.xmin = std::max(pr.core.xmin, cfg.clip.xmin);
+            pr.core.ymin = std::max(pr.core.ymin, cfg.clip.ymin);
+            pr.core.xmax = std::min(pr.core.xmax, cfg.clip.xmax);
+            pr.core.ymax = std::min(pr.core.ymax, cfg.clip.ymax);
+        }
+        pr.halo.xmin = pr.core.xmin - cfg.margin;
+        pr.halo.ymin = pr.core.ymin - cfg.margin;
+        pr.halo.xmax = pr.core.xmax + cfg.margin;
+        pr.halo.ymax = pr.core.ymax + cfg.margin;
+        if (pr.core.xmax > pr.core.xmin && pr.core.ymax > pr.core.ymin) out.push_back(pr);
+        return;
+    }
+
+    bool split_x = nx > 1;
+    bool split_y = ny > 1;
+    int mx = ix0 + (nx + 1) / 2;
+    int my = iy0 + (ny + 1) / 2;
+
+    if (split_x && split_y) {
+        // Quadtree split into four children.
+        recurse(ix0, iy0, mx, my, depth + 1, cfg, out);
+        recurse(mx, iy0, ix1, my, depth + 1, cfg, out);
+        recurse(ix0, my, mx, iy1, depth + 1, cfg, out);
+        recurse(mx, my, ix1, iy1, depth + 1, cfg, out);
+    } else if (split_x) {
+        recurse(ix0, iy0, mx, iy1, depth + 1, cfg, out);
+        recurse(mx, iy0, ix1, iy1, depth + 1, cfg, out);
+    } else {
+        recurse(ix0, iy0, ix1, my, depth + 1, cfg, out);
+        recurse(ix0, my, ix1, iy1, depth + 1, cfg, out);
+    }
+}
+
+}  // namespace
+
+std::vector<PartitionRect> partition_recursive(const BBox& area, const PartitionConfig& cfg) {
+    std::vector<PartitionRect> out;
+    if (!area.valid() || cfg.align_x <= 0 || cfg.align_y <= 0) return out;
+
+    dbu w = area.xmax - cfg.anchor_x;
+    dbu h = area.ymax - cfg.anchor_y;
+    int nX = static_cast<int>((w + cfg.align_x - 1) / cfg.align_x);
+    int nY = static_cast<int>((h + cfg.align_y - 1) / cfg.align_y);
+    nX = std::max(nX, 1);
+    nY = std::max(nY, 1);
+    recurse(0, 0, nX, nY, 0, cfg, out);
+    return out;
+}
+
+}  // namespace mf
diff --git a/gpu_metal_fill/src/raster.cpp b/gpu_metal_fill/src/raster.cpp
new file mode 100644
index 0000000..ab9a1ef
--- /dev/null
+++ b/gpu_metal_fill/src/raster.cpp
@@ -0,0 +1,80 @@
+#include "metalfill/raster.hpp"
+
+#include <algorithm>
+#include <cmath>
+
+namespace mf {
+
+Grid make_grid(const BBox& area, dbu cell) {
+    Grid g;
+    if (cell <= 0) cell = 1;
+    dbu ox = area.valid() ? area.xmin : 0;
+    dbu oy = area.valid() ? area.ymin : 0;
+    dbu w = area.valid() ? area.width() : 0;
+    dbu h = area.valid() ? area.height() : 0;
+    // Enough cells to fully cover the area (ceil).
+    int nx = static_cast<int>((w + cell - 1) / cell);
+    int ny = static_cast<int>((h + cell - 1) / cell);
+    g.resize(std::max(nx, 1), std::max(ny, 1), cell, ox, oy);
+    return g;
+}
+
+void stamp_rect(Grid& grid, dbu x0, dbu y0, dbu x1, dbu y1, uint8_t value) {
+    if (x1 < x0) std::swap(x0, x1);
+    if (y1 < y0) std::swap(y1, y0);
+    // Cells whose center lies within [x0,x1) x [y0,y1).
+    int ix0 = static_cast<int>(std::floor((static_cast<double>(x0 - grid.ox)) / grid.cell));
+    int iy0 = static_cast<int>(std::floor((static_cast<double>(y0 - grid.oy)) / grid.cell));
+    int ix1 = static_cast<int>(std::ceil((static_cast<double>(x1 - grid.ox)) / grid.cell));
+    int iy1 = static_cast<int>(std::ceil((static_cast<double>(y1 - grid.oy)) / grid.cell));
+    ix0 = std::max(ix0, 0);
+    iy0 = std::max(iy0, 0);
+    ix1 = std::min(ix1, grid.nx);
+    iy1 = std::min(iy1, grid.ny);
+    for (int iy = iy0; iy < iy1; ++iy)
+        for (int ix = ix0; ix < ix1; ++ix) grid.set(ix, iy, value);
+}
+
+void rasterize(Grid& grid, const std::vector<Polygon>& polys, int layer, int datatype) {
+    for (const auto& poly : polys) {
+        if (poly.layer != layer || poly.datatype != datatype) continue;
+        if (poly.pts.size() < 3) continue;
+
+        BBox b = poly.bbox();
+        int row0 = static_cast<int>(std::floor((static_cast<double>(b.ymin - grid.oy)) / grid.cell));
+        int row1 = static_cast<int>(std::ceil((static_cast<double>(b.ymax - grid.oy)) / grid.cell));
+        row0 = std::max(row0, 0);
+        row1 = std::min(row1, grid.ny);
+
+        const size_t n = poly.pts.size();
+        std::vector<double> xs;
+        for (int iy = row0; iy < row1; ++iy) {
+            // Scanline at the vertical center of the cell row.
+            double yc = grid.oy + (iy + 0.5) * grid.cell;
+            xs.clear();
+            for (size_t i = 0; i < n; ++i) {
+                const Point& a = poly.pts[i];
+                const Point& c = poly.pts[(i + 1) % n];
+                double ay = a.y, cy = c.y;
+                if ((ay <= yc && cy > yc) || (cy <= yc && ay > yc)) {
+                    double t = (yc - ay) / (cy - ay);
+                    xs.push_back(a.x + t * (c.x - a.x));
+                }
+            }
+            std::sort(xs.begin(), xs.end());
+            for (size_t k = 0; k + 1 < xs.size(); k += 2) {
+                double xl = xs[k], xr = xs[k + 1];
+                // Cells whose center is inside [xl, xr).
+                int cx0 = static_cast<int>(
+                    std::ceil((xl - grid.ox) / grid.cell - 0.5));
+                int cx1 = static_cast<int>(
+                    std::floor((xr - grid.ox) / grid.cell - 0.5));
+                cx0 = std::max(cx0, 0);
+                cx1 = std::min(cx1, grid.nx - 1);
+                for (int ix = cx0; ix <= cx1; ++ix) grid.set(ix, iy, 1);
+            }
+        }
+    }
+}
+
+}  // namespace mf
diff --git a/gpu_metal_fill/tests/test_main.cpp b/gpu_metal_fill/tests/test_main.cpp
new file mode 100644
index 0000000..2ec65a0
--- /dev/null
+++ b/gpu_metal_fill/tests/test_main.cpp
@@ -0,0 +1,241 @@
+// Lightweight, dependency-free test harness for the metal-fill library.
+#include <cmath>
+#include <cstdio>
+#include <cstdlib>
+#include <string>
+#include <vector>
+
+#include "metalfill/density.hpp"
+#include "metalfill/drc.hpp"
+#include "metalfill/engine.hpp"
+#include "metalfill/fill.hpp"
+#include "metalfill/gdsii.hpp"
+#include "metalfill/layermap.hpp"
+#include "metalfill/partition.hpp"
+#include "metalfill/raster.hpp"
+
+using namespace mf;
+
+static int g_fails = 0;
+static int g_checks = 0;
+#define CHECK(cond, msg)                                                       \
+    do {                                                                       \
+        ++g_checks;                                                            \
+        if (!(cond)) {                                                         \
+            ++g_fails;                                                         \
+            std::printf("  FAIL: %s (%s:%d)\n", msg, __FILE__, __LINE__);      \
+        }                                                                      \
+    } while (0)
+
+static void test_gds_roundtrip() {
+    std::printf("[gds roundtrip]\n");
+    Layout L;
+    L.cell_name = "TOP";
+    L.polygons.push_back(make_rect(10, 0, 0, 0, 5000, 3000));
+    L.polygons.push_back(make_rect(23, 5, -2000, -1000, 4000, 7000));
+    const std::string path = "build/_test_rt.gds";
+    write_gds(path, L);
+    Layout R = read_gds(path);
+    CHECK(R.polygons.size() == 2, "polygon count preserved");
+    CHECK(R.cell_name == "TOP", "cell name preserved");
+    CHECK(std::llround(R.dbu_per_um()) == 1000, "units preserved (1000 dbu/um)");
+    bool found = false;
+    for (const auto& p : R.polygons)
+        if (p.layer == 23 && p.datatype == 5) {
+            BBox b = p.bbox();
+            found = (b.xmin == -2000 && b.ymin == -1000 && b.xmax == 4000 && b.ymax == 7000);
+        }
+    CHECK(found, "layer/datatype/coords preserved");
+}
+
+static void test_sat_density() {
+    std::printf("[sat + density]\n");
+    Grid g;
+    g.resize(4, 4, 100, 0, 0);
+    // Fill a 2x2 block in the corner.
+    for (int y = 0; y < 2; ++y)
+        for (int x = 0; x < 2; ++x) g.set(x, y, 1);
+    SummedAreaTable sat = build_sat(g);
+    CHECK(sat.area(0, 0, 4, 4) == 4, "total occupied = 4");
+    CHECK(sat.area(0, 0, 2, 2) == 4, "block area = 4");
+    CHECK(sat.area(2, 2, 4, 4) == 0, "empty quadrant = 0");
+    DensityMap dm = compute_density(sat, 2, 2);  // four 2x2 windows
+    CHECK(dm.ntiles_x == 2 && dm.ntiles_y == 2, "2x2 tiling");
+    CHECK(std::fabs(dm.at(0, 0).density - 1.0) < 1e-9, "corner window density 1.0");
+    CHECK(std::fabs(dm.at(1, 1).density - 0.0) < 1e-9, "far window density 0.0");
+    CHECK(std::fabs(dm.mean_density() - 0.25) < 1e-9, "mean density 0.25");
+}
+
+static void test_rasterize_rect() {
+    std::printf("[rasterize]\n");
+    Grid g = make_grid(BBox{0, 0, 1000, 1000}, 100);  // 10x10 cells (11x11 with +1)
+    std::vector<Polygon> polys = {make_rect(1, 0, 200, 200, 700, 700)};
+    rasterize(g, polys, 1, 0);
+    // 500x500 dbu rect at 100 dbu/cell -> ~5x5 = 25 cells.
+    CHECK(g.count() == 25, "rasterized rect area = 25 cells");
+}
+
+static void test_keepout() {
+    std::printf("[keepout dilation]\n");
+    Grid g;
+    g.resize(7, 7, 100, 0, 0);
+    g.set(3, 3, 1);  // single occupied cell in the middle
+    Grid b = compute_keepout(g, 1);
+    CHECK(b.at(3, 3) == 1, "center blocked");
+    CHECK(b.at(2, 2) == 1 && b.at(4, 4) == 1, "diagonal neighbors blocked (r=1)");
+    CHECK(b.at(1, 1) == 0, "distance-2 cell not blocked");
+    Grid b0 = compute_keepout(g, 0);
+    CHECK(b0.count() == 1, "r=0 keepout equals occupancy");
+}
+
+static void test_place_fill() {
+    std::printf("[place fill]\n");
+    Grid occ;
+    occ.resize(40, 40, 100, 0, 0);  // empty layer
+    Grid blocked = compute_keepout(occ, 0);
+    FillParams fp;
+    fp.fill_w = 1;
+    fp.fill_h = 1;
+    fp.pitch_x = 2;
+    fp.pitch_y = 2;
+    fp.win_cells = 40;
+    fp.step_cells = 40;
+    fp.target_density = 0.20;
+    fp.max_density = 0.80;
+    Grid fill_occ;
+    auto shapes = place_fill(occ, blocked, fp, fill_occ);
+    double dens = double(fill_occ.count()) / (40.0 * 40.0);
+    CHECK(!shapes.empty(), "fill placed");
+    CHECK(dens >= 0.20 - 1e-6, "reached target density");
+    CHECK(dens <= 0.80 + 1e-6, "did not exceed max density");
+
+    // With keepout covering everything, no fill can be placed.
+    Grid occ2;
+    occ2.resize(10, 10, 100, 0, 0);
+    for (auto& v : occ2.occ) v = 1;
+    Grid blocked2 = compute_keepout(occ2, 0);
+    Grid fo2;
+    auto s2 = place_fill(occ2, blocked2, fp, fo2);
+    CHECK(s2.empty(), "no fill where fully blocked");
+}
+
+static void test_partition_cover() {
+    std::printf("[partition cover]\n");
+    BBox area{0, 0, 100000, 100000};  // 100um die
+    PartitionConfig pc;
+    pc.anchor_x = 0;
+    pc.anchor_y = 0;
+    pc.align_x = 20000;  // 20um
+    pc.align_y = 20000;
+    pc.margin = 20000;
+    pc.max_leaf = 2;
+    pc.clip = area;
+    auto parts = partition_recursive(area, pc);
+    CHECK(parts.size() >= 4, "die split into multiple leaves");
+
+    // Cores must be disjoint and exactly cover the die area (by summed area).
+    long long covered = 0;
+    bool aligned = true;
+    for (const auto& p : parts) {
+        covered += (long long)p.core.width() * p.core.height();
+        if ((p.core.xmin % pc.align_x) != 0 || (p.core.ymin % pc.align_y) != 0) aligned = false;
+        CHECK(p.halo.xmin <= p.core.xmin && p.halo.xmax >= p.core.xmax, "halo contains core");
+    }
+    CHECK(aligned, "core edges aligned to lattice");
+    CHECK(covered == (long long)area.width() * area.height(), "cores tile the die exactly");
+}
+
+// The whole point of partitioning: the merged result must match a single-shot
+// (global) fill. We run the engine with a tiny leaf size and with a huge leaf
+// size (one partition) and require identical post-fill density.
+static void test_partition_equals_global() {
+    std::printf("[partition == global]\n");
+    Layout base;
+    base.cell_name = "TOP";
+    const dbu win = 20000;
+    for (int j = 0; j < 5; ++j)
+        for (int i = 0; i < 5; ++i) {
+            double f = 0.05 + 0.05 * ((i + j) % 7);
+            double s = std::sqrt(f) * win;
+            dbu x0 = i * win + dbu((win - s) / 2);
+            dbu y0 = j * win + dbu((win - s) / 2);
+            base.polygons.push_back(make_rect(metal_layer(3), 0, x0, y0, x0 + dbu(s), y0 + dbu(s)));
+        }
+
+    std::vector<FillRule> rules;
+    for (const auto& r : default_layermap())
+        if (r.name == "M3") rules.push_back(r);
+
+    Layout a = base, b = base;
+    EngineConfig ca;
+    ca.max_leaf = 1;      // finely partitioned
+    ca.max_iterations = 3;
+    EngineConfig cb;
+    cb.max_leaf = 100000; // effectively one partition (global)
+    cb.max_iterations = 3;
+
+    FillSummary sa = run_fill(a, rules, ca);
+    FillSummary sb = run_fill(b, rules, cb);
+
+    CHECK(sa.layers[0].partitions > 1, "fine config actually partitioned");
+    CHECK(sb.layers[0].partitions == 1, "global config is a single partition");
+    CHECK(sa.layers[0].fill_shapes == sb.layers[0].fill_shapes,
+          "same number of fill shapes partitioned vs global");
+    CHECK(std::fabs(sa.layers[0].density_after_mean - sb.layers[0].density_after_mean) < 1e-6,
+          "same post-fill mean density");
+    CHECK(std::fabs(sa.layers[0].density_after_min - sb.layers[0].density_after_min) < 1e-6,
+          "same post-fill min density");
+}
+
+static void test_engine_endtoend() {
+    std::printf("[engine end-to-end]\n");
+    Layout L;
+    L.cell_name = "TOP";
+    const dbu win = 20000;
+    auto rules = default_layermap();
+    int li = 0;
+    for (const auto& r : rules) {
+        for (int j = 0; j < 5; ++j)
+            for (int i = 0; i < 5; ++i) {
+                double f = 0.04 + 0.06 * ((i + j + li) % 9);
+                double s = std::sqrt(f) * win;
+                dbu x0 = i * win + dbu((win - s) / 2);
+                dbu y0 = j * win + dbu((win - s) / 2);
+                L.polygons.push_back(make_rect(r.layer, r.datatype, x0, y0, x0 + dbu(s), y0 + dbu(s)));
+            }
+        ++li;
+    }
+    size_t before = L.polygons.size();
+    EngineConfig cfg;
+    cfg.max_iterations = 3;
+    FillSummary sum = run_fill(L, rules, cfg);
+
+    CHECK(L.polygons.size() > before, "fill polygons appended");
+    CHECK(sum.total_fill_shapes() > 0, "fill produced");
+    // No spacing or max-density violations may ever be produced by the tool.
+    int spacing = 0, maxd = 0;
+    for (const auto& l : sum.layers) {
+        spacing += l.drc.count(ViolationType::Spacing);
+        maxd += l.drc.count(ViolationType::MaxDensity);
+    }
+    CHECK(spacing == 0, "no fill-in-keepout spacing violations");
+    CHECK(maxd == 0, "no max-density violations");
+    for (const auto& l : sum.layers) {
+        CHECK(l.density_after_min >= l.density_before_min - 1e-9, "min density did not decrease");
+        CHECK(l.density_after_mean >= l.density_before_mean - 1e-9, "mean density did not decrease");
+    }
+}
+
+int main() {
+    std::printf("== metalfill tests ==\n");
+    test_gds_roundtrip();
+    test_sat_density();
+    test_rasterize_rect();
+    test_keepout();
+    test_place_fill();
+    test_partition_cover();
+    test_partition_equals_global();
+    test_engine_endtoend();
+    std::printf("\n%d checks, %d failures\n", g_checks, g_fails);
+    return g_fails == 0 ? 0 : 1;
+}
diff --git a/gpu_metal_fill/tools/gdsinfo.cpp b/gpu_metal_fill/tools/gdsinfo.cpp
new file mode 100644
index 0000000..639d0b4
--- /dev/null
+++ b/gpu_metal_fill/tools/gdsinfo.cpp
@@ -0,0 +1,62 @@
+// Prints a per-(layer,datatype) polygon-count summary for one or more GDSII
+// files. Fill emitted by this project uses datatype 10, so it is flagged as
+// "FILL" to make it easy to confirm that fill happened.
+#include <cstdio>
+#include <map>
+#include <string>
+#include <utility>
+
+#include "metalfill/gdsii.hpp"
+#include "metalfill/layermap.hpp"
+
+using namespace mf;
+
+namespace {
+
+// Best-effort human-readable name for the known layer numbers.
+std::string layer_name(int layer) {
+    if (layer == kOdLayer) return "OD";
+    if (layer == kPoLayer) return "PO";
+    if (layer > kMetalBaseLayer && layer <= kMetalBaseLayer + kNumMetals)
+        return "M" + std::to_string(layer - kMetalBaseLayer);
+    return "?";
+}
+
+void dump(const char* path) {
+    Layout L = read_gds(path);
+    std::map<std::pair<int, int>, long> counts;
+    for (const auto& p : L.polygons) counts[{p.layer, p.datatype}]++;
+
+    BBox b = L.bbox();
+    double dppu = L.dbu_per_um();
+    std::printf("=== %s ===\n", path);
+    std::printf("  cell '%s', %zu polygons, bbox %.2f x %.2f um\n", L.cell_name.c_str(),
+                L.polygons.size(), b.width() / dppu, b.height() / dppu);
+    std::printf("  %-6s %-6s %10s   %s\n", "layer", "L/DT", "polygons", "kind");
+    std::printf("  ------------------------------------------------\n");
+
+    long design = 0, fill = 0;
+    for (const auto& kv : counts) {
+        int layer = kv.first.first, dt = kv.first.second;
+        bool is_fill = (dt == kFillDatatype);
+        if (is_fill) fill += kv.second;
+        else design += kv.second;
+        char ld[32];
+        std::snprintf(ld, sizeof(ld), "%d/%d", layer, dt);
+        std::printf("  %-6s %-6s %10ld   %s\n", layer_name(layer).c_str(), ld, kv.second,
+                    is_fill ? "FILL" : "design");
+    }
+    std::printf("  ------------------------------------------------\n");
+    std::printf("  design=%ld  fill=%ld  total=%ld\n\n", design, fill, design + fill);
+}
+
+}  // namespace
+
+int main(int argc, char** argv) {
+    if (argc < 2) {
+        std::printf("usage: gdsinfo file1.gds [file2.gds ...]\n");
+        return 2;
+    }
+    for (int i = 1; i < argc; ++i) dump(argv[i]);
+    return 0;
+}
diff --git a/gpu_metal_fill/tools/make_dummy_gds.cpp b/gpu_metal_fill/tools/make_dummy_gds.cpp
new file mode 100644
index 0000000..44dc7f2
--- /dev/null
+++ b/gpu_metal_fill/tools/make_dummy_gds.cpp
@@ -0,0 +1,67 @@
+// Generates a dummy GDSII with FEOL base layers (OD, PO) and BEOL metals
+// (M1..M14). Each 20um window is given a deterministic, spatially varying
+// existing density so the fill engine has real work to do (some windows below
+// the min-density target, some already dense, with gradients across the die).
+#include <cmath>
+#include <cstdlib>
+#include <iostream>
+#include <string>
+
+#include "metalfill/gdsii.hpp"
+#include "metalfill/layermap.hpp"
+
+using namespace mf;
+
+int main(int argc, char** argv) {
+    std::string out = "dummy.gds";
+    int windows = 5;  // 5x5 windows -> 100um die
+    for (int i = 1; i < argc; ++i) {
+        std::string a = argv[i];
+        if (a == "-o" && i + 1 < argc)
+            out = argv[++i];
+        else if (a == "-n" && i + 1 < argc)
+            windows = std::atoi(argv[++i]);
+        else if (a == "-h" || a == "--help") {
+            std::cout << "usage: make_dummy_gds [-o out.gds] [-n windows_per_side]\n";
+            return 0;
+        }
+    }
+
+    Layout L;
+    L.lib_name = "METALFILL";
+    L.cell_name = "TOP";
+    const double dppu = L.dbu_per_um();       // 1000 dbu / um
+    const dbu win = static_cast<dbu>(20.0 * dppu);  // 20um window
+
+    auto rules = default_layermap();
+    int layer_index = 0;
+    for (const auto& r : rules) {
+        for (int j = 0; j < windows; ++j) {
+            for (int i = 0; i < windows; ++i) {
+                // Spatially varying density, offset per layer so stacks differ.
+                int phase = (i + j + layer_index) % 9;
+                double f = 0.04 + 0.06 * phase;  // 0.04 .. 0.52
+                if (f < 0.0) f = 0.0;
+                if (f > 0.85) f = 0.85;
+                double s = std::sqrt(f) * win;  // side of a centered square
+                dbu wx0 = static_cast<dbu>(i) * win;
+                dbu wy0 = static_cast<dbu>(j) * win;
+                dbu off = static_cast<dbu>((win - s) / 2.0);
+                dbu x0 = wx0 + off;
+                dbu y0 = wy0 + off;
+                dbu x1 = x0 + static_cast<dbu>(s);
+                dbu y1 = y0 + static_cast<dbu>(s);
+                if (x1 > x0 && y1 > y0)
+                    L.polygons.push_back(make_rect(r.layer, r.datatype, x0, y0, x1, y1));
+            }
+        }
+        ++layer_index;
+    }
+
+    write_gds(out, L);
+    BBox b = L.bbox();
+    std::cout << "wrote " << out << ": " << L.polygons.size() << " polygons, die "
+              << (b.width() / dppu) << " x " << (b.height() / dppu) << " um, " << rules.size()
+              << " layers\n";
+    return 0;
+}
diff --git a/gpu_metal_fill/tools/make_gpu_block.cpp b/gpu_metal_fill/tools/make_gpu_block.cpp
new file mode 100644
index 0000000..105c9e1
--- /dev/null
+++ b/gpu_metal_fill/tools/make_gpu_block.cpp
@@ -0,0 +1,161 @@
+// Generates a synthetic but realistic "GPU block" GDSII to exercise the fill
+// engine. It is NOT a real design -- it just reproduces the *floorplan-level
+// density structure* a GPU block would have so that fill has varied work:
+//
+//   * A grid of SM (streaming-multiprocessor) tiles. Each tile has two SRAM
+//     macros (register file + shared memory) and a standard-cell logic region.
+//   * SRAM macros: very dense on OD/PO/M1..M3 (bitcell arrays), open above.
+//   * Logic regions: medium-density cell rows (OD/PO) + local routing (M1..M3),
+//     with per-tile variation.
+//   * Routing channels between tiles: bus routing on mid metals (M4..M8).
+//   * Global clock/spine routes on M9/M10.
+//   * A regular power grid (wide straps) on the top metals (M11..M14), which is
+//     sparse (~15-20%) and therefore needs fill up to the min-density target.
+//
+// GDS layers follow metalfill/layermap.hpp: OD=1, PO=2, M1..M14=10..23, dt 0.
+#include <cmath>
+#include <cstdint>
+#include <cstdlib>
+#include <iostream>
+#include <string>
+
+#include "metalfill/gdsii.hpp"
+#include "metalfill/layermap.hpp"
+
+using namespace mf;
+
+namespace {
+
+struct Gen {
+    Layout L;
+    double dppu;
+    uint64_t s = 0x9e3779b97f4a7c15ULL;
+
+    dbu um(double v) { return static_cast<dbu>(std::llround(v * dppu)); }
+    uint32_t rnd() {
+        s = s * 6364136223846793005ULL + 1442695040888963407ULL;
+        return static_cast<uint32_t>(s >> 33);
+    }
+    double frand() { return (rnd() % 100000) / 100000.0; }
+
+    void rect(int layer, dbu x0, dbu y0, dbu x1, dbu y1) {
+        if (x1 > x0 && y1 > y0) L.polygons.push_back(make_rect(layer, 0, x0, y0, x1, y1));
+    }
+    // Parallel stripes filling a rectangle. density ~= width/pitch.
+    void stripes(int layer, dbu x0, dbu y0, dbu x1, dbu y1, dbu pitch, dbu width, bool vertical) {
+        if (pitch <= 0) return;
+        if (vertical) {
+            for (dbu x = x0; x + width <= x1; x += pitch) rect(layer, x, y0, x + width, y1);
+        } else {
+            for (dbu y = y0; y + width <= y1; y += pitch) rect(layer, x0, y, x1, y + width);
+        }
+    }
+    // A dense macro: bitcell-like stripes on OD/PO/M1..M3. A guard ring of a
+    // couple of microns is left inside so density tapers toward the edge (as in
+    // a real macro placement), which keeps the density gradient fillable.
+    void macro(dbu x0, dbu y0, dbu x1, dbu y1) {
+        stripes(kOdLayer, x0, y0, x1, y1, um(0.5), um(0.28), false);   // ~56%
+        stripes(kPoLayer, x0, y0, x1, y1, um(0.5), um(0.24), true);    // ~48%
+        stripes(metal_layer(1), x0, y0, x1, y1, um(0.4), um(0.20), true);   // ~50%
+        stripes(metal_layer(2), x0, y0, x1, y1, um(0.4), um(0.18), false);  // ~45%
+        stripes(metal_layer(3), x0, y0, x1, y1, um(0.6), um(0.24), true);   // ~40%
+    }
+    // Standard-cell logic region: rows of active + poly gates, local routing.
+    void logic(dbu x0, dbu y0, dbu x1, dbu y1, double lo) {
+        // Cell rows: OD islands per row, poly gates crossing them.
+        dbu row = um(1.2);
+        for (dbu y = y0; y + um(0.7) <= y1; y += row) {
+            for (dbu x = x0; x + um(1.0) <= x1; x += um(1.4)) {
+                if (frand() < lo + 0.35) rect(kOdLayer, x, y, x + um(0.9), y + um(0.6));
+            }
+        }
+        stripes(kPoLayer, x0, y0, x1, y1, um(0.9), um(0.18), true);  // poly gates
+        // Local routing, density varies per region.
+        stripes(metal_layer(1), x0, y0, x1, y1, um(0.6), um(0.2 + 0.15 * lo), true);
+        stripes(metal_layer(2), x0, y0, x1, y1, um(0.8), um(0.24 + 0.2 * lo), false);
+        stripes(metal_layer(3), x0, y0, x1, y1, um(1.0), um(0.24 + 0.2 * lo), true);
+    }
+};
+
+}  // namespace
+
+int main(int argc, char** argv) {
+    std::string out = "gpu_block.gds";
+    int sm = 3;  // SM grid is sm x sm
+    for (int i = 1; i < argc; ++i) {
+        std::string a = argv[i];
+        if (a == "-o" && i + 1 < argc) out = argv[++i];
+        else if (a == "-n" && i + 1 < argc) sm = std::atoi(argv[++i]);
+        else if (a == "-h" || a == "--help") {
+            std::cout << "usage: make_gpu_block [-o out.gds] [-n sm_grid]\n";
+            return 0;
+        }
+    }
+
+    Gen g;
+    g.dppu = g.L.dbu_per_um();
+    g.L.lib_name = "GPU_BLOCK";
+    g.L.cell_name = "GPU_TOP";
+
+    const dbu margin = g.um(8);
+    const dbu tile = g.um(56);
+    const dbu chan = g.um(6);
+    const dbu step = tile + chan;
+    const dbu die = margin * 2 + sm * tile + (sm - 1) * chan;
+
+    // ---- SM tiles: macros + logic ----------------------------------------
+    for (int ty = 0; ty < sm; ++ty) {
+        for (int tx = 0; tx < sm; ++tx) {
+            dbu ox = margin + tx * step;
+            dbu oy = margin + ty * step;
+            // Two SRAM macros in the top half (register file + shared memory).
+            dbu mac_h = g.um(24);
+            dbu half = tile / 2;
+            g.macro(ox + g.um(1), oy + tile - mac_h - g.um(1), ox + half - g.um(1), oy + tile - g.um(1));
+            g.macro(ox + half + g.um(1), oy + tile - mac_h - g.um(1), ox + tile - g.um(1),
+                    oy + tile - g.um(1));
+            // Logic occupies the bottom part; density varies per tile.
+            double lo = 0.30 + 0.4 * g.frand();
+            g.logic(ox + g.um(1), oy + g.um(1), ox + tile - g.um(1), oy + tile - mac_h - g.um(2), lo);
+        }
+    }
+
+    // ---- Routing channels: bus routing on mid metals ---------------------
+    for (int t = 1; t < sm; ++t) {
+        dbu cx = margin + t * step - chan;
+        dbu cy = margin + t * step - chan;
+        for (int m = 4; m <= 8; ++m) {
+            bool vert = (m % 2 == 0);
+            // vertical channel (routes run vertically) and horizontal channel
+            g.stripes(metal_layer(m), cx, margin, cx + chan, die - margin, g.um(0.8), g.um(0.4),
+                      vert);
+            g.stripes(metal_layer(m), margin, cy, die - margin, cy + chan, g.um(0.8), g.um(0.4),
+                      !vert);
+        }
+    }
+
+    // ---- Global clock / spine on M9, M10 ---------------------------------
+    for (int k = 0; k < sm; ++k) {
+        dbu c = margin + k * step + tile / 2;
+        g.rect(metal_layer(9), c - g.um(1), margin, c + g.um(1), die - margin);
+        g.rect(metal_layer(10), margin, c - g.um(1), die - margin, c + g.um(1));
+    }
+
+    // ---- Power grid on top metals (wide, sparse straps) ------------------
+    // Vertical on M12/M14, horizontal on M11/M13; ~4um straps on a 24um pitch.
+    for (dbu x = margin; x + g.um(4) <= die - margin; x += g.um(24)) {
+        g.rect(metal_layer(12), x, margin, x + g.um(4), die - margin);
+        g.rect(metal_layer(14), x, margin, x + g.um(4), die - margin);
+    }
+    for (dbu y = margin; y + g.um(4) <= die - margin; y += g.um(24)) {
+        g.rect(metal_layer(11), margin, y, die - margin, y + g.um(4));
+        g.rect(metal_layer(13), margin, y, die - margin, y + g.um(4));
+    }
+
+    write_gds(out, g.L);
+    BBox b = g.L.bbox();
+    std::cout << "wrote " << out << ": " << g.L.polygons.size() << " polygons, GPU block "
+              << (b.width() / g.dppu) << " x " << (b.height() / g.dppu) << " um, " << (sm * sm)
+              << " SM tiles\n";
+    return 0;
+}
diff --git a/gpu_metal_fill/tools/render_layer.cpp b/gpu_metal_fill/tools/render_layer.cpp
new file mode 100644
index 0000000..466ede4
--- /dev/null
+++ b/gpu_metal_fill/tools/render_layer.cpp
@@ -0,0 +1,116 @@
+// Renders one layer of a GDSII to a PPM image: existing geometry, dummy fill,
+// and the recursive-partition (quadtree) core boundaries. Useful for eyeballing
+// what the fill engine produced.
+#include <algorithm>
+#include <cstdint>
+#include <cstdio>
+#include <cstdlib>
+#include <fstream>
+#include <string>
+#include <vector>
+
+#include "metalfill/engine.hpp"
+#include "metalfill/gdsii.hpp"
+#include "metalfill/layermap.hpp"
+#include "metalfill/partition.hpp"
+#include "metalfill/raster.hpp"
+
+using namespace mf;
+
+struct Img {
+    int w, h;
+    std::vector<uint8_t> px;  // RGB
+    Img(int w_, int h_) : w(w_), h(h_), px(size_t(w_) * h_ * 3, 255) {}
+    void set(int x, int y, uint8_t r, uint8_t g, uint8_t b) {
+        if (x < 0 || y < 0 || x >= w || y >= h) return;
+        size_t i = (size_t(y) * w + x) * 3;
+        px[i] = r;
+        px[i + 1] = g;
+        px[i + 2] = b;
+    }
+    void hline(int x0, int x1, int y, uint8_t r, uint8_t g, uint8_t b) {
+        for (int x = x0; x <= x1; ++x) set(x, y, r, g, b);
+    }
+    void vline(int x, int y0, int y1, uint8_t r, uint8_t g, uint8_t b) {
+        for (int y = y0; y <= y1; ++y) set(x, y, r, g, b);
+    }
+    void write_ppm(const std::string& path) {
+        std::ofstream os(path, std::ios::binary);
+        os << "P6\n" << w << " " << h << "\n255\n";
+        os.write(reinterpret_cast<char*>(px.data()), px.size());
+    }
+};
+
+int main(int argc, char** argv) {
+    std::string in, out = "layer.ppm";
+    int target_layer = metal_layer(3);  // M3 by default
+    int px_target = 900;
+    for (int i = 1; i < argc; ++i) {
+        std::string a = argv[i];
+        if (a == "-i" && i + 1 < argc) in = argv[++i];
+        else if (a == "-o" && i + 1 < argc) out = argv[++i];
+        else if (a == "-l" && i + 1 < argc) target_layer = std::atoi(argv[++i]);
+        else if (a == "-p" && i + 1 < argc) px_target = std::atoi(argv[++i]);
+    }
+    if (in.empty()) { std::fprintf(stderr, "usage: render_layer -i in.gds -o out.ppm [-l layer]\n"); return 2; }
+
+    Layout L = read_gds(in);
+    BBox area = L.bbox();
+    if (!area.valid()) { std::fprintf(stderr, "empty layout\n"); return 1; }
+
+    dbu cell = std::max<dbu>(area.width() / px_target, 1);
+    Grid ge = make_grid(area, cell), gf = make_grid(area, cell);
+    rasterize(ge, L.polygons, target_layer, 0);
+    rasterize(gf, L.polygons, target_layer, kFillDatatype);
+
+    int W = ge.nx, H = ge.ny;
+    Img img(W, H);
+
+    auto to_px = [&](dbu x, dbu y, int& ix, int& iy) {
+        ix = int((x - area.xmin) / cell);
+        iy = H - 1 - int((y - area.ymin) / cell);
+    };
+
+    // Light window grid every 20um.
+    dbu win = static_cast<dbu>(20.0 * L.dbu_per_um());
+    for (dbu x = area.xmin; x <= area.xmax; x += win) {
+        int ix, iy;
+        to_px(x, area.ymin, ix, iy);
+        img.vline(ix, 0, H - 1, 235, 235, 235);
+    }
+    for (dbu y = area.ymin; y <= area.ymax; y += win) {
+        int ix, iy;
+        to_px(area.xmin, y, ix, iy);
+        img.hline(0, W - 1, iy, 235, 235, 235);
+    }
+
+    // Existing geometry (navy) and fill (orange).
+    for (int y = 0; y < H; ++y)
+        for (int x = 0; x < W; ++x) {
+            int row = H - 1 - y;
+            if (ge.at(x, y)) img.set(x, row, 30, 60, 150);
+            else if (gf.at(x, y)) img.set(x, row, 250, 160, 60);
+        }
+
+    // Recursive-partition (quadtree) core boundaries in red.
+    dbu align = lcm_dbu(win, win);
+    PartitionConfig pc;
+    pc.anchor_x = area.xmin; pc.anchor_y = area.ymin;
+    pc.align_x = align; pc.align_y = align;
+    pc.margin = align; pc.max_leaf = 2; pc.clip = area;
+    for (const auto& part : partition_recursive(area, pc)) {
+        int x0, y0, x1, y1;
+        to_px(part.core.xmin, part.core.ymin, x0, y0);
+        to_px(part.core.xmax, part.core.ymax, x1, y1);
+        if (x1 < x0) std::swap(x0, x1);
+        if (y1 < y0) std::swap(y0, y1);
+        img.hline(x0, x1, y0, 220, 30, 30);
+        img.hline(x0, x1, y1, 220, 30, 30);
+        img.vline(x0, y0, y1, 220, 30, 30);
+        img.vline(x1, y0, y1, 220, 30, 30);
+    }
+
+    img.write_ppm(out);
+    std::printf("wrote %s (%dx%d) layer=%d\n", out.c_str(), W, H, target_layer);
+    return 0;
+}
diff --git a/gpu_metal_fill/tools/run_fill.cpp b/gpu_metal_fill/tools/run_fill.cpp
new file mode 100644
index 0000000..7a01470
--- /dev/null
+++ b/gpu_metal_fill/tools/run_fill.cpp
@@ -0,0 +1,58 @@
+// Reads a GDSII, runs recursive-partitioned FEOL/BEOL metal fill, writes the
+// filled GDSII and a text report.
+#include <cstdlib>
+#include <fstream>
+#include <iostream>
+#include <string>
+
+#include "metalfill/engine.hpp"
+#include "metalfill/gdsii.hpp"
+#include "metalfill/layermap.hpp"
+
+using namespace mf;
+
+int main(int argc, char** argv) {
+    std::string in, out = "filled.gds", report;
+    EngineConfig cfg;
+    for (int i = 1; i < argc; ++i) {
+        std::string a = argv[i];
+        if ((a == "-i" || a == "--in") && i + 1 < argc)
+            in = argv[++i];
+        else if ((a == "-o" || a == "--out") && i + 1 < argc)
+            out = argv[++i];
+        else if ((a == "-r" || a == "--report") && i + 1 < argc)
+            report = argv[++i];
+        else if (a == "--iters" && i + 1 < argc)
+            cfg.max_iterations = std::atoi(argv[++i]);
+        else if (a == "--gpu")
+            cfg.prefer_gpu = true;
+        else if (a == "-h" || a == "--help") {
+            std::cout << "usage: run_fill -i in.gds [-o filled.gds] [-r report.txt] "
+                         "[--iters N] [--gpu]\n";
+            return 0;
+        } else if (in.empty()) {
+            in = a;
+        }
+    }
+    if (in.empty()) {
+        std::cerr << "error: no input GDS (use -i in.gds)\n";
+        return 2;
+    }
+
+    Layout layout = read_gds(in);
+    std::cout << "read " << in << ": " << layout.polygons.size() << " polygons\n";
+
+    auto rules = default_layermap();
+    FillSummary summary = run_fill(layout, rules, cfg);
+
+    write_gds(out, layout);
+    std::string rpt = format_report(summary);
+    std::cout << rpt;
+    std::cout << "wrote filled GDS: " << out << " (" << layout.polygons.size() << " polygons)\n";
+    if (!report.empty()) {
+        std::ofstream f(report);
+        f << rpt;
+        std::cout << "wrote report: " << report << "\n";
+    }
+    return summary.total_violations() == 0 ? 0 : 1;
+}