[petsc-users] Error with KSPSetUp and MatNest

Thu Apr 11 05:07:23 CDT 2019

On Thu, Apr 11, 2019 at 5:57 AM Manuel Colera Rico <m.colera at upm.es> wrote:

> MatCreateLocalRef() seems to be useful to extract submatrices from a
> global matrix. But how to construct A from the beginning from its
> submatrices? I.e., with MatNest, I can just do:
>
>     Mat A;
>
>     Mat matA[] = {A00, A01, NULL, A11};
>
>     MatCreateNest(PETSC_COMM_SELF, 2, NULL, 2, NULL, matA, &A);
>
> without worrying about memory preallocations for A. Can I do the same with
> MatCreateLocalRef(), or do I first have to preallocate memory for A?
>
> It is a good point. You would have to preallocate the large matrix.
However, the preallocation is just the sum for the
preallocations of the small matrices, which you are already specifying. It
would be no problem to give a similar
LocalRef interface to preallocation to do that sum automatically for the
user.

Moreover, we now have a MatPreallocator interface. You can run your matrix
assembly twice. Once with the
MatPreallocator matrix, which can then preallocate you matrix exactly, and
then again with the real matrix. This
is more expensive of course, but in many applications, this expense is
negligible compared to the solver.

  Thanks,

     Matt

