Re: [petsc-users] Question about rank of matrix

[Stefano Zampini]

On Fri, Feb 17, 2023, 10:43 user_gong Kim  wrote:

> Hello,
>
> I have a question about rank of matrix.
> At the problem
> Au = b,
>
> In my case, sometimes global matrix A is not full rank.
> In this case, the global matrix A is more likely to be singular, and if it
> becomes singular, the problem cannot be solved even in the case of the
> direct solver.
> I haven't solved the problem with an iterative solver yet, but I would
> like to ask someone who has experienced this kind of problem.
>
> 1. If it is not full rank, is there a numerical technique to solve it by
> catching rows and columns with empty ranks in advance?
>
> 2.If anyone has solved it in a different way than the above numerical
> analysis method, please tell me your experience.
>
> Thanks,
> Hyung Kim
>

My experience with this is usually associated to reading a book and find
the solution I'm looking for.

>
>
>

[petsc-users] Question about rank of matrix

[user_gong Kim]

Hello,

I have a question about rank of matrix.
At the problem
Au = b,

In my case, sometimes global matrix A is not full rank.
In this case, the global matrix A is more likely to be singular, and if it
becomes singular, the problem cannot be solved even in the case of the
direct solver.
I haven't solved the problem with an iterative solver yet, but I would like
to ask someone who has experienced this kind of problem.

1. If it is not full rank, is there a numerical technique to solve it by
catching rows and columns with empty ranks in advance?

2.If anyone has solved it in a different way than the above numerical
analysis method, please tell me your experience.

Thanks,
Hyung Kim

Re: [petsc-users] PetscViewer with 64bit

[Jed Brown]

Okay, this works now. I'm pretty sure I tested this long ago with connectivity 
using Int64 and found that didn't work. That may have been ancient history, but 
I'm hesitant to revamp to match PetscInt. If doing that, it would require 
changing the signature of DMPlexGetVTKConnectivity to use PetscInt instead of 
PetscVTKInt. I'm already underwater and don't have the stamina to test it, but 
this MR will get you going for problems in which individual parts don't have 
more than 2B dofs.

https://gitlab.com/petsc/petsc/-/merge_requests/6081/diffs?commit_id=27ba695b8b62ee2bef0e5776c33883276a2a1735

Mike Michell  writes:

> Jed,
>
> It does not work for me even with the reproducer case with the small 2D
> square mesh. Can you run the case that I sent and open the created
> "sol.vtu" file with paraview?
>
> Thanks,
>
>
>> Thanks, Dave.
>>
>> Mike, can you test that this branch works with your large problems? I
>> tested that .vtu works in parallel for small problems, where works = loads
>> correctly in Paraview and VisIt.
>>
>> https://gitlab.com/petsc/petsc/-/merge_requests/6081
>>
>> Dave May  writes:
>>
>> > On Tue 14. Feb 2023 at 21:27, Jed Brown  wrote:
>> >
>> >> Dave May  writes:
>> >>
>> >> > On Tue 14. Feb 2023 at 17:17, Jed Brown  wrote:
>> >> >
>> >> >> Can you share a reproducer? I think I recall the format requiring
>> >> certain
>> >> >> things to be Int32.
>> >> >
>> >> >
>> >> > By default, the byte offset used with the appended data format is
>> >> UInt32. I
>> >> > believe that’s where the sizeof(int) is coming from. This default is
>> >> > annoying as it limits the total size of your appended data to be < 3
>> GB.
>> >> > That said, in the opening of the paraview file you can add this
>> attribute
>> >> >
>> >> > header_type="UInt64"
>> >>
>> >> You mean in the header of the .vtu?
>> >
>> >
>> > Yeah, within the open VTKFile tag.
>> > Like this
>> > < VTKFile type=“xxx”,  byte_order="LittleEndian" header_type="UInt64" >
>> >
>> > Do you happen to have an example or pointers to docs describing this
>> >> feature?
>> >
>> >
>> > Example yes - will send it to you tomorrow. Docs… not really. Only stuff
>> > like this
>> >
>> >
>> https://kitware.github.io/paraview-docs/latest/python/paraview.simple.XMLPStructuredGridWriter.html
>> >
>> >
>> >
>> https://kitware.github.io/paraview-docs/v5.8.0/python/paraview.simple.XMLMultiBlockDataWriter.html
>> >
>> > All the writers seem to support it.
>> >
>> >
>> > Can we always do this?
>> >
>> >
>> > Yep!
>> >
>> >
>> > It isn't mentioned in these:
>> >>
>> >> https://vtk.org/wp-content/uploads/2015/04/file-formats.pdf   (PDF was
>> >> created in 2003)
>> >>
>> >>
>> https://kitware.github.io/vtk-examples/site/VTKFileFormats/#xml-file-formats
>> >>
>> >
>> > Yes I know. I’ve tied myself in knots for years because the of the
>> > assumption that the offset had to be an int.
>> >
>> > Credit for the discovery goes to Carsten Uphoff. He showed me this.
>> >
>> > Cheers,
>> > Dave
>> >
>> >
>> >
>> >> > then the size of the offset is now UInt64 and now large files can be
>> >> > finally written.
>> >> >
>> >> >
>> >> > Cheers,
>> >> > Dave
>> >> >
>> >> >
>> >> >
>> >> >
>> >> >>
>> >> >> Mike Michell  writes:
>> >> >>
>> >> >> > Thanks for the note.
>> >> >> > I understood that PETSc calculates the offsets for me through
>> >> "boffset"
>> >> >> > variable in plexvtu.c file. Please correct me if it is wrong.
>> >> >> >
>> >> >> > If plexvtu.c has a bug, it could be around "write file header"
>> part in
>> >> >> > which the boffset is also computed. Is this correct? I am not using
>> >> >> complex
>> >> >> > number.
>> >> >> > There are several mixed parts among "Int32, UInt8, PetscInt_FMT,
>> >> >> > PetscInt64_FMT" in writing the header.
>> >> >> >
>> >> >> > Which combination of those flags is correct for 64bit indices? I am
>> >> gonna
>> >> >> > modify plexvtu.c file with "#if defined(PETSC_USE_64BIT_INDICES)"
>> >> >> > statement, but I do not know what is the correct form of the header
>> >> flag
>> >> >> > for 64bit indices.
>> >> >> >
>> >> >> > It is also confusing to me:
>> >> >> > boffset += gpiece[r].ncells * sizeof(PetscInt) + sizeof(int);
>> >> >> > How is sizeof(PetscInt) different from sizeof(int)?
>> >> >> >
>> >> >> > Thanks,
>> >> >> > Mike
>> >> >> >
>> >> >> >
>> >> >> >> On Tue, Feb 14, 2023 at 11:45 AM Mike Michell <
>> mi.mike1...@gmail.com
>> >> >
>> >> >> >> wrote:
>> >> >> >>
>> >> >> >>> I was trying to modify the header flags from "Int32" to "Int64",
>> but
>> >> >> the
>> >> >> >>> problem was not resolved. Could I get any additional comments?
>> >> >> >>>
>> >> >> >>
>> >> >> >> The calculated offsets are not correct I think.
>> >> >> >>
>> >> >> >>   Matt
>> >> >> >>
>> >> >> >>
>> >> >> >>> Thanks,
>> >> >> >>> Mike
>> >> >> >>>
>> >> >> >>>
>> >> >>  Thanks for the comments.
>> >> >>  To be precise on the question, the entire part of the header of
>> the
>> >> >> .vtu
>> >> >>  file is attached:
>> >> >> 
>> 

