原文链接
两小时入门MPI与并行计算(五):对等模式(实现Jacobi迭代的并行计算) - 知乎
因为原文只有fortran,而且评论链接里的代码没有考虑c++是row major的储存格式,在这里提供对应的C++代码。
原文中对应的2D case
在这里为了和原文对应,运行时只可以使用“mysize”个processor来运算,即当mysize=4时,运行
mpirun -n 4 ./jacobi
否则会提示printf("Please use %d processors!\n", mysize);
除此之外,为了更好理解mpi initialization的过程,我添加了 calcu_jacobi的变量,如果将其comment掉,就不会进行jacobi iteration的运算,而是直接输出initialization的结果
/* ************************************************************************
> File Name: jacobi.cpp
> Author: Roy
> Created Time: Tuesday, October 17, 2023 AM11:27:18
> Description:
************************************************************************/
#include "mpi.h"
#include <iostream>
#include <iomanip>
#include <stdio.h>
#define calcu_jacobi // A flag whether calculate jacobi iiterations or not
int main(int argc, char *argv[])
{
int steps = 10; //Number of iterations
const int Acolsize = 16; //Number of colomns
const int mysize = 4; //Number of processors
const int Arowsize = mysize + 2; //2 additional rows for buffer
double A[Arowsize][Acolsize] = {0};// processors' data matrix
double B[Arowsize][Acolsize] = {0};// processors' buffer matrix
int begin_row = 1, end_row = Arowsize - 2; // Indices for begin col and end col
int myid = 0, numprocs = 0; //MPI rank and total numprocessors, but here we have already set it 4
int up = 0, down = 0; //rank ID for up and down processors
int tag1 = 3, tag2 = 4; //MPI send and recv tags
MPI_Status status;
MPI_Init(&argc, &argv);
MPI_Comm_rank(
MPI_COMM_WORLD,
&myid
);
MPI_Comm_size(
MPI_COMM_WORLD,
&numprocs
);
if (numprocs!=mysize){
printf("Please use %d processors!\n", mysize);
MPI_Finalize();
return 0;
}
printf("Process %d of %d is alive\n", myid, numprocs); // print MPI proc information
// assign values to A so that the gathered matrix are all 0 except for boundaries (assign to 8.0)
for (int i = 0; i < Arowsize; i++){
for (int j = 0; j < Acolsize; j++){
A[i][j] = 0.0;
}
A[i][0] = 8.0;
A[i][Acolsize-1] = 8.0;
}
// Assume rank 0 is the uppermost part and 3 is the lowest part of the gathered matrix
// proc0 additional assignment
if (myid == 0){
for (int j = 0; j < Acolsize; j++){
A[1][j] = 8.0;
}
}
// proc3 additional assignment
if (myid == 3){
for (int j = 0; j < Acolsize; j++){
A[Arowsize-2][j] = 8.0;
}
}
#ifdef calcu_jacobi
// Calculate nearby ranks
// 0 is the uppermost and 3 is the lowest
if (myid > 0){
up = myid - 1;
}
else
up = MPI_PROC_NULL;
if (myid < 3){
down = myid + 1;
}
else
down = MPI_PROC_NULL;
// Jacobi iterations
for (int i = 0; i < steps; i++){
// up to down transfor
MPI_Sendrecv(
&A[Arowsize - 2][0], // onst void *sendbuf; lowest data row (excluding the buffer row)
Acolsize, // int sendcount; length of data sended
MPI_DOUBLE, // MPI_Datatype;
down, // int dest; rank of receiver;
tag1, // int sending,
//----------------------------------
&A[0][0], // onst void *recvbuf; uppermost buffer row
Acolsize, // int recvcount; length of data received
MPI_DOUBLE, // MPI_Datatype;
up, // int source; rank of sender;
tag1, // int sendtag,
MPI_COMM_WORLD, // MPI_Comm comm;
&status
);
// down to up transfor
MPI_Sendrecv(
&A[1][0], // onst void *sendbuf; uppermost data row (excluding the buffer row)
Acolsize, // int sendcount; length of data sended
MPI_DOUBLE, // MPI_Datatype;
up, // int dest; rank of receiver;
tag2, // int sendtag,
//----------------------------------
&A[Arowsize - 1][0], // onst void *recvbuf; lowest buffer row
Acolsize, // int recvcount; length of data received
MPI_DOUBLE, // MPI_Datatype;
down, // int source; rank of sender;
tag2, // int sending,
MPI_COMM_WORLD, // MPI_Comm comm;
&status
);
}
if (myid == 0){
begin_row = 2; // For 0, not modify the upper boundary
}
if (myid == 3){
end_row = Arowsize - 3; // For 3, not modify the lower boundary
}
// Iterate from the starting row to end rwo, and excluding the buffer colomns
for (int i = begin_row; i <= end_row; i++){
for (int j = 1; j <= Acolsize - 2; j++){
B[i][j] = (A[i][j + 1] + A[i][j - 1] + A[i + 1][j] + A[i - 1][j]) * 0.25;
}
}
// Update A from B
for (int i = begin_row; i <= end_row; i++){
for (int j = 1; j <= Acolsize - 2; j++){
A[i][j] = B[i][j];
}
}
#endif
// Print result
for (int row = 1; row <= Arowsize - 2; row++){
std::cout << "proc (" << myid << "): ";
for(int col = 0; col <= Acolsize - 1; col++){
std::cout << std::setiosflags(std::ios_base::left)
<< std::setw(2) <<A[row][col]
<< std::resetiosflags(std::ios_base::left);
}
std::cout << std::endl;
}
MPI_Finalize();
}
拓展一下:3D jacobi iterator
和上面代码设置一样,唯一不同是这个是计算3D case
/* ************************************************************************
> File Name: 3Djacobi.cpp
> Author: Roy
> Created Time: Thursday, October 19, 2023 PM04:25:03
> Description: mpi jacobi calculator for 3D matrix operations
************************************************************************/
#include "mpi.h"
#include <iostream>
#include <iomanip>
#include <stdio.h>
#include <cmath>
#define calcu_jacobi // A flag whether calculate jacobi iiterations or not
int main(int argc, char *argv[])
{
int steps = 1; //Number of iterations
// We want a 3D matrix of size P*R*C, or more precisely [mysize*4][5][16], whose index=(p*R+r)*C+c
// Two additional planes should be allocated for buffer purpose, so we need to have A as A[mysize+2][5][16] with 4 copies opposed by 4 processors
const int Acolsize = 16; //Number of colomns
const int mysize = 4; //Number of processors
const int Arowsize = 5; //Number of rows
const int Aplanesize = mysize + 2; //2 additional planes for buffer
double A[Aplanesize][Arowsize][Acolsize] = {0};// processors' data matrix
double B[Aplanesize][Arowsize][Acolsize] = {0};// processors' buffer matrix
int begin_p = 1, end_p = Aplanesize - 2; // Indices for begin planes and end planes
int myid = 0, numprocs = 0; //MPI rank and total numprocessors, but here we have already set it 4
int up = 0, down = 0; //rank ID for up and down processors
int tag1 = 3, tag2 = 4; //MPI send and recv tags
MPI_Status status;
MPI_Init(&argc, &argv);
MPI_Comm_rank(
MPI_COMM_WORLD,
&myid
);
MPI_Comm_size(
MPI_COMM_WORLD,
&numprocs
);
if (numprocs!=mysize){
printf("Please use %d processors!\n", mysize);
MPI_Finalize();
return 0;
}
printf("Process %d of %d is alive\n", myid, numprocs); // print MPI proc information
// assign values to A so that the gathered matrix are all 0 except for boundaries (assign to 8.0)
for (int i = 0; i < Aplanesize; i++){
for (int j = 1; j < Arowsize - 1; j++){
for (int k = 1; k < Acolsize - 1; k++){
A[i][j][k] = 0.0;
}
}
for (int k = 0; k < Acolsize; k++){
A[i][0][k] = 8.0;
A[i][Arowsize - 1][k] = 8.0;
}
for (int j = 0; j < Arowsize; j++){
A[i][j][0] = 8.0;
A[i][j][Acolsize - 1] = 8.0;
}
}
// Assume rank 0 is the uppermost plane and 3 is the lowest plane of the gathered matrix
// proc0 additional assignment
if (myid == 0){
for (int j = 0; j < Arowsize; j++){
for (int k = 0; k < Acolsize; k++){
A[1][j][k] = 8.0;
}
}
}
// proc3 additional assignment
if (myid == 3){
for (int j = 0; j < Arowsize; j++){
for (int k = 0; k < Acolsize; k++){
A[Aplanesize - 2][j][k] = 8.0;
}
}
}
#ifdef calcu_jacobi
// Calculate nearby ranks
// 0 is the uppermost and 3 is the lowest
if (myid > 0){
up = myid - 1;
}
else
up = MPI_PROC_NULL;
if (myid < 3){
down = myid + 1;
}
else
down = MPI_PROC_NULL;
// Jacobi iterations
for (int i = 0; i < steps; i++){
// up to down transfor
MPI_Sendrecv(
&A[Aplanesize - 2][0][0], // onst void *sendbuf; lowest data plane (excluding the buffer plane)
Arowsize*Acolsize, // int sendcount; length of data sended
MPI_DOUBLE, // MPI_Datatype;
down, // int dest; rank of receiver;
tag1, // int sending,
//----------------------------------
&A[0][0][0], // onst void *recvbuf; uppermost buffer plane
Arowsize*Acolsize, // int recvcount; length of data received
MPI_DOUBLE, // MPI_Datatype;
up, // int source; rank of sender;
tag1, // int sendtag,
MPI_COMM_WORLD, // MPI_Comm comm;
&status
);
// down to up transfor
MPI_Sendrecv(
&A[1][0][0], // onst void *sendbuf; uppermost data row (excluding the buffer row)
Arowsize*Acolsize, // int sendcount; length of data sended
MPI_DOUBLE, // MPI_Datatype;
up, // int dest; rank of receiver;
tag2, // int sendtag,
//----------------------------------
&A[Aplanesize - 1][0][0], // onst void *recvbuf; lowest buffer plane
Arowsize*Acolsize, // int recvcount; length of data received
MPI_DOUBLE, // MPI_Datatype;
down, // int source; rank of sender;
tag2, // int sending,
MPI_COMM_WORLD, // MPI_Comm comm;
&status
);
}
if (myid == 0){
begin_p = 2; // For 0, not modify the upper boundary
}
if (myid == 3){
end_p = Aplanesize - 3; // For 3, not modify the lower boundary
}
// Iterate from the starting plane to end plane, and excluding the buffer elements
for (int i = begin_p; i <= end_p; i++){
for (int j = 1; j <= Arowsize - 2; j++){
for (int k = 1; k <= Acolsize - 2; k++){
B[i][j][k] = (A[i][j][k + 1] + A[i][j][k - 1] + A[i][j + 1][k] + A[i][j - 1][k] + A[i + 1][j][k] + A[i - 1][j][k]) / 6.0;
}
}
}
// Update A from B
for (int i = begin_p; i <= end_p; i++){
for (int j = 1; j <= Arowsize - 2; j++){
for (int k = 1; k <= Acolsize - 2; k++){
A[i][j][k] = B[i][j][k];
}
}
}
#endif
// Print result
int printid = 3;
if (myid == printid)
{
for (int plane = 1; plane <= Aplanesize - 2; plane++){
for(int row = 0; row <= Arowsize - 1; row++){
printf("proc (\" %d \") @ A(%d, %d, k) = ", myid, plane, row);
for(int col = 0; col <= Acolsize - 1; col++){
std::cout << std::setiosflags(std::ios_base::left)
<< std::setw(5) << std::ceil(A[plane][row][col] * 100.0) / 100.0 //rounding to 2 decimal places
<< std::resetiosflags(std::ios_base::left);
}
std::cout << std::endl;
}
}
}
MPI_Finalize();
}