LU Dense Solver

This project investigates the computational cost of solving large dense linear systems of the form Ax = b, where A is a square coefficient matrix, b is a known right-hand-side vector, and x is the unknown solution. Such systems arise in numerical modeling and engineering calculations in which many variables interact. Because a dense matrix stores entries throughout its dimensions, both arithmetic work and memory requirements grow rapidly with problem size.

The study compares a CPU baseline with naive and optimized CUDA implementations to examine how much of this workload can benefit from GPU parallelism. The objective is to accelerate the complete solution process while understanding the relationship between algorithmic dependencies, memory movement, and measured execution time.

LU decomposition expresses the coefficient matrix as the product of a lower triangular matrix L and an upper triangular matrix U. Once these factors are available, the original system is solved in two stages: forward substitution solves Ly = b for an intermediate vector y, and backward substitution solves Ux = y for x. This separates the expensive factorization from the subsequent triangular solves.

During elimination, a pivot is used to construct multipliers that remove entries below the diagonal. Those multipliers then update the remaining matrix entries. Dense LU factorization requires work that grows approximately with the cube of the matrix dimension, whereas each triangular solve grows quadratically. The comparison therefore centers on accelerating elimination while accounting for the additional work required to return a solution.

The CUDA implementation distributes independent matrix updates across GPU threads. At a given elimination stage, many entries in the remaining submatrix can be updated concurrently once the required pivot information and multipliers are available. However, successive elimination stages are dependent: a later stage must use values produced by earlier updates. The available parallelism therefore lies primarily within each stage rather than across the entire sequence of pivots.

The optimized version uses shared-memory tiling and improved memory access to reduce the cost of repeatedly accessing matrix values. Tiling allows threads in a block to cooperate on portions of the matrix and reuse data closer to the processors. These optimizations address data movement as well as arithmetic throughput; their value depends on whether the reduction in memory traffic outweighs coordination and execution overhead.

For the reported 5000 × 5000 benchmark, the CPU implementation required approximately 120.4 seconds, compared with 4.68 seconds for the naive CUDA version and approximately 3.0 seconds for the optimized version. Using these rounded measurements, the naive implementation achieved about 25.7× speedup over the CPU baseline, while the optimized implementation achieved about 40.1×, consistent with the project's approximate 40–41× reported result.

Relative to naive CUDA, the optimized version reduced total runtime by approximately 36%, corresponding to an additional 1.56× speedup. These measurements demonstrate the benefit of both GPU execution and implementation refinement for the evaluated workload. They characterize this benchmark and its implementations; performance at other matrix sizes or on different hardware requires separate measurement.

A central challenge was distinguishing improvements in the LU kernel from improvements in the complete solver. Total runtime includes work outside factorization, such as allocation, memory transfers, and other solution stages. A substantially faster kernel can therefore produce a smaller end-to-end improvement when these remaining costs occupy a significant portion of execution.

The project reinforced the importance of identifying which operations are independent, placing reusable data efficiently, and measuring performance at a scope that matches the intended claim. It also separates speed from numerical reliability: execution time alone does not establish solution accuracy or robustness for poorly conditioned matrices. Numerical validation and a clearly documented treatment of pivots remain distinct considerations when assessing a solver's suitability for scientific use.