Strided access with array views

[ivan-pi]

I noticed that all the interpolation classes make copies of the x-, y-, and z-arrays (when applicable). This is of course a safe and pragmatic design choice, after-all memory is cheap today.

For large applications, one might not be able to afford the copies. In this case it might be preferable to have a network_view class, which only stores a reference to the data:

real, pointer :: x(:) => null(), y(:) => null()

What I've realized today is that such a class could access the coordinates in their natural storage layout:

subroutine create_mesh_2d( this, x, y )
    class(network_view_2d), intent(inout)  :: this
    real, intent(in), target          :: x(:)
    real, intent(in), target          :: y(:)

This would mean an application could support both Fortran- and C-like layouts,i.e.

    real, target :: x(4), y(4)

    x = [ 0.0, 0.0, 1.0, 1.0]
    y = [ 0.0, 1.0, 0.0, 1.0]

    call network%create_mesh( x, y )

or

    real, target :: xy(2,4)

    xy(1,:) = [ 0.0, 0.0, 1.0, 1.0]
    xy(2,:) = [ 0.0, 1.0, 0.0, 1.0]

    call network%create_mesh( xy(1,:), xy(2,:) )

I think the strided array access is a very cool Fortran feature. I haven't looked deeper into the original routines from Renka. Presumably those expect (contiguous) assumed-size arrays. Those would need to be modified to take assumed-shape arrays for the above to work

Note I have no urgent need for this, but the idea seemed worth documenting at least.

[arjenmarkus]

Definitely worth thinking about. The Renka code uses the classic assumed-size feature (the spag utility does ot take care of that), so that a copy would be made in case of non-contiguous arrays, Modifying the routines to take advantage of the assumed-shape feature is trivial. This would be a next step in modernising the code.

[arjenmarkus]

Note on spag: it leaves the routines as independent, globally visible routines, not included in modules. This is one of the manual changes I did make.

[ivan-pi]

I learned later that what I proposed above,

subroutine create_mesh_2d( this, x, y )
    class(network_view_2d), intent(inout)  :: this
    real, intent(in), target          :: x(:)
    real, intent(in), target          :: y(:)

is not safe in general. The pointer attribute has to be used here instead:

subroutine create_mesh_2d( this, x, y )
    class(network_view_2d), intent(inout)  :: this
    real, intent(in), pointer  :: x(:), y(:)

(Note the intent only applies to the pointer target, the user would still be allowed to change the data, which may break internal data structures, meaning that they would have to be recomputed.)

[arjenmarkus]

Indeed, any solution that does not copy the data is unsafe. For efficiency, an alternative interface could be provided with all those caveats.

Op do 22 aug 2024 om 10:55 schreef Ivan Pribec @.***>:

I learned later that what I proposed above,

subroutine create_mesh_2d( this, x, y ) class(network_view_2d), intent(inout) :: this real, intent(in), target :: x(:) real, intent(in), target :: y(:)

is not safe in general. The pointer attribute has to be used here instead:

subroutine create_mesh_2d( this, x, y ) class(network_view_2d), intent(inout) :: this real, intent(in), pointer :: x(:), y(:)

(Note the intent only applies to the pointer target, the user would still be allowed to change the data, which may break internal data structures, meaning that they would have to be recomputed.)

— Reply to this email directly, view it on GitHub https://github.com/arjenmarkus/interpolation2d3d/issues/1#issuecomment-2304133492, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAN6YR5DYLQAKFTKVJ626ETZSWRP5AVCNFSM6AAAAABM5VL55SVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDGMBUGEZTGNBZGI . You are receiving this because you commented.Message ID: @.***>