dbm_miniapp.c

Go to the documentation of this file.

/*----------------------------------------------------------------------------*/
/*  CP2K: A general program to perform molecular dynamics simulations         */
/*  Copyright 2000-2026 CP2K developers group <https://cp2k.org>              */
/*                                                                            */
/*  SPDX-License-Identifier: BSD-3-Clause                                     */
/*----------------------------------------------------------------------------*/
#include "dbm_library.h"
#include "dbm_matrix.h"
 
#include "../mpiwrap/cp_mpi.h"
#include "../offload/offload_library.h"
#include "../offload/offload_mempool.h"
 
#include <assert.h>
#include <limits.h>
#include <omp.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
 
/*******************************************************************************
 * \brief Wrapper for printf, passed to dbm_library_print_stats.
 * \author Ole Schuett
 ******************************************************************************/
static void print_func(const char *msg, int msglen, int output_unit) {
  (void)msglen;           // mark used
  if (output_unit == 0) { // i.e. my_rank == 0
    printf("%s", msg);
  }
}
 
/*******************************************************************************
 * \brief Returns the smaller of the two integers (missing from the C standard).
 * \author Ole Schuett
 ******************************************************************************/
static inline int imin(int x, int y) { return (x < y ? x : y); }
 
/*******************************************************************************
 * \brief Private routine for creating a distribution.
 * \author Ole Schuett
 ******************************************************************************/
static dbm_distribution_t *create_dist(const int nrows, const int ncols,
                                       const cp_mpi_comm_t comm) {
  int cart_dims[2], cart_periods[2], cart_coords[2];
  cp_mpi_cart_get(comm, 2, cart_dims, cart_periods, cart_coords);
 
  // Create distribution.
  assert(0 < nrows && 0 < ncols);
  int *row_dist = malloc(nrows * sizeof(int));
  int *col_dist = malloc(ncols * sizeof(int));
  assert(row_dist != NULL && col_dist != NULL);
  for (int i = 0; i < nrows; i++) {
    row_dist[i] = i % cart_dims[0];
  }
  for (int i = 0; i < ncols; i++) {
    col_dist[i] = i % cart_dims[1];
  }
  const int fortran_comm = cp_mpi_comm_c2f(comm);
  dbm_distribution_t *dist = NULL;
  dbm_distribution_new(&dist, fortran_comm, nrows, ncols, row_dist, col_dist);
  free(row_dist);
  free(col_dist);
  return dist;
}
 
/*******************************************************************************
 * \brief Private routine for drawing block sizes from [size_min, size_max].
 * \author Ole Schuett
 ******************************************************************************/
static void set_some_sizes(const int nsizes, const int size_min,
                           const int size_max, int sizes[nsizes]) {
  assert(0 < size_min && size_min <= size_max);
  const int range = size_max - size_min + 1;
  for (int i = 0; i < nsizes; i++) {
    sizes[i] = (1 < range ? size_min + rand() % range : size_max);
  }
}
 
/*******************************************************************************
 * \brief Private routine for creating a distribution and an empty matrix.
 * \author Ole Schuett
 ******************************************************************************/
static dbm_matrix_t *
create_some_matrix(const int nrows, const int ncols, const int nrows_min,
                   const int nrows_max, const int ncols_min,
                   const int ncols_max, const cp_mpi_comm_t comm) {
  // Create distribution.
  dbm_distribution_t *dist = create_dist(nrows, ncols, comm);
 
  // Create matrix.
  assert(0 < nrows && 0 < ncols);
  int *row_sizes = malloc(nrows * sizeof(int));
  int *col_sizes = malloc(ncols * sizeof(int));
  assert(row_sizes != NULL && col_sizes != NULL);
  set_some_sizes(nrows, nrows_min, nrows_max, row_sizes);
  set_some_sizes(ncols, ncols_min, ncols_max, col_sizes);
  dbm_matrix_t *matrix = NULL;
  dbm_create(&matrix, dist, "some name", nrows, ncols, row_sizes, col_sizes);
  dbm_distribution_release(dist);
  free(row_sizes);
  free(col_sizes);
  return matrix;
}
 
