1. 項目概述為什么數(shù)組是Fortran的“王牌”如果你剛開始接觸Fortran可能會覺得這門語言有點“老古董”。但當(dāng)你真正用它來處理科學(xué)計算、數(shù)值模擬或者大規(guī)模數(shù)據(jù)處理時你會立刻發(fā)現(xiàn)它的魅力所在。而數(shù)組就是Fortran這門“老將”手中最鋒利、最趁手的武器。我最初接觸Fortran是為了處理一個氣候模型的海量網(wǎng)格數(shù)據(jù)當(dāng)時用其他語言寫的腳本跑起來慢得像蝸牛內(nèi)存占用還高得嚇人。后來切換到Fortran僅僅是利用了其原生的、高效的數(shù)組操作性能就提升了不止一個數(shù)量級。這讓我深刻體會到在Fortran的世界里不懂?dāng)?shù)組就等于沒入門。簡單來說數(shù)組就是一組具有相同類型的數(shù)據(jù)元素的集合通過一個名字和下標(biāo)來訪問。這聽起來和C、Python里的數(shù)組或列表差不多對吧但Fortran的數(shù)組遠不止于此。它從語言設(shè)計之初就是為了高效處理大型數(shù)值數(shù)組而生的。這種“基因”層面的優(yōu)化體現(xiàn)在內(nèi)存布局、向量化操作、內(nèi)置函數(shù)等方方面面。無論是處理一個簡單的向量還是一個龐大的四維時空數(shù)據(jù)場比如溫度、壓力、速度在三維空間加一維時間上的分布Fortran都能提供極其簡潔和高效的語法支持。所以這篇教程的目標(biāo)很明確不僅要讓你知道Fortran數(shù)組怎么聲明、怎么用更要讓你理解它為什么這么設(shè)計以及如何利用這些特性寫出既快又穩(wěn)的代碼。無論你是要用Fortran做有限元分析比如配合Abaqus、計算流體力學(xué)還是處理實驗數(shù)據(jù)數(shù)組都是你必須跨過去的第一道也是最重要的一道坎。2. 數(shù)組的核心概念與聲明從一維到多維2.1 數(shù)組的聲明與維度在Fortran中聲明數(shù)組的核心是明確三要素數(shù)組名、數(shù)據(jù)類型和形狀維度與范圍。最經(jīng)典的聲明方式是在變量名后加上用括號括起來的維度信息。! 聲明一個長度為10的一維實型數(shù)組 real :: temperature(10) ! 聲明一個3行4列的二維整型數(shù)組 integer :: matrix(3, 4) ! 聲明一個三維數(shù)組例如表示一個空間體x, y, z上的物理量 real :: density(100, 100, 50)這里需要注意下標(biāo)范圍的默認(rèn)約定。在早期的Fortran如Fortran 77中數(shù)組下標(biāo)默認(rèn)從1開始。這是Fortran的一個傳統(tǒng)也符合很多數(shù)學(xué)和物理公式的書寫習(xí)慣第一個元素是第1個。在現(xiàn)代FortranFortran 90及以后中你可以顯式地指定下標(biāo)范圍這帶來了更大的靈活性。! 顯式聲明下標(biāo)從0開始到9結(jié)束共10個元素 real :: vector(0:9) ! 下標(biāo)從-5到5共11個元素常用于對稱區(qū)間處理 integer :: symmetric_index(-5:5) ! 二維數(shù)組行從1到3列從0到2 real :: grid(1:3, 0:2)為什么下標(biāo)從1開始是傳統(tǒng)這很大程度上源于其科學(xué)計算背景。在數(shù)學(xué)中向量和矩陣的第一個元素通常標(biāo)記為1如a?, M??。Fortran遵循這一慣例使得代碼公式和數(shù)學(xué)表達式之間的轉(zhuǎn)換更為直觀減少了因下標(biāo)偏移如C語言的0起始導(dǎo)致的“差一錯誤”off-by-one error。2.2 動態(tài)數(shù)組運行時決定大小靜態(tài)數(shù)組的大小在編譯時就必須確定。但在實際科研中我們經(jīng)常需要處理大小在程序運行時才能確定的數(shù)據(jù)比如從文件讀取的未知行數(shù)數(shù)據(jù)集。這時就需要動態(tài)數(shù)組Allocatable Arrays。聲明動態(tài)數(shù)組使用allocatable屬性并在聲明時用冒號:來標(biāo)記未指定的維度。數(shù)組的實際大小通過allocate語句在運行時分配使用完畢后應(yīng)使用deallocate語句釋放內(nèi)存這是一個好習(xí)慣。program dynamic_array_demo implicit none real, allocatable :: data(:) ! 一維動態(tài)數(shù)組 real, allocatable :: image(:, :) ! 二維動態(tài)數(shù)組 integer :: n, m, i, stat ! 從用戶輸入或文件讀取數(shù)組大小 print *, 請輸入一維數(shù)組的長度 read *, n ! 分配內(nèi)存并檢查是否成功 allocate(data(n), statstat) if (stat / 0) then print *, 內(nèi)存分配失敗 stop end if ! 使用數(shù)組... data 3.14 ! 分配二維數(shù)組 m 5 allocate(image(n, m)) ! ... 使用 image 數(shù)組進行計算 ! 計算結(jié)束后釋放內(nèi)存 deallocate(data, image) end program dynamic_array_demo實操心得務(wù)必檢查分配狀態(tài)。allocate語句的stat參數(shù)至關(guān)重要特別是在處理可能非常大的數(shù)組時。內(nèi)存不足是運行時常見錯誤通過檢查stat程序可以優(yōu)雅地失敗并給出提示而不是直接崩潰。養(yǎng)成allocate后必跟if (stat/0)檢查的習(xí)慣。2.3 數(shù)組的整體操作與片段操作這是Fortran數(shù)組語法中最優(yōu)雅、最高效的部分之一。你不需要寫循環(huán)就能對整個數(shù)組或數(shù)組的一部分進行操作。整體賦值與運算real :: A(10, 10), B(10, 10), C(10, 10) ! 將數(shù)組A所有元素賦值為1.0 A 1.0 ! 數(shù)組對應(yīng)元素相加結(jié)果存入C。注意形狀必須完全一致 C A B ! 數(shù)組與標(biāo)量運算所有元素乘以2.0 B 2.0 * A ! 內(nèi)置函數(shù)作用在整個數(shù)組上返回一個同形狀的數(shù)組 C sin(A) ! 計算A中每個元素的正弦值數(shù)組片段Array Sections你可以使用冒號:操作符來選取數(shù)組的一個連續(xù)子集也就是片段。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 ! 使用步長選取奇數(shù)索引的元素 vec(1:100:2) -1 ! 第1,3,5,...,99個元素被賦值為-1為什么這很強大這種語法不僅讓代碼更簡潔、更接近數(shù)學(xué)表達比如C A B更重要的是它給了編譯器極大的優(yōu)化空間。編譯器可以識別出這是對連續(xù)內(nèi)存塊的操作從而更容易利用現(xiàn)代CPU的SIMD單指令多數(shù)據(jù)流指令進行向量化大幅提升計算性能。你手寫一個循環(huán)編譯器可能還需要分析才能決定是否向量化而這種整體操作幾乎是直接告訴編譯器“這里可以向量化”3. 數(shù)組在內(nèi)存中的布局與高效訪問理解數(shù)組在內(nèi)存中是如何存儲的是寫出高性能Fortran代碼的關(guān)鍵。Fortran使用“列優(yōu)先”存儲順序。對于一個二維數(shù)組A(m, n)它在內(nèi)存中是按列連續(xù)存放的先存儲第一列的所有行元素A(1,1), A(2,1), ..., A(m,1)然后是第二列A(1,2), A(2,2), ..., A(m,2)依此類推。! 假設(shè)有數(shù)組 A(3, 2) ! 內(nèi)存布局順序是 ! A(1,1) - A(2,1) - A(3,1) - A(1,2) - A(2,2) - A(3,2)這與C/C、Python (NumPy默認(rèn)) 的“行優(yōu)先”順序正好相反。這個差異在混合編程如用C調(diào)用Fortran庫時至關(guān)重要。對性能的影響訪問局部性現(xiàn)代CPU從內(nèi)存中讀取數(shù)據(jù)時并不是一次只拿一個數(shù)而是會一次性抓取一個“緩存行”通常64字節(jié)的數(shù)據(jù)到高速緩存中。如果程序訪問的數(shù)據(jù)在內(nèi)存中是連續(xù)的那么緩存命中率就高速度飛快如果程序跳躍式訪問就會導(dǎo)致大量的“緩存未命中”CPU不得不等待慢速的內(nèi)存性能急劇下降。因此在Fortran中遍歷數(shù)組時為了獲得最佳性能內(nèi)層循環(huán)應(yīng)該是第一個下標(biāo)行索引。! 高效的遍歷方式內(nèi)層循環(huán)變行i real :: arr(1000, 1000) integer :: i, j do j 1, 1000 ! 外層循環(huán)列 do i 1, 1000 ! 內(nèi)層循環(huán)行 arr(i, j) arr(i, j) * 2.0 end do end do ! 低效的遍歷方式內(nèi)層循環(huán)變列j do i 1, 1000 ! 外層循環(huán)行 do j 1, 1000 ! 內(nèi)層循環(huán)列 arr(i, j) arr(i, j) * 2.0 ! 每次訪問都跳過了大量內(nèi)存 end do end do在低效的方式中當(dāng)i1時我們訪問arr(1,1)下一個訪問arr(1,2)。但在內(nèi)存中arr(1,2)離arr(1,1)很遠中間隔著第1列的其他999個元素這破壞了訪問的局部性。注意這是一個非常經(jīng)典的性能陷阱。很多從C/C轉(zhuǎn)過來的程序員會不自覺地寫出低效的循環(huán)順序。記住口訣“Fortran列優(yōu)先內(nèi)層循環(huán)變行索引”。4. 強大的內(nèi)置數(shù)組函數(shù)Fortran提供了一系列內(nèi)置函數(shù)來處理數(shù)組這些函數(shù)很多都能自動進行向量化并行計算用起來非常方便。4.1 形狀查詢與操作函數(shù)size(array [, dim]): 返回數(shù)組的總元素數(shù)或指定維度dim上的大小。integer :: A(5, 10) print *, size(A) ! 輸出 50 print *, size(A, dim1) ! 輸出 5 (第一維大小) print *, size(A, dim2) ! 輸出 10 (第二維大小)shape(array): 返回一個一維數(shù)組包含各維度的大小。print *, shape(A) ! 輸出 [5, 10]lbound(array [, dim])/ubound(array [, dim]): 返回數(shù)組在指定維度dim上的下標(biāo)下界和上界。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 數(shù)學(xué)統(tǒng)計函數(shù)這些函數(shù)通?？梢詭б粋€可選的dim參數(shù)。如果不指定dim則對整個數(shù)組操作如果指定了dim則沿著該維度操作結(jié)果數(shù)組的維度會減少一維。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]): 返回最大值/最小值所在的位置下標(biāo)。返回的是一個一維數(shù)組長度等于原數(shù)組的秩維度數(shù)。integer :: pos(2) pos maxloc(mat) ! 找出mat中最大值的位置 print *, pos ! 輸出 [2, 2] (即第2行第2列值為4.0)dot_product(vector_a, vector_b): 計算兩個一維數(shù)組向量的點積。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進行線性代數(shù)計算的核心函數(shù)之一高度優(yōu)化。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參數(shù)。很多統(tǒng)計函數(shù)如sum,maxval支持mask參數(shù)這是一個邏輯型數(shù)組形狀與原數(shù)組相同。只有對應(yīng)mask為.true.的元素才會被納入計算。這在處理帶有無效值如NaN或需要條件統(tǒng)計的數(shù)據(jù)時非常有用。real :: data(100) logical :: valid(100) ! ... 給 data 賦值并設(shè)置 valid 掩碼例如標(biāo)記大于0的數(shù)據(jù)為有效 valid (data 0.0) print *, 有效數(shù)據(jù)的平均值, sum(data, maskvalid) / count(valid)5. 數(shù)組的輸入輸出與格式化將數(shù)組數(shù)據(jù)讀寫到文件或屏幕是常見操作。Fortran提供了靈活的輸入輸出方式。5.1 隱式Do循環(huán)輸出這是一種簡潔的輸出數(shù)組所有元素的方式。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個整數(shù)每個整數(shù)占5個字符寬度。5.2 讀寫文件將數(shù)組寫入文件或從文件讀入同樣方便。program io_demo implicit none real :: data(100, 50) integer :: i, j, unit ! 生成一些示例數(shù)據(jù) 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) ! 將整個二維數(shù)組寫入文件按列優(yōu)先順序 write(unit, *) data close(unit) ! 從文件讀回數(shù)據(jù)到另一個數(shù)組 real, allocatable :: data_in(:, :) open(newunitunit, fileoutput_data.dat, statusold, actionread) ! 動態(tài)分配數(shù)組大小需要先知道大小這里假設(shè)已知 allocate(data_in(100, 50)) read(unit, *) data_in close(unit) end program io_demo注意使用write(unit, *)這種列表導(dǎo)向格式輸出時數(shù)據(jù)是以空格分隔的純文本。讀取時read(unit, *)會按順序讀入足夠的數(shù)據(jù)填滿data_in。這種方式簡單但文件可能很大且不保留維度信息。對于需要保存維度或更復(fù)雜結(jié)構(gòu)的數(shù)據(jù)可以考慮使用未格式化二進制文件或NetCDF、HDF5等科學(xué)數(shù)據(jù)格式。5.3 格式化控制對于需要精細控制輸出格式的場景格式描述符是關(guān)鍵。real :: results(5) [3.14159, 1.41421, 2.71828, 1.61803, 0.57721] ! 以表格形式輸出每個數(shù)占12個字符寬度保留6位小數(shù) write(*, (5F12.6)) results ! 輸出類似 ! 3.141590 1.414210 2.718280 1.618030 0.577210 ! 每個數(shù)單獨一行帶說明文字 do i 1, 5 write(*, (A, I1, A, F10.5)) Result(, i, ) , results(i) end do6. 數(shù)組實戰(zhàn)一個簡單的數(shù)值微分示例讓我們用一個完整的例子把上面講的知識點串起來。假設(shè)我們有一個由函數(shù)f(x) sin(x)在區(qū)間[0, 2π]上離散采樣得到的一維數(shù)組我們要用中心差分法計算其導(dǎo)數(shù)的近似值。中心差分公式f(x_i) ≈ (f(x_{i1}) - f(x_{i-1})) / (2 * h)其中h是采樣間距。program numerical_derivative implicit none integer, parameter :: n 100 ! 采樣點數(shù)量 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. 生成等間距網(wǎng)格點和函數(shù)值 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. 使用中心差分法計算內(nèi)部點的導(dǎo)數(shù) ! 注意邊界點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. 輸出結(jié)果前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范數(shù) ! 忽略邊界點計算內(nèi)部點的誤差均方根 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這個例子涵蓋了數(shù)組的聲明、循環(huán)填充、基于下標(biāo)的計算、格式化輸出以及簡單的誤差分析。你可以嘗試修改函數(shù)f(x)或者使用更復(fù)雜的差分格式。7. 常見問題與避坑指南在實際使用Fortran數(shù)組時下面這些“坑”我?guī)缀醵疾冗^希望你能避開。7.1 數(shù)組越界Index out of bounds這是最常見的運行時錯誤。訪問了數(shù)組聲明范圍之外的下標(biāo)。integer :: arr(10) arr(11) 5 ! 編譯可能通過但運行時會崩潰或產(chǎn)生不可預(yù)知結(jié)果。排查與預(yù)防開啟編譯器邊界檢查。在編譯時加上調(diào)試選項如gfortran的-fcheckboundsIntel Fortran的/check:bounds。這會在運行時捕獲越界訪問并給出明確的錯誤信息。仔細檢查循環(huán)的起始和結(jié)束條件。特別是當(dāng)循環(huán)變量依賴于其他計算時。使用lbound和ubound函數(shù)來獲取數(shù)組邊界而不是硬編碼數(shù)字。do i lbound(arr, 1), ubound(arr, 1) ! 安全操作 end do7.2 形狀不匹配Shape mismatch在進行數(shù)組整體賦值或運算時等號兩邊的數(shù)組形狀必須一致。real :: A(5), B(10) A B ! 錯誤形狀 (5) 和 (10) 不匹配。排查編譯器通常能捕獲這類錯誤。確保進行操作的數(shù)組維度、各維度大小完全相同。使用shape()函數(shù)打印數(shù)組形狀進行調(diào)試。7.3 動態(tài)數(shù)組未分配或已釋放訪問一個尚未分配allocate或已經(jīng)釋放deallocate的動態(tài)數(shù)組會導(dǎo)致嚴(yán)重錯誤。real, allocatable :: temp(:) temp(1) 10.0 ! 錯誤temp尚未分配內(nèi)存。 allocate(temp(100)) ! ... 使用 temp deallocate(temp) temp(50) 20.0 ! 錯誤temp已被釋放。預(yù)防在使用前用allocated()函數(shù)檢查數(shù)組是否已分配。if (.not. allocated(temp)) then allocate(temp(100)) end if釋放數(shù)組后最好將其顯式置為空nullify不適用于普通可分配數(shù)組但可以將其指向一個空狀態(tài)或者更簡單地在釋放后避免再次使用。更好的做法是在子程序或作用域結(jié)束時才釋放數(shù)組。7.4 性能陷阱低效的循環(huán)與臨時數(shù)組除了前面提到的內(nèi)存訪問順序問題過度使用臨時數(shù)組也會拖慢程序。! 低效創(chuàng)建了臨時數(shù)組保存中間結(jié)果 real, allocatable :: tmp(:) tmp A B C tmp * D ! 高效合并操作減少內(nèi)存分配和拷貝 C (A B) * D現(xiàn)代Fortran編譯器很智能像(AB)*D這樣的表達式編譯器通常會生成優(yōu)化的代碼避免生成完整的臨時數(shù)組AB而是采用融合循環(huán)的方式計算。但如果你顯式地分配了tmp就強制了中間結(jié)果的存儲。7.5 與C語言交互時的順序問題當(dāng)你用C語言調(diào)用Fortran子程序或者反過來時必須特別注意數(shù)組的存儲順序。假設(shè)Fortran有一個子程序! Fortran side: sub.f90 subroutine process_matrix(mat, m, n) integer, intent(in) :: m, n real, intent(inout) :: mat(m, n) ! 列優(yōu)先 ! ... 處理矩陣 end subroutine process_matrix在C語言中調(diào)用時傳入的數(shù)組必須是轉(zhuǎn)置后的或者C代碼需要按照列優(yōu)先的順序來理解和填充數(shù)據(jù)。/* C side: main.c */ extern void process_matrix_(float *mat, int *m, int *n); /* 注意名稱可能加了下劃線 */ int main() { int m3, n2; float mat[3][2]; /* 在C中是行優(yōu)先: mat[行][列] */ /* 填充數(shù)據(jù)... */ /* 但Fortran期望的是列優(yōu)先所以這里直接傳遞可能導(dǎo)致數(shù)據(jù)解釋錯誤 */ process_matrix_(mat[0][0], m, n); return 0; }解決方案要么在C端將數(shù)據(jù)按列優(yōu)先存儲要么在接口層進行轉(zhuǎn)置。對于復(fù)雜的交互建議使用ISO_C_BINDING模塊它提供了更安全、標(biāo)準(zhǔn)的互操作方式。數(shù)組是Fortran的靈魂掌握它你就掌握了用Fortran進行高效數(shù)值計算的鑰匙。從簡單的聲明到高效的內(nèi)存訪問從強大的內(nèi)置函數(shù)到實際的應(yīng)用案例每一步都需要理解和練習(xí)。剛開始可能會覺得有些概念如下標(biāo)順序、片段操作比較別扭但一旦習(xí)慣你就會愛上這種直接、高效的表達方式。我個人的體會是在處理真正的多維網(wǎng)格數(shù)據(jù)時Fortran數(shù)組的簡潔性和性能是其他語言難以比擬的。多寫多試多調(diào)優(yōu)尤其是多用編譯器提供的診斷和性能分析工具你會進步飛快。