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>