下面的mexed c代码只是调用magma来反转矩阵。独立的C代码(也被发布)可以工作,但MEX代码崩溃。
我已经三次检查了文档,验证了其他magma函数是否按预期工作,以及
posted on the Magma forum and was told my code is fine
(这篇文章是来自Magma论坛的交叉帖子)。这意味着问题出在墨西哥。我想知道是什么导致了MEX代码的SEG错误,以及如何让它按预期运行。
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-lmagma-lcudart ccode.c-o ccode.o
