sgetrf matlab mex file - crashing...

[Boxed Cylon]

Below is a test matlab mex routine for calling sgetrf. Insofar as I can tell (which isn't saying much...) the routine is correct, yet matlab just crashes with a Segmentation violation. I am at a loss for what might be wrong, other than the obvious "a memory allocation problem". Anyone have any suggestions for what I am missing? Thx. (I don't believe one can get away with calling cudaMallocHost within a mex routine.)

Code: Select all

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

#include "mex.h"
#include "cuda.h"
#include "cuda_runtime_api.h"
#include "cublas.h"
#include "magma.h"

#include "sys/time.h"

void mexFunction( int nlhs, mxArray *plhs[],
        int nrhs, const mxArray *prhs[])

{
      int L;
      int Lc;
      int dims0[2];

      // INPUT VARIABLES   %%%%%%%%%%%%%%%%%%%%%%%%%
      // A is dimensioned LXL
      float *A;

      // OUTPUT VARIABLE, X factorized   %%%%%%%%%%%%%%%%%%
      float *X;
 
      // CUDA/GPU VARIABLES %%%%%%%%%%%%%%%%%%%%%%%%
      float *ga;
      float *h_work;
      int* ipiv = 0;
      int info[1];

      if (nrhs != 1) {
         mexErrMsgTxt("gpu_sgetrf_magma requires 1 input arguments");
      } else if (nlhs != 1) {
         mexErrMsgTxt("gpu_sgesv_magma requires 1 output argument");
      }

      if ( !mxIsSingle(prhs[0]) ) {
           mexErrMsgTxt("Input array must be single precision.");
      }

 
// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// Single-precision input array */
// Dimensions, and then array data
      L = mxGetN(prhs[0]);
     // printf("L = %i\n",L);
      A =   (float*) mxGetData(prhs[0]);

// Left hand side matrix set up    (the solution) 
      dims0[0]=L;
      plhs[0] = mxCreateNumericArray(2,dims0,mxSINGLE_CLASS,mxREAL);
      X = (float*) mxGetData(plhs[0]);

      cuInit( 0 );
      cublasInit();

      int maxnb = magma_get_sgetrf_nb(L);
      int lwork = L*maxnb;

      printf("Get maxnb...\n");

      cublasAlloc((L+32)*(L+32) + 32*maxnb + lwork+2*maxnb*maxnb,sizeof(float), (void**)&ga);
      cudaMemset(ga,0,((L+32)*(L+32) + 32*maxnb + lwork+2*maxnb*maxnb)*4);

      Lc = (L/32)*32;
      if (Lc<L) Lc+=32;
      printf("%i, %i, %i, %i\n",L,Lc,maxnb,lwork);

      cublasSetMatrix (L, L, sizeof(float), A, L, ga, Lc);
    //  cublasSetMatrix (L, L, sizeof(float), A, L, (void *)ga, Lc);

      printf("Set A...\n");

      ipiv = (int *) mxCalloc (Lc,sizeof (int));

    //  cudaMallocHost( (void**)&h_work, (lwork+32*maxnb)*sizeof(float) );
      h_work = (float *) mxCalloc((lwork+32*maxnb),sizeof(float));

      printf("mxCalloc...\n");

    // Ready to go...
      magma_sgetrf_gpu(&L, &L, ga, &Lc, ipiv, h_work, info);
  //  magma_sgetrf_gpu(&N,  &N,  d_A, &lda, ipiv, h_work, info);
 // *****  CRASHES AT THIS CALL  TO MAGMA *****
      printf("sgetrf...\n");

    // Get the solution off the GPU
      cublasGetMatrix (L, L, sizeof(float), ga, Lc, X, L);

      printf("Solution...\n");
    // X has the solution we need; now back to matlab after a bit of clean up.

    // Print the first three elements of the first row (debugging)
      printf("X-top = %e %e %e\n",X[0],X[L],X[L+L]); 
    // Print the  last three elements of the  last row (debugging)
      printf("X-bottom = %e %e %e\n",X[L*(L-2)-1],X[L*(L-1)-1],X[L*L-1]); 

    // Clear the variables to avoid GPU memory leak (and GPU crash!)
      free(ipiv);
      free(h_work);
      cublasFree (ga);
      cublasShutdown(); 

}


[rnath]

Can you tell what did you use for the input matrix dimension? N X N
I mean , what is the value for N ?

[Boxed Cylon]

Can you tell what did you use for the input matrix dimension? N X N
I mean , what is the value for N ?

In this test case, I see that A was 131X131.

[rnath]

