为什么这个MEXed C / magma代码在独立的C代码工作时会出错？

Question

以下MEXed C代码只是调用magma来反转矩阵。独立的C代码（也发布）工作，但mex代码崩溃。

我已经三重检查了文档，验证了其他岩浆功能按预期工作，并且posted on the Magma forum and was told my code is fine（这篇文章是Magma论坛的一篇文章）。这意味着问题出在mex上。我想知道是什么导致mex代码出现seg-fault以及如何让它按预期运行。

Mexed代码：

#include <mex.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <stddef.h>
#include <magma_v2.h>
#include <cuda_runtime.h>

void mat2magma(magmaDoubleComplex* p, double* pr, double* pi,int numElements)
{
    int j=0;
    for(j=0;j<numElements;j++){
        p[j].x=pr[j];
        p[j].y=pi[j];
    }
}

void magma2mat(magmaDoubleComplex* p, double* pr, double* pi,int numElements)
{
    int j=0;
    for(j=0;j<numElements;j++){
        pr[j]= p[j].x;
        pi[j]= p[j].y;
    }
}

/*gateway function*/
void mexFunction( int nlhs, mxArray *plhs[],
        int nrhs, const mxArray *prhs[]) {

    /*initialize magma*/
    magma_init();
    magma_queue_t queue = NULL;
    magma_device_t dev;
    magma_getdevice(&dev);
    magma_queue_create(dev,&queue );

    magma_int_t m,ldwork,info;
    magma_int_t *piv;
    magmaDoubleComplex *a,*da,*dwork;

    /* Matlab -> Host */
    m=mxGetM(prhs[0]);
    piv=(magma_int_t*) malloc(m*sizeof(magma_int_t));
    magma_zmalloc_cpu(&a,m*m);
    mat2magma(a,mxGetPr(prhs[0]),mxGetPi(prhs[0]),m*m);
    ldwork = m*magma_get_zgetri_nb(m);

    /* Host -> GPU */
    magma_zmalloc(&dwork,ldwork);
    magma_zmalloc(&da,m*m);
    magma_zsetmatrix(m,m,a,m,da,m,queue);

    /*LU and Inverse */
    magma_zgetrf_gpu(m,m,da,m,piv,&info);
    magma_zgetri_gpu(m,da,m,piv,dwork,ldwork,&info);

    /*GPU -> Host */
    magma_zgetmatrix(m,m,da,m,a,m,queue);

    /*Host -> Matlab*/
    plhs[0] = mxCreateDoubleMatrix(m,m,mxCOMPLEX);
    magma2mat(a,mxGetPr(plhs[0]),mxGetPi(plhs[0]),m*m);
    free(a);
    free(piv);
    magma_free(dwork);
    magma_free(da);
    magma_queue_destroy(queue);
    magma_finalize();
}

我用mex CC=gcc LDFLAGS="-lmagma -lcudart -lcublas" magmaZinv.c然后从matlab编译它，我运行：

a=magic(3)+magic(3)*1i;
magmaZinv(a)

独立C代码：

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <stddef.h>
#include <magma_v2.h>
#include <cuda_runtime.h>
#include <sys/time.h>
#include <time.h>

/*gateway function*/
int main() {

    /*initialize magma*/
    magma_init();
    magma_queue_t queue = NULL;
    magma_device_t dev;
    magma_getdevice(&dev);
    magma_queue_create(dev,&queue );

    int m,ldwork,info;
    int *piv;
    magmaDoubleComplex *a,*da,*dwork;

    /* allocate and initialize a = magic(3)+magic(3)*1i; */
    m=3;
    piv=(int*) malloc(m*sizeof(int));
    ldwork = m*magma_get_zgetri_nb(m);
    magma_zmalloc_cpu(&a,m*m);
    a[0].x=8;a[0].y=8;
    a[1].x=3;a[1].y=3;
    a[2].x=4;a[2].y=4;
    a[3].x=1;a[3].y=1;
    a[4].x=5;a[4].y=5;
    a[5].x=9;a[5].y=9;
    a[6].x=6;a[6].y=6;
    a[7].x=7;a[7].y=7;
    a[8].x=2;a[8].y=2;

    /* Host -> GPU */
    magma_zmalloc(&dwork,ldwork);
    magma_zmalloc(&da,m*m);
    magma_zsetmatrix(m,m,a,m,da,m,queue);

    /*LU and Inverse */
    magma_zgetrf_gpu(m,m,da,m,piv,&info);
    magma_zgetri_gpu(m,da,m,piv,dwork,ldwork,&info);

    /*GPU -> Host */
    magma_zgetmatrix(m,m,da,m,a,m,queue);

    /* display inv(a) */
    for (int i=0;i<(m*m);i++){
        printf("%f +%fi\n",a[i].x,a[i].y);
    }

    /* free memory */
    free(a);
    free(piv);
    magma_free(dwork);
    magma_free(da);
    magma_queue_destroy(queue);
    magma_finalize();

    return 0;
}

我用以下代码编译：gcc -magazine -lcudart Ccode.c -o Ccode.o

Answer 1

我的系统管理员已经找到了为什么独立的C代码工作，而mexed C代码没有。我会在这里发布原因，因为在Matlab中使用Magma时，任何遇到相同问题的人都会有所帮助。

1. 我使用的Matlab版本是2014a。这个版本的supported compiler是4.7.x.我使用更高版本的gcc来编译代码。我使用不同版本的GCC与matlab一直没有问题，尽管它给出了警告，但对于上面的代码它确实很重要。
2. Compile with the MKL_ilp64 flag when using Magma with Matlab确保magma_int_t为int64。

有了这两个建议，可以将Magma混合到matlab中，没有任何问题。