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

[Timothée Nicolas]

Small erratum, in the declaration for v, it should be

PetscScalar :: v(n)

where n is the same as for col (6 in the example, not 8 which I copied from
my particular case)

2015-08-24 18:54 GMT+09:00 Timothée Nicolas <timothee.nicolas at gmail.com>:

> 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>:
>
>> 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/20150824/8cdf0d98/attachment-0001.html>