code: def get_preconditioner():
      """Compute the preconditioner M"""
      diags_x = zeros((3, nx))
      diags_x[0, :] = 1 / hx / hx
      diags_x[1, :] = -2 / hx / hx
      diags_x[2, :] = 1 / hx / hx
      Lx = spdiags(diags_x, [-1, 0, 1], nx, nx)
  
      diags_y = zeros((3, ny))
      diags_y[0, :] = 1 / hy / hy
      diags_y[1, :] = -2 / hy / hy
      diags_y[2, :] = 1 / hy / hy
      Ly = spdiags(diags_y, [-1, 0, 1], ny, ny)
  
      J1 = kron(Lx, eye(ny)) + kron(eye(nx), Ly)
  
      # Now we have the matrix `J_1`. We need to find its inverse `M` --
      # however, since an approximate inverse is enough, we can use
      # the *incomplete LU* decomposition
  
      J1_ilu = spilu(J1)
  
      # This returns an object with a method .solve() that evaluates
      # the corresponding matrix-vector product. We need to wrap it into
      # a LinearOperator before it can be passed to the Krylov methods:
  
      M = LinearOperator(shape=(nx * ny, nx * ny), matvec=J1_ilu.solve)
      return M
  
imports: from scipy.optimize import root
  from scipy.sparse import spdiags, kron
  from scipy.sparse.linalg import spilu, LinearOperator
  from numpy import cosh, zeros_like, mgrid, zeros, eye
label_desc: Use image node_id 1 for calling a global function with example usage: get_preconditioner() and returns: M