跳轉到內容

並行譜數值方法/Fortran 程式和 Windows 入門

Add links
來自華夏公益教科書，開放書籍，開放世界
< 並行譜數值方法
此頁面可能需要稽核質量。

Fortran 程式和 Windows 入門

[編輯 | 編輯原始碼]

示例程式

[編輯 | 編輯原始碼]

要使用 OpenMP 或 MPI（訊息傳遞介面）進行並行程式設計，我們通常需要使用比 Matlab 更低階的語言，例如 Fortran。另一種可能的語言選擇是 C，但是 Fortran 與 C 相比具有更強大的陣列處理能力，並且語法類似於 Matlab，因此通常更容易用於對規則陣列進行大量使用的科學計算。因此，在我們開始學習如何建立並行程式之前，在 Fortran 中介紹一些簡單的程式是有用的。關於 Fortran 的一個好的最新參考是 Metcalf、Reid 和 Cohen[1]。我們認識到，大多數人會不熟悉 Fortran，可能更熟悉 Matlab[2]、C 或 C++，但我們預計示例程式碼將使任何具有入門程式設計背景的人更容易使用。最近的一本指南描述瞭如何編寫高效的並行 Fortran 程式碼，該指南是 Levesque 和 Wagenbreth[3]。我們的程式被編寫為在密歇根大學的 Flux 叢集上執行。有關此叢集的更多資訊，請訪問 http://cac.engin.umich.edu/resources/systems/flux/ 和 http://cac.engin.umich.edu/started/index.html。以下是執行此程式所需的四個檔案。

  1. 一個 Makefile 用於在 Flux 上編譯 Fortran 程式碼，如列表 A 中所示。這應該儲存為 makefile。在使用 Makefile 編譯程式碼之前，您需要輸入
    module load fftw/3.2.1-intel
    在登入到 Flux 後的命令列提示符處。然後將 Makefile 和 heat.f90 放置在同一個目錄中，以下示例檔案假設此目錄為
    {${HOME}/ParallelMethods/Heat}
    並輸入
    make
    以編譯檔案。檔案編譯完成後，輸入
    qsub fluxsubscript
    以讓叢集執行您的程式，然後輸出結果。以下程式使用 FFTW 庫來執行快速傅立葉變換。有關此庫的更多資訊，請訪問 http://www.fftw.org/。

     

     

     

     
    ( A)
    用於編譯傅立葉譜 Fortran 熱方程程式的示例 Makefile 程式碼下載
    #define the complier
COMPILER = mpif90
# compilation settings, optimization, precision, parallelization
FLAGS = -O0

# libraries
LIBS = -L${FFTW_LINK} -lfftw3 -lm
# source list for main program
SOURCES = heat.f90

test: $(SOURCES)
${COMPILER} -o heat $(FLAGS) $(SOURCES) $(LIBS)

clean:
rm *.o

  2. 列表 B 中的 Fortran 程式 - 這應該儲存為 heat.f90

     

     

     

     
    ( B)
    使用向後尤拉時間步進法求解熱方程的 Fortran 傅立葉譜程式 程式碼下載
    !--------------------------------------------------------------------