Re: [petsc-users] PetscViewer with 64bit

[Mike Michell]

Jed,

It does not work for me even with the reproducer case with the small 2D
square mesh. Can you run the case that I sent and open the created
"sol.vtu" file with paraview?

Thanks,


> Thanks, Dave.
>
> Mike, can you test that this branch works with your large problems? I
> tested that .vtu works in parallel for small problems, where works = loads
> correctly in Paraview and VisIt.
>
> https://gitlab.com/petsc/petsc/-/merge_requests/6081
>
> Dave May  writes:
>
> > On Tue 14. Feb 2023 at 21:27, Jed Brown  wrote:
> >
> >> Dave May  writes:
> >>
> >> > On Tue 14. Feb 2023 at 17:17, Jed Brown  wrote:
> >> >
> >> >> Can you share a reproducer? I think I recall the format requiring
> >> certain
> >> >> things to be Int32.
> >> >
> >> >
> >> > By default, the byte offset used with the appended data format is
> >> UInt32. I
> >> > believe that’s where the sizeof(int) is coming from. This default is
> >> > annoying as it limits the total size of your appended data to be < 3
> GB.
> >> > That said, in the opening of the paraview file you can add this
> attribute
> >> >
> >> > header_type="UInt64"
> >>
> >> You mean in the header of the .vtu?
> >
> >
> > Yeah, within the open VTKFile tag.
> > Like this
> > < VTKFile type=“xxx”,  byte_order="LittleEndian" header_type="UInt64" >
> >
> > Do you happen to have an example or pointers to docs describing this
> >> feature?
> >
> >
> > Example yes - will send it to you tomorrow. Docs… not really. Only stuff
> > like this
> >
> >
> https://kitware.github.io/paraview-docs/latest/python/paraview.simple.XMLPStructuredGridWriter.html
> >
> >
> >
> https://kitware.github.io/paraview-docs/v5.8.0/python/paraview.simple.XMLMultiBlockDataWriter.html
> >
> > All the writers seem to support it.
> >
> >
> > Can we always do this?
> >
> >
> > Yep!
> >
> >
> > It isn't mentioned in these:
> >>
> >> https://vtk.org/wp-content/uploads/2015/04/file-formats.pdf   (PDF was
> >> created in 2003)
> >>
> >>
> https://kitware.github.io/vtk-examples/site/VTKFileFormats/#xml-file-formats
> >>
> >
> > Yes I know. I’ve tied myself in knots for years because the of the
> > assumption that the offset had to be an int.
> >
> > Credit for the discovery goes to Carsten Uphoff. He showed me this.
> >
> > Cheers,
> > Dave
> >
> >
> >
> >> > then the size of the offset is now UInt64 and now large files can be
> >> > finally written.
> >> >
> >> >
> >> > Cheers,
> >> > Dave
> >> >
> >> >
> >> >
> >> >
> >> >>
> >> >> Mike Michell  writes:
> >> >>
> >> >> > Thanks for the note.
> >> >> > I understood that PETSc calculates the offsets for me through
> >> "boffset"
> >> >> > variable in plexvtu.c file. Please correct me if it is wrong.
> >> >> >
> >> >> > If plexvtu.c has a bug, it could be around "write file header"
> part in
> >> >> > which the boffset is also computed. Is this correct? I am not using
> >> >> complex
> >> >> > number.
> >> >> > There are several mixed parts among "Int32, UInt8, PetscInt_FMT,
> >> >> > PetscInt64_FMT" in writing the header.
> >> >> >
> >> >> > Which combination of those flags is correct for 64bit indices? I am
> >> gonna
> >> >> > modify plexvtu.c file with "#if defined(PETSC_USE_64BIT_INDICES)"
> >> >> > statement, but I do not know what is the correct form of the header
> >> flag
> >> >> > for 64bit indices.
> >> >> >
> >> >> > It is also confusing to me:
> >> >> > boffset += gpiece[r].ncells * sizeof(PetscInt) + sizeof(int);
> >> >> > How is sizeof(PetscInt) different from sizeof(int)?
> >> >> >
> >> >> > Thanks,
> >> >> > Mike
> >> >> >
> >> >> >
> >> >> >> On Tue, Feb 14, 2023 at 11:45 AM Mike Michell <
> mi.mike1...@gmail.com
> >> >
> >> >> >> wrote:
> >> >> >>
> >> >> >>> I was trying to modify the header flags from "Int32" to "Int64",
> but
> >> >> the
> >> >> >>> problem was not resolved. Could I get any additional comments?
> >> >> >>>
> >> >> >>
> >> >> >> The calculated offsets are not correct I think.
> >> >> >>
> >> >> >>   Matt
> >> >> >>
> >> >> >>
> >> >> >>> Thanks,
> >> >> >>> Mike
> >> >> >>>
> >> >> >>>
> >> >>  Thanks for the comments.
> >> >>  To be precise on the question, the entire part of the header of
> the
> >> >> .vtu
> >> >>  file is attached:
> >> >> 
> >> >>  
> >> >>   >> >> byte_order="LittleEndian">
> >> >>    
> >> >>  
> >> >>    
> >> >>   >> >> NumberOfComponents="3"
> >> >>  format="appended" offset="0" />
> >> >>    
> >> >>    
> >> >>   >> >>  NumberOfComponents="1" format="appended" offset="116932" />
> >> >>   >> >>   NumberOfComponents="1" format="appended" offset="372936" />
> >> >>   >> >>   NumberOfComponents="1" format="appended" offset="404940" />
> >> >>    
> >> >>    
> >> >>   NumberOfComponents="1"
> >> >>  format="appended" offset="408944" />
> >> >>    
> >> >>    
> >> >>   Name="Vec_0x37c89c0_4Field_0.0"
> >> >>  

