Jacobi迭代並行演算法

    Jacobi迭代是一種常見的迭代方法，迭代得到的新值是原來舊值點相鄰資料點的平均。串列程式片段如下：

    並行化方法之一，可以考慮按列劃分，邊界點新值的計算需要相鄰邊界其它塊的資料，所以在劃分後，每一個資料區塊的兩邊各增加一列，用於存放通訊得到的資料。如：

program mainimplicit noneinclude 'mpif.h'! define the size of array, 100 x 100integer, parameter::totalsize = 100, steps = 50real(kind=8) time_begin, time_endinteger n, myid, numprocs, i, j, rcinteger dp(0:totalsize), deltareal a(totalsize, totalsize), b(totalsize, totalsize)integer begin_col, end_col, ierrinteger left, right, tag1, tag2, taginteger status(MPI_STATUS_SIZE),REQ(16)call MPI_INIT(ierr)call MPI_COMM_RANK(MPI_COMM_WORLD, myid, ierr)call MPI_COMM_SIZE(MPI_COMM_WORLD, numprocs, ierr)!print *, 'Process', myid, 'of', numprocs, 'is alive'if(myid.eq.numprocs-1) thentime_begin = MPI_WTIME(ierr)endif! boundary values initialized as 8.0, others as 0do j = 1, totalsizedo i = 1, totalsizea(i, j) = 0.0enddoenddodo i = 1, totalsizea(1,i) = 8.0a(totalsize, i) = 8.0a(i, 1) = 8.0a(i, totalsize) = 8.0enddo!partition the taskdelta =ceiling(dfloat(totalsize) / numprocs)dp(0) = 2do i = 1, numprocsdp(i) = dp(i-1) + deltaenddodp(numprocs) = totalsizebegin_col = dp(myid)end_col = dp(myid+1) - 1!print the computing domain of each processcall MPI_Barrier(MPI_COMM_WORLD,IERR)!call sleep((myid+1)*5)write(*,*) '[',begin_col,'-----------------------',end_col,']'tag1 = 3tag2 = 4!define left process and right processif(myid.gt.0) thenleft = myid -1elseleft = MPI_PROC_NULLendifif(myid.lt.(numprocs-1)) thenright = myid + 1elseright = MPI_PROC_NULLendif!Jacobi iteration begins, number of  iterations: stepsdo n = 1, steps! data exchange before iteration begins!shift from left to rightcall MPI_SENDRECV(a(1, end_col), totalsize, MPI_REAL,     &right, tag1,a(1,begin_col-1),totalsize, MPI_REAL, left, tag1,     &MPI_COMM_WORLD, status, ierr)!shift from right to leftcall MPI_SENDRECV(a(1,begin_col), totalsize, MPI_REAL,     &  left, tag2,a(1,end_col+1),totalsize, MPI_REAL, right, tag2,     &  MPI_COMM_WORLD, status, ierr)do j = begin_col, end_coldo i = 2, totalsize-1b(i,j) = (a(i,j+1)+a(i,j-1)+a(i+1,j)+a(i-1,j))*0.25enddoenddodo j = begin_col, end_col        do i = 2, totalsize-1a(i,j) = b(i,j)enddoenddoenddoif(myid.eq.numprocs-1) thentime_end = MPI_WTIME(ierr)write(*,*)time_end - time_beginendif!each process print the data in its computing domain!do i = 2, totalsize-1!print *,myid,(a(i,j),j=begin_col, end_col)!enddo!the last process (process numprocs-1) collects data. process 0 ~ numprocs-2 sends data to !process  numprocs-1.!call MPI_Barrier(MPI_COMM_WORLD,IERR)if(myid.ne.(numprocs-1)) thentag = myid +10call MPI_Isend(a(1,begin_col),totalsize*delta,     &MPI_DOUBLE_PRECISION ,numprocs-1,tag,     &MPI_COMM_WORLD,REQ(myid+1),IERR)endifif(myid.eq.(numprocs-1)) thendo i = 0, numprocs-2tag = i + 10call MPI_recv(a(1,dp(i)),totalsize*delta,MPI_DOUBLE_PRECISION,     &i,tag,MPI_COMM_WORLD,REQ(i+1),IERR)enddoendifcall MPI_FINALIZE(rc)end

這是運行結果：

[root@c0108 ~]# mpiexec -n  5 ./jacobi [           2 -----------------------          21 ] [          22 -----------------------          41 ] [          42 -----------------------          61 ] [          82 -----------------------          99 ]  9.968996047973633E-003 [          62 -----------------------          81 ][root@c0108 ~]# mpiexec -n 10 ./jacobi [           2 -----------------------          11 ] [          12 -----------------------          21 ] [          22 -----------------------          31 ] [          42 -----------------------          51 ] [          32 -----------------------          41 ] [          82 -----------------------          91 ] [          72 -----------------------          81 ] [          92 -----------------------          99 ]  1.573109626770020E-002 [          62 -----------------------          71 ] [          52 -----------------------          61 ][root@c0108 ~]# mpiexec -n 3 ./jacobi [           2 -----------------------          35 ] [          70 -----------------------          99 ]  1.249790191650391E-002 [          36 -----------------------          69 ][root@c0108 ~]# ls

    上述代碼中，使用了虛擬進程。虛擬進程(MPI_PROC_NULL)是不存在的假想進程，在MPI中的主要作用是充當真實進程通訊的目的或源，引入虛擬進程的目的是在某些情況下，編寫通訊語句的方便。當一個真實進程向一個虛擬進程發送資料或從一個虛擬進程接收資料時，該真實進程會立即正確返回，如同執行了一個空操作。引入虛擬進程不僅可以大大簡化處理邊界的代碼，而且使程式顯得簡潔易懂。在捆綁發送接收操作中經常這種通訊手段。