[firedrake] Matrix multiplication in bilinear form

[Anna Kalogirou]

Hi Lawrence,

Can I manually configure/multiply those matrices in Python and then 
solve the linear system Ax=b in UFL by writing solve(A, x, b), where b 
can be assembled using other known functions?

Anna.


On 26/10/15 11:18, Lawrence Mitchell wrote:
> -----BEGIN PGP SIGNED MESSAGE-----
> Hash: SHA1
>
> Hi Anna,
>
> On 23/10/15 18:54, Anna Kalogirou wrote:
>> Hi Lawrence,
>>
>> Could you please help me with finding the inverse of a matrix and
>> also performing matrix multiplications using PETSc? I tried (but
>> failed) to do it.
>>
>> I know computing the inverse is costly (I only need to do it once
>> in the code), but it is crucial in the wave-ship problem I am
>> considering to have a LHS matrix of the form
>>
>> (M_b*M^{-1}) * A * (M^{-1}*M_b)
>>
>> instead of just A_b, where _b denotes a block matrix assembled
>> locally in one part of the domain by using a heavyside function.
>> The mentioned matrices have the following forms: M_form = u*v*dx
>> Mb_form = step_b*u*v*dx A_form = inner(grad(u),grad(v))*dx Ab_form
>> = step_b*inner(grad(u),grad(v))*dx.
> Do you really need this matrix to be formed, or do you just need its
> action?  If you really need it to be formed, we can do it, but it
> might need more effort to work in parallel.  This PETSc FAQ gives you
> the idea:
>
> http://www.mcs.anl.gov/petsc/documentation/faq.html#invertmatrix
>
> effectively, you do something like:
>
> from firedrake.petsc import PETSc
>
> M = assemble(u*v*dx).M.handle
>
> tmp = PETSc.Mat().createDense(M.getSizes())I'm slightly confused
> though, because looking at the equation, you have terms:
>
> M^{-1} M_b
>
> But M_b is just M (except where it's masked out). So surely this is
> just equivalent to the identity matrix where the heaviside function is
> 1 and 0 otherwise?
>
> tmp.setUp()
> diag = tmp.createVecLeft()
> diag.set(1.0)
> tmp.setDiagonal(diag, addv=PETSc.InsertMode.INSERT_VALUES)
>
> Minv = tmp.duplicate()
>
> M.factorLU()
>
> M.matSolve(tmp, Minv)
>
> M_b = assemble(Mb_form).M.handle
>
> MinvMb = Minv.matMult(M_b)
>
> A = assemble(A_form).M.handle
>
> AMinvMb = A.matMult(MinvMb)
>
> ...
>
>
> If you only need the action, then we can build a shell matrix that
> applies this operation and you can use that in a linear solver.
> You'll then want to think about a spectrally equivalent operator that
> you /can/ invert to use as a preconditioner (maybe using a
> mass-lumped, or diagonal, inverse of the mass matrix).
>
>
> Lawrence
> -----BEGIN PGP SIGNATURE-----
> Version: GnuPG v2.0.22 (GNU/Linux)
>
> iQEcBAEBAgAGBQJWLgwJAAoJECOc1kQ8PEYv27oIAJ+1kxzR8VFSkBeVZkMLUoej
> V0OTjfzUaVTb0rMqHWedoKnkucPfNN78GFW+2xDf6WgbddTaFT5OK3bT1tzDplT8
> SE/us3oyLot5NGKfJT9zfqjT+x1s1zHw7WeXK5hAWnZf+JFVi1ZoSyFjthEV3I06
> vY222YNEe9yoZQPlcNiNeqUOFSBeW8IP2imH6EwdbsYqK/5EIccPUruv4inSCJ8F
> r18SvXDX7D7jh566fOCCzBsDkv5WfnKKk3YGjd04lkS+q1Lg2aBGw4VJ8hrAKp26
> 0vScLVy5lTqqz3uSDfPNxNieob0V+kKZWPZcFHDJGNzrmzOMfT0DD9pTPbV9Mxc=
> =+3ID
> -----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/