Re: [petsc-users] PetscViewer with 64bit

[Jed Brown]

Thanks, Dave.

Mike, can you test that this branch works with your large problems? I tested 
that .vtu works in parallel for small problems, where works = loads correctly 
in Paraview and VisIt.

https://gitlab.com/petsc/petsc/-/merge_requests/6081

Dave May  writes:

> On Tue 14. Feb 2023 at 21:27, Jed Brown  wrote:
>
>> Dave May  writes:
>>
>> > On Tue 14. Feb 2023 at 17:17, Jed Brown  wrote:
>> >
>> >> Can you share a reproducer? I think I recall the format requiring
>> certain
>> >> things to be Int32.
>> >
>> >
>> > By default, the byte offset used with the appended data format is
>> UInt32. I
>> > believe that’s where the sizeof(int) is coming from. This default is
>> > annoying as it limits the total size of your appended data to be < 3 GB.
>> > That said, in the opening of the paraview file you can add this attribute
>> >
>> > header_type="UInt64"
>>
>> You mean in the header of the .vtu?
>
>
> Yeah, within the open VTKFile tag.
> Like this
> < VTKFile type=“xxx”,  byte_order="LittleEndian" header_type="UInt64" >
>
> Do you happen to have an example or pointers to docs describing this
>> feature?
>
>
> Example yes - will send it to you tomorrow. Docs… not really. Only stuff
> like this
>
> https://kitware.github.io/paraview-docs/latest/python/paraview.simple.XMLPStructuredGridWriter.html
>
>
> https://kitware.github.io/paraview-docs/v5.8.0/python/paraview.simple.XMLMultiBlockDataWriter.html
>
> All the writers seem to support it.
>
>
> Can we always do this?
>
>
> Yep!
>
>
> It isn't mentioned in these:
>>
>> https://vtk.org/wp-content/uploads/2015/04/file-formats.pdf   (PDF was
>> created in 2003)
>>
>> https://kitware.github.io/vtk-examples/site/VTKFileFormats/#xml-file-formats
>>
>
> Yes I know. I’ve tied myself in knots for years because the of the
> assumption that the offset had to be an int.
>
> Credit for the discovery goes to Carsten Uphoff. He showed me this.
>
> Cheers,
> Dave
>
>
>
>> > then the size of the offset is now UInt64 and now large files can be
>> > finally written.
>> >
>> >
>> > Cheers,
>> > Dave
>> >
>> >
>> >
>> >
>> >>
>> >> Mike Michell  writes:
>> >>
>> >> > Thanks for the note.
>> >> > I understood that PETSc calculates the offsets for me through
>> "boffset"
>> >> > variable in plexvtu.c file. Please correct me if it is wrong.
>> >> >
>> >> > If plexvtu.c has a bug, it could be around "write file header" part in
>> >> > which the boffset is also computed. Is this correct? I am not using
>> >> complex
>> >> > number.
>> >> > There are several mixed parts among "Int32, UInt8, PetscInt_FMT,
>> >> > PetscInt64_FMT" in writing the header.
>> >> >
>> >> > Which combination of those flags is correct for 64bit indices? I am
>> gonna
>> >> > modify plexvtu.c file with "#if defined(PETSC_USE_64BIT_INDICES)"
>> >> > statement, but I do not know what is the correct form of the header
>> flag
>> >> > for 64bit indices.
>> >> >
>> >> > It is also confusing to me:
>> >> > boffset += gpiece[r].ncells * sizeof(PetscInt) + sizeof(int);
>> >> > How is sizeof(PetscInt) different from sizeof(int)?
>> >> >
>> >> > Thanks,
>> >> > Mike
>> >> >
>> >> >
>> >> >> On Tue, Feb 14, 2023 at 11:45 AM Mike Michell > >
>> >> >> wrote:
>> >> >>
>> >> >>> I was trying to modify the header flags from "Int32" to "Int64", but
>> >> the
>> >> >>> problem was not resolved. Could I get any additional comments?
>> >> >>>
>> >> >>
>> >> >> The calculated offsets are not correct I think.
>> >> >>
>> >> >>   Matt
>> >> >>
>> >> >>
>> >> >>> Thanks,
>> >> >>> Mike
>> >> >>>
>> >> >>>
>> >>  Thanks for the comments.
>> >>  To be precise on the question, the entire part of the header of the
>> >> .vtu
>> >>  file is attached:
>> >> 
>> >>  
>> >>  > >> byte_order="LittleEndian">
>> >>    
>> >>  
>> >>    
>> >>  > >> NumberOfComponents="3"
>> >>  format="appended" offset="0" />
>> >>    
>> >>    
>> >>  > >>  NumberOfComponents="1" format="appended" offset="116932" />
>> >>  > >>   NumberOfComponents="1" format="appended" offset="372936" />
>> >>  > >>   NumberOfComponents="1" format="appended" offset="404940" />
>> >>    
>> >>    
>> >>  > >>  format="appended" offset="408944" />
>> >>    
>> >>    
>> >>  > >>  NumberOfComponents="1" format="appended" offset="424948" />
>> >>    
>> >>  
>> >>  
>> >>    
>> >>  > >> NumberOfComponents="3"
>> >>  format="appended" offset="463928" />
>> >>    
>> >>    
>> >>  > >>  NumberOfComponents="1" format="appended" offset="580860" />
>> >>  > >>   NumberOfComponents="1" format="appended" offset="836864" />
>> >>  > >>   NumberOfComponents="1" format="appended" offset="868868" />
>> >>    
>> >>    
>> >>  > >>  format="appended" 

