[petsc-users] unexpected ordering when VecSetValues set multiples values
Satish Balay
balay at mcs.anl.gov
Fri Jun 7 15:14:22 CDT 2013
On Fri, 7 Jun 2013, Barry Smith wrote:
>
> How large is your integer on the Fortran side and the C side?
>
> Normally a PetscInt is a C int and a Fortran integer*4 which are each 4 bytes long.
>
> If PETSc is configured with --with-64-bit-indices then PetscInt is a C long long int and Fortran integer*8 both 8 bytes long.
>
> It looks like you are using 8 byte Fortran integers but have not compiled PETSc with --with-64-bit-indices Possibly you are using come fortran compile time flag that "promotes" float to double and integer*4 to integer*8 in Fortran.
>
> 1) we recommend never using any fortran compiler flag that promotes types automatically.
If Petsc datatypes are used in fortran code [i.e declare variable as
PetscInt - and not integer for functions that expect PetscInt] - then
the code [from PETSc part] would be resiliant to -i8 type options.
Satish
>
> 2) likely you just want to use 4 byte integers in fortran. If you want to use 8 bytes then configure PETSc with --with-64-bit-indices
>
> Barry
>
>
>
>
> On Jun 7, 2013, at 2:57 PM, Matteo Parsani <parsani.matteo at gmail.com> wrote:
>
> > 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
>
>
More information about the petsc-users
mailing list