维度变大时出现矩阵计算错误

Question

我是运行我只是创建 2 个矩阵的代码：一个矩阵的维度为 arows x nsame，另一个矩阵的维度为 nsame x bcols。结果是维度数组 arows x bcols。这使用 BLAS 实现起来相当简单，当将以下主从模型与 OpenMPI 一起使用时，以下代码似乎按预期工作：`

#include <iostream>
#include <stdio.h>
#include <iostream>
#include <cmath>
#include <mpi.h>
#include <gsl/gsl_blas.h>
using namespace std;`

int main(int argc, char** argv){
    int noprocs, nid;
    MPI_Status status;
    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &nid);
    MPI_Comm_size(MPI_COMM_WORLD, &noprocs);
    int master = 0;

    const int nsame = 500; //must be same if matrices multiplied together = acols = brows
    const int arows = 500;
    const int bcols = 527; //works for 500 x 500 x 527 and 6000 x 100 x 36
    int rowsent;
    double buff[nsame];
    double b[nsame*bcols];
    double c[arows][bcols];
    double CC[1*bcols]; //here ncols corresponds to numbers of rows for matrix b
    for (int i = 0; i < bcols; i++){
                CC[i] = 0.;
    }; 
    // Master part
    if (nid == master ) { 

        double a [arows][nsame]; //creating identity matrix of dimensions arows x nsame (it is I if arows = nsame)
        for (int i = 0; i < arows; i++){
            for (int j = 0; j < nsame; j++){
                if (i == j)
                    a[i][j] = 1.;
                else
                    a[i][j] = 0.;
            }
        }
        double b[nsame*bcols];//here ncols corresponds to numbers of rows for matrix b
            for (int i = 0; i < (nsame*bcols); i++){
                b[i] = (10.*i + 3.)/(3.*i - 2.) ;
            }; 
        MPI_Bcast(b,nsame*bcols, MPI_DOUBLE_PRECISION, master, MPI_COMM_WORLD);  
        rowsent=0;
        for (int i=1; i < (noprocs); i++) {  
            // Note A is a 2D array so A[rowsent]=&A[rowsent][0]
            MPI_Send(a[rowsent], nsame, MPI_DOUBLE_PRECISION,i,rowsent+1,MPI_COMM_WORLD);
            rowsent++; 
        }

        for (int i=0; i<arows; i++) { 
            MPI_Recv(CC, bcols, MPI_DOUBLE_PRECISION, MPI_ANY_SOURCE, MPI_ANY_TAG,
                     MPI_COMM_WORLD, &status); 
            int sender = status.MPI_SOURCE;
            int anstype = status.MPI_TAG;            //row number+1
            int IND_I = 0;
            while (IND_I < bcols){
                c[anstype - 1][IND_I] = CC[IND_I]; 
                IND_I++;
            }
            if (rowsent < arows) {
                MPI_Send(a[rowsent], nsame,MPI_DOUBLE_PRECISION,sender,rowsent+1,MPI_COMM_WORLD);
                rowsent++; 
            }
            else {       // tell sender no more work to do via a 0 TAG
                MPI_Send(MPI_BOTTOM,0,MPI_DOUBLE_PRECISION,sender,0,MPI_COMM_WORLD);
            }
        }
    }

    // Slave part
    else { 
        MPI_Bcast(b,nsame*bcols, MPI_DOUBLE_PRECISION, master, MPI_COMM_WORLD); 
        MPI_Recv(buff,nsame,MPI_DOUBLE_PRECISION,master,MPI_ANY_TAG,MPI_COMM_WORLD,&status); 
        while(status.MPI_TAG != 0) {
            int crow = status.MPI_TAG; 
            gsl_matrix_view AAAA = gsl_matrix_view_array(buff, 1, nsame);
            gsl_matrix_view BBBB = gsl_matrix_view_array(b, nsame, bcols);
            gsl_matrix_view CCCC = gsl_matrix_view_array(CC, 1, bcols);

            /* Compute C = A B */
            gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, &AAAA.matrix, &BBBB.matrix,
                            0.0, &CCCC.matrix); 

            MPI_Send(CC,bcols,MPI_DOUBLE_PRECISION, master, crow, MPI_COMM_WORLD);
            MPI_Recv(buff,nsame,MPI_DOUBLE_PRECISION,master,MPI_ANY_TAG,MPI_COMM_WORLD,&status); 
        }
    }

    // output c here on master node //uncomment the below lines if I wish to see the output
    //    if (nid == master){
//        if (rowsent == arows){
//            //            cout << rowsent;
//            int IND_F = 0;
//            while (IND_F < arows){
//                int IND_K = 0;
//                while (IND_K < bcols){
//                    cout << "[" << IND_F << "]" << "[" << IND_K << "] = " << c[IND_F][IND_K] << " ";
//                    IND_K++;
//                }
//                cout << "\n";
//                IND_F++;
//            }
//        }
//    }
    MPI_Finalize();
    //free any allocated space here
    return 0;
};

现在看起来奇怪的是，当我增加矩阵的大小时（例如，从 nsame = 500 到 nsame = 501），代码不再有效。我收到以下错误：

mpirun noticed that process rank 0 with PID 0 on node Users-MacBook-Air exited on signal 11 (Segmentation fault: 11).

我已经尝试过使用矩阵的其他大小组合，并且矩阵本身的大小似乎总是有上限（这似乎根据我如何改变不同维度本身而有所不同）。我也尝试过修改矩阵本身的值，尽管这似乎没有任何改变。我意识到在我的示例中有其他方法可以初始化矩阵（例如使用向量），但我只是想知道为什么我当前的任意大小矩阵相乘方案似乎只能在一定程度上起作用。

Answer 1

您声明了过多的大局部变量，这会导致堆栈 space 相关问题。 a，尤其是 500x500 双精度数（250000 个 8 字节元素，或 200 万字节）。 b更大。

您需要为这些数组中的部分或全部动态分配 space。

可能有一个编译器选项可以增加初始堆栈 space 但这不是一个好的长期解决方案。

维度变大时出现矩阵计算错误

Matrix calculation error appears when dimensions become large

c++

arrays

matrix

openmpi

blas