2013-02-26 88 views
0

我試圖通過jacobi迭代解決ax = b,我的串行代碼工作正常,但MPI版本甚至不會運行。誰能幫我?Fortran中的Jacobi迭代和MPI

串行

program jacobis 

implicit none 

integer, parameter :: n=10 
integer :: i,j,k,ni,s,seed 
double precision :: tol,t1,t2,sig 
double precision, dimension(0:n-1,0:n-1) :: A 
double precision, dimension(0:n-1) :: B, x, xb, buff 

ni=1000 

seed=time() 
call srand(seed) 

do i=0, n-1 
    do j=0, n-1 
    A(i,j)=rand(0) 
    B(i)=rand(0) 
    end do 
end do 

do i = 0, n-1 
A(i,i) = sum(A(i,:)) + 1 
enddo 

!do i=0,n-1 
!A(i,i)=4 
!end do 

print *, "a", A 
print *, "b", B 

x=B 
call cpu_time(t1) 
do k=1,ni 
xb=x 
do i=0,n-1 
    s=0 
    do j=0,n-1 
    if (j/=i) then 
     s=s+A(i,j)*xb(j) 
     endif 
    end do 
    x(i)=(B(i)-s)/A(i,i) 

    sig=(x(i)-xb(i))*(x(i)-xb(i)) 
    tol=tol+sig 
    tol=sqrt(tol) 
end do 


print *, "x", x 

!print *, "tol=", tol 

print *, "iter =",k 

if (tol<1.000001) EXIT 
if (k==(ni-1)) then 
    print *, "Numero Maximo de Iteracoes" 
    EXIT 
endif 
end do 

call cpu_time(t2) 
print *, "t=",t2-t1 


end 

MPI版本

program jacobis 

use mpi 
implicit none 

integer, parameter :: n=2 
integer :: i_local,i_global,j,k,ni,s,m 
double precision :: tol,t,t2,sig 
double precision, dimension(:,:), ALLOCATABLE :: A_local 
double precision, dimension(:), ALLOCATABLE :: B_local, x_local, x_temp1,x_temp2,x_old,x_new, buff 
INTEGER, DIMENSION (MPI_STATUS_SIZE) :: STATUS 
integer :: rank,procs,tag,ierror 


CALL MPI_INIT(ierror) 
CALL MPI_COMM_RANK(MPI_COMM_WORLD,rank,ierror) 
CALL MPI_COMM_SIZE(MPI_COMM_WORLD,procs,ierror) 

ni=100 
m=n/procs 

ALLOCATE (A_local(0:n-1,0:n-1)) 
ALLOCATE (B_local(0:m-1)) 
ALLOCATE (x_temp1(0:m-1)) 
ALLOCATE (x_temp2(0:m-1)) 

A_local=0 
B_local=2 

do i_global=0,n-1 
A_local(i_global,i_global)=2 
end do 

CALL MPI_ALLGATHER(B_local, m, MPI_DOUBLE, x_temp1, m, MPI_DOUBLE, MPI_COMM_WORLD,ierror) 

x_new=x_temp1 
x_old=x_temp2 

print *, "a", A_local 
print *, "b", B_local 


t=mpi_wtime() 
do k=1,ni 
x_old=x_new 
do i_local=0,m-1 
    i_global=i_local+rank*m 
    !x_local(i_local)=b_local(i_local) 
    s=0 
    do j=0,n-1 
    if (j/=i_local) then 
     s=s+A_local(i_local,j)*x_old(j) 
     endif 
    end do 
    x_local(i_local)=(B_local(i_local)-s)/A_local(i_local,i_global) 

end do 
CALL MPI_ALLGATHER(x_local,m, MPI_DOUBLE, x_new, m, MPI_DOUBLE, MPI_COMM_WORLD,ierror) 

do i_global=0,n-1 
    sig=(x_new(i_global)-x_old(i_global))*(x_new(i_global)-x_old(i_global)) 
    tol=tol+sig 
    tol=sqrt(tol) 
end do 

print *, "x", x_local 

print *, "tol=", tol 

print *, "iter =",k 

if (tol<1.000001) EXIT 
if (k==(ni-1)) then 
    print *, "Numero Maximo de Iteracoes" 
    EXIT 
endif 
end do 

t2=mpi_wtime()-t; 
print *, "t=",t2 

CALL MPI_FINALIZE(ierror) 
end 

任何人都可以指出我在做什麼錯?這是索引問題嗎?請我真的需要今天解決這個問題,否則我會忽略這個過程。我已經花了無數小時,並且無法完成工作。

好吧,你是對的!現在我有一個分段錯誤,但找不到它!用新版本取代了代碼

+0

你仍然在賦值表達式的LHS上有非分配數組。你使用什麼編譯器?使用'-O0 -g -C'標誌編譯以獲取錯誤消息,並提供有關錯誤的更多提示。因爲這是一個類的分配,所以我們不能調試你的代碼,但只能提供什麼是錯誤的以及如何繼續的提示。祝你好運。 – milancurcic 2013-02-27 19:18:48

