mkl with matlab mex interface

mkl with matlab mex interface

Dear users:

I am having trouble interfacing my intel mkl code with Matlab using the mex interface. I have a code that does a certain factorization based on the pivoted QR decomposition. It works fine when I compile the code by itself. However, I would like to make a mex code for other users. I have done this before and have been able to compile simple mex files with icc and mkl. However, now I have a problem. Perhaps someone here has ideas for me to try. I am trying to use the pivoted QR decomposition function and this is where my code fails:

m = nrows; n = ncols;

data = (double*)mxCalloc(nrows*ncols,sizeof(double));

Iarr = (lapack_int*)mxCalloc(n,sizeof(lapack_int));
tauarr = (double*)mxCalloc(min(m,n),sizeof(double));

LAPACKE_dgeqp3(LAPACK_COL_MAJOR, (lapack_int)nrows, (lapack_int)ncols, data, (lapack_int)nrows, Iarr, tauarr);

sometimes it simply segfaults (when mex file is run inside matlab). Other times it complains that parameter 8 to dgeqp3 is incorrect.

Intel MKL ERROR: Parameter 8 was incorrect on entry to DGEQP3.

Parameter 8 is not present in the lapacke function, it is the info parameter in the fortran routines. Notice that I initialize space using mxCalloc vs regular C calloc..

https://software.intel.com/sites/products/documentation/hpc/mkl/mklman/G...

Here is how I compile my mex file:

#!/bin/bash

icc -c -DMX_COMPAT_32   -D_GNU_SOURCE -DMATLAB_MEX_FILE  -I"/apps/matlab2014a/extern/include" -I"/apps/matlab2014a/simulink/include" -fexceptions -fPIC -fno-omit-frame-pointer -openmp -shared -O -DNDEBUG mex_code1.c  -o mex_code1.o

icc -openmp -shared -L"/opt/intel/mkl/lib/intel64/" -I"/opt/intel/mkl/include/" /opt/intel/mkl/lib/intel64/libmkl_intel_lp64.a  /opt/intel/mkl/lib/intel64/libmkl_intel_thread.a /opt/intel/mkl/lib/intel64/libmkl_core.a /opt/intel/mkl/lib/intel64/libmkl_blas95_lp64.a  /opt/intel/mkl/lib/intel64/libmkl_lapack95_lp64.a -lmkl_rt -lmkl_core -lmkl_gnu_thread -lgomp -lpthread  -Wl,--no-undefined -Wl,-rpath-link,apps/matlab2014a/bin/glnxa64 -O -Wl,--version-script,"apps/matlab2014a/extern/lib/glnxa64/mexFunction.map" mex_code1.o   -L"apps/matlab2014a/bin/glnxa64" -lmx -lmex -lmat -lm -lstdc++ -o mex_code1.mexa64

rm -f mex_code1.o

I suspect the issue maybe in the compilation step (which compiles fine). As i mentioned, the code runs fine outside of the mex interface, when I replace back

all the mxCalloc calls by regular calloc.

thanks for your help

 

46 posts / 0 new
Last post
For more complete information about compiler optimizations, see our Optimization Notice.

Hi 

we have one article talking about Using Intel® MKL in MATLAB Executable (MEX) Files for your reference

in https://software.intel.com/en-us/articles/intel-math-kernel-library-inte...

i saw there are mixed mkl library in link line, You may remove some of them according to the MKL link adviser: 

https://software.intel.com/en-us/articles/intel-mkl-link-line-advisor/

The compile option and linked library is like below

 -openmp -I$(MKLROOT)/include

-Wl,--start-group $(MKLROOT)/lib/intel64/libmkl_intel_lp64.a $(MKLROOT)/lib/intel64/libmkl_core.a $(MKLROOT)/lib/intel64/libmkl_intel_thread.a -Wl,--end-group -L$(CompilerPath) -liomp5 -lpthread -lm

let's see any change if with correct library. 

Best Regards,

Ying 

Thanks for your suggestions. I've tried a few things and now the mex file does not segfault, but still unable to execute the QR function. Here is the code I use in the mex file:

    double *data = (double*)calloc(m*n,sizeof(double));
    for(i=0; i<(m*n); i++){
        data[i] = 1.0;
    }
    Iarr = (lapack_int*)calloc(n,sizeof(lapack_int));
    tauarr = (double*)calloc(min(m,n),sizeof(double));
    LAPACKE_dgeqp3(LAPACK_COL_MAJOR, m, n, data, n, Iarr, tauarr);
 

Here is the error reported in matlab:

Intel MKL ERROR: Parameter 4 was incorrect on entry to DGEQP3

Here is how I compile the mex file:

#!/bin/bash

icc -c -DMKL_ILP64 -fPIC -fno-omit-frame-pointer -DMX_COMPAT_32 -I"/apps/matlab2014a/extern/include" -I"/apps/matlab2014a/simulink/include" -I"/opt/intel/mkl/include" -fexceptions -openmp -shared -O -DNDEBUG mex_code.c  -o mex_code.o

icc -shared -openmp -I"/opt/intel/mkl/include/" -I"/opt/intel/include/" -L"/opt/intel/mkl/lib/intel64/"  -L"/opt/intel/lib/intel64/" -lmkl_rt -lmkl_core -lblas -lmkl_gnu_thread -lgomp  -liomp5 -lpthread -Wl,--no-undefined -Wl,-rpath-link,/apps/matlab2014a/bin/glnxa64 -O -Wl,--version-script,"/apps/matlab2014a/extern/lib/glnxa64/mexFunction.map" mex_code.o   -L"/apps/matlab2014a/bin/glnxa64" -lm -lmx -lmex -lmat -lstdc++ -o mex_code.mexa64

rm -f mex_code.o

I couldn't get the linking/compiling to work with the way Ying suggested (it kept complaining about undefined references to cblas functions in my code etc), but the above seems to not link duplicate libraries. Notice that before I was trying to do QR on data passed from matlab (allocated via mxCalloc), now I am simply allocating data inside the mex file and passing it to the QR function. For some reason it still doesn't like it (though now it complains about parameter 4 the data vs parameter 8 which doesn't exist in this function).

Thanks

 

 

 

Hi 

It seem the parameter n is not correct when you have LAPACK_COL_MAJOR.   lda INTEGER. The leading dimension of a; at least max(1, m)

How about change to m? 

Regarding your link line, it look too complex to understand :).  

If you need ILP 64 , while all Interge are 64bit ,   the mkl_rt need do run-time setting to align with ILP 64 interface.  So please make sure LP64 or ILP64. by default, mkl_rt use LP64. 

if you link mkl_rt (single Dynamic library) you don't need other mkl library. the -lmkl_core -lblas -lmkl_gnu_thread -lgomp  -liomp5 is don't need. 

If you'd like ILP64 dynamic link, 

 -L$(MKLROOT)/lib/intel64 -lmkl_intel_ilp64 -lmkl_core -lmkl_intel_thread -lpthread -lm -openmp  should be enough. 

