Dear all, Is it allowed to work with an index set (IS) with size 0? In particular, using such an index set in VecScatter and MatGetSubMatrix, I am creating the following objects: v : a vector of size n>0, w : a vector of size 0, A : a matrix of size n\times0. Then I do MatMultAdd (A, w, v, v); This actually seems to change the value of v, although I supposed it not to do so. Are such index sets allowed or not? If not, what about index sets having local size 0 on some processors? Best Regards, Tim -- Dr. Tim Kroeger CeVis -- Center of Complex Systems and Visualization University of Bremen [email protected] Universitaetsallee 29 [email protected] D-28359 Bremen Phone +49-421-218-7710 Germany Fax +49-421-218-4236
On Tue, Nov 9, 2010 at 12:11, Tim Kroeger <[email protected]>wrote:
Is it allowed to work with an index set (IS) with size 0?
Yes. In particular, using such an index set in VecScatter and MatGetSubMatrix, I
am creating the following objects:
v : a vector of size n>0, w : a vector of size 0, A : a matrix of size n\times0.
Then I do
MatMultAdd (A, w, v, v);
This actually seems to change the value of v, although I supposed it not to do so.
It sounds like you have a simple test case for this, could you send it? Jed
On Tue, 9 Nov 2010, Jed Brown wrote:
On Tue, Nov 9, 2010 at 12:11, Tim Kroeger <[email protected]> wrote: Is it allowed to work with an index set (IS) with size 0?
Yes.
Good to hear, thank you.
In particular, using such an index set in VecScatter and MatGetSubMatrix, I am creating the following objects:
v : a vector of size n>0, w : a vector of size 0, A : a matrix of size n\times0.
Then I do
MatMultAdd (A, w, v, v);
This actually seems to change the value of v, although I supposed it not to do so.
It sounds like you have a simple test case for this, could you send it?
Well, as you know, I'm using PETSc only via libMesh. Creating a standalone PETSc test case is non-trivial for me. But I'll keep to this and will either find out that I did something wrong or produce a test case for you. It might be imporant to know that the original matrix is a ShellMatrix. Best Regards, Tim -- Dr. Tim Kroeger CeVis -- Center of Complex Systems and Visualization University of Bremen [email protected] Universitaetsallee 29 [email protected] D-28359 Bremen Phone +49-421-218-7710 Germany Fax +49-421-218-4236
On Tue, Nov 9, 2010 at 12:49, Tim Kroeger <[email protected]>wrote:
Well, as you know, I'm using PETSc only via libMesh. Creating a standalone PETSc test case is non-trivial for me. But I'll keep to this and will either find out that I did something wrong or produce a test case for you.
It might be imporant to know that the original matrix is a ShellMatrix.
Yes, this is important. Check the implementation of MatMultAdd for the Shell. Jed
On Tue, 9 Nov 2010, Jed Brown wrote:
On Tue, Nov 9, 2010 at 12:49, Tim Kroeger <[email protected]> wrote:
It might be imporant to know that the original matrix is a ShellMatrix.
Yes, this is important. Check the implementation of MatMultAdd for the Shell.
Thank you, I found it now: The shell matrix does not implement MatMultAdd, only MatMult. Of course, it can't work then. Thank you very much! Best Regards, Tim -- Dr. Tim Kroeger CeVis -- Center of Complex Systems and Visualization University of Bremen [email protected] Universitaetsallee 29 [email protected] D-28359 Bremen Phone +49-421-218-7710 Germany Fax +49-421-218-4236
On Tue, Nov 9, 2010 at 12:56, Tim Kroeger <[email protected]>wrote:
Thank you, I found it now: The shell matrix does not implement MatMultAdd, only MatMult. Of course, it can't work then. Thank you very much!
Okay, but then it will error instead of silently modifying the vector. Jed
On Tue, 9 Nov 2010, Jed Brown wrote:
On Tue, Nov 9, 2010 at 12:56, Tim Kroeger <[email protected]> wrote: Thank you, I found it now: The shell matrix does not implement MatMultAdd, only MatMult. Of course, it can't work then. Thank you very much!
Okay, but then it will error instead of silently modifying the vector.
It didn't. And I did not forget CHKERRABORT() after the call to MatMultAdd(). Also, implementing MatMultAdd() for the shell did not resolve the problem. I had a look into src/mat/impls/submat/submat.c (of petsc-3.1-p4) and found that MatMultAdd_SubMatrix() internally calls MatMult() (line 155), not MatMultAdd(). Later in that method, some VecAYPX() is called (line 159), which might do the required addition. Can you please confirm whether (a) this is correct but explains why leaving MatMultAdd() of the shell undefined does not error; or (b) everything remains as you said before; or (c) this is a bug in that PETSc version? Best Regards, Tim -- Dr. Tim Kroeger CeVis -- Center of Complex Systems and Visualization University of Bremen [email protected] Universitaetsallee 29 [email protected] D-28359 Bremen Phone +49-421-218-7710 Germany Fax +49-421-218-4236
On Tue, Nov 9, 2010 at 14:41, Tim Kroeger <[email protected]>wrote:
I had a look into src/mat/impls/submat/submat.c (of petsc-3.1-p4) and found that MatMultAdd_SubMatrix() internally calls MatMult() (line 155), not MatMultAdd(). Later in that method, some VecAYPX() is called (line 159), which might do the required addition. Can you please confirm whether
Ah, I missed that you called MatGetSubMatrix on a MatShell in which you hadn't implemented MatGetSubMatrix, therefore you had a MatSubMatrix. Looking at MatMultAdd_SubMatrix, I see that the code is wrong if v2==v3. We don't have "weighted scatters", so I think it needs to get another work vector for a correct answer with "left scaling" and v2==v3. I'll fix this in 3.1 and -dev, but I might not get to it today. Jed
On Tue, 9 Nov 2010, Jed Brown wrote:
On Tue, Nov 9, 2010 at 14:41, Tim Kroeger <[email protected]> wrote: I had a look into src/mat/impls/submat/submat.c (of petsc-3.1-p4) and found that MatMultAdd_SubMatrix() internally calls MatMult() (line 155), not MatMultAdd(). Later in that method, some VecAYPX() is called (line 159), which might do the required addition. Can you please confirm whether
Ah, I missed that you called MatGetSubMatrix on a MatShell in which you hadn't implemented MatGetSubMatrix, therefore you had a MatSubMatrix.
Sorry for not pointing this our more clearly.
Looking at MatMultAdd_SubMatrix, I see that the code is wrong if v2==v3. We don't have "weighted scatters", so I think it needs to get another work vector for a correct answer with "left scaling" and v2==v3. I'll fix this in 3.1 and -dev, but I might not get to it today.
Okay, I implemented a workaround on user side (using a temporary vector and MatMult()). This seems to resolve the problem now. I'll leave it like this in the libMesh part, so that libMesh remains compatible with PETSc-3.1-p4. Nevertheless, if you fix this within a few days and launch PETSc-3.1-p5 with this fix (or just send me a patch), let me know, I can then test it. Best Regards, Tim -- Dr. Tim Kroeger CeVis -- Center of Complex Systems and Visualization University of Bremen [email protected] Universitaetsallee 29 [email protected] D-28359 Bremen Phone +49-421-218-7710 Germany Fax +49-421-218-4236
On Tue, Nov 9, 2010 at 15:43, Tim Kroeger <[email protected]>wrote:
Okay, I implemented a workaround on user side (using a temporary vector and MatMult()). This seems to resolve the problem now. I'll leave it like this in the libMesh part, so that libMesh remains compatible with PETSc-3.1-p4. Nevertheless, if you fix this within a few days and launch PETSc-3.1-p5 with this fix (or just send me a patch), let me know, I can then test it.
I fixed this, and also a bug with nontrivial MatDiagonalScale after MatGetSubMatrix. You can get it from hg, or in the next patch level (I don't know when that will be). Jed
On Wed, 10 Nov 2010, Jed Brown wrote:
On Tue, Nov 9, 2010 at 15:43, Tim Kroeger <[email protected]> wrote: Okay, I implemented a workaround on user side (using a temporary vector and MatMult()). This seems to resolve the problem now. I'll leave it like this in the libMesh part, so that libMesh remains compatible with PETSc-3.1-p4. Nevertheless, if you fix this within a few days and launch PETSc-3.1-p5 with this fix (or just send me a patch), let me know, I can then test it.
I fixed this, and also a bug with nontrivial MatDiagonalScale after MatGetSubMatrix. You can get it from hg, or in the next patch level (I don't know when that will be).
Yes, this fixes it. Thank you very much again. Best Regards, Tim -- Dr. Tim Kroeger CeVis -- Center of Complex Systems and Visualization University of Bremen [email protected] Universitaetsallee 29 [email protected] D-28359 Bremen Phone +49-421-218-7710 Germany Fax +49-421-218-4236
On Thu, 11 Nov 2010, Tim Kroeger wrote:
On Wed, 10 Nov 2010, Jed Brown wrote:
On Tue, Nov 9, 2010 at 15:43, Tim Kroeger <[email protected]> wrote: Okay, I implemented a workaround on user side (using a temporary vector and MatMult()). This seems to resolve the problem now. I'll leave it like this in the libMesh part, so that libMesh remains compatible with PETSc-3.1-p4. Nevertheless, if you fix this within a few days and launch PETSc-3.1-p5 with this fix (or just send me a patch), let me know, I can then test it.
I fixed this, and also a bug with nontrivial MatDiagonalScale after MatGetSubMatrix. You can get it from hg, or in the next patch level (I don't know when that will be).
Yes, this fixes it. Thank you very much again.
One more question: Does this bux and fix also apply to a submatrix of an MPIAIJ matrix? That is, can I safely use MatMultAdd() having v2==v3 for such a matrix in 3.1.0-p4 or not? (I'm having trouble with such a matrix now, but the error might be on my side.) Best Regards, Tim -- Dr. Tim Kroeger CeVis -- Center of Complex Systems and Visualization University of Bremen [email protected] Universitaetsallee 29 [email protected] D-28359 Bremen Phone +49-421-218-7710 Germany Fax +49-421-218-4236
On Thu, Nov 11, 2010 at 17:03, Tim Kroeger <[email protected]>wrote:
One more question: Does this bux and fix also apply to a submatrix of an MPIAIJ matrix? That is, can I safely use MatMultAdd() having v2==v3 for such a matrix in 3.1.0-p4 or not? (I'm having trouble with such a matrix now, but the error might be on my side.)
MPIAIJ implements MatGetSubMatrix, so it returns a new, efficient MatMPIAIJ matrix. This patch only affects matrix types that do not implement MatGetSubMatrix, and thus fall back on creating a MatSubMatrix (which just holds a reference to the original matrix and implements multiplication through scatters). So calling MatMultAdd using a matrix of any type other than MatSubMatrix is unrelated to my patch. Jed
participants (2)
-
Jed Brown -
Tim Kroeger