/*******************************************************************************
 * \brief Private routine for reserving all blocks of the given matrix.
 * \author Ole Schuett
 ******************************************************************************/
static void reserve_all_blocks(dbm_matrix_t *matrix) {
  int nrows, ncols;
  const int *row_sizes, *col_sizes;
  dbm_get_row_sizes(matrix, &nrows, &row_sizes);
  dbm_get_col_sizes(matrix, &ncols, &col_sizes);
 
#pragma omp parallel
  {
    int nblocks = 0;
#pragma omp for collapse(2)
    for (int row = 0; row < nrows; row++) {
      for (int col = 0; col < ncols; col++) {
        if (dbm_get_stored_coordinates(matrix, row, col) ==
            matrix->dist->my_rank) {
          ++nblocks;
        }
      }
    }
    int *reserve_row = NULL, *reserve_col = NULL;
    if (0 < nblocks) {
      reserve_row = malloc(nblocks * sizeof(int));
      reserve_col = malloc(nblocks * sizeof(int));
      assert(reserve_row != NULL && reserve_col != NULL);
    }
    int iblock = 0;
#pragma omp for collapse(2)
    for (int row = 0; row < nrows; row++) {
      for (int col = 0; col < ncols; col++) {
        if (dbm_get_stored_coordinates(matrix, row, col) ==
            matrix->dist->my_rank) {
          reserve_row[iblock] = row;
          reserve_col[iblock] = col;
          ++iblock;
        }
      }
    }
    assert(iblock == nblocks);
    dbm_reserve_blocks(matrix, nblocks, reserve_row, reserve_col);
    free(reserve_row);
    free(reserve_col);
  }
}
 
/*******************************************************************************
 * \brief Private routine for setting all blocks.
 * \author Ole Schuett
 ******************************************************************************/
static void set_all_blocks(dbm_matrix_t *matrix) {
#pragma omp parallel
  {
    dbm_iterator_t *iter = NULL;
    dbm_iterator_start(&iter, matrix);
    while (dbm_iterator_blocks_left(iter)) {
      int row, col, row_size, col_size;
      double *block;
      dbm_iterator_next_block(iter, &row, &col, &block, &row_size, &col_size);
      // Blocks of the same size differ, which exposes mixed up blocks.
      const int shift = (31 * row + 17 * col) % 101;
      const int block_size = row_size * col_size;
      for (int i = 0; i < block_size; i++) {
        block[i] = 1.0 / (i + 1 + shift);
      }
    }
    dbm_iterator_stop(iter);
  }
}
/*******************************************************************************
 * \brief Run a benchmark of dbm_multiply with given block sizes. The blocks
 *        of A have 1 to m rows and those of B have 1 to n columns, unless
 *        DBM_MINIAPP_FIXED is set (non-zero), e.g., for homogeneous batches.
 * \author Ole Schuett
 ******************************************************************************/