You mentioned the static link : -Wl,--start-group $(MKLROOT)/lib/intel64/libmkl_intel_lp64.a $(MKLROOT)/lib/intel64/libmkl_core.a $(MKLROOT)/lib/intel64/libmkl_intel_thread.a -Wl,--end-group -lpthread -lm  -openmp

some cblas  symbols missing. could you copy the error detail or try write the library two/three  times in below order. 

$(MKLROOT)/lib/intel64/libmkl_intel_ilp64.a $(MKLROOT)/lib/intel64/libmkl_intel_thread.a  $(MKLROOT)/lib/intel64/libmkl_core.a $(MKLROOT)/lib/intel64/libmkl_intel_thread.a  $(MKLROOT)/lib/intel64/libmkl_core.a $(MKLROOT)/lib/intel64/libmkl_intel_thread.a  $(MKLROOT)/lib/intel64/libmkl_core.a  lpthread -lm  -openmp and let me know how it works? 

I may suggest to don't mix the dynamic and static library, intel thread and gnu thread library. Here is one article about libraries and linkage model of mkl for your reference 

https://software.intel.com/en-us/articles/intel-math-kernel-library-intel-mkl-linkage-and-distribution-quick-reference-guide

Best Regards,

Ying

 

 

 

Numerix1: once you get the MKL calls corrected, be sure to pass a meaningful matrix to the GEQP3 routine. A matrix with all elements set to 1 is rank-deficient and, if you do not remember this, many mysterious errors can be encountered later in your program.

Thanks for your suggestions. To start with, the m vs n is not the problem. I have just tested the following code:

 

    srand(time(NULL));
    double *data = (double*)calloc(m*n,sizeof(double));
    for(i=0; i<(m*n); i++){
        data[i] = ((double) rand() / (RAND_MAX));
        printf("data[%d] = %f\n", i, data[i]);
    }
    Iarr = (lapack_int*)calloc(n,sizeof(lapack_int));
    tauarr = (double*)calloc(min(m,n),sizeof(double));
    LAPACKE_dgeqp3(LAPACK_COL_MAJOR, m, n, data, m, Iarr, tauarr);
 

and it still returns:

Intel MKL ERROR: Parameter 4 was incorrect on entry to DGEQP3. It is true that it should be m since it is column major format and that the matrix of all ones I used is rank-deficient thanks for pointing this out, I use now a random matrix and m, but it still gives the same error. I will try now to compile differently as suggested by Ying.

 

Dear all,

Here is some more details on my problem. I've simplified my code to try to narrow down the problem (which I think is probably in the mex file compilation step). Here is my code:

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include "mkl.h"
#include "mkl_lapacke.h"
#include "matrix.h"
#include "mex.h"



void mexFunction( int nlhs, mxArray *plhs[],
             int nrhs, const mxArray *prhs[] )
{
    int i, j, m, n;
    time_t start_time, end_time;
    double  *pr;
    lapack_int *Iarr;
    double *tauarr;

    /* Check for proper number of arguments. */
    if (nrhs != 1) 
        mexErrMsgTxt("One input required, the matrix.\n");

    srand(time(NULL));

    /* get input parameters */
    n = mxGetN(prhs[0]);
    m = mxGetM(prhs[0]);
    pr = mxGetPr(prhs[0]);
    double *data = (double*)calloc(m*n,sizeof(double));
    for(i=0; i<(m*n); i++){
        // neither of below options work (matlab data or in mex initialized data)
        //data[i] = pr[i];
        data[i] = ((double) rand() / (RAND_MAX));
    }

    printf("the matrix is of size %d x %d\n", m, n);

    printf("call QR\n");
    Iarr = (lapack_int*)calloc(n,sizeof(lapack_int));
    tauarr = (double*)calloc(min(m,n),sizeof(double));
    LAPACKE_dgeqp3(LAPACK_COL_MAJOR, m, n, data, m, Iarr, tauarr);
    printf("end calling QR\n");

   
    free(data);
    free(Iarr);
    free(tauarr); 
    return;
}

Here the script that compiles this:

#!/bin/bash

icc -c -DMKL_ILP64 -fPIC -fno-omit-frame-pointer -DMX_COMPAT_32 -I"/apps/matlab2014a/extern/include" -I"/apps/matlab2014a/simulink/include" -I"/opt/intel/mkl/include" -fexceptions -openmp -shared -O -DNDEBUG mex_code.c  -o mex_code.o

icc -shared -openmp -I"/opt/intel/mkl/include/" -I"/opt/intel/include/" -L"/opt/intel/mkl/lib/intel64/"  -L"/opt/intel/lib/intel64/" -lmkl_rt -lpthread -Wl,--no-undefined -Wl,-rpath-link,/apps/matlab2014a/bin/glnxa64 -O -Wl,--version-script,"/apps/matlab2014a/extern/lib/glnxa64/mexFunction.map" id_code_mex1.o   -L"/apps/matlab2014a/bin/glnxa64" -lm -lmx -lmex -lmat -lstdc++ -o mex_code.mexa64

rm -f mex_code.o

This compiles fine. Here the error in matlab:

>> M = rand(2,3);
>> mex_code(M)
the matrix is of size 2 x 3
call QR

Intel MKL ERROR: Parameter 8 was incorrect on entry to DGEQP3.
end calling QR
>>

When I try Ying's compilation suggestion I use the script:

#!/bin/bash

icc -c -DMKL_ILP64 -fPIC -fno-omit-frame-pointer -DMX_COMPAT_32 -I"/apps/matlab2014a/extern/include" -I"/apps/matlab2014a/simulink/include" -I"/opt/intel/mkl/include" -fexceptions -openmp -shared -O -DNDEBUG mex_code.c  -o mex_code.o

icc -openmp -I"/opt/intel/mkl/include/" -I"/opt/intel/include/" -L"/opt/intel/mkl/lib/intel64/"  -L"/opt/intel/lib/intel64/" -Wl,--start-group /opt/intel/mkl/lib/intel64/libmkl_intel_lp64.a /opt/intel/mkl/lib/intel64/libmkl_core.a /opt/intel/mkl/lib/intel64/libmkl_intel_thread.a  /opt/intel/mkl/lib/intel64/libmkl_intel_lp64.a /opt/intel/mkl/lib/intel64/libmkl_blas95_lp64.a /opt/intel/mkl/lib/intel64/libmkl_lapack95_lp64.a -Wl,--end-group  -liomp5 -lpthread -lm -Wl,--no-undefined -Wl,-rpath-link,/apps/matlab2014a/bin/glnxa64 -O -Wl,--version-script,"/apps/matlab2014a/extern/lib/glnxa64/mexFunction.map" mex_code.o   -L"/apps/matlab2014a/bin/glnxa64" -lmx -lmex -lmat -lstdc++ -o mex_code.mexa64

rm -f mex_code.o

when I compile I get undefined reference to lapacke function:

$./compile_mex.sh
/usr/lib64/gcc/x86_64-suse-linux/4.8/../../../../lib64/crt1.o: In function `_start':
/home/abuild/rpmbuild/BUILD/glibc-2.18/csu/../sysdeps/x86_64/start.S:118: undefined reference to `main'
mex_code.o: In function `mexFunction':
mex_code.c:(.text+0x139): undefined reference to `LAPACKE_dgeqp3'

