<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Tue, May 20, 2014 at 7:12 AM, Christophe Ortiz <span dir="ltr"><<a href="mailto:christophe.ortiz@ciemat.es" target="_blank">christophe.ortiz@ciemat.es</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 class="gmail_extra">I found another problem when using two-dimensional arrays defined using pointers of pointers.</div>
<div class="gmail_extra"><br></div><div class="gmail_extra">When I use a "classical" two-dimensional array defined by</div>
<div class="gmail_extra"><br></div><div class="gmail_extra">PetscScalar array[dof][dof];</div></div></blockquote><div><br></div><div>This declaration will use contiguous memory since its on the stack.</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 class="gmail_extra">and then build the Jacobian using </div><div class="gmail_extra"><br>
</div><div class="gmail_extra">ierr = MatSetValuesBlocked(*Jpre,1,&row,1,&col,&array[0][0],INSERT_VALUES);<br></div><div class="gmail_extra"><br></div><div class="gmail_extra">It works fine.</div><div class="gmail_extra">
<br></div><div class="gmail_extra">The problem comes when I define the two-dimensional array as follows:</div><div class="gmail_extra"><br></div><div class="gmail_extra">PetscScalar **array;</div></div></blockquote><div>
<br></div><div>This one uses non-contiguous memory since its on the heap.</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 class="gmail_extra">
<div class="gmail_extra"> array = (PetscScalar**)malloc(sizeof(PetscScalar*) * dof);</div><div class="gmail_extra"> for (k = 0; k < dof; k++){</div><div class="gmail_extra"> array[k] = (PetscScalar*)malloc(sizeof(PetscScalar) * dof);</div>
<div class="gmail_extra"> }</div><div class="gmail_extra"><br></div><div class="gmail_extra">When I pass it to MatSetValuesBlocked() there is a problem. Either Petsc complains because I am not passing it the right way or when it accepts it, wrong data is passed because the solution is not correct. Maybe Petsc expect dof*dof values and only dof are passed ?</div>
</div></div></blockquote><div><br></div><div>You can only pass contiguous memory to MatSetValues().</div><div><br></div><div> Matt</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 class="gmail_extra"><div class="gmail_extra">How a two-dimensional array declared with pointers of pointers should be passed to MatSetValuesBlocked() ?</div><div class="gmail_extra"><br></div><div class="gmail_extra">
Many thanks in advance.</div><div class="gmail_extra">Christophe</div><div class="gmail_extra"><br></div><div class="gmail_extra"> </div></div><div class="gmail_extra">
<br><br><div class="gmail_quote">On Thu, May 15, 2014 at 2:08 AM, Jed Brown <span dir="ltr"><<a href="mailto:jed@jedbrown.org" target="_blank">jed@jedbrown.org</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
<div>Christophe Ortiz <<a href="mailto:christophe.ortiz@ciemat.es" target="_blank">christophe.ortiz@ciemat.es</a>> writes:<br>
<br>
> Hi all,<br>
><br>
> I am experiencing some problems of memory corruption with PetscMemzero().<br>
><br>
> I set the values of the Jacobian by blocks using MatSetValuesBlocked(). To<br>
> do so, I use some temporary two-dimensional arrays[dof][dof] that I must<br>
> reset at each loop.<br>
><br>
> Inside FormIJacobian, for instance, I declare the following two-dimensional<br>
> array:<br>
><br>
> PetscScalar diag[dof][dof];<br>
><br>
> and then, to zero the array diag[][] I do<br>
><br>
> ierr = PetscMemzero(diag,dof*dof*sizeof(PetscScalar));<br>
<br>
</div>Note that this can also be spelled<br>
<br>
PetscMemzero(diag,sizeof diag);<br>
<div><br>
> Then, inside main(), once dof is determined, I allocate memory for diag as<br>
> follows:<br>
><br>
> diag = (PetscScalar**)malloc(sizeof(PetscScalar*) * dof);<br>
><br>
> for (k = 0; k < dof; k++){<br>
> diag[k] = (PetscScalar*)malloc(sizeof(PetscScalar) * dof);<br>
> }<br>
> That is, the classical way to allocate memory using the pointer notation.<br>
<br>
</div>Note that you can do a contiguous allocation by creating a Vec, then use<br>
VecGetArray2D to get 2D indexing of it.<br>
</blockquote></div><br></div></div>
</blockquote></div><br><br clear="all"><div><br></div>-- <br>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
-- Norbert Wiener
</div></div>