Re: [petsc-users] Valgrind Errors with KSPSolve

2018-10-26 Thread Phil Tooley
Thanks Matt,

I got it figured out.  As with a lot of uninitialised value bugs it
started a long way away (in a scalar value) and propagated through. 
What is slightly bizarre is that valgrind didn't pick it up where I
define and first use it, even with --track-origins=yes.

Thank you for your help.

Phil


On 26/10/2018 15:01, Matthew Knepley wrote:
> On Fri, Oct 26, 2018 at 9:54 AM Phil Tooley
> mailto:phil.too...@sheffield.ac.uk>> wrote:
>
> Hi All,
>
> Running valgrind over my code reveals a huge number (it stops
> collecting
> after 1e6) of uninitialised value errors occuring the KSPSolve
> routine. 
> These don't occur with the ksp example scripts but I can't figure out
> what I am doing wrong in my code.  My code works as expected.  The
> same
> errors occur with both version 3.9.4 and 3.10.0.
>
> I am developing in c++ and using smart pointers to allow
> autodestruction
> of Petsc objects, but I don't think this should have any side
> effects. 
>
>
> From casual inspection, it sure seems like the matrix you are passing
> KSP has unitialized
> values.
>
> I used smart pointers for years, and work on
> packages that had smart
> pointers. In almost every case, we ended up ripping out the smart
> pointers and going with
> explicit destruction, the amount of work they save never outweighed
> the problems with
> detecting errors
>
>   Thanks,
>
>      Matt
>
>   Thanks,
>
>     Matt
>  
>
> I do the ksp setup as:
>
> KSP_unique m_ksp = create_unique_ksp();
>
> perr = KSPCreate(m_comm, m_ksp.get());CHKERRABORT(m_comm, perr);
>
> perr = KSPSetOperators(*m_ksp, *normmat, *normmat);CHKERRABORT(m_comm,
> perr);
> perr = KSPSetUp(*m_ksp);CHKERRABORT(m_comm, perr);
>
> perr = KSPSetFromOptions(*m_ksp);CHKERRABORT(m_comm, perr);
> perr = KSPSetUp(*m_ksp);CHKERRABORT(m_comm, perr);
>
> perr = KSPSolve(*m_ksp, *m_rhs, *m_delta);CHKERRABORT(m_comm, perr);
>
> Where normmat is a smart pointer to a previously created Mat
> object and
> m_rhs and m_delta are smart pointers to Vec objects which are
> initialised using MatCreateVecs on normmat.
>
> I am using MPICH to provide a (hopefully) valgrind clean MPI
> implementation.
>
> Any insights are appreciated.
>
> Many Thanks
>
> Phil Tooley
>
> Typical errors are as follows:
>
> ==7368== Conditional jump or move depends on uninitialised value(s)
> ==7368==    at 0x58894EE: PetscAbsScalar(double) (petscmath.h:389)
> ==7368==    by 0x5891702: MatLUFactorNumeric_SeqAIJ (aijfact.c:553)
> ==7368==    by 0x5406B83: MatLUFactorNumeric (matrix.c:3031)
> ==7368==    by 0x609F317: PCSetUp_ILU(_p_PC*) (ilu.c:176)
> ==7368==    by 0x6192D4F: PCSetUp (precon.c:923)
> ==7368==    by 0x62EAE56: KSPSetUp (itfunc.c:381)
> ==7368==    by 0x60190BB: PCSetUpOnBlocks_BJacobi_Singleblock(_p_PC*)
> (bjacobi.c:621)
> ==7368==    by 0x619352C: PCSetUpOnBlocks (precon.c:954)
> ==7368==    by 0x62E8F0A: KSPSetUpOnBlocks (itfunc.c:213)
> ==7368==    by 0x62ECCB1: KSPSolve (itfunc.c:613)
> ==7368==    by 0x10FCDD: Elastic::innerstep(double) (in
> /home/telemin/repos/pfire_petsc/bin/pfire)
> ==7368==    by 0x110600: Elastic::innerloop(int) (in
> /home/telemin/repos/pfire_petsc/bin/pfire)
> ==7368==
> ==7368== Conditional jump or move depends on uninitialised value(s)
> ==7368==    at 0x58894EE: PetscAbsScalar(double) (petscmath.h:389)
> ==7368==    by 0x588A75F: MatPivotCheck_none(_p_Mat*, _p_Mat*,
> MatFactorInfo const*, FactorShiftCtx*, int) (matimpl.h:704)
> ==7368==    by 0x588AD6E: MatPivotCheck(_p_Mat*, _p_Mat*,
> MatFactorInfo
> const*, FactorShiftCtx*, int) (matimpl.h:727)
> ==7368==    by 0x589177D: MatLUFactorNumeric_SeqAIJ (aijfact.c:558)
> ==7368==    by 0x5406B83: MatLUFactorNumeric (matrix.c:3031)
> ==7368==    by 0x609F317: PCSetUp_ILU(_p_PC*) (ilu.c:176)
> ==7368==    by 0x6192D4F: PCSetUp (precon.c:923)
> ==7368==    by 0x62EAE56: KSPSetUp (itfunc.c:381)
> ==7368==    by 0x60190BB: PCSetUpOnBlocks_BJacobi_Singleblock(_p_PC*)
> (bjacobi.c:621)
> ==7368==    by 0x619352C: PCSetUpOnBlocks (precon.c:954)
> ==7368==    by 0x62E8F0A: KSPSetUpOnBlocks (itfunc.c:213)
> ==7368==    by 0x62ECCB1: KSPSolve (itfunc.c:613)
> ==7368==
> ==7368== Conditional jump or move depends on uninitialised value(s)
> ==7368==    at 0x588A76D: MatPivotCheck_none(_p_Mat*, _p_Mat*,
> MatFactorInfo const*, FactorShiftCtx*, int) (matimpl.h:704)
> ==7368==    by 0x588AD6E: MatPivotCheck(_p_Mat*, 

[petsc-users] Valgrind Errors with KSPSolve

2018-10-26 Thread Phil Tooley
Hi All,

Running valgrind over my code reveals a huge number (it stops collecting
after 1e6) of uninitialised value errors occuring the KSPSolve routine. 
These don't occur with the ksp example scripts but I can't figure out
what I am doing wrong in my code.  My code works as expected.  The same
errors occur with both version 3.9.4 and 3.10.0.

I am developing in c++ and using smart pointers to allow autodestruction
of Petsc objects, but I don't think this should have any side effects. 

I do the ksp setup as:

KSP_unique m_ksp = create_unique_ksp();

perr = KSPCreate(m_comm, m_ksp.get());CHKERRABORT(m_comm, perr);

perr = KSPSetOperators(*m_ksp, *normmat, *normmat);CHKERRABORT(m_comm,
perr);
perr = KSPSetUp(*m_ksp);CHKERRABORT(m_comm, perr);

perr = KSPSetFromOptions(*m_ksp);CHKERRABORT(m_comm, perr);
perr = KSPSetUp(*m_ksp);CHKERRABORT(m_comm, perr);

perr = KSPSolve(*m_ksp, *m_rhs, *m_delta);CHKERRABORT(m_comm, perr);

Where normmat is a smart pointer to a previously created Mat object and
m_rhs and m_delta are smart pointers to Vec objects which are
initialised using MatCreateVecs on normmat.

I am using MPICH to provide a (hopefully) valgrind clean MPI implementation.

Any insights are appreciated.

Many Thanks

Phil Tooley

Typical errors are as follows:

==7368== Conditional jump or move depends on uninitialised value(s)
==7368==    at 0x58894EE: PetscAbsScalar(double) (petscmath.h:389)
==7368==    by 0x5891702: MatLUFactorNumeric_SeqAIJ (aijfact.c:553)
==7368==    by 0x5406B83: MatLUFactorNumeric (matrix.c:3031)
==7368==    by 0x609F317: PCSetUp_ILU(_p_PC*) (ilu.c:176)
==7368==    by 0x6192D4F: PCSetUp (precon.c:923)
==7368==    by 0x62EAE56: KSPSetUp (itfunc.c:381)
==7368==    by 0x60190BB: PCSetUpOnBlocks_BJacobi_Singleblock(_p_PC*)
(bjacobi.c:621)
==7368==    by 0x619352C: PCSetUpOnBlocks (precon.c:954)
==7368==    by 0x62E8F0A: KSPSetUpOnBlocks (itfunc.c:213)
==7368==    by 0x62ECCB1: KSPSolve (itfunc.c:613)
==7368==    by 0x10FCDD: Elastic::innerstep(double) (in
/home/telemin/repos/pfire_petsc/bin/pfire)
==7368==    by 0x110600: Elastic::innerloop(int) (in
/home/telemin/repos/pfire_petsc/bin/pfire)
==7368==
==7368== Conditional jump or move depends on uninitialised value(s)
==7368==    at 0x58894EE: PetscAbsScalar(double) (petscmath.h:389)
==7368==    by 0x588A75F: MatPivotCheck_none(_p_Mat*, _p_Mat*,
MatFactorInfo const*, FactorShiftCtx*, int) (matimpl.h:704)
==7368==    by 0x588AD6E: MatPivotCheck(_p_Mat*, _p_Mat*, MatFactorInfo
const*, FactorShiftCtx*, int) (matimpl.h:727)
==7368==    by 0x589177D: MatLUFactorNumeric_SeqAIJ (aijfact.c:558)
==7368==    by 0x5406B83: MatLUFactorNumeric (matrix.c:3031)
==7368==    by 0x609F317: PCSetUp_ILU(_p_PC*) (ilu.c:176)
==7368==    by 0x6192D4F: PCSetUp (precon.c:923)
==7368==    by 0x62EAE56: KSPSetUp (itfunc.c:381)
==7368==    by 0x60190BB: PCSetUpOnBlocks_BJacobi_Singleblock(_p_PC*)
(bjacobi.c:621)
==7368==    by 0x619352C: PCSetUpOnBlocks (precon.c:954)
==7368==    by 0x62E8F0A: KSPSetUpOnBlocks (itfunc.c:213)
==7368==    by 0x62ECCB1: KSPSolve (itfunc.c:613)
==7368==
==7368== Conditional jump or move depends on uninitialised value(s)
==7368==    at 0x588A76D: MatPivotCheck_none(_p_Mat*, _p_Mat*,
MatFactorInfo const*, FactorShiftCtx*, int) (matimpl.h:704)
==7368==    by 0x588AD6E: MatPivotCheck(_p_Mat*, _p_Mat*, MatFactorInfo
const*, FactorShiftCtx*, int) (matimpl.h:727)
==7368==    by 0x589177D: MatLUFactorNumeric_SeqAIJ (aijfact.c:558)
==7368==    by 0x5406B83: MatLUFactorNumeric (matrix.c:3031)
==7368==    by 0x609F317: PCSetUp_ILU(_p_PC*) (ilu.c:176)
==7368==    by 0x6192D4F: PCSetUp (precon.c:923)
==7368==    by 0x62EAE56: KSPSetUp (itfunc.c:381)
==7368==    by 0x60190BB: PCSetUpOnBlocks_BJacobi_Singleblock(_p_PC*)
(bjacobi.c:621)
==7368==    by 0x619352C: PCSetUpOnBlocks (precon.c:954)
==7368==    by 0x62E8F0A: KSPSetUpOnBlocks (itfunc.c:213)
==7368==    by 0x62ECCB1: KSPSolve (itfunc.c:613)
==7368==    by 0x10FCDD: Elastic::innerstep(double) (in
/home/telemin/repos/pfire_petsc/bin/pfire)
==7368==
==7368== Conditional jump or move depends on uninitialised value(s)
==7368==    at 0x58894EE: PetscAbsScalar(double) (petscmath.h:389)
==7368==    by 0x5889528: PetscIsNanScalar(double) (petscmath.h:674)
==7368==    by 0x588A784: MatPivotCheck_none(_p_Mat*, _p_Mat*,
MatFactorInfo const*, FactorShiftCtx*, int) (matimpl.h:704)
==7368==    by 0x588AD6E: MatPivotCheck(_p_Mat*, _p_Mat*, MatFactorInfo
const*, FactorShiftCtx*, int) (matimpl.h:727)
==7368==    by 0x589177D: MatLUFactorNumeric_SeqAIJ (aijfact.c:558)
==7368==    by 0x5406B83: MatLUFactorNumeric (matrix.c:3031)
==7368==    by 0x609F317: PCSetUp_ILU(_p_PC*) (ilu.c:176)
==7368==    by 0x6192D4F: PCSetUp (precon.c:923)
==7368==    by 0x62EAE56: KSPSetUp (itfunc.c:381)
==7368==    by 0x60190BB: PCSetUpOnBlocks_BJacobi_Singleblock(_p_PC*)
(bjacobi.c:621)
==7368==    by 0x619352C: PCSetUpOnBlocks (precon.c:954)
==7368==    by 0x62E8F0A: KSPSetUpOnBlocks (itfunc.c:213)
==7368==
==7368== Conditional

[petsc-users] Better understanding VecScatter

2018-10-04 Thread Phil Tooley
Hi all,

I am trying to understand how to create a custom scatter from petsc
ordering in a local vector to natural ordering in a global vector.

I have a 3D DMDA and local vector containing field data and am
calculating the x y and z gradients into their own local vectors.  I
then need to scatter these gradients to different parts of a single
vector so they are arranged [x_1 .. x_n y_1 .. y_n z_1 .. z_n] with
natural ordering.  (This is for use in a custom numerical scheme that I
have been asked to parallelise)

I can of course do this using multiple scatters but I need to perform
this operation regularly and would like to express it as a single
scatter operation if possible.

I know I need to get the relevant AO from my DMDA and use it to
construct appropriate ISs to build the scatter context but I am unsure
how exactly to go about this.

Can someone point me in the right direction please?

Many Thanks

Phil

-- 
Phil Tooley
Research Software Engineering
University of Sheffield



Re: [petsc-users] DMDA with dimension of size 1

2018-10-01 Thread Phil Tooley
Hi Patrick,

You are spot on, I was copying in a vector and not taking care to pad
with 1s to length 3.

Thanks

Phil

On 01/10/2018 15:57, Patrick Sanan wrote:
> Whoops, sent that patch too fast (forgot to update SETERRQ3 to
> SETERRQ4). Updated patch for maint attached.
>
> Am Mo., 1. Okt. 2018 um 16:55 Uhr schrieb Patrick Sanan
> mailto:patrick.sa...@gmail.com>>:
>
> Meshes of size 1 should work.
>
> Looks like there is a bug in the code to produce this error
> message (patch for maint attached). It's not outputting the "P"
> (size in the 3rd dimension) component.
>
> This is just speculation without a full error message or code to
> reproduce, but perhaps the size in the third dimension is an
> uninitialized value which is triggering this warning with a 10 x
> 10 x (huge garbage) x 1 (dof) mesh.
>
> Am Mo., 1. Okt. 2018 um 15:59 Uhr schrieb Phil Tooley
> mailto:phil.too...@sheffield.ac.uk>>:
>
> Hi all,
>
> Is it valid to have a DMDA with one of the dimensions of size
> 1.  I was
> hoping to avoid having to write explicit logic to handle the
> case that I
> am working on a 2D rather than a 3D image for my current
> application.  
> Unfortunately when I try to construct such a DMDA I get an error:
>
> [0]PETSC ERROR: - Error Message
> --
> [0]PETSC ERROR: Overflow in integer operation:
> [0]PETSC ERROR: Mesh of 10 by 10 by 1 (dof) is too large for
> 32 bit indices
>
>     Is there a workaround other than "write everything twice"?
>
>     Thanks
>
> Phil
>
> -- 
> Phil Tooley
> Research Software Engineering
> University of Sheffield
>
-- 
Phil Tooley
Research Software Engineering
University of Sheffield



[petsc-users] DMDA with dimension of size 1

2018-10-01 Thread Phil Tooley
Hi all,

Is it valid to have a DMDA with one of the dimensions of size 1.  I was
hoping to avoid having to write explicit logic to handle the case that I
am working on a 2D rather than a 3D image for my current application.  
Unfortunately when I try to construct such a DMDA I get an error:

[0]PETSC ERROR: - Error Message
--
[0]PETSC ERROR: Overflow in integer operation:
[0]PETSC ERROR: Mesh of 10 by 10 by 1 (dof) is too large for 32 bit indices

Is there a workaround other than "write everything twice"?

Thanks

Phil

-- 
Phil Tooley
Research Software Engineering
University of Sheffield



Re: [petsc-users] Checking if a vector is a localvector of a given DMDA

2018-09-25 Thread Phil Tooley
Thanks both,

I now have what I need.  For now I am checking that the vector I am
passed has the same local size, global size, and Comm as the vector
provided by DMGetLocalVector, mostly because I already have a
compatibility check function written.  (I assume this requires a malloc
and free behind the scenes)

At some point I will likely change to explicitly checking for comm size
of one and appropriate global and local sizes based on the DMDA
properties instead, for now I want to get to an alpha version I can let
people play with.

Phil


On 25/09/18 13:07, Dave May wrote:
>
>
> On Tue, 25 Sep 2018 at 13:20, Matthew Knepley  <mailto:knep...@gmail.com>> wrote:
>
> On Tue, Sep 25, 2018 at 7:03 AM Dave May  <mailto:dave.mayhe...@gmail.com>> wrote:
>
> On Tue, 25 Sep 2018 at 11:49, Phil Tooley
>  <mailto:phil.too...@sheffield.ac.uk>> wrote:
>
> Hi all,
>
> Given a vector I know I can get an associated DM (if there
> is one) by
> calling VecGetDM, but I need to also be able to check that
>
> a) the vector is the localvector of that DM rather than
> the global
>
>
> Given the vector, you can check the communicator size via
> PetscObjectGetComm()
>
>  
> 
> https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/PetscObjectGetComm.html
> and then MPI_Comm_size()
> If the comm size 1, it is local vector.
>
>
> In serial, both vectors have comm size 1.
>
>
> Right - and the local and global sizes are the same.
>
>  My point was to check the comm size first. If it's 1 then you have a
> candidate for a local vector. Then you'd check the vec global size
> matches the dmda local size. If the commsize is anything other than 1
> then it cannot be a local vector 
>
>
>    Matt
>  
>
> You can check the size matches your local DMDA space by using
> DMDAGetGhostCorners()
>
> 
> https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMDA/DMDAGetGhostCorners.html
>
> and return the quantities m, n, and p.
>
> You also need to use  DMDAGetInfo()
>
> 
> https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMDA/DMDAGetInfo.html
>
> The important quantity you want returned is "dof"
>
> If m x n x p x dof matches the number returned by VecGetSize()
> (assuming you know the vector is sequential) then you know the
> local space will fit within your vector.
>
>  
>
>
> b) the DM is a DMDA rather than some other subclass
>
>
> See Matt's answer 
>  
>
>
>     I can't immediately see routines that do what I need, but
> I am likely
> missing something obvious. Is there a way to achieve the
> above?
>
> Thanks
>
> Phil
>
> -- 
> Phil Tooley
> Research Software Engineering
> University of Sheffield
>
>
>
> -- 
>     What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to
> which their experiments lead.
> -- Norbert Wiener
>
> https://www.cse.buffalo.edu/~knepley/
> <http://www.cse.buffalo.edu/%7Eknepley/>
>

