代码实现的是基于MPI的并行计算，代码如下：
#include <stdio.h> // needed for printing
#include <math.h> // needed for tanh, used in init function
#include "params.h" // model & simulation parameters
#include <mpi.h>
static int workerRank;
static int numWorkers;
static int rowsRange[501];
static MPI_Status status;
static int rangeBegin;
static int rangeEnd;
void init(double u[N][N], double v[N][N]){
double uhi, ulo, vhi, vlo;
uhi = 0.5; ulo = -0.5; vhi = 0.1; vlo = -0.1;
for (int i=rangeBegin; i < rangeEnd; i++){
for (int j=0; j < N; j++){
u[i][j] = ulo + (uhi-ulo)*0.5*(1.0 + tanh((i-N/2)/16.0));
v[i][j] = vlo + (vhi-vlo)*0.5*(1.0 + tanh((j-N/2)/16.0));
}
}
}
void dxdt(double du[N][N], double dv[N][N], double u[N][N], double v[N][N]){
double lapu, lapv;
int up, down, left, right;
for (int i = rangeBegin; i < rangeEnd; i++){
for (int j = 0; j < N; j++){
if (i == 0){
down = i;
}
else{
down = i-1;
}
if (i == N-1){
up = i;
}
else{
up = i+1;
}
if (j == 0){
left = j;
}
else{
left = j-1;
}
if (j == N-1){
right = j;
}
else{
right = j+1;
}
lapu = u[up][j] + u[down][j] + u[i][left] + u[i][right] + -4.0*u[i][j];
lapv = v[up][j] + v[down][j] + v[i][left] + v[i][right] + -4.0*v[i][j];
du[i][j] = DD*lapu + u[i][j]*(1.0 - u[i][j])*(u[i][j]-b) - v[i][j];
dv[i][j] = d*DD*lapv + c*(a*u[i][j] - v[i][j]);
}
}
}
void step(double du[N][N], double dv[N][N], double u[N][N], double v[N][N]){
for (int i = rangeBegin; i < rangeEnd; i++){
for (int j = 0; j < N; j++){
u[i][j] += dt*du[i][j];
v[i][j] += dt*dv[i][j];
}
}
}
double norm(double x[N][N]){
double nrmx = 0.0;
for (int i = rangeBegin; i < rangeEnd; i++){
for (int j = 0; j < N; j++){
nrmx += x[i][j]*x[i][j];
}
}
/*
集合通信函数
规约函数 MPI_Reduce()，将通信子内各进程的同一个变量参与规约计算，并向指定的进程输出计算结果
int MPI_Reduce(
void *input_data, //指向发送消息的内存块的指针
void *output_data, //指向接收（输出）消息的内存块的指针
int count，//数据量
MPI_Datatype datatype,//数据类型
MPI_Op operator,//规约操作
int dest，//要接收（输出）消息的进程的进程号
MPI_Comm comm);//通信器，指定通信范围
*/
double result;
MPI_Reduce(&nrmx, &result, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
return result;
}
// Calculate the range of rows for different processes
void calcRowsRange(){
int numIndividual = N / numWorkers;
int moreRows = N % numWorkers;
for(int idx = 0; idx < numWorkers; ++idx){
if(idx < moreRows){
rowsRange[idx + 1] = rowsRange[idx] + numIndividual + 1;
} else {
rowsRange[idx + 1] = rowsRange[idx] + numIndividual;
}
}
// Set the range of current worker
rangeBegin = rowsRange[workerRank];
rangeEnd = rowsRange[workerRank + 1];
}
// Swap local rows between adjacent worker
void swapRows(double u[N][N], double v[N][N]){
// Swap local rows with next worker
int nextWorkerRank = workerRank + 1;
int prevWorkerRank = workerRank - 1;
if(workerRank < numWorkers - 1){
double *addr = &u[rowsRange[nextWorkerRank]][0];
int n = (rowsRange[nextWorkerRank + 1] - rowsRange[nextWorkerRank]) * N;
/*
接收操作，在收到匹配消息之前不会返回
addr:指向包含要发送的数据的缓冲区的指针
n:缓冲区中元素的数目。 如果消息的数据部分为空，请将 count 参数设置为 0
MPI_DOUBLE:缓冲区数组中元素的数据类型
nextWorkerRank:指定通信器内发送进程的排名。 指定 MPI_ANY_SOURCE 常量，以指定任何源都是可接受的
标签(0):用于区分不同类型的消息的消息标记。 指定 MPI_ANY_TAG 常量以指示任何标记都是可接受的
MPI_COMM_WORLD:通信器的句柄
status:返回时，包含指向 MPI_Status 结构的指针，其中存储了有关接收的消息的信息
*/
MPI_Recv(addr, n, MPI_DOUBLE, nextWorkerRank, 0, MPI_COMM_WORLD, &status);
addr = &u[rowsRange[workerRank]][0];
n = (rowsRange[workerRank + 1] - rowsRange[workerRank]) * N;
/*
执行标准模式发送操作，并在可以安全重用发送缓冲区时返回
addr:指向包含要发送的数据的缓冲区的指针
n:缓冲区中元素的数目。 如果消息的数据部分为空，请将 count 参数设置为 0
MPI_DOUBLE:缓冲区数组中元素的数据类型
nextWorkerRank:指定通信器内发送进程的排名
标签(0):用于区分不同类型的消息的消息标记
MPI_COMM_WORLD:通信器的句柄
*/
MPI_Send(addr, n, MPI_DOUBLE, nextWorkerRank, 0, MPI_COMM_WORLD);
}
// Swap local rows with previous worker
if(workerRank > 0){
double *addr = &u[rowsRange[workerRank]][0];
int n = (rowsRange[workerRank + 1] - rowsRange[workerRank]) * N;
MPI_Send(addr, n, MPI_DOUBLE, prevWorkerRank, 0, MPI_COMM_WORLD);
addr = u[rowsRange[prevWorkerRank]];
n = (rowsRange[prevWorkerRank + 1] - rowsRange[prevWorkerRank]) * N;
MPI_Recv(addr, n, MPI_DOUBLE, prevWorkerRank, 0, MPI_COMM_WORLD, &status);
}
}
int main(int argc, char** argv){
double t = 0.0, nrmu, nrmv;
double u[N][N], v[N][N], du[N][N], dv[N][N];
// Initialize the MPI program
/*
每一个被MPI程序调用的第一个MPI函数都是MPI_Init,
要在调用任何MPI函数之前调用
*/
MPI_Init(&argc, &argv);
/*
当MPI初始化后，每一个活动进程变成了一个叫MPI_COMM_WORLD的通信域中的成员。
一个通信域是一个不透明的对象，提供了在进程之间传递消息的环境。
在一个通信域内的进程是有序的，一个进程的序号便是它在整个排序中的位置。
*/
MPI_Comm_rank(MPI_COMM_WORLD, &workerRank);
/*
进程可以通过调用函数MPI_Comm_size来确定一个通信域中的进程总数
*/
MPI_Comm_size(MPI_COMM_WORLD, &numWorkers);
// Calculate the range of rows for different processes
calcRowsRange();
FILE *fptr = fopen("part2.dat", "w");
fprintf(fptr, "#t\t\tnrmu\t\tnrmv\n");
// initialize the state
init(u, v);
// Swap local rows between adjacent worker
swapRows(u, v);
// time-loop
for (int k=0; k < M; k++){
// track the time
t = dt*k;
// evaluate the PDE
dxdt(du, dv, u, v);
// update the state variables u,v
step(du, dv, u, v);
// Swap local rows between adjacent worker
swapRows(u, v);
if (k%m == 0){
// calculate the norms
nrmu = norm(u);
nrmv = norm(v);
if(workerRank == 0){
printf("t = %2.1f\tu-norm = %2.5f\tv-norm = %2.5f\n", t, nrmu, nrmv);
fprintf(fptr, "%f\t%f\t%f\n", t, nrmu, nrmv);
}
}
}
if(workerRank == 0){
fclose(fptr);
}
/*
并行结束
*/
MPI_Finalize();
return 0;
}
)