1. Do you get the same segmentation fault
using
cudaMallocHost( (void**)&h_work, (lwork+32*maxnb)*sizeof(float) );
instead of
h_work = (float *) mxCalloc((lwork+32*maxnb),sizeof(float));

2. What about using problem size 128 X 128?

[Stan Tomov]

I commented out the matlab related functionality in your code and ran it without problem for various Ls, including 131. Is it possible that the matrix A does not get a "good pointer", e.g. you have

Code: Select all

A =   (float*) mxGetData(prhs[0]);


This should point to L*L floats at least (if lda is L). I only replaced this (and the pointer for X) by allocating directly memory for A and initializing it

Code: Select all

 //  A =  (float*) mxGetData(prhs[0]);
 A =  (float*)malloc(L*L * sizeof(float));
 for(int j = 0; j < L*L; j++)
    A[j] = rand() / (float)RAND_MAX;


and the code runs.

[Boxed Cylon]

1. Do you get the same segmentation fault
using
cudaMallocHost( (void**)&h_work, (lwork+32*maxnb)*sizeof(float) );
instead of
h_work = (float *) mxCalloc((lwork+32*maxnb),sizeof(float));

2. What about using problem size 128 X 128?

I get the same fault using cudaMallocHost - that routine has never worked in a mex file for me, however.
I can try the smaller size. But I am off to MLS right now...Beckam is supposed to bend it.

[Boxed Cylon]

I commented out the matlab related functionality in your code and ran it without problem for various Ls, including 131. Is it possible that the matrix A does not get a "good pointer", e.g. you have

Code: Select all

A =   (float*) mxGetData(prhs[0]);


This should point to L*L floats at least (if lda is L). I only replaced this (and the pointer for X) by allocating directly memory for A and initializing it

Code: Select all

 //  A =  (float*) mxGetData(prhs[0]);
 A =  (float*)malloc(L*L * sizeof(float));
 for(int j = 0; j < L*L; j++)
    A[j] = rand() / (float)RAND_MAX;


and the code runs.

I can take a look at the mxGetData again. But this is a tried and true procedure I've used elsewhere. The issue is matlab's extreme memory sensitivity - things have to be done just right. mxGetData is filling A with the input from the call from matlab.

Thanks for taking a look at the issue for me. It is helpful that it runs without the matlab calls.

[Boxed Cylon]

I've wrestled with the problem some more. If I comment out the sgetrf call the routine runs fine. You tell me if you comment out the matlab stuff it runs fine (but not run from matlab as a mex routine I presume?). My present thinking is that the sgetrf calls BLAS routines/MKL that run into trouble with matlab causing the crash. A google search finds lots of material on this problem. One is supposed to link in the BLAS library statically to resolve the issue - I've tried to do that to no avail. This seems to me a matlab/mex issue, but I am at a loss for a solution. I am giving up for now; no magma sgetrf/matlab mex for me for now!

It has been my impression that Matlab mex files seem to allocate a specific area of memory, which is meshed in with the running matlab. Routines called from the mex file that allocate memory without telling matlab about it seem to cause the most problem - the memory is allocated over mexfile memory or something. Its a pain...

As an aside, it would be nice to have an equivalent routine to cudaMallocHost that could allocate pinned memory within a mex file. As I've said, that routine just causes matlab to crash.

(It was a subpar MLS cup game with a less-than-tight play on both sides. It was not raining, however; a small miracle there.)

[jpeinado]

Please please help!!!

Stan:

Could you answer to this question? : Did you get to run the Boxed Cyclon code as a mex file only by commenting out the sgetrf_gpu call....?

I tried to run a simmilar code using the shared magma libraries, but I cannot run it even in the case of commenting out the line sgetrf_gpu. If I link with the magma shared libraries (-lmagmablas -lmagma), my code crashes, even in the case of not use any magma call....

Thanks

jpeinado

I commented out the matlab related functionality in your code and ran it without problem for various Ls, including 131. Is it possible that the matrix A does not get a "good pointer", e.g. you have

Code: Select all

A =   (float*) mxGetData(prhs[0]);


This should point to L*L floats at least (if lda is L). I only replaced this (and the pointer for X) by allocating directly memory for A and initializing it

Code: Select all

 //  A =  (float*) mxGetData(prhs[0]);
 A =  (float*)malloc(L*L * sizeof(float));
 for(int j = 0; j < L*L; j++)
    A[j] = rand() / (float)RAND_MAX;


and the code runs.

[Stan Tomov]

Jpeinado,
I only commented out getting the matrix A using mxGetData. I directly allocated A and used the rest of the Boxed Cyclon code as a simple cpp code (not as a mex file), linking it with the shared libraries posted in the other thread ("Matlab/nvmex - failure to compile").
Stan