1. 项目概述为什么数组是Fortran的“王牌”如果你刚开始接触Fortran可能会觉得这门语言有点“老古董”。但当你真正用它来处理科学计算、数值模拟或者大规模数据处理时你会立刻发现它的魅力所在。而数组就是Fortran这门“老将”手中最锋利、最趁手的武器。我最初接触Fortran是为了处理一个气候模型的海量网格数据当时用其他语言写的脚本跑起来慢得像蜗牛内存占用还高得吓人。后来切换到Fortran仅仅是利用了其原生的、高效的数组操作性能就提升了不止一个数量级。这让我深刻体会到在Fortran的世界里不懂数组就等于没入门。简单来说数组就是一组具有相同类型的数据元素的集合通过一个名字和下标来访问。这听起来和C、Python里的数组或列表差不多对吧但Fortran的数组远不止于此。它从语言设计之初就是为了高效处理大型数值数组而生的。这种“基因”层面的优化体现在内存布局、向量化操作、内置函数等方方面面。无论是处理一个简单的向量还是一个庞大的四维时空数据场比如温度、压力、速度在三维空间加一维时间上的分布Fortran都能提供极其简洁和高效的语法支持。所以这篇教程的目标很明确不仅要让你知道Fortran数组怎么声明、怎么用更要让你理解它为什么这么设计以及如何利用这些特性写出既快又稳的代码。无论你是要用Fortran做有限元分析比如配合Abaqus、计算流体力学还是处理实验数据数组都是你必须跨过去的第一道也是最重要的一道坎。2. 数组的核心概念与声明从一维到多维2.1 数组的声明与维度在Fortran中声明数组的核心是明确三要素数组名、数据类型和形状维度与范围。最经典的声明方式是在变量名后加上用括号括起来的维度信息。! 声明一个长度为10的一维实型数组 real :: temperature(10) ! 声明一个3行4列的二维整型数组 integer :: matrix(3, 4) ! 声明一个三维数组例如表示一个空间体x, y, z上的物理量 real :: density(100, 100, 50)这里需要注意下标范围的默认约定。在早期的Fortran如Fortran 77中数组下标默认从1开始。这是Fortran的一个传统也符合很多数学和物理公式的书写习惯第一个元素是第1个。在现代FortranFortran 90及以后中你可以显式地指定下标范围这带来了更大的灵活性。! 显式声明下标从0开始到9结束共10个元素 real :: vector(0:9) ! 下标从-5到5共11个元素常用于对称区间处理 integer :: symmetric_index(-5:5) ! 二维数组行从1到3列从0到2 real :: grid(1:3, 0:2)为什么下标从1开始是传统这很大程度上源于其科学计算背景。在数学中向量和矩阵的第一个元素通常标记为1如a₁, M₁₁。Fortran遵循这一惯例使得代码公式和数学表达式之间的转换更为直观减少了因下标偏移如C语言的0起始导致的“差一错误”off-by-one error。2.2 动态数组运行时决定大小静态数组的大小在编译时就必须确定。但在实际科研中我们经常需要处理大小在程序运行时才能确定的数据比如从文件读取的未知行数数据集。这时就需要动态数组Allocatable Arrays。声明动态数组使用allocatable属性并在声明时用冒号:来标记未指定的维度。数组的实际大小通过allocate语句在运行时分配使用完毕后应使用deallocate语句释放内存这是一个好习惯。program dynamic_array_demo implicit none real, allocatable :: data(:) ! 一维动态数组 real, allocatable :: image(:, :) ! 二维动态数组 integer :: n, m, i, stat ! 从用户输入或文件读取数组大小 print *, 请输入一维数组的长度 read *, n ! 分配内存并检查是否成功 allocate(data(n), statstat) if (stat / 0) then print *, 内存分配失败 stop end if ! 使用数组... data 3.14 ! 分配二维数组 m 5 allocate(image(n, m)) ! ... 使用 image 数组进行计算 ! 计算结束后释放内存 deallocate(data, image) end program dynamic_array_demo实操心得务必检查分配状态。allocate语句的stat参数至关重要特别是在处理可能非常大的数组时。内存不足是运行时常见错误通过检查stat程序可以优雅地失败并给出提示而不是直接崩溃。养成allocate后必跟if (stat/0)检查的习惯。2.3 数组的整体操作与片段操作这是Fortran数组语法中最优雅、最高效的部分之一。你不需要写循环就能对整个数组或数组的一部分进行操作。整体赋值与运算real :: A(10, 10), B(10, 10), C(10, 10) ! 将数组A所有元素赋值为1.0 A 1.0 ! 数组对应元素相加结果存入C。注意形状必须完全一致 C A B ! 数组与标量运算所有元素乘以2.0 B 2.0 * A ! 内置函数作用在整个数组上返回一个同形状的数组 C sin(A) ! 计算A中每个元素的正弦值数组片段Array Sections你可以使用冒号:操作符来选取数组的一个连续子集也就是片段。integer :: vec(100) integer :: mat(10, 10) ! 选取第5到第15个元素 vec(5:15) 0 ! 选取所有行用冒号:表示第3列 mat(:, 3) 1 ! 将整个第三列赋值为1 ! 选取第2到第5行第4到第7列 mat(2:5, 4:7) 9 ! 使用步长选取奇数索引的元素 vec(1:100:2) -1 ! 第1,3,5,...,99个元素被赋值为-1为什么这很强大这种语法不仅让代码更简洁、更接近数学表达比如C A B更重要的是它给了编译器极大的优化空间。编译器可以识别出这是对连续内存块的操作从而更容易利用现代CPU的SIMD单指令多数据流指令进行向量化大幅提升计算性能。你手写一个循环编译器可能还需要分析才能决定是否向量化而这种整体操作几乎是直接告诉编译器“这里可以向量化”3. 数组在内存中的布局与高效访问理解数组在内存中是如何存储的是写出高性能Fortran代码的关键。Fortran使用“列优先”存储顺序。对于一个二维数组A(m, n)它在内存中是按列连续存放的先存储第一列的所有行元素A(1,1), A(2,1), ..., A(m,1)然后是第二列A(1,2), A(2,2), ..., A(m,2)依此类推。! 假设有数组 A(3, 2) ! 内存布局顺序是 ! A(1,1) - A(2,1) - A(3,1) - A(1,2) - A(2,2) - A(3,2)这与C/C、Python (NumPy默认) 的“行优先”顺序正好相反。这个差异在混合编程如用C调用Fortran库时至关重要。对性能的影响访问局部性现代CPU从内存中读取数据时并不是一次只拿一个数而是会一次性抓取一个“缓存行”通常64字节的数据到高速缓存中。如果程序访问的数据在内存中是连续的那么缓存命中率就高速度飞快如果程序跳跃式访问就会导致大量的“缓存未命中”CPU不得不等待慢速的内存性能急剧下降。因此在Fortran中遍历数组时为了获得最佳性能内层循环应该是第一个下标行索引。! 高效的遍历方式内层循环变行i real :: arr(1000, 1000) integer :: i, j do j 1, 1000 ! 外层循环列 do i 1, 1000 ! 内层循环行 arr(i, j) arr(i, j) * 2.0 end do end do ! 低效的遍历方式内层循环变列j do i 1, 1000 ! 外层循环行 do j 1, 1000 ! 内层循环列 arr(i, j) arr(i, j) * 2.0 ! 每次访问都跳过了大量内存 end do end do在低效的方式中当i1时我们访问arr(1,1)下一个访问arr(1,2)。但在内存中arr(1,2)离arr(1,1)很远中间隔着第1列的其他999个元素这破坏了访问的局部性。注意这是一个非常经典的性能陷阱。很多从C/C转过来的程序员会不自觉地写出低效的循环顺序。记住口诀“Fortran列优先内层循环变行索引”。4. 强大的内置数组函数Fortran提供了一系列内置函数来处理数组这些函数很多都能自动进行向量化并行计算用起来非常方便。4.1 形状查询与操作函数size(array [, dim]): 返回数组的总元素数或指定维度dim上的大小。integer :: A(5, 10) print *, size(A) ! 输出 50 print *, size(A, dim1) ! 输出 5 (第一维大小) print *, size(A, dim2) ! 输出 10 (第二维大小)shape(array): 返回一个一维数组包含各维度的大小。print *, shape(A) ! 输出 [5, 10]lbound(array [, dim])/ubound(array [, dim]): 返回数组在指定维度dim上的下标下界和上界。real :: B(0:4, -1:8) print *, lbound(B, dim1), ubound(B, dim1) ! 输出 0 4 print *, lbound(B, dim2), ubound(B, dim2) ! 输出 -1 84.2 数学统计函数这些函数通常可以带一个可选的dim参数。如果不指定dim则对整个数组操作如果指定了dim则沿着该维度操作结果数组的维度会减少一维。sum(array [, dim, mask]): 求和。real :: vals(3) [1.0, 2.0, 3.0] print *, sum(vals) ! 输出 6.0 real :: mat(2,2) reshape([1.,2.,3.,4.], [2,2]) ! 矩阵 [1.0 3.0] ! [2.0 4.0] print *, sum(mat, dim1) ! 按列求和输出 [3.0, 7.0] print *, sum(mat, dim2) ! 按行求和输出 [4.0, 6.0]product(array [, dim]): 求乘积。maxval(array [, dim, mask])/minval(array [, dim, mask]): 求最大值/最小值。print *, maxval(vals) ! 输出 3.0 print *, minval(mat) ! 输出 1.0maxloc(array [, mask])/minloc(array [, mask]): 返回最大值/最小值所在的位置下标。返回的是一个一维数组长度等于原数组的秩维度数。integer :: pos(2) pos maxloc(mat) ! 找出mat中最大值的位置 print *, pos ! 输出 [2, 2] (即第2行第2列值为4.0)dot_product(vector_a, vector_b): 计算两个一维数组向量的点积。real :: v1(3)[1.,2.,3.], v2(3)[4.,5.,6.] print *, dot_product(v1, v2) ! 输出 32.0 (1*42*53*6)matmul(matrix_a, matrix_b): 矩阵乘法。这是Fortran进行线性代数计算的核心函数之一高度优化。real :: A(2,3), B(3,2), C(2,2) A reshape([1,2,3,4,5,6], [2,3]) B reshape([7,8,9,10,11,12], [3,2]) C matmul(A, B) ! C 将是 A 和 B 的矩阵乘积实操心得善用mask参数。很多统计函数如sum,maxval支持mask参数这是一个逻辑型数组形状与原数组相同。只有对应mask为.true.的元素才会被纳入计算。这在处理带有无效值如NaN或需要条件统计的数据时非常有用。real :: data(100) logical :: valid(100) ! ... 给 data 赋值并设置 valid 掩码例如标记大于0的数据为有效 valid (data 0.0) print *, 有效数据的平均值, sum(data, maskvalid) / count(valid)5. 数组的输入输出与格式化将数组数据读写到文件或屏幕是常见操作。Fortran提供了灵活的输入输出方式。5.1 隐式Do循环输出这是一种简洁的输出数组所有元素的方式。integer :: arr(5) [1, 2, 3, 4, 5] ! 输出所有元素用空格分隔 print *, arr ! 输出 1 2 3 4 5 ! 更可控的格式化输出每行3个元素 write(*, (3I5)) arr ! 输出 ! 1 2 3 ! 4 5(3I5)是格式描述符表示每行输出3个整数每个整数占5个字符宽度。5.2 读写文件将数组写入文件或从文件读入同样方便。program io_demo implicit none real :: data(100, 50) integer :: i, j, unit ! 生成一些示例数据 do j 1, 50 do i 1, 100 data(i, j) sin(0.1*i) * cos(0.05*j) end do end do ! 打开文件进行写入 open(newunitunit, fileoutput_data.dat, statusreplace, actionwrite) ! 将整个二维数组写入文件按列优先顺序 write(unit, *) data close(unit) ! 从文件读回数据到另一个数组 real, allocatable :: data_in(:, :) open(newunitunit, fileoutput_data.dat, statusold, actionread) ! 动态分配数组大小需要先知道大小这里假设已知 allocate(data_in(100, 50)) read(unit, *) data_in close(unit) end program io_demo注意使用write(unit, *)这种列表导向格式输出时数据是以空格分隔的纯文本。读取时read(unit, *)会按顺序读入足够的数据填满data_in。这种方式简单但文件可能很大且不保留维度信息。对于需要保存维度或更复杂结构的数据可以考虑使用未格式化二进制文件或NetCDF、HDF5等科学数据格式。5.3 格式化控制对于需要精细控制输出格式的场景格式描述符是关键。real :: results(5) [3.14159, 1.41421, 2.71828, 1.61803, 0.57721] ! 以表格形式输出每个数占12个字符宽度保留6位小数 write(*, (5F12.6)) results ! 输出类似 ! 3.141590 1.414210 2.718280 1.618030 0.577210 ! 每个数单独一行带说明文字 do i 1, 5 write(*, (A, I1, A, F10.5)) Result(, i, ) , results(i) end do6. 数组实战一个简单的数值微分示例让我们用一个完整的例子把上面讲的知识点串起来。假设我们有一个由函数f(x) sin(x)在区间[0, 2π]上离散采样得到的一维数组我们要用中心差分法计算其导数的近似值。中心差分公式f(x_i) ≈ (f(x_{i1}) - f(x_{i-1})) / (2 * h)其中h是采样间距。program numerical_derivative implicit none integer, parameter :: n 100 ! 采样点数量 real, parameter :: pi 3.141592653589793 real, parameter :: x_start 0.0, x_end 2.0 * pi real :: x(n), f(n), dfdx_approx(n) real :: h integer :: i ! 1. 生成等间距网格点和函数值 h (x_end - x_start) / real(n - 1) do i 1, n x(i) x_start real(i-1) * h f(i) sin(x(i)) end do ! 2. 使用中心差分法计算内部点的导数 ! 注意边界点i1和in无法用中心差分这里简单处理 dfdx_approx(1) 0.0 ! 前向或后向差分此处赋0 dfdx_approx(n) 0.0 do i 2, n-1 dfdx_approx(i) (f(i1) - f(i-1)) / (2.0 * h) end do ! 3. 输出结果前10个点 print *, x sin(x) approx cos(x) exact cos(x) print *, -------------------------------------------------------- do i 1, 10 write(*, (4F12.6)) x(i), f(i), dfdx_approx(i), cos(x(i)) end do ! 4. 计算并输出整体误差L2范数 ! 忽略边界点计算内部点的误差均方根 real :: error, rms_error error 0.0 do i 2, n-1 error error (dfdx_approx(i) - cos(x(i)))**2 end do rms_error sqrt(error / real(n-2)) print *, print *, RMS error (internal points): , rms_error end program numerical_derivative这个例子涵盖了数组的声明、循环填充、基于下标的计算、格式化输出以及简单的误差分析。你可以尝试修改函数f(x)或者使用更复杂的差分格式。7. 常见问题与避坑指南在实际使用Fortran数组时下面这些“坑”我几乎都踩过希望你能避开。7.1 数组越界Index out of bounds这是最常见的运行时错误。访问了数组声明范围之外的下标。integer :: arr(10) arr(11) 5 ! 编译可能通过但运行时会崩溃或产生不可预知结果。排查与预防开启编译器边界检查。在编译时加上调试选项如gfortran的-fcheckboundsIntel Fortran的/check:bounds。这会在运行时捕获越界访问并给出明确的错误信息。仔细检查循环的起始和结束条件。特别是当循环变量依赖于其他计算时。使用lbound和ubound函数来获取数组边界而不是硬编码数字。do i lbound(arr, 1), ubound(arr, 1) ! 安全操作 end do7.2 形状不匹配Shape mismatch在进行数组整体赋值或运算时等号两边的数组形状必须一致。real :: A(5), B(10) A B ! 错误形状 (5) 和 (10) 不匹配。排查编译器通常能捕获这类错误。确保进行操作的数组维度、各维度大小完全相同。使用shape()函数打印数组形状进行调试。7.3 动态数组未分配或已释放访问一个尚未分配allocate或已经释放deallocate的动态数组会导致严重错误。real, allocatable :: temp(:) temp(1) 10.0 ! 错误temp尚未分配内存。 allocate(temp(100)) ! ... 使用 temp deallocate(temp) temp(50) 20.0 ! 错误temp已被释放。预防在使用前用allocated()函数检查数组是否已分配。if (.not. allocated(temp)) then allocate(temp(100)) end if释放数组后最好将其显式置为空nullify不适用于普通可分配数组但可以将其指向一个空状态或者更简单地在释放后避免再次使用。更好的做法是在子程序或作用域结束时才释放数组。7.4 性能陷阱低效的循环与临时数组除了前面提到的内存访问顺序问题过度使用临时数组也会拖慢程序。! 低效创建了临时数组保存中间结果 real, allocatable :: tmp(:) tmp A B C tmp * D ! 高效合并操作减少内存分配和拷贝 C (A B) * D现代Fortran编译器很智能像(AB)*D这样的表达式编译器通常会生成优化的代码避免生成完整的临时数组AB而是采用融合循环的方式计算。但如果你显式地分配了tmp就强制了中间结果的存储。7.5 与C语言交互时的顺序问题当你用C语言调用Fortran子程序或者反过来时必须特别注意数组的存储顺序。假设Fortran有一个子程序! Fortran side: sub.f90 subroutine process_matrix(mat, m, n) integer, intent(in) :: m, n real, intent(inout) :: mat(m, n) ! 列优先 ! ... 处理矩阵 end subroutine process_matrix在C语言中调用时传入的数组必须是转置后的或者C代码需要按照列优先的顺序来理解和填充数据。/* C side: main.c */ extern void process_matrix_(float *mat, int *m, int *n); /* 注意名称可能加了下划线 */ int main() { int m3, n2; float mat[3][2]; /* 在C中是行优先: mat[行][列] */ /* 填充数据... */ /* 但Fortran期望的是列优先所以这里直接传递可能导致数据解释错误 */ process_matrix_(mat[0][0], m, n); return 0; }解决方案要么在C端将数据按列优先存储要么在接口层进行转置。对于复杂的交互建议使用ISO_C_BINDING模块它提供了更安全、标准的互操作方式。数组是Fortran的灵魂掌握它你就掌握了用Fortran进行高效数值计算的钥匙。从简单的声明到高效的内存访问从强大的内置函数到实际的应用案例每一步都需要理解和练习。刚开始可能会觉得有些概念如下标顺序、片段操作比较别扭但一旦习惯你就会爱上这种直接、高效的表达方式。我个人的体会是在处理真正的多维网格数据时Fortran数组的简洁性和性能是其他语言难以比拟的。多写多试多调优尤其是多用编译器提供的诊断和性能分析工具你会进步飞快。
