[petsc-dev] MetSetBlockSize issue

Wed Apr 25 16:56:03 CDT 2012

On Apr 25, 2012, at 5:41 PM, Barry Smith wrote:

> 
> On Apr 25, 2012, at 4:33 PM, Dmitry Karpeev wrote:
> 
>> 
>> 
>> On Wed, Apr 25, 2012 at 4:09 PM, Mark F. Adams <mark.adams at columbia.edu> wrote:
>> 
>> On Apr 25, 2012, at 3:48 PM, Barry Smith wrote:
>> 
>>> or is it a more vague thing that defined by "well if you compute Pt A P the result has a blocksize of nbs".  In that case we would provide a MatSetVirtualColumnBlockSize()
>>> 
>> 
>> All I really need is for the coarse grid matrix to have the block size, because that is used to create a graph (for "vertices" not equations).  The method that constructs the prolongator knows the block size, so I attach it to P.  I could actually attach the block size to the coarse matrix in the way that I do with P, or pass it around as a parameter internally, I suppose.
>> 
>> Actually, in my opinion, this is the cleanest solution: attach the column block size to P and then to Ac = PtAP.
>> All of this currently happens inside GAMG, which has access to the composed datum's id.
>> Otherwise we are expanding the API with MatSetVirtualColumnSize(), which basically has no meaning outside GAMG.
>> I think the data composition mechanism exists exactly for these specific cases.
>> 
>> Dmitry.
> 
>    Dmitry, 
> 
>      The problem with that approach is that the routine that forms PtAP is a generic routine, (nothing to do with GAMG) and so it cannot be getting GAMG specific composed datuums out of the P.  
> 
>      I do not want the block size passed as an argument to PtAP() because that is a real ugly interface and the block size information arises from the matrix P and A, not from some arbitrary value one could pass into PtAp().
> 

Also, Dmitry,

The coarse grid matrices do have this block data structure.  So I can survive by doing it all in GAMG as I said but it is bad because then Mat methods, that can be faster with block size information, will not be able to exploit these optimizations.

Mark

