Solution (source code)

= Solution

Extend a grid vector $v$ by zero to the <Dirichlet boundary condition>[Dirichlet boundary]. Pairing contributions along undirected nearest-neighbour edges gives the discrete energy identity
$$
\boxed{v^TAv=-\sum_{\{r,s\}\in E}(v_r-v_s)^2,}
$$
where $E$ includes edges from an interior node to a boundary node. This is the three-dimensional version of <summation by parts> for the <seven-point Dirichlet Laplacian>.

The right side is nonpositive. If it vanishes, every pair of neighbouring values agrees; connectivity of the grid and the zero boundary values then imply $v=0$. Thus $v^TAv<0$ for every nonzero $v$, so $A$ is <negative-definite matrix>[negative definite]. In particular, zero is not an <eigenvalue>, and therefore
$$
\boxed{A\text{ is nonsingular}.}
$$