-- 
Phil Tooley
Research Software Engineering
University of Sheffield



Re: [petsc-users] Checking if a vector is a localvector of a given DMDA

2018-09-25 Thread Phil Tooley
Hi Matt,

Thanks for the rapid response.  The reason I want to ensure I have a
DMDA is that this is a custom gradient function, the results of which
are meaningless for anything but a regular grid. I would rather raise an
exception upfront with the user than allow the program to continue with
garbage data.

Phil

On 25/09/18 11:55, Matthew Knepley wrote:
> On Tue, Sep 25, 2018 at 6:49 AM Phil Tooley
> mailto:phil.too...@sheffield.ac.uk>> wrote:
>
> Hi all,
>
> Given a vector I know I can get an associated DM (if there is one) by
> calling VecGetDM, but I need to also be able to check that
>
> a) the vector is the localvector of that DM rather than the global
>
>
> Get a local vector from the DM and check the sizes.
>  
>
> b) the DM is a DMDA rather than some other subclass
>
>
> You can do this (PetscObjectTypeCompare), but it is not recommended.
> Why would
> you want to do that?
>
>   Thanks,
>
>      Matt
>  
>
> I can't immediately see routines that do what I need, but I am likely
> missing something obvious. Is there a way to achieve the above?
>
> Thanks
>
> Phil
>
> -- 
> Phil Tooley
> Research Software Engineering
> University of Sheffield
>
>
>
> -- 
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which
> their experiments lead.
> -- Norbert Wiener
>
> https://www.cse.buffalo.edu/~knepley/
> <http://www.cse.buffalo.edu/%7Eknepley/>

