Mathematics of arrays

The second lecture (and the rest) in the series will be held on Zoom:

Lenore Mullin is inviting you to a scheduled Zoom meeting.
Topic: MoA Lecture 1
Time: Sep 3, 2021 02:30 PM Eastern Time (US and Canada)
Every 2 weeks on Fri, until Nov 12, 2021, 6 occurrence(s)
Sep 3, 2021 02:30 PM
Sep 17, 2021 02:30 PM
Oct 1, 2021 02:30 PM
Oct 15, 2021 02:30 PM
Oct 29, 2021 02:30 PM
Nov 12, 2021 02:30 PM
Please download and import the following iCalendar (.ics) files to your calendar system.
Weekly: https://albany.zoom.us/meeting/tJwkfumuqT8sGdVCPn9GavN8j88zWjS399rF/ics?icsToken=98tyKuCgpzIqHNORthqGRow-HYr4c-7wtildjadrvy_rWgdSdC2uPLoaKIV1I4uJ
Join Zoom Meeting
Launch Meeting - Zoom
Meeting ID: 989 9067 4631
Passcode: 122039

2 Likes

Sorry I couldn’t make the meeting this time, too busy at work. If you could post the slides, that would be very helpful for me.

Hi @certik, yes I will post the slides sometime this weekend after some corrections (mainly typos). Also, since this time it is hosted on Zoom, we recorded it, which I will also post once it’s ready!

1 Like

@wyphan perfect, thank you! Much appreciated.

Hi everyone,

The slides for the second lecture can be accessed here:

And here’s the link for the recording:

Passcode: P2=fQfK2

2 Likes

@Lenore ,

The math with MoA is most interesting, there is a lot there to be studied by compiler writers, those working on compiler optimizations, and certain library authors.

Some questions re: Fortran and its leaning toward an “array language” since its inception and all the work you have done with MoA:

  • if you can do a quick review (or perhaps you’ve done it already) of current Fortran standard published in 2018, what comes across to you that the Fortran base language itself can do better to be a reasonably full-fledged “array language”?
  • Are there any improvements in terms of semantics or syntax or intrinsic facilities (e.g., functions, etc.) that can be introduced in the Fortran standard to enable better practice of MoA using Fortran?

A summary or a paper (a short one or a very detailed) on recommended Fortran language standard enhancements, should there be any, by you and those working with you on MoA will be highly valuable for the language.

1 Like

I would love to work with you and others on this.

When I get home I’ll respond more on Fortran discourse.

Looking at the slides, it’s not clear on p23 how to express “catenation” in Fortran. Another issue, which should be solvable, is that MoA assumes array indexing starting at 0, while Fortran arrays start at 1 by default. One can declare a lower bound other than 1, but when passing such an array to a procedure the lower bound is typically lost (example below).

Thanks to Lenore for her lectures.

module m
contains
function take(ivec,itake) result(jvec)
integer, intent(in) :: ivec(:),itake(:)
integer             :: jvec(size(itake))
jvec = ivec(itake)
end function take
end module m

program main
use m
integer :: ivec(0:2)
ivec = [10,20,30]
print*,take(ivec,[1,2]) ! gives [10,20]
end program main

@Beliavsky Yes, I encountered the 1-based indexing in Fortran to be an issue when testing the Fortran analog for the “take” operation (p.19). Thankfully the fix is pretty simple, i.e. manually add +1.
For instance, here’s the “take” operation implemented as a pointer function for the default real data type.

PROGRAM taketest
  IMPLICIT NONE

  INTEGER, PARAMETER :: n = 10
  INTEGER :: i, s
  REAL, TARGET :: v(0:n-1)
  REAL, POINTER :: p(:)

  ! Init v
  v = [( 0.1 * REAL(i), i = 0, n-1 )]
  PRINT *, 'vvec = '
  PRINT *, ( v(i), ' ' , i = 0, n-1 )

  ! Positive take
  s = 2
  p => take( s, v )
  PRINT *, s, ' take vvec = '
  PRINT *, ( p(i), ' ' , i = 1, s )
  NULLIFY(p)

  ! Negative take
  s = -2
  p => take( s, v )
  PRINT *, s, ' take vvec = '
  PRINT *, ( p(i), ' ' , i = 1, ABS(s) )
  NULLIFY(p)

  STOP

CONTAINS

  FUNCTION take( sigma, vec )
    IMPLICIT NONE

    ! Arguments
    INTEGER, INTENT(IN) :: sigma
    REAL, INTENT(IN), TARGET :: vec(:)
    REAL, POINTER :: take(:)

    ! Internal variables
    INTEGER :: n

    n = SIZE( vec, 1 )

    IF( ABS(sigma) <= n ) THEN

      IF( sigma > 0 ) THEN

        ! Positive take
        take => vec( :sigma )

      ELSE

        ! Negative take
        ! Note how +1 is needed here to account for 1-based Fortran indexing
        take => vec( (n+sigma+1): )

      END IF

    ELSE

      PRINT *, 'Error[take]: out of bounds ', sigma, ' of ', n

    END IF

  END FUNCTION take

END PROGRAM taketest

@wyphan , @Beliavsky , and especially anyone interested in library development:

I would “debate” instead the flexibility provided in the Fortran standard with assumed-shape arrays is “good enough” for the needs here.

The “take” procedure truly belongs in a library and under the circumstances, that library can have named constant for the desired lower bound (zero-based or whatever) and work with assumed-shaped arrays accordingly.

integer, parameter :: LB = 0 ! Desired lower bound
..
function take (..)
..
real, .., intent(in) :: vec(LB:) ! Set the lower bound
..

While the above may fail to tickle some tummies, it is hardly the painful barriers in 0-based indexing languages (Djikstra-inspired ones) when it comes to a library author wanting to use a lower bound other than that hardwired into the language.

With Fortran language standard and MoA, the question instead is whether something deeper is lacking.

1 Like

@certik @gak and other Fortran experts who know the standard : we could do this together and publish to help others get on board.
BTW index origin is arbitrary in a language. The theory is behind the scenes.

@Beliavsky Sorry I missed the question on catenation. Yes, I don’t have any idea yet on how to make it work properly for cases such as catenating a vector with a scalar.

1 Like

I am picking up this challenge :slight_smile:

1 Like

@Arjen Please DM me your GitHub username, so I can add you to a common repo that I made for the MoA Fortran Library

Okay, it took me some time to acutally sit down for it, but I have a (slightly disfunctional) example of how you could do the catenate operation (something is still wrong, need to figure that out). Here is the code:

! moa_cat.f90 --
!     First version of a cat operation
!
module moa_operations
    implicit none

    !
    ! Type for defining "views" on arrays - to allow for a#b#c: pointers to other views
    !
    type view_type
        integer, dimension(:), pointer :: left_array  => null()
        integer, dimension(:), pointer :: right_array => null()
        type(view_type), pointer       :: left_view   => null()
        type(view_type), pointer       :: right_view  => null()
    contains
    !    procedure, nopass :: cataa => catenate_array_array
    !    procedure, nopass :: catav => catenate_array_view
    !    procedure, nopass :: catva => catenate_view_array
    !    procedure, nopass :: catvv => catenate_view_view
         procedure         :: elem  => get_elem
    !    generic   :: cat => cataa, catav, catva, catvv
    end type view_type

    interface operator(//)
        module procedure catenate_array_array, catenate_array_view, catenate_view_array, catenate_view_view
    end interface
contains
function catenate_array_array( array1, array2 ) result(new_view)
    integer, dimension(:), target, intent(in) :: array1
    integer, dimension(:), target, intent(in) :: array2
    class(view_type), allocatable             :: new_view

    allocate( new_view )

    new_view%left_array  => array1
    new_view%right_array => array2
end function catenate_array_array

function catenate_view_array( view1, array2 ) result(new_view)
    class(view_type), target, intent(in)      :: view1
    integer, dimension(:), target, intent(in) :: array2
    class(view_type), allocatable             :: new_view

    allocate( new_view )

    new_view%left_view   => view1
    new_view%right_array => array2
end function catenate_view_array

function catenate_array_view( array1, view2 ) result(new_view)
    integer, dimension(:), target, intent(in) :: array1
    class(view_type), target, intent(in)      :: view2
    class(view_type), allocatable             :: new_view

    allocate( new_view )

    new_view%left_array => array1
    new_view%right_view => view2
end function catenate_array_view

function catenate_view_view( view1, view2 ) result(new_view)
    class(view_type), target, intent(in) :: view1
    class(view_type), target, intent(in) :: view2
    class(view_type), allocatable        :: new_view

    allocate( new_view )

    new_view%left_view  => view1
    new_view%right_view => view2
end function catenate_view_view

recursive function get_elem(view, i) result(elem)
    class(view_type), intent(inout) :: view
    integer, intent(in)             :: i
    integer, pointer                :: elem

    integer                         :: sz

    if ( associated(view%left_array) ) then
        sz = size(view%left_array)
        if ( i <= sz ) then
            elem => view%left_array(i)
        else
            if ( associated(view%right_array) ) then
                elem => view%right_array(i-sz)
            else
                elem => view%right_view%elem(i-sz)
            endif
        endif
    else
        elem => view%left_view%elem(i)
    endif
end function
end module moa_operations

program test_cat
    use moa_operations
    implicit none

    type(view_type)        :: x, y
    integer, dimension(10) :: a, b, c
    integer                :: i, j
    !     1   2   3   4    5    6    7    8    9   10
    a = [ 2,  3,  5,  7,  11,  13,  17,  19,  23,  29]
    b = [-2, -3, -5, -7, -11, -13, -17, -19, -23, -29]
    c = [31, 37, 41, 43,  47,  53,  59,  61,  67,  71]

    x = a // b
    y = x // c

    do i = 1,30
        write(*,*) i, y%elem(i)
    enddo
end program test_cat

In this branch you still need to determine if i <= size of the left view. I think a recursive procedure to calculate the size of a view is needed.

if (associated(view%left_array) then
  if (associated(view%right_array) then
    view_size = size(view%left_array) + size(view%right_array)
  else
    view_size = size(view%left_array) + view%right_view%size()
  end if
else
  if (associated(view%right_array) then
    view_size = view%left_view%size() + size(view%right_array)
  else
    view_size = view%left_view%size() + view%right_view%size()
  end if
end if

An additional thought, should you check for index out of bounds and error stop? Or should this behave like normal arrays and only crash if you’ve got bounds checking on?

@everythingfunctional , yes, I found that out myself. Trying to wrap my poor brains around this, while als o doing my regular job :wink:

1 Like

Okay, it took a bit serious thinking and debugging, but here is an improved version (not extensively tested, mind you):

! moa_cat.f90 --
!     First version of a cat operation
!
module moa_operations
    implicit none

    !
    ! Type for defining "views" on arrays - to allow for a#b#c: pointers to other views
    !
    type view_type
        integer, dimension(:), pointer :: left_array  => null()
        integer, dimension(:), pointer :: right_array => null()
        type(view_type), pointer       :: left_view   => null()
        type(view_type), pointer       :: right_view  => null()
    contains
    !    procedure, nopass :: cataa => catenate_array_array
    !    procedure, nopass :: catav => catenate_array_view
    !    procedure, nopass :: catva => catenate_view_array
    !    procedure, nopass :: catvv => catenate_view_view
         procedure         :: elem  => get_elem
    !    generic   :: cat => cataa, catav, catva, catvv
    end type view_type

    interface operator(//)
        module procedure catenate_array_array, catenate_array_view, catenate_view_array, catenate_view_view
    end interface
contains
function catenate_array_array( array1, array2 ) result(new_view)
    integer, dimension(:), target, intent(in) :: array1
    integer, dimension(:), target, intent(in) :: array2
    class(view_type), allocatable             :: new_view

    allocate( new_view )

    new_view%left_array  => array1
    new_view%right_array => array2
end function catenate_array_array

function catenate_view_array( view1, array2 ) result(new_view)
    class(view_type), target, intent(in)      :: view1
    integer, dimension(:), target, intent(in) :: array2
    class(view_type), allocatable             :: new_view

    allocate( new_view )

    new_view%left_view   => view1
    new_view%right_array => array2
end function catenate_view_array

function catenate_array_view( array1, view2 ) result(new_view)
    integer, dimension(:), target, intent(in) :: array1
    class(view_type), target, intent(in)      :: view2
    class(view_type), allocatable             :: new_view

    allocate( new_view )

    new_view%left_array => array1
    new_view%right_view => view2
end function catenate_array_view

function catenate_view_view( view1, view2 ) result(new_view)
    class(view_type), target, intent(in) :: view1
    class(view_type), target, intent(in) :: view2
    class(view_type), allocatable        :: new_view

    allocate( new_view )

    new_view%left_view  => view1
    new_view%right_view => view2
end function catenate_view_view

recursive subroutine get_pointer(view, i, elem, found)
    class(view_type), intent(inout) :: view
    integer, intent(inout)          :: i
    integer, pointer                :: elem
    logical, intent(out)            :: found

    integer                         :: sz


    found = .false.

    if ( associated(view%left_array) ) then
        sz = size(view%left_array)
        if ( i <= sz ) then
            found =  .true.
            elem  => view%left_array(i)
        else
            i = i - sz
            if ( associated(view%right_array) ) then
                sz = size(view%right_array)
                if ( i <= sz ) then
                    found =  .true.
                    elem  => view%right_array(i)
                else
                    i = i - sz
                    return
                endif
            else
                call get_pointer( view%right_view, i, elem, found )
            endif
        endif
    else
        call get_pointer(view%left_view, i, elem, found )

        if ( .not. found ) then
            if ( associated(view%right_array) ) then
                sz = size(view%right_array)
                if ( i <= sz ) then
                    found =  .true.
                    elem  => view%right_array(i)
                else
                    error stop
                endif
            else
                call get_pointer( view%right_view, i, elem, found )
            endif
        endif
    endif
end subroutine get_pointer

function get_elem(view, i) result(elem)
    class(view_type), intent(inout) :: view
    integer, intent(in)             :: i
    integer, pointer                :: elem

    integer                         :: inew
    logical                         :: found


    inew = i
    call get_pointer( view, inew, elem, found )
end function
end module moa_operations

program test_cat
    use moa_operations
    implicit none

    type(view_type)        :: x, y
    integer, dimension(10) :: a, b, c
    integer                :: i, j
    !     1   2   3   4    5    6    7    8    9   10
    a = [ 2,  3,  5,  7,  11,  13,  17,  19,  23,  29]
    b = [-2, -3, -5, -7, -11, -13, -17, -19, -23, -29]
    c = [31, 37, 41, 43,  47,  53,  59,  61,  67,  71]

    x = a // b
    y = x // c

    do i = 1,30
        write(*,*) i, y%elem(i)
    enddo
end program test_cat
2 Likes

Google’s Dex language looks interesting. Is it somehow related to MoA?