我正在尝试使用 Fortran、MPI 和 Jacobi 方法实现方程组的解。请告诉我这里的错误是什么以及如何修复代码以获得正确的答案
program jacobi_test
! Подключение MPI
use mpi
implicit none
integer :: nob
real :: sbr, en = 0
! p_t_s - старт, параллельное
! p_t_f - финиш, параллельное
real :: p_t_s, p_t_f
! Переменные цикла
integer :: i, j, iter
! Максимальное количество
! итераций в м. Якоби
integer :: maxIters = 100000
! Количество строк и
! столбцов матрицы
integer, parameter :: n = 10
! A - матрица коэффициентов
! B - правая часть системы
! XE - точное решение системы
! X - параллельное решение
real(kind=8) :: A(n, n), B(n), XE(n), X(n)
! Переменные для MPI
integer :: size, rank, ierr
! Определение "кусков" матрицы
! для MPI решения
real(kind=8), dimension(:, :), allocatable :: A_local
real(kind=8), dimension(:), allocatable :: S_local
real(kind=8), dimension(:), allocatable :: B_local
! Инициализация MPI
call MPI_Init(ierr)
call MPI_COMM_SIZE(MPI_COMM_WORLD, size, ierr)
call MPI_COMM_RANK(MPI_COMM_WORLD, rank, ierr)
if (rank == 0) then
! Сгенерируем с датчика
! случайных чисел матрицу A
call RANDOM_NUMBER(A)
! Сгегерируем с датчика
! случайных чисел матриц XE
call RANDOM_NUMBER(XE)
! Диагональное преобладание
do i = 1, n
sbr = 0
do j = 1, n
sbr = sbr + ABS(A(i, j))
end do
A(i, i) = sbr
end do
! Определим вектор B
B = MATMUL(A, XE)
end if
! Время начала работы MPI
call MPI_Barrier(MPI_COMM_WORLD, ierr)
p_t_s = MPI_Wtime()
nob = n / size
! Выделение памяти на куски MPI
ALLOCATE(B_local(nob))
ALLOCATE(S_local(nob))
ALLOCATE(A_local(nob, n))
! Отправляем матрицу A и вектор
! B по всем процессам в MPI
CALL MPI_Scatter(A, nob * n, MPI_REAL8, A_local, nob * n, MPI_REAL8, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Scatter(B, nob, MPI_REAL8, B_local, nob, MPI_REAL8, 0, MPI_COMM_WORLD, ierr)
!
! Прописываем алгоритм для
! работы с каждым кусочком
! матрицы тут
!
X = 0.0
print *, A_local, "rank=", rank
print *, ""
print *, B_local, "rank=", rank
print *, ""
print *, "nob=", nob
do iter = 1, maxIters
S_local = 0.0
do i = 1, nob
S_local(i) = B_local(i)
do j = 1, n
if (i /= j) then
S_local(i) = S_local(i) - A_local(i, j) * X(j + nob * rank)
end if
end do
S_local(i) = S_local(i) / A_local(i, i)
X(i + nob * rank) = S_local(i)
end do
! Обмен полученных
! значений среди всех
! парал. процессов
call MPI_Allgather(X, nob, MPI_REAL8, X, nob, MPI_REAL8, MPI_COMM_WORLD, ierr)
end do
!
! Завершение работы с матрицами
!
! Время завершения работы MPI
call MPI_Barrier(MPI_COMM_WORLD, ierr)
p_t_f = MPI_Wtime()
if (rank == 0) then
print *, ""
print *, "Matrix XE:"
do i = 1, n
print *, "x(", i, ") = ", XE(i)
end do
print *, ""
print *, "Parallel solution: (", p_t_f - p_t_s, "s):"
do i = 1, n
print *, "x(", i, ") = ", X(i)
end do
print *, ""
print *, "Parallel solution time: (", p_t_f - p_t_s, "s)"
print *, ""
print *, "Euclidean norm:"
do i = 1, n
en = en + (XE(i) - X(i)) ** 2
end do
print *, SQRT(en)
end if
! Завершение работы MPI
call MPI_Finalize(ierr)
end program jacobi_test
我正在尝试使用 MPI_SCATTER，但我不明白这个问题。我有 MPI_Bcast 解决方案，效果很好
! mpif90 makarov_lab_0_s.f90 makarov_lab_0_test.f90 -o makarov_lab_0 && mpirun -np 4 ./makarov_lab_0
program jacobi_test
! Подключение MPI
use mpi
implicit none
real :: sbr, en = 0
! p_t_s - старт, параллельное
! p_t_f - финиш, параллельное
real :: p_t_s, p_t_f
! Переменные цикла
integer :: i, j, iter
! Максимальное количество
! итераций в м. Якоби
integer :: maxIters = 100000
! Количество строк и
! столбцов матрицы
integer, parameter :: n = 100
! A - матрица коэффициентов
! B - правая часть системы
! XE - точное решение системы
! X - параллельное решение
real(kind=8) :: A(n, n), B(n), XE(n), X(n)
! Переменные для MPI
integer :: size, rank, ierr
! (1) sRow - индекс стартовой
! полосы для матрицы
! и вектора X
! (2) eRow - индекс последней
! полосы для матрицы
! и вектора X
! (3) lRng - строки, стобцы
integer :: sRow, eRow, lRng
! Определение "кусков" матрицы
! для MPI решения
real(kind=8), dimension(:, :), allocatable :: A_local
real(kind=8), dimension(:), allocatable :: S_local
! Инициализация MPI
call MPI_Init(ierr)
call MPI_COMM_SIZE(MPI_COMM_WORLD, size, ierr)
call MPI_COMM_RANK(MPI_COMM_WORLD, rank, ierr)
if (rank == 0) then
! Сгенерируем с датчика
! случайных чисел матрицу A
call RANDOM_NUMBER(A)
! Сгегерируем с датчика
! случайных чисел матриц XE
call RANDOM_NUMBER(XE)
! Диагональное преобладание
do i = 1, n
sbr = 0
do j = 1, n
sbr = sbr + ABS(A(i, j))
end do
A(i, i) = sbr
end do
! Определим вектор B
B = MATMUL(A, XE)
end if
! Время начала работы MPI
call MPI_Barrier(MPI_COMM_WORLD, ierr)
p_t_s = MPI_Wtime()
! (1) - строки, стобцы; индексы
! (2) - // -- //
! (3) - строки, стобцы; кол-во
sRow = 1 + rank * (n / size)
eRow = sRow + n / size - 1
lRng = eRow - sRow + 1
! Выделение памяти на куски MPI
allocate(S_local(lRng))
allocate(A_local(lRng, n))
! Отправляем матрицу A и вектор
! B по всем процессам в MPI
call MPI_Bcast(A, n * n, MPI_REAL8, 0, MPI_COMM_WORLD, ierr)
call MPI_Bcast(B, n, MPI_REAL8, 0, MPI_COMM_WORLD, ierr)
! Копирование частей матрицы A
! и вектора B в локальные перем
A_local = A(sRow:eRow, :)
! Вычисление методом Якоби
do iter = 1, maxIters
S_local = 0.0
do i = 1, lRng
S_local(i) = B(sRow + i - 1)
do j = 1, n
if (sRow + i - 1 /= j) then
S_local(i) = S_local(i) - A_local(i, j) * X(j)
end if
end do
S_local(i) = S_local(i) / A_local(i, sRow + i - 1)
end do
! Обмен полученных
! значений среди всех
! парал. процессов
call MPI_Allgather(S_local, lRng, MPI_REAL8, X, lRng, MPI_REAL8, MPI_COMM_WORLD, ierr)
end do
! Время завершения работы MPI
call MPI_Barrier(MPI_COMM_WORLD, ierr)
p_t_f = MPI_Wtime()
if (rank == 0) then
print *, ""
print *, "Matrix XE:"
do i = 1, n
print *, "x(", i, ") = ", XE(i)
end do
print *, ""
print *, "Parallel solution: (", p_t_f - p_t_s, "s):"
do i = 1, n
print *, "x(", i, ") = ", X(i)
end do
print *, ""
print *, "Parallel solution time: (", p_t_f - p_t_s, "s)"
print *, ""
print *, "Euclidean norm:"
do i = 1, n
en = en + (XE(i) - X(i)) ** 2
end do
print *, SQRT(en)
end if
! Завершение работы MPI
call MPI_Finalize(ierr)
end program jacobi_test
现在我得到了这样的输出
Matrix XE:
x( 1 ) = 0.89133570357258263
x( 2 ) = 0.63776384763443716
x( 3 ) = 5.9840042244841740E-002
x( 4 ) = 0.33879246967332022
x( 5 ) = 0.76695216356640361
x( 6 ) = 0.83136160974869977
x( 7 ) = 0.16494638062655809
x( 8 ) = 0.98882870791022592
x( 9 ) = 0.84931545759316152
x( 10 ) = 0.21269066751687060
Parallel solution: ( 0.103964999 s):
x( 1 ) = NaN
x( 2 ) = NaN
x( 3 ) = NaN
x( 4 ) = NaN
x( 5 ) = NaN
x( 6 ) = NaN
x( 7 ) = NaN
x( 8 ) = NaN
x( 9 ) = NaN
x( 10 ) = NaN
您对 MPI 的工作原理存在根本性误解。假设您有一个 30 行的矩阵，并且有 5 个进程。然后 scatter 调用将 6 行放入每个进程中。因此每个进程只需要一个 6 行的数组，分散的结果位于第 1--6 行 (Fortran) 或 0--5 (C) 中。
但是，您在每个进程中分配 30 行，并且您似乎假设这些行分散到适当的位置。事实并非如此。
使用MPI不仅仅是分配工作，还可以分配数据！每个进程只需要分配总数据的1/P部分。 1. 这意味着使用 MPI 可以解决比共享内存大得多的问题。 2. 这意味着分散方法不好：你想让每个进程在本地create它的一部分。试想一下，我有一台拥有 10 万核的超级计算机。对于您的程序，其中一个内核应该能够访问其他内核 10 万倍的内存。不会发生吧？