Fortran已知首尾值、步长构造数组及循环包含末值方法

问题描述 #

• 编写do循环对变量domega的每个取值执行运算时，原有代码存在缺陷：会遗漏domega的末值（该值应当等于tau2）。
• 期望实现类似Matlab中domega = first_value:step:last_value的逻辑：构造从首值到末值、固定步长的序列后，遍历所有取值执行循环。

原有缺陷代码 #

domega  = tau1                       ! 序列首值
   
  if(tau2 < 0.d0) tau2 = tau2 + twopi  ! 序列末值tau2
      
  delta    =  4.617207817725326D-003  ! 遍历步长
  n_domega = 1 + ((tau2 - tau1)/delta)     ! 区间内取值总数计算


  do i_om = 1,n_domega    ! 循环开始
      
     
 !  ( 此处执行业务计算逻辑 )
   
   
      
      domega = tau1 + i_om*delta
  enddo  ! Delta_omega循环结束

修复方案 #

通过实现类似Matlab linspace功能的子程序，预先生成包含首尾端点的等距取值数组，再遍历数组完成运算，从根源避免循环内递推计算带来的索引偏移、浮点误差导致的末值遗漏问题。

等距序列生成通用子程序 #

!***********************************************************************
    subroutine linspace_real(from, to, array)
    !***********************************************************************
    ! 生成从from到to的等距数值序列，结果包含首尾两个端点
    !
    ! 输入参数:
    ! -------
    ! from, to : 待生成序列的上下边界
    !
    ! 输出参数:
    ! -------
    ! array : 存储等距序列的数组
    !***********************************************************************
        implicit none
    
        ! 参数定义
        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
    !*********************************************************************** 

主程序调用逻辑 #

不再在循环执行过程中递增domega，而是先结合给定步长delta计算序列总长度，构造从tau1到tau2的完整domega数组后再执行遍历：

! 变量声明
real(pr), dimension(:), allocatable  :: d_omega_vec 

 
! 遍历Delta_omega区间（已知步长，先计算区间内取值总数）
domega  = tau1
                     
if(tau2 < 0.d0) tau2 = tau2 + twopi
                      
delta    =  4.617207817725326D-003  ! 遍历步长
n_domega = 1 + ((tau2 - tau1)/delta)     ! 区间内取值总数
              
n_domega_vec = n_domega + 1 
                    
allocate(d_omega_vec(n_domega_vec))
call linspace_real(tau1, tau2, d_omega_vec)

后续可直接遍历d_omega_vec数组中的所有元素执行业务计算，确保不会遗漏末值tau2。

内容的提问来源于stack exchange，提问作者g_don