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

[Barry Smith]

  Shouldn't matter that it is called from Fortran we do it all the time.

  Does it work if the final m is not very large?

   You may need to run in the debugger and follow the values through the code to see why they don't get to where they belong.


   Barry

   Could also be a buggy fortran compiler.


On Apr 26, 2013, at 2:06 PM, Matteo Parsani <parsani.matteo at gmail.com> wrote:

> Hello Barry,
> sorry I modified few things just before to write the mail.
> The correct loop with the correct name of the variables is the following
> 
>         ! Number of DOFs owned by this process
>         len_local = (elem_high-elem_low+1)*nodesperelem*nequations
>         
>         ! Allocate memory for x_vec_local and index_list
>         allocate(x_vec_local(len_local))
>         allocate(index_list(len_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
>                     ! Add element to x_vec_local
>                     x_vec_local(m) = ug(ieq,inode,ielem)
>                     ! Add element to index list
>                     index_list(m) = (elem_low-1)*nodesperelem*nequations+m-1      
>                     ! Update m index
>                     m = m+1
>                 end do
>             end do
>         end do
> 
>         ! HERE I HAVE PRINTED x_vec_local, ug and index_list
> 
>         ! Set values in the portion of the vector owned by the process
>         call VecSetValues(x_vec_in,m-1,index_list,x_vec_local,INSERT_VALUES,& 
>                          & ierr_local)
> 
>         ! Assemble initial guess
>         call VecAssemblyBegin(x_vec_in,ierr_local)
>         call VecAssemblyEnd(x_vec_in,ierr_local)
> 
> 
> 
> I have printed the values and the indices I want to pass to VecSetValues() and they are correct.
> 
> I also printed the values after VecSetValues() has been called and they are wrong.
> 
> The attachment shows that.
> 
> 
> Could it be a problem of VecSetValues() + F90 when more than 1 elements is set?
> 
> I order to debug my the code I am running just with 1 processor. Thus the process owns all the DOFs.
> 
> Thank you.  
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> On Fri, Apr 26, 2013 at 2:39 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
> 
> 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>
> 
> 
> 
> 
> -- 
> Matteo
> <comparison.txt>