在mex文件matlab中使用magma_dysevd

我尝试在matlab中编写使用magma库,所以我基本上编写了一个mex函数,它使用magma函数合并了c代码,然后将这个mex函数编译成mexa64文件,因此我可以在matlab中使用。

mexfunction或source c代码如下:(称为eig_magma)

#include  #include  #include  #include  #include  #include  // includes, project #include "flops.h" #include "magma.h" #include "magma_lapack.h" #include "testings.h" #include  #include "mex.h" #define A(i,j) A[i + j*lda] extern "C" void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { #define L_OUT plhs[0] #define A_IN prhs[0] #define S_OUT plhs[1] magma_init(); real_Double_t gpu_perf, gpu_time; double *h_A, *h_R, *h_work, *L,*S; double *w; magma_int_t *iwork; magma_int_t N, n2, info, lwork, liwork, lda, aux_iwork[1]; double aux_work[1]; magma_vec_t jobz = MagmaVec; magma_uplo_t uplo = MagmaLower; magma_int_t i,j; N = mxGetM(A_IN); lda = N; n2 = lda*N; // query for workspace sizes magma_dsyevd( jobz, uplo, N, NULL, lda, NULL, aux_work, -1, aux_iwork, -1, &info ); lwork = (magma_int_t) aux_work[0]; liwork = aux_iwork[0]; TESTING_MALLOC_CPU( h_A, double, n2); h_A = (double *)mxGetData(A_IN); L_OUT = mxCreateDoubleMatrix(N,N,mxREAL); L = mxGetPr(L_OUT); S_OUT = mxCreateDoubleMatrix(N,1,mxREAL); S = mxGetPr(S_OUT); //print and check printf("print out h_A\n"); for(i=0;i<N;i++) { for(j=0;j<N;j++) printf("%f\t",h_A[i+j*lda]); printf("\n"); } /* Allocate host memory for the matrix */ TESTING_MALLOC_CPU( w, double, N ); TESTING_MALLOC_CPU( iwork, magma_int_t, liwork ); TESTING_MALLOC_PIN( h_R, double, N*lda ); TESTING_MALLOC_PIN( h_work, double, lwork ); printf("allocate memory on cpu with pin!\n"); printf("copy h_A to h_R, and then print h_R\n"); memcpy( h_R, h_A, sizeof(double)*n2); // lapackf77_dlacpy( MagmaUpperLowerStr, &N, &N, h_A, &lda,h_R, &lda ); // print and check for(i=0;i<N;i++) { for(j=0;j<N;j++) { printf("%f\t",h_R[i+j*lda]); } printf("\n"); } /* Performs operation using MAGA */ gpu_time = magma_wtime(); printf("start eig\n"); magma_dsyevd( jobz, uplo, N, h_R, lda, w, h_work, lwork, iwork, liwork, &info ); gpu_time = magma_wtime() - gpu_time; printf("%d\n",info); if (info != 0) printf("magma_dsyevd returned error %d: %s.\n", (int) info, magma_strerror( info )); printf("convey the output\n"); memcpy(L,h_R,sizeof(double)*n2); memcpy(S,w,sizeof(double)*N); TESTING_FREE_CPU( w ); TESTING_FREE_CPU( iwork ); TESTING_FREE_PIN( h_R ); TESTING_FREE_PIN( h_work); TESTING_FINALIZE(); } 

因为使用岩浆需要cuda和magma lib我编写makefile来将代码编译成mex文件:

 # Paths where MAGMA, CUDA, and OpenBLAS are installed MAGMADIR = /home/yuxin/magma-1.5.0 CUDADIR = /usr/local/cuda MATLABHOME = /opt/share/MATLAB/R2012b.app/ CC = icpc LD = icpc CFLAGS = -Wall LDFLAGS = -Wall -mkl=parallel MEX_CFLAGS = -fPIC -ansi -pthread -DMX_COMPAT_32 \ -DMATLAB_MEX_FILE MAGMA_CFLAGS := -DADD_ -DHAVE_CUBLAS -I$(MAGMADIR)/include -I$(CUDADIR)/include -I$(MAGMADIR)/testing MAGMA_LIBS := -L$(MAGMADIR)/lib -L$(CUDADIR)/lib64 \ -lmagma -lcublas -lcudart MEXLIBS := -L/opt/share/MATLAB/R2012b.app/bin/glnxa64 -lmx -lmex -lmat -lm -lstdc++ MEX_INCLU := -I$(MATLABHOME)/extern/include -I$(MATLABHOME)/simulink/include MEXFLAGS := -shared OBJECT = eig_magma.o EXECUTABLE = eig_magma $(EXECUTABLE): $(OBJECT) $(CC) $(LDFLAGS) $(MEXFLAGS) $(OBJECT) -o $@ $(MAGMA_LIBS) $(MEXLIBS) eig_magma.o: eig_magma.c $(CC) $(CFLAGS) $(MAGMA_CFLAGS) $(MEX_INCLU) $(MEX_CFLAGS) -c $< -o $@ clean: rm -rf eig_magma *.o eig_magma.mexa64 

