SciPy 3-D Laplacian: spsolve 118 s, symmetric splu 41 s, cg 0.4 s
Measured on a 7-point 3-D conduction matrix (SPD, N=157k, nnz=1.0M, SciPy 1.18.1). spsolve with the default COLAMD ordering took 118 s. splu(A.tocsc(), permc_spec='MMD_AT_PLUS_A', options=dict(SymmetricMode=True)) took 41 s, because a symmetric ordering on A+A^T suits SPD matrices and COLAMD does not. At N=512k, splu was killed for memory. cg(A, b, rtol=1e-9) took 0.37 s (315 iterations) and was within 5e-5 K of the direct solution. A Jacobi preconditioner gave no gain on a single-material grid (324 iterations), since the diagonal is nearly constant. pyamg smoothed aggregation cut iterations to 8-9 at similar wall time. Use the direct solver only when one factorisation serves many right-hand sides on small grids. Pin OPENBLAS_NUM_THREADS before timing cg, because a busy box made it 13x slower.
Author re-run (same script, N=157464, nnz=1017532, SciPy 1.18.1, OPENBLAS_NUM_THREADS=1, load average 0.9): cg 0.39 s and 315 iterations (twice, identical), Jacobi-cg 0.46 s and 324 iterations, symmetric splu 33.7 s (41 s in the original run on a busier box; I did not separate load from run-to-run variance), max |cg - direct| 5.2e-5 K. The conclusions stand. Iteration counts are deterministic and reproduce exactly, so prefer them over seconds when comparing solvers across machines.