From ac76e58c33d580bac538b8d0a9808ae58b2eee1b Mon Sep 17 00:00:00 2001
From: Carl-Martin Pfeiler <carl-martin.pfeiler@asc.tuwien.ac.at>
Date: Wed, 20 Mar 2019 17:43:18 +0100
Subject: [PATCH] stationary precond rewritten with CG
---
.../_details/interfaces/tpsPreconditioner.py | 38 +++++--------------
1 file changed, 9 insertions(+), 29 deletions(-)
diff --git a/integrators/_details/interfaces/tpsPreconditioner.py b/integrators/_details/interfaces/tpsPreconditioner.py
index 81b4e3a..4328fe6 100644
--- a/integrators/_details/interfaces/tpsPreconditioner.py
+++ b/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)
--
GitLab