On Tue, 22 Jun 2010 09:41:15 -0400, Valerio Grazioso <graziosov at me.queensu.ca> wrote: > Hi Jed, > > On 2010-06-22, at 6:03 AM, Jed Brown wrote: > > > On Tue, 22 Jun 2010 04:02:29 -0400, Valerio Grazioso <graziosov at > > me.queensu.ca> wrote: > >> and I've got a matrix A ordered in a Petsc ordering: rows and variables > >> have a sequential numbering relative to each processor. > > > > This is always the order used internally, but your application code > > doesn't have to "see" this ordering. I suggest using > > MatSetValuesStencil() which allows you to effectively insert values in > > the natural ordering (as coordinates k,j,i) which will get mapped > > appropriately to the internal storage. > > Yes I've used MatSetValuesStencil() to insert the values and it worked > perfectly. > When I say that I've got a matrix in Petsc ordering is > because when I see it with -mat_view (at run time) that is what is printed at > screen.
Yeah, that's PETSc's ordering (the one used internally). There isn't much point outputting a matrix in a different ordering, it's just nice to assemble it using natural indexing (MatSetValuesStencil). > Also in this case, when I say that I get a naturally ordered vector is > because when I use > > > ..... > > call DALocalToGlobal(da,qloc,INSERT_VALUES,q,err) > call vecView(q,PETSC_VIEWER_STDOUT_WORLD,err) > > ..... > > that what is printed at screen (a naturally ordered vector). This is because it's often convenient to have a vector in the natural format (because you might read it in on a different number of processes). This is only done for DA vectors, you can PetscViewerPushFormat(viewer,PETSC_VIEWER_NATIVE) to write it in PETSc's ordering (consistent with the matrix). > But if I'm understanding well what you suggest, at the end, no matter > what I see printed at screen (for the matrix), I should be already > working in natural ordering for both the matrix and the rhs? This is easiest, unstructured codes usually require you to deal with both a local and global ordering, and there usually isn't any concept of a "natural ordering", just a concatenation of the owned local spaces (equivalent to "PETSc ordering" for structured grids). > > > >> cont=cont+1 > >> qloc_a(cont)=.... > > > > Can use qloc_a(i,j,k) here. > > I didn't manage to do this. I get compiling errors (I'm using ifort) if I > define a qloc_a(:,:,:) pointer array: PetscScalar, pointer :: qloc_a(:,:,:) Define the qloc_a(:) as usual, but index it with i,j,k (looks weird, I know), see src/snes/examples/tutorials/ex5f90.F. Jed
