[petsc-users] unexpected ordering when VecSetValues set multiples values
Matteo Parsani
parsani.matteo at gmail.com
Fri Jun 7 14:57:16 CDT 2013
Hello Barry,
I have tried to pass to VectSetValues two values
*call
VecSetValues(x_vec_in,2,index_list(1:2),x_vec_local(1:2),INSERT_VALUES,ierr_local)
*
By looking at the values of index_list and x_vec_local with gdb I see the
correct numbers:
(gdb) print x_vec_local(1)
$1 = 0.9394877605834675
(gdb) print x_vec_local(2)
$2 = 0.55294328167143081
(gdb) print index_list(1)
$3 = 0
(gdb) print index_list(2)
$4 = 1
Now, when I ask gdb to print the values of ni, ix and y in
*PetscErrorCode VecSetValues_Seq(Vec xin,PetscInt ni,const PetscInt
ix[],const PetscScalar y[],InsertMode m)*
I get
(gdb) print ni
$5 = 2
(gdb) print ix[0]
$6 = 0
(gdb) print ix[1]
$7 = 0
print y[0]
$13 = 0.9394877605834675
(gdb) print y[1]
$14 = 0.55294328167143081
The values of y match with the values I am passing in, i.e. x_vec_local.
However the index ix is wrong because the first two components of it are
zero. Thus, when i = 1 in the VecSetValues_Seq for loop, the first element
that was placed correctly is replaced by the second element. This is also
waht I see if I use
call VecView(x_vec_in,PETSC_VIEWER_STDOUT_WORLD,ierr_local)
Now if I keep printing ix in gdb I see the following pattern:
(gdb) print ix[0]
$18 = 0
(gdb) print ix[1]
$19 = 0
(gdb) print ix[2]
$20 = 1
(gdb) print ix[3]
$21 = 0
(gdb) print ix[4]
$22 = 2
(gdb) print ix[5]
$23 = 0
(gdb) print ix[6]
$24 = 3
which I think is strange since I ix should be a vector of length 2 because
I am calling VecSetValues as
*call
VecSetValues(x_vec_in,2,index_list(1:2),x_vec_local(1:2),INSERT_VALUES,ierr_local)
*
Am I doing something wrong?
Thank you in advance.
On Fri, Jun 7, 2013 at 3:03 PM, Matteo Parsani <parsani.matteo at gmail.com>wrote:
> Yes, they are correct. When asked, the debugger prints the right values
> for bot the vector of the index and the vector of the values.
>
> Sorry, how can I set a break point in VecSetValues_Seq() and
> VecSetValues_MPI()?
> If I type
>
> break VecSetValues_Seq()
>
> or
>
> break VecSetValues_MPI()
>
> or
>
> /home0/pmatteo/research/lib_src/petsc/src/vec/vec/impls/seq/bvec2.c:1050
>
> I get
>
> Function "vecsetvalues_seq" not defined.
>
> or
>
> Function "vecsetvalues_mpi" not defined.
>
>
> Thanks
>
>
>
> On Fri, Jun 7, 2013 at 2:32 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
>
>>
>> Are, seemingly the correct values getting into the vecsetvalues_ call?
>>
>> Put a break point also in VecSetValues_Seq() and VecSetValues_MPI()
>> when it gets into either of those routines you can do step in the debugger
>> and watch it take the values passed in and attempt to put them in the
>> appropriate places in the vector's array.
>>
>> Barry
>>
>> On Jun 7, 2013, at 1:20 PM, Matteo Parsani <parsani.matteo at gmail.com>
>> wrote:
>>
>> > Hello Barry,
>> > I have tried to follow the debugger but I got lost.
>> >
>> > I have run my code in this way:
>> >
>> > ./NSE -start_in_debugger noxterm
>> >
>> > then after getting some output I get
>> >
>> > gdb()
>> >
>> > Then I type break vecsetvalues_
>> >
>> >
>> > and a break point is generated
>> >
>> > Breakpoint 1 at 0x2ae97eda5827: file
>> /scratch/home0/pmatteo/research/lib_src/petsc/src/vec/vec/interface/ftn-auto/rvectorf.c,
>> line 239.
>> >
>> >
>> > I keep typing n and it goes on. It shows me every time vectsetvalues is
>> called but it never goes inside vectsetvalues.
>> >
>> > I have tested the VecSetValues with just two elements and it does not
>> work. It sets correctly the values both in serial and in parallel if I set
>> one value at the time.
>> >
>> > Thanks
>> >
>> >
>> >
>> > On Tue, Apr 30, 2013 at 8:35 AM, Barry Smith <bsmith at mcs.anl.gov>
>> wrote:
>> >
>> > 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
>> >
>> >
>> >
>> >
>> > --
>> > Matteo
>>
>>
>
>
> --
> Matteo
>
--
Matteo
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20130607/45d21bce/attachment-0001.html>
More information about the petsc-users
mailing list