void benchmark_multiply(const int M, const int N, const int K, const int m,
                        const int n, const int k, const cp_mpi_comm_t comm) {
  const char *const fixed_env = getenv("DBM_MINIAPP_FIXED");
  const bool fixed = (NULL != fixed_env && 0 != atoi(fixed_env));
 
  srand(25071975); // same matrices regardless of preceding benchmarks
  dbm_matrix_t *matrix_a =
      create_some_matrix(M, K, fixed ? m : 1, m, k, k, comm);
  dbm_matrix_t *matrix_b =
      create_some_matrix(K, N, k, k, fixed ? n : 1, n, comm);
  dbm_distribution_t *dist_c = create_dist(M, N, comm);
  dbm_matrix_t *matrix_c = NULL, *matrix_d = NULL;
  dbm_create(&matrix_c, dist_c, "result", M, N, matrix_a->row_sizes,
             matrix_b->col_sizes);
  dbm_distribution_release(dist_c);
 
  reserve_all_blocks(matrix_a);
  reserve_all_blocks(matrix_b);
  set_all_blocks(matrix_a);
  set_all_blocks(matrix_b);
 
  const char *const verify_env = getenv("DBM_MULTIPLY_VERIFY");
  const int skip_verify = (NULL == verify_env ? 0 : (atoi(verify_env) + 1));
 
  if (0 == skip_verify) {
    dbm_distribution_t *const dist_shared = matrix_c->dist;
    dbm_create(&matrix_d, dist_shared, matrix_c->name, matrix_c->nrows,
               matrix_c->ncols, matrix_c->row_sizes, matrix_c->col_sizes);
    dbm_copy(matrix_d, matrix_c);
  }
 
  int64_t flop = 0;
  const double time_start_multiply = omp_get_wtime();
  dbm_multiply(false, false, 1.0, matrix_a, matrix_b, 1.0, matrix_c, false,
               1e-8, &flop);
  const double time_end_multiply = omp_get_wtime();
 
  if (cp_mpi_comm_rank(comm) == 0) {
    printf("%5i x %5i x %5i  with  %3i x %3i x %3i blocks: ", M, N, K, m, n, k);
  }
 
  if (NULL != matrix_d) { // Calculate result on the host for validation.
    dbm_multiply(false, false, 1.0, matrix_a, matrix_b, 1.0, matrix_d, false,
                 1e-8, NULL);
 
    const double maxeps = 1E-5, epsilon = dbm_maxeps(matrix_d, matrix_c);
    if (maxeps < epsilon) {
      printf("ERROR\n");
      fprintf(stderr, "Failed validation (epsilon=%f).\n", epsilon);
      exit(1);
    }
    dbm_release(matrix_d);
  }
 
  cp_mpi_sum_int64(&flop, 1, comm);
  if (cp_mpi_comm_rank(comm) == 0) {
    const double duration = time_end_multiply - time_start_multiply;
    printf("%6.3f s =>  %6.1f GFLOP/s\n", duration, 1e-9 * flop / duration);
    fflush(stdout);
  }
 
  dbm_release(matrix_a);
  dbm_release(matrix_b);
  dbm_release(matrix_c);
}
 
/*******************************************************************************
 * \brief Private routine for parsing a triplet like MxNxK, where omitted
 *        components default to the first one. Returns false if invalid.
 * \author Hans Pabst
 ******************************************************************************/
static bool parse_triplet(const char *arg, int triplet[3]) {
  const char *const delims = "x,;:|/";
  int ntriplet = 0;
  while (ntriplet < 3 && '\0' != *arg) {
    char *end = NULL;
    const long value = strtol(arg, &end, 10);
    if (end == arg || value <= 0 || INT_MAX < value) {
      break;
    }
    triplet[ntriplet++] = (int)value;
    for (arg = end; '\0' != *arg && NULL != strchr(delims, *arg); ++arg) {
    }
  }
  for (int i = ntriplet; 0 < ntriplet && i < 3; i++) {
    triplet[i] = triplet[0];
  }
  return (0 < ntriplet && '\0' == *arg);
}
 
/*******************************************************************************
 * \brief Private routine for running benchmarks given as pairs of triplets,
 *        i.e., block sizes followed by the matrix size in blocks (default:
 *        128x128x128 if the last pair is incomplete).
 * \author Hans Pabst
 ******************************************************************************/
static bool run_benchmarks(const int nargs, char *args[],
                           const cp_mpi_comm_t comm) {
  bool valid = (0 < nargs);
  for (int i = 0; valid && i < nargs; i += 2) {
    int mnk[3], MNK[3] = {128, 128, 128};
    valid = parse_triplet(args[i], mnk) &&
            (nargs <= i + 1 || parse_triplet(args[i + 1], MNK));
    if (valid) {
      benchmark_multiply(MNK[0], MNK[1], MNK[2], mnk[0], mnk[1], mnk[2], comm);
    }
  }
  return valid;
}
 
/*******************************************************************************
 * \brief Stand-alone miniapp for smoke-testing and benchmarking dbm_multiply.
 * \author Ole Schuett
 ******************************************************************************/
