[petsc-users] unexpected ordering when VecSetValues set multiples values

[Barry Smith]

On Apr 26, 2013, at 10:20 AM, Matteo Parsani <parsani.matteo at gmail.com> wrote:

> Hello,
> I have some problem when I try to set multiple values to a PETSc vector that I will use later on with SNES. I am using Fortran 90.
> Here the problem and two fixes that however are not so good for performances reasons. The code is very simple.
> 
> Standard approach that does not work correctly: (I am probably doing something wrong)
> 
>         m = 1
>         ! Loop over all elements
>         do ielem = elem_low, elem_high
>             ! Loop over all nodes in the element
>             do inode = 1, nodesperelem
>                 !Loop over all equations
>                 do ieq = 1, nequations
>                     ! Add element to x_vec_local
>                     x_vec_local(m) = ug(ieq,inode,ielem)
>                     ! Add element to index list
>                     ind(m) = (elem_low-1)*nodesperelem*nequations+m-1
>                     ! Update m index
>                     m = m+1
>                 end do
>             end do
>         end do
> 
>        ! Set values in the portion of the vector owned by the process
>         call VecSetValues(x_vec_in,len_local,index_list,x_vec_local,INSERT_VALUES,& 

    What is len_local and index_list? They do not appear in the loop above. Shouldn't you be passing m-1 for the length and ind for the indices?

    I would first print out the all the values in your input to VecSetValues() and make sure they are correct.

    Barry

>                                    & ierr_local)
> 
>         ! Assemble initial guess
>         call VecAssemblyBegin(x_vec_in,ierr_local)
>         call VecAssemblyEnd(x_vec_in,ierr_local)
> 
> Then I print my expected values and the values contained in the PETSc vector to a file. See attachment. I am running in serial for the moment BUT strangely if you look at the file I have attached the first 79 DOFs values have a wrong ordering and the remaining 80 are zero.
> 
> 
> 1st approach: set just one value at the time inside the loop.
>         m = 1
>         ! Loop over all elements
>         do ielem = elem_low, elem_high
>             ! Loop over all nodes in the element
>             do inode = 1, nodesperelem
>                 !Loop over all equations
>                 do ieq = 1, nequations
>                     ! Add element to x_vec_local
>                     value = ug(ieq,inode,ielem)
>                     ! Add element to index list
>                     ind = (elem_low-1)*nodesperelem*nequations+m-1
>                     call VecSetValues(x_vec_in,1,ind,value,INSERT_VALUES,& 
>                                        & ierr_local)
>                     ! Update m index
>                     m = m+1
>                 end do
>             end do
>         end do
>  
> 
> This works fine. As you can see I am using the same expression used in the previous loop to compute the index of the element that I have to add in the x_vec_in, i.e.
> ind = (elem_low-1)*nodesperelem*nequations+m-1
> 
> Thus I cannot see which is the problem. 
> 
> 2nd approach: get the pointer to the local part of the global vector and use it to set the values in the global vector
> 
>         m = 1
>         ! Loop over all elements
>         do ielem = elem_low, elem_high
>             ! Loop over all nodes in the element
>             do inode = 1, nodesperelem
>                 !Loop over all equations
>                 do ieq = 1, nequations
>                     ! Add element to x_vec_local
>                     tmp(m) = ug(ieq,inode,ielem)
>                     ! Update m index
>                     m = m+1
>                 end do
>             end do
>         end do
> 
> 
> This works fine too.
> 
> 
> Jut to be complete. I use the following two approaches to view the vector:
> 
> call VecView(x_vec_in,PETSC_VIEWER_STDOUT_WORLD,ierr_local)
> 
> 
> and 
> 
> call VecGetArrayF90(x_vec_in,tmp,ierr_local)
> 
> 
>         m = 1
>         ! Loop over all elements
>         do ielem = elem_low, elem_high
>             ! Loop over all nodes in the element
>             do inode = 1, nodesperelem
>                 !Loop over all equations
>                 do ieq = 1, nequations
>                     write(*,*) m,index_list(m),x_vec_local(m),tmp(m)                
>                     ! Update m index
>                     m = m+1
>                 end do
>             end do
>         end do
> 
> 
> Thank you.
> 
> 
> -- 
> Matteo
> <values.txt>