!
!
! PURPOSE
!
! This program solves heat equation in 1 dimension
! u_t=\alpha*u_{xx}
! using a the backward Euler method for x\in[0,2\pi]
!
! The boundary conditions are u(0)=u(2\pi)
! The initial condition is u=sin(x)
!
! .. Parameters ..
! Nx = number of modes in x - power of 2 for FFT
! Nt = number of timesteps to take
! Tmax = maximum simulation time
! plotgap = number of timesteps between plots
! FFTW_IN_PLACE = value for FFTW input
! FFTW_MEASURE = value for FFTW input
! FFTW_EXHAUSTIVE = value for FFTW input
! FFTW_PATIENT = value for FFTW input
! FFTW_ESTIMATE = value for FFTW input
! FFTW_FORWARD = value for FFTW input
! FFTW_BACKWARD = value for FFTW input
! pi = 3.14159265358979323846264338327950288419716939937510d0
! L = width of box
! alpha = heat conductivity
! .. Scalars ..
! i = loop counter in x direction
! n = loop counter for timesteps direction
! allocatestatus = error indicator during allocation
! start = variable to record start time of program
! finish = variable to record end time of program
! count_rate = variable for clock count rate
! planfx = Forward 1d fft plan in x
! planbx = Backward 1d fft plan in x
! dt = timestep
! .. Arrays ..
! u = approximate REAL solution
! v = Fourier transform of approximate solution
! vna = temporary field
! .. Vectors ..
! kx = fourier frequencies in x direction
! x = x locations
! time = times at which save data
! name_config = array to store filename for data to be saved
!
! REFERENCES
!
! ACKNOWLEDGEMENTS
!
! ACCURACY
!
! ERROR INDICATORS AND WARNINGS
!
! FURTHER COMMENTS
! Check that the initial iterate is consistent with the
! boundary conditions for the domain specified
!--------------------------------------------------------------------
! External routines required
!
! External libraries required
! FFTW3 -- Fast Fourier Transform in the West Library
! (http://www.fftw.org/)

PROGRAM main

! Declare variables
IMPLICIT NONE
INTEGER(kind=4),	PARAMETER :: Nx=64
INTEGER(kind=4),	PARAMETER	:: Nt=20
REAL(kind=8),	PARAMETER	&
:: pi=3.14159265358979323846264338327950288419716939937510d0
REAL(kind=8),	PARAMETER	:: L=5.0d0	
REAL(kind=8),	PARAMETER	:: alpha=0.50d0	
REAL(kind=8)	:: dt=0.2d0/REAL(Nt,kind(0d0))	
COMPLEX(KIND=8), DIMENSION(:),ALLOCATABLE	:: kx	
REAL(kind=8), DIMENSION(:),ALLOCATABLE	:: x	
COMPLEX(KIND=8), DIMENSION(:,:),ALLOCATABLE	:: u,v	
REAL(kind=8), DIMENSION(:),ALLOCATABLE	:: time
COMPLEX(KIND=8), DIMENSION(:),ALLOCATABLE	:: vna	
INTEGER(kind=4)	:: i,j,k,n
INTEGER(kind=4)	:: start, finish, count_rate, AllocateStatus
INTEGER(kind=4), PARAMETER	:: FFTW_IN_PLACE = 8, FFTW_MEASURE = 0, &
FFTW_EXHAUSTIVE = 8, FFTW_PATIENT = 32, FFTW_ESTIMATE = 64
INTEGER(kind=4), PARAMETER :: FFTW_FORWARD = -1, FFTW_BACKWARD=1	
COMPLEX(KIND=8), DIMENSION(:),ALLOCATABLE :: fftfx,fftbx
INTEGER(kind=8)	:: planfx,planbx
CHARACTER*100	:: name_config

CALL system_clock(start,count_rate)	
ALLOCATE(kx(1:Nx),x(1:Nx),u(1:Nx,1:1+Nt),v(1:Nx,1:1+Nt),&
    time(1:1+Nt),vna(1:Nx),fftfx(1:Nx),fftbx(1:Nx),&
    stat=AllocateStatus)
    IF (AllocateStatus .ne. 0) STOP

! set up ffts
CALL dfftw_plan_dft_1d(planfx,Nx,fftfx(1:Nx),fftbx(1:Nx),&
FFTW_FORWARD,FFTW_ESTIMATE)
CALL dfftw_plan_dft_1d(planbx,Nx,fftbx(1:Nx),fftfx(1:Nx),&
FFTW_BACKWARD,FFTW_ESTIMATE)

PRINT *,'Setup FFTs'

! setup fourier frequencies
DO i=1,1+Nx/2
kx(i)= cmplx(0.0d0,1.0d0)*REAL(i-1,kind(0d0))/L
END DO
kx(1+Nx/2)=0.00d0
DO i = 1,Nx/2 -1
kx(i+1+Nx/2)=-kx(1-i+Nx/2)
END DO
DO i=1,Nx
x(i)=(-1.00d0 + 2.00d0*REAL(i-1,kind(0d0))/REAL(Nx,KIND(0d0)))*pi*L
END DO
PRINT *,'Setup grid and fourier frequencies and splitting coefficients'

u(1:Nx,1)=sin(x(1:Nx))
! transform initial data
CALL dfftw_execute_dft_(planfx,u(1:Nx,1),v(1:Nx,1))
PRINT *,'Got initial data, starting timestepping'
time(1)=0.0d0

vna(1:Nx)=v(1:Nx,1)
PRINT *,'Starting timestepping'
DO n=1,Nt
DO i=1,Nx
vna(i)=vna(i)/(1-dt*kx(i)*kx(i))
END DO
PRINT *,'storing plot data ',n
time(n+1)=time(n)+dt
v(1:Nx,n+1)=vna(1:Nx)
CALL dfftw_execute_dft_(planbx,v(1:Nx,n+1),u(1:Nx,n+1))
u(1:Nx,n+1)=u(1:Nx,n+1)/REAL(Nx,KIND(0d0))	! normalize
END DO
PRINT *,'Finished time stepping'
CALL system_clock(finish,count_rate)
PRINT*,'Program took ',REAL(finish-start)/REAL(count_rate),'for execution'

! Write data out to disk

name_config = 'u.dat'
OPEN(unit=11,FILE=name_config,status="UNKNOWN")
REWIND(11)
DO j=1,1+Nt
DO i=1,Nx
WRITE(11,*) REAL(u(i,j))
END DO
END DO
CLOSE(11)

name_config = 'tdata.dat'
OPEN(unit=11,FILE=name_config,status="UNKNOWN")
REWIND(11)
DO j=1,1+Nt
WRITE(11,*) time(j)
END DO
CLOSE(11)

name_config = 'xcoord.dat'
OPEN(unit=11,FILE=name_config,status="UNKNOWN")
REWIND(11)
DO i=1,Nx
WRITE(11,*) x(i)
END DO
CLOSE(11)

PRINT *,'Saved data'
DEALLOCATE(kx,x,u,v,&
    time,vna,fftfx,fftbx,&
    stat=AllocateStatus)
    IF (AllocateStatus .ne. 0) STOP
CALL dfftw_destroy_plan(planbx)
CALL dfftw_destroy_plan(planfx)
CALL dfftw_cleanup()
PRINT *,'Program execution complete'
END PROGRAM main


  3. 列表 C 中用於在叢集上使用的示例提交指令碼 - 這應該儲存為 fluxsubscript。更多示例可以在 http://cac.engin.umich.edu/resources/software/pbs.html 找到。要使用它，請將電子郵件地址從 your_uniqname@umich.edu 更改為您可以接收作業開始和完成通知的電子郵件地址。

     

     

     

     
    ( C)
    在 Flux 上使用的示例提交指令碼 程式碼下載
    #!/bin/bash
#PBS -N heatequation
#PBS -l nodes=1,walltime=00:10:00
#PBS -l qos=math471f11_flux
#PBS -A math471f11_flux
#PBS -q flux
#PBS -M your_uniqname@umich.edu
#PBS -m abe
#PBS -V
# Create a local directory to run and copy your files to local.
# Let PBS handle your output
mkdir /tmp/${PBS_JOBID}
cp ${HOME}/ParallelMethods/Heat/heatequation /tmp/${PBS_JOBID}/heatequation
cd /tmp/${PBS_JOBID}
./heatequation
#Clean up your files
cd
cd ParallelMethods/Heat
# Retrieve your output
cp /tmp/${PBS_JOBID}/u.dat ${HOME}/ParallelMethods/Heat/u.dat
cp /tmp/${PBS_JOBID}/xcoord.dat ${HOME}/ParallelMethods/Heat/xcoord.dat
cp /tmp/${PBS_JOBID}/tdata.dat ${HOME}/ParallelMethods/Heat/tdata.dat

/bin/rm -rf /tmp/${PBS_JOBID}

  4. 一個 Matlab 繪圖指令碼[4] 用於生成圖 i 在列表 D 中。

     

     

     

     
    ( D)
    一個 Matlab 程式，用於繪製計算結果 程式碼下載
    % A Matlab program to plot the computed results

clear all; format compact, format short,
set(0,'defaultaxesfontsize',18,'defaultaxeslinewidth',.9,...
    'defaultlinelinewidth',3.5,'defaultpatchlinewidth',5.5);

% Load data
load('./u.dat');
load('./tdata.dat');
load('./xcoord.dat');
Tsteps = length(tdata);

Nx = length(xcoord); Nt = length(tdata);
 
u = reshape(u,Nx,Nt);

% Plot data
figure(3); clf; mesh(tdata,xcoord,u); xlabel t; ylabel x; zlabel('u');


     

     

     

     
    ( i)
    由 Fortran 程式在 B 中計算，並由 Matlab 程式在 D 中後處理的熱方程解。

Windows 入門

[編輯 | 編輯原始碼]

要在 Windows 上執行這些程式，您需要安裝 Cygwin。Cygwin 在 Windows 上模擬 bash shell。您還需要為 Cygwin 安裝以下軟體包：make、mpi、fftw3（列表不完整）。以下是修改後的 Makefile，其他檔案保持不變。

一個 Makefile 用於編譯列表 E 中的 Fortran 程式碼。這應該儲存為 makefile 而不是 ( A )。然後將 Makefile 和 heat.f90 ( B ) 放置在同一個目錄中。輸入

make

to compile the file. Once the file is compiled type

./heat.exe

to get your computer to run the program.

使用上面的 Matlab 繪圖指令碼獲取圖片。

 

 

 

 
( E)
用於編譯傅立葉譜 Fortran 熱方程程式的示例 Makefile 程式碼下載
#define the complier
COMPILER = mpif90
# compilation settings, optimization, precision, parallelization
FLAGS = -O0

# libraries
LIBS = -lfftw3 -lm
# source list for main program
SOURCES = heat.f90

test: $(SOURCES)
                ${COMPILER} -o heat $(FLAGS) $(SOURCES) $(LIBS)

clean:
       rm *.o 
       rm heat

練習

[編輯 | 編輯原始碼]
  1. 請閱讀網頁 http://cac.engin.umich.edu/started/index.html 上的資源，瞭解如何使用 Flux 叢集。
  2. 修改一維熱方程的 Fortran 程式，以求解 Allen-Cahn 方程，並選擇時間步進方案。建立執行輸出的繪圖。在您的解決方案中包含原始碼和繪圖。
  3. 修改一維熱方程的 Fortran 程式，以求解二維熱方程，並選擇時間步進方案。您的程式應該在每個時間步儲存欄位，而不是將所有欄位放在一個大型陣列中。建立執行的初始狀態和最終狀態的繪圖。在您的解決方案中包含原始碼和繪圖。

註釋

[編輯 | 編輯原始碼]
  1. ↑ Metcalf、Reid 和 Cohen (2009)
  2. ↑ 雖然 Matlab 是用 C 編寫的，但它最初是用 Fortran 編寫的，因此風格與 Fortran 相似。
  3. ↑ Levesque 和 Wagenbreth (2009)
  4. ↑

    對於許多計算問題，人們可以使用比生成結果所需的計算能力低 10-100 倍的計算能力來視覺化結果，因此對於規模不大的問題，使用 Matlab 等高階語言來後處理資料要容易得多。

參考文獻

[編輯 | 編輯原始碼]

Levesque, J.; Wagenbreth, G. (2011). 高效能計算：程式設計和應用. CRC 出版社。

Metcalf, M.; Reid, J.; Cohen, M. (2011). 現代 Fortran 解釋. 牛津大學出版社。

檢索自 "https://wikibook.tw/wiki/Parallel_Spectral_Numerical_Methods/Fortran_Programs_and_Getting_Started_on_Windows"
華夏公益教科書