[petsc-users] use of VecGetArray() with dof>1 (block size > 1)

[Ed Bueler]

Dave --

Perhaps you should re-read my questions.  My points were missed.  But you
motivated me to look slightly deeper; see below.

> I would say the reason why the last arg to VecGetArray() is not void* is
because
> it is intended to give you direct access to the pointer associated with
the entries
> within the vector - these are also PetscScalar's

Yes, that is fine.  But DMDAVecGetArray() does it the other way, using
"void *" instead.  The question was, why does the PETSc API have two
different approaches.

Jed's "more information" answer about unstructured grids is a bit vague,
but it must be the one that motivates the difference.

> I don't believe it is ever recommended anywhere to do the following:
>  VecGetArray(DM da,Vec local,(PetscScalar**)&u)

I do not recommend it; I am annoyed that I need it!

My point was that this awkward construct allowed the "recommended
approach"---which is a quote from the PETSc manual and not my own
recommendation---to proceed with VecGetArray().  Whereas this awkward
construct is not needed in user code which uses DMDAVecGetArray().

> Trying to trick the compile with such a cast is just begging for memory
> corruption to occur.

Maybe so, but you do it all the time with PETSc structured grids.  In fact,
the third line of the implementation of DMDAVecGetArray() is

" VecGetArray(vec,(PetscScalar**)array); "

That is, in the "recommended approach", the cast occurs inside
DMDAVecGetArray().  The same cast is simply exposed to the user if they
want to get a struct-valued pointer from VecGetArray().

Ed


On Fri, May 27, 2016 at 12:06 PM, Dave May <dave.mayhem23 at gmail.com> wrote:

>
>
> On 27 May 2016 at 20:34, Ed Bueler <elbueler at alaska.edu> wrote:
>
>> Dear PETSc --
>>
>> This is an "am I using it correctly" question.  Probably the API has the
>> current design because of something I am missing.
>>
>> First, a quote from the PETSc manual which I fully understand; it works
>> great and gives literate code (to the extent possible...):
>>
>> """
>> The recommended approach for multi-component PDEs is to declare a struct
>> representing the fields defined
>> at each node of the grid, e.g.
>>
>> typedef struct {
>> PetscScalar u,v,omega,temperature;
>> } Node;
>>
>> and write residual evaluation using
>>
>> Node **f,**u;
>> DMDAVecGetArray(DM da,Vec local,&u);
>> DMDAVecGetArray(DM da,Vec global,&f);
>> ...
>> f[i][j].omega = ...
>>
>
> Note that here the indexing should be
>   f[ *j* ][ *i* ].omega
>
>
>
>> ...
>> DMDAVecRestoreArray(DM da,Vec local,&u);
>> DMDAVecRestoreArray(DM da,Vec global,&f);
>> """
>>
>> Now the three questions:
>>
>> 1)  The third argument to DMDAVec{Get,Restore}Array() is of type "void
>> *".  It makes the above convenient.  But the third argument of the
>> unstructured version Vec{Get,Restore}Array() is of type "PetscScalar **",
>> which means that in an unstructured case, with the same Node struct, I
>> would write
>> "VecGetArray(DM da,Vec local,(PetscScalar **)&u);"
>> to get the same functionality.  Why is it this way?  More specifically,
>> why not have the argument to VecGetArray() be of type "void *"?
>>
>
> I would say the reason why the last arg to VecGetArray() is not void* is
> because it is intended to give you direct access to the pointer associated
> with the entries within the vector - these are also PetscScalar's
>
>
>>
>> 2) Given that the "recommended approach"
>>
>
> I don't believe it is ever recommended anywhere to do the following:
>   VecGetArray(DM da,Vec local,(PetscScalar**)&u)
>
> Trying to trick the compile with such a cast is just begging for memory
> corruption to occur.
>
> above works just fine, why do DMDAVec{Get,Restore}ArrayDOF() exist?  (I.e.
>> is there something I am missing about C indexing?)
>>
>
> As an additional point, DMDAVec{Get,Restore}ArrayDOF() return
>
> void *array
>
> so that the same API will work for 1D, 2D and 3D DMDA's which would
> require PetscScalar **data, PetscScalar ***data, PetscScalar ****data
> respectively.
>
>
> Cheers,
>   Dave
>
>
>> 3) There are parts of the PETSc API that refer to "dof" and parts that
>> refer to "block size".  Is this a systematic distinction with an underlying
>> reason?  It seems "block size" is more generic, but also it seems that it
>> could replace "dof" everywhere.
>>
>> Thanks for addressing silly questions.
>>
>> Ed
>>
>>
>> --
>> Ed Bueler
>> Dept of Math and Stat and Geophysical Institute
>> University of Alaska Fairbanks
>> Fairbanks, AK 99775-6660
>> 301C Chapman and 410D Elvey
>> 907 474-7693 and 907 474-7199  (fax 907 474-5394)
>>
>
>


-- 
Ed Bueler
Dept of Math and Stat and Geophysical Institute
University of Alaska Fairbanks
Fairbanks, AK 99775-6660
301C Chapman and 410D Elvey
907 474-7693 and 907 474-7199  (fax 907 474-5394)
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20160527/45765ee5/attachment-0001.html>