> Manuel
>
> ---
>
> On 4/11/19 11:30 AM, Matthew Knepley wrote:
>
> On Thu, Apr 11, 2019 at 2:07 AM Smith, Barry F. <bsmith at mcs.anl.gov>
> wrote:
>
>>
>>    Matt,
>>
>>      You can't have some functionality in the public API for a library
>> and then criticize someone for using it. I looked at the manual page for
>> MatCreateNest() and it doesn't say in big letters "don't use this", nor
>> does the compiler tell the user "don't use this". Either MatCreateNest()
>> and/or MATNEST is so horrible it needs to be stripped from the face of the
>> earth or it needs to be supported. It can't be in this nebulous state of
>> "yeah it exists but don't ever use it".
>>
>
> I agree this is a documentation problem. That is why I try to provide
> complete explanations for my recommendations.
>
>
>>    Barry
>>
>>    In order to support easy testing of MatCreateNest() based user code
>> perhaps we could add a little bit of code that allowed users to call either
>> MatCreateNest() or MatNestSetSubMats() but with an option to have it
>> generate a regular AIJ matrix instead. For example -mat_nest_use_aij (and
>> MatNestUseAIJ()) This would allow  users to always write and think about
>> their code as nested matrices but have it "fall back" to AIJ for
>> testing/debugging purposes.
>>
>
> I think we should remove MatCreateNest() completely. We should have 1
> mechanism for getting submatrices for creation, not 2, so we retain
> MatCreateLocalRef(). Then if -mat_type nest was specified, it gives you
> back an actually submatrix, otherwise it gives a view. This would do
> everything that we currently do without this horrible MatNest interface
> bubbling to the top.
>
>    Matt
>
>
>> > On Apr 10, 2019, at 11:49 AM, Manuel Colera Rico via petsc-users <
>> petsc-users at mcs.anl.gov> wrote:
>> >
>> > Thank you for your answer, Matt. In the MWE example attached before,
>> both Nest vectors (the r.h.s. of the system and the vector of unknowns) are
>> composed of the same number of blocks (2). Indeed, PETSc is able to solve
>> the system if KSPSetUp() is not called, so the system/MatNest/MatVec's must
>> not incompatible at all. Therefore, I wonder if I have missed to called
>> something before this routine or if this is a KSPSetUp's bug.
>> >
>> > Of course one can always directly define a single matrix and a single
>> vector, but I find it easier to work with Nest matrices and vectors.
>> Moreover, I think that the moment to use them is from the beginning... once
>> all the code is developed, it is very hard to switch matrices types.
>> >
>> > Regards,
>> >
>> > Manuel
>> >
>> > ---
>> >
>> >
>> >
>> > On 4/10/19 5:41 PM, Matthew Knepley wrote:
>> >> On Wed, Apr 10, 2019 at 11:29 AM Manuel Colera Rico via petsc-users <
>> petsc-users at mcs.anl.gov> wrote:
>> >> Hello,
>> >>
>> >> I am trying to solve a system whose matrix is of type MatNest. If I
>> >> don't use KSPSetUp(), everything is fine. However, if I use that
>> >> routine, I get the following error:
>> >>
>> >> 0]PETSC ERROR: --------------------- Error Message
>> >> --------------------------------------------------------------
>> >> [0]PETSC ERROR: Invalid argument
>> >> [0]PETSC ERROR: Nest vector arguments 1 and 2 have different numbers
>> of
>> >> blocks.
>> >>
>> >> This seems self-explanatory. Nest vectors must have the same number of
>> blocks to be compatible.
>> >>
>> >> More broadly, there should be no reason to use Nest vectors or
>> matrices. It is an optimization to
>> >> be used at the very end, only after you have profiled the code and
>> seen that its important. You can
>> >> do everything you want to do without ever touching Nest, and it looks
>> like the Nest interface is a
>> >> problem for your code right now.
>> >>
>> >>   Thanks,
>> >>
>> >>     Matt
>> >>
>> >> [0]PETSC ERROR: See
>> http://www.mcs.anl.gov/petsc/documentation/faq.html
>> >> for trouble shooting.
>> >> [0]PETSC ERROR: Petsc Release Version 3.11.0, unknown
>> >> [0]PETSC ERROR:
>> >>
>> /home/manu/Documents/FEM-fluids/C-codes/CLG2-ConvectionDiffusion/Debug/CLG2-ConvectionDiffusion
>>
>> >> on a mcr_20190405 named mancolric by Unknown Wed Apr 10 17:20:16 2019
>> >> [0]PETSC ERROR: Configure options --with-cc=gcc --with-cxx=g++
>> >> --with-fc=gfortran COPTFLAGS="-O3 -march=native -mtune=native"
>> >> CXXOPTFLAGS="-O3 -march=native -mtune=native" FOPTFLAGS="-O3
>> >> -march=native -mtune=native" --with-debugging=0 --download-fblaslapack
>> >> --download--f2cblaslapack --download-mpich --download--hypre
>> >> --download-scalapack --download-mumps --download-suitesparse
>> >> --download-ptscotch --download-pastix --with-matlab --with-openmp
>> >> [0]PETSC ERROR: #1 VecCopy_Nest() line 68 in
>> >> /opt/PETSc_library/petsc-3.11.0/src/vec/vec/impls/nest/vecnest.c
>> >> [0]PETSC ERROR: #2 VecCopy() line 1614 in
>> >> /opt/PETSc_library/petsc-3.11.0/src/vec/vec/interface/vector.c
>> >> [0]PETSC ERROR: #3 KSPInitialResidual() line 63 in
>> >> /opt/PETSc_library/petsc-3.11.0/src/ksp/ksp/interface/itres.c
>> >> [0]PETSC ERROR: #4 KSPSolve_GMRES() line 236 in
>> >> /opt/PETSc_library/petsc-3.11.0/src/ksp/ksp/impls/gmres/gmres.c
>> >> [0]PETSC ERROR: #5 KSPSolve() line 782 in
>> >> /opt/PETSc_library/petsc-3.11.0/src/ksp/ksp/interface/itfunc.c
>> >> [0]PETSC ERROR: #6 mwe() line 55 in ../Tests/tests.c
>> >>
>> >> Please find attached a MWE (it is a slight modification of that of the
>> >> post opened by Ce Qin,
>> >>
>> https://lists.mcs.anl.gov/pipermail/petsc-users/2015-February/024230.html,
>>
>> >> whose answer I have not found).
>> >>
>> >> By the way, with the newest version of PETSc, Eclipse marks as errors
>> >> the commands PetscFree, CHKERRQ, PETSC_COMM_SELF,... although it
>> >> compiles and executes well. Perhaps it is a problem related to
>> Eclipse,
>> >> but this did not happen with the older versions of PETSc.
>> >>
>> >> Thanks and regards,
>> >>
>> >> Manuel
>> >>
>> >> ---
>> >>
>> >>
>> >>
>> >> --
>> >> 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/
>>
>>
>
> --
> 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/~knepley/>
>
>

-- 
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/~knepley/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20190411/9008045e/attachment-0001.html>