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

TAY wee-beng zonexo at gmail.com
Tue Aug 25 23:12:49 CDT 2015

```Hi,

I have wrote the routine for my Poisson eqn. I have only 1 DOF, which is
for pressure. The center cell is coupled with 6 other cells (north,
south, east, west, front, back), so together 7 couplings.

size x/y/z = 4/8/10

*/MatStencil  :: row(4,1),col(4,7)/**/
/**/
/**/PetscScalar :: value_insert(7)/**/
/**/
/**/PetscInt :: ione,iseven/**/
/**/
/**/ione = 1;   iseven = 7/**/
/**/
/**/do k=ksta,kend/**/
/**/
/**/        do j = jsta,jend/**/
/**/
/**/            do i=1,size_x/**/
/**//**/
/**/                row(MatStencil_i,1) = i - 1/**/
/**//**/
/**/                row(MatStencil_j,1) = j - 1/**/
/**//**/
/**/                row(MatStencil_k,1) = k - 1/**/
/**//**/
/**/                row(MatStencil_c,1) = 0 ! 1 - 1/**/
/**//**/
/**/                value_insert = 0.d0/**/
/**//**/
/**/                if (i /= size_x) then/**/
/**//**/
/**/                    col(MatStencil_i,3) = i + 1 - 1 !east/**/
/**//**/
/**/                    col(MatStencil_j,3) = j - 1/**/
/**//**/
/**/                    col(MatStencil_k,3) = k - 1/**/
/**//**/
/**/                    col(MatStencil_c,3) = 0/**/
/**//**/
/**/                    value_insert(3) =
(cp_yz(j,k)%fc_E)/(cp_x(i)%pd_E+cp_x(i+1)%pd_W)/**/
/**//**/
/**/                end if/**/
/**//**/
/**/                if (i /= 1) then/**/
/**//**/
/**/                    col(MatStencil_i,5) = i - 1 - 1 !west/**/
/**//**/
/**/                    col(MatStencil_j,5) = j - 1/**/
/**//**/
/**/                    col(MatStencil_k,5) = k - 1/**/
/**//**/
/**/                    col(MatStencil_c,5) = 0/**/
/**//**/
/**/                    value_insert(5) =
(cp_yz(j,k)%fc_E)/(cp_x(i)%pd_W+cp_x(i-1)%pd_E)/**/
/**//**/
/**/                end if/**/
/**//**/
/**/                if (j /= size_y) then/**/
/**//**/
/**/                    col(MatStencil_i,2) = i - 1 !north/**/
/**//**/
/**/                    col(MatStencil_j,2) = j + 1 - 1/**/
/**//**/
/**/                    col(MatStencil_k,2) = k - 1/**/
/**//**/
/**/                    col(MatStencil_c,2) = 0/**/
/**//**/
/**/                    value_insert(2) =
(cp_zx(i,k)%fc_N)/(cp_y(j)%pd_N+cp_y(j+1)%pd_S)/**/
/**//**/
/**/                end if/**/
/**//**/
/**/                .../**/
/**//**/
/**/                col(MatStencil_i,1) = i - 1/**/
/**//**/
/**/                col(MatStencil_j,1) = j - 1/**/
/**//**/
/**/                col(MatStencil_k,1) = k - 1/**/
/**//**/
/**/                col(MatStencil_c,1) = 0/**/
/**//**/
/**/                value_insert(1) = -value_insert(2) - value_insert(3)
- value_insert(4) - value_insert(5) - value_insert(6) - value_insert(7)/**/
/**//**/
/**/                call
MatSetValuesStencil(A_mat,ione,row,iseven,col,value_insert,INSERT_VALUES,ierr)/**/
/**//**/
/**/            end do/**/
/**//**/
/**/        end do/**/
/**//**/
/**/    end do/*

but I got the error :

[0]PETSC ERROR: Argument out of range
[0]PETSC ERROR: Inserting a new nonzero at (3,0) in the matrix.

The error happens at i = 4, j = 1, k = 1. So I guess it has something to
do with the boundary condition. However, I can't figure out what's
wrong. Can someone help?

Thank you

Yours sincerely,

TAY wee-beng

On 24/8/2015 5:54 PM, Timothée Nicolas 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
>
>
> 2015-08-24 18:09 GMT+09:00 Wee-Beng Tay <zonexo at gmail.com
> <mailto: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/20150826/e41bc9ac/attachment.html>
```

More information about the petsc-users mailing list