回答

0

我已經解決了這個問題,現在它可以正確地計算迭代,並通過使用相同矩陣的串行程序進行驗證。這是一個分配和索引問題。感謝之前的回答,非常有幫助。

program jacobis 

use mpi 
implicit none 

integer, parameter :: n=1000 
integer :: i_local,i_global,j,k,ni,s,m,seed 
double precision :: tol,t,t2,sig 
double precision, dimension(:,:), ALLOCATABLE :: A_local 
double precision, dimension(:), ALLOCATABLE :: B_local, x_local, x_temp1,x_old,x_new, buff 
INTEGER, DIMENSION (MPI_STATUS_SIZE) :: STATUS 
integer :: rank,procs,tag,ierror 


CALL MPI_INIT(ierror) 
CALL MPI_COMM_RANK(MPI_COMM_WORLD,rank,ierror) 
CALL MPI_COMM_SIZE(MPI_COMM_WORLD,procs,ierror) 

ni=1000 
m=n/procs 

ALLOCATE (A_local(0:n-1,0:n-1)) 
ALLOCATE (B_local(0:n-1)) 
ALLOCATE (x_local(0:n-1)) 
ALLOCATE (x_temp1(0:n-1)) 
ALLOCATE (x_new(0:n-1)) 

!A_local=23 
!B_local=47 

seed=time() 
call srand(seed) 

do k=0, n-1 
do j=0, n-1 
    A_local(k,j)=rand(0) 
    B_local(k)=rand(0) 
end do 
end do 

do i_global = 0, m-1 
A_local(i_global,i_global) = sum(A_local(i_global,:)) + n 
enddo 

CALL MPI_ALLGATHER(B_local, m, MPI_DOUBLE, x_temp1, m, MPI_DOUBLE, MPI_COMM_WORLD,ierror) 

x_new=x_temp1 

print *, "a", A_local 
print *, "b", B_local 


t=mpi_wtime() 
do k=1,ni 
x_old=x_new 
do i_local=0,m-1 
    i_global=i_local+rank*m 
    !x_local(i_local)=b_local(i_local) 
    s=0 
    do j=0,n-1 
    if (j/=i_local) then 
     s=s+A_local(i_local,j)*x_old(j) 
     endif 
    end do 
    x_local(i_local)=(B_local(i_local)-s)/A_local(i_local,i_global) 
end do 
CALL MPI_ALLGATHER(x_local,m, MPI_DOUBLE, x_new, m, MPI_DOUBLE, MPI_COMM_WORLD,ierror) 
do j=0,n-1 
    sig=(x_new(j)-x_old(j))*(x_new(j)-x_old(j)) 
    tol=tol+sig 
    tol=sqrt(tol) 
end do 

print *, "x", x_local 

print *, "tol=", tol 

print *, "iter =",k 

if (tol<1.01) EXIT 
if (k==(ni-1)) then 
    print *, "Numero Maximo de Iteracoes" 
    EXIT 
endif 
end do 

t2=mpi_wtime()-t; 
print *, "t=",t2 

CALL MPI_FINALIZE(ierror) 
end 
+0

很高興你解決了你的問題。請注意,您可以點擊複選標記以接受任一答案。 – milancurcic 2013-02-27 21:15:33

+1

考慮點擊@ IRO-bot答案上的複選標記,它幫助你解決問題(不是你自己的解決方案是基於他的建議)。 – bcumming 2013-12-13 11:39:18

2

您的程序有幾個我可以看到的問題。你包括該錯誤消息表示在此調用未分配的接收緩衝區:

CALL MPI_ALLGATHER(B_local, m, MPI_DOUBLE, x_temp1, m, MPI_DOUBLE, MPI_COMM_WORLD) 

陣列x_temp1,接收緩衝區,需要使用在這方面之前進行分配。

解決這個問題只會讓你得到更多的信息,並且你將得到更少的信息分割錯誤。在您的MPI實施中查找MPI_AllGather的正確用法將很有用。大多數MPI例程在最後都有一個整數錯誤狀態參數:

MPI_ALLGATHER(SENDBUF, SENDCOUNT, SENDTYPE, RECVBUF, RECVCOUNT, 
     RECVTYPE, COMM, IERROR) 
    <type> SENDBUF (*), RECVBUF (*) 
    INTEGER SENDCOUNT, SENDTYPE, RECVCOUNT, RECVTYPE, COMM, 
    INTEGER IERROR 

這應該讓你順利完成任務。請確保分配您使用的所有allocatable陣列,併爲您的MPI實現和編譯器手冊使用適當的文檔。

0

您的程序存在嚴重問題,您可能會得到錯誤的結果。變量s被聲明爲整數,而它被分配了非整數值。以雙精度重新聲明它以獲得正確的結果。 (發佈給那些曾經複製過這段代碼的人)