[petsc-users] Some questions regarding VecNest type.

[Vijay S. Mahadevan]

> vecnest.c:639 ?

Oops. Must have had a typo in my search.

> What do you mean by multiple indices? The union of multiple constituent
> index sets? You could make a hierarchical VecNest if you really want that
> layout, otherwise we need to find an inexpensive way to handle exotic index
> sets without losing structure.

Yes. I meant the union of indexes. i.e, if parent_vec = {child1,
child2, .. childn}, then VecGetSubVector([1,2]) would then return
subparent_vec={child1, child2} ? Not sure if I would call that exotic
but I feel this is natural since I think of the vector as truly
nested. If you are expecting the IS to contain [1..len(child1),
len(child1)+1..len(child1)+len(child2)], then the parent_vec isn't
truly nested but rather a regular parallel vector. Aren't these the
distinctions between -vec_type nest and otherwise ? This is also what
I gather from your explanation below. Please correct me if I am wrong.

> Since you can write your
> software so that this choice is a run-time option, I don't see a good
> rationale for choosing an interface that only does one thing.

Agreed. As I mentioned earlier, I originally did have both the options
but only compile time since I was using petsc-ext. Perhaps this is the
most optimal way forward.

I have just started playing around with the examples and hopefully the
tests will clarify some of my issues. Thanks.

On Wed, Oct 12, 2011 at 7:58 AM, Jed Brown <jedbrown at mcs.anl.gov> wrote:
> On Tue, Oct 11, 2011 at 15:53, Vijay S. Mahadevan <vijay.m at gmail.com> wrote:
>>
>> > It would make sense to have VecNestSetSubVec(), but it is not
>> > implemented
>> > yet.
>>
>> In src/vec/vec/examples/tests/ex37.c, I see a call to this function
>> but it is part of dead code.
>>
>> I don't see a specialization VecGetSubVector/VecRestoreSubVector in
>> VecNest
>
> vecnest.c:639 ?
>
>>
>> and so my related question is if I call VecGetSubVector with
>> multiple indices,
>
> What do you mean by multiple indices? The union of multiple constituent
> index sets? You could make a hierarchical VecNest if you really want that
> layout, otherwise we need to find an inexpensive way to handle exotic index
> sets without losing structure.
>
>>
>> do I expect a nested vector back ? And in the case
>> when the index set has only one entry, perhaps the returned vector is
>> in the parallel layout of the corresponding block entry. Or does the
>> returned Vector have a different parallel layout in these cases ?
>
> Same parallel layout, the returned vector has the distribution specified by
> the IS. The main difference between using VecNest and standard Vec is that
> the former can have a non-contiguous index set, but store it contiguously in
> memory, and provide no-copy access to that piece. With a normal Vec, you
> only get no-copy access if the index space is contiguous. (But usually you
> choose the ordering of the index space so that related things are nearby, so
> it should be a good ordering.)
>
>>
>> > You get a lot more consistent support if you
>> > use a standard contiguous Vec type and call VecGetSubVector() or a
>> > VecScatter when you need to get part of it.
>>
>> I used to take this route until I found the block extensions. I feel
>> it is elegant to handle different physics blocks using the block or
>> Nest interface specifically since dofs for each physics can be
>> distributed independently based on its own mesh and other
>> requirements. IMO, the blocking also makes the solve/preconditioning
>> for multi-physics problems natural (yes the sparsity pattern for
>> off-diagonal blocks are tricky). I will also investigate the Nest
>> independent interface to generalize this access.
>
> With the Nest-independent interface, you access pieces using an IS instead
> of an integer. I find it quite valuable when comparing algorithms to switch
> between monolithic formats with Schwarz/coupled multigrid/direct solvers and
> the Nest format with FieldSplit preconditioners. Since you can write your
> software so that this choice is a run-time option, I don't see a good
> rationale for choosing an interface that only does one thing.