[petsc-users] unexpected ordering when VecSetValues set multiples values
Barry Smith
bsmith at mcs.anl.gov
Tue Apr 30 07:35:14 CDT 2013
On Apr 30, 2013, at 7:21 AM, Matteo Parsani <parsani.matteo at gmail.com> wrote:
> Hello Barry,
> it does not work also if m is not so large.
Then you need to run in the debugger. You can use the command line options
-start_in_debugger noxterm
in the debugger put a break point in vecsetvalues_ (yes lower case followed by an underscore) then when it
gets to that point print the size of the array and the values in the integer and double array to see if they match, then step into the VecSetValues() function that is called by vecsetvalues_() and look at the array values ago, continue to step along and it will go into VecSetValues_Seq() and start actually putting the values into the Vec array.
Barry
>
> 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
More information about the petsc-users
mailing list