如何提升Fortran代码效率？素数求解程序优化及效率方法问询

背景与现有实现 #

我为工作自学Fortran，目前是新手。完成了一道练习题：开发Fortran 90程序计算并输出2到1000之间的素数数组，结果正确，但认为存在更高效的实现方式。

初始解决方案 #

subroutine locate_primes(prime_array,count)
    integer(4), parameter               :: n = 1000
    integer                             :: i,j
    integer, allocatable, intent(out)   :: prime_array(:)
    integer, intent(out)                :: count

    j = 2
    i = 3
    count = 1
    
    do while (i<n+1) !统计素数数量
        do while (j<i)
            if (mod(i,j) == 0 ) then
                go to 100
            end if
            j = j + 1
            if (j == i) then
                count = count + 1
            end if
        end do
100     i = i+1
        j = 2
    end do
    
    allocate(prime_array(count)) !分配对应大小的数组
    count = 1
    prime_array(1) = 2
    j = 2
    i = 3

    do while (i<n+1) !生成素数数组
        do while (j<i)
            if (mod(i,j) == 0 ) then
                go to 101
            end if
            j = j + 1
            if (j == i) then
                count = count + 1
                prime_array(count)=i
            end if
        end do
101     i = i+1
        j = 2
    end do
end subroutine

问题分析：需要重复遍历3到1000的数字两次，处理效率低；且需要先统计数量再分配数组，内存利用流程冗余。

改进后的第二个方案 #

subroutine locate_primes_2(prime_array_final,count)
    integer(4), parameter               :: n = 1000
    integer                             :: i,j
    integer, allocatable, intent(out)   :: prime_array_final(:)
    integer, allocatable                :: prime_array_500(:)
    integer, intent(out)                :: count

    allocate(prime_array_500(n/2)) !创建临时数组存储素数
    j = 2
    i = 3
    count = 1
    prime_array_500(1) = 2
    
    do while (i<n+1) !统计并存储素数
        do while (j<i)
            if (mod(i,j) == 0 ) then
                go to 102
            end if
            j = j + 1
            if (j == i) then
                count = count + 1
                prime_array_500(count) = i
            end if
        end do
102     i = i+1
        j = 2
    end do
    
    allocate(prime_array_final(count))
    prime_array_final = prime_array_500(1:count) !复制有效素数到结果数组
    deallocate(prime_array_500) !释放临时数组
    
end subroutine

问题分析：无需重复遍历数字，但临时数组按n/2预先分配，存在内存冗余问题。

请教问题 #

• 此类问题的更优实现方案；
• 如何判断新方案是否优于旧方案；
• 处理效率与内存效率优化的方法及参考资料；
• 结合Fortran指针的实现思路。

问题解答 #

1. 更优实现方案 #

方案一：埃拉托斯特尼筛法（Sieve of Eratosthenes）

这是素数计算的经典高效算法，时间复杂度远低于逐个因子检查的方案。核心逻辑是先标记所有数为素数，再从最小素数开始，标记其所有倍数为非素数，最后收集未被标记的数。

subroutine sieve_of_eratosthenes(prime_array, count)
    integer, parameter :: n = 1000
    logical, dimension(n) :: is_prime
    integer, allocatable, intent(out) :: prime_array(:)
    integer, intent(out) :: count
    integer :: i, j

    ! 初始化：1不是素数，其余默认标记为素数
    is_prime = .true.
    is_prime(1) = .false.

    ! 筛法核心：标记非素数
    do i = 2, int(sqrt(real(n)))
        if (is_prime(i)) then
            ! 从i*i开始标记倍数（更小的倍数已被之前的素数标记）
            do j = i*i, n, i
                is_prime(j) = .false.
            end do
        end if
    end do

    ! 统计素数数量并分配数组
    count = count(is_prime)
    allocate(prime_array(count))

    ! 收集素数到结果数组
    j = 1
    do i = 2, n
        if (is_prime(i)) then
            prime_array(j) = i
            j = j + 1
        end if
    end do
end subroutine sieve_of_eratosthenes

该方案仅需两次遍历（标记+收集），内存上使用逻辑数组，比整数数组更节省空间。

方案二：动态扩容的素数收集

针对你想动态调整数组的需求，可通过move_alloc实现数组扩容，避免预先分配大内存：