> 
>    Barry
> 
>  I admit MatSetVirtualColumnBlockSize() is clunky and would be happy to hear alternatives.
> 
> 
> 
>> 
>> Mark
>> 
>>>   Thanks
>>> 
>>>  Barry
>>> 
>>>> 
>>>> Mark
>>>> 
>>>>> Hong
>>>>> 
>>>>> And MatGetSubMatrix has to inherit block size, which I assume is trial.
>>>>> 
>>>>> Mark
>>>>> 
>>>>> On Apr 24, 2012, at 11:14 PM, Barry Smith wrote:
>>>>> 
>>>>>> 
>>>>>> On Apr 24, 2012, at 8:14 PM, Mark F. Adams wrote:
>>>>>> 
>>>>>>> 
>>>>>>> On Apr 24, 2012, at 5:01 PM, Barry Smith wrote:
>>>>>>> 
>>>>>>>> 
>>>>>>>> On Apr 24, 2012, at 2:38 PM, Mark F. Adams wrote:
>>>>>>>> 
>>>>>>>>> 
>>>>>>>>> On Apr 24, 2012, at 3:28 PM, Hong Zhang wrote:
>>>>>>>>> 
>>>>>>>>>> Mark :
>>>>>>>>>> Shall C=PtAP inherit the block size of A?
>>>>>>>>>> Currently, MatPtAP() is only implemented with aij bs=1.
>>>>>>>>> 
>>>>>>>>> No, as I said, algebraically it should inherit the column block size of P, but that is not available.
>>>>>>>> 
>>>>>>>> What is the column block size of P and how do you know what it is?
>>>>>>> 
>>>>>>> I now attach it to the P matrix instead of adding it as a return parameter of the create P method.  So I keep track of it in the code.  Its size is the number of null space vectors.
>>>>>> 
>>>>>> Ok if P knows its block size then the code that computes PtAP can set the block size of the resulting matrix as it creates it using the information in P. So is the issue on calling MatSetBlockSize() now resolved and everyone happy?
>>>>>> 
>>>>>> Barry
>>>>>> 
>>>>>>> 
>>>>>>> Mark
>>>>>>> 
>>>>>>>> 
>>>>>>>> Barry
>>>>>>>> 
>>>>>>>>> 
>>>>>>>>> Perhaps PtAP should take a parameter for the block size ...
>>>>>>>>> 
>>>>>>>>> Mark
>>>>>>>>> 
>>>>>>>>>> Hong
>>>>>>>>>> 
>>>>>>>>>> On Apr 24, 2012, at 2:55 PM, Barry Smith wrote:
>>>>>>>>>> 
>>>>>>>>>>> 
>>>>>>>>>>> On Apr 24, 2012, at 1:39 PM, Mark F. Adams wrote:
>>>>>>>>>>> 
>>>>>>>>>>>> Now I'm getting this error with code like this:
>>>>>>>>>>>> 
>>>>>>>>>>>> ierr = MatGetSubMatrix( Cmat, new_eq_indices, new_eq_indices, MAT_INITIAL_MATRIX, &mat );
>>>>>>>>>>>> CHKERRQ(ierr);
>>>>>>>>>>>> ierr = MatSetBlockSize( mat, cbs );      CHKERRQ(ierr);
>>>>>>>>>>>> 
>>>>>>>>>>>> and like this:
>>>>>>>>>>>> 
>>>>>>>>>>>> ierr = MatPtAP( Amat_fine, Pold, MAT_INITIAL_MATRIX, 2.0, &Cmat ); CHKERRQ(ierr);
>>>>>>>>>>>> ierr = MatSetBlockSize( Cmat, cbs );      CHKERRQ(ierr);
>>>>>>>>>>> 
>>>>>>>>>>> Right you cannot do this. The matrix already exists so you cannot now set its size.
>>>>>>>>>>> 
>>>>>>>>>>> Does Cmat have a block size of cbs
>>>>>>>>>>> 
>>>>>>>>>>> What about Amat?   By default these routines should create the new matrix with the correct blocksize.
>>>>>>>>>>> 
>>>>>>>>>>> The harder part is if the blocksize of mat would be different than Cmat?
>>>>>>>>>>> 
>>>>>>>>>> 
>>>>>>>>>> It seems like MatGetSubMatrix should just inherit the block size but PtAP is harder.  The block size of the Cmat is the column block size of P.  But I can not set a column block size (!= row block size) so I don't set block size on P at all.  Here cbs is this column block size of P, or what it should be.
>>>>>>>>>> 
>>>>>>>>>> Mark
>>>>>>>>>> 
>>>>>>>>>>> Barry
>>>>>>>>>>> 
>>>>>>>>>>>> 
>>>>>>>>>>>> Mark
>>>>>>>>>>>> 
>>>>>>>>>>>> On Apr 23, 2012, at 9:09 PM, Barry Smith wrote:
>>>>>>>>>>>> 
>>>>>>>>>>>>> 
>>>>>>>>>>>>> Yes, look, for example how ex32 is handled at the bottom of the makefile.
>>>>>>>>>>>>> 
>>>>>>>>>>>>> Thanks
>>>>>>>>>>>>> 
>>>>>>>>>>>>> Barry
>>>>>>>>>>>>> 
>>>>>>>>>>>>> On Apr 23, 2012, at 5:50 PM, Mark F. Adams wrote:
>>>>>>>>>>>>> 
>>>>>>>>>>>>>> 
>>>>>>>>>>>>>> On Apr 23, 2012, at 5:59 PM, Barry Smith wrote:
>>>>>>>>>>>>>> 
>>>>>>>>>>>>>>> 
>>>>>>>>>>>>>>> It was changed a while ago that MatSetBlockSize() couldn't be set after the matrix was full instantiated. I guess that example did not get fixed because it is not listed in the makefile to run in the makeall.
>>>>>>>>>>>>>> 
>>>>>>>>>>>>>> May I add it?
>>>>>>>>>>>>>> 
>>>>>>>>>>>>>>> 
>>>>>>>>>>>>>>> I have fixed the example to call MatSetBlockSize() at a safe point.
>>>>>>>>>>>>>>> 
>>>>>>>>>>>>>>> Barry
>>>>>>>>>>>>>>> 
>>>>>>>>>>>>>>> On Apr 23, 2012, at 4:29 PM, Mark F. Adams wrote:
>>>>>>>>>>>>>>> 
>>>>>>>>>>>>>>>> ex55.c in ksp is failing with:
>>>>>>>>>>>>>>>> 
>>>>>>>>>>>>>>>> [0]PETSC ERROR: --------------------- Error Message ------------------------------------
>>>>>>>>>>>>>>>> [0]PETSC ERROR: Arguments are incompatible!
>>>>>>>>>>>>>>>> [0]PETSC ERROR: Cannot change block size 1 to 2!
>>>>>>>>>>>>>>>> [0]PETSC ERROR: ------------------------------------------------------------------------
>>>>>>>>>>>>>>>> 
>>>>>>>>>>>>>>>> on this line 57 of ex55.c:
>>>>>>>>>>>>>>>> 
>>>>>>>>>>>>>>>> ierr = MatSetBlockSize(Amat,2);      CHKERRQ(ierr);
>>>>>>>>>>>>>>>> 
>>>>>>>>>>>>>>>> Any idea what happened here?
>>>>>>>>>>>>>>>> 
>>>>>>>>>>>>>>>> Mark
>>>>>>>>>>>>>>> 
>>>>>>>>>>>>>>> 
>>>>>>>>>>>>>> 
>>>>>>>>>>>>> 
>>>>>>>>>>>>> 
>>>>>>>>>>>> 
>>>>>>>>>>> 
>>>>>>>>>>> 
>>>>>>>>>> 
>>>>>>>>>> 
>>>>>>>>> 
>>>>>>>> 
>>>>>>>> 
>>>>>>> 
>>>>>> 
>>>>>> 
>>>>> 
>>>>> 
>>>> 
>>> 
>>> 
>> 
>> 
> 
> 