Kurseinheit 7: Multi-GPU Computing

/* Allokation von Host und GPU-Speicher, Initialisierung etc */
if (rank == 0) {
    cudaMemcpy(s_buf_h, s_buf_d, size, cudaMemcpyDeviceToHost);
    MPI_Send(s_buf_h, size, MPI_BYTE, 1, tag, MPI_COMM_WORLD);
}
if (rank == 1) {
    MPI_Recv(r_buf_h, size, MPI_BYTE, 0, tag, MPI_COMM_WORLD, &stat);
    cudaMemcpy(r_buf_d, r_buf_h, size, cudaMemcpyHostToDevice);
}

/* Allokation von GPU-Speicher, Initialisierung etc */
if (rank == 0) {
    MPI_Send(s_buf_d, size, MPI_BYTE, 1, tag, MPI_COMM_WORLD);
}
if (rank == 1) {
    MPI_Recv(r_buf_d, size, MPI_BYTE, 0, tag, MPI_COMM_WORLD, &stat);
}

int provided;
MPI_Init_thread(&argc, &argv, MPI_THREAD_MULTIPLE, &provided);
if (provided < MPI_THREAD_MULTIPLE) {
    fprintf(stderr, "MPI_THREAD_MULTIPLE nicht unterstützt!\n");
    MPI_Abort(MPI_COMM_WORLD, 1);
}

for (int iy = 1; iy < ny-1; iy++)
    for (int ix = 1; ix < nx-1; ix++)
        a_new[iy*nx+ix] = 0.25 * (
            a[iy*nx+(ix+1)] + a[iy*nx+ix-1] +
            a[(iy-1)*nx+ix] + a[(iy+1)*nx+ix]
        );

int iy = blockIdx.y * blockDim.y + threadIdx.y;
int ix = blockIdx.x * blockDim.x + threadIdx.x;
if (iy < ny-1 && ix < nx-1) {
    float new_val = 0.25f * (
        a[iy * nx + ix + 1] + a[iy * nx + ix - 1] +
        a[(iy + 1) * nx + ix] + a[(iy - 1) * nx + ix]
    );
    a_new[iy * nx + ix] = new_val;
}

Solange nicht konvergiert:
    1. Jacobi-Schritt (Kernel)
    2. Randbedingungen setzen
    3. a_new und a vertauschen
    4. Nächste Iteration

// Lokale L2-Norm auf der GPU berechnen (als Teil des Jacobi-Kernels)
// l2_norm_d enthält die lokale Summe der quadrierten Residuen

// GPU-Ergebnis auf Host kopieren
double l2_norm_h = 0.0;
cudaMemcpy(&l2_norm_h, l2_norm_d, sizeof(double), cudaMemcpyDeviceToHost);

// Globale L2-Norm über alle Ranks summieren
double global_l2_norm = 0.0;
MPI_Allreduce(&l2_norm_h, &global_l2_norm, 1, MPI_DOUBLE,
              MPI_SUM, MPI_COMM_WORLD);

global_l2_norm = sqrt(global_l2_norm);

// Konvergenzbedingung prüfen
if (global_l2_norm < tolerance) break;

// L2-Norm direkt GPU-seitig reduzieren (bleibt auf GPU)
ncclAllReduce(l2_norm_d, l2_norm_global_d, 1,
              ncclDouble, ncclSum, nccl_comm, compute_stream);

┌──────────────────┐
│   GPU 0 (Rank 0) │  ← oberer Rand → GPU 1
├──────────────────┤
│   GPU 1 (Rank 1) │  ← oberer Rand → GPU 0
│                  │  ← unterer Rand → GPU 2
├──────────────────┤
│   GPU 2 (Rank 2) │
└──────────────────┘

// Sende letzte Zeile an unteren Nachbar, empfange obere Halo-Zeile
MPI_Sendrecv(a_new_d + offset_last_row,  m-2, MPI_DOUBLE, b_nb, 1,
             a_new_d + offset_top_boundary, m-2, MPI_DOUBLE, t_nb, 1,
             MPI_COMM_WORLD, MPI_STATUS_IGNORE);

// Sende erste Zeile an oberen Nachbar, empfange untere Halo-Zeile
MPI_Sendrecv(a_new_d + offset_first_row, m-2, MPI_DOUBLE, t_nb, 0,
             a_new_d + offset_bottom_boundary, m-2, MPI_DOUBLE, b_nb, 0,
             MPI_COMM_WORLD, MPI_STATUS_IGNORE);

