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>