subroutine dynamic_prime_collection(prime_array, count)
    integer, parameter :: n = 1000
    integer, allocatable, intent(out) :: prime_array(:)
    integer, intent(out) :: count
    integer :: i, j, current_size
    logical :: is_p
    integer, allocatable :: temp(:)

    count = 0
    current_size = 10  ! 初始分配小容量
    allocate(prime_array(current_size))

    ! 处理唯一的偶素数2
    count = count + 1
    prime_array(count) = 2

    ! 只检查奇数，减少计算量
    do i = 3, n, 2
        is_p = .true.
        ! 仅检查到当前数的平方根，无需遍历所有小于它的数
        do j = 1, count
            if (prime_array(j) > sqrt(real(i))) exit
            if (mod(i, prime_array(j)) == 0) then
                is_p = .false.
                exit
            end if
        end do

        if (is_p) then
            count = count + 1
            ! 数组满则扩容（翻倍）
            if (count > current_size) then
                current_size = current_size * 2
                allocate(temp(current_size))
                temp(1:count-1) = prime_array(1:count-1)
                call move_alloc(temp, prime_array)
            end if
            prime_array(count) = i
        end if
    end do

    ! 最后调整数组到实际大小
    allocate(temp(count))
    temp = prime_array(1:count)
    call move_alloc(temp, prime_array)
end subroutine dynamic_prime_collection

2. 判断方案优劣的方法 #

• 时间效率：用system_clock函数统计程序运行时间，对比不同方案的执行时长。比如在调用子routine前后分别获取时钟计数，计算差值。
• 内存效率：分析代码中数组的分配大小，对比峰值内存占用；也可通过工具查看程序运行时的内存使用情况。
• 可读性与可维护性：优先选择逻辑清晰、避免冗余跳转（如goto）的代码，便于后续修改和调试。
• 正确性：验证新方案输出的素数数组与旧方案完全一致，确保效率优化不牺牲正确性。

3. 效率优化方法与参考资料 #

处理效率优化

• 选对算法：算法的时间复杂度是效率的核心影响因素，比如筛法的时间复杂度为O(n log log n)，远优于逐个检查的O(n²)。
• 减少冗余计算：比如只检查奇数、仅遍历到目标数的平方根、筛法中从i*i开始标记倍数。
• 编译器优化：开启编译器优化选项，如Intel Fortran的-O3、Gfortran的-O3，编译器会自动进行循环展开、指令级并行等优化。
• 循环优化：将循环内的不变计算（如sqrt(n)）提到循环外；遵循Fortran列优先的内存布局，尽量按列顺序操作数组。

内存效率优化

• 按需分配：用可分配数组代替固定大小数组，只分配实际需要的内存；使用move_alloc高效转移数组内容，避免不必要的数据拷贝。
• 合适的数据类型：标记素数时用逻辑数组（logical），每个元素仅占1字节，比整数数组节省内存。
• 避免临时数组冗余：用动态扩容代替预先分配大的临时数组，或直接在原数组上操作。

参考资料

• Fortran 90/95标准文档：详细介绍可分配数组、指针等核心特性。
• 《Fortran 90程序设计》：经典入门书籍，涵盖数组操作、性能优化的实用内容。
• 编译器官方手册：Intel Fortran Compiler Documentation、Gfortran手册，提供编译器优化的详细说明。

4. Fortran指针的实现思路 #

Fortran指针可用于更灵活的内存管理，比如动态扩容时通过指针关联新分配的内存：

subroutine pointer_based_primes(prime_ptr, count)
    integer, parameter :: n = 1000
    integer, pointer, intent(out) :: prime_ptr(:)
    integer, intent(out) :: count
    integer, pointer :: temp_ptr(:)
    integer :: i, j, current_size
    logical :: is_p

    count = 0
    current_size = 10
    allocate(prime_ptr(current_size))

    ! 添加2
    count = count + 1
    prime_ptr(count) = 2

    do i = 3, n, 2
        is_p = .true.
        do j = 1, count
            if (prime_ptr(j) > sqrt(real(i))) exit
            if (mod(i, prime_ptr(j)) == 0) then
                is_p = .false.
                exit
            end if
        end do

        if (is_p) then
            count = count + 1
            if (count > current_size) then
                current_size = current_size * 2
                allocate(temp_ptr(current_size))
                temp_ptr(1:count-1) = prime_ptr(1:count-1)
                deallocate(prime_ptr)
                prime_ptr => temp_ptr  ! 更新指针指向新内存
            end if
            prime_ptr(count) = i
        end if
    end do

    ! 调整到实际大小
    allocate(temp_ptr(count))
    temp_ptr = prime_ptr(1:count)
    deallocate(prime_ptr)
    prime_ptr => temp_ptr
end subroutine pointer_based_primes

注意：Fortran指针必须指向已分配的内存，需通过=>显式关联。使用指针可灵活管理内存，但要避免悬空指针和内存泄漏，需谨慎操作。

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