[petsc-users] DA matrices and vectors ordering

[Jed Brown]

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