How to construct an array knowing first value, last value and step size in Fortran?

Question

In the following algorithm I set up a do loop in such a way to perform some operations for each value of the variable domega. The problem of my code is that it misses the last value of domega, that should be equal to tau2.

Can you tell me how to solve this issue? In matlab I would write an array of domega like: domega = first_value:step:last_value and then I would perform a for loop on the values of domega array.

Can you show how to modify my code to get a do loop on domega, by starting from first value up to last value with the indicated step size?

  domega  = tau1                       ! first value 
   
  if(tau2 < 0.d0) tau2 = tau2 + twopi  ! tau2 is the last value
      
  delta    =  4.617207817725326D-003  ! Interval step size
  n_domega = 1 + ((tau2 - tau1)/delta)     ! numbers of values within interval


  do i_om = 1,n_domega    ! Loop 
      
     
 !  ( COMPUTATIONS...)
   
   
      
      domega = tau1 + i_om*delta
  enddo  ! end loop on Delta_omega

PROBLEM SOLVED! After reading your comments and some web searches, I managed to solve the problem by using the following linspace subroutine:

!***********************************************************************
    subroutine linspace_real(from, to, array)
    !***********************************************************************
    ! Generates evenly spaced numbers from `from` to `to` (inclusive).
    !
    ! Inputs:
    ! -------
    !
    ! from, to : the lower and upper boundaries of the numbers to generate
    !
    ! Outputs:
    ! -------
    !
    ! array : Array of evenly spaced numbers
    !***********************************************************************
        implicit none
    
        !Arguments
        real(pr), intent(in)  :: from, to
        real(pr), intent(out) :: array(:)
        
        real(pr) :: range
        integer :: n, i
    
        n = size(array)
        range = to - from

        if (n == 0) return

        if (n == 1) then
            array(1) = from
            return
        end if


        do i=1, n
            array(i) = from + range * (i - 1) / (n - 1)
        end do
        
        return
        
    end subroutine linspace_real
    !***********************************************************************

In the main program I wrote:

!Declarations
real(pr), dimension(:), allocatable  :: d_omega_vec 

 
! Scan interval of Delta_omega (we know the interval step, so we have to compute the number of values within the interval)
domega  = tau1
                     
if(tau2 < 0.d0) tau2 = tau2 + twopi
                      
delta    =  4.617207817725326D-003  ! Interval step size
n_domega = 1 + ((tau2 - tau1)/delta)     ! numbers of values within interval
              
n_domega_vec = n_domega + 1 
                    
allocate(d_omega_vec(n_domega_vec))
call linspace_real(tau1, tau2, d_omega_vec)

So, instead to increment domega during the loop, I created an array of domega from tau1 to tau2 with number of elements properly computed by exploiting my delta.

Please show the values you do see. Note that you can't expect, thanks to how floating point works, that in general the final value after adding numbers successively will exactly equal the value you desire. — francescalus, Jun 16 '22 at 13:26
See [this other question](https://stackoverflow.com/q/30707668/3157076) for ideas of the implications of doing this. — francescalus, Jun 16 '22 at 13:31
n_domega = 1191; tau1 = 0.392464539129341; tau2 = 5.89072076805024; delta= 4.617207817725326D-003. Hi @francescalus, I know the problem about do loop with real numbers, that's the reason for why I'd like to use something like the matlab command I indicated in my question, any ideas about the latter one? — g_don, Jun 16 '22 at 13:39
Last value of domega before program intreruption: domega = 5.88694184222248 — g_don, Jun 16 '22 at 13:44
Can you show that `tau1 + n_domega*delta` (exactly) equals `tau2`? (For the best answer to you, we need to understand well what you may already know.) — francescalus, Jun 16 '22 at 13:58
I think it is not possible to exactly get that equality because they are real numbers and so we can work, for example, with the difference of the two quantities and set it lower/greater a sort of "numeric zero like": `small = 1.0d-15` — g_don, Jun 16 '22 at 14:07
Problem solved! You can find my solution at the bottom of the original question. — g_don, Jun 17 '22 at 08:56

score 1 · Answer 1 · answered Jun 16 '22 at 22:18

I'm not sure if delta = 0.0026 (as in the code) or delta= 4.617207817725326D-003 (as in the comment), but here is a way to give an initial value and loop adding a step each time:

program test
  double precision, parameter :: tau1 = 0.392464539129341, tau2 = 5.89072076805024, delta = 0.0026 
  double precision :: domega

  domega  = tau1                       ! first value 
  do
    if(domega > tau2) exit
    write(*,"('domega = ',es16.9)") domega
    domega = domega + delta
  enddo  ! end loop on Delta_omega 

  write(*,"('next domega would be = ',es16.9)") domega + delta
end program test

End of output:

(...)
domega =  5.883664563E+00
domega =  5.886264563E+00
domega =  5.888864563E+00
next domega would be =  5.894064563E+00

The last line is printed outside of the loop, and is already above tau2 = 5.89072076805024.

score 1 · Answer 2 · answered Jun 17 '22 at 08:50

From what I can tell what you are looking for is basically Python's numpy.linspace. You can write effectively the same in Fortran with an implied do-loop:

array = start_val + val_spacing * [(i, i=0,N-1)]

Sample

program main
   use :: iso_fortran_env, only: dp=>real64
   implicit none

   real(dp), allocatable :: arr(:)

   arr =  linspace(1.0_dp, 2.0_dp, 11)
   print*, arr

   contains

   function linspace(a,b,N) result(lin_arr)
        !! Equivalent to Numpy's numpy.linspace

        real(dp), intent(in) :: a,b
        integer, intent(in) :: N
        real(dp), allocatable :: lin_arr(:)

        ! Local
        integer :: i
        real(dp) :: dx

        dx = (b-a) / (N-1)
        lin_arr = a + dx * [(i, i=0, N-1)]
   end function

end program main

Hi @nikjohn, thanks for your solution. I think it is very similar to mine that I added at the bottom of my original question. — g_don, Jun 17 '22 at 08:55

How to construct an array knowing first value, last value and step size in Fortran?

2 Answers2

Sample