[petsc-users] Insert values into matrix using MatSetValuesStencil or MatSetValuesLocal
TAY wee-beng
zonexo at gmail.com
Wed Aug 26 00:02:33 CDT 2015
Hi Timothee,
Yes, I only parallelized in j and k. ksta,jsta are the starting k and j
values. kend,jend are the ending k and j values.
However, now I am using only 1 procs.
I was going to resend you my code but then I realised my mistake. I used:
*/call
MatSetValuesStencil(A_mat,ione,row,iseven,col,value_insert,INSERT_VALUES,ierr)/**/
/*
for all pts, including those at the boundary. Hence, those coupling
outside the boundary is also included.
I changed to:
*/call
MatSetValuesStencil(A_mat,ione,row,ione,col(:,7),value_insert(7),INSERT_VALUES,ierr)/*
so I am now entering values individually.
Is there anyway I can use the 1st option to enter all the values
together even those some pts are invalid. I think it should be faster.
Can I somehow tell PETSc to ignore them?
Thank you
Yours sincerely,
TAY wee-beng
On 26/8/2015 12:24 PM, Timothée Nicolas wrote:
> What is the definition of ksta, kend, jsta, jend ? Etc ? You are
> parallelized only in j and k ?
>
> What I said about the "-1" holds only if you have translated the start
> and end points to FORTRAN numbering after getting the corners and
> ghost corners from the DMDA (see ex ex5f90.F from snes)
>
> Would you mind sending the complete routine with the complete
> definitions of ksta,kend,jsta,jend,and size_x ?
>
> Timothee
>
> 2015-08-26 13:12 GMT+09:00 TAY wee-beng <zonexo at gmail.com
> <mailto:zonexo at gmail.com>>:
>
> 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/4deb536d/attachment.html>
More information about the petsc-users
mailing list