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. Thanks a lot and enjoy your weekend! Best, Anna. On 21/10/15 15:03, David Ham wrote:
On Wed, 21 Oct 2015 at 14:59 Mitchell, Lawrence <lawrence.mitchell@imperial.ac.uk <mailto:lawrence.mitchell@imperial.ac.uk>> wrote:
-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1
On 21/10/15 14:32, Anna Kalogirou wrote: > But the matrix A has the same form as C, namely it has nonzero > variables in a block submatrix only.
Ah, oops, I see.
So because you're using an iterative solver without any preconditioning, the operator is never being inverted directly (hence you never get zero pivot errors in an LU factorisation). Your operator just has a gigantic null space. If the right hand side is consistent you might get lucky in that the Krylov solver can just explore the space orthogonal to the null space (and you end up with a solution in that subspace).
I guess what you really want to do is only solve the problem in the region where the heaviside function is nonzero, but we don't currently have proper support for this. I guess if you want to be able to invert the system somehow, we can add the identity matrix to the zero block and (presumably) zero the right hand side in the appropriate rows.
Note that this is equivalent to prescribing a dirichlet BC everywhere in the 0 block.
Lawrence
> > On 21/10/15 14:26, Lawrence Mitchell wrote: On 21/10/15 14:14, Anna > Kalogirou wrote: >>>> Dear Lawrence, >>>> >>>> I believe the code (attached) sets up the >>>> matrices/vectrors/operators correctly. However, I have one >>>> question: Due to the use of a heavyside function, the >>>> operator B on the LHS is nonzero in a block NxN, say, in the >>>> bottom right corner. For that reason, when I solved the >>>> linear system in Matlab I only considered these nonzero >>>> values in the block, because otherwise the LHS would not be >>>> invertible. How come and it doesn't complain here? > So it's clear that the matrix C = Q1*Q2^T is not invertible, > however, presumably you're solving for: > > B = A + C > > where A is full rank and invertible. So I don't immediately see > why there would be a problem in trying to invert B which is what > the solver does. > > Lawrence > >> >> _______________________________________________ firedrake mailing >> list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> >> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >
-----BEGIN PGP SIGNATURE----- Version: GnuPG v2.0.22 (GNU/Linux)
iQEcBAEBAgAGBQJWJ5o0AAoJECOc1kQ8PEYvxpAH/AmnOnhCRU+gVH5Yxjf76jhE evR42Oie1wLnZPjgP1CeYpbb/QaMk9i6X3TMWN7fgTIheoIIQMusu7a1rl7oMCdN RzUcVyXACNd4lcKT/SpSm2AcQxHkf1rLGzXsU1R8OoCqj+gFPWuLkCX0QX1gzWWe OGd47c1v+QbFzJgTUbI/2ivd3ACYnBGkHHuDapFYQYJvCyCF0aAdDrsDWu5EruIJ Xt7wUCbtcauEFih3WnDKDO1e7GwERyTbFhrnHWwqqcjIbDRFl6ENr6mrkduN1r/3 78d+4YdUabyx1Yq9C8P0srAsOAshpYYcTViVAemLcMC8LoCqILW3o/n9GnW7ntU= =L2Fs -----END PGP SIGNATURE-----
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake