Re: [petsc-dev] BlockGetIndices and GetBlockIndices

[Pierre Jolivet via petsc-dev]

> On 21 Oct 2019, at 7:52 AM, Smith, Barry F.  wrote:
> 
> 
> 
>> On Oct 21, 2019, at 12:23 AM, Pierre Jolivet  
>> wrote:
>> 
>> 
>> 
>>> On 21 Oct 2019, at 7:11 AM, Smith, Barry F.  wrote:
>>> 
>>> 
>>> 
 On Oct 20, 2019, at 11:52 PM, Pierre Jolivet  
 wrote:
 
 
 
> On 21 Oct 2019, at 6:42 AM, Smith, Barry F.  wrote:
> 
> Could you provide a use case where you want to access/have a block size 
> of a IS that is not an ISBlock? 
 
 In the end, all I really want is get access to the underlying 
 is->data->idx without having to worry about the subclass of is.
 I don’t have such a use case, but I don’t think this is really related to 
 what I want to achieve (or maybe it is…).
>>> 
>>> ISGetIndices()
>> 
>> Not true for ISBlock with bs > 1.
> 
>  Certainly suppose to be, is there a bug?
> 
> static PetscErrorCode ISGetIndices_Block(IS in,const PetscInt *idx[])
> {
>  IS_Block   *sub = (IS_Block*)in->data;
>  PetscErrorCode ierr;
>  PetscInt   i,j,k,bs,n,*ii,*jj;
> 
>  PetscFunctionBegin;
>  ierr = PetscLayoutGetBlockSize(in->map, &bs);CHKERRQ(ierr);
> 
> Dang, there is that stupid layout stuff again. Who put this crap in. 
> 
>  ierr = PetscLayoutGetLocalSize(in->map, &n);CHKERRQ(ierr);
>  n   /= bs;
>  if (bs == 1) *idx = sub->idx;
>  else {

There it is, I don’t want this if/else. ISGetBlockIndices would have been a 
function always returning sub->idx.

Thanks,
Pierre

>ierr = PetscMalloc1(bs*n,&jj);CHKERRQ(ierr);
>*idx = jj;
>k= 0;
>ii   = sub->idx;
>for (i=0; i  for (j=0; jjj[k++] = bs*ii[i] + j;
>  }
>  PetscFunctionReturn(0);
> }
> 
> 
>> 
 Again, my goal is just to avoid having to write something like: 
 https://www.mcs.anl.gov/petsc/petsc-dev/src/vec/is/utils/isltog.c.html#line322
 If things were done as in ISLocalToGlobalMapping (which I know is not the 
 same at all as IS/ISBlock), I think this could just read:
 ISGetBlockIndices(is,&indices)
 ISLocalToGlobalMappingCreate(comm,bs,n/bs,indices,PETSC_COPY_VALUES,mapping)
>>> 
>>> If you want your code to take advantage of the underlying structure of the 
>>> IS (be it strided, blocked) for optimizations then you HAVE to go through 
>>> the subclasses by design and nature (for example VecScatterCreate() takes 
>>> advantage of the structure of the IS as does 
>>> ISLocalToGlobalMappingCreate()). There is no universal concept of a block 
>>> size in IS (despite the crap that other people put in using PetscLayout 
>>> which has nothing to do with ISBlock and is completely buggy.) nor should 
>>> there be or is there a need for. 
>> 
>> OK then.
>> 
>> Thanks,
>> Pierre
>> 
>>> Barry
>>> 
>>> 
>>> 
>>> 
 
 Thanks,
 Pierre
 
> 
>> On Oct 16, 2019, at 2:50 AM, Pierre Jolivet via petsc-dev 
>>  wrote:
>> 
>> Hello,
>> I’m trying to understand what is the rationale for naming a function 
>> ISBlockGetIndices and another ISLocalToGlobalMappingGetBlockIndices 
>> (BlockGet vs. GetBlock).
>> Also, it looks to me like the implementation of ISBlockGetIndices is 
>> somehow less versatile than ISLocalToGlobalMappingGetBlockIndices.
>> Indeed, I may call ISLocalToGlobalMappingGetBlockIndices with an 
>> underlying non-block IS, while I can’t for ISBlockGetIndices, so a check 
>> must be performed by the user, e.g., 
>> https://www.mcs.anl.gov/petsc/petsc-current/src/vec/is/is/examples/tutorials/ex3.c.html#line58
>>  or this (IMHO) ugly code duplication 
>> https://www.mcs.anl.gov/petsc/petsc-dev/src/vec/is/utils/isltog.c.html#line322.
>> 
>> Thoughts and/or comments? Would it make sense to add an 
>> ISGetBlockIndices/ISRestoreBlockIndices or would that be too confusing 
>> for the user?
>> 
>> Thanks,
>> Pierre

Re: [petsc-dev] ksp_error_if_not_converged in multilevel solvers

[Pierre Jolivet via petsc-dev]

On Oct 20, 2019, at 6:07 PM, "Smith, Barry F."  wrote:

> 
>   The reason the code works this way is that normally 
> -ksp_error_if_not_converged is propagated into the inner (and innerer) solves 
> and normally it is desirable that these inner solves do not error simply 
> because they reach the maximum number of iterations since for nested 
> iterative methods generally we don't need or care if the inner solves 
> "converge". 

I fully agree with you on the last part of the above sentence. Thus, this makes 
me question the first part (which I wasn't aware of): why is 
error_if_not_converged being propagated to inner solves?
I'm sure there are good usages, but if one cares that ksp_1 (which depends on 
ksp_2) converges, why should an error be thrown if ksp_2 does not converge as 
long as ksp_1 does (I guess this goes along your last paragraph)?

Thanks,
Pierre

>Of course, in your particular use case this backfires. 
> 
>I'm not sure what the best design fix is. It seems anything will 
> complicate the code making it harder to maintain. 
> 
>Possible changes:
>   * if the error_if_not_converged is set directly on a KSP then also 
> cause it to error on maximum iterations, but not if it is propagated into a 
> KSP. This would require tracking how the error if not converged got set
>   * make error_if_not_converged an enum with, for example, the value 
> __always__ causing the error to be generated even on maximum iterations. This 
> is easier on the code then above 
> 
>A larger more far-reaching change would be to change the current model and 
> have the "inner" KSP use CONVERGED_ITERATING instead of DIVERGED_ITERATING 
> when it doesn't matter if they converge or not: 
> KSPSetReachMaximumIterationsConverged(), then the outer KSP would call this 
> on the inner KSP at construction time. This is appealing from code clarity 
> point of view since now the inner solves return a negative (diverged) reason 
> when that is acceptable and that is confusing, to return a diverged reason 
> while everything is fine. For the given use case one would need to do 
> -pc_bddc_coarse_ksp_error_if_not_converged  
> -pc_bddc_coarse_ksp_reach_maximum_iterations_converged false but this is 
> cumbersome, who would remember the second option.
> 
>   I'm not really excited by any of my proposed solutions.
> 
>  Thoughts?
> 
> 
>  Barry
> 
> 
> 
> 
>> On Oct 20, 2019, at 7:55 AM, Pierre Jolivet via petsc-dev 
>>  wrote:
>> 
>> Hello,
>> I’m trying to get multilevel solvers to error out when coarse levels are not 
>> converging, but I’m failing…
>> Could someone either tell me if this is not possible to do so, or help me 
>> find what the problem in my options is, please?
>> (in src/ksp/ksp/examples/tutorials)
>> $ mpirun -n 8 ./ex71 -pde_type Elasticity -cells 7,9,8 -dim 3 
>> -pc_bddc_monolithic -pc_bddc_use_faces -pc_bddc_coarse_pc_type none 
>> -pc_bddc_coarse_ksp_type gmres -pc_bddc_coarse_ksp_converged_reason 
>> -pc_bddc_coarse_ksp_error_if_not_converged -ksp_converged_reason
>>  Linear pc_bddc_coarse_ solve did not converge due to DIVERGED_ITS 
>> iterations 1 <- I want to error out here
>> […]
>>  Linear pc_bddc_coarse_ solve did not converge due to DIVERGED_ITS 
>> iterations 1
>> Linear solve converged due to CONVERGED_RTOL iterations 31
>> $ mpirun -np 1 ./ex1 -ksp_type gmres -pc_type gamg -ksp_converged_reason 
>> -mg_coarse_pc_type lu -mg_coarse_ksp_max_it 5 -mg_coarse_ksp_type gmres 
>> -ksp_type fgmres -mg_levels_1_ksp_type fgmres 
>> -mg_coarse_ksp_error_if_not_converged -mg_coarse_ksp_converged_reason
>>Linear mg_coarse_ solve did not converge due to DIVERGED_ITS iterations 5 
>>  <- I want to error out here
>> […]
>>Linear mg_coarse_ solve did not converge due to DIVERGED_ITS iterations 5
>> Linear solve converged due to CONVERGED_RTOL iterations 3
>> 
>> This is probably linked with this SETERRQ1 
>> https://www.mcs.anl.gov/petsc/petsc-dev/src/ksp/ksp/interface/itfunc.c.html#line836
>>  and the condition ksp->reason != KSP_DIVERGED_ITS.
>> If one just wants to run a fixed number of iterations, not checking for 
>> convergence, why would one set ksp->errorifnotconverged to true?
>> 
>> Thanks,
>> Pierre
> 

Re: [petsc-dev] BlockGetIndices and GetBlockIndices

[Smith, Barry F. via petsc-dev]

> On Oct 21, 2019, at 12:23 AM, Pierre Jolivet  
> wrote:
> 
> 
> 
>> On 21 Oct 2019, at 7:11 AM, Smith, Barry F.  wrote:
>> 
>> 
>> 
>>> On Oct 20, 2019, at 11:52 PM, Pierre Jolivet  
>>> wrote:
>>> 
>>> 
>>> 
 On 21 Oct 2019, at 6:42 AM, Smith, Barry F.  wrote:
 
 Could you provide a use case where you want to access/have a block size of 
 a IS that is not an ISBlock? 
>>> 
>>> In the end, all I really want is get access to the underlying is->data->idx 
>>> without having to worry about the subclass of is.
>>> I don’t have such a use case, but I don’t think this is really related to 
>>> what I want to achieve (or maybe it is…).
>> 
>>  ISGetIndices()
> 
> Not true for ISBlock with bs > 1.

  Certainly suppose to be, is there a bug?

static PetscErrorCode ISGetIndices_Block(IS in,const PetscInt *idx[])
{
  IS_Block   *sub = (IS_Block*)in->data;
  PetscErrorCode ierr;
  PetscInt   i,j,k,bs,n,*ii,*jj;

  PetscFunctionBegin;
  ierr = PetscLayoutGetBlockSize(in->map, &bs);CHKERRQ(ierr);

Dang, there is that stupid layout stuff again. Who put this crap in. 

  ierr = PetscLayoutGetLocalSize(in->map, &n);CHKERRQ(ierr);
  n   /= bs;
  if (bs == 1) *idx = sub->idx;
  else {
ierr = PetscMalloc1(bs*n,&jj);CHKERRQ(ierr);
*idx = jj;
k= 0;
ii   = sub->idx;
for (i=0; i 
>>> Again, my goal is just to avoid having to write something like: 
>>> https://www.mcs.anl.gov/petsc/petsc-dev/src/vec/is/utils/isltog.c.html#line322
>>> If things were done as in ISLocalToGlobalMapping (which I know is not the 
>>> same at all as IS/ISBlock), I think this could just read:
>>> ISGetBlockIndices(is,&indices)
>>> ISLocalToGlobalMappingCreate(comm,bs,n/bs,indices,PETSC_COPY_VALUES,mapping)
>> 
>> If you want your code to take advantage of the underlying structure of the 
>> IS (be it strided, blocked) for optimizations then you HAVE to go through 
>> the subclasses by design and nature (for example VecScatterCreate() takes 
>> advantage of the structure of the IS as does 
>> ISLocalToGlobalMappingCreate()). There is no universal concept of a block 
>> size in IS (despite the crap that other people put in using PetscLayout 
>> which has nothing to do with ISBlock and is completely buggy.) nor should 
>> there be or is there a need for. 
> 
> OK then.
> 
> Thanks,
> Pierre
> 
>> Barry
>> 
>> 
>> 
>> 
>>> 
>>> Thanks,
>>> Pierre
>>> 
 
> On Oct 16, 2019, at 2:50 AM, Pierre Jolivet via petsc-dev 
>  wrote:
> 
> Hello,
> I’m trying to understand what is the rationale for naming a function 
> ISBlockGetIndices and another ISLocalToGlobalMappingGetBlockIndices 
> (BlockGet vs. GetBlock).
> Also, it looks to me like the implementation of ISBlockGetIndices is 
> somehow less versatile than ISLocalToGlobalMappingGetBlockIndices.
> Indeed, I may call ISLocalToGlobalMappingGetBlockIndices with an 
> underlying non-block IS, while I can’t for ISBlockGetIndices, so a check 
> must be performed by the user, e.g., 
> https://www.mcs.anl.gov/petsc/petsc-current/src/vec/is/is/examples/tutorials/ex3.c.html#line58
>  or this (IMHO) ugly code duplication 
> https://www.mcs.anl.gov/petsc/petsc-dev/src/vec/is/utils/isltog.c.html#line322.
> 
> Thoughts and/or comments? Would it make sense to add an 
> ISGetBlockIndices/ISRestoreBlockIndices or would that be too confusing 
> for the user?
> 
> Thanks,
> Pierre

Re: [petsc-dev] BlockGetIndices and GetBlockIndices

[Pierre Jolivet via petsc-dev]

> On 21 Oct 2019, at 7:11 AM, Smith, Barry F.  wrote:
> 
> 
> 
>> On Oct 20, 2019, at 11:52 PM, Pierre Jolivet  
>> wrote:
>> 
>> 
>> 
>>> On 21 Oct 2019, at 6:42 AM, Smith, Barry F.  wrote:
>>> 
>>> Could you provide a use case where you want to access/have a block size of 
>>> a IS that is not an ISBlock? 
>> 
>> In the end, all I really want is get access to the underlying is->data->idx 
>> without having to worry about the subclass of is.
>> I don’t have such a use case, but I don’t think this is really related to 
>> what I want to achieve (or maybe it is…).
> 
>   ISGetIndices()

Not true for ISBlock with bs > 1.

>> Again, my goal is just to avoid having to write something like: 
>> https://www.mcs.anl.gov/petsc/petsc-dev/src/vec/is/utils/isltog.c.html#line322
>> If things were done as in ISLocalToGlobalMapping (which I know is not the 
>> same at all as IS/ISBlock), I think this could just read:
>> ISGetBlockIndices(is,&indices)
>> ISLocalToGlobalMappingCreate(comm,bs,n/bs,indices,PETSC_COPY_VALUES,mapping)
> 
>  If you want your code to take advantage of the underlying structure of the 
> IS (be it strided, blocked) for optimizations then you HAVE to go through the 
> subclasses by design and nature (for example VecScatterCreate() takes 
> advantage of the structure of the IS as does ISLocalToGlobalMappingCreate()). 
> There is no universal concept of a block size in IS (despite the crap that 
> other people put in using PetscLayout which has nothing to do with ISBlock 
> and is completely buggy.) nor should there be or is there a need for. 

OK then.

Thanks,
Pierre

>  Barry
> 
> 
> 
> 
>> 
>> Thanks,
>> Pierre
>> 
>>> 
 On Oct 16, 2019, at 2:50 AM, Pierre Jolivet via petsc-dev 
  wrote:
 
 Hello,
 I’m trying to understand what is the rationale for naming a function 
 ISBlockGetIndices and another ISLocalToGlobalMappingGetBlockIndices 
 (BlockGet vs. GetBlock).
 Also, it looks to me like the implementation of ISBlockGetIndices is 
 somehow less versatile than ISLocalToGlobalMappingGetBlockIndices.
 Indeed, I may call ISLocalToGlobalMappingGetBlockIndices with an 
 underlying non-block IS, while I can’t for ISBlockGetIndices, so a check 
 must be performed by the user, e.g., 
 https://www.mcs.anl.gov/petsc/petsc-current/src/vec/is/is/examples/tutorials/ex3.c.html#line58
  or this (IMHO) ugly code duplication 
 https://www.mcs.anl.gov/petsc/petsc-dev/src/vec/is/utils/isltog.c.html#line322.
 
 Thoughts and/or comments? Would it make sense to add an 
 ISGetBlockIndices/ISRestoreBlockIndices or would that be too confusing for 
 the user?
 
 Thanks,
 Pierre
>>> 
>> 
> 

Re: [petsc-dev] BlockGetIndices and GetBlockIndices

[Smith, Barry F. via petsc-dev]

> On Oct 20, 2019, at 11:52 PM, Pierre Jolivet  
> wrote:
> 
> 
> 
>> On 21 Oct 2019, at 6:42 AM, Smith, Barry F.  wrote:
>> 
>>  Could you provide a use case where you want to access/have a block size of 
>> a IS that is not an ISBlock? 
> 
> In the end, all I really want is get access to the underlying is->data->idx 
> without having to worry about the subclass of is.
> I don’t have such a use case, but I don’t think this is really related to 
> what I want to achieve (or maybe it is…).

   ISGetIndices()

> Again, my goal is just to avoid having to write something like: 
> https://www.mcs.anl.gov/petsc/petsc-dev/src/vec/is/utils/isltog.c.html#line322
> If things were done as in ISLocalToGlobalMapping (which I know is not the 
> same at all as IS/ISBlock), I think this could just read:
> ISGetBlockIndices(is,&indices)
> ISLocalToGlobalMappingCreate(comm,bs,n/bs,indices,PETSC_COPY_VALUES,mapping)

  If you want your code to take advantage of the underlying structure of the IS 
(be it strided, blocked) for optimizations then you HAVE to go through the 
subclasses by design and nature (for example VecScatterCreate() takes advantage 
of the structure of the IS as does ISLocalToGlobalMappingCreate()). There is no 
universal concept of a block size in IS (despite the crap that other people put 
in using PetscLayout which has nothing to do with ISBlock and is completely 
buggy.) nor should there be or is there a need for. 

  Barry




> 
> Thanks,
> Pierre
> 
>> 
>>> On Oct 16, 2019, at 2:50 AM, Pierre Jolivet via petsc-dev 
>>>  wrote:
>>> 
>>> Hello,
>>> I’m trying to understand what is the rationale for naming a function 
>>> ISBlockGetIndices and another ISLocalToGlobalMappingGetBlockIndices 
>>> (BlockGet vs. GetBlock).
>>> Also, it looks to me like the implementation of ISBlockGetIndices is 
>>> somehow less versatile than ISLocalToGlobalMappingGetBlockIndices.
>>> Indeed, I may call ISLocalToGlobalMappingGetBlockIndices with an underlying 
>>> non-block IS, while I can’t for ISBlockGetIndices, so a check must be 
>>> performed by the user, e.g., 
>>> https://www.mcs.anl.gov/petsc/petsc-current/src/vec/is/is/examples/tutorials/ex3.c.html#line58
>>>  or this (IMHO) ugly code duplication 
>>> https://www.mcs.anl.gov/petsc/petsc-dev/src/vec/is/utils/isltog.c.html#line322.
>>> 
>>> Thoughts and/or comments? Would it make sense to add an 
>>> ISGetBlockIndices/ISRestoreBlockIndices or would that be too confusing for 
>>> the user?
>>> 
>>> Thanks,
>>> Pierre
>> 
> 

Re: [petsc-dev] BlockGetIndices and GetBlockIndices

[Pierre Jolivet via petsc-dev]

> On 21 Oct 2019, at 6:42 AM, Smith, Barry F.  wrote:
> 
>  Could you provide a use case where you want to access/have a block size of a 
> IS that is not an ISBlock? 

In the end, all I really want is get access to the underlying is->data->idx 
without having to worry about the subclass of is.
I don’t have such a use case, but I don’t think this is really related to what 
I want to achieve (or maybe it is…).
Again, my goal is just to avoid having to write something like: 
https://www.mcs.anl.gov/petsc/petsc-dev/src/vec/is/utils/isltog.c.html#line322 

If things were done as in ISLocalToGlobalMapping (which I know is not the same 
at all as IS/ISBlock), I think this could just read:
ISGetBlockIndices(is,&indices)
ISLocalToGlobalMappingCreate(comm,bs,n/bs,indices,PETSC_COPY_VALUES,mapping)

Thanks,
Pierre

> 
>> On Oct 16, 2019, at 2:50 AM, Pierre Jolivet via petsc-dev 
>>  wrote:
>> 
>> Hello,
>> I’m trying to understand what is the rationale for naming a function 
>> ISBlockGetIndices and another ISLocalToGlobalMappingGetBlockIndices 
>> (BlockGet vs. GetBlock).
>> Also, it looks to me like the implementation of ISBlockGetIndices is somehow 
>> less versatile than ISLocalToGlobalMappingGetBlockIndices.
>> Indeed, I may call ISLocalToGlobalMappingGetBlockIndices with an underlying 
>> non-block IS, while I can’t for ISBlockGetIndices, so a check must be 
>> performed by the user, e.g., 
>> https://www.mcs.anl.gov/petsc/petsc-current/src/vec/is/is/examples/tutorials/ex3.c.html#line58
>>  or this (IMHO) ugly code duplication 
>> https://www.mcs.anl.gov/petsc/petsc-dev/src/vec/is/utils/isltog.c.html#line322.
>> 
>> Thoughts and/or comments? Would it make sense to add an 
>> ISGetBlockIndices/ISRestoreBlockIndices or would that be too confusing for 
>> the user?
>> 
>> Thanks,
>> Pierre
> 

Re: [petsc-dev] BlockGetIndices and GetBlockIndices

[Smith, Barry F. via petsc-dev]

  Could you provide a use case where you want to access/have a block size of a 
IS that is not an ISBlock? 

> On Oct 16, 2019, at 2:50 AM, Pierre Jolivet via petsc-dev 
>  wrote:
> 
> Hello,
> I’m trying to understand what is the rationale for naming a function 
> ISBlockGetIndices and another ISLocalToGlobalMappingGetBlockIndices (BlockGet 
> vs. GetBlock).
> Also, it looks to me like the implementation of ISBlockGetIndices is somehow 
> less versatile than ISLocalToGlobalMappingGetBlockIndices.
> Indeed, I may call ISLocalToGlobalMappingGetBlockIndices with an underlying 
> non-block IS, while I can’t for ISBlockGetIndices, so a check must be 
> performed by the user, e.g., 
> https://www.mcs.anl.gov/petsc/petsc-current/src/vec/is/is/examples/tutorials/ex3.c.html#line58
>  or this (IMHO) ugly code duplication 
> https://www.mcs.anl.gov/petsc/petsc-dev/src/vec/is/utils/isltog.c.html#line322.
> 
> Thoughts and/or comments? Would it make sense to add an 
> ISGetBlockIndices/ISRestoreBlockIndices or would that be too confusing for 
> the user?
> 
> Thanks,
> Pierre

Re: [petsc-dev] BlockGetIndices and GetBlockIndices

[Smith, Barry F. via petsc-dev]

> On Oct 16, 2019, at 9:41 AM, Stefano Zampini via petsc-dev 
>  wrote:
> 
> I just took a look at the ISGENERAL code. ISSetBlockSize_General just sets 
> the block size of the layout (??)
> ISGetIndices always return the data->idx memory.
> So, a more profound question is: what is the model behind setting the block 
> size on a ISGENERAL? And on a IS in general?

   I still don't understand why IS even has a layout so I was lost long ago.

   ISSetBlockSize_*() appear to be complete gibberish for ISGeneral, ISBlock, 
and ISStride! 

> 
>> On Oct 16, 2019, at 4:37 PM, Jed Brown  wrote:
>> 
>> Stefano Zampini via petsc-dev  writes:
>> 
 Thoughts and/or comments? Would it make sense to add an 
 ISGetBlockIndices/ISRestoreBlockIndices or would that be too confusing for 
 the user?

   Seems insane, ISGetBlockIndices/ISRestoreBlockIndices() should not only not 
exist the entire absurd SSetBlockSize() concept should be removed. 

   Barry


>>> 
>>> That would be more general and I think it makes sense, and should pair with 
>>>  ISGetBlockSize
>> 
>> What happens if you call ISGetBlockIndices on an ISGeneral for which
>> blocks are not contiguous?
> 

Re: [petsc-dev] BlockGetIndices and GetBlockIndices

[Smith, Barry F. via petsc-dev]

> On Oct 16, 2019, at 2:50 AM, Pierre Jolivet via petsc-dev 
>  wrote:
> 
> Hello,
> I’m trying to understand what is the rationale for naming a function 
> ISBlockGetIndices and another ISLocalToGlobalMappingGetBlockIndices (BlockGet 
> vs. GetBlock).

   ISBlockGetIndices returns the indices for each block; it is method that only 
applies too the  ISBlock subclass. Given its purpose its name is perfect.

   ISLocalToGlobalMapping is a class that has direct support for a block size 
(unlike IS which has no concept of a block size, it is an abstraction of any 
collection of integers).

   ISLocalToGlobalMappingGetBlockIndices is then a method on 
ISLocalToGlobalMapping that gives you the block indices (one index per block). 
Given its purpose its name is also perfect; the method is getblockindices()   
Note there is also ISLocalToGlobalMappingGetIndices(). 

  Though IS and ISLocalToGlobalMapping both begin with IS they are really not 
closely related.

   

> Also, it looks to me like the implementation of ISBlockGetIndices is somehow 
> less versatile than ISLocalToGlobalMappingGetBlockIndices.
> Indeed, I may call ISLocalToGlobalMappingGetBlockIndices with an underlying 
> non-block IS, while I can’t for ISBlockGetIndices, so a check must be 
> performed by the user, e.g., 
> https://www.mcs.anl.gov/petsc/petsc-current/src/vec/is/is/examples/tutorials/ex3.c.html#line58
>  or this (IMHO) ugly code duplication 
> https://www.mcs.anl.gov/petsc/petsc-dev/src/vec/is/utils/isltog.c.html#line322.
> 
> Thoughts and/or comments? Would it make sense to add an 
> ISGetBlockIndices/ISRestoreBlockIndices or would that be too confusing for 
> the user?

   For non ISBlock subclass of IS I don't know what ISGetBlockIndices is 
suppose to even mean. 

> 
> Thanks,
> Pierre

Re: [petsc-dev] ksp_error_if_not_converged in multilevel solvers

[Smith, Barry F. via petsc-dev]

   The reason the code works this way is that normally 
-ksp_error_if_not_converged is propagated into the inner (and innerer) solves 
and normally it is desirable that these inner solves do not error simply 
because they reach the maximum number of iterations since for nested iterative 
methods generally we don't need or care if the inner solves "converge". 

Of course, in your particular use case this backfires. 

I'm not sure what the best design fix is. It seems anything will complicate 
the code making it harder to maintain. 

Possible changes:
   * if the error_if_not_converged is set directly on a KSP then also cause 
it to error on maximum iterations, but not if it is propagated into a KSP. This 
would require tracking how the error if not converged got set
   * make error_if_not_converged an enum with, for example, the value 
__always__ causing the error to be generated even on maximum iterations. This 
is easier on the code then above 

A larger more far-reaching change would be to change the current model and 
have the "inner" KSP use CONVERGED_ITERATING instead of DIVERGED_ITERATING when 
it doesn't matter if they converge or not: 
KSPSetReachMaximumIterationsConverged(), then the outer KSP would call this on 
the inner KSP at construction time. This is appealing from code clarity point 
of view since now the inner solves return a negative (diverged) reason when 
that is acceptable and that is confusing, to return a diverged reason while 
everything is fine. For the given use case one would need to do 
-pc_bddc_coarse_ksp_error_if_not_converged  
-pc_bddc_coarse_ksp_reach_maximum_iterations_converged false but this is 
cumbersome, who would remember the second option.

   I'm not really excited by any of my proposed solutions.

  Thoughts?


  Barry




> On Oct 20, 2019, at 7:55 AM, Pierre Jolivet via petsc-dev 
>  wrote:
> 
> Hello,
> I’m trying to get multilevel solvers to error out when coarse levels are not 
> converging, but I’m failing…
> Could someone either tell me if this is not possible to do so, or help me 
> find what the problem in my options is, please?
> (in src/ksp/ksp/examples/tutorials)
> $ mpirun -n 8 ./ex71 -pde_type Elasticity -cells 7,9,8 -dim 3 
> -pc_bddc_monolithic -pc_bddc_use_faces -pc_bddc_coarse_pc_type none 
> -pc_bddc_coarse_ksp_type gmres -pc_bddc_coarse_ksp_converged_reason 
> -pc_bddc_coarse_ksp_error_if_not_converged -ksp_converged_reason
>   Linear pc_bddc_coarse_ solve did not converge due to DIVERGED_ITS 
> iterations 1 <- I want to error out here
> […]
>   Linear pc_bddc_coarse_ solve did not converge due to DIVERGED_ITS 
> iterations 1
> Linear solve converged due to CONVERGED_RTOL iterations 31
> $ mpirun -np 1 ./ex1 -ksp_type gmres -pc_type gamg -ksp_converged_reason 
> -mg_coarse_pc_type lu -mg_coarse_ksp_max_it 5 -mg_coarse_ksp_type gmres 
> -ksp_type fgmres -mg_levels_1_ksp_type fgmres 
> -mg_coarse_ksp_error_if_not_converged -mg_coarse_ksp_converged_reason
> Linear mg_coarse_ solve did not converge due to DIVERGED_ITS iterations 5 
>  <- I want to error out here
> […]
> Linear mg_coarse_ solve did not converge due to DIVERGED_ITS iterations 5
> Linear solve converged due to CONVERGED_RTOL iterations 3
> 
> This is probably linked with this SETERRQ1 
> https://www.mcs.anl.gov/petsc/petsc-dev/src/ksp/ksp/interface/itfunc.c.html#line836
>  and the condition ksp->reason != KSP_DIVERGED_ITS.
> If one just wants to run a fixed number of iterations, not checking for 
> convergence, why would one set ksp->errorifnotconverged to true?
> 
> Thanks,
> Pierre
> 

Re: [petsc-dev] ksp_error_if_not_converged in multilevel solvers

[Mark Adams via petsc-dev]

> If one just wants to run a fixed number of iterations, not checking for
> convergence, why would one set ksp->errorifnotconverged to true?
>
>
Good question. I can see not worrying too much about convergence on the
coarse grids, but to not allow it ... and now that I think about it, it
seems like we might want to error out with in indefinite PC. Maybe make
ksp->errorifnotconverged an int <1>:
0: no error
1: error if indefinite only
2: error if any error



> Thanks,
> Pierre
>
>

[petsc-dev] ksp_error_if_not_converged in multilevel solvers

[Pierre Jolivet via petsc-dev]

Hello,
I’m trying to get multilevel solvers to error out when coarse levels are not 
converging, but I’m failing…
Could someone either tell me if this is not possible to do so, or help me find 
what the problem in my options is, please?
(in src/ksp/ksp/examples/tutorials)
$ mpirun -n 8 ./ex71 -pde_type Elasticity -cells 7,9,8 -dim 3 
-pc_bddc_monolithic -pc_bddc_use_faces -pc_bddc_coarse_pc_type none 
-pc_bddc_coarse_ksp_type gmres -pc_bddc_coarse_ksp_converged_reason 
-pc_bddc_coarse_ksp_error_if_not_converged -ksp_converged_reason
  Linear pc_bddc_coarse_ solve did not converge due to DIVERGED_ITS iterations 
1 <- I want to error out here
[…]
  Linear pc_bddc_coarse_ solve did not converge due to DIVERGED_ITS iterations 1
Linear solve converged due to CONVERGED_RTOL iterations 31
$ mpirun -np 1 ./ex1 -ksp_type gmres -pc_type gamg -ksp_converged_reason 
-mg_coarse_pc_type lu -mg_coarse_ksp_max_it 5 -mg_coarse_ksp_type gmres 
-ksp_type fgmres -mg_levels_1_ksp_type fgmres 
-mg_coarse_ksp_error_if_not_converged -mg_coarse_ksp_converged_reason
Linear mg_coarse_ solve did not converge due to DIVERGED_ITS iterations 5  
<- I want to error out here
[…]
Linear mg_coarse_ solve did not converge due to DIVERGED_ITS iterations 5
Linear solve converged due to CONVERGED_RTOL iterations 3

This is probably linked with this SETERRQ1 
https://www.mcs.anl.gov/petsc/petsc-dev/src/ksp/ksp/interface/itfunc.c.html#line836
 

 and the condition ksp->reason != KSP_DIVERGED_ITS.
If one just wants to run a fixed number of iterations, not checking for 
convergence, why would one set ksp->errorifnotconverged to true?

Thanks,
Pierre