[petsc-users] unexpected ordering when VecSetValues set multiples values
Matteo Parsani
parsani.matteo at gmail.com
Tue Apr 30 07:21:31 CDT 2013
Hello Barry,
it does not work also if m is not so large.
Thanks again.
On Fri, Apr 26, 2013 at 3:23 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
>
> 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>
>
>
--
Matteo
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20130430/f797188c/attachment-0001.html>
More information about the petsc-users
mailing list