Fortran高级数组技术探讨 - 摩柯技术社区

Fortran 数组的基础回顾

在深入探讨 Fortran 的高级数组技术之前，我们先来回顾一下 Fortran 数组的基础知识。Fortran 中的数组是一种有序的数据集合，所有元素具有相同的数据类型。

数组声明

声明数组的基本语法如下：

type, dimension([lower_bound: ]upper_bound [, [lower_bound: ]upper_bound...]) :: array_name

例如，声明一个一维整数数组 a，长度为 10：

integer, dimension(10) :: a

这里默认下限为 1。如果想要指定下限，例如从 0 开始，可这样声明：

integer, dimension(0:9) :: a

对于二维数组，比如一个 3×4 的实数数组 b：

real, dimension(3, 4) :: b

第一个维度表示行，第二个维度表示列。

数组初始化

静态初始化 在声明数组时可以同时进行初始化，例如：

integer, dimension(3) :: a = [1, 2, 3]

对于二维数组：

integer, dimension(2, 2) :: b = reshape([1, 2, 3, 4], [2, 2])

这里使用 reshape 函数将一维数组重新塑造为二维数组。

动态初始化 在程序执行过程中对数组进行赋值，例如：

integer, dimension(5) :: a
do i = 1, 5
    a(i) = i * 2
end do

高级数组索引

数组切片

数组切片允许我们访问数组的一部分，而不是单个元素。在 Fortran 中，切片的语法非常直观。

一维数组切片 对于一维数组 a，获取其第 3 到第 5 个元素：

integer, dimension(10) :: a = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
integer, dimension(3) :: sub_a
sub_a = a(3:5)

这里 a(3:5) 就是对 a 数组的切片操作，返回一个包含 a(3)、a(4) 和 a(5) 的新数组。

多维数组切片 对于二维数组 b，获取第 2 行的所有元素：

integer, dimension(3, 4) :: b = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12], [3, 4])
integer, dimension(4) :: sub_b
sub_b = b(2, :)

这里 b(2, :) 表示选取 b 数组第 2 行的所有列元素。同样，b(:, 3) 表示选取所有行的第 3 列元素。

步长切片

在切片操作中，我们还可以指定步长。例如，对于一维数组 a，获取从第 1 个元素开始，每隔 2 个元素的子数组：

integer, dimension(10) :: a = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
integer, dimension(5) :: sub_a
sub_a = a(1:10:2)

这里 a(1:10:2) 中，第三个参数 2 就是步长，表示每隔一个元素取一个值。

对于二维数组 b，获取第 1 行，从第 2 列开始，每隔 1 列的元素：

integer, dimension(3, 4) :: b = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12], [3, 4])
integer, dimension(2) :: sub_b
sub_b = b(1, 2:4:1)

数组运算

数组与标量的运算

在 Fortran 中，数组可以与标量进行各种算术运算。例如，将数组 a 的每个元素都乘以 2：

integer, dimension(5) :: a = [1, 2, 3, 4, 5]
a = a * 2

这里 a * 2 会将数组 a 的每个元素与标量 2 相乘，并将结果赋回给 a。同样，a + 3、a - 1、a / 2 等运算也是合法的。

数组与数组的运算

对应元素运算 两个形状相同的数组可以进行对应元素的算术运算。例如，将数组 a 和 b 的对应元素相加：

integer, dimension(5) :: a = [1, 2, 3, 4, 5]
integer, dimension(5) :: b = [5, 4, 3, 2, 1]
integer, dimension(5) :: c
c = a + b

这里 c 数组的元素 c(i) 等于 a(i) + b(i)，i = 1, 2, 3, 4, 5。

矩阵运算（二维数组） 对于二维数组（可视为矩阵），Fortran 支持矩阵乘法运算。例如，矩阵 A 和 B 相乘：

real, dimension(2, 3) :: A = reshape([1.0, 2.0, 3.0, 4.0, 5.0, 6.0], [2, 3])
real, dimension(3, 2) :: B = reshape([7.0, 8.0, 9.0, 10.0, 11.0, 12.0], [3, 2])
real, dimension(2, 2) :: C
C = matmul(A, B)

这里 matmul 函数用于矩阵乘法，A 是 2×3 的矩阵，B 是 3×2 的矩阵，结果 C 是 2×2 的矩阵。

数组聚合函数

基本聚合函数

求和函数 sum 计算数组所有元素的总和。例如，对于数组 a：

integer, dimension(5) :: a = [1, 2, 3, 4, 5]
integer :: total
total = sum(a)

这里 total 的值为 1 + 2 + 3 + 4 + 5 = 15。

求最大值函数 maxval 和求最小值函数 minval 获取数组中的最大值和最小值。例如：

integer, dimension(5) :: a = [1, 2, 3, 4, 5]
integer :: max_num, min_num
max_num = maxval(a)
min_num = minval(a)

这里 max_num 为 5，min_num 为 1。

按维度聚合

对于多维数组，我们可以按特定维度进行聚合。例如，对于二维数组 b，按行求和：