int main(int argc, char *argv[]) {
  int result = EXIT_SUCCESS;
 
  cp_mpi_init(&argc, &argv);
  dbm_library_init();
 
  const cp_mpi_comm_t world_comm = cp_mpi_get_comm_world();
  const int nranks = cp_mpi_comm_size(world_comm);
  const int my_rank = cp_mpi_comm_rank(world_comm);
 
  if (offload_get_device_count() > 0) {
    offload_set_chosen_device(my_rank % offload_get_device_count());
  }
 
  // Create 2D cart.
  int dims[2] = {0, 0};
  cp_mpi_dims_create(nranks, 2, dims);
  const int periods[2] = {true, true};
  cp_mpi_comm_t comm = cp_mpi_cart_create(world_comm, 2, dims, periods, false);
 
  if (my_rank == 0) {
    printf("OpenMP-threads: %i  GPUs: %i", omp_get_max_threads(),
           imin(offload_get_device_count(), nranks));
#if defined(__parallel)
    printf("  MPI-ranks: %i  MPI-cart: %i x %i", nranks, dims[0], dims[1]);
#else
    printf("  MPI: n/a");
#endif
    printf("\n\n");
    fflush(stdout);
  }
 
  if (1 >= argc) {
    benchmark_multiply(16384, 128, 128, 4, 4, 4, comm);
    benchmark_multiply(128, 16384, 128, 4, 4, 4, comm);
    benchmark_multiply(128, 128, 16384, 4, 4, 4, comm);
    benchmark_multiply(645, 645, 645, 4, 4, 4, comm);
    if (my_rank == 0)
      printf("\n");
 
    benchmark_multiply(60, 500, 500, 128, 4, 4, comm);
    benchmark_multiply(500, 60, 500, 4, 128, 4, comm);
    benchmark_multiply(500, 500, 60, 4, 4, 128, comm);
    if (my_rank == 0)
      printf("\n");
 
    benchmark_multiply(500, 60, 60, 4, 128, 128, comm);
    benchmark_multiply(60, 500, 60, 128, 4, 128, comm);
    benchmark_multiply(60, 60, 500, 128, 128, 4, comm);
    if (my_rank == 0)
      printf("\n");
 
    benchmark_multiply(350, 350, 350, 23, 23, 23, comm);
    benchmark_multiply(250, 250, 250, 32, 32, 32, comm);
    benchmark_multiply(60, 60, 60, 128, 128, 128, comm);
  } else { // Pairs of triplets from the command line or per line of a file.
    FILE *const file = fopen(argv[1], "r");
    bool valid = true;
    if (NULL != file) {
      char buffer[1024];
      while (valid && NULL != fgets(buffer, sizeof(buffer), file)) {
        char *args[2] = {NULL, NULL};
        int nargs = 0;
        for (char *arg = strtok(buffer, " \t\r\n"); NULL != arg;
             arg = strtok(NULL, " \t\r\n")) {
          valid = (nargs < 2);
          if (!valid) {
            break;
          }
          args[nargs++] = arg;
        }
        if (valid && 0 < nargs) { // skip empty lines
          valid = run_benchmarks(nargs, args, comm);
        }
      }
      fclose(file);
    } else {
      valid = run_benchmarks(argc - 1, argv + 1, comm);
    }
    if (!valid) {
      if (my_rank == 0) {
        fprintf(stderr, "Usage: %s [file | mxnxk [MxNxK] ...]\n", argv[0]);
      }
      result = EXIT_FAILURE;
    }
  }
 
  if (EXIT_SUCCESS == result) {
    const int fortran_comm = cp_mpi_comm_c2f(comm);
    dbm_library_print_stats(fortran_comm, &print_func, my_rank);
    offload_mempool_stats_print(fortran_comm, &print_func, my_rank);
  }
  dbm_library_finalize();
  cp_mpi_comm_free(&comm);
  cp_mpi_finalize();
  return result;
}
 
// EOF