[petsc-users] VecPow clarification

Barry Smith bsmith at petsc.dev
Thu Nov 14 09:29:17 CST 2024


  I see that currently VecPow is used only in a small number of places including:

  ksp/ksp/utils/lmvm/diagbrdn/diagbrdn.c:          PetscCall(VecPow(ldb->U, ldb->beta - 1));


  I am unsure if  the usage here requires the special handling of negative numbers.

  I was wrong and it is not used in the semi-smooth code, that access the vector elements directly.

  If could be we can strip out all the special infinity cases completely.


  Barry


> On Nov 14, 2024, at 9:56 AM, Stefano Zampini <stefano.zampini at gmail.com> wrote:
> 
> That is a very old bug! Can you make an MR to just call PetscPowScalar in a loop here https://urldefense.us/v3/__https://gitlab.com/petsc/petsc/-/blob/main/src/vec/vec/utils/projection.c*L1022__;Iw!!G_uCfscf7eWS!b2Nfr7i01ZnyfVjweWqi87it8hcCDv0s6MtIsQhmSU8s4dq4jBPi-Ca87RRTwS20Srwyh9wQdhXcsbR6MVC8WSU$  <https://urldefense.us/v3/__https://gitlab.com/petsc/petsc/-/blob/main/src/vec/vec/utils/projection.c*L1022__;Iw!!G_uCfscf7eWS!bbSdSMnU5KpH03jHI7aV5j4WLGQ3yPvxWzR6Lwr14QLy_7EJ2MNT-qhL6J1x6z3vpF6M5GQk9lMQkMoLJ_JNSn1mqwqict4$> ?
> 
> Thanks
> 
> Il giorno gio 14 nov 2024 alle ore 17:39 Peder Jørgensgaard Olesen via petsc-users <petsc-users at mcs.anl.gov <mailto:petsc-users at mcs.anl.gov>> ha scritto:
>> Given a vector containing roots of unity, v[i] = exp(i*k*x[i]) I wanted to compute the vector u[i]=exp(i*n*k*x[i]), for some real number n. From the face of it this should be easily achieved with VecPow, as u[i] = v[i]^n.
>> 
>> That didn't work as expected, though I got around it using VecGetArray() and a loop with PetscPowComplex(). The source designated in the docs (src/vec/vec/utils/projection.c) reveals that VecPow() maps v[i] to PETSC_INFINITY when the PetscRealPart(v[i]) < 0, unless the power is any of 0, ±0.5, ±1 or ±2. Even in the simple case of a purely real vector (with negative entries) raised to any other integer power, the results would not be what one might  reasonably expect from the description of VecPow().
>> 
>> While I do have a solution suiting my need, I'm left wondering what might be the rationale for VecPow working the way it does.
>> 
>> Best,
>> Peder
> 
> 
> 
> --
> Stefano

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20241114/971492ff/attachment.html>


More information about the petsc-users mailing list