串行程序如何改并行

qinliang_1985

下面是求二阶duffing方程幅频特性的程序，主要调用了imsl库里的ivmrk进行积分计算，先是选定频率，然后调用瞬态计算子程序使响应稳定，接着稳态计算一定周期取最大振幅，程序串行执行没有问题，为了加快计算速度，在频率的外层do循环加了并行，但始终报错，错误出现在ivmrk的内部计算中，求助大神们指点！

提示：
*** TERMINAL ERROR 38 from Di2mrk.  I4MRK has not been called so I7MRK
***          cannot be used.
     Here is a traceback of subprogram calls in reverse order:
     Routine name                    Error type  Error code
     ------------                    ----------  ----------
     Di2mrk                               5          38    (Called internally)
     Divmrk                               0           0
     USER                                 0           0

[Fortran] 纯文本查看 复制代码

module shared_vars
    implicit none
    integer, parameter :: dp = kind(0.0d0)
    ! 全局变量
    integer, parameter :: step = 200  ! 每周期步数
    integer, parameter :: transient_period = 500  ! 瞬态仿真周期 (稳定后提前终止)
    integer, parameter :: stable_period = 10  ! 稳态仿真周期
    ! 预先分配稳态仿真所需的位移数组
    integer, parameter :: steady_state_steps = stable_period * step
    real(dp), dimension(steady_state_steps) :: x_values
    ! 外部激励频率和模型参数
    real(dp), parameter :: alpha = 0.4_dp, beta = 0.5_dp, delta = -4.0_dp, F = 3.0_dp
    real(dp) :: p
end module shared_vars

program duffingfreqanalysis
    include 'link_fnl_shared.h'
    use shared_vars
    use ivmrk_int
    implicit none     
    integer, parameter :: n = 2  ! 方程数量 (Duffing模型中为2个方程)    
    real(dp) :: t, y(n)
    real(dp) :: period ! 周期
    real(dp) :: max_displacement
    real(dp) :: p_start, p_end, dp_p  ! 频率p的范围和步长
    real(dp) :: start_time, end_time, elapsed_time  ! 记录时间的变量
    external :: fcn  ! 外部函数
    integer :: i, num_steps

    ! 输出文件的文件名
    open(unit=10, file='amplitude_vs_frequency.txt', status='unknown')            

    ! 设置频率p的范围和步长
    p_start = 0.5_dp
    p_end = 4.0_dp
    dp_p = 0.01_dp
    num_steps = int((p_end - p_start) / dp_p) + 1

    ! 开始计时
    call cpu_time(start_time)

    !$omp parallel do private(i, p, period, t, y) shared(p_start, dp_p, num_steps)
    do i = 1, num_steps
        p = p_start + (i - 1) * dp_p   

        ! 计算当前频率对应的周期
        period = 2 * 3.141592653589793_dp / p       
        ! 调用瞬态计算子程序
        call transient_simulation(t, y, period)        
        ! 调用稳态计算子程序
        call steady_state_simulation(t, y, period, max_displacement)
        ! 输出初始结果到文件
        !$omp critical
        write (10, '(2F12.6)') p, max_displacement    
        !$omp end critical
    end do
    !$omp end parallel do

    ! 关闭文件
    close (unit=10)  

    ! 结束计时
    call cpu_time(end_time)
    ! 计算运行时间
    elapsed_time = end_time - start_time
    ! 输出运行时间到控制台
    print *, '运行时间 (秒):', elapsed_time

contains

    subroutine transient_simulation(t, y, period)
        use shared_vars
        implicit none 
        real(dp), intent(in) :: period
        real(dp) :: t, tend, y(n), yprime(n)
        integer :: ido, istep

        ! 初始条件设定
        t = 0.0
        tend = 0.0
        y(1) = 0.0   ! 初始位移
        y(2) = 0.0   ! 初始速度

        ! 初始调用ivmrk
        ido = 1        
        ! 瞬态仿真循环
        do istep = 1, transient_period            
            ! 更新 tend
            tend = tend + period  ! 更新结束时间           
            ! 调用ivmrk函数进行积分
            call ivmrk(ido, fcn, t, tend, y, yprime)           

        end do  

        ! 释放空间       
        ido = 3     
        call ivmrk(ido, fcn, t, tend, y, yprime)

    end subroutine transient_simulation

    subroutine steady_state_simulation(t, y, period, max_displacement)
        use shared_vars
        implicit none 
        real(dp), intent(in) :: period
        real(dp) :: t, tend, y(n), yprime(n), h
        integer :: ido, istep
        real(dp), intent(out) :: max_displacement

        ! 初始条件设定
        tend = t
        h = period / step

        ! 初始调用ivmrk
        ido = 1        
        ! 循环调用ivmrk，计算每个时间步的结果
        do istep = 1, steady_state_steps            
            ! 更新 tend
            tend = tend + h  ! 更新当前时间        
            ! 调用ivmrk函数进行积分
            call ivmrk(ido, fcn, t, tend, y, yprime)           
            ! 存储位移到数组
            x_values(istep) = y(1)  ! y(1)表示位移   
        end do       
        ! 释放空间
        ido = 3
        call ivmrk(ido, fcn, t, tend, y, yprime)  

        ! 计算最大位移的绝对值
        max_displacement = maxval(abs(x_values))

    end subroutine steady_state_simulation

end program duffingfreqanalysis

subroutine fcn(n, t, y, yprime)
    use shared_vars
    implicit none 
    integer, intent(in) :: n
    real(dp), intent(in) :: t
    real(dp), intent(in) :: y(n)
    real(dp), intent(out) :: yprime(n)

    ! 计算 Duffing 方程的导数
    yprime(1) = y(2)
    yprime(2) = -alpha * y(2) - beta * y(1)**3 + delta * y(1) + F * cos(p * t)

end subroutine fcn