• Hacking the pointer (Re: Mathematics of arrays

    From Jos Bergervoet@21:1/5 to Jos Bergervoet on Fri Aug 20 17:46:22 2021
    [One extra comment line in the program is needed: ]

    On 21/08/20 5:27 PM, Jos Bergervoet wrote:
    On 21/08/20 1:54 AM, gah4 wrote:
    On Thursday, August 19, 2021 at 2:08:04 PM UTC-7, Jos Bergervoet wrote:
    On 21/08/16 4:24 AM, Robin Vowels wrote:

    (snip)
    It does.
    Try LBOUND and UBOUND.

    That's only part of the array descriptor metadata (the stride still
    missing..) And it gives the user only read access.

    Yes, it would be interesting to have an extension, or alternative, to
    C_F_POINTER that would allow one to specify the strides separarately,

    What we want is a transposed pointer assignment, like in:

    !-- Warning, dangerous hacking, tested only for gfortran 9.1.0
    program swap_strides
      implicit none
      integer, target   :: A(4,4), B(4,4), C(4,4)
      integer, pointer  :: P(:,:)

      A(1,:) = [1 ,2 ,3 ,4]
      A(2,:) = [0 ,5 ,6 ,7]
      A(3,:) = [0 ,0 ,8 ,9]
      A(4,:) = [0 ,0 ,0 ,10]

      call printmat(A)   ! The original matrix

    !-- Different ways to get transpose without copying:
      P => transptr(A);  call printmat(P)

      P => transptr(B); P = A;  call printmat(B)


    !-- Now with copying, but without intrinsic 'transpose':

      call printmat( transptr(A) )

      transptr(C) = A;   call printmat(C)

    !-- More complicated mirroring and transposing:
      P => transptr( A( :, 4:1:-1 ) );  call printmat(P)

      call printmat(transptr( A( :, 4:1:-1 ) ))

    !-- This one does not work: (optimization dependent?)
      B = transptr(A);  call printmat(B)  ! Error !!!


    contains

      subroutine printmat(m)
        integer, intent(in)  :: m(:,:)
        integer              :: i

        print *; do i=1,ubound(m,1); print "(19i4)", m(i,:); end do
      end


    !-- Pointer hacking function using descriptor wrapper type:

      function transptr(m)       result(p)
        integer, target, intent(in)  :: m(:,:)
        integer, pointer             :: p(:,:)
        integer                      :: ar(22)

        type data; integer, pointer :: p(:,:); end type;  type(data) :: d

        d%p => m;   ar = transfer(d, ar)     ! fetch descriptor in ar(:)

        ar(11:22) = [ar(17:22), ar(11:16)]   ! swap index information

        d = transfer(ar, d)                  ! put back in the descriptor

        p => d%p                             ! copy pointer from wrapper
      end

    end program


    such that one could arrange for all the different ways of addressing
    array elements.

    The above could perhaps be extended, but it is of course a very
    unreliable way of hacking, so the real challenge is to propose
    new syntax and make it part of the language standard (after which
    it will of course be perfectly reliable, I mean!)

    As shown in the above, gfortran apparently uses 22 integers in
    the descriptor of a matrix. At least for gcc 9.1.0.


    --- SoupGate-Win32 v1.05
    * Origin: fsxNet Usenet Gateway (21:1/5)
  • From Jos Bergervoet@21:1/5 to All on Fri Aug 20 17:27:44 2021
    On 21/08/20 1:54 AM, gah4 wrote:
    On Thursday, August 19, 2021 at 2:08:04 PM UTC-7, Jos Bergervoet wrote:
    On 21/08/16 4:24 AM, Robin Vowels wrote:

    (snip)
    It does.
    Try LBOUND and UBOUND.

    That's only part of the array descriptor metadata (the stride still
    missing..) And it gives the user only read access.

    Yes, it would be interesting to have an extension, or alternative, to C_F_POINTER that would allow one to specify the strides separarately,

    What we want is a transposed pointer assignment, like in:

    !-- Warning, dangerous hacking, tested only for gfortran 9.1.0
    program swap_strides
    implicit none
    integer, target :: A(4,4), B(4,4), C(4,4)
    integer, pointer :: P(:,:)

    A(1,:) = [1 ,2 ,3 ,4]
    A(2,:) = [0 ,5 ,6 ,7]
    A(3,:) = [0 ,0 ,8 ,9]
    A(4,:) = [0 ,0 ,0 ,10]

    call printmat(A) ! The original matrix

    !-- Different ways to get transpose without copying:
    P => transptr(A); call printmat(P)

    P => transptr(B); P = A; call printmat(B)

    call printmat( transptr(A) )

    transptr(C) = A; call printmat(C)

    !-- More complicated mirroring and transposing:
    P => transptr( A( :, 4:1:-1 ) ); call printmat(P)

    call printmat(transptr( A( :, 4:1:-1 ) ))

    !-- This one does not work: (optimization dependent?)
    B = transptr(A); call printmat(B) ! Error !!!


    contains

    subroutine printmat(m)
    integer, intent(in) :: m(:,:)
    integer :: i

    print *; do i=1,ubound(m,1); print "(19i4)", m(i,:); end do
    end


    !-- Pointer hacking function using descriptor wrapper type:

    function transptr(m) result(p)
    integer, target, intent(in) :: m(:,:)
    integer, pointer :: p(:,:)
    integer :: ar(22)

    type data; integer, pointer :: p(:,:); end type; type(data) :: d

    d%p => m; ar = transfer(d, ar) ! fetch descriptor in ar(:)

    ar(11:22) = [ar(17:22), ar(11:16)] ! swap index information

    d = transfer(ar, d) ! put back in the descriptor

    p => d%p ! copy pointer from wrapper
    end

    end program


    such that one could arrange for all the different ways of addressing
    array elements.

    The above could perhaps be extended, but it is of course a very
    unreliable way of hacking, so the real challenge is to propose
    new syntax and make it part of the language standard (after which
    it will of course be perfectly reliable, I mean!)

    As shown in the above, gfortran apparently uses 22 integers in
    the descriptor of a matrix. At least for gcc 9.1.0.

    --
    Jos

    --- SoupGate-Win32 v1.05
    * Origin: fsxNet Usenet Gateway (21:1/5)
  • From gah4@21:1/5 to All on Fri Aug 20 18:30:07 2021
    On Friday, August 20, 2021 at 8:28:04 AM UTC-7, Jos Bergervoet wrote:

    (snip, I wrote)
    Yes, it would be interesting to have an extension, or alternative, to C_F_POINTER that would allow one to specify the strides separarately,

    What we want is a transposed pointer assignment, like in:

    (snip)

    The above could perhaps be extended, but it is of course a very
    unreliable way of hacking, so the real challenge is to propose
    new syntax and make it part of the language standard (after which
    it will of course be perfectly reliable, I mean!)

    There could be other ways, but C_F_POINTER is one way to create a Fortran pointer to existing data. It would seem a not so big extension to allow making pointers, such as for a transposed array or triangular matrix.

    --- SoupGate-Win32 v1.05
    * Origin: fsxNet Usenet Gateway (21:1/5)
  • From Jos Bergervoet@21:1/5 to All on Sat Aug 21 23:48:10 2021
    On 21/08/21 3:30 AM, gah4 wrote:
    On Friday, August 20, 2021 at 8:28:04 AM UTC-7, Jos Bergervoet wrote:

    (snip, I wrote)
    Yes, it would be interesting to have an extension, or alternative, to
    C_F_POINTER that would allow one to specify the strides separarately,

    What we want is a transposed pointer assignment, like in:

    (snip)

    The above could perhaps be extended, but it is of course a very
    unreliable way of hacking, so the real challenge is to propose
    new syntax and make it part of the language standard (after which
    it will of course be perfectly reliable, I mean!)

    There could be other ways, but C_F_POINTER is one way to create a Fortran pointer to existing data.

    It cannot only make Fortran pointers to contiguous arrays, so it
    is not even capable of making all possible pointers that can exist
    in Fortran as it is.

    ... It would seem a not so big extension to allow making
    pointers, such as for a transposed array

    You would have to add more optional arguments, to avoid being
    restricted to column-major ordering. (And probably remove the
    contiguous restriction as well, to allow all pointers.) So
    what additional arguments would you propose?

    or triangular matrix.

    That's not possible with only a "not so big extension" to
    C_F_POINTER, because such pointers do not exist in Fortran,
    and the compilers cannot handle them. There is no intrinsic
    sparse-array handling (of which triangular matrices would
    just be one possible case, ragged arrays or more complicated
    cases could then also be considered..)

    New Fortran syntax for declaration of such arrays and for
    handling them in allocations and array-expressions would
    be needed. (Whether C_F_POINTER would be of much help with
    that is questionable..)

    --
    Jos

    --- SoupGate-Win32 v1.05
    * Origin: fsxNet Usenet Gateway (21:1/5)