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