[firedrake] Matrix multiplication in bilinear form

[Anna Kalogirou]

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 at 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/