Start from line 192 at "curvature_perFace.compute", the decomposition is LU, not LDLT. I'm unsure about the stability for LU on a normal equation.