它可以成功编译,但是当我在matlab中运行eig_magma时,以及当matlab执行时

magma_dsyevd(jobz,uplo,N,h_R,lda,w,h_work,lwork,iwork,liwork,&info);

matlab崩溃…….我也试过只在c中编写eig_magma,不使用mex函数,它编译成功并且工作正常。

任何人都知道如何解决这个问题?

谢谢

雨欣

让我举个例子。 由于我之前从未使用过MAGMA,因此我只使用常规LAPACK显示相同的代码(代码改编自我 之前给出的答案 )。 将其转换为使用MAGMA等效函数应该不难。

请注意,MATLAB已经附带了BLAS / LAPACK库以及必要的头文件供您在MEX文件中使用。 实际上这是英特尔MKL库,因此它是一个优化良好的实现。

eig_lapack.cpp

 #include "mex.h" #include "lapack.h" /* // the following prototype is defined in "lapack.h" header extern void dsyevd(char *jobz, char *uplo, ptrdiff_t *n, double *a, ptrdiff_t *lda, double *w, double *work, ptrdiff_t *lwork, ptrdiff_t *iwork, ptrdiff_t *liwork, ptrdiff_t *info); */ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { // validate input if (nrhs != 1 || nlhs > 2) mexErrMsgIdAndTxt("mex:error", "Wrong number of arguments."); if (!mxIsDouble(prhs[0]) || mxIsComplex(prhs[0])) mexErrMsgIdAndTxt("mex:error", "Input is not real double matrix."); mwSignedIndex m = mxGetM(prhs[0]), n = mxGetN(prhs[0]); if (m != n) mexErrMsgIdAndTxt("mex:error", "Input is not square symmetric matrix."); // allocate output plhs[0] = mxDuplicateArray(prhs[0]); plhs[1] = mxCreateDoubleMatrix(1, n, mxREAL); double *A = mxGetPr(plhs[0]); double *w = mxGetPr(plhs[1]); // query and allocate the optimal workspace size double workopt = 0; mwSignedIndex iworkopt = 0; mwSignedIndex lwork = -1, liwork = -1, info = 0; dsyevd("Vectors", "Upper", &n, A, &n, w, &workopt, &lwork, &iworkopt, &liwork, &info); lwork = (mwSignedIndex) workopt; liwork = (mwSignedIndex) iworkopt; double *work = (double*) mxMalloc(lwork * sizeof(double)); mwSignedIndex *iwork = (mwSignedIndex*) mxMalloc(liwork * sizeof(mwSignedIndex)); // compute eigenvectors/eigenvalues dsyevd("Vectors", "Upper", &n, A, &n, w, work, &lwork, iwork, &liwork, &info); mxFree(work); mxFree(iwork); // check successful computation if (info < 0) mexErrMsgIdAndTxt("mex:error", "Illegal values in arguments."); else if (info > 0) mexWarnMsgIdAndTxt("mex:warn", "Algorithm failed to converge."); } 

首先我们编译MEX文件(我在Windows上使用VS2013编译器):

 >> mex -largeArrayDims eig_lapack.cpp libmwlapack.lib 

现在我与内置的eig函数进行比较:

 >> A = [1 2 3 4 5; 2 2 3 4 5; 3 3 3 4 5; 4 4 4 4 5; 5 5 5 5 5] A = 1 2 3 4 5 2 2 3 4 5 3 3 3 4 5 4 4 4 4 5 5 5 5 5 5 >> [V,D] = eig(A) V = 0.6420 -0.5234 0.3778 -0.1982 0.3631 0.4404 0.1746 -0.6017 0.5176 0.3816 0.1006 0.6398 -0.0213 -0.6356 0.4196 -0.2708 0.2518 0.6143 0.5063 0.4791 -0.5572 -0.4720 -0.3427 -0.1801 0.5629 D = -3.1851 0 0 0 0 0 -0.7499 0 0 0 0 0 -0.3857 0 0 0 0 0 -0.2769 0 0 0 0 0 19.5976 >> [VV,DD] = eig_lapack(A) VV = 0.6420 -0.5234 0.3778 -0.1982 0.3631 0.4404 0.1746 -0.6017 0.5176 0.3816 0.1006 0.6398 -0.0213 -0.6356 0.4196 -0.2708 0.2518 0.6143 0.5063 0.4791 -0.5572 -0.4720 -0.3427 -0.1801 0.5629 DD = -3.1851 -0.7499 -0.3857 -0.2769 19.5976 

结果是相同的(虽然不能保证,看到特征向量是按比例定义的 ):

 >> V-VV, D-diag(DD) ans = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ans = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0