// Pack-Kernel: kopiert Spalten-Daten in zusammenhängenden Puffer
pack<<<gs, bs, 0, s>>>(to_left_d, u_new_d, n, m);
cudaStreamSynchronize(s);

MPI_Sendrecv(to_left_d,   n-2, MPI_DOUBLE, l_nb, 0,
             from_right_d, n-2, MPI_DOUBLE, r_nb, 0,
             MPI_COMM_WORLD, MPI_STATUS_IGNORE);

// Unpack-Kernel: schreibt empfangene Daten zurück ins Grid
unpack<<<gs, bs, 0, s>>>(u_new_d, from_right_d, n, m);
cudaStreamSynchronize(s);

// Schritt 1: Kommunikator nur für Prozesse auf dem gleichen Knoten
MPI_Comm local_comm;
MPI_Comm_split_type(MPI_COMM_WORLD, MPI_COMM_TYPE_SHARED, rank,
                    MPI_INFO_NULL, &local_comm);

// Schritt 2: Lokalen Rank auf dem Knoten abfragen
int local_rank = -1;
MPI_Comm_rank(local_comm, &local_rank);
MPI_Comm_free(&local_comm);

// Schritt 3: GPU entsprechend dem lokalen Rank auswählen
int num_devs = 0;
cudaGetDeviceCount(&num_devs);
cudaSetDevice(local_rank % num_devs);

Adressraum: 0x0000 ──────────────────────────── 0xFFFF
                   │ System Memory  │ GPU Memory │

GPU A ──NVLink / PCIe Switch── GPU B
              (direkt)

int can_access = 0;
cudaDeviceCanAccessPeer(&can_access, src_device, dst_device);

if (can_access) {
    cudaSetDevice(src_device);
    cudaDeviceEnablePeerAccess(dst_device, 0);  // 0 = flags (reserviert)

    cudaSetDevice(dst_device);
    cudaDeviceEnablePeerAccess(src_device, 0);
} else {
    // Fallback: über Host-Speicher kopieren
}

// Direkte GPU-zu-GPU-Kopie (kein Host-Puffer nötig)
cudaMemcpyPeerAsync(dst_ptr, dst_device,
                    src_ptr, src_device,
                    size, stream);

GPU ──PCIe Switch── HCA ──Netzwerk── HCA ──PCIe Switch── GPU
                  (direkt, kein Host-Durchlauf)

// API-Aufruf ist in allen drei Fällen gleich:
MPI_Send(s_buf_d, size, MPI_BYTE, 1, tag, MPI_COMM_WORLD);
MPI_Recv(r_buf_d, size, MPI_BYTE, 0, tag, MPI_COMM_WORLD, &stat);

[──── GPU Jacobi ────][─MPI─][──── GPU Jacobi ────][─MPI─]

[Rand][──── Jacobi Innen ────][Rand][──── Jacobi Innen ────]
      [─── MPI ───]                 [─── MPI ───]

// 1. Rand-Kernel auf push_top_stream (oberste Zeile)
launch_jacobi_kernel(a_new, a, l2_norm_d,
    iy_start, (iy_start + 1), nx, push_top_stream);

// 2. Rand-Kernel auf push_bottom_stream (unterste Zeile)
launch_jacobi_kernel(a_new, a, l2_norm_d,
    (iy_end - 1), iy_end, nx, push_bottom_stream);

// 3. Innen-Kernel auf compute_stream (parallel zu MPI)
launch_jacobi_kernel(a_new, a, l2_norm_d,
    (iy_start + 1), (iy_end - 1), nx, compute_stream);

// 4. Warte auf oberen Rand, dann Austausch mit oberem Nachbar
const int top    = rank > 0 ? rank - 1 : (size - 1);
const int bottom = (rank + 1) % size;

cudaStreamSynchronize(push_top_stream);
MPI_Sendrecv(a_new + iy_start * nx,  nx, MPI_REAL_TYPE, top, 0,
             a_new + (iy_end * nx),  nx, MPI_REAL_TYPE, bottom, 0,
             MPI_COMM_WORLD, MPI_STATUS_IGNORE);

// 5. Warte auf unteren Rand, dann Austausch mit unterem Nachbar
cudaStreamSynchronize(push_bottom_stream);
MPI_Sendrecv(a_new + (iy_end - 1) * nx, nx, MPI_REAL_TYPE, bottom, 0,
             a_new,                      nx, MPI_REAL_TYPE, top, 0,
             MPI_COMM_WORLD, MPI_STATUS_IGNORE);

