Stdlib: linspace and logspace

Created on 19 Dec 2019  路  8Comments  路  Source: fortran-lang/stdlib

The linspace API is implemented here: https://github.com/certik/fortran-utils/blob/b43bd24cd421509a5bc6d3b9c3eeae8ce856ed88/src/mesh.f90#L157

Matlab's linspace.

The logspace is similar, but I don't have it implemented yet -- historically I have used a function called meshexp, which is more general --- it allows you to change the gradation of the mesh, which the Matlab's logspace does not allow. The NumPy's logspace allows to set different base which allows to change gradation. So I think my meshexp can be implemented using NumPy's logspace. NumPy also has geomspace where you can specify the end points directly (just like in my meshexp) but it does not allow to change gradation. So I think there is room for meshexp, perhaps we should change the name somehow to be consistent with the other functions.

All 8 comments

Personally, for linspace I like the extended API of the numpy linspace; ported to Fortran this is:

function linspace(start,end,num,endpoint,step) result(samples)

I have created a gist with my linspace version here.

A minor downside of mirroring the numpy API, is that since the step is an intent(out) argument, meaning the function can not have the pure attribute.

Have you ever used the step argument? I've never used it.

Along the lines of "Let's use pure Fortran whenever we can" from #20, I suggest we aim for pure implementations whenever we can.

If we have an intent(out) parameter here, I'd argue that linspace should be a subroutine (a function shouldn't modify state elsewhere), which I'm not in favor of.

I have used the step argument before. A simple use case is setting up a simple finite difference method:

x = linspace(0.0_dp,1.0_dp,101,step=dx) 
! call linspace(0.0_dp,1.0_dp,101,x,dx) ! subroutine version

! ... assemble the tridiagonal matrix and rhs, dx appears in the matrix entries ... 
! ... solve system for field  u  using Thomas algorithm or gttrf and gttrs from LAPACK ...

call print_file("result.txt",x,u)

But I can agree with the pure argument of @milancurcic . If the step size is added to the function prototype, then linspace should be a subroutine.

In the end, the user can always recover the step size as follows

x = linspace(0.0_dp,1.0_dp,11)
dx = x(2)-x(1)

which I suppose is easier to remember and more intuitive than a subroutine version...

Indeed, I feel the second option:

x = linspace(0.0_dp, 1.0_dp, 11)
dx = x(2)-x(1)

is perhaps even clearer what it is doing, and it allows us to stay "pure". Matlab's linspace does not have this dx either.

How would this module be called? stdlib_numerical and stdlib_experimental_numerical?

meshgrid() from #18 would belong here as well.

In fortran-utils I call it a mesh.f90. NumPy has this directly in the numpy namespace. Matlab also seems to have this directly in the built-in namespace. For us, here are some options:

  • stdlib_mesh.f90
  • stdlib_grid.f90
  • stdlib_numerical.f90
  • stdlib_space.f90
  • stdlib_mesh_utilities.f90

I feel "numerical" is too general. I would prefer something more concrete. I like "mesh" or "grid" so far the most.

Was this page helpful?
0 / 5 - 0 ratings

Related issues

leonfoks picture leonfoks  路  9Comments

certik picture certik  路  6Comments

certik picture certik  路  7Comments

vmagnin picture vmagnin  路  10Comments

zbeekman picture zbeekman  路  4Comments