[petsc-users] Create a nest not aligned by processors
Matthew Knepley
knepley at gmail.com
Mon Mar 20 09:51:11 CDT 2023
On Mon, Mar 20, 2023 at 10:09 AM Berger Clement <clement.berger at ens-lyon.fr>
wrote:
> Ok so this means that if I define my vectors via VecNest (organized as the
> matrices), everything will be correctly ordered ?
>
Yes.
> How does that behave with MatConvert ? In the sense that if I convert a
> MatNest to MatAIJ via MatConvert, will the vector as built in the example I
> showed you work properly ?
>
> No. MatConvert just changes the storage format, not the ordering.
Thanks,
Matt
> Thank you !
> ---
> Clément BERGER
> ENS de Lyon
>
>
> Le 2023-03-20 14:58, Matthew Knepley a écrit :
>
> On Mon, Mar 20, 2023 at 9:35 AM Berger Clement <clement.berger at ens-lyon.fr>
> wrote:
>
>> 1) Yes in my program I use PETSC_DETERMINE, but I don't see what is the
>> issue there. From what I understand, it just lets PETSc set the total size
>> from the local sizes provided, am I mistaken ?
>>
>> 2) I attached a small script, when I run it with 1 proc the output vector
>> is not the same as if I run it with 2 procs, I don't know what I should do
>> to make them match.
>>
>> PS : I precise that I am not trying to point out a bug here, I realize
>> that my usage is wrong somehow, I just can't determine why, sorry if I gave
>> you the wrong impression !
>>
>>
>> I think I can now explain this clearly. Thank you for the nice simple
> example. I attach my slightly changed version (I think better in C). Here
> is running on one process:
>
> master *:~/Downloads/tmp/Berger$ /PETSc3/petsc/apple/bin/mpiexec -n 1
> ./nestTest -left_view -right_view -nest_view -full_view
> Mat Object: 1 MPI process
> type: nest
> Matrix object:
> type=nest, rows=2, cols=2
> MatNest structure:
> (0,0) : type=constantdiagonal, rows=4, cols=4
> (0,1) : NULL
> (1,0) : NULL
> (1,1) : type=constantdiagonal, rows=4, cols=4
> Mat Object: 1 MPI process
> type: seqaij
> row 0: (0, 2.)
> row 1: (1, 2.)
> row 2: (2, 2.)
> row 3: (3, 2.)
> row 4: (4, 1.)
> row 5: (5, 1.)
> row 6: (6, 1.)
> row 7: (7, 1.)
> Vec Object: 1 MPI process
> type: seq
> 0.
> 1.
> 2.
> 3.
> 4.
> 5.
> 6.
> 7.
> Vec Object: 1 MPI process
> type: seq
> 0.
> 2.
> 4.
> 6.
> 4.
> 5.
> 6.
> 7.
>
> This looks like what you expect. Doubling the first four rows and
> reproducing the last four. Now let's run on two processes:
>
> master *:~/Downloads/tmp/Berger$ /PETSc3/petsc/apple/bin/mpiexec -n 2
> ./nestTest -left_view -right_view -nest_view -full_view
> Mat Object: 2 MPI processes
> type: nest
> Matrix object:
> type=nest, rows=2, cols=2
> MatNest structure:
> (0,0) : type=constantdiagonal, rows=4, cols=4
> (0,1) : NULL
> (1,0) : NULL
> (1,1) : type=constantdiagonal, rows=4, cols=4
> Mat Object: 2 MPI processes
> type: mpiaij
> row 0: (0, 2.)
> row 1: (1, 2.)
> row 2: (2, 1.)
> row 3: (3, 1.)
> row 4: (4, 2.)
> row 5: (5, 2.)
> row 6: (6, 1.)
> row 7: (7, 1.)
> Vec Object: 2 MPI processes
> type: mpi
> Process [0]
> 0.
> 1.
> 2.
> 3.
> Process [1]
> 4.
> 5.
> 6.
> 7.
> Vec Object: 2 MPI processes
> type: mpi
> Process [0]
> 0.
> 2.
> 2.
> 3.
> Process [1]
> 8.
> 10.
> 6.
> 7.
>
> Let me describe what has changed. The matrices A and B are parallel, so
> each has two rows on process 0 and two rows on process 1. In the MatNest
> they are interleaved because we asked for contiguous numbering (by giving
> NULL for the IS of global row numbers). If we want to reproduce the same
> output, we would need to produce our input vector with the same interleaved
> numbering.
>
> Thanks,
>
> Matt
>
>> Thank you,
>> ---
>> Clément BERGER
>> ENS de Lyon
>>
>>
>> Le 2023-03-20 12:53, Matthew Knepley a écrit :
>>
>> On Mon, Mar 20, 2023 at 6:18 AM Berger Clement <
>> clement.berger at ens-lyon.fr> wrote:
>>
>>> I simplified the problem with the initial test I talked about because I
>>> thought I identified the issue, so I will walk you through my whole problem
>>> :
>>>
>>> - first the solve doesn't produce the same results as mentioned
>>>
>>> - I noticed that the duration of the factorization step of the matrix
>>> was also not consistent with the number of processors used (it is longer
>>> with 3 processes than with 1), I didn't think much of it but I now realize
>>> that for instance with 4 processes, MUMPS crashes when factorizing
>>>
>>> - I thought my matrices were wrong, but it's hard for me to use MatView
>>> to compare them with 1 or 2 proc because I work with a quite specific
>>> geometry, so in order not to fall into some weird particular case I need to
>>> use at least roughly 100 points, so looking at 100x100 matrices is not
>>> really nice...Instead I tried to multiply them by a vector full of one
>>> (after I used the vector v such that v(i)=i). I tried it on two matrices,
>>> and the results didn't depend on the number of procs, but when I tried to
>>> multiply against the nest of these two matrices (a 2x2 block diagonal
>>> nest), the result changed depending on the number of processors used
>>>
>>> - that's why I tried the toy problem I wrote to you in the first place
>>>
>>> I hope it's clearer now.
>>>
>>
>> Unfortunately, it is not clear to me. There is nothing attached to this
>> email. I will try to describe things from my end.
>>
>> 1) There are lots of tests. Internally, Nest does not depend on the
>> number of processes unless you make it so. This leads
>> me to believe that your construction of the matrix changes with the
>> number of processes. For example, using PETSC_DETERMINE
>> for sizes will do this.
>>
>> 2) In order to understand what you want to achieve, we need to have
>> something running in two cases, one with "correct" output and one
>> with something different. It sounds like you have such a small
>> example, but I have missed it.
>>
>> Can you attach this example? Then I can run it, look at the matrices, and
>> see what is different.
>>
>> Thanks,
>>
>> Matt
>>
>>
>>> Thank you
>>> ---
>>> Clément BERGER
>>> ENS de Lyon
>>>
>>>
>>> Le 2023-03-17 21:57, Barry Smith a écrit :
>>>
>>>
>>> Yes, you would benefit from a VecConvert() to produce a standard
>>> vector. But you should be able to use VecGetArray() on the nest array and
>>> on the standard array and copy the values between the arrays any way you
>>> like. You don't need to do any reordering when you copy. Is that not
>>> working and what are the symptoms (more than just the answers to the linear
>>> solve are different)? Again you can run on one and two MPI processes with a
>>> tiny problem to see if things are not in the correct order in the vectors
>>> and matrices.
>>>
>>> Barry
>>>
>>>
>>> On Mar 17, 2023, at 3:22 PM, Berger Clement <clement.berger at ens-lyon.fr>
>>> wrote:
>>>
>>> To use MUMPS I need to convert my matrix in MATAIJ format (or at least
>>> not MATNEST), after that if I use a VECNEST for the left and right hanside,
>>> I get an error during the solve procedure, it is removed if I copy my data
>>> in a vector with standard format, I couldn't find any other way
>>> ---
>>> Clément BERGER
>>> ENS de Lyon
>>>
>>>
>>> Le 2023-03-17 19:53, Matthew Knepley a écrit :
>>>
>>> On Fri, Mar 17, 2023 at 2:53 PM Berger Clement <
>>> clement.berger at ens-lyon.fr> wrote:
>>>
>>>> But this is to properly fill up the VecNest am I right ? Because this
>>>> one is correct, but I can't directly use it in the KSPSolve, I need to copy
>>>> it into a standard vector
>>>>
>>>>
>>> I do not understand what you mean here. You can definitely use a
>>> VecNest in a KSP.
>>>
>>> Thanks,
>>>
>>> Matt
>>>
>>>
>>>
>>>> ---
>>>> Clément BERGER
>>>> ENS de Lyon
>>>>
>>>>
>>>> Le 2023-03-17 19:39, Barry Smith a écrit :
>>>>
>>>>
>>>> I think the intention is that you use VecNestGetSubVecs()
>>>> or VecNestGetSubVec() and fill up the sub-vectors in the same style as the
>>>> matrices; this decreases the change of a reordering mistake in trying to do
>>>> it by hand in your code.
>>>>
>>>>
>>>>
>>>> On Mar 17, 2023, at 2:35 PM, Berger Clement <clement.berger at ens-lyon.fr>
>>>> wrote:
>>>>
>>>> That might be it, I didn't find the equivalent of MatConvert for the
>>>> vectors, so when I need to solve my linear system, with my righthandside
>>>> properly computed in nest format, I create a new vector using VecDuplicate,
>>>> and then I copy into it my data using VecGetArrayF90 and copiing each
>>>> element by hand. Does it create an incorrect ordering ? If so how can I get
>>>> the correct one ?
>>>> ---
>>>> Clément BERGER
>>>> ENS de Lyon
>>>>
>>>>
>>>> Le 2023-03-17 19:27, Barry Smith a écrit :
>>>>
>>>>
>>>> I would run your code with small sizes on 1, 2, 3 MPI ranks and use
>>>> MatView() to examine the matrices. They will definitely be ordered
>>>> differently but should otherwise be the same. My guess is that the right
>>>> hand side may not have the correct ordering with respect to the matrix
>>>> ordering in parallel. Note also that when the right hand side does have the
>>>> correct ordering the solution will have a different ordering for each
>>>> different number of MPI ranks when printed (but changing the ordering
>>>> should give the same results up to machine precision.
>>>>
>>>> Barry
>>>>
>>>>
>>>> On Mar 17, 2023, at 2:23 PM, Berger Clement <clement.berger at ens-lyon.fr>
>>>> wrote:
>>>>
>>>> My issue is that it seems to improperly with some step of my process,
>>>> the solve step doesn't provide the same result depending on the number of
>>>> processors I use. I manually tried to multiply one the matrices I defined
>>>> as a nest against a vector, and the result is not the same with e.g. 1 and
>>>> 3 processors. That's why I tried the toy program I wrote in the first
>>>> place, which highlights the misplacement of elements.
>>>> ---
>>>> Clément BERGER
>>>> ENS de Lyon
>>>>
>>>>
>>>> Le 2023-03-17 19:14, Barry Smith a écrit :
>>>>
>>>>
>>>> This sounds like a fine use of MATNEST. Now back to the original
>>>> question
>>>>
>>>>
>>>> I want to construct a matrix by blocs, each block having different
>>>> sizes and partially stored by multiple processors. If I am not mistaken,
>>>> the right way to do so is by using the MATNEST type. However, the following
>>>> code
>>>>
>>>> Call
>>>> MatCreateConstantDiagonal(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,4,4,2.0E0_wp,A,ierr)
>>>> Call
>>>> MatCreateConstantDiagonal(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,4,4,1.0E0_wp,B,ierr)
>>>> Call
>>>> MatCreateNest(PETSC_COMM_WORLD,2,PETSC_NULL_INTEGER,2,PETSC_NULL_INTEGER,(/A,PETSC_NULL_MAT,PETSC_NULL_MAT,B/),C,ierr)
>>>>
>>>> does not generate the same matrix depending on the number of
>>>> processors. It seems that it starts by everything owned by the first proc
>>>> for A and B, then goes on to the second proc and so on (I hope I am being
>>>> clear).
>>>>
>>>> Is it possible to change that ?
>>>>
>>>> If I understand correctly it is behaving as expected. It is the same
>>>> matrix on 1 and 2 MPI processes, the only difference is the ordering of the
>>>> rows and columns.
>>>>
>>>> Both matrix blocks are split among the two MPI processes. This is how
>>>> MATNEST works and likely what you want in practice.
>>>>
>>>> On Mar 17, 2023, at 1:19 PM, Berger Clement <clement.berger at ens-lyon.fr>
>>>> wrote:
>>>>
>>>> I have a matrix with four different blocks (2rows - 2columns). The
>>>> block sizes differ from one another, because they correspond to a different
>>>> physical variable. One of the block has the particularity that it has to be
>>>> updated at each iteration. This update is performed by replacing it with a
>>>> product of multiple matrices that depend on the result of the previous
>>>> iteration. Note that these intermediate matrices are not square (because
>>>> they also correspond to other types of variables), and that they must be
>>>> completely refilled by hand (i.e. they are not the result of some simple
>>>> linear operations). Finally, I use this final block matrix to solve
>>>> multiple linear systems (with different righthand sides), so for now I use
>>>> MUMPS as only the first solve takes time (but I might change it).
>>>>
>>>> Considering this setting, I created each type of variable separately,
>>>> filled the different matrices, and created different nests of vectors /
>>>> matrices for my operations. When the time comes to use KSPSolve, I use
>>>> MatConvert on my matrix to get a MATAIJ compatible with MUMPS, I also copy
>>>> the few vector data I need from my nests in a regular Vector, I solve, I
>>>> get back my data in my nest and carry on with the operations needed for my
>>>> updates.
>>>>
>>>> Is that clear ? I don't know if I provided too many or not enough
>>>> details.
>>>>
>>>> Thank you
>>>> ---
>>>> Clément BERGER
>>>> ENS de Lyon
>>>>
>>>>
>>>> Le 2023-03-17 17:34, Barry Smith a écrit :
>>>>
>>>>
>>>> Perhaps if you provide a brief summary of what you would like to do
>>>> and we may have ideas on how to achieve it.
>>>>
>>>> Barry
>>>>
>>>> Note: that MATNEST does require that all matrices live on all the MPI
>>>> processes within the original communicator. That is if the original
>>>> communicator has ranks 0,1, and 2 you cannot have a matrix inside MATNEST
>>>> that only lives on ranks 1,2 but you could have it have 0 rows on rank zero
>>>> so effectively it lives only on rank 1 and 2 (though its communicator is
>>>> all three ranks).
>>>>
>>>> On Mar 17, 2023, at 12:14 PM, Berger Clement <
>>>> clement.berger at ens-lyon.fr> wrote:
>>>>
>>>> It would be possible in the case I showed you but in mine that would
>>>> actually be quite complicated, isn't there any other workaround ? I precise
>>>> that I am not entitled to utilizing the MATNEST format, it's just that I
>>>> think the other ones wouldn't work.
>>>> ---
>>>> Clément BERGER
>>>> ENS de Lyon
>>>>
>>>>
>>>> Le 2023-03-17 15:48, Barry Smith a écrit :
>>>>
>>>>
>>>> You may be able to mimic what you want by not using PETSC_DECIDE but
>>>> instead computing up front how many rows of each matrix you want stored on
>>>> each MPI process. You can use 0 for on certain MPI processes for certain
>>>> matrices if you don't want any rows of that particular matrix stored on
>>>> that particular MPI process.
>>>>
>>>> Barry
>>>>
>>>>
>>>> On Mar 17, 2023, at 10:10 AM, Berger Clement <
>>>> clement.berger at ens-lyon.fr> wrote:
>>>>
>>>> Dear all,
>>>>
>>>> I want to construct a matrix by blocs, each block having different
>>>> sizes and partially stored by multiple processors. If I am not mistaken,
>>>> the right way to do so is by using the MATNEST type. However, the following
>>>> code
>>>>
>>>> Call
>>>> MatCreateConstantDiagonal(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,4,4,2.0E0_wp,A,ierr)
>>>> Call
>>>> MatCreateConstantDiagonal(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,4,4,1.0E0_wp,B,ierr)
>>>> Call
>>>> MatCreateNest(PETSC_COMM_WORLD,2,PETSC_NULL_INTEGER,2,PETSC_NULL_INTEGER,(/A,PETSC_NULL_MAT,PETSC_NULL_MAT,B/),C,ierr)
>>>>
>>>> does not generate the same matrix depending on the number of
>>>> processors. It seems that it starts by everything owned by the first proc
>>>> for A and B, then goes on to the second proc and so on (I hope I am being
>>>> clear).
>>>>
>>>> Is it possible to change that ?
>>>>
>>>> Note that I am coding in fortran if that has ay consequence.
>>>>
>>>> Thank you,
>>>>
>>>> Sincerely,
>>>> --
>>>> Clément BERGER
>>>> ENS de Lyon
>>>>
>>>>
>>>
>>> --
>>> 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/>
>>
>>
>
> --
> 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/20230320/51bb65a1/attachment-0001.html>
More information about the petsc-users
mailing list