cudaStreamSynchronize(push_top_stream);  // blockiert CPU
MPI_Sendrecv(...);                        // blockiert CPU weiter

MPI_Init(&argc, &argv);
MPI_Comm_size(MPI_COMM_WORLD, &size);
MPI_Comm_rank(MPI_COMM_WORLD, &rank);

// Unique ID auf Rank 0 erzeugen und an alle verteilen
ncclUniqueId nccl_uid;
if (rank == 0) ncclGetUniqueId(&nccl_uid);
MPI_Bcast(&nccl_uid, sizeof(ncclUniqueId), MPI_BYTE, 0, MPI_COMM_WORLD);

// NCCL-Kommunikator initialisieren
ncclComm_t nccl_comm;
ncclCommInitRank(&nccl_comm, size, nccl_uid, rank);

/* ... Programmcode ... */

ncclCommDestroy(nccl_comm);
MPI_Finalize();

ncclSend(void* sbuff, size_t count, ncclDataType_t type,
         int peer, ncclComm_t comm, cudaStream_t stream);

ncclRecv(void* rbuff, size_t count, ncclDataType_t type,
         int peer, ncclComm_t comm, cudaStream_t stream);

ncclAllReduce(sbuff, rbuff, count, type, op, comm, stream);
ncclBroadcast(sbuff, rbuff, count, type, root, comm, stream);
ncclReduce(sbuff, rbuff, count, type, op, root, comm, stream);
ncclReduceScatter(sbuff, rbuff, count, type, op, comm, stream);
ncclAllGather(sbuff, rbuff, count, type, comm, stream);

// SendRecv zwischen zwei Ranks
ncclGroupStart();
ncclSend(sendbuff, sendcount, sendtype, peer, comm, stream);
ncclRecv(recvbuff, recvcount, recvtype, peer, comm, stream);
ncclGroupEnd();

// Broadcast (manuell mit Send/Recv implementiert)
ncclGroupStart();
if (rank == root) {
    for (int r = 0; r < nranks; r++)
        ncclSend(sendbuff[r], size, type, r, comm, stream);
}
ncclRecv(recvbuff, size, type, root, comm, stream);
ncclGroupEnd();

// Nachbar-Austausch in mehrdimensionalem Gitter
ncclGroupStart();
for (int d = 0; d < ndims; d++) {
    ncclSend(sendbuff[d], sendcount, sendtype, next[d], comm, stream);
    ncclRecv(recvbuff[d], recvcount, recvtype, prev[d], comm, stream);
}
ncclGroupEnd();

// Gesamten Jacobi-Schritt auf compute_stream berechnen
launch_jacobi_kernel(a_new, a, l2_norm_d,
    iy_start, iy_end, nx, compute_stream);

// Halo-Austausch auf demselben Stream (kein Synchronisationspunkt nötig!)
ncclGroupStart();
ncclRecv(a_new,                    nx, NCCL_REAL_TYPE, top, nccl_comm, compute_stream);
ncclSend(a_new + (iy_end-1)*nx,   nx, NCCL_REAL_TYPE, btm, nccl_comm, compute_stream);
ncclRecv(a_new + (iy_end*nx),     nx, NCCL_REAL_TYPE, btm, nccl_comm, compute_stream);
ncclSend(a_new + iy_start * nx,   nx, NCCL_REAL_TYPE, top, nccl_comm, compute_stream);
ncclGroupEnd();

// Prioritätsbereich abfragen
int leastPriority   = 0;
int greatestPriority = leastPriority;
cudaDeviceGetStreamPriorityRange(&leastPriority, &greatestPriority);

// Streams mit unterschiedlichen Prioritäten erstellen
cudaStream_t compute_stream;
cudaStream_t push_stream;
cudaStreamCreateWithPriority(&compute_stream, cudaStreamDefault, leastPriority);
cudaStreamCreateWithPriority(&push_stream,    cudaStreamDefault, greatestPriority);

// Rand-Kernel auf push_stream (höchste Priorität)
launch_jacobi_kernel(a_new, a, l2_norm_d,
    iy_start, iy_start + 1,         nx, push_stream);
launch_jacobi_kernel(a_new, a, l2_norm_d,
    iy_end - 1,  iy_end,            nx, push_stream);

