// mri_exp_mult_mex.c // Sangwoo Lee, University of Michigan, Jun, 2004 // 2004-6-21 modified by JF #include "mex.h" #include #include #define Usage \ "Usage: D = function(A, u, v);\n\ in:\n\ A [N L] complex matrix\n\ u [N 1] vector\n\ v [M 1] vector\n\ one (and only one) of u and v must be complex!\n\ out:\n\ D [L M] complex vector, D = A' * exp(-u * v.')\n\ D_lm = sum_n A_nl^* exp(-u_n v_m)\n\ \n\ This function is a memory efficient and fast implementation\n\ of AA matrix computation in int_tim_seg.m function in NUFFT package.\n\ \n" #ifndef Fail # define Fail(msg) { \ (void)fprintf(stderr, "FAIL %s %d: %s\n", __FILE__, __LINE__, msg); \ (void)fflush(stderr); \ return 0; } #endif #define dtype double static bool mri_exp_mult_double #include "mri_exp_mult_template.c" #undef dtype #define dtype float static bool mri_exp_mult_float #include "mri_exp_mult_template.c" #undef dtype // mri_exp_mult_mex() // secondary gateway static bool mri_exp_mult_mex( mxArray *plhs[], const mxArray *mx_A, const mxArray *mx_u, const mxArray *mx_v) { const int is_single = mxIsSingle(mx_A); double *work_r, *work_i; if (is_single ^ mxIsSingle(mx_u)) Fail("A & u type differ") if (is_single ^ mxIsSingle(mx_v)) Fail("A & v type differ") if (!mxIsComplex(mx_A)) Fail("A must be complex") if (!(mxIsComplex(mx_u) ^ mxIsComplex(mx_v))) Fail("need exactly one of u and v to be complex") // input sizes const int NN = mxGetM(mx_A); // N const int LL = mxGetN(mx_A); // L // check input size if (NN != (int) mxGetM(mx_u)) Fail("size of 2nd input != 1st dim of 1st input") if (1 != mxGetN(mx_u)) Fail("2nd input argument must be a column vector!") const int MM = mxGetM(mx_v); if (1 != mxGetN(mx_v)) Fail("3rd input argument must be a column vector!") // create output matrix [L M] if (!(plhs[0] = mxCreateNumericMatrix (LL, MM, is_single ? mxSINGLE_CLASS : mxDOUBLE_CLASS, mxCOMPLEX))) Fail("allocate output") if (!(work_r = (double *) calloc(NN, sizeof(*work_r)))) Fail("alloc") if (!(work_i = (double *) calloc(NN, sizeof(*work_i)))) Fail("alloc") if (is_single) { #define dtype float const dtype *Ar = (const dtype *) mxGetData(mx_A); const dtype *Ai = (const dtype *) mxGetImagData(mx_A); const dtype *ur = (const dtype *) mxGetData(mx_u); const dtype *ui = (const dtype *) mxGetImagData(mx_u); const dtype *vr = (const dtype *) mxGetData(mx_v); const dtype *vi = (const dtype *) mxGetImagData(mx_v); dtype *Dr = (dtype *) mxGetData(plhs[0]); dtype *Di = (dtype *) mxGetImagData(plhs[0]); #undef dtype if (!mri_exp_mult_float(Dr, Di, work_r, work_i, Ar, Ai, ur, ui, vr, vi, LL, MM, NN)) Fail("mri_exp_mult") } else { #define dtype double const dtype *Ar = (const dtype *) mxGetData(mx_A); const dtype *Ai = (const dtype *) mxGetImagData(mx_A); const dtype *ur = (const dtype *) mxGetData(mx_u); const dtype *ui = (const dtype *) mxGetImagData(mx_u); const dtype *vr = (const dtype *) mxGetData(mx_v); const dtype *vi = (const dtype *) mxGetImagData(mx_v); dtype *Dr = (dtype *) mxGetData(plhs[0]); dtype *Di = (dtype *) mxGetImagData(plhs[0]); #undef dtype if (!mri_exp_mult_double(Dr, Di, work_r, work_i, Ar, Ai, ur, ui, vr, vi, LL, MM, NN)) Fail("mri_exp_mult") } free(work_r); free(work_i); return 1; } // gateway routine - interface with Matlab void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { if (!nlhs && !nrhs) { #if 0 // defined(Date) printf("compiled %s\n", Date); #endif printf(Usage); return; } if (nlhs == 0 && nrhs == 1 && mxIsChar(prhs[0])) // check return; if (nlhs > 1 || nrhs != 3) fprintf(stderr, Usage); if (nlhs > 1) mexErrMsgTxt("One output argument required."); if (nrhs != 3) mexErrMsgTxt("Three input arguments required"); if (!mri_exp_mult_mex(plhs, prhs[0], prhs[1], prhs[2])) mexErrMsgTxt("mri_exp_mult_mex() failed"); } #undef Usage