Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.

mkl with matlab mex interface

numerix1
Beginner
2,813 Views

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/GUID-4326D19C-246A-4A7B-9728-03D371FB39AD.htm

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

 

0 Kudos
45 Replies
mecej4
Honored Contributor III
834 Views

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.

0 Kudos
numerix1
Beginner
834 Views

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.

0 Kudos
mecej4
Honored Contributor III
834 Views

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 = ((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,tauarr); 
       for(j=0; j<n; j++)printf(" %12.4f",pr[i+j*m]);
       }

    // call with random data
    for(i=0; i<n; i++)Iarr=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,tauarr); 
       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

 

0 Kudos
numerix1
Beginner
834 Views

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.

0 Kudos
mecej4
Honored Contributor III
834 Views

Let us try something simpler. Please build the tiny example code in https://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,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,A[i+3],A[i+6],jpiv,tau);
}

 

0 Kudos
numerix1
Beginner
834 Views

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 = ((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,tauarr); 
       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
$

0 Kudos
numerix1
Beginner
834 Views

For your program in 26,

icc -mkl -shared -fpic test_prog2.c

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

0 Kudos
mecej4
Honored Contributor III
834 Views

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.

0 Kudos
VipinKumar_E_Intel
834 Views

 

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.

 

0 Kudos
numerix1
Beginner
834 Views

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.

0 Kudos
mecej4
Honored Contributor III
834 Views

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,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,A[i+3],A[i+6],jpiv,tau);
}

 

0 Kudos
numerix1
Beginner
834 Views
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,tauarr); 
       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.

0 Kudos
Ying_H_Intel
Employee
834 Views

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

0 Kudos
numerix1
Beginner
834 Views

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 = ((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,tauarr); 
       for(j=0; j<n; j++)printf(" %12.4f",pr[i+j*m]);
       }
    printf("\n"); 
   
    free(data);
    free(Iarr);
    free(tauarr); 
    return;
}

 

0 Kudos
mecej4
Honored Contributor III
834 Views

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];

 

0 Kudos
numerix1
Beginner
834 Views

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 = ((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,tauarr); 
       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.

0 Kudos
mecej4
Honored Contributor III
834 Views

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.

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.

0 Kudos
numerix1
Beginner
834 Views

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 
>> 

 

0 Kudos
mecej4
Honored Contributor III
834 Views

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,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,A[i+3],A[i+6],jpiv,tau);
}

 

0 Kudos
numerix1
Beginner
814 Views

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 = ((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,tauarr); 
       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
$
0 Kudos
mecej4
Honored Contributor III
814 Views

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.

0 Kudos
Reply