mpi-matrix-inv.c

/*
 * MIT License
 *
 * Copyright (c) 2018 Paderborn Center for Parallel Computing
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */

#include <fcntl.h>
#include <mkl.h>
#include <mpi.h>
#include <omp.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/mman.h>
#include <unistd.h>

#define PATHLEN 255

struct properties {
  int size;
  int density;
  int condition;
};

MKL_INT find_elem(MKL_INT needle, MKL_INT *haystack, MKL_INT size) {
  MKL_INT l, r, c;
  l = 0;
  r = size;
  do {
    c = (l + r) / 2;
    if (haystack[c] == needle) {
      return c;
    }
    if (haystack[c] < needle) {
      l = c + 1;
    } else {
      r = c;
    }
  } while (l != r);
  return -1;
}

lapack_int invert_matrix(double *matrix, lapack_int size) {
  // First we need to compute the LU factorization using ?getrf
  lapack_int *ipiv, ret;
  ipiv = (lapack_int *)mkl_calloc(size, sizeof(lapack_int), 64);

  ret = LAPACKE_dgetrf(LAPACK_COL_MAJOR, size, size, matrix, size, ipiv);
  if (ret) {
    mkl_free(ipiv);
    return ret;
  }

  // And now we calculate the inverse using the LU factorization
  ret = LAPACKE_dgetri(LAPACK_COL_MAJOR, size, matrix, size, ipiv);
  mkl_free(ipiv);
  return ret;
}

void print_matrix(double *matrix, MKL_INT size) {
  MKL_INT i, j;
  for (i = 0; i < size; i++) {
    for (j = 0; j < size; j++) {
      printf("%.2f\t", matrix[i * size + j]);
    }
    printf("\n");
  }
  printf("\n");
}

void invert_submatrix(double *values, MKL_INT *row_ind, MKL_INT *col_ptr,
                      double *values_inv, int i, double *locDurBuild,
                      double *locDurCalc) {

  MKL_INT nnz, k, l, kcal, lcal, idx;
  lapack_int ret;
  double *submatrix;
  double tStart, tEnd;

  // 代表了第i列非零元的数量，同样就等于第i个子矩阵的维数
  nnz = col_ptr[i + 1] - col_ptr[i];
  submatrix = (double *)mkl_calloc(nnz * nnz, sizeof(double), 64);

  tStart = omp_get_wtime();
  for (k = 0; k < nnz; k++) {
    for (l = 0; l < nnz; l++) {
      // kcal、lcal是原始矩阵的行列坐标
      kcal = row_ind[col_ptr[i] + k];
      lcal = row_ind[col_ptr[i] + l];
      // We now have to copy M[kcal][lcal] to submatrix[k][l]
      // How to access M[kcal][lcal]? Calculate idx
      // find_elem是从原始矩阵的第lcal列中找到第kcal行的非零值，由于只有CSC格式，所以得用二分查找
      idx = find_elem(kcal, &(row_ind[col_ptr[lcal]]),
                      col_ptr[lcal + 1] - col_ptr[lcal]);
      if (idx != -1) {
        submatrix[k * nnz + l] = values[col_ptr[lcal] + idx];
      }
    }
  }
  tEnd = omp_get_wtime();
  *locDurBuild = (tEnd - tStart);

  tStart = omp_get_wtime();
  ret = invert_matrix(submatrix, nnz);
  tEnd = omp_get_wtime();
  *locDurCalc = (tEnd - tStart);
  if (ret) {
    fprintf(stderr, "Inverting submatrix failed\n");
  }

  // 大矩阵和子矩阵都是对称矩阵，所以读取某一行就是在读取对应列
  // tStart = omp_get_wtime();
  memcpy(values_inv,
         &(submatrix[find_elem(i, &(row_ind[col_ptr[i]]), nnz) * nnz]),
         nnz * sizeof(double));
  // tEnd = omp_get_wtime();
  // *locDurBuild += (tEnd - tStart);

  mkl_free(submatrix);
}

