<div dir="ltr">On Sat, Dec 29, 2012 at 11:15 PM, amlan barua <span dir="ltr"><<a href="mailto:abarua@iit.edu" target="_blank">abarua@iit.edu</a>></span> wrote:<br><div class="gmail_extra"><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div dir="ltr"><div><div><div>Hi,<br></div>I am actually trying to implement a 'parallel' ordinary differential equation solver. <br>For proper functioning of the algorithm (the name is parareal), I need to implement a simple recurrence relation of the form x[i+1] = f(x[i]), f known, depends on quadrature one would like to use. <br>
</div>What is the best way to implement a sequential operation on a parallel structure?<br></div><div>Somehow I need to keep all but one process idle which.</div></div></blockquote><div><br></div><div style>It's not very parallel when all but one process is idle. ;-D</div>
<div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div> So I wrote a loop over all processes and within the loop I forced only one process to update its part. Say when j=0 ideally only the processor 0 should work. But others will update their local j value before 0th processor much faster. Thus I am concerned about the safety of the operation. Will it be okay if I modify my code as following, please advice? <br>
</div></div></blockquote><div><br></div><div style>Yes, this works, but the performance won't be great when using DMGlobalToLocal despite only really wanting to send the update to the next process in the sequence. (Maybe you don't care about performance yet, or this particular part won't be performance-sensitive.)</div>
<div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div>
PetscInitialize(&argc,&argv,(char *)0,help);<br><div dir="ltr"><div class="im"> ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size); CHKERRQ(ierr); <br> ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr); <br>
ierr = DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,size*5,1,1,PETSC_NULL,&da); CHKERRQ(ierr);<br> ierr = DMCreateGlobalVector(da,&vec); CHKERRQ(ierr);<br> ierr = VecSet(vec,1.00);<br> ierr = DMCreateLocalVector(da,&local);<br>
ierr = DMDAGetLocalInfo(da,&info);<br></div><div class="im"> temp = 1;<br> for (j=0;j<size;j++) {<br></div><div class="im"> ierr = DMGlobalToLocalBegin(da,vec,INSERT_VALUES,local); CHKERRQ(ierr);<br>
ierr = DMGlobalToLocalEnd(da,vec,INSERT_VALUES,local); CHKERRQ(ierr);<br></div><div class="im"> ierr = DMDAVecGetArray(da,vec,&arr);<br> ierr = DMDAVecGetArray(da,local,&array);<br></div><div class="im">
if (rank==j) {<br> for (i=info.xs;i<info.xs+info.xm;i++) {<br>
if ((!i)==0) {<br></div></div></div></div></blockquote><div><br></div><div style>I would write if (i) unless I was intentionally trying to confuse the reader.</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div dir="ltr"><div><div dir="ltr"><div class="im"> array[i] = array[i] + array[i-1];<br></div><div class="im">
arr[i] = array[i];<br> }<br> }<br> }<br> ierr = DMDAVecRestoreArray(da,local,&array);<br><div dir="ltr"> ierr = DMDAVecRestoreArray(da,vec,&arr); <br> }<br></div> ierr = VecView(vec,PETSC_VIEWER_STDOUT_WORLD);<br>
PetscFinalize();<br> return 0;</div></div></div><span class="HOEnZb"><font color="#888888">Amlan<br></font></span></div><div class="HOEnZb"><div class="h5"><div class="gmail_extra"><br><br><div class="gmail_quote">
On Sat, Dec 29, 2012 at 10:19 AM, Jed Brown <span dir="ltr"><<a href="mailto:jedbrown@mcs.anl.gov" target="_blank">jedbrown@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"><div dir="ltr"><div>On Sat, Dec 29, 2012 at 1:40 AM, amlan barua <span dir="ltr"><<a href="mailto:abarua@iit.edu" target="_blank">abarua@iit.edu</a>></span> wrote:<br>
</div><div class="gmail_extra"><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div dir="ltr"><div><div>Hi Barry,<br></div><div>I wrote the following piece according to your suggestions. Currently it does nothing but creates a vector with 1 at 1th position, 2 at 2th and so on. But I made it serial, i.e. (n+1)th place is computed using the value of nth place. My question, did I do it correctly, i.e. is it safe or results may change depending on problem size? This is much faster than VecSetValues, I believe the communication is minimum here because I take the advantage of ghost points.<br>
</div></div><div>Amlan<br><br>PetscInitialize(&argc,&argv,(char *)0,help);<br> ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size); CHKERRQ(ierr); <br> ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr); <br>
ierr = DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,size*5,1,1,PETSC_NULL,&da); CHKERRQ(ierr);<br> ierr = DMCreateGlobalVector(da,&vec); CHKERRQ(ierr);<br> ierr = VecSet(vec,1.00);<br> ierr = DMCreateLocalVector(da,&local);<br>
ierr = DMDAGetLocalInfo(da,&info);<br> ierr = DMDAVecGetArray(da,vec,&arr);<br> ierr = DMDAVecGetArray(da,local,&array);<br> temp = 1;<br> for (j=0;j<size;j++) {<br></div></div></blockquote>
<div>
<br></div><div>This is needlessly sequential (or should be).</div><div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"> ierr = DMGlobalToLocalBegin(da,vec,INSERT_VALUES,local); CHKERRQ(ierr);<br>
ierr = DMGlobalToLocalEnd(da,vec,INSERT_VALUES,local); CHKERRQ(ierr);<br></div></blockquote></div><div><br>You should never use a communication routine while you have access to the array (*VecGetArray()).<br> </div>
<div>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"> if (rank==j) {<br> for (i=info.xs;i<info.xs+info.xm;i++) {<br> if ((!i)==0) {<br>
array[i] = array[i] + array[i-1];<br></div></blockquote><div><br></div></div><div>What sort of recurrence do you actually want to implement. When possible, it's much better to reorganize so that you can do local work followed by an MPI_Scan followed by more local work. MPI_Scan is fast (logarithmic).</div>
<div><div>
<div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr">
arr[i] = array[i];<br> }<br> }<br> }<br> }<br> ierr = DMDAVecRestoreArray(da,local,&array);<br> ierr = DMDAVecRestoreArray(da,vec,&arr);<br> ierr = VecView(vec,PETSC_VIEWER_STDOUT_WORLD);<br>
PetscFinalize();<br> return 0;<br><br>
</div><div><div><div class="gmail_extra"><br><br><div class="gmail_quote">On Thu, Dec 27, 2012 at 10:40 AM, 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"><div><br>
On Dec 27, 2012, at 10:34 AM, amlan barua <<a href="mailto:abarua@iit.edu" target="_blank">abarua@iit.edu</a>> wrote:<br>
<br>
> I think I can use VecSetValues, is that right?<br>
<br>
</div> Yes you could do that. But since you are using a DMDA you could also use DMGetLocalVector(), DMGlobalToLocalBegin/End() followed by DMDAVecGetArray() to access the ghost values.<br>
<span><font color="#888888"><br>
Barry<br>
</font></span><div><div><br>
> Amlan<br>
><br>
><br>
> On Thu, Dec 27, 2012 at 9:04 AM, amlan barua <<a href="mailto:abarua@iit.edu" target="_blank">abarua@iit.edu</a>> wrote:<br>
> Hi Barry,<br>
> Is this scattering a very costly operation? I have to compute x[i] = f(x[i-1]) where f is known. Since this operation is strictly sequential, I thought of gathering the entire vector on processor 0, do the sequential operation there and scatter the result back. However this is unnecessary because I only need the bordering x[i] values. What can be a better way?<br>
> Amlan<br>
><br>
><br>
> On Thu, Dec 27, 2012 at 8:18 AM, Barry Smith <<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>> wrote:<br>
><br>
> ierr = DMDACreateNaturalVector(da,&natural);CHKERRQ(ierr);<br>
> ierr = DMDAGlobalToNaturalBegin(da,xin,INSERT_VALUES,natural);CHKERRQ(ierr);<br>
> ierr = DMDAGlobalToNaturalEnd(da,xin,INSERT_VALUES,natural);CHKERRQ(ierr);<br>
><br>
> Now do VecScatterCreateToZero() from natural and the vector will be in the natural ordering on process zero with the dof interlaced.<br>
><br>
><br>
> Barry<br>
><br>
> On Dec 27, 2012, at 12:22 AM, amlan barua <<a href="mailto:abarua@iit.edu" target="_blank">abarua@iit.edu</a>> wrote:<br>
><br>
> > Hi,<br>
> > Is there an analogue of VecScatterCreateToZero for DA vectors? The DMDA object has more than one degrees of freedom.<br>
> > If there isn't any, should I use an IS object to do the scattering?<br>
> > Amlan<br>
><br>
><br>
><br>
<br>
</div></div></blockquote></div><br></div>
</div></div></blockquote></div></div></div><br></div></div>
</blockquote></div><br></div>
</div></div></blockquote></div><br></div></div>