matlab version 8.3.0.532 (R2014a)

icc version 14.0.3 (gcc version 4.8.0 compatibility)

I would be happy to try any suggestions. Thanks.

 

Oh and btw the undefined reference to main is b/c I missed the shared flag. Doing:

icc -shared -openmp -I"/opt/intel/mkl/include/" ... (rest same as before) just gives me undefined reference to the lapacke routine:

mex_code.o: In function `mexFunction':
/mex_code.c:(.text+0x139): undefined reference to `LAPACKE_dgeqp3'

 

You probably guaranteed run-time failure by compiling with -DMKL_ILP64, which is for 8-byte integers, and then linking with LP64 libraries, which expect 4-byte integer arguments.

Your second attempt failed with "undefined reference to `main'" because you left out the -shared option when linking the MEX file.

Yes thank you the main undefined is due to shared option. But when I change the flag to -DMKL_LP64 (from ILP64) I still get the error

Intel MKL ERROR: Parameter 8 was incorrect on entry to DGEQP3.

 

Note also that the problem is specific to lapacke functions. For example, doing cblas_dgemm seems to work ok.

Not many people here may use a system with the same OS, C compiler, MKL and Matlab as you, but the following may help. I have a 32-bit Matlab (2007b) on Windows 8.1-x64. I built a Mex file (DLL) with the following modification of your program:

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include "mkl.h"
#include "mkl_lapacke.h"
#include "matrix.h"
#include "mex.h"



void mexFunction( int nlhs, mxArray *plhs[],
             int nrhs, const mxArray *prhs[] )
{
    int i, j, m, n;
    time_t start_time, end_time;
    double  *pr;
    lapack_int *Iarr,info;
    double *tau,*pl;

    /* Check for proper number of arguments. */
    if (nrhs != 1)
        mexErrMsgTxt("One input required, the matrix.\n");

    srand(time(NULL));

    /* get input parameters */
    n = mxGetN(prhs[0]);
    m = mxGetM(prhs[0]);
    pr = mxGetPr(prhs[0]);

    printf("the matrix is of size %d x %d\n", m, n);

    printf("call QR\n");
    Iarr = (lapack_int*)calloc(n,sizeof(lapack_int));
    tau = (double*)calloc(min(m,n),sizeof(double));
    info=LAPACKE_dgeqp3(LAPACK_COL_MAJOR, m, n, pr, m, Iarr, tau);
    printf("end calling QR, info = %d\n",info);
    for(i=0; i<m; i++)printf("%12.3e %12.3e %12.3e %2d %12.3e\n",
       pr[i],pr[i+3],pr[i+6],Iarr[i],tau[i]);

    free(Iarr);
    free(tau);
    return;
}

Invoking the Mex function from the Matlab command window gave:

>> A=[  4.100e+000  -1.900e+000  -1.300e+000
  1.200e+000   5.100e+000  -2.400e+000
 -2.200e+000   3.500e+000   7.700e+000]

A =

    4.1000   -1.9000   -1.3000
    1.2000    5.1000   -2.4000
   -2.2000    3.5000    7.7000

>> dgeqp3(A)
the matrix is of size 3 x 3
call QR
end calling QR, info = 0
  8.169e+000   2.103e+000  -3.079e+000  3   1.159e+000
  2.534e-001  -6.119e+000   4.732e-001  2   1.999e+000
 -8.131e-001   2.003e-002   3.659e+000  1   0.000e+000
>> 

Perhaps, you can get this example running on your system, and then modify the source code to run your problem. Of course, you would also have to write code to get the results from the Mex context back to Matlab instead of printing to the console.

Thank you a lot mecej4.  What is your compile option?

HI numeric, 

