stationary precond rewritten with CG

integrators/_details/interfaces/tpsPreconditioner.py

@@ -115,38 +115,18 @@ class TPSPreconditioner(_preconditionerInterface._PreconditionerInterface):
       A = self._A_stat_precond[:self._A_stat_precond.shape[0]//3 \
         , :self._A_stat_precond.shape[0]//3]
 
-    Ainv = np.zeros(A.shape, dtype=float)
-
-    e = np.zeros(A.shape[0], dtype=float)
-    e[0] = 1.0
-
-    print("Inverting matrix for practical preconditioner")
-    sol, succ = scipy.sparse.linalg.cg(A, e \
-      , tol=self._solvetol, maxiter=4000)
-    idx = (abs(sol) > self._solvetol)
-    Ainv[idx,0] = sol[idx]
-    for j in range(1, A.shape[0]):
-      e[j-1] = 0.0
-      e[j] = 1.0
-      sol, succ = scipy.sparse.linalg.cg(A, e, tol=self._solvetol, maxiter=4000)
-      idx = (abs(sol) > self._solvetol)
-      Ainv[idx, j] = sol[idx]
-      if (j % 100 == 0):
-        print("Finished ", j, "out of ", A.shape[0] \
-          , "cols of practical preconditioner")
-
-    N = A.shape[0]
-    practical3D = np.zeros((N*3, N*3), dtype=float)
-    practical3D[:N, :N] = Ainv
-    practical3D[N:2*N, N:2*N] = Ainv
-    practical3D[2*N:3*N, 2*N:3*N] = Ainv
-
-    self._prePractical = sp.csr_matrix(practical3D)
-    self._prePractical.eliminate_zeros()
     
     self._Compute_Q()
+    metaPre = 1.0 / A.diagonal()
+    metaP = lambda x: metaPre*x
+    P = scipy.sparse.linalg.LinearOperator(A.shape, matvec=metaP)
+    self._prePractical = lambda y: np.concatenate(( \
+      scipy.sparse.linalg.cg(A, y[:y.size//3], tol=self._solvetol/10, M=P)[0], \
+      scipy.sparse.linalg.cg(A, y[y.size//3:2*y.size//3], tol=self._solvetol/10, M=P)[0], \
+      scipy.sparse.linalg.cg(A, y[2*y.size//3:], tol=self._solvetol/10, M=P)[0]))
+
     practical2D = lambda x: \
-      self._Q.dot(self._prePractical.dot(self._Q.transpose().dot(x)))
+      self._Q.dot(self._prePractical(self._Q.transpose().dot(x)))
 
     self._preconditioner = scipy.sparse.linalg.LinearOperator( \
       (self._Q.shape[0], self._Q.shape[0]), matvec=practical2D)