LoginSignup
3
3

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 3 years have passed since last update.

@ishigaki

JuliaでMPIによる並列計算

Last updated at Posted at 2021-10-25

お詫び

Gatherv!を使ったコードにバグがあったので，修正版のコードを載せました。

MPIによる並列計算

Juliaのコードで並列計算したいとおもって，MPIによる並列計算をテスト。
（Juliaのプロセス並列も試してみたが，オーバーヘッドが大きいのか並列効果が出なくてMPIをためそうとしてする状態です）

上記の記事を参考にさせてもらいました。

MPIの準備

MPIを使えるようにしておきます。その状態でJuliaのREPLで

add MPI

しておきます。

行列積の計算

コード

ベタですが，行列積の計算。自分がマスターのころ（15年ほど前）TAやってたMPIの講義の資料を引っ張り出して，後輩たちのコードを見つつ作成。

using MPI

function main()
    MPI.Init() # MPI初期化
    comm = MPI.COMM_WORLD
    nprocs = MPI.Comm_size(comm)
    myrank = MPI.Comm_rank(comm)
    println("myrank/nprocs=", myrank, "/", nprocs)

    l=3000
    m=3000
    n=3000
    A = zeros(l, n)
    B = ones(l, m)
    C = 0.5 .* ones(m, n)
    Atemp = zeros(l, n)
    
    ista, iend = parallel_range(1, m, nprocs, myrank)
    println("myrank=", myrank, ", ista, iend=", ista, ", ", iend)
    
    for k=ista:iend
        for j=1:n
            for i=1:l
                Atemp[i,j] = Atemp[i,j]+B[i,k]C[k,j]
            end
        end
    end
    
    A = MPI.Reduce(Atemp, MPI.SUM, 0, comm)
    
    MPI.Barrier(comm)
    
    if myrank == 0 
        println("A[1:5,1:5]=", A[1:5, 1:5])
        println("A[3000-5:3000,3000-5:3000]=", A[3000-5:3000,3000-5:3000])
    end

    MPI.Finalize() #MPI終了
end

function parallel_range(n1, n2, nprocs, irank)
    nb = div(n2-n1+1, nprocs)
    nm = (n2-n1+1)%nprocs
    ista = irank*nb+n1 + min(irank, nm)
    iend = ista+nb-1 
    if nm > irank
        iend = iend+1
    end
    return ista, iend
end

@time main()

各プロセスで計算した計算結果をAtempに保存して，MPI.ReduceでAtempの結果をAに集約（合算）します。

処理する要素を分割するparallel_rangeですが，当時のテキストをコピーしました。
可能な限り，同程度の要素数になるように分割します。要素数がプロセス数で割り切れなくても
うまく分割できます。

実行

ソースコードをMPI-matrixmul01.jlとして保存したとして，以下で実行。環境はLinuxです。
Macだとmpirunだと動かないかも。出てくるコメントに応じて適宜変更すると動きます。

8並列くらいまでは並列性能が出るのを確かめました。

mpirun -np 9 julia MPI-matrixmul01.jl

MPI.Gatherを使った場合

データの集約にMPI.Gatherを使った場合もやってみました。MPI.Gatherの返り値に計算結果が入りますが，なぜか1次元配列になります。なので，あとからreshapeしてます。
理由ご存知だったら教えていただけるとうれしいです。

MPI.Gatherだと送信バッファのサイズがすべて等しいと仮定するので，ランクごとに送信バッファのサイズが変わるとエラーします。

using MPI

function main()
    MPI.Init() # MPI初期化
    comm = MPI.COMM_WORLD
    nprocs = MPI.Comm_size(comm)
    myrank = MPI.Comm_rank(comm)
    println("myrank/nprocs=", myrank, "/", nprocs)

    l=3000
    m=2079
    n=3000
    A = zeros(l, n)
    B = ones(l, m)
    C = 0.5 .* ones(m, n)
       
    jsta, jend = parallel_range(1, n, nprocs, myrank)
    println("myrank=", myrank, ", jsta, jend=", jsta, ", ", jend)
    
    for j=jsta:jend
        for k=1:m
            for i=1:l
                A[i,j] = A[i,j]+B[i,k]C[k,j]
            end
        end
    end
    
     A = MPI.Gather(A[:,jsta:jend], 0, comm)
    
    MPI.Barrier(comm)
    
    if myrank == 0 
        A=reshape(A, 3000, 3000)
        println("A[1:5,1:5]=", A[1:5, 1:5])
        println("A[3000-5:3000,3000-5:3000]=", A[3000-5:3000,3000-5:3000])
    end

    MPI.Finalize() #MPI終了
end

function parallel_range(n1, n2, nprocs, irank)
    nb = div(n2-n1+1, nprocs)
    nm = (n2-n1+1)%nprocs
    ista = irank*nb+n1 + min(irank, nm)
    iend = ista+nb-1 
    if nm > irank
        iend = iend+1
    end
    return ista, iend
end

@time main()

MPI.Send, MPI.Recvを使った場合

send, recvを使った場合。受け取るrecv側の配列にはA[:,jsta:jend]のようなスライスが使えなかったので，別途配列を作って受け取ってます（tempの部分）。tempのデータを元の配列に戻す処理が必要になってくるので若干手間が増えます。

using MPI

function main()
    MPI.Init() # MPI初期化
    comm = MPI.COMM_WORLD
    nprocs = MPI.Comm_size(comm)
    myrank = MPI.Comm_rank(comm)
    println("myrank/nprocs=", myrank, "/", nprocs)

    l=3000
    m=2080
    n=3000
    #l=5
    #m=4
    #n=5
    A = zeros(l, n)
    B = ones(l, m)
    C = 0.5 .* ones(m, n)
    
    temp = zeros(l,n)
       
    jsta, jend = parallel_range(1, n, nprocs, myrank)
    println("myrank=", myrank, ", jsta, jend=", jsta, ", ", jend)
    
    for j=jsta:jend
        for k=1:m
            for i=1:l
                A[i,j] = A[i,j]+B[i,k]C[k,j]
            end
        end
    end
    
    if myrank == 0
        for rank=1:nprocs-1
            MPI.Recv!(temp, rank, 0, comm)
            println("temp",temp[1:2,1])
            jsta, jend = parallel_range(1, n, nprocs, rank)
            println("jsta,jend=", jsta, ", ", jend)
            A[:,jsta:jend] .= temp[:,jsta-jsta+1:jend-jsta+1]
        end
    else 
        MPI.Send(A[:,jsta:jend], 0, 0, comm)
    end
    
    if myrank == 0 
        println("A[1:5,1:5]=", A[1:5, 1:5])
        println("A[3000-5:3000,3000-5:3000]=", A[l-5:l,l-5:l])
    end

    MPI.Finalize() #MPI終了
end

function parallel_range(n1, n2, nprocs, irank)
    nb = div(n2-n1+1, nprocs)
    nm = (n2-n1+1)%nprocs
    ista = irank*nb+n1 + min(irank, nm)
    iend = ista+nb-1 
    if nm > irank
        iend = iend+1
    end
    return ista, iend
end

@time main()

MPI.Gatherv!を使った場合

Gatherv!を使えば，各ランクで要素数が違っても対応できます。そのために，各ランクで送受信の対象となるデータの要素数を先に数えておく必要があります。それがsizesという配列の計算の部分です。

Gatherv!では受信バッファにはVBufferを作ってやります。VBuffer(A, sizes)でAにデータが集約されます。

配列データを送受信するとき，受信する側では配列の列に沿った方向に順番にデータが格納されます。そのため，配列は列に沿った方向に分割する必要があります。

修正版のコード

using MPI

function main()
    MPI.Init() # MPI初期化
    comm = MPI.COMM_WORLD
    nprocs = MPI.Comm_size(comm)
    myrank = MPI.Comm_rank(comm)
    println("myrank/nprocs=", myrank, "/", nprocs)

    l=3000
    m=2079
    n=3000
    A = zeros(l, n)
    B = ones(l, m)
    C = 0.5 .* ones(m, n)
    Atemp = similar(A)
    
    sizes = zeros(Int64, nprocs)
    
    if myrank == 0
        for i=0:nprocs-1
            jsta, jend = parallel_range(1, n, nprocs, i)
            sizes[i+1] = l*(jend-jsta+1)
        end
        println("sizes=", sizes, ", ", sum(sizes))
    end
    
    jsta, jend = parallel_range(1, n, nprocs, myrank)
    println("myrank=", myrank, ", jsta, jend=", jsta, ", ", jend)

    for j=jsta:jend
        for k=1:m
            for i=1:l
                Atemp[i,j] = Atemp[i,j]+B[i,k]C[k,j]
            end
        end
    end
    
    if myrank == 0
        MPI.Gatherv!(Atemp[:,jsta:jend], VBuffer(A, sizes), 0, comm)
    else
        MPI.Gatherv!(Atemp[:,jsta:jend], VBuffer(nothing), 0, comm)
    end

    #MPI.Barrier(comm)
    
    if myrank == 0 
        println(findall(x->(x!=1039.5), A)) #値のチェック
        println("A[1:5,1:5]=", A[1:5, 1:5])
        println("A[3000-5:3000,3000-5:3000]=", A[3000-5:3000,3000-5:3000])
    end

    MPI.Finalize() #MPI終了
end

function parallel_range(n1, n2, nprocs, irank)
    nb = div(n2-n1+1, nprocs)
    nm = (n2-n1+1)%nprocs
    ista = irank*nb+n1 + min(irank, nm)
    iend = ista+nb-1 
    if nm > irank
        iend = iend+1
    end
    return ista, iend
end

@time main()

お詫び

下記のコードだと，jsta, jendの範囲が正しく計算されず，正しい結果が得られません。お詫びします。
上記に修正したコードを載せます。

修正前コード

using MPI

function main()
    MPI.Init() # MPI初期化
    comm = MPI.COMM_WORLD
    nprocs = MPI.Comm_size(comm)
    myrank = MPI.Comm_rank(comm)
    println("myrank/nprocs=", myrank, "/", nprocs)

    l=3000
    m=2079
    n=3000
    A = zeros(l, n)
    B = ones(l, m)
    C = 0.5 .* ones(m, n)
    Atemp = similar(A)
    
    sizes = zeros(Int64, nprocs)

    jsta, jend = parallel_range(1, n, nprocs, myrank)
    println("myrank=", myrank, ", jsta, jend=", jsta, ", ", jend)
    
    if myrank == 0
        for i=0:nprocs-1
            jsta, jend = parallel_range(1, n, nprocs, i)
            sizes[i+1] = l*(jend-jsta+1)
        end
        println("sizes=", sizes, ", ", sum(sizes))
    end

    for j=jsta:jend
        for k=1:m
            for i=1:l
                Atemp[i,j] = Atemp[i,j]+B[i,k]C[k,j]
            end
        end
    end
    
    if myrank == 0
        MPI.Gatherv!(Atemp[:,jsta:jend], VBuffer(A, sizes), 0, comm)
    else
        MPI.Gatherv!(Atemp[:,jsta:jend], VBuffer(nothing), 0, comm)
    end

    #MPI.Barrier(comm)
    
    if myrank == 0 
        A=reshape(A, 3000, 3000)
        println("A[1:5,1:5]=", A[1:5, 1:5])
        println("A[3000-5:3000,3000-5:3000]=", A[3000-5:3000,3000-5:3000])
    end

    MPI.Finalize() #MPI終了
end

function parallel_range(n1, n2, nprocs, irank)
    nb = div(n2-n1+1, nprocs)
    nm = (n2-n1+1)%nprocs
    ista = irank*nb+n1 + min(irank, nm)
    iend = ista+nb-1 
    if nm > irank
        iend = iend+1
    end
    return ista, iend
end

@time main()

計算速度

MPI.Gatherv!使った場合の結果の抜粋。

1コア
myrank/nprocs=0/1
myrank=0, jsta, jend=1, 3000
sizes=[9000000], 9000000
A[1:5,1:5]=[1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5]
A[3000-5:3000,3000-5:3000]=[1039.5 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5 1039.5]
 49.078051 seconds (936.48 k allocations: 402.841 MiB, 0.49% gc time, 1.13% compilation time)

5コア
[1800000, 1800000, 1800000, 1800000, 1800000], 9000000
A[1:5,1:5]=[1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5]
A[3000-5:3000,3000-5:3000]=[1039.5 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5 1039.5]
  8.579884 seconds (936.51 k allocations: 347.910 MiB, 1.61% gc time, 5.38% compilation time)

9コア
sizes=[1002000, 1002000, 1002000, 999000, 999000, 999000, 999000, 999000, 999000], 9000000
A[1:5,1:5]=[1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5]
A[3000-5:3000,3000-5:3000]=[1039.5 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5 1039.5]
  5.390722 seconds (936.55 k allocations: 341.800 MiB, 3.50% gc time, 8.65% compilation time)  5.410653 seconds (79.32 k allocations: 292.609 MiB, 3.54% gc time, 1.78% compilation time)

10コア
sizes=[900000, 900000, 900000, 900000, 900000, 900000, 900000, 900000, 900000, 900000], 9000000
A[1:5,1:5]=[1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5]
A[3000-5:3000,3000-5:3000]=[1039.5 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5 1039.5]
  6.697126 seconds (79.31 k allocations: 291.854 MiB, 2.35% gc time, 1.66% compilation time)
  6.704198 seconds (936.56 k allocations: 341.045 MiB, 2.35% gc time, 8.01% compilation time)

20コア
sizes=[450000, 450000, 450000, 450000, 450000, 450000, 450000, 450000, 450000, 450000, 450000, 450000, 450000, 450000, 450000, 450000, 450000, 450000, 450000, 450000], 9000000
A[1:5,1:5]=[1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5]
A[3000-5:3000,3000-5:3000]=[1039.5 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5 1039.5]
  5.459388 seconds (79.31 k allocations: 288.420 MiB, 3.14% gc time, 2.05% compilation time)  5.461391 seconds (79.31 k allocations: 288.420 MiB, 3.16% gc time, 2.05% compilation time)

40コア
sizes=[225000, 225000, 225000, 225000, 225000, 225000, 225000, 225000, 225000, 225000, 225000, 225000, 225000, 225000, 225000, 225000, 225000, 225000, 225000, 225000, 225000, 225000, 225000, 225000, 225000, 225000, 225000, 225000, 225000, 225000, 225000, 225000, 225000, 225000, 225000, 225000, 225000, 225000, 225000, 225000], 9000000
A[1:5,1:5]=[1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5]
A[3000-5:3000,3000-5:3000]=[1039.5 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5 1039.5; 1039.5 1039.5 1039.5 1039.5 1039.5 1039.5]
  6.470760 seconds (79.31 k allocations: 286.704 MiB, 3.23% gc time, 1.76% compilation time)
3
3
2

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
3
3

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?