Could you please try to call lapacke_sgeqp3  or dgeqp3 direclty (include the work space query, i double parameter 8 is from lapacke_dgeqp3_work()) and see if any issues? (sorry, I haven't matlab at hand)

Best Regards,

Ying

 

Quote:

Ying H (Intel) wrote:
What is your compile option?

I add the Matlab directories to %INCLUDE% and %LIB% within an IFort/Icl command window, then build the Mex file with the command

     icl /Qmkl /LD dgeqp3.c libmx.lib libmex.lib /link /export:mexFunction

 

Hi mecej4 and ying thanks a lot for your suggestions. First I tried mecej's approach- it seems to be same code but calling QR function on the matlab data so equivalent to:

double *data = (double*)mxCalloc(m*n,sizeof(double));
    for(i=0; i<(m*n); i++){
        data[i] = pr[i];
    }

    Iarr = (lapack_int*)calloc(n,sizeof(lapack_int));
    tauarr = (double*)calloc(min(m,n),sizeof(double));
    LAPACKE_dgeqp3(LAPACK_COL_MAJOR, m, n, data, m, Iarr, tauarr);

// should behave same as this:
LAPACKE_dgeqp3(LAPACK_COL_MAJOR, m, n, pr, m, Iarr, tauarr);

but this segfaults for me.. I am using Linux (I quote matlab and icc versions above, both very recent). Ying, could you please

explain what you mean by "Could you please try to call lapacke_sgeqp3  or dgeqp3 direclty " what do you mean directly? I have

tried calling it on data that I've initialized inside the mex file, but it segfaults or reports parameter 8.

 

 

 

 

I mean call dgeqp3(  9 parameter) directly instead of lapacke c interface. 

dgeqp3(&m, &n, a, &lda, jpvt, tau, work, &lwork, &info)

Best Regards,

Ying

 

Hi Ying,

Thanks. I tried

  Iarr = (lapack_int*)calloc(n,sizeof(lapack_int));
    tauarr = (double*)calloc(min(m,n),sizeof(double));

    double *lwork = (double*)malloc(2*(3*n+1)*sizeof(double));
    LAPACKE_dgeqp3_work(LAPACK_COL_MAJOR, m, n,
                                data, m, Iarr,
                                tauarr, lwork, 2*(3*n+1) );

But this fails. Does dgeqp3 exist as a function in C? I thought it is only in fortran according to docs. I will look more carefully in includes/

 

Numerix1, with you posting different snippets of code and with us offering comments based on those snippets, we could go in circles for eons.

Please post complete source code and instructions to allow us to reproduce the problem. Avoid using random numbers for filling input matrices, since they are -- random and not reproducible.

The names dgeqp3, etc., are Fortran subroutines, but they can be called from C if the Fortran calling conventions are obeyed. Secondly, routines such as LAPACKE_dgeqp3 are simply wrapper routines that provide a natural C interface to the user and map to the matching Fortran routine. The LAPACKE routines do not actually do any linear algebra computations.

I tried your 2 X 3 example using the same Mex file as in #12, and obtained normal output:

>> M=rand(2,3)

M =

    0.8147    0.1270    0.6324
    0.9058    0.9134    0.0975

>> dgeqp3(M)
the matrix is of size 2 x 3
call QR
end calling QR, info = 0
 -1.218e+000   5.164e-001   2.156e-314  1   1.669e+000
  4.455e-001  -4.954e-001   1.004e-314  2   0.000e+000

 

Hi, let me repost the code again (similar to #7 above) to clarify, sorry for any confusion.

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include "mkl.h"
#include "mkl_lapacke.h"
#include "matrix.h"
#include "mex.h"


void mexFunction( int nlhs, mxArray *plhs[],
             int nrhs, const mxArray *prhs[] )
{
    int i, j, m, n, k, solve_type;;
    time_t start_time, end_time;
    double  *pr;
    lapack_int *Iarr;
    double *tauarr;

    /* Check for proper number of arguments. */
    if (nrhs != 1) 
        mexErrMsgTxt("One input required, the matrix.\n");

    srand(time(NULL));

    /* get input parameters */
    n = mxGetN(prhs[0]);
    m = mxGetM(prhs[0]);
    pr = mxGetPr(prhs[0]);
    
    // can use mxCalloc here too..
    double *data = (double*)calloc(m*n,sizeof(double));
    
    for(i=0; i<(m*n); i++){
        data[i] = ((double) rand() / (RAND_MAX));
    }

    printf("the matrix is of size %d x %d\n", m, n);

    printf("call QR\n");
    Iarr = (lapack_int*)calloc(n,sizeof(lapack_int));
    tauarr = (double*)calloc(min(m,n),sizeof(double));
    
    // this doesn't work (calling on supplied matlab data)
    // matlab segfaults and exits
    LAPACKE_dgeqp3(LAPACK_COL_MAJOR, m, n, pr, m, Iarr, tauarr);
    
    // and this doesn't work either (calling on data declared in program)
    // no matter if calloc or mxCalloc was used, mxCalloc on data segfaults
    // calloc produces the error
    // Intel MKL ERROR: Parameter 8 was incorrect on entry to DGEQP3
    // mxCalloc of data results in segfault 
    LAPACKE_dgeqp3(LAPACK_COL_MAJOR, m, n, data, m, Iarr, tauarr);

    printf("after call to QR\n");

   
    free(data);
    free(Iarr);
    free(tauarr); 
    return;
}

the compilation of the mex file is via the script:

#!/bin/bash

icc -c -DMKL_LP64  -fPIC -fno-omit-frame-pointer -I"/apps/matlab2014a/extern/include" -I"/apps/matlab2014a/simulink/include" -I"/opt/intel/mkl/include" -fexceptions -openmp -shared -O -DNDEBUG mex_code.c  -o mex_code.o

icc -shared -openmp -I"/opt/intel/mkl/include/" -I"/opt/intel/include/" -L"/opt/intel/mkl/lib/intel64/"  -L"/opt/intel/lib/intel64/" -lmkl_rt -lpthread -Wl,--no-undefined -Wl,-rpath-link,/apps/matlab2014a/bin/glnxa64 -O -Wl,--version-script,"/apps/matlab2014a/extern/lib/glnxa64/mexFunction.map" mex_code.o   -L"/apps/matlab2014a/bin/glnxa64" -lm -lmx -lmex -lmat -lstdc++ -o mex_code.mexa64

rm -f mex_code.o

 

The error I get is either complain about parameter 8 if lapacke_dgeqp3 is called on data calloc'd in the program or segfault if called on data that's mxCalloc'd (or passed from matlab - pr).

 

I tried also (code snipped for full code above)

double *lwork = (double*)malloc(2*(3*n+1)*sizeof(double));
LAPACKE_dgeqp3_work(LAPACK_COL_MAJOR, m, n, data, m, Iarr, tauarr, lwork, 2*(3*n+1) );

which always segfaults (even if data is calloc'd  or mxCalloc'd or data is replaced by pr).

fails also if I try to query the workspace size

    /* Query optimal working size */
    double work_query;
    lapack_int lwork = -1;
    printf("work query..\n");
    LAPACKE_dgeqp3_work( LAPACK_COL_MAJOR, m, n, data, m, Iarr, tauarr, &work_query, lwork );
    printf("work query = %f\n", work_query);

I suspect the issue is with compilation. tried newest matlab prerelease 2014 b, same result. will attempt to call dgeqp3 directly.

 

I found that your program in #19 works fine (see #12 for information on my system) provided that the array Iarr is reset to zero before the second call to LAPACKE_dgeqp3. I tried with several different input matrices, and the execution was normal with all of them.

Thanks mecej4. Yeah Iarr should be set to zero before each call. But for me the first call to dgeqp3 fails, I don't make it past that. It is an issue with system (linux vs windows), matlab version, and compilation.

Here is a version that runs on openSuse Linux 12.3 X64 with 64-bit Matlab 2007b and Intel C 13.0.1.

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include "mkl.h"
#include "mkl_lapacke.h"
#include "matrix.h"
#include "mex.h"

#define MIN(a,b) ((a) < (b) ? (a) : (b))

void mexFunction( int nlhs, mxArray *plhs[],
             int nrhs, const mxArray *prhs[] )
{
    int i, j, m, n, k, info, solve_type;
    time_t start_time, end_time;
    double  *pr;
    lapack_int *Iarr;
    double *tauarr,*data;

    /* Check for proper number of arguments. */
    if (nrhs != 1) 
        mexErrMsgTxt("One input required, the matrix.\n");

    srand(time(NULL));

    /* get input parameters */
    n = mxGetN(prhs[0]);
    m = mxGetM(prhs[0]);
    pr = mxGetPr(prhs[0]);
    
    data = (double*)calloc(m*n,sizeof(double));
    
    for(i=0; i<(m*n); i++) data[i] = ((double) rand() / (RAND_MAX));
    
    printf("the matrix is of size %d x %d\n", m, n);

    printf("call QR\n");
    Iarr = (lapack_int*)calloc(n,sizeof(lapack_int));
    tauarr = (double*)calloc(MIN(m,n),sizeof(double));
    
    // call with data from Matlab
    info=LAPACKE_dgeqp3(LAPACK_COL_MAJOR, m, n, pr, m, Iarr, tauarr);
    
    printf("After first call to QR, info = %d\n",info);
    for(i=0; i<m; i++){
       printf("\n%2d %12.5f ",Iarr[i],tauarr[i]); 
       for(j=0; j<n; j++)printf(" %12.4f",pr[i+j*m]);
       }

    // call with random data
    for(i=0; i<n; i++)Iarr[i]=0;
    info=LAPACKE_dgeqp3(LAPACK_COL_MAJOR, m, n, data, m, Iarr, tauarr);

    printf("\n\nAfter second call to QR, info = %d\n",info);
    for(i=0; i<m; i++){
       printf("\n%2d %12.5f ",Iarr[i],tauarr[i]); 
       for(j=0; j<n; j++)printf(" %12.4f",data[i+j*m]);
       }
    printf("\nAbout to release arrays\n");

   
    free(data);
    free(Iarr);
    free(tauarr); 
    return;
}

I built the Mex file with the command

icc -mkl -fpic -shared -I ~/MIRROR/ML2007b/extern/include/ dgeqp3.c -L ~/MIRROR/ML2007b/bin/glnxa64/ -lmex -lmx -o dgeqp3.mexa64

 

Hi mecej4, thanks a lot but sorry I've had to edit my post... I realize it actually still doesn't work as it should. If I use your version of the code or mines, I get as output something like this:

>> test_mex1(M)
the matrix is of size 2 x 3
call QR

Intel MKL ERROR: Parameter 8 was incorrect on entry to DGEQP3.
After first call to QR, info = -9

 0      0.00000        0.2785       0.9575       0.1576
 0      0.00000        0.5469       0.9649       0.9706
Intel MKL ERROR: Parameter 8 was incorrect on entry to DGEQP3.

After second call to QR, info = -9

 0      0.00000        0.9832       0.5922       0.3149
 0      0.00000        0.1101       0.8079       0.2358
About to release arrays
>>

The compilation commands I tried were:

icc -mkl -fpic -shared -I"/apps/matlab2014b_pre/extern/include/" test_mex1.c -L"/apps/matlab2014b_pre/bin/glnxa64" -lmex -lmx -o test_mex1.mexa64
icc -openmp -mkl -fpic -shared -I "/apps/matlab2014b_pre/extern/include" test_mex1.c -L"/apps/matlab2014b_pre/bin/glnxa64" -L"/opt/intel/mkl/lib/intel64" -lmkl_rt -lmex -lmx -o test_mex1.mexa64

with matlab 2014b : it doesn't segfault but it still complains about parameter 8..

with matlab 2014a: it segfaults. os is opensuse 13.1.

Let us try something simpler. Please build the tiny example code ihttps://software.intel.com/en-us/forums/topic/394013#comment-1740656 (sorry, that is in Fortran) below. Use the same ICC options and the same version of Intel C as you have been using. Does this example also give the same error message (parameter 8 to DGEQP3 wrong)? Note that Matlab is out of the picture for this test.

#include <stdio.h>
#include <mkl.h>

main(){
double A[]={4.1,1.2,-2.2, -1.9,5.1,3.5, -1.3, -2.4, 7.7};
lapack_int jpiv[]={0,0,0};
lapack_int m,n,lda,info;
double tau[3]; int i;

m=n=lda=3;
for(i=0; i<m; i++)printf("%12.3e %12.3e %12.3e\n",
   A[i],A[i+3],A[i+6]);
printf("\n");
info=LAPACKE_dgeqp3(LAPACK_COL_MAJOR,m,n,A, lda, jpiv, tau);
printf("DGEQP3 returned info = %d\n",info);
for(i=0; i<m; i++)printf("%12.3e %12.3e %12.3e %2d %12.3e\n",
   A[i],A[i+3],A[i+6],jpiv[i],tau[i]);
}

 

Hi mecej4,

You mean just run the mex file as a standalone C program? Here is what I tried:

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include "mkl.h"
#include "mkl_lapacke.h"

#define MIN(a,b) ((a) < (b) ? (a) : (b))

void main()
{
    int i, j, m, n, info;
    time_t start_time, end_time;
    double  *pr;
    lapack_int *Iarr;
    double *tauarr,*data;

    srand(time(NULL));

    /* set input parameters */
    m = 2;
    n = 3;
    
    data = (double*)calloc(m*n,sizeof(double));
    
    for(i=0; i<(m*n); i++){
        data[i] = ((double) rand() / (RAND_MAX));
    }    

    printf("the matrix is of size %d x %d\n", m, n);

    printf("call QR\n");
    Iarr = (lapack_int*)calloc(n,sizeof(lapack_int));
    tauarr = (double*)calloc(MIN(m,n),sizeof(double));
    
    info=LAPACKE_dgeqp3(LAPACK_COL_MAJOR, m, n, data, m, Iarr, tauarr);
    
    printf("After first call to QR, info = %d\n",info);
    for(i=0; i<m; i++){
       printf("\n%2d %12.5f ",Iarr[i],tauarr[i]); 
       for(j=0; j<n; j++)printf(" %12.4f",data[i+j*m]);
       }



    printf("\nAbout to release arrays\n");   
    free(data);
    free(Iarr);
    free(tauarr); 
    return;
}

first compile via:

icc -mkl test_prog.c

and all is ok!

$ ./a.out

the matrix is of size 2 x 3
call QR
After first call to QR, info = 0

 2      1.95839       -1.0207      -0.4832      -0.3024
 1      0.00000        0.1458       0.5430       0.3673
About to release arrays

now compile via:

icc -openmp -mkl -fpic -shared -I "/apps/matlab2014a/extern/include" test_prog.c -L"/apps/matlab2014a/bin/glnxa64" -L"/opt/intel/mkl/lib/intel64" -lmkl_rt -lmex -lmx -o test_prog.mexa64

gives segementation fault after running executable. Not sure if this is what you meant. This is not relevant I guess outside matlab,

BUT compiling like this:

icc -openmp -mkl -fpic -shared -L"/opt/intel/mkl/lib/intel64" -lmkl_rt test_prog.c

gives segmentation fault when running a.out. However, when the shared option is removed all is ok!

icc -openmp -mkl -fpic -L"/opt/intel/mkl/lib/intel64" -lmkl_rt test_prog.c

$./a.out

the matrix is of size 2 x 3
call QR
After first call to QR, info = 0

 2      1.95212       -1.0207      -0.6132      -0.5600
 1      0.00000        0.1566       0.5782       0.4417
About to release arrays
$

For your program in 26,

icc -mkl -shared -fpic test_prog2.c

gives segmentation fault. When shared is removed all is ok.

The "-shared" option is not suitable when producing an executable. On Linux, some (esp. GNU origin-) shared libraries, when invoked as commands, simply run a stub that prints information about the shared library and then return to the shell. If the shared library does not contain such a stub, it may well seg-fault if run as a command.

What we need at this point is a stand-alone (i.e., not related to Matlab) program that uses MKL and results in a "Intel MKL ERROR: Parameter 8 was incorrect on entry to DGEQP3" error. Such a "reproducer" is what would give the MKL team something to work with.

 

Could you attach a standalone test case as mecej4 mentioned?  You may submit the issue to premier.intel.com as well with the test case.

 

Sorry I am unable to produce a standalone program with the error. If I make a program, it is similar to the one I post in # 27. compiling it simply with icc -mkl option causes no problems. The problems arise when compiling it against the matlab libraries (and using the -shared option). I have tested with matlab version 2014a and 2014b (pre-release). Sometimes they give different results (either segfault and exit upon running mexfile or report error about parameter 8). But the mex file never runs error free. But it is only in the mex file that I run into problems, not in standalone codes.

How can I call dgeqp3 directly in C without using LAPACKE_dgeqp3 or LAPACKE_dgeqp3_work routines? Maybe if I am able to call the routine directly I can offer more insight into the problems. Thanks.

Quote:

numerix1 wrote:
How can I call dgeqp3 directly in C without using LAPACKE_dgeqp3 or LAPACKE_dgeqp3_work routines? Maybe if I am able to call the routine directly I can offer more insight into the problems.
You can call Fortran subroutines from C either by (i) understanding and delivering the subroutine arguments using the conventions of the Fortran compiler used to compile MKL, or (ii) following the Fortran-C interoperability conventions and writing a Fortran wrapper subroutine that calls MKL and which you call from C..

 

I doubt that this line of inquiry will be of much help, but here is the version of the program in #26, modified to call the Fortran subroutine DGEQP3 directly, using Intel C and MKL on Windows 32-bit. Note the name decoration used for DGEQP3 (all caps, no additional underscores). On Linux, the decorated name for use with GCC is "dgeqp3_".

#include <stdio.h>
#include <mkl.h>

main(){
double A[]={4.1,1.2,-2.2, -1.9,5.1,3.5, -1.3, -2.4, 7.7};
lapack_int jpiv[]={0,0,0};
lapack_int m,n,lda,info;
double tau[3]; int i;
double work[200]; lapack_int lwork=200;

m=n=lda=3;
for(i=0; i<m; i++)printf("%12.3e %12.3e %12.3e\n",
   A[i],A[i+3],A[i+6]);
printf("\n");
DGEQP3(&m,&n,A,&lda,jpiv,tau,work,&lwork,&info);
printf("DGEQP3 returned info = %d\n",info);
for(i=0; i<m; i++)printf("%12.3e %12.3e %12.3e %2d %12.3e\n",
   A[i],A[i+3],A[i+6],jpiv[i],tau[i]);
}

 

Thanks a lot mecej4, I tried running the code below. But same result:

Intel MKL ERROR: Parameter 8 was incorrect on entry to DGEQP3.
after QR
After first call to QR, info = -8
Here is the code:

#define min(x,y) (((x) < (y)) ? (x) : (y))
#define max(x,y) (((x) > (y)) ? (x) : (y))

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include "mkl.h"
#include "mkl_lapacke.h"
#include "matrix.h"
#include "mex.h"


void mexFunction( int nlhs, mxArray *plhs[],
             int nrhs, const mxArray *prhs[] )
{
    int i, j, m, n, k, solve_type, info;
    time_t start_time, end_time;
    double  *pr;
    lapack_int *Iarr;
    double *tauarr, *data;

    /* Check for proper number of arguments. */
    if (nrhs != 1) 
        mexErrMsgTxt("One input required, the matrix.\n");

    srand(time(NULL));

    /* get input parameters */
    n = mxGetN(prhs[0]);
    m = mxGetM(prhs[0]);
    pr = mxGetPr(prhs[0]);

    printf("the matrix is of size %d x %d\n", m, n);

    printf("call QR\n");
    Iarr = (lapack_int*)calloc(n,sizeof(lapack_int));
    tauarr = (double*)calloc(min(m,n),sizeof(double));

    lapack_int lwork = 2*min(m,n) - 1;
    double *work = (double*)malloc(lwork*sizeof(double));
    DGEQP3(&m,&n,pr,&m,Iarr,tauarr,work,&lwork,&info);
    printf("after QR\n");
    
    printf("After first call to QR, info = %d\n",info);
    for(i=0; i<m; i++){
       printf("\n%2d %12.5f ",Iarr[i],tauarr[i]); 
       for(j=0; j<n; j++)printf(" %12.4f",pr[i+j*m]);
       }
    printf("\n"); 
   
    free(Iarr);
    free(tauarr); 
    return;
}

compilation command:

icc -openmp -mkl -fpic -shared -I "/apps/matlab2014b_pre/extern/include" test_mex.c -L"/apps/matlab2014b_pre/bin/glnxa64" -L"/opt/intel/mkl/lib/intel64" -lmkl_rt -lmex -lmx -o test_mex.mexa64

Something is fishy with this icc/matlab combo that I am using (icc 14.03; matlab 2014a/b; opensuse 13.1). The fact that you can run this code fine with previous matlab versions suggests compatibility issues with these versions. p.s. I've tried various values for  lwork and also running dgeqp3 on data declared in program instead of pr and it gives same result.

Hi numerix1

1.  could you please try a bigger work size and see if any change ? 

From manual, it mentioned:  lwork: INTEGER. The size of the workarray; must be at least max(1, 3*n+1)for

real flavors, and at least max(1, n+1)for complex flavor

2. could you please try query the workspace as below code also, 

       lapack_int lwork = -1;
        DGEQP3(&m,&n,pr,&m,Iarr,tauarr,work,&lwork,&info);

printf("After first call to QR, info = %d\n",info);

        lwork = (MKL_INT) work;
        work = (double*) malloc( lwork * sizeof(double) );
        
        DGEQP3(&m,&n,pr,&m,Iarr,tauarr,work,&lwork,&info);

printf("After Second call to QR, info = %d\n",info);

Thanks

Ying

Hi Ying,

I tried changing lwork  to 3*n+1 , 1000*(3*n+1), etc. All give same result (parameter 8 error). When trying your suggestion #2, during compilation I get:

test_mex.c(57): warning #810: conversion from "double *" to "int" may lose significant bits
      lwork = (MKL_INT) work;

(as expected). When I run the mex file, it now simply segfaults upon reaching the first dgeqp3 call

printf("first call to QR..\n");
    DGEQP3(&m,&n,pr,&m,Iarr,tauarr,work,&lwork,&info);

and matlab quits.

Here is whole code which I ran for reference:

#define min(x,y) (((x) < (y)) ? (x) : (y))
#define max(x,y) (((x) > (y)) ? (x) : (y))

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include "mkl.h"
#include "mkl_lapacke.h"
#include "matrix.h"
#include "mex.h"



void mexFunction( int nlhs, mxArray *plhs[],
             int nrhs, const mxArray *prhs[] )
{
    int i, j, m, n, k, solve_type, info;
    time_t start_time, end_time;
    double  *pr;
    lapack_int *Iarr;
    double *tauarr, *data;

    /* Check for proper number of arguments. */
    if (nrhs != 1) 
        mexErrMsgTxt("One input required, the matrix.\n");

    srand(time(NULL));

    /* get input parameters */
    n = mxGetN(prhs[0]);
    m = mxGetM(prhs[0]);
    pr = mxGetPr(prhs[0]);
    
    data = (double*)calloc(m*n,sizeof(double));
    for(i=0; i<(m*n); i++){
        data[i] = ((double) rand() / (RAND_MAX));
    }

    printf("the matrix is of size %d x %d\n", m, n);

    printf("call QR\n");
    Iarr = (lapack_int*)calloc(n,sizeof(lapack_int));
    tauarr = (double*)calloc(min(m,n),sizeof(double));

    lapack_int lwork = -1;
    double *work = (double*)malloc(lwork*sizeof(double));

    printf("first call to QR..\n");
    DGEQP3(&m,&n,pr,&m,Iarr,tauarr,work,&lwork,&info);

printf("After first call to QR, info = %d\n",info);

    lwork = (MKL_INT) work;

    work = (double*) malloc( lwork * sizeof(double) );


    DGEQP3(&m,&n,pr,&m,Iarr,tauarr,work,&lwork,&info);

    printf("After Second call to QR, info = %d\n",info);


    printf("After first call to QR, info = %d\n",info);
    for(i=0; i<m; i++){
       printf("\n%2d %12.5f ",Iarr[i],tauarr[i]); 
       for(j=0; j<n; j++)printf(" %12.4f",pr[i+j*m]);
       }
    printf("\n"); 
   
    free(data);
    free(Iarr);
    free(tauarr); 
    return;
}

 

I think that the problem with the code segments in the last two posts is that, when a workspace optimal size query is made by setting lwork=-1, in order to have the results returned in work[0], work must be a valid pointer. Therefore, you should allocate the array to a minimum size of 1, or you should declare a double precision scalar variable and pass the address of that in the workspace query call. In other words,

double dlw,*work; int lwork;

m=n=lda=3;
lwork=-1; 
DGEQP3(&m,&n,A,&lda,jpiv,tau,&dlw,&lwork,&info);  // workspace query
lwork=dlw+0.5; work=malloc(lwork*sizeof(double)); // allocate optimum size
DGEQP3(&m,&n,A,&lda,jpiv,tau,work,&lwork,&info);  // workhorse query

 

The preceding posts suggest that the array work[] was unallocated at the time of the workspace query call, which would understandably cause an access violation. In particular, line-46 of the code in #35 calls malloc() with a negative size request, and the NULL that is returned by malloc() for such a meaningless request is then dereferenced to obtain the value of lwork. There is also an error on line-53; it should read

lwork = (MKL_INT) work[0];

 

Oops yes thanks mecej4 the previous code indeed had some problems. Sorry about that. Here is what I just ran:

#define min(x,y) (((x) < (y)) ? (x) : (y))
#define max(x,y) (((x) > (y)) ? (x) : (y))

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include "mkl.h"
#include "mkl_lapacke.h"
#include "matrix.h"
#include "mex.h"



void mexFunction( int nlhs, mxArray *plhs[],
             int nrhs, const mxArray *prhs[] )
{
    int i, j, m, n, k, solve_type, info;
    time_t start_time, end_time;
    double  *pr;
    lapack_int *Iarr;
    double *tauarr, *data;

    /* Check for proper number of arguments. */
    if (nrhs != 1) 
        mexErrMsgTxt("One input required, the matrix.\n");

    srand(time(NULL));

    /* get input parameters */
    n = mxGetN(prhs[0]);
    m = mxGetM(prhs[0]);
    pr = mxGetPr(prhs[0]);
    
    data = (double*)calloc(m*n,sizeof(double));
    for(i=0; i<(m*n); i++){
        data[i] = ((double) rand() / (RAND_MAX));
    }

    printf("the matrix is of size %d x %d\n", m, n);

    printf("call QR\n");
    Iarr = (lapack_int*)calloc(n,sizeof(lapack_int));
    tauarr = (double*)calloc(min(m,n),sizeof(double));

    lapack_int lwork = -1;
    double *work = (double*)malloc(sizeof(double));

    printf("first call to QR..\n");
    DGEQP3(&m,&n,pr,&m,Iarr,tauarr,work,&lwork,&info);

    printf("After first call to QR, info = %d and lwork = %d\n",info,lwork);

    lwork = (MKL_INT) work[0];
    printf("lwork from work array = %d\n", lwork);

    if(lwork > 0){
        work = (double*) malloc( lwork * sizeof(double) );
    }else{
        lwork = 1000;
        work = (double*) malloc( lwork * sizeof(double) );
    }

    DGEQP3(&m,&n,pr,&m,Iarr,tauarr,work,&lwork,&info);

    printf("After Second call to QR, info = %d\n",info);

    printf("printing result:\n");
    for(i=0; i<m; i++){
       printf("\n%2d %12.5f ",Iarr[i],tauarr[i]); 
       for(j=0; j<n; j++)printf(" %12.4f",pr[i+j*m]);
       }
    printf("\n"); 
   
    free(data);
    free(Iarr);
    free(tauarr); 
    return;
}

and the output:

>> A = randn(3,5);
>> test_mex4(A) 
the matrix is of size 3 x 5
call QR
first call to QR..

Intel MKL ERROR: Parameter 1 was incorrect on entry to DGEQP3.
After first call to QR, info = -1 and lwork = -1
lwork from work array = 0

Intel MKL ERROR: Parameter 1 was incorrect on entry to DGEQP3.
After Second call to QR, info = -1
printing result:

 0      0.00000 
 0      0.00000 
 0      0.00000 
>> 

It complains about parameter 1? That seems strange. Also, is there a default ordering that's assumed when calling dgeqp3 without the lapacke interface? there is no parameter specifying ordering, I assume it thinks input is column major by default? thanks.

I should expect the decorated name for calling MKL DGEQP3 from Intel C on 64 bit Linux to be dgeqp3_, not DGEQP3, as stated in #32.

Quote:

Also, is there a default ordering that's assumed when calling dgeqp3 without the lapacke interface? there is no parameter specifying ordering, I assume it thinks input is column major by default? 

The entry point DGEQP3 is the name of a Fortran-callable subroutine in MKL, and so uses Fortran calling conventions. Such entries have no awareness of being called from other languages, other Fortran compilers or even with the same compiler as the one used to build MKL but with different options and directives.

Sorry forgot, just tested it. I replaced both DGEQP3 calls by dgeqp3_ i.e.:

//DGEQP3(&m,&n,pr,&m,Iarr,tauarr,work,&lwork,&info); replaced by:
dgeqp3_(&m,&n,pr,&m,Iarr,tauarr,work,&lwork,&info);

it compiles the same and the output is exactly the same as for #37.

>> test_mex4(A)
the matrix is of size 3 x 4
call QR
first call to QR..

Intel MKL ERROR: Parameter 1 was incorrect on entry to DGEQP3.
After first call to QR, info = -1 and lwork = -1
lwork from work array = 0

Intel MKL ERROR: Parameter 1 was incorrect on entry to DGEQP3.
After Second call to QR, info = -1
printing result:

 0      0.00000 
 0      0.00000 
 0      0.00000 
>> 

 

Numerix1, I think that your clinging to the Matlab-coupled example imposes an unnecessary impediment to debugging your code. As far as the issues concerning the arguments passed to DGEQP3 are concerned, Matlab has no role that one can discern, since the dimensions of the matrix are being correctly passed from Matlab to the Mex process,

I really urge that you attempt to build and run the stand-alone example in #32, modified to make an optimum workspace size query as discussed above (in #36). The modified code is given below for your convenience. If you are able to get this example to run correctly on your system, it is quite likely that with the code put into a Mex function you will have better chances of building a functioning combination with Matlab..

#include <stdio.h>
#include <mkl.h>

main(){
double A[]={4.1,1.2,-2.2, -1.9,5.1,3.5, -1.3, -2.4, 7.7};
lapack_int jpiv[]={0,0,0};
lapack_int m,n,lda,info;
double tau[3]; int i;
double dlw,*work; int lwork;

m=n=lda=3;
for(i=0; i<m; i++)printf("%12.3e %12.3e %12.3e\n",
   A[i],A[i+3],A[i+6]);
printf("\n");
lwork=-1; dgeqp3_(&m,&n,A,&lda,jpiv,tau,&dlw,&lwork,&info);
lwork=dlw+0.5; printf("\n dlw = %.2f, lwork=%d\n",dlw,lwork);
work=malloc(lwork*sizeof(double));
dgeqp3_(&m,&n,A,&lda,jpiv,tau,work,&lwork,&info);
printf("DGEQP3 returned info = %d\n",info);
for(i=0; i<m; i++)printf("%12.3e %12.3e %12.3e %2d %12.3e\n",
   A[i],A[i+3],A[i+6],jpiv[i],tau[i]);
}

 

Hi mecej4, well the standalone C code always works, like in #27 and etc. I tried the code below with query of workspace and it works fine, just like all previous examples with lapacke_work function, etc. The issue seems to be in linking with matlab libraries and it seems to be version dependent since I understood you are able to run all these mex examples with your older setup.

Here is the standalone code I ran:

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include "mkl.h"
#include "mkl_lapacke.h"

#define MIN(a,b) ((a) < (b) ? (a) : (b))


void main()
{
    int i, j, m, n, info;
    time_t start_time, end_time;
    double  *pr;
    lapack_int *Iarr;
    double *tauarr,*data;

    srand(time(NULL));

    /* set input parameters */
    m = 20;
    n = 30;
    
    data = (double*)calloc(m*n,sizeof(double));
    
    for(i=0; i<(m*n); i++){
        data[i] = ((double) rand() / (RAND_MAX));
    }    

    printf("the matrix is of size %d x %d\n", m, n);

    printf("call QR\n");
    Iarr = (lapack_int*)calloc(n,sizeof(lapack_int));
    tauarr = (double*)calloc(MIN(m,n),sizeof(double));

    lapack_int lwork = -1;
    double *work = (double*)malloc(sizeof(double));

    printf("first call to QR to get workspace size..\n");
    dgeqp3_(&m,&n,data,&m,Iarr,tauarr,work,&lwork,&info);

    printf("After first call to QR, info = %d and lwork = %d\n",info,lwork);

    lwork = (MKL_INT) work[0];
    printf("lwork from work array = %d\n", lwork);

    if(lwork > 0){
        work = (double*) malloc( lwork * sizeof(double) );
    }else{
        lwork = 1000;
        work = (double*) malloc( lwork * sizeof(double) );
    }

    dgeqp3_(&m,&n,data,&m,Iarr,tauarr,work,&lwork,&info);

    printf("After Second call to QR, info = %d\n",info);

    dgeqp3_(&m,&n,data,&m,Iarr,tauarr,work,&lwork,&info);
    
    printf("After call to QR, info = %d\n",info);
    /*for(i=0; i<m; i++){
       printf("\n%2d %12.5f ",Iarr[i],tauarr[i]); 
       for(j=0; j<n; j++)printf(" %12.4f",data[i+j*m]);
       }*/



    printf("\nAbout to release arrays\n");   
    free(data);
    free(Iarr);
    free(tauarr); 
    return;
}

I compiled via: icc -mkl test_prog.c -o test_prog. When executing test_prog I get as expected:

$ ./test_prog
the matrix is of size 20 x 30
call QR
first call to QR to get workspace size..
After first call to QR, info = 0 and lwork = -1
lwork from work array = 1052
After Second call to QR, info = 0
After call to QR, info = 0

About to release arrays
$

Quote:

The issue seems to be in linking with matlab libraries and it seems to be version dependent since I understood you are able to run all these mex examples with your older setup.
Correct. I have run out of ideas on this problem.

Numerix1, have you reached any resolution of the problem discussed in this thread?

I recently gained access to a PC with Matlab 2014a - X64 installed (Windows 8.1 X64 OS), and ran your example in #37. 

The output from Matlab:

the matrix is of size 3 x 5
call QR
first call to QR..
After first call to QR, info = 0 and lwork = -1
lwork from work array = 16
After Second call to QR, info = 0
printing result:

 4      1.64038       -4.3247      -2.1266       1.8133       0.4651      -0.9858
 3      1.22685       -0.1903      -2.9305       1.3075       1.3496      -0.0427
 1      0.00000        0.4278       0.7938      -1.9383      -0.7192      -0.2596

 

I never got it to work properly with Linux based installations (I tried a few different Linux boxes I had access to with different Matlab versions up to 2014b and a couple of different icc versions). At the end I wrote some code to pass data back and forth through Linux based sockets within a mex file (i.e. to avoid passing data through disk). This way I can call inside mex files external library functions (such as from mkl) but do not need to link mkl when compiling the mex file. The mkl code then runs completely separately from the mex code and passes results back to mex via a simple socket server when finished (so no issues with openmp, etc). You can find my code @MatlabMexSockets http://tinyurl.com/nk34t83

This is what I currently use in my project to interface mkl code and Matlab. It does obviously have some overhead for data i/o but that can be fairly fast if executing on local machine and if socket transfer settings are set properly. I plan to write another version later using another IPC approach instead of sockets.

Thanks for your help with this.

Hi  Numerix1,

I checked what happens if create mex-file with MKL LAPACKE routine.

During Matlab execution LAPACKE routine is invoked from MKL library but base routine (dgeqp3_ for LAPACKE_dgeqp3) is invoked from Matlab library libmwlapack.so. libmwlapack.so supports 8-byte integers. So call of LAPACKE routine with LP64 interface leads to crash.

There are 3 options:

1) Use ILP64 libraries when compile mex file. Please be careful and ensure that:

·         integer variables passed to MKL routine is declared as lapack_int or MKL_INT

·         compile with –DMKL_ILP64

·         link with ILP64 libraries

Link line for this case:

icc -g -DMKL_ILP64 -fpic -shared -I${MKLROOT}/include -I "${matlab_root}/extern/include/" test_mex5.c -L${MKLROOT}/lib/intel64 -lmkl_intel_ilp64 -lmkl_intel_thread -lmkl_core -L"${matlab_root}/bin/glnxa64/" -lmex -lmx -o test_mex5.mexa64

 

2) Say matlab to use MKL libraries

export BLAS_VERSION=${MKLROOT}/lib/intel64/libmkl_rt.so

export LAPACK_VERSION=${MKLROOT}/lib/intel64/libmkl_rt.so

In this case all LAPACK routines will be invoked from MKL with specified interface (default is LP64,  export MKL_INTERFACE_LAYER=ILP64 to use ILP64 interface).

 

3) Also there is option to invoke LAPACK routine from libmwlapack.so. Great if there is no standalone MKL.

icc -fpic -shared -I "${matlab_root}/extern/include/" test_mex5.c -L"${matlab_root}/bin/glnxa64/" -lmwlapack -lmex -lmx -o test_mex5.mexa64

In this case you call dgeqp3_ from matlab library and need to pass 8-byte integers to routine.

Best Regards,

        Vladimir

 

I would add that I made these checks with 2013b Matlab version.

3.1) corrected version of option 3)

icc -DMKL_ILP64 -fpic -shared -I "${matlab_root}/extern/include/" test_mex5.c "${matlab_root}/bin/gln
xa64/"mkl.so -L"${matlab_root}/bin/glnxa64/" -lmex -lmx -o test_mex5.mexa64

Login to leave a comment.