Re: [firedrake] Matrix multiplication in bilinear form
On Wed, 21 Oct 2015 at 14:59 Mitchell, Lawrence < 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 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 https://mailman.ic.ac.uk/mailman/listinfo/firedrake
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
-----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-----
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@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/
-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1 On 11/11/15 11:49, Anna Kalogirou wrote:
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?
No, A must be a firedrake Matrix for solve(A, x, b) to work. If it is a PETSc Mat, you will need to get the PETSc Vecs from x and b and use a PETSc KSP to solve the system. If you're working in serial and have A as a numpy matrix, you'll need to use numpy.linalg.solve having obtained the data arrays defining x and b. Lawrence -----BEGIN PGP SIGNATURE----- Version: GnuPG v2.0.22 (GNU/Linux) iQEcBAEBAgAGBQJWQy1IAAoJECOc1kQ8PEYvV9EIAI2/gnMn/nfeNIvVZYTBqFgG wGJUhfhdGoe1aHPzU47BUtYLKusj9ai1Zz9hNBUa6E/nFAHO1aqvVj15EWGqt6Yf K7qp5hsMw+vvhCPd17vlHHyAtimFduLegtEzRQTHeHVJqFoX5e4+cY/AhMUEVdO7 MzCIndac7QSTIHarJXIPgNZTQK5BOaZarj6liwC8UOETkGgGcDuiY1QamzyXqSwn nBd80EcQsFlVquE1WV3lsR1QTA2Ft3JwVL/SDbFUhzuzIDS00Y4vN4ZcFWaZ2AbY gXfAGPleomblmm7cujeEiaSsQFKM2w2tJFD5sb0GtzRF6HcJHFcgX+t7N7i+H3c= =9gWg -----END PGP SIGNATURE-----
Also, I don't really understand what the difference is between forming the matrix or taking its action... 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@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/
-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1 On 11/11/15 12:02, Anna Kalogirou wrote:
Also, I don't really understand what the difference is between forming the matrix or taking its action...
So if you form the matrix, then you have A, and you can do things like explicitly compute its inverse (using LU, say). If you only give me the action, you don't have A, you just give me a recipe for computing Ax if I give you x. In this case, I don't know what the entries of A are (so I can't explicitly compute its inverse). Lawrence -----BEGIN PGP SIGNATURE----- Version: GnuPG v2.0.22 (GNU/Linux) iQEcBAEBAgAGBQJWQzB5AAoJECOc1kQ8PEYvugEIAPEyk/4lVpt1f0G/NSnoPOoX DK4FD+Kjydg+h2sKWx8U0wyIJEvh442vQ3wgtWj+iNVrt9MG8i2dMADUZugvYnc0 IUqLLE8AvHBkLY0++8qe2MPw4tCeX/r09OB3hMXo/39R//KjK3RnRGd0WNe7xncn aJC/1VrfJqxmJvyzYFps5bBIloTA+Gz6C3XZyr44k0xpjQRFnPBgIGp0AXBtVmei cXA1k/3F4TUrA1jFjNC1cNODtrI4V1+Hi9Amgs6XWj5Prk/370dC3sGDGeAO0x2i UH7LVPgsdFo6vsAF0TjOa8zQWEhvw7Iuhx//CtArIO3/dGen3gVQ0CMhVALntSY= =OLyo -----END PGP SIGNATURE-----
participants (3)
-
Anna Kalogirou
-
David Ham
-
Lawrence Mitchell