我正在尝试使用 Blac 的 pdgemm 执行矩阵乘法。我正在使用的矩阵乘法的确切子程序可以在这里找到：http: //www.netlib.org/scalapack/html/pblas_qref.html#PvGEMM
但是，我的代码在 pdgemm 调用上返回“无法为线程本地数据分配内存：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 没有太多经验，但是从我看过的示例中，我不确定为什么会出现分段错误，我们将不胜感激。