<div dir="ltr"><div><div>Hello Barry,<br></div>it does not work also if m is not so large.<br><br></div>Thanks again.<br><div><br> <br><br><br></div></div><div class="gmail_extra"><br><br><div class="gmail_quote">On Fri, Apr 26, 2013 at 3:23 PM, Barry Smith <span dir="ltr"><<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><br>
Shouldn't matter that it is called from Fortran we do it all the time.<br>
<br>
Does it work if the final m is not very large?<br>
<br>
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.<br>
<br>
<br>
Barry<br>
<br>
Could also be a buggy fortran compiler.<br>
<div><div class="h5"><br>
<br>
On Apr 26, 2013, at 2:06 PM, Matteo Parsani <<a href="mailto:parsani.matteo@gmail.com">parsani.matteo@gmail.com</a>> wrote:<br>
<br>
> Hello Barry,<br>
> sorry I modified few things just before to write the mail.<br>
> The correct loop with the correct name of the variables is the following<br>
><br>
> ! Number of DOFs owned by this process<br>
> len_local = (elem_high-elem_low+1)*nodesperelem*nequations<br>
><br>
> ! Allocate memory for x_vec_local and index_list<br>
> allocate(x_vec_local(len_local))<br>
> allocate(index_list(len_local))<br>
><br>
><br>
> m = 1<br>
> ! Loop over all elements<br>
> do ielem = elem_low, elem_high<br>
> ! Loop over all nodes in the element<br>
> do inode = 1, nodesperelem<br>
> !Loop over all equations<br>
> do ieq = 1, nequations<br>
> ! Add element to x_vec_local<br>
> x_vec_local(m) = ug(ieq,inode,ielem)<br>
> ! Add element to index list<br>
> index_list(m) = (elem_low-1)*nodesperelem*nequations+m-1<br>
> ! Update m index<br>
> m = m+1<br>
> end do<br>
> end do<br>
> end do<br>
><br>
> ! HERE I HAVE PRINTED x_vec_local, ug and index_list<br>
><br>
> ! Set values in the portion of the vector owned by the process<br>
> call VecSetValues(x_vec_in,m-1,index_list,x_vec_local,INSERT_VALUES,&<br>
> & ierr_local)<br>
><br>
> ! Assemble initial guess<br>
> call VecAssemblyBegin(x_vec_in,ierr_local)<br>
> call VecAssemblyEnd(x_vec_in,ierr_local)<br>
><br>
><br>
><br>
> I have printed the values and the indices I want to pass to VecSetValues() and they are correct.<br>
><br>
> I also printed the values after VecSetValues() has been called and they are wrong.<br>
><br>
> The attachment shows that.<br>
><br>
><br>
> Could it be a problem of VecSetValues() + F90 when more than 1 elements is set?<br>
><br>
> I order to debug my the code I am running just with 1 processor. Thus the process owns all the DOFs.<br>
><br>
> Thank you.<br>
><br>
><br>
><br>
><br>
><br>
><br>
><br>
><br>
><br>
><br>
><br>
> On Fri, Apr 26, 2013 at 2:39 PM, Barry Smith <<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>> wrote:<br>
><br>
> On Apr 26, 2013, at 10:20 AM, Matteo Parsani <<a href="mailto:parsani.matteo@gmail.com">parsani.matteo@gmail.com</a>> wrote:<br>
><br>
> > Hello,<br>
> > 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.<br>
> > Here the problem and two fixes that however are not so good for performances reasons. The code is very simple.<br>
> ><br>
> > Standard approach that does not work correctly: (I am probably doing something wrong)<br>
> ><br>
> > m = 1<br>
> > ! Loop over all elements<br>
> > do ielem = elem_low, elem_high<br>
> > ! Loop over all nodes in the element<br>
> > do inode = 1, nodesperelem<br>
> > !Loop over all equations<br>
> > do ieq = 1, nequations<br>
> > ! Add element to x_vec_local<br>
> > x_vec_local(m) = ug(ieq,inode,ielem)<br>
> > ! Add element to index list<br>
> > ind(m) = (elem_low-1)*nodesperelem*nequations+m-1<br>
> > ! Update m index<br>
> > m = m+1<br>
> > end do<br>
> > end do<br>
> > end do<br>
> ><br>
> > ! Set values in the portion of the vector owned by the process<br>
> > call VecSetValues(x_vec_in,len_local,index_list,x_vec_local,INSERT_VALUES,&<br>
><br>
> 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?<br>
><br>
> I would first print out the all the values in your input to VecSetValues() and make sure they are correct.<br>
><br>
> Barry<br>
><br>
> > & ierr_local)<br>
> ><br>
> > ! Assemble initial guess<br>
> > call VecAssemblyBegin(x_vec_in,ierr_local)<br>
> > call VecAssemblyEnd(x_vec_in,ierr_local)<br>
> ><br>
> > 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.<br>
> ><br>
> ><br>
> > 1st approach: set just one value at the time inside the loop.<br>
> > m = 1<br>
> > ! Loop over all elements<br>
> > do ielem = elem_low, elem_high<br>
> > ! Loop over all nodes in the element<br>
> > do inode = 1, nodesperelem<br>
> > !Loop over all equations<br>
> > do ieq = 1, nequations<br>
> > ! Add element to x_vec_local<br>
> > value = ug(ieq,inode,ielem)<br>
> > ! Add element to index list<br>
> > ind = (elem_low-1)*nodesperelem*nequations+m-1<br>
> > call VecSetValues(x_vec_in,1,ind,value,INSERT_VALUES,&<br>
> > & ierr_local)<br>
> > ! Update m index<br>
> > m = m+1<br>
> > end do<br>
> > end do<br>
> > end do<br>
> ><br>
> ><br>
> > 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.<br>
> > ind = (elem_low-1)*nodesperelem*nequations+m-1<br>
> ><br>
> > Thus I cannot see which is the problem.<br>
> ><br>
> > 2nd approach: get the pointer to the local part of the global vector and use it to set the values in the global vector<br>
> ><br>
> > m = 1<br>
> > ! Loop over all elements<br>
> > do ielem = elem_low, elem_high<br>
> > ! Loop over all nodes in the element<br>
> > do inode = 1, nodesperelem<br>
> > !Loop over all equations<br>
> > do ieq = 1, nequations<br>
> > ! Add element to x_vec_local<br>
> > tmp(m) = ug(ieq,inode,ielem)<br>
> > ! Update m index<br>
> > m = m+1<br>
> > end do<br>
> > end do<br>
> > end do<br>
> ><br>
> ><br>
> > This works fine too.<br>
> ><br>
> ><br>
> > Jut to be complete. I use the following two approaches to view the vector:<br>
> ><br>
> > call VecView(x_vec_in,PETSC_VIEWER_STDOUT_WORLD,ierr_local)<br>
> ><br>
> ><br>
> > and<br>
> ><br>
> > call VecGetArrayF90(x_vec_in,tmp,ierr_local)<br>
> ><br>
> ><br>
> > m = 1<br>
> > ! Loop over all elements<br>
> > do ielem = elem_low, elem_high<br>
> > ! Loop over all nodes in the element<br>
> > do inode = 1, nodesperelem<br>
> > !Loop over all equations<br>
> > do ieq = 1, nequations<br>
> > write(*,*) m,index_list(m),x_vec_local(m),tmp(m)<br>
> > ! Update m index<br>
> > m = m+1<br>
> > end do<br>
> > end do<br>
> > end do<br>
> ><br>
> ><br>
> > Thank you.<br>
> ><br>
> ><br>
> > --<br>
> > Matteo<br>
> > <values.txt><br>
><br>
><br>
><br>
><br>
> --<br>
> Matteo<br>
</div></div>> <comparison.txt><br>
<br>
</blockquote></div><br><br clear="all"><br>-- <br><div dir="ltr">Matteo<br></div>
</div>