invert diag

Former-commit-id: df9b25d7636dd2a27c02efcdfa841f8db4986a43

include/igl/invert_diag.cpp

@@ -0,0 +1,39 @@
+#include "invert_diag.h"
+#include "diag.h"
+
+template <typename T>
+IGL_INLINE void igl::invert_diag(
+  const Eigen::SparseMatrix<T>& X, 
+  Eigen::SparseMatrix<T>& Y)
+{
+#ifndef NDEBUG
+  typename Eigen::SparseVector<T> dX;
+  igl::diag(X,dX);
+  // Check that there are no zeros along the diagonal
+  assert(dX.nonZeros() == dX.size());
+#endif
+  // http://www.alecjacobson.com/weblog/?p=2552
+  if(&Y != &X)
+  {
+    Y = X;
+  }
+  // Iterate over outside
+  for(int k=0; k<Y.outerSize(); ++k)
+  {
+    // Iterate over inside
+    for(typename Eigen::SparseMatrix<T>::InnerIterator it (Y,k); it; ++it)
+    {
+      if(it.col() == it.row())
+      {
+        T v = it.value();
+        assert(v != 0);
+        v = ((T)1.0)/v;
+        Y.coeffRef(it.row(),it.col()) = v;
+      }
+    }
+  }
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+#endif

include/igl/invert_diag.h

@@ -0,0 +1,29 @@
+#ifndef IGL_INVERT_DIAG_H
+#define IGL_INVERT_DIAG_H
+#include "igl_inline.h"
+#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
+#include <Eigen/Sparse>
+
+namespace igl
+{
+  // Invert the diagonal entries of a matrix (if the matrix is a diagonal
+  // matrix then this amounts to inverting the matrix)
+
+  // Templates:
+  //   T  should be a eigen sparse matrix primitive type like int or double
+  // Inputs:
+  //   X  an m by n sparse matrix
+  // Outputs:
+  //   Y  an m by n sparse matrix
+  template <typename T>
+  IGL_INLINE void invert_diag(
+    const Eigen::SparseMatrix<T>& X, 
+    Eigen::SparseMatrix<T>& Y);
+}
+
+#ifdef IGL_HEADER_ONLY
+#  include "invert_diag.cpp"
+#endif
+
+#endif
+