-- 
Phil Tooley
Research Software Engineering
University of Sheffield



[petsc-users] Checking if a vector is a localvector of a given DMDA

2018-09-25 Thread Phil Tooley
Hi all,

Given a vector I know I can get an associated DM (if there is one) by
calling VecGetDM, but I need to also be able to check that

a) the vector is the localvector of that DM rather than the global

b) the DM is a DMDA rather than some other subclass

I can't immediately see routines that do what I need, but I am likely
missing something obvious. Is there a way to achieve the above?

Thanks

Phil

-- 
Phil Tooley
Research Software Engineering
University of Sheffield



[petsc-users] DMCreateXXXVector vs VecDuplicate

2018-09-24 Thread Phil Tooley
Hi all,

Can I get come clarification regarding  DMCreateXXXVector vs
VecDuplicate. The manual says that I can obtain additional vectors with
the latter, does this mean that I *must* only call DMCreateXXXVector
once each for global and local?

I ask because I want multiple instances of an image class all with the
same distributed layout. It makes life much tidier if I can just pass
the DMDA rather than also having to pass vectors to be copied.

Cheers


Phil

-- 
Phil Tooley
Research Software Engineering
University of Sheffield



Re: [petsc-users] Indexing/using a 3D DMDA as a 1D vector

2018-09-12 Thread Phil Tooley
Ok, So I have figured this out a bit better now, facepalm moment as I
finally realised that this is exactly what the DMDA exists to do...  I
need to create a vector using the DMCreate*Vector functions and can then
insert values using the assembly functions. 

