Hi Lawrence, In the email below you explained to me how to define a matrix free operator B. Is there a way to inspect the values of the matrix behind this operator? Best, Anna. On 20/10/15 12:30, Lawrence Mitchell wrote:
-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1
On 19/10/15 19:43, Anna Kalogirou wrote:
I am not familiar with this expression |Q1> <Q2|. What I really want is Q1*Q2^T, where Q1 and Q2 and vectors defined in my previous email. Actually I think it's easier if I just provide you with the actual equation I need to solve, please see attached. OK, so we're on the same page.
Note that the operator Q1*Q2^T is dense, so you don't want to form it. However, its action is cheap, since if you have the assembled vectors you just do:
Q1*Q2^T x == Q1 dot(Q2, x)
For the purposes of solving the system you want, I think what you'll want to do is use firedrake to assemble the pieces and then define a matrix-free operator to pass to PETSc (via petsc4py) to actually solve the problem. Finally, you'll want to think about how to precondition the system: B == A + Q1*Q2^T.
Fortunately, the perturbuation is rank-1 and so the Sherman-Morrison-Woodbury formula tells you how to compute a good preconditioner for B, given one for A:
B^{-1} = A^{-1} + (A^{-1} Q1 * Q2^T A^{-1}) / (1 + Q2^T A^{-1} Q1)
To hook this up, you'll need to know some petsc4py, as well as firedrake.
Assume the forms are already defined elsewhere.
First I'll build the operator B:
from firedrake.petsc import PETSc
class MatrixFreeB(object): def __init__(self, A, Q1, Q2): self.A = A self.Q1 = Q1 self.Q2 = Q2
def mult(self, mat, x, y): "Compute y = B*x" # y = Ax self.A.mult(x, y) # alpha = Q2^T x alpha = self.Q2.dot(x) # y = alpha Q1 + y y.axpy(alpha, self.Q1)
# Get reference to PETSc Matrix A = assemble(A_form).M.handle Q1f = assemble(Q1_form) Q2f = assemble(Q2_form)
# Get reference to PETSc Vectors for Q1 and Q2 # I take a copy here, so if Q1 and Q2 change, # you'll need to do something difference with Q1f.dat.vec_ro as v: Q1 = v.copy() A^{- with Q2f.dat.vec_ro as v: Q2 = v.copy()
B = PETSc.Mat().create()
B.setSizes(*A.getSizes()) B.setType(B.Type.PYTHON) # Indicate that B should use the MatrixFreeB class to # compute mult. B.setPythonContext(MatrixFreeB(A, Q1, Q2)) B.setUp()
# For small problems, we don't need a preconditioner at all, so let's # check this works:
solver = PETSc.KSP().create()
solver.setOperators(B)
opts = PETSc.Options() opts["pc_type"] = "none"
solver.setUp() solver.setFromOptions()
# Now let's solve the system.
rhs = assemble(rhs_form) solution = Function(result_space)
with rhs.dat.vec_ro as b: with solution.dat.vec as x: solver.solve(b, x)
Let's try and get this working on a small problem first and then move on to constructing the preconditioner as well. Note that if you have boundary conditions in your operator, we'll have to do a little more work (but not much).
Cheers,
Lawrence -----BEGIN PGP SIGNATURE----- Version: GnuPG v2.0.22 (GNU/Linux)
iQEcBAEBAgAGBQJWJiXEAAoJECOc1kQ8PEYvSGsH/2mDSfh9VkHb+JwJpMt6IEhs +Wx3SfqyfaztbEKay1QPfd50qU7Ojfx1Gv2j7bKFxNh1ppHBqRWo/awZmgbsm2wI aPbVxUL9v5Pj6294DiWVoB8MKFxdQvDcMxm5FNXKRof9M8aLHOOFuAX2fK+aTae7 5a4u/a7yduOh/2qZ19tylaCvDlDx1pBMMU8/mGI28ecmsBaGsBmLR5ObvOnn/LPH HSe/uvBPwautxD5VssgmveF6C5Nwqa/LkYRh5a8hKZhAE3r1PqQhVCiyLbFtEw34 +RdQhldLKKwGenFhF5fuii5v1F4hzctfrotoCkAojuf75UIFXlmdsOKD3eP7lPM= =fus0 -----END PGP SIGNATURE-----
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
-- Dr Anna Kalogirou Research Fellow School of Mathematics University of Leeds http://www1.maths.leeds.ac.uk/~matak/