<div dir="ltr"><div><div><div><div><div>Hi,<br><br></div>ex5 of snes can give you an example of the two routines.<br><br></div>The C version ex5.c uses MatSetValuesStencil whereas the Fortran90 version ex5f90.F uses MatSetValuesLocal.<br><br></div>However, I use MatSetValuesStencil also in Fortran, there is no problem, and no need to mess around with DMDAGetAO, I think.<br><br></div>To input values in the matrix, you need to do the following :<br><br></div>! Declare the matstencils for matrix columns and rows<br><div>MatStencil :: row(4,1),col(4,n)<br></div><div>! Declare the quantity which will store the actual matrix elements<br></div><div>PetscScalar :: v(8)<br><br></div><div>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<br><br></div><div>Then you define the row number by naturally doing the following, inside a local loop :<br><br>row(MatStencil_i,1) = i -1<br>row(MatStencil_j,1) = j -1<br>row(MatStencil_c,1) = 1 -1<br><br></div><div>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.<br><br></div><div>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)<br><br>col(MatStencil_i,1) = i -1<br>col(MatStencil_j,1) = j -1<br>col(MatStencil_c,1) = 1 -1<br>v(1) = whatever_it_is<br><br>col(MatStencil_i,2) = i-1 -1<br>col(MatStencil_j,2) = j -1<br>col(MatStencil_c,2) = 1 -1<br>v(2) = whatever_it_is<br><br>col(MatStencil_i,3) = i -1<br>col(MatStencil_j,3) = j -1<br>col(MatStencil_c,3) = 2 -1<br>v(3) = whatever_it_is<br><br>...<br>...<br>..<br><br>...<br>...<br>...<br><br></div><div>Note that the index of the degree of freedom (or what field you are coupling to), is indicated by MatStencil_c<br><br></div><div><br></div><div>Finally use MatSetValuesStencil<br></div><div><br></div><div>ione = 1<br></div><div>isix = 6<br></div><div>call MatSetValuesStencil(Matrix,ione,row,isix,col,v,INSERT_VALUES,ierr)<br><br></div><div>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.<br><br></div><div>Best<br><br></div><div>Timothee<br><br></div></div><div class="gmail_extra"><br><div class="gmail_quote">2015-08-24 18:09 GMT+09:00 Wee-Beng Tay <span dir="ltr"><<a href="mailto:zonexo@gmail.com" target="_blank">zonexo@gmail.com</a>></span>:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr">Hi,<div><br></div><div>I'm modifying my 3d fortran code from MPI along 1 direction (z) to MPI along 2 directions (y,z)</div><div><br></div><div>Previously I was using MatSetValues with global indices. However, now I'm using DM and global indices is much more difficult.</div><div><br></div><div>I come across MatSetValuesStencil or MatSetValuesLocal.</div><div><br></div><div>So what's the difference bet the one since they both seem to work locally?</div><div><br></div><div>Which is a simpler/better option?</div><div><br></div><div>Is there an example in Fortran for MatSetValuesStencil?</div><div><br></div><div>Do I also need to use DMDAGetAO together with MatSetValuesStencil or MatSetValuesLocal?</div><div><br></div><div>Thanks!</div></div>
</blockquote></div><br></div>