Fix underflow issue when computing normalization

include/igl/heat_geodesics.cpp

@@ -84,8 +84,7 @@ IGL_INLINE bool igl::heat_geodesics_precompute(
         return false;
       }
     }
-    const DerivedV M_diag_tr = M.diagonal().transpose();
-    const Eigen::SparseMatrix<Scalar> Aeq = M_diag_tr.sparseView();
+    const Eigen::SparseMatrix<double> Aeq = M.diagonal().transpose().sparseView();
     L *= -0.5;
     if(!igl::min_quad_with_fixed_precompute(
       L,Eigen::VectorXi(),Aeq,true,data.Poisson))
@@ -127,13 +126,22 @@ IGL_INLINE void igl::heat_geodesics_solve(
   const int m = data.Grad.rows()/data.ng;
   for(int i = 0;i<m;i++)
   {
+    // It is very important to use a stable norm calculation here. If the
+    // triangle is far from a source, then the floating point values in the
+    // gradient can be _very_ small (e.g., 1e-300). The standard/naive norm
+    // calculation will suffer from underflow. Dividing by the max value is more
+    // stable. (Eigen implements this as stableNorm or blueNorm).
     Scalar norm = 0;
+    Scalar ma = 0;
+    for(int d = 0;d<data.ng;d++) {ma = std::max(ma,std::fabs(grad_u(d*m+i)));}
     for(int d = 0;d<data.ng;d++)
     {
-      norm += grad_u(d*m+i)*grad_u(d*m+i);
+      const Scalar gui = grad_u(d*m+i) / ma;
+      norm += gui*gui;
     }
-    norm = sqrt(norm);
-    if(norm == 0)
+    norm = ma*sqrt(norm);
+    // These are probably over kill; ma==0 should be enough
+    if(ma == 0 || norm == 0 || norm!=norm)
     {
       for(int d = 0;d<data.ng;d++) { grad_u(d*m+i) = 0; }
     }else