[petsc-users] Insert values into matrix using MatSetValuesStencil or MatSetValuesLocal
Timothée Nicolas
timothee.nicolas at gmail.com
Wed Aug 26 00:15:10 CDT 2015
I don't really understand what you say, but it does not sound right.
You can enter the boundary points separately and then the points outside
the boundary on separate calls, like this :
do j=user%ys,user%ye
do i=user%xs,user%xe
if (i.eq.1 .or. i.eq.user%mx .or. j .eq. 1 .or. j .eq. user%my) then
! boundary point
row(MatStencil_i,1) = i -1
row(MatStencil_j,1) = j -1
row(MatStencil_c,1) = 1 -1
col(MatStencil_i,1) = i -1
col(MatStencil_j,1) = j -1
col(MatStencil_c,1) = 1 -1
v(1) = one
call MatSetValuesStencil(jac_prec,ione,row,ione,col,v, &
& INSERT_VALUES,ierr)
else
row(MatStencil_i,1) = i -1
row(MatStencil_j,1) = j -1
row(MatStencil_c,1) = 1 -1
col(MatStencil_i,1) = i -1
col(MatStencil_j,1) = j -1
col(MatStencil_c,1) = 1 -1
v(1) = undemi*dxm1*(vx_ip1j-vx_im1j) +
two*user%nu*(dxm1**2+dym1**2)
col(MatStencil_i,2) = i+1 -1
col(MatStencil_j,2) = j -1
col(MatStencil_c,2) = 1 -1
v(2) = undemi*dxm1*(vx_ij-vx_ip1j) - user%nu*dxm1**2
col(MatStencil_i,3) = i-1 -1
col(MatStencil_j,3) = j -1
col(MatStencil_c,3) = 1 -1
v(3) = -undemi*dxm1*(vx_ij-vx_im1j) - user%nu*dxm1**2
col(MatStencil_i,4) = i -1
col(MatStencil_j,4) = j+1 -1
col(MatStencil_c,4) = 1 -1
v(4) = undemi*dym1*vy_ij - user%nu*dym1**2
col(MatStencil_i,5) = i -1
col(MatStencil_j,5) = j-1 -1
col(MatStencil_c,5) = 1 -1
v(5) = -undemi*dym1*vy_ij - user%nu*dym1**2
col(MatStencil_i,6) = i -1
col(MatStencil_j,6) = j -1
col(MatStencil_c,6) = 2 -1
v(6) = undemi*dym1*(vx_ijp1-vx_ijm1)
col(MatStencil_i,7) = i+1 -1
col(MatStencil_j,7) = j -1
col(MatStencil_c,7) = 2 -1
v(7) = -undemi*dxm1*vy_ip1j
col(MatStencil_i,8) = i-1 -1
col(MatStencil_j,8) = j -1
col(MatStencil_c,8) = 2 -1
v(8) = undemi*dxm1*vy_im1j
call MatSetValuesStencil(jac_prec,ione,row,ieight,col,v, &
& INSERT_VALUES,ierr)
endif
enddo
enddo
Timothee
2015-08-26 14:02 GMT+09:00 TAY wee-beng <zonexo at gmail.com>:
> 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>:
>
>> 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>
>> 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/0301a30e/attachment-0001.html>
More information about the petsc-users
mailing list