[petsc-dev] MetSetBlockSize issue

Dmitry Karpeev karpeev at mcs.anl.gov
Wed Apr 25 16:33:56 CDT 2012 

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:
>
> >
> > On Apr 25, 2012, at 2:44 PM, Mark F. Adams wrote:
> >
> >>
> >> On Apr 25, 2012, at 12:22 PM, Hong Zhang wrote:
> >>
> >>> Mark,
> >>>
> >>> I attach it with this code:
> >>>
> >>> /* attach block size of columns */
> >>> if( pc_gamg->col_bs_id == -1 ) {
> >>>   ierr = PetscObjectComposedDataRegister( &pc_gamg->col_bs_id );
> assert(pc_gamg->col_bs_id != -1 );
> >>> }
> >>> ierr = PetscObjectComposedDataSetInt( (PetscObject)Prol,
> pc_gamg->col_bs_id, data_cols ); CHKERRQ(ierr);
> >>>
> >>> which Hong can not see w/o my ID.  so we need a mechanism for me to
> get this to PtAP.
> >>>
> >>> Sorry, I do not understand what do you mean here.
> >>> Do you want add above into the implementation of PtAP or petsc?
> >>
> >> This is an architectural/API decision.  I will do what you and Barry
> want.  You could add a parameter to PtAP, or add a field to Mat (column
> block size), or whatever.
> >
> >   Mark,
> >
> >     I'm leaning towards setting a column block size,
> MatSetColumnBlockSize(). But what exactly is a "column block size"?  Is it
> really a column blocksize in the sense that nbs columns next to each other
> always have the exact same nonzero pattern: for example
> >
> > xx  00   xx 00
> > yy  00   yy 00
> > zz  00  zz  00
>
> No.  The blocks are dense for smoothed aggregation but for geometric MG
> and classical AMG the blocks are scaled identity.
>
> >
> > 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.

>
> 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
> >>>>>>>>>>>>>
> >>>>>>>>>>>>>
> >>>>>>>>>>>>
> >>>>>>>>>>>
> >>>>>>>>>>>
> >>>>>>>>>>
> >>>>>>>>>
> >>>>>>>>>
> >>>>>>>>
> >>>>>>>>
> >>>>>>>
> >>>>>>
> >>>>>>
> >>>>>
> >>>>
> >>>>
> >>>
> >>>
> >>
> >
> >
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20120425/eb2beac5/attachment.html>

More information about the petsc-dev mailing list