allow passing use_cholmod flag (#8)

lapy/Conformal.py

@@ -368,12 +368,11 @@ def sparse_symmetric_solve(A, b, use_cholmod=True):
         except ImportError:
             use_cholmod = False
     if use_cholmod:
-        print("Solver: cholesky decomp - performance optimal ...")
+        print("Solver: Cholesky decomposition from scikit-sparse cholmod ...")
         chol = cholesky(A)
         x = chol(b)
     else:
-        print("Package scikit-sparse not found (Cholesky decomp)")
-        print("Solver: spsolve (LU decomp) - performance not optimal ...")
+        print("Solver: spsolve (LU decomposition) ...")
         lu = splu(A)
         x = lu.solve(b)
     return x    

lapy/DiffGeo.py

@@ -248,7 +248,7 @@ def tria_compute_rotated_f(tria, vfunc):
     return vf
 
 
-def tria_mean_curvature_flow(tria, max_iter=30, stop_eps=1e-13, step=1.0):
+def tria_mean_curvature_flow(tria, max_iter=30, stop_eps=1e-13, step=1.0, use_cholmod=True):
     """
     mean_curvature_flow iteratively flows a triangle mesh along mean curvature
     normal (non-singular, see Kazhdan 2012)
@@ -266,14 +266,16 @@ def tria_mean_curvature_flow(tria, max_iter=30, stop_eps=1e-13, step=1.0):
     steps. It will normalize surface area of the mesh and translate the barycenter
     to the origin. Closed meshes will map to the unit sphere. 
     """
-    use_cholmod = True
-    try:
-        from sksparse.cholmod import cholesky
-    except ImportError:
-        use_cholmod = False
-        from scipy.sparse.linalg import spsolve
+    if use_cholmod:
+        try:
+            from sksparse.cholmod import cholesky
+        except ImportError:
+            use_cholmod = False
+
+    if not use_cholmod:
         # Note, it would be better to do sparse Cholesky (CHOLMOD)
         # as it can be 5-6 times faster
+        from scipy.sparse.linalg import spsolve
     # pre-normalize
     trianorm = TriaMesh(tria.v, tria.t)
     trianorm.normalize_()
@@ -282,9 +284,9 @@ def tria_mean_curvature_flow(tria, max_iter=30, stop_eps=1e-13, step=1.0):
     fem = Solver(trianorm, lump)
     a_mat = fem.stiffness
     if use_cholmod:
-        print("Solver: cholesky decomp - performance optimal ...")
+        print("Solver: Cholesky decomposition from scikit-sparse cholmod ...")
     else:
-        print("Solver: spsolve (LU decomp) - performance not optimal ...")
+        print("Solver: spsolve (LU decomposition) ...")
     for x in range(max_iter):
         # store last position (for delta computation below)
         vlast = trianorm.v

lapy/Heat.py

@@ -57,7 +57,7 @@ def diffusion(geometry, vids, m=1.0, aniso=None, use_cholmod=True):
       vids          vertex index or indices where initial heat is applied
       m             factor (default 1) to compute time of heat evolution:
                     t = m * avg_edge_length^2
-      use_cholmod   (default True), if Cholmod not installed automatically
+      use_cholmod   (default True), if Cholmod is not found
                     revert to LU decomposition (slower)
 
     :return:
@@ -83,12 +83,11 @@ def diffusion(geometry, vids, m=1.0, aniso=None, use_cholmod=True):
     # solve H x = b0
     print("Matrix Format now:  " + hmat.getformat())
     if use_cholmod:
-        print("Solver: cholesky decomp - performance optimal ...")
+        print("Solver: Cholesky decomposition from scikit-sparse cholmod ...")
         chol = cholesky(hmat)
         vfunc = chol(b0)
     else:
-        print("Package scikit-sparse not found (Cholesky decomp)")
-        print("Solver: spsolve (LU decomp) - performance not optimal ...")
+        print("Solver: spsolve (LU decomposition) ...")
         lu = splu(hmat)
         vfunc = lu.solve(np.float32(b0))
     return vfunc

lapy/Solver.py

@@ -16,22 +16,27 @@ class Solver:
      for external function that do not need stiffness.
      """
 
-    def __init__(self, geometry, lump=False, aniso=None, aniso_smooth=10):
+    def __init__(self, geometry, lump=False, aniso=None, aniso_smooth=10, use_cholmod=True):
         """
         Construct the Solver class. Computes linear Laplace FEM stiffness and
         mass matrix for TriaMesh or TetMesh input geometries. For TriaMesh it can also
         construct the anisotropic Laplace.
 
         Inputs:
-            geometry  : is a geometry class, currently either TriaMesh or TetMesh
-            aniso     : float, anisotropy for curvature based anisotopic Laplace
-                        can also be tuple (a_min, a_max) to differentially affect
-                        the min and max curvature directions. E.g. (0,50) will set
-                        scaling to 1 into min curv direction even if the max curvature
-                        is large in those regions (= isotropic in regions with 
-                        large max curv and min curv close to zero= concave cylinder)
-            lump      : whether to lump the mass matrix (diagonal), default False
+            geometry    : is a geometry class, currently either TriaMesh or TetMesh
+            aniso       : float, anisotropy for curvature based anisotopic Laplace
+                          can also be tuple (a_min, a_max) to differentially affect
+                          the min and max curvature directions. E.g. (0,50) will set
+                          scaling to 1 into min curv direction even if the max curvature
+                          is large in those regions (= isotropic in regions with 
+                          large max curv and min curv close to zero= concave cylinder)
+            lump        : whether to lump the mass matrix (diagonal), default False
+            use_cholmod : try to use the Cholesky decomposition from the cholmod
+                          library for improved speed. This requires skikit sparse to
+                          be installed. If it cannot be found, we fallback to LU
+                          decomposition. 
         """
+        self.use_cholmod = use_cholmod
         if type(geometry).__name__ == "TriaMesh":
             if aniso is not None:
                 # anisotropic Laplace
@@ -456,7 +461,7 @@ def _fem_voxels(vox, lump=False):
         b = sparse.csc_matrix((local_b, (i, j)))
         return a, b
 
-    def eigs(self, k=10, use_cholmod=True):
+    def eigs(self, k=10):
         """
         Compute linear finite-element method Laplace-Beltrami spectrum
 
@@ -466,18 +471,19 @@ def eigs(self, k=10, use_cholmod=True):
                     eigenfunctions  array: (N x k) with k eigenfunctions (in the columns)
         """
         from scipy.sparse.linalg import LinearOperator, eigsh, splu
-        try:
-            from sksparse.cholmod import cholesky
-        except ImportError:
-            use_cholmod = False
-        if use_cholmod:
-            print("Solver: cholesky decomp - performance optimal ...")
+        if self.use_cholmod:
+            try:
+                from sksparse.cholmod import cholesky
+            except ImportError:
+                self.use_cholmod = False
+
+        if self.use_cholmod:
+            print("Solver: Cholesky decomposition from scikit-sparse cholmod ...")
         else:
-            print("Package scikit-sparse not found (Cholesky decomp)")
-            print("Solver: spsolve (LU decomp) - performance not optimal ...")
-        # turns out it is much faster to use cholesky and pass operator
+            print("Solver: spsolve (LU decomposition) ...")
+            # turns out it is much faster to use cholesky and pass operator
         sigma = -0.01
-        if use_cholmod:
+        if self.use_cholmod:
             chol = cholesky(self.stiffness - sigma * self.mass)
             op_inv = LinearOperator(matvec=chol, shape=self.stiffness.shape,
                                     dtype=self.stiffness.dtype)
@@ -488,7 +494,7 @@ def eigs(self, k=10, use_cholmod=True):
         eigenvalues, eigenvectors = eigsh(self.stiffness, k, self.mass, sigma=sigma, OPinv=op_inv)
         return eigenvalues, eigenvectors
 
-    def poisson(self, h=0.0, dtup=(), ntup=(), use_cholmod=True):
+    def poisson(self, h=0.0, dtup=(), ntup=()):
         """
         poissonSolver solves the poisson equation with boundary conditions
                based on the A and B Laplace matrices:  A x = B h
@@ -505,17 +511,15 @@ def poisson(self, h=0.0, dtup=(), ntup=(), use_cholmod=True):
                   ntup - Neumann boundary condition as a tuple.
                          Tuple contains index and data arrays of same length.
                          Default: Neumann on all boundaries
-                  useCholmod - default True: try to find Cholesky decomp from sksparse
-                               falls back to scipy sparse splu if not found.
 
         Outputs:  x - array with vertex values of solution
         """
         from scipy.sparse.linalg import splu
-        if use_cholmod:
+        if self.use_cholmod:
             try:
                 from sksparse.cholmod import cholesky
             except ImportError:
-                use_cholmod = False
+                self.use_cholmod = False
         # check matrices
         dim = self.stiffness.shape[0]
         if self.stiffness.shape != self.mass.shape or self.stiffness.shape[1] != dim:
@@ -575,13 +579,12 @@ def poisson(self, h=0.0, dtup=(), ntup=(), use_cholmod=True):
             a = self.stiffness
         # solve A x = b
         print("Matrix Format now: " + a.getformat())
-        if use_cholmod:
-            print("Solver: cholesky decomp - performance optimal ...")
+        if self.use_cholmod:
+            print("Solver: Cholesky decomposition from scikit-sparse cholmod ...")
             chol = cholesky(a)
             x = chol(b)
         else:
-            print("Package scikit-sparse not found (Cholesky decomp)")
-            print("Solver: spsolve (LU decomp) - performance not optimal ...")
+            print("Solver: spsolve (LU decomposition) ...")
             lu = splu(a)
             x = lu.solve(b.astype(np.float32))
         x = np.squeeze(np.array(x))