Re: [petsc-users] dmplex overlap questions

[Matthew Knepley]

On Thu, Feb 16, 2023 at 1:09 PM Lawrence Mitchell  wrote:

> On Thu, 16 Feb 2023 at 16:43, Matthew Knepley  wrote:
> >
> > On Thu, Feb 16, 2023 at 10:54 AM Lawrence Mitchell  wrote:
> >>
> >> Hi Blaise,
> >>
> >> On Thu, 16 Feb 2023 at 15:17, Blaise Bourdin 
> wrote:
> >> >
> >> > Hi,
> >> >
> >> > I am trying to implement a non-local finite elements reconstruction
> operator in parallel.
> >> >
> >> > Given a dmplex distributed with an overlap, is there a way to figure
> out which cells are in the overlap and which are not?
> >>
> >> Yes. Get the pointSF of the DM, and the cell chart
> >>
> >> DMPlexGetPointSF(dm, );
> >> DMPlexGetHeightStratum(dm, 0, , );
> >>
> >> Now get the graph (specifically ilocal of the sf):
> >>
> >> PetscSFGetGraph(sf, NULL, , ,  NULL);
> >>
> >> Now any value of ilocal that lies in [cstart, cend) is a cell which is
> >> not owned by this process (i.e. in the overlap). Note that ilocal can
> >> be NULL which just means it is the identity map [0, ..., nleaves), so
> >> you just intersect [cstart, cend) with [0, nleaves) in that case to
> >> find the overlap cells.
> >>
> >> But that is very unlikely to be true, so:
> >
> >
> > Note that you can use
> >
> >   PetscFindInt(nleaves, ilocal, cell, );
>
> Modulo argument order, is it not PetscFindInt(cell, nleaves, ilocal, )?
>

You are right.


> I guess one should probably document that PetscSFSetGraph ensures that
> ilocal is sorted.
>

Yes, I thought it was already there, but does not seem so.

  Thanks,

Matt


> Lawrence

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

Re: [petsc-users] dmplex overlap questions

[Lawrence Mitchell]

On Thu, 16 Feb 2023 at 16:43, Matthew Knepley  wrote:
>
> On Thu, Feb 16, 2023 at 10:54 AM Lawrence Mitchell  wrote:
>>
>> Hi Blaise,
>>
>> On Thu, 16 Feb 2023 at 15:17, Blaise Bourdin  wrote:
>> >
>> > Hi,
>> >
>> > I am trying to implement a non-local finite elements reconstruction 
>> > operator in parallel.
>> >
>> > Given a dmplex distributed with an overlap, is there a way to figure out 
>> > which cells are in the overlap and which are not?
>>
>> Yes. Get the pointSF of the DM, and the cell chart
>>
>> DMPlexGetPointSF(dm, );
>> DMPlexGetHeightStratum(dm, 0, , );
>>
>> Now get the graph (specifically ilocal of the sf):
>>
>> PetscSFGetGraph(sf, NULL, , ,  NULL);
>>
>> Now any value of ilocal that lies in [cstart, cend) is a cell which is
>> not owned by this process (i.e. in the overlap). Note that ilocal can
>> be NULL which just means it is the identity map [0, ..., nleaves), so
>> you just intersect [cstart, cend) with [0, nleaves) in that case to
>> find the overlap cells.
>>
>> But that is very unlikely to be true, so:
>
>
> Note that you can use
>
>   PetscFindInt(nleaves, ilocal, cell, );