int main(int argc, char *argv[]) {

  int threadsupport;
  MPI_Init_thread(NULL, NULL, MPI_THREAD_FUNNELED, &threadsupport);
  if (threadsupport < MPI_THREAD_FUNNELED) {
    fprintf(stderr, "Could not initialize thread support.");
    omp_set_num_threads(1);
  }

  struct properties prop;
  int fd, world_rank, world_size, *displs, *recvcounts, mkl_threads;
  char fn_in_val[PATHLEN], fn_in_ri[PATHLEN], fn_in_cp[PATHLEN],
      fn_out_val[PATHLEN];
  MKL_INT *col_ptr, *row_ind, total_nnz, i, submatrices_per_worker, total_elem,
      my_first_col, next_first_col, submatrices_for_me;
  double *values, *values_inv, tStart, tEnd;
  FILE *fp;

  MPI_Comm_rank(MPI_COMM_WORLD, &world_rank);
  MPI_Comm_size(MPI_COMM_WORLD, &world_size);

  // printf("%d: I'm alive\n", world_rank);

  if (world_rank == 0) {

    /**************
     * MPI Rank 0 *
     **************/

    // We are the boot process. Decide on what to do, broadcast and gather.
    if (argc != 4) {
      fprintf(
          stderr,
          "%d: Main process needs to be called with parameters size density "
          "condition\n",
          world_rank);

      // printf("%d: Shutting down workers...\n", world_rank);
      prop.size = 0;
      MPI_Bcast(&prop, 3, MPI_INT, 0, MPI_COMM_WORLD);

      exit(EXIT_FAILURE);
    }
    prop.size = strtol(argv[1], NULL, 10);
    prop.density = strtol(argv[2], NULL, 10);
    prop.condition = strtol(argv[3], NULL, 10);

    submatrices_per_worker = prop.size / (world_size - 1);
    printf("%d: Each of the %d workers will solve %d submatrices.\n",
           world_rank, (world_size - 1), submatrices_per_worker);
    if (prop.size % (world_size - 1) != 0) {
      fprintf(stderr,
              "%d: WARNING: Load imbalanced. Last worker will have to "
              "solve %d additional submatrices\n",
              world_rank, prop.size % (world_size - 1));
    }

    /* Main evaluation loop */

    int evalRep, evalChoice;

    for (evalRep = 0; evalRep < 5; evalRep++) {
      for (evalChoice = 1; evalChoice > 0; evalChoice--) {

        snprintf(fn_in_cp, PATHLEN, "sprandsym-s%d-d%d-c%d-n%d.cp", prop.size,
                 prop.density, prop.condition, evalChoice);
        snprintf(fn_out_val, PATHLEN, "sprandsym-s%d-d%d-c%d-n%d.inv.val",
                 prop.size, prop.density, prop.condition, evalChoice);
        snprintf(fn_in_val, PATHLEN, "sprandsym-s%d-d%d-c%d-n%d.val", prop.size,
                 prop.density, prop.condition, evalChoice);
        snprintf(fn_in_ri, PATHLEN, "sprandsym-s%d-d%d-c%d-n%d.ri", prop.size,
                 prop.density, prop.condition, evalChoice);
        snprintf(fn_out_val, PATHLEN, "sprandsym-s%d-d%d-c%d-n%d.out.val",
                 prop.size, prop.density, prop.condition, evalChoice);

        // col_ptr最后一个元素存的是总非0元素数量，这样就能计算出最后一列的非0元数量
        col_ptr = (MKL_INT *)calloc(prop.size + 1, sizeof(MKL_INT));
        fp = fopen(fn_in_cp, "rb");
        fread(col_ptr, sizeof(MKL_INT), prop.size + 1, fp);
        fclose(fp);

        total_nnz = col_ptr[prop.size];
        row_ind = (MKL_INT *)calloc(total_nnz, sizeof(MKL_INT));
        fp = fopen(fn_in_ri, "rb");
        fread(row_ind, sizeof(MKL_INT), total_nnz, fp);
        fclose(fp);

        values = (double *)calloc(total_nnz, sizeof(double));
        fp = fopen(fn_in_val, "rb");
        fread(values, sizeof(double), total_nnz, fp);
        fclose(fp);
        /*
         fp = fopen(fn_out_val, "wb");
         fseek(fp, total_nnz*sizeof(double)-1, SEEK_SET);
         fputc('\0', fp);
         fclose(fp);
         fd = open(fn_out_val, O_RDWR);

         values_inv = (double*) mmap(NULL, total_nnz*sizeof(double),
                                     PROT_READ|PROT_WRITE, MAP_SHARED, fd, 0);
         close(fd);
        */
        values_inv = (double *)calloc(total_nnz, sizeof(double));

        displs = (int *)calloc(world_size, sizeof(int));
        recvcounts = (int *)calloc(world_size, sizeof(int));
        for (i = 1; i < world_size; i++) {
          displs[i] = col_ptr[(i - 1) * submatrices_per_worker];
          recvcounts[i] = col_ptr[i * submatrices_per_worker] -
                          col_ptr[(i - 1) * submatrices_per_worker];
        }
        // Allow last worker to send all remaining results
        recvcounts[world_size - 1] =
            col_ptr[prop.size] -
            col_ptr[(world_size - 2) * submatrices_per_worker];

        tStart = MPI_Wtime();
        // printf("%d: Broadcasting information to all workers...\n",
        // world_rank);
        MPI_Bcast(&prop, 3, MPI_INT, 0, MPI_COMM_WORLD);
        // printf("%d: ... done\n", world_rank);

#ifndef USE_BEEGFS
        // Send data to all workers
        MPI_Bcast(col_ptr, prop.size + 1, MPI_INT, 0, MPI_COMM_WORLD);
        MPI_Bcast(row_ind, total_nnz, MPI_INT, 0, MPI_COMM_WORLD);
        MPI_Bcast(values, total_nnz, MPI_DOUBLE, 0, MPI_COMM_WORLD);
#endif
        tEnd = MPI_Wtime();

        printf("%d: Wall time elapsed for Bcast: %dms\n", world_rank,
               (int)((tEnd - tStart) * 1000));

        tStart = MPI_Wtime();
        // printf("%d: Waiting for results...\n", world_rank);
        MPI_Gatherv(NULL, 0, MPI_DOUBLE, values_inv, recvcounts, displs,
                    MPI_DOUBLE, 0, MPI_COMM_WORLD);
        // printf("%d: ... done\n", world_rank);
        tEnd = MPI_Wtime();

        printf("%d: Wall time elapsed for Gatherv: %dms\n", world_rank,
               (int)((tEnd - tStart) * 1000));

        fp = fopen(fn_out_val, "wb");
        fwrite(values_inv, sizeof(double), total_nnz, fp);
        fclose(fp);

        free(row_ind);
        free(values);
        free(recvcounts);
        free(displs);
        free(col_ptr);
        free(values_inv);
        // munmap(values_inv, total_nnz*sizeof(double));
      }
    }

    /* End of main evaluation loop */

    // printf("%d: Shutting down workers...\n", world_rank);
    prop.size = 0;
    MPI_Bcast(&prop, 3, MPI_INT, 0, MPI_COMM_WORLD);

  } else {

    /***************
     * Worker code *
     ***************/

    // We are one of the workers. Run in a loop and wait for jobs.
    while (1) {
      printf("%d: Waiting for matrix properties...\n", world_rank);
      MPI_Bcast(&prop, 3, MPI_INT, 0, MPI_COMM_WORLD);
      printf("%d: ... received\n", world_rank);
      if (prop.size == 0) {
        printf("%d: Received signal to halt.\n", world_rank);
        break;
      }

      snprintf(fn_in_cp, PATHLEN, "sprandsym-s%d-d%d-c%d-n1.cp", prop.size,
               prop.density, prop.condition);
      snprintf(fn_out_val, PATHLEN, "sprandsym-s%d-d%d-c%d-n1.inv.val",
               prop.size, prop.density, prop.condition);
      snprintf(fn_in_val, PATHLEN, "sprandsym-s%d-d%d-c%d-n1.val", prop.size,
               prop.density, prop.condition);
      snprintf(fn_in_ri, PATHLEN, "sprandsym-s%d-d%d-c%d-n1.ri", prop.size,
               prop.density, prop.condition);

#ifdef USE_BEEGFS
#ifdef USE_MMAP
      fd = open(fn_in_cp, O_RDONLY);
      col_ptr = (MKL_INT *)mmap(NULL, (prop.size + 1) * sizeof(MKL_INT),
                                PROT_READ, MAP_SHARED, fd, 0);
      close(fd);
#else
      fp = fopen(fn_in_cp, "rb");
      col_ptr = (MKL_INT *)calloc(prop.size + 1, sizeof(MKL_INT));
      fread(col_ptr, sizeof(MKL_INT), prop.size + 1, fp);
      fclose(fp);
#endif
#else
      col_ptr = (MKL_INT *)calloc(prop.size + 1, sizeof(MKL_INT));
      MPI_Bcast(col_ptr, prop.size + 1, MPI_INT, 0, MPI_COMM_WORLD);
#endif

      total_nnz = col_ptr[prop.size];

#ifdef USE_BEEGFS
#ifdef USE_MMAP
      fd = open(fn_in_ri, O_RDONLY);
      row_ind = (MKL_INT *)mmap(NULL, total_nnz * sizeof(MKL_INT), PROT_READ,
                                MAP_SHARED, fd, 0);
      close(fd);
      fd = open(fn_in_val, O_RDONLY);
      values = (double *)mmap(NULL, total_nnz * sizeof(double), PROT_READ,
                              MAP_SHARED, fd, 0);
      close(fd);
#else
      row_ind = (MKL_INT *)calloc(total_nnz, sizeof(MKL_INT));
      values = (double *)calloc(total_nnz, sizeof(double));
      fp = fopen(fn_in_ri, "rb");
      fread(row_ind, sizeof(MKL_INT), total_nnz, fp);
      fclose(fp);
      fp = fopen(fn_in_val, "rb");
      fread(values, sizeof(double), total_nnz, fp);
      fclose(fp);
#endif
#else
      row_ind = (MKL_INT *)calloc(total_nnz, sizeof(MKL_INT));
      values = (double *)calloc(total_nnz, sizeof(double));
      MPI_Bcast(row_ind, total_nnz, MPI_INT, 0, MPI_COMM_WORLD);
      MPI_Bcast(values, total_nnz, MPI_DOUBLE, 0, MPI_COMM_WORLD);
#endif

      submatrices_per_worker = prop.size / (world_size - 1);
      my_first_col = (world_rank - 1) * submatrices_per_worker;
      if (world_rank == world_size - 1) {
        // We are the last worker. Maybe we have to do additional work...
        submatrices_for_me =
            submatrices_per_worker + (prop.size % (world_size - 1));
        next_first_col = prop.size;
      } else {
        submatrices_for_me = submatrices_per_worker;
        next_first_col = world_rank * submatrices_per_worker;
      }
      total_elem = col_ptr[next_first_col] - col_ptr[my_first_col];
      values_inv = (double *)calloc(total_elem, sizeof(double));

      /* Optimize threading: We should do as much submatrices as possible in
       * parallel. If threads are left, leave them for MKL's internal
       * parallelism. */
      mkl_threads = omp_get_max_threads() / submatrices_for_me;
      if (mkl_threads < 1) {
        mkl_threads = 1;
      }
      mkl_set_num_threads(mkl_threads);

      printf("%d: We have %d thread(s) to solve %d submatrices. Give %d "
             "thread(s) to MKL for each submatrix operation.\n",
             world_rank, omp_get_max_threads(), submatrices_for_me,
             mkl_threads);

      double durationBuild = .0;
      double durationCalc = .0;
      tStart = MPI_Wtime();
// printf("%d: Starting the number crunching\n", world_rank);
#pragma omp parallel for schedule(dynamic)                                     \
    reduction(+ : durationBuild, durationCalc)
      for (i = 0; i < submatrices_for_me; i++) {
        double locDurBuild, locDurCalc;
        // printf("%d: Inverting submatrix %d in thread %d.\n", world_rank,
        //       (world_rank - 1) * submatrices_for_me + i,
        //       omp_get_thread_num());
        invert_submatrix(
            values, row_ind, col_ptr,
            &(values_inv[col_ptr[my_first_col + i] - col_ptr[my_first_col]]),
            my_first_col + i, &locDurBuild, &locDurCalc);
        durationBuild += locDurBuild;
        durationCalc += locDurCalc;
      }
      tEnd = MPI_Wtime();

      printf("%d: Wall time elapsed: %dms\n", world_rank,
             (int)((tEnd - tStart) * 1000));
      printf("%d: CPU time sm build: %dms\n", world_rank,
             (int)(durationBuild * 1000));
      printf("%d: CPU time sm calc: %dms\n", world_rank,
             (int)(durationCalc * 1000));

      // printf("%d: Send results to root\n", world_rank);
      MPI_Gatherv(values_inv, total_elem, MPI_DOUBLE, NULL, NULL, NULL,
                  MPI_DOUBLE, 0, MPI_COMM_WORLD);
      // printf("%d: ... done\n", world_rank);

      memset(values_inv, 0, total_elem * sizeof(double));
      free(values_inv);
#if defined USE_BEEGFS && defined USE_MMAP
      munmap(values, total_nnz * sizeof(double));
      munmap(row_ind, total_nnz * sizeof(MKL_INT));
      munmap(col_ptr, (prop.size + 1) * sizeof(MKL_INT));
#else
      memset(values, 0, total_nnz * sizeof(double));
      memset(row_ind, 0, total_nnz * sizeof(MKL_INT));
      memset(col_ptr, 0, (prop.size + 1) * sizeof(MKL_INT));
      free(values);
      free(row_ind);
      free(col_ptr);
#endif
    }
  }

  MPI_Barrier(MPI_COMM_WORLD);
  MPI_Finalize();
  exit(EXIT_SUCCESS);
}