// Innen-Kernel auf compute_stream (läuft parallel zur Kommunikation)
launch_jacobi_kernel(a_new, a, l2_norm_d,
    iy_start + 1, iy_end - 1,       nx, compute_stream);

// Kommunikation auf push_stream (nach den Rand-Kerneln, parallel zum Innen-Kernel)
ncclGroupStart();
ncclRecv(a_new,                    nx, NCCL_REAL_TYPE, top, nccl_comm, push_stream);
ncclSend(a_new + (iy_end-1)*nx,   nx, NCCL_REAL_TYPE, btm, nccl_comm, push_stream);
ncclRecv(a_new + (iy_end*nx),     nx, NCCL_REAL_TYPE, btm, nccl_comm, push_stream);
ncclSend(a_new + iy_start * nx,   nx, NCCL_REAL_TYPE, top, nccl_comm, push_stream);
ncclGroupEnd();

#include <nccl.h>

MPICXX_FLAGS = -I$(CUDA_HOME)/include -I$(NCCL_HOME)/include
LD_FLAGS     = -L$(CUDA_HOME)/lib64 -lcudart -lnccl

# CUDA-Kernels kompilieren
$(NVCC) $(NVCC_FLAGS) jacobi_kernels.cu -c -o jacobi_kernels.o

# MPI-Code kompilieren und linken
$(MPICXX) $(MPICXX_FLAGS) jacobi.cpp jacobi_kernels.o $(LD_FLAGS) -o jacobi

┌────────────────────────────────────────────────────────────────┐
│               Partitioned Global Address Space                  │
│  ┌─────────┐       ┌─────────┐             ┌─────────┐        │
│  │GPU 0    │       │GPU 1    │             │GPU N    │        │
│  │(PE 0)   │       │(PE 1)   │             │(PE N)   │        │
│  │┌───────┐│       │┌───────┐│             │┌───────┐│        │
│  ││Symm.  ││       ││Symm.  ││  ←────────→ ││Symm.  ││        │
│  ││Heap   ││       ││Heap   ││             ││Heap   ││        │
│  │└───────┘│       │└───────┘│             │└───────┘│        │
│  │┌───────┐│       │┌───────┐│             │┌───────┐│        │
│  ││Private││       ││Private││             ││Private││        │
│  │└───────┘│       │└───────┘│             │└───────┘│        │
│  └─────────┘       └─────────┘             └─────────┘        │
└────────────────────────────────────────────────────────────────┘

// Symmetrische Objekte: kollektiv auf allen PEs allokiert
// (muss auf ALLEN PEs mit gleicher Größe aufgerufen werden!)
void* shared_data = nvshmem_malloc(shared_size);

// Privater Speicher: wie gewohnt
void* private_data;
cudaMalloc(&private_data, private_size);

// Daten von diesem PE zu einem anderen PE schicken (Put)
// nvshmem_TYPE_put(dest_ptr, src_ptr, nelems, pe)
__global__ void kernel_put(double* data, int my_pe, int neighbor_pe, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < n) {
        // Schreibe data[idx] direkt in den symmetrischen Heap von neighbor_pe
        nvshmem_double_p(&data[idx], data[idx], neighbor_pe);
    }
    // Sicherstellen, dass alle Puts abgeschlossen sind
    nvshmem_quiet();
}

// Daten von einem anderen PE lesen (Get)
__global__ void kernel_get(double* local, double* remote_sym, int src_pe, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < n) {
        // Lese Wert aus dem symmetrischen Heap von src_pe
        local[idx] = nvshmem_double_g(&remote_sym[idx], src_pe);
    }
}

// Allokation im symmetrischen Heap (kollektiv auf allen PEs)
double* a = (double*) nvshmem_malloc(nx * ny * sizeof(double));

// Im Kernel: Rand direkt in den Heap des Nachbarn schreiben
__global__ void jacobi_with_halo(double* a_new, double* a,
                                  int iy_start, int iy_end, int nx,
                                  int top_pe, int btm_pe) {
    // ... Jacobi-Berechnung ...

    // Halo-Austausch GPU-seitig (nur ein Thread pro Reihe nötig)
    if (threadIdx.x == 0 && iy == iy_start) {
        // Sende oberste Zeile an top_pe
        nvshmem_double_put(a_new,                       // Ziel: Anfang von a_new auf top_pe
                           a_new + iy_start * nx,       // Quelle: lokale erste Zeile
                           nx, top_pe);
    }
    nvshmem_quiet();  // warte bis alle Puts fertig
}