Rustで円周率を計算する: Chudnovskyの公式をプロセス並列 

今回はMPIを用いてChudnovskyの公式による円周率計算をマルチノード並列化できるようにしていきます。今回使用するのはMPI bindings for Rustです。(rustネイティブな分散並列用のクレートを知らないので、馴染みのあるMPIを選択しました。)実行環境は以下のとおりです。

Chip: Apple M3 Pro
Total Number of Cores: 11 (5 performance and 6 efficiency)
Memory: 36 GB
OS: macOS Sonoma 14.4.1

準備

OpenMPIのインストール

手元環境はmacなので、brewを用いてOpenMPIをインストールします。

brew install openmpi

mpiクレートの追加

cargo add mpi

実装

※rsmpiのexampleを参考に、単純なsendとreceiveを用いて実装しました。当方MPI初心者(数年前に数回触った程度)なのでベストプラクティスではないかもしれません。

スレッド並列版とは異なり、マスタープロセス(ランク0)が各プロセスの計算結果を配列で受け取ってから各プロセスの計算結果を統合し、円周率を計算します。rugのIntegerを他プロセスへ送信するため、IntegerString型に変換後、Vec<u8>に変換します。最終的に得られるVec<u8>の長さはIntegerで表現されていた元の数値の桁数に依存します。そのため、マスタープロセスからすると、データを受け取るために用意すべき配列の大きさを事前に知ることはできません。そこで、ワーカープロセスは、まず送信したいVec<u8>の長さをマスタープロセスに送信します。マスタープロセスはこの情報を用いて必要な長さのVec<u8>を用意します。その後、ワーカープロセスから送信されたVec<u8>を用意した配列で受け取るという手順です。

ワーカープロセスから集めたデータから最終的な答えを復元するため、マスタープロセスはVec<u8>をString型に変換した後、Integer型に戻します。以降の手順はThread並列版と同様です。

use mpi::{traits::*, topology::Rank};

pub fn chudnovsky_openmpi(n: u32, prec: u32, n_threads: u64) -> Float{
    let __n: u64 = (n / 14).into(); 
    let _n = __n + __n % 2;
    if n_threads == 1 {
        return chudnovsky_serial(n, prec);
    }
    let universe = mpi::initialize().unwrap();
    let world = universe.world();
    let world_size = world.size() as usize;
    let rank = world.rank();
    let root_rank = 0;
    let size = _n / (world_size as u64);
    let left = size * (rank as u64) + _n % (world_size as u64);
    let right = std::cmp::min(left + size, _n);

    let mut pi = Float::with_val(prec, 0);
    
    if rank == root_rank {
        let (xp, yp, zp) = compute_xyz(0, right, _n);
        let mut x = xp;
        let mut y = yp;
        let mut z = zp;
        for r in 1..world_size {
            let (msg, status) = world.process_at_rank(r as i32).receive_vec::<Rank>();
            let mut recv_buf_x = vec![0u8; msg[0] as usize];
            let mut recv_buf_y = vec![0u8; msg[2] as usize];
            let mut recv_buf_z = vec![0u8; msg[4] as usize];
            world.process_at_rank(r as i32).receive_into(&mut recv_buf_x[..]);
            world.process_at_rank(r as i32).receive_into(&mut recv_buf_y[..]);
            world.process_at_rank(r as i32).receive_into(&mut recv_buf_z[..]);
            // println!("Process {} got x {:?} from {}.", rank, recv_buf_x, status.source_rank());
            let mut xt = Integer::new();
            let mut yt = Integer::new();
            let mut zt = Integer::new();
            xt.assign(Integer::parse(String::from_utf8(recv_buf_x).unwrap()).unwrap());
            yt.assign(Integer::parse(String::from_utf8(recv_buf_y).unwrap()).unwrap());
            zt.assign(Integer::parse(String::from_utf8(recv_buf_z).unwrap()).unwrap());
            // println!("From rank{}: l={}, r={}, _n={}\n x: {}, y: {}, z: {}", r, size * (status.source_rank() as u64) + _n % (world_size as u64), size * (status.source_rank() as u64) + size + _n % (world_size as u64), _n, xt, yt, zt);
            x *= &xt;
            y = &xt * y + &yt * &z;
            z *= &zt;
        }
        pi = Float::with_val(prec, C_CUBED).sqrt() * x / 12 / y;
    } else {
        let (xp, yp, zp) = compute_xyz(left, right, _n);
        let x_sv =  xp.to_string().into_bytes();
        let y_sv =  yp.to_string().into_bytes();
        let z_sv =  zp.to_string().into_bytes();
        let msg = [x_sv.len() as u64, y_sv.len() as u64, z_sv.len() as u64];
        // println!("Process {} sent msg {:?} to {}.", rank, msg, root_rank);
        world.process_at_rank(root_rank).send(&msg[..]);
        world.process_at_rank(root_rank).send(&x_sv[..]);
        world.process_at_rank(root_rank).send(&y_sv[..]);
        world.process_at_rank(root_rank).send(&z_sv[..]);
        // println!("Process {} sent x {:?} to {}.", rank, x_sv, root_rank);
    }
    pi
}

今後

ここまで行ってきた実装で、円周率1億2800万桁程度までは計算できるようになりました。しかし、それ以上はオーバーフローが発生するのか、正しく計算できなくなります。したがって、これ以上を目指す場合は何らかの工夫が必要です。

この記事は役に立ちましたか?

もし参考になりましたら、下記のボタンで教えてください。

関連記事

コメント

この記事へのコメントはありません。

CAPTCHA