Then, if I understand correctly, in order to have the state I need to
multiply with my operator I should make use of the
DMDAGlobalToNaturalBegin() and End() functions to ensure the correct
ordering.  Do I then need to ensure I have reverted it to global
ordering before I try to apply finite difference methods again?

Thanks

P


On 12/09/18 14:27, Phil Tooley wrote:
> I will preface this by saying I am new to PETSc and am still trying to
> get my head around all of the layout mapping that is done.  That means I
> may well have fundamentally misunderstood something, but hopefully
> someone will be able to to put me right.
>
> In my application I have some 3D pixel data which I want to manipulate
> using finite difference methods and then transform by viewing as a 1-D
> vector and multiplying by a large sparse matrix operator.
>
> I would assume that the correct way to do this is by creating a DMDA to
> hold the image data and ghosting appropriately to apply my finite
> difference operations.  Then I had hoped that I could use some form of
> application ordering to allow viewing the data as a vector that can be
> multiplied with my operator matrix.  This is where I have come unstuck,
> I may just be missing something obivous but I can't figure out how to do
> this.  Can anyone point me in the correct direction please?
>
> Many Thanks
>

-- 
Phil Tooley
Research Software Engineering
University of Sheffield



[petsc-users] Indexing/using a 3D DMDA as a 1D vector

2018-09-12 Thread Phil Tooley
I will preface this by saying I am new to PETSc and am still trying to
get my head around all of the layout mapping that is done.  That means I
may well have fundamentally misunderstood something, but hopefully
someone will be able to to put me right.

In my application I have some 3D pixel data which I want to manipulate
using finite difference methods and then transform by viewing as a 1-D
vector and multiplying by a large sparse matrix operator.

I would assume that the correct way to do this is by creating a DMDA to
hold the image data and ghosting appropriately to apply my finite
difference operations.  Then I had hoped that I could use some form of
application ordering to allow viewing the data as a vector that can be
multiplied with my operator matrix.  This is where I have come unstuck,
I may just be missing something obivous but I can't figure out how to do
this.  Can anyone point me in the correct direction please?

Many Thanks

-- 
Phil Tooley
Research Software Engineering
University of Sheffield