Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Aliasing in UVINIT #7

Open
ivan-pi opened this issue Mar 10, 2023 · 1 comment
Open

Aliasing in UVINIT #7

ivan-pi opened this issue Mar 10, 2023 · 1 comment

Comments

@ivan-pi
Copy link
Owner

ivan-pi commented Mar 10, 2023

The V argument is passed as an array of size 1.

pdecheb/pdecheb/cset.f

Lines 1 to 11 in f429435

SUBROUTINE CSET(NPDE,NPTS,U,X,OMEGA,DU,XBK,NEL,NPTL,XC,CCR,XBH,
* IBK,DUTEM,V,NV)
C***********************************************************************
C FORTRAN FUNCTIONS USED: SIN COS .
C***********************************************************************
C .. Scalar Arguments ..
INTEGER IBK, NEL, NPDE, NPTL, NPTS, NV
C .. Array Arguments ..
DOUBLE PRECISIONCCR(NPTL), DU(NPTL,NPTL), DUTEM(NPTL,NPTL),
* OMEGA(NPTL,NPTL), U(NPDE,NPTS), V(1), X(NPTS),
* XBH(IBK), XBK(IBK), XC(NPTL)

...

CALL UVINIT(NPDE,NPTS,X,U,NV,V)

This can have some surprising effects, if NV = 0, but the solution array would be declared as Y(NPDE*NPTS + 1) in the calling program. Here's the initialization routine that exposes the issue:

SUBROUTINE UVINIT( NPDE, NPTS, X, U, NV, V)
IMPLICIT NONE
INTEGER :: NPDE, NPTS, NV
DOUBLE PRECISION :: X(NPTS), U(NPDE,NPTS), TIME, V(NV)
U(1,:) = 0
U(1,NPTS) = 1
V(1) = 0       ! <-- eliminates effect of previous line
END SUBROUTINE

In practice this should be flagged as an out-of-bounds violation (or produce a segfault).

@ivan-pi
Copy link
Owner Author

ivan-pi commented Mar 10, 2023

The problem is more obvious in INICHB where CSET gets called:

      NVST = NPDE*NPTS + 1
C ...
      IV = NPDE*NPTS
      IF (NV.GT.0) THEN
         IV = NVST
C ...
      END IF
      CALL CSET(NPDE,NPTS,U,WK(I6),WK,WK(I2),WK(I5),NEL,NPTL,WK(I4),
     *          WK(I12),XBK,IBK,WK(I3),U(IV),NV)

If NV = 0, the element U(IV) will be the alias for U(NPDE*NPTS) or U(1,NPTS) in the UVINIT routine.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant