[stdlib_linalg] Add empty function.

[zoziha]

• add empty function. (see numpy/empty)

Tasks

• ready to allocatable solution. (see [stdlib_linalg] Add empty function.  #477 (comment))

More routines to do, not this PR

see #476.

[awvwgk]

The main purpose of numpy.empty would be to create an allocation on the heap and return its reference. For a similar mechanism in Fortran we would have to use the pointer attribute, also we cannot assign in such case:

implicit none
real, pointer :: array(:)
integer :: n = 6000000

array => empty_ptr(n)
! ... do something with array
deallocate(array)  ! explicit free required
contains
   function empty_ptr(n) result(array)
      integer, intent(in) :: n
      real, pointer :: array(:)
      allocate(array(n))
   end function
end

The current implementation uses an automatic array to trigger an automatic LHS (re)allocation. This has several drawback, first the automatic array will most likely be created on the stack and cannot be moved to an allocatable array on the heap, but must be copied to the newly allocated array. This will be considerably worse than just using allocate. Also large stack arrays allocated by this mean are at risk of overflowing the stack, which leads to hard to debug stack overflows.

Preferably, to match the original intent of the functionality provided by numpy.empty to provide a functional access to the allocate procedure without introducing additional overhead, the implementation would be along the lines of:

function empty(n) result(array)
  integer, intent(in) :: n
  real, allocatable :: array(:)
  allocate(array(n))
end function empty

In the best case the compiler might move the allocation produced by the function to the LHS array instead of separately allocating and copying the uninitialized array.

[zoziha]

In the best case the compiler might move the allocation produced by the function to the LHS array instead of separately allocating and copying the uninitialized array.

Take a look at these two examples @awvwgk , I see from examples here, the current code implementation is more efficient in gfortran (not ifort), and the code implementation method you mentioned can also be (re)allocated.

Example 1

In the case of small stack usage, gfortran is more efficient and ifort is less efficient for the current implementation. (Um..)

program test

    real :: stime, etime
    real(kind=8), allocatable :: A(:,:), B(:,:)

    call cpu_time(stime)
    A = empty1(100,100)
    A = empty1(200,200)
    call cpu_time(etime)
    print *, "etime - stime (seconds) : ", etime - stime !! 3.4e-4 (gfortran); 1.5e-5 (ifx); 5.5e-4 (ifort)

    call cpu_time(stime)
    B = empty2(100,100)
    B = empty2(200,200)
    call cpu_time(etime)
    print *, "etime - stime (seconds) : ", etime - stime !! 8.0e-6 (gfortran); 4.3e-4 (ifx); 4.5e-4 (ifort)

contains

    pure function empty1(ndim1, ndim2) result(result)
        implicit none
        integer, intent(in) :: ndim1, ndim2
        real(kind=8), allocatable :: result(:,:)
        allocate(result(ndim1, ndim2))
    end function empty1

    pure function empty2(ndim1, ndim2) result(result)
        implicit none
        integer, intent(in) :: ndim1, ndim2
        real(kind=8) :: result(ndim1, ndim2)
    end function empty2

end program test

Example 2

In the case of large stack usage, gfortran is more efficient as well and ifort stack overflows for the current implementation. (Um.. +1)

program test

    real :: stime, etime
    real(kind=8), allocatable :: A(:,:), B(:,:)

    call cpu_time(stime)
    A = empty1(1000,1000)
    A = empty1(2000,2000)
    call cpu_time(etime)
    print *, "etime - stime (seconds) : ", etime - stime !! 3.3e-2 (gfortran); 2.2e-5 (ifx); 5.3e-2 (ifort)

    call cpu_time(stime)
    B = empty2(1000,1000)
    B = empty2(2000,2000)
    call cpu_time(etime)
    print *, "etime - stime (seconds) : ", etime - stime !! 1.6e-5 (gfortran); stack overflow (ifx/ifort) 

contains

    pure function empty1(ndim1, ndim2) result(result)
        implicit none
        integer, intent(in) :: ndim1, ndim2
        real(kind=8), allocatable :: result(:,:)
        allocate(result(ndim1, ndim2))
    end function empty1

    pure function empty2(ndim1, ndim2) result(result)
        implicit none
        integer, intent(in) :: ndim1, ndim2
        real(kind=8) :: result(ndim1, ndim2)
    end function empty2

end program test

[zoziha]

I think based on the above examples, to be compatible with ifort's compiler and maintain the robustness of empty, we use the allocatable solution seems to be a better choice.

[milancurcic]

Suggested change

	program demo_linlag_empty_1
	use stdlib_linlag, only: empty
	program demo_linalg_empty_1
	use stdlib_linalg, only: empty

[milancurcic]

Like with eye, I don't think that the allocatable result is a better choice here over the automatic array and the examples show it. However, the difference is probably low impact because this is likely to be used as a convenience function rather than a high-performance one, and in my view UX > performance. I think this can go forward as is.

[awvwgk]

The only use case I see for this function is to trigger an automatic LHS allocation, which I think is rather limited compared to the intrinsic function, because it will always require a statement to make use of this.

allocate(array(n, n))
array = empty(n, n)

Using it in an expression would only work when multiplying by zero, because the initial entries are undefined, but for this purpose we could use zeros instead. Also, in this context the result could change depending on compiler options, e.g. when requesting the compiler to initialize all reals with signalling NaNs.

I can see why this function is required in numpy to return a pointer to a memory allocation, but I don't think there is need for it in Fortran.

[zoziha]

I think the purpose of empty is indeed not clear at present, and I agree to cancel this PR.
We can just use allocate(array(m, n)) in Fortran.