[petsc-users] Insert values into matrix using MatSetValuesStencil or MatSetValuesLocal

[Wee Beng Tay]

Sent using CloudMagic [https://cloudmagic.com/k/d/mailapp?ct=pa&cv=7.2.9&pv=5.0.2] On Mon, Aug 24, 2015 at 5:54 PM, Timothée Nicolas < timothee.nicolas at gmail.com [timothee.nicolas at gmail.com] > wrote:
Hi,

ex5 of snes can give you an example of the two routines.

The C version ex5.c uses MatSetValuesStencil whereas the Fortran90 version
ex5f90.F uses MatSetValuesLocal.

However, I use MatSetValuesStencil also in Fortran, there is no problem, and no
need to mess around with DMDAGetAO, I think.

To input values in the matrix, you need to do the following :

! Declare the matstencils for matrix columns and rows
MatStencil :: row(4,1),col(4,n)
! Declare the quantity which will store the actual matrix elements
PetscScalar :: v(8)

The first dimension in row and col is 4 to allow for 3 spatial dimensions (even
if you use only 2) plus one degree of freedom if you have several fields in your
DMDA. The second dimension is 1 for row (you input one row at a time) and n for
col, where n is the number of columns that you input. For instance, if at node
(1,i,j) (1 is the index of the degree of freedom), you have, say, 6 couplings,
with nodes (1,i,j), (1,i+1,j), (1,i-1,j), (1,i,j-1), (1,i,j+1), (2,i,j) for
example, then you need to set n=6

Then you define the row number by naturally doing the following, inside a local
loop :

row(MatStencil_i,1) = i -1
row(MatStencil_j,1) = j -1
row(MatStencil_c,1) = 1 -1

the -1 are here because FORTRAN indexing is different from the native C
indexing. I put them on the right to make this more apparent.

Then the column information. For instance to declare the coupling with node
(1,i,j), (1,i-1,j) and (2,i,j) (you can make up for the rest) you will have to
write (still within the same local loop on i and j)

col(MatStencil_i,1) = i -1
col(MatStencil_j,1) = j -1
col(MatStencil_c,1) = 1 -1
v(1) = whatever_it_is

col(MatStencil_i,2) = i-1 -1
col(MatStencil_j,2) = j -1
col(MatStencil_c,2) = 1 -1
v(2) = whatever_it_is

col(MatStencil_i,3) = i -1
col(MatStencil_j,3) = j -1
col(MatStencil_c,3) = 2 -1
v(3) = whatever_it_is

...
...
..

...
...
...

Note that the index of the degree of freedom (or what field you are coupling
to), is indicated by MatStencil_c


Finally use MatSetValuesStencil

ione = 1
isix = 6
call MatSetValuesStencil(Matrix,ione,row,isix,col,v,INSERT_VALUES,ierr)

If it is not clear don't hesitate to ask more details. For me it worked that
way, I succesfully computed a Jacobian that way. It is very sensitive. If you
slightly depart from the right jacobian, you will see a huge difference compared
to using matrix free with -snes_mf, so you can hardly make a mistake because you
would see it. That's how I finally got it to work.

Best

Timothee


Hi Timothee,
Thanks for the help. So for boundary pts I will just leave blank for non
existent locations?
Also, can I use PETSc multigrid to solve this problem? This is a poisson eqn.
2015-08-24 18:09 GMT+09:00 Wee-Beng Tay < zonexo at gmail.com [zonexo at gmail.com] > :
Hi,
I'm modifying my 3d fortran code from MPI along 1 direction (z) to MPI along 2
directions (y,z)
Previously I was using MatSetValues with global indices. However, now I'm using
DM and global indices is much more difficult.
I come across MatSetValuesStencil or MatSetValuesLocal.
So what's the difference bet the one since they both seem to work locally?
Which is a simpler/better option?
Is there an example in Fortran for MatSetValuesStencil?
Do I also need to use DMDAGetAO together with MatSetValuesStencil or
MatSetValuesLocal?
Thanks!
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20150825/c4f34d87/attachment.html>