integer, dimension(3, 4) :: b = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12], [3, 4])
integer, dimension(3) :: row_sum
row_sum = sum(b, dim = 2)

这里 sum(b, dim = 2) 表示对 b 数组按第二维度（列）进行求和，结果 row_sum 是一个长度为 3 的一维数组，分别对应 b 数组每一行的元素和。

同样，sum(b, dim = 1) 表示按第一维度（行）进行求和，结果是一个长度为 4 的一维数组，对应 b 数组每一列的元素和。

数组重塑与转置

数组重塑

reshape 函数不仅可以用于初始化数组，还可以在程序运行过程中改变数组的形状。例如，将一个一维数组 a 重塑为二维数组：

integer, dimension(6) :: a = [1, 2, 3, 4, 5, 6]
integer, dimension(2, 3) :: b
b = reshape(a, [2, 3])

这里将长度为 6 的一维数组 a 重塑为 2×3 的二维数组 b。

数组转置

对于二维数组（矩阵），可以使用 transpose 函数进行转置。例如，对于矩阵 A：

real, dimension(2, 3) :: A = reshape([1.0, 2.0, 3.0, 4.0, 5.0, 6.0], [2, 3])
real, dimension(3, 2) :: At
At = transpose(A)

这里 At 就是 A 的转置矩阵，A 是 2×3 的矩阵，转置后 At 是 3×2 的矩阵。

数组掩码操作

掩码数组的概念

掩码数组是一种根据掩码（一个布尔数组）来选择或操作数组元素的技术。例如，我们有一个数组 a 和一个掩码数组 mask：

integer, dimension(5) :: a = [1, 2, 3, 4, 5]
logical, dimension(5) :: mask = [.true., .false., .true., .false., .true.]

掩码数组操作

选择元素 我们可以根据掩码选择 a 数组中的元素。例如，将掩码为 .true. 的元素提取出来：

integer, dimension(count(mask)) :: selected_a
selected_a = a(mask)

这里 count(mask) 用于计算掩码数组中 .true. 的个数，a(mask) 就是根据掩码选择的 a 数组中的元素组成的新数组。

条件赋值 根据掩码对数组元素进行条件赋值。例如，将掩码为 .true. 的 a 数组元素乘以 2：

a(mask) = a(mask) * 2

这样，a 数组中对应掩码为 .true. 的元素就会被乘以 2，而掩码为 .false. 的元素保持不变。

数组的内存管理与性能优化

数组存储布局

Fortran 数组在内存中是按列优先存储的。以二维数组 b 为例，b(1, 1)、b(2, 1)、b(3, 1) 等列元素在内存中是连续存储的。这种存储布局对于某些数值计算算法有重要影响，例如矩阵乘法，如果按照列优先的方式访问数组元素，可以提高缓存命中率，从而提升性能。

动态数组与内存分配

在 Fortran 中，可以使用动态数组来根据程序运行时的需求分配内存。例如：

integer, dimension(:), allocatable :: a
integer :: n
n = 10
allocate(a(n))
a = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
deallocate(a)

这里先声明一个可分配的一维数组 a，然后根据变量 n 的值动态分配内存，使用完后再释放内存。合理使用动态数组可以避免内存浪费，提高程序的灵活性。

性能优化技巧

减少数组边界检查 在编译时使用相关编译选项（如 -O2 等优化选项）可以减少数组边界检查，从而提高程序运行速度。但需要注意的是，关闭边界检查可能会导致程序在访问越界数组元素时出现未定义行为，所以在确保数组访问安全的前提下使用该方法。
向量化操作 尽可能使用数组运算代替循环操作，因为现代 Fortran 编译器能够对数组运算进行向量化优化，利用 CPU 的向量指令集提高运算速度。例如，将数组每个元素乘以 2，使用数组运算 a = a * 2 通常比使用循环逐个元素赋值要快。

案例分析：数值计算中的数组应用

线性方程组求解

考虑一个简单的线性方程组 Ax = b，其中 A 是系数矩阵，x 是未知数向量，b 是常数向量。我们可以使用 Fortran 数组来实现高斯消元法求解该方程组。

program gauss_elimination
    implicit none
    integer, parameter :: n = 3
    real, dimension(n, n) :: A = reshape([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 10.0], [n, n])
    real, dimension(n) :: b = [1.0, 2.0, 3.0]
    real, dimension(n) :: x
    integer :: i, j, k
    real :: pivot

    do k = 1, n - 1
        pivot = A(k, k)
        do i = k + 1, n
            do j = k, n
                A(i, j) = A(i, j) - A(k, j) * A(i, k) / pivot
            end do
            b(i) = b(i) - b(k) * A(i, k) / pivot
        end do
    end do

    x(n) = b(n) / A(n, n)
    do i = n - 1, 1, -1
        x(i) = (b(i) - dot_product(A(i, i + 1:n), x(i + 1:n))) / A(i, i)
    end do

    write(*, *) 'Solution vector x:'
    do i = 1, n
        write(*, *) x(i)
    end do
end program gauss_elimination

在这个程序中，我们使用二维数组 A 表示系数矩阵，一维数组 b 和 x 分别表示常数向量和未知数向量。通过高斯消元法对 A 和 b 进行操作，最终求解出 x。

