在C中使用mpi和lapack时出现分段错误

Question

我正在尝试使用Blac的pdgemm来实现矩阵乘法。我正在使用的矩阵乘法的确切的子例程可以在这里找到：http://www.netlib.org/scalapack/html/pblas_qref.html#PvGEMM

但是，我的代码在pdgemm调用时返回"cannot allocate memory for thread-local data: ABORT“。在我的代码中，我将一个矩阵乘以它自己，它是一个方阵，所以得到的矩阵是相同的维数。下面是有问题的代码

#include <stdlib.h>
#include <stdio.h>
#include "mpi.h"
#include <math.h>

#define gridSize 10

int main(int argc, char* argv[]) {
  int i,j,k, np, myid;
  int bufsize;
  double *buf;          /* Buffer for reading */

  MPI_Offset filesize;
  MPI_File myfile;    /* Shared file */
  MPI_Status status;  /* Status returned from read */

  /* Initialize MPI */
  MPI_Init(&argc, &argv);
  MPI_Comm_rank(MPI_COMM_WORLD, &myid);
  MPI_Comm_size(MPI_COMM_WORLD, &np);

  double *A = (double*) malloc(gridSize*gridSize*sizeof(double));

  /*MPI-IO Code removed for clarity including buf and bufsize assignments
  .
  .
  .
  .
  */
    //I use mpi-IO to read in a matrix from a file, each processor reads in a row and that row is store in the array called buf
    for (j = 0; j <bufsize;j++){
    A[myid*bufsize+j] = buf[j];
  }

  //BLACS portion
  int ictxt, nprow, npcol, myrow, mycol, nb;
  int info,itemp;
  int ZERO = 0, ONE = 1;
  nprow = 2; npcol = 2; nb = 2;

  Cblacs_pinfo(&myid,&np);
  Cblacs_get(-1,0,&ictxt);
  Cblacs_gridinit(&ictxt,"Row",nprow,npcol);
  Cblacs_gridinfo(ictxt,&nprow,&npcol,&myrow,&mycol);

  int M = gridSize;

  int descA[9], descx[9], descy[9];
  int mA = numroc_(&M, &nb, &myrow, &ZERO, &nprow);
  int nA = numroc_(&M, &nb, &mycol, &ZERO, &npcol);
  int nx = numroc_(&M, &nb, &myrow, &ZERO, &nprow);
  int my = numroc_(&M ,&nb, &myrow, &ZERO, &nprow);

  descinit_(descA,&M,&M,&nb,&nb,&ZERO,&ZERO,&ictxt,&mA,&info);
  descinit_(descx,&M,&M,&nb,&nb,&ZERO,&ZERO,&ictxt,&nx,&info);

  descinit_(descy,&M,&M,&nb,&nb,&ZERO,&ZERO,&ictxt,&my,&info);


  double* matrixA = (double*)malloc(mA*nA*sizeof(double));
  double* matrixB = (double*)malloc(mA*nA*sizeof(double));
  double* matrixC = (double*)calloc(mA*nA,sizeof(double));
  int sat,sut;


  for(i=0;i<mA;i++){
    for(j=0;j<nA;j++){
        sat = (myrow*nb)+i+(i/nb)*nb;
        sut = (mycol*nb)+j+(j/nb)*nb;
        matrixA[j*mA+i] = A[sat*M+sut];
        matrixB[j*mA+i] = A[sat*M+sut];
    }
  }

  double alpha = 1.0; double beta = 0.0;

  //call where seg fault happens
  pdgemm_("N","N",&M,&M,&M,&alpha,matrixA,&ONE,&ONE,descA,matrixB,&ONE,&ONE,descx,
      &beta,matrixC,&ONE,&ONE,descy);

  Cblacs_barrier(ictxt,"A");

  Cblacs_gridexit(0);

  /* Close the file */
  MPI_File_close(&myfile);

  if (myid==0) {
    printf("Done\n");

  }
  MPI_Finalize();

  exit(0);

}

我没有太多使用ScaLapacs的经验，但从我看过的例子中，我不确定为什么会出现分段错误，任何帮助都将不胜感激。

Answer 1

我正确地把"int“改成了"long”，把"double“改成了"long double”

我尝试过的另一种有效的方法是静态链接库。

mpicc -w -o a a.c -L$MKLPATH -I$IMKLPATH -Wl，--开始组$MKLPATH/libmkl_scalapack.a $MKLPATH/libmkl_blacs_openmpi_lp64.a $MKLPATH/libmkl_intel_lp64.a $MKLPATH/libmkl_intel_thread.a $MKLPATH/libmkl_ -static_mpi -Wl，--结束组-lpthread -lm -openmp