Modulo argument order, is it not PetscFindInt(cell, nleaves, ilocal, )?

I guess one should probably document that PetscSFSetGraph ensures that
ilocal is sorted.

Lawrence

Re: [petsc-users] Question for Petsc

[Stefano Zampini]

For bddc, you can also take a look at
https://gitlab.com/petsc/petsc/-/blob/main/src/ksp/ksp/tutorials/ex71.c

On Thu, Feb 16, 2023, 19:41 Matthew Knepley  wrote:

> On Thu, Feb 16, 2023 at 9:14 AM ziming xiong 
> wrote:
>
>> Hello，
>> I want to use Petsc to implement high performance computing, and I mainly
>> want to apply DDM methods to parallel computing. I have implemented some of
>> the DDM methods (such as ASM, Bjacobi, etc.), but I don't understand the
>> PCBDDC method. The official example (src/ksp/ksp/tutorials/ex59.c.html) is
>> too complicated, so I have not been able to figure out the setting process.
>> I would like to ask you if you have other simple and clearer examples for
>> reference.
>>
>
> You could look at the paper:
> https://epubs.siam.org/doi/abs/10.1137/15M1025785
>
>
>> Secondly, I would like to apply mklPardiso to Petsc. But not work, can u
>> help me figure out the problem? i use oneAPI for the mklpardiso, and when i
>> configure, i give the blaslapack lib.
>>
>
> You should reconfigure with --download-mkl_pardiso
>
>   Thanks,
>
>  Matt
>
>
>> there are the errors:
>>
>> [0]PETSC ERROR: See
>> https://petsc.org/release/overview/linear_solve_table/ for possible LU
>> and Cholesky solvers
>> [0]PETSC ERROR: Could not locate solver type mkl_pardiso for
>> factorization type LU and matrix type seqaij. Perhaps you must ./configure
>> with --download-mkl_pardiso
>> [0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting.
>> [0]PETSC ERROR: Petsc Release Version 3.18.3, Dec 28, 2022
>> [0]PETSC ERROR:
>> C:\Users\XiongZiming\Desktop\test_petsc_FEM\test_petsc_fem\x64\Release\test_petsc_fem.exe
>> on a arch-mswin-c-debug named lmeep-329 by XiongZiming Thu Feb 16 15:05:14
>> 2023
>> [0]PETSC ERROR: Configure options --with-cc="win32fe cl" --with-fc=0
>> --with-cxx="win32fe cl" --with-shared-libraries=0
>> --with-mpi-include="[/cygdrive/c/PROGRA~2/Intel/MPI/Include,/cygdrive/c/PROGRA~2/Intel/MPI/Include/x64]"
>> --with-mpi-lib="-L/cygdrive/c/PROGRA~2/Intel/MPI/lib/x64 msmpifec.lib
>> msmpi.lib" --with-mpiexec=/cygdrive/c/PROGRA~1/Microsoft_MPI/Bin/mpiexec
>> --with-blaslapack-lib="-L/cygdrive/c/PROGRA~2/Intel/oneAPI/mkl/2023.0.0/lib/intel64
>> mkl_intel_lp64_dll.lib mkl_sequential_dll.lib mkl_core_dll.lib"
>>
>> Thanks,
>> Ziming XIONG
>>
>
>
> --
> 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/
> 
>

Re: [petsc-users] dmplex overlap questions

[Matthew Knepley]

On Thu, Feb 16, 2023 at 10:54 AM Lawrence Mitchell  wrote:

> Hi Blaise,
>
> On Thu, 16 Feb 2023 at 15:17, Blaise Bourdin  wrote:
> >
> > Hi,
> >
> > I am trying to implement a non-local finite elements reconstruction
> operator in parallel.
> >
> > Given a dmplex distributed with an overlap, is there a way to figure out
> which cells are in the overlap and which are not?
>
> Yes. Get the pointSF of the DM, and the cell chart
>
> DMPlexGetPointSF(dm, );
> DMPlexGetHeightStratum(dm, 0, , );
>
> Now get the graph (specifically ilocal of the sf):
>
> PetscSFGetGraph(sf, NULL, , ,  NULL);
>
> Now any value of ilocal that lies in [cstart, cend) is a cell which is
> not owned by this process (i.e. in the overlap). Note that ilocal can
> be NULL which just means it is the identity map [0, ..., nleaves), so
> you just intersect [cstart, cend) with [0, nleaves) in that case to
> find the overlap cells.
>
> But that is very unlikely to be true, so:
>

Note that you can use

  PetscFindInt(nleaves, ilocal, cell, );

as well. I do this a lot in the library.

  Thanks,

 Matt


> for (PetscInt i = 0; i < nleaves; i++) {
> if (cstart <= ilocal[i] && ilocal[i] < cend) {
>// i is an overlap cell
> }
> }
> > Alternatively, suppose that I distribute the same DM with and without an
> overlap. Is there any warranty that the distributions are compatible (i.e.
> coincide when the overlap is ignored)? If this is the case, can I assume
> that the non-overlap cells are numbered first in the overlapped dm?
>
> If you do:
>
> DMPlexDistribute(dm, 0, , );
> DMPlexDistributeOverlap(paralleldm, 1, , );
>
> Then paralleldm and overlapdm will be compatible, and I think it is
> still the case that the overlap cells are numbered last (and
> contiguously).
>
> If you just do DMPlexDistribute(dm, 1, , ) then
> you obviously don't have the non-overlapped one to compare, but it is
> in this case still true that the overlap cells are numbered last.
>
> Thanks,
>
> Lawrence
>


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

Re: [petsc-users] Question for Petsc

[Matthew Knepley]

On Thu, Feb 16, 2023 at 9:14 AM ziming xiong 
wrote:

> Hello，
> I want to use Petsc to implement high performance computing, and I mainly
> want to apply DDM methods to parallel computing. I have implemented some of
> the DDM methods (such as ASM, Bjacobi, etc.), but I don't understand the
> PCBDDC method. The official example (src/ksp/ksp/tutorials/ex59.c.html) is
> too complicated, so I have not been able to figure out the setting process.
> I would like to ask you if you have other simple and clearer examples for
> reference.
>

You could look at the paper:
https://epubs.siam.org/doi/abs/10.1137/15M1025785


> Secondly, I would like to apply mklPardiso to Petsc. But not work, can u
> help me figure out the problem? i use oneAPI for the mklpardiso, and when i
> configure, i give the blaslapack lib.
>

You should reconfigure with --download-mkl_pardiso

  Thanks,

 Matt


> there are the errors:
>
> [0]PETSC ERROR: See https://petsc.org/release/overview/linear_solve_table/
> for possible LU and Cholesky solvers
> [0]PETSC ERROR: Could not locate solver type mkl_pardiso for factorization
> type LU and matrix type seqaij. Perhaps you must ./configure with
> --download-mkl_pardiso
> [0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting.
> [0]PETSC ERROR: Petsc Release Version 3.18.3, Dec 28, 2022
> [0]PETSC ERROR:
> C:\Users\XiongZiming\Desktop\test_petsc_FEM\test_petsc_fem\x64\Release\test_petsc_fem.exe
> on a arch-mswin-c-debug named lmeep-329 by XiongZiming Thu Feb 16 15:05:14
> 2023
> [0]PETSC ERROR: Configure options --with-cc="win32fe cl" --with-fc=0
> --with-cxx="win32fe cl" --with-shared-libraries=0
> --with-mpi-include="[/cygdrive/c/PROGRA~2/Intel/MPI/Include,/cygdrive/c/PROGRA~2/Intel/MPI/Include/x64]"
> --with-mpi-lib="-L/cygdrive/c/PROGRA~2/Intel/MPI/lib/x64 msmpifec.lib
> msmpi.lib" --with-mpiexec=/cygdrive/c/PROGRA~1/Microsoft_MPI/Bin/mpiexec
> --with-blaslapack-lib="-L/cygdrive/c/PROGRA~2/Intel/oneAPI/mkl/2023.0.0/lib/intel64
> mkl_intel_lp64_dll.lib mkl_sequential_dll.lib mkl_core_dll.lib"
>
> Thanks,
> Ziming XIONG
>


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

Re: [petsc-users] dmplex overlap questions

[Blaise Bourdin]

On Feb 16, 2023, at 10:54 AM, Lawrence Mitchell  wrote:


Hi Blaise,

On Thu, 16 Feb 2023 at 15:17, Blaise Bourdin  wrote:

Hi,

I am trying to implement a non-local finite elements reconstruction operator in parallel.

Given a dmplex distributed with an overlap, is there a way to figure out which cells are in the overlap and which are not?


Yes. Get the pointSF of the DM, and the cell chart

DMPlexGetPointSF(dm, );
DMPlexGetHeightStratum(dm, 0, , );

Now get the graph (specifically ilocal of the sf):

PetscSFGetGraph(sf, NULL, , ,  NULL);

Now any value of ilocal that lies in [cstart, cend) is a cell which is
not owned by this process (i.e. in the overlap). Note that ilocal can
be NULL which just means it is the identity map [0, ..., nleaves), so
you just intersect [cstart, cend) with [0, nleaves) in that case to
find the overlap cells.

But that is very unlikely to be true, so:

for (PetscInt i = 0; i < nleaves; i++) {
   if (cstart <= ilocal[i] && ilocal[i] < cend) {
  // i is an overlap cell
   }
}
Alternatively, suppose that I distribute the same DM with and without an overlap. Is there any warranty that the distributions are compatible (i.e. coincide when the overlap is ignored)? If this is the case, can I assume that the non-overlap
 cells are numbered first in the overlapped dm?


If you do:

DMPlexDistribute(dm, 0, , );
DMPlexDistributeOverlap(paralleldm, 1, , );

Then paralleldm and overlapdm will be compatible, and I think it is
still the case that the overlap cells are numbered last (and
contiguously).

If you just do DMPlexDistribute(dm, 1, , ) then
you obviously don't have the non-overlapped one to compare, but it is
in this case still true that the overlap cells are numbered last.







Excellent. That is what I was hoping.


Regards,
Blaise






Thanks,

Lawrence





















— 
Canada Research Chair in Mathematical and Computational Aspects of Solid Mechanics (Tier 1)
Professor, Department of Mathematics & Statistics
Hamilton Hall room 409A, McMaster University
1280 Main Street West, Hamilton, Ontario L8S 4K1, Canada 
https://www.math.mcmaster.ca/bourdin | +1 (905) 525 9140 ext. 27243

Re: [petsc-users] dmplex overlap questions

[Lawrence Mitchell]

Hi Blaise,

On Thu, 16 Feb 2023 at 15:17, Blaise Bourdin  wrote:
>
> Hi,
>
> I am trying to implement a non-local finite elements reconstruction operator 
> in parallel.
>
> Given a dmplex distributed with an overlap, is there a way to figure out 
> which cells are in the overlap and which are not?

Yes. Get the pointSF of the DM, and the cell chart

DMPlexGetPointSF(dm, );
DMPlexGetHeightStratum(dm, 0, , );

Now get the graph (specifically ilocal of the sf):

PetscSFGetGraph(sf, NULL, , ,  NULL);

Now any value of ilocal that lies in [cstart, cend) is a cell which is
not owned by this process (i.e. in the overlap). Note that ilocal can
be NULL which just means it is the identity map [0, ..., nleaves), so
you just intersect [cstart, cend) with [0, nleaves) in that case to
find the overlap cells.

But that is very unlikely to be true, so:

for (PetscInt i = 0; i < nleaves; i++) {
if (cstart <= ilocal[i] && ilocal[i] < cend) {
   // i is an overlap cell
}
}
> Alternatively, suppose that I distribute the same DM with and without an 
> overlap. Is there any warranty that the distributions are compatible (i.e. 
> coincide when the overlap is ignored)? If this is the case, can I assume that 
> the non-overlap cells are numbered first in the overlapped dm?

If you do:

DMPlexDistribute(dm, 0, , );
DMPlexDistributeOverlap(paralleldm, 1, , );

Then paralleldm and overlapdm will be compatible, and I think it is
still the case that the overlap cells are numbered last (and
contiguously).

If you just do DMPlexDistribute(dm, 1, , ) then
you obviously don't have the non-overlapped one to compare, but it is
in this case still true that the overlap cells are numbered last.

Thanks,

Lawrence

Re: [petsc-users] Question about preconditioner

[Barry Smith]

   If your  matrix has the form 

(   A   B )
(   C   0 )

then often PCFIELDSPLIT can be a useful preconditioner with the option 
-pc_fieldsplit_detect_saddle_point



> On Feb 16, 2023, at 2:42 AM, user_gong Kim  wrote:
> 
>  
> Hello,
> 
>  
> There are some questions about some preconditioners.
> 
> The questions are from problem Au=b. The global matrix A has zero value 
> diagonal terms.
> 
> 1. Which preconditioner is preferred for matrix A which has zero value in 
> diagonal terms?
> The most frequently used basic 2 preconditioners are jacobi and SOR (gauss 
> seidel). As people knows both methods should have non zero diagonal terms. 
> Although the improved method is applied in PETSc, jacobi can also solve the 
> case with zero diagonal term, but I ask because I know that it is not 
> recommended.
> 
> 2. Second question is about running code with the two command options 
> below in a single process.
> 1st command : -ksp_type gmres -pc_type bjacobi -sub_pc_type jacobi
> 2nd command : -ksp_type gmres -pc_type hpddm -sub_pc_type jacobi
> When domain decomposition methods such as bjacobi or hpddm are parallel, the 
> global matrix is divided for each process. As far as I know, running it in a 
> single process should eventually produce the same result if the sub pc type 
> is the same. However, in the second option, ksp did not converge.
> In this case, I wonder how to analyze the situation.
> How can I monitor and see the difference between the two?
> 
>  
>  
> Thanks,
> 
> Hyung Kim
> 

[petsc-users] dmplex overlap questions

[Blaise Bourdin]

Hi,

I am trying to implement a non-local finite elements reconstruction operator in 
parallel.

Given a dmplex distributed with an overlap, is there a way to figure out which 
cells are in the overlap and which are not?
Alternatively, suppose that I distribute the same DM with and without an 
overlap. Is there any warranty that the distributions are compatible (i.e. 
coincide when the overlap is ignored)? If this is the case, can I assume that 
the non-overlap cells are numbered first in the overlapped dm?


Regards,
Blaise


— 
Canada Research Chair in Mathematical and Computational Aspects of Solid 
Mechanics (Tier 1)
Professor, Department of Mathematics & Statistics
Hamilton Hall room 409A, McMaster University
1280 Main Street West, Hamilton, Ontario L8S 4K1, Canada 
https://www.math.mcmaster.ca/bourdin | +1 (905) 525 9140 ext. 27243

[petsc-users] Question for Petsc

[ziming xiong]

Hello，
I want to use Petsc to implement high performance computing, and I mainly
want to apply DDM methods to parallel computing. I have implemented some of
the DDM methods (such as ASM, Bjacobi, etc.), but I don't understand the
PCBDDC method. The official example (src/ksp/ksp/tutorials/ex59.c.html) is
too complicated, so I have not been able to figure out the setting process.
I would like to ask you if you have other simple and clearer examples for
reference.

Secondly, I would like to apply mklPardiso to Petsc. But not work, can u
help me figure out the problem? i use oneAPI for the mklpardiso, and when i
configure, i give the blaslapack lib.

there are the errors:

[0]PETSC ERROR: See https://petsc.org/release/overview/linear_solve_table/
for possible LU and Cholesky solvers
[0]PETSC ERROR: Could not locate solver type mkl_pardiso for factorization
type LU and matrix type seqaij. Perhaps you must ./configure with
--download-mkl_pardiso
[0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.18.3, Dec 28, 2022
[0]PETSC ERROR:
C:\Users\XiongZiming\Desktop\test_petsc_FEM\test_petsc_fem\x64\Release\test_petsc_fem.exe
on a arch-mswin-c-debug named lmeep-329 by XiongZiming Thu Feb 16 15:05:14
2023
[0]PETSC ERROR: Configure options --with-cc="win32fe cl" --with-fc=0
--with-cxx="win32fe cl" --with-shared-libraries=0
--with-mpi-include="[/cygdrive/c/PROGRA~2/Intel/MPI/Include,/cygdrive/c/PROGRA~2/Intel/MPI/Include/x64]"
--with-mpi-lib="-L/cygdrive/c/PROGRA~2/Intel/MPI/lib/x64 msmpifec.lib
msmpi.lib" --with-mpiexec=/cygdrive/c/PROGRA~1/Microsoft_MPI/Bin/mpiexec
--with-blaslapack-lib="-L/cygdrive/c/PROGRA~2/Intel/oneAPI/mkl/2023.0.0/lib/intel64
mkl_intel_lp64_dll.lib mkl_sequential_dll.lib mkl_core_dll.lib"

Thanks,
Ziming XIONG

Re: [petsc-users] Question about preconditioner

[Matthew Knepley]

On Thu, Feb 16, 2023 at 2:43 AM user_gong Kim  wrote:

>
>
> Hello,
>
>
>
> There are some questions about some preconditioners.
>
> The questions are from problem Au=b. The global matrix A has zero value
> diagonal terms.
>
> 1. Which preconditioner is preferred for matrix A which has zero
> value in diagonal terms?
> The most frequently used basic 2 preconditioners are jacobi and SOR (gauss
> seidel). As people knows both methods should have non zero diagonal terms.
> Although the improved method is applied in PETSc, jacobi can also solve the
> case with zero diagonal term, but I ask because I know that it is not
> recommended.
>
> 2. Second question is about running code with the two command options
> below in a single process.
> 1st command : -ksp_type gmres -pc_type bjacobi -sub_pc_type jacobi
> 2nd command : -ksp_type gmres -pc_type hpddm -sub_pc_type jacobi
> When domain decomposition methods such as bjacobi or hpddm are parallel,
> the global matrix is divided for each process. As far as I know, running it
> in a single process should eventually produce the same result if the sub pc
> type is the same. However, in the second option, ksp did not converge.
> In this case, I wonder how to analyze the situation.
> How can I monitor and see the difference between the two?
>
> Pierre is correct. I will just note that the best way to use PETSc is
generally to find the preconditioner you need
by looking in literature for what has been successful for other people, and
then reproducing it in PETSc, which
should be easy.

  Thanks,

 Matt


>
>
> Thanks,
>
> Hyung Kim
>


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

Re: [petsc-users] Question about preconditioner

[Pierre Jolivet]

> On 16 Feb 2023, at 8:43 AM, user_gong Kim  wrote:
> 
> 
>  
> 
> Hello,
> 
>  
> 
> There are some questions about some preconditioners.
> 
> The questions are from problem Au=b. The global matrix A has zero value 
> diagonal terms.
> 
> 1. Which preconditioner is preferred for matrix A which has zero value in 
> diagonal terms?
> 

This question has not a single answer. It all depends on where your A and b are 
coming from.

> The most frequently used basic 2 preconditioners are jacobi and SOR (gauss 
> seidel).
> 
They are not the most frequently used. And rightfully so, as they very often 
can’t handle non-trivial systems.

> As people knows both methods should have non zero diagonal terms. Although 
> the improved method is applied in PETSc, jacobi can also solve the case with 
> zero diagonal term, but I ask because I know that it is not recommended.
> 
> 2. Second question is about running code with the two command options 
> below in a single process.
> 1st command : -ksp_type gmres -pc_type bjacobi -sub_pc_type jacobi
> 2nd command : -ksp_type gmres -pc_type hpddm -sub_pc_type jacobi
> When domain decomposition methods such as bjacobi or hpddm are parallel, the 
> global matrix is divided for each process. As far as I know, running it in a 
> single process should eventually produce the same result if the sub pc type 
> is the same. However, in the second option, ksp did not converge.
> 
1st command: it’s pointless to couple PCBJACOBI with PCJABOCI, it’s equivalent 
to only using PCJABOBI.
2nd command: it’s pointless to use PCHPDDM if you don’t specify in some way how 
to coarsen your problem (either algebraically or via an auxiliary operator). 
You just have a single level (equivalent to PCBJACOBI), but its options are 
prefixed by -pc_hpddm_coarse_ instead of -sub_
Again, both sets of options do not make sense.
If you want, you could share your A and b (or tell us what you are 
discretizing) and we will be able to provide a better feedback.

Thanks,
Pierre

> In this case, I wonder how to analyze the situation.
> How can I monitor and see the difference between the two?
> 
>  
> 
>  
> 
> Thanks,
> 
> Hyung Kim