数值积分

以梯形积分法为例，计算函数 f(x) = x^2 在区间 [0, 1] 上的积分。

program trapezoidal_integration
    implicit none
    real, parameter :: a = 0.0, b = 1.0
    integer, parameter :: n = 100
    real :: h, integral
    real, dimension(n + 1) :: x
    real, dimension(n + 1) :: f
    integer :: i

    h = (b - a) / n
    do i = 0, n
        x(i + 1) = a + i * h
        f(i + 1) = x(i + 1) ** 2
    end do

    integral = 0.5 * h * (f(1) + f(n + 1))
    do i = 2, n
        integral = integral + h * f(i)
    end do

    write(*, *) 'Integral value:', integral
end program trapezoidal_integration

这里我们使用一维数组 x 存储积分区间的离散点，f 数组存储对应点的函数值。通过梯形积分公式计算积分值，数组操作在这个过程中起到了关键作用，方便地存储和处理离散数据。

与其他语言的数组交互

Fortran 与 C 的数组交互

在一些混合编程场景中，可能需要 Fortran 与 C 语言进行数组交互。Fortran 数组按列优先存储，而 C 数组按行优先存储，这是需要注意的关键差异。

Fortran 调用 C 函数传递数组 假设 C 语言中有一个函数 add_arrays 用于将两个数组对应元素相加：

#include <stdio.h>

void add_arrays(double *a, double *b, double *c, int n) {
    for (int i = 0; i < n; i++) {
        c[i] = a[i] + b[i];
    }
}

在 Fortran 中调用该函数：

program call_c_function
    implicit none
    interface
        subroutine add_arrays(a, b, c, n) bind(C, name = "add_arrays")
            use iso_c_binding
            real(c_double), intent(in) :: a(*)
            real(c_double), intent(in) :: b(*)
            real(c_double), intent(out) :: c(*)
            integer(c_int), intent(in) :: n
        end subroutine add_arrays
    end interface

    real(8), dimension(5) :: a = [1.0d0, 2.0d0, 3.0d0, 4.0d0, 5.0d0]
    real(8), dimension(5) :: b = [5.0d0, 4.0d0, 3.0d0, 2.0d0, 1.0d0]
    real(8), dimension(5) :: c
    integer :: n

    n = 5
    call add_arrays(a, b, c, n)

    write(*, *) 'Result array c:'
    do i = 1, 5
        write(*, *) c(i)
    end do
end program call_c_function

这里通过 bind(C) 声明与 C 函数的接口，并注意数据类型的对应（real(c_double) 对应 C 语言的 double，integer(c_int) 对应 C 语言的 int）。

C 调用 Fortran 函数传递数组 假设 Fortran 中有一个函数 multiply_array 用于将数组每个元素乘以一个标量：

subroutine multiply_array(a, scalar, n) bind(C)
    use iso_c_binding
    real(c_double), intent(inout) :: a(*)
    real(c_double), intent(in) :: scalar
    integer(c_int), intent(in) :: n
    integer :: i
    do i = 1, n
        a(i) = a(i) * scalar
    end do
end subroutine multiply_array

在 C 语言中调用该函数：

#include <stdio.h>
#include <stdlib.h>

extern void multiply_array(double *a, double scalar, int n);

int main() {
    double a[] = {1.0, 2.0, 3.0, 4.0, 5.0};
    double scalar = 2.0;
    int n = 5;

    multiply_array(a, scalar, n);

    for (int i = 0; i < n; i++) {
        printf("%lf ", a[i]);
    }
    printf("\n");

    return 0;
}

在这个过程中，要注意 Fortran 和 C 语言数组存储顺序的差异，在实际应用中可能需要进行适当的转换。

Fortran 与 Python 的数组交互

在科学计算领域，Python 是常用的语言，有时也需要与 Fortran 进行数组交互。可以通过 f2py 工具来实现 Fortran 与 Python 的接口。

使用 f2py 生成 Python 接口 假设 Fortran 有一个简单的函数 sum_array 用于计算数组元素和：

function sum_array(a, n) result(sum_val)
    implicit none
    integer, intent(in) :: n
    real, intent(in) :: a(n)
    real :: sum_val
    integer :: i
    sum_val = 0.0
    do i = 1, n
        sum_val = sum_val + a(i)
    end do
end function sum_array

使用 f2py 生成 Python 接口：

f2py -c -m sum_module sum_array.f90

这会生成一个 Python 模块 sum_module，可以在 Python 中使用：

import sum_module

a = [1.0, 2.0, 3.0, 4.0, 5.0]
sum_val = sum_module.sum_array(a)
print("Sum of array:", sum_val)

通过 f2py，Fortran 数组可以方便地与 Python 数组（如 numpy 数组）进行交互，充分发挥两种语言的优势。

通过以上对 Fortran 高级数组技术的探讨，我们可以看到 Fortran 在数组处理方面具有强大的功能和灵活性，无论是在数值计算、科学模拟还是与其他语言的交互中，都能为开发者提供有效的工具和方法。