Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I'm trying link one of Lapack's function to Mathematica through LibraryLink, and the Library function works great at the first execution but crashes at the second time.

Here is the code (modified from here)

MMASrc="
#include <WolframLibrary.h>
#include <WolframCompileLibrary.h>

DLLEXPORT mint WolframLibrary_getVersion(){
  return WolframLibraryVersion;
}
DLLEXPORT int WolframLibrary_initialize(WolframLibraryData libData){
  return 0;
}

EXTERN_C void zspsv_(char* UPLO, const mint* N, mint* NRHS, mcomplex* A, mint* IPIV, mcomplex* B, mint* LDB, mint* INFO);


EXTERN_C DLLEXPORT int lpkSolveAXBS(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res){

  MTensor ta=MArgument_getMTensor(Args[0]);
  MTensor tb=MArgument_getMTensor(Args[1]);
  mint N=*MTensor_getDimensionsMacro(tb);
  mcomplex* A=MTensor_getComplexDataMacro(ta);
  mcomplex* B=MTensor_getComplexDataMacro(tb);

  char U='U';
  mint IPIV[N];
  mint NRHS=1;
  mint LDA=N,LDB=N;
  mint INFO;
  mint i,j;

  mcomplex AP[N*(N+1)/2];

  for(j=0;j<N;j++)
    for(i=0;i<=j;i++)
      AP[i + (j+1)*j/2] = A[i*N+j];

  zspsv_(&U,&N,&NRHS,AP,IPIV,B,&LDB,&INFO);

  MArgument_setMTensor(Res,tb);
  return LIBRARY_NO_ERROR;
}


";
Needs["CCompilerDriver`"];

CreateLibrary[MMASrc,"lpk","Debug"->True,"TargetDirectory"->"/tmp",
        "CompileOptions"->"-llapack"];

A={{1. +0. I,0. +0. I,0. -6.94908*^-13 I,0. -6.94908*^-13 I},
           {0. +0. I,1. +0.0352595 I,0. -4.51893*^-11 I,0. -4.51893*^-11 I},
           {0. -6.94908*^-13 I,0. -4.51893*^-11 I,1. +0.0376938 I,0. +0. I},
           {0. -6.94908*^-13 I,0. -4.51893*^-11 I,0. +0. I,1. +0.0378932 I}};

B={1. +0. I,0. +0. I,0. +6.94908*^-13 I,0. +6.94908*^-13 I};

lpkSolve=LibraryFunctionLoad["/tmp/lpk.dylib","lpkSolveAXBS",
        {{Complex,2},{Complex,1}},{Complex,1}]
(*LibraryFunction[<>,lpkSolveAXBS,{{Complex,2},{Complex,1}},{Complex,1}]*)

Abs[LinearSolve[A,B]-lpkSolve[A,B]]//Max
(*2.02047*10^-28*)

Table[LinearSolve[A,B],{40000}];//AbsoluteTiming
(*{0.470767,Null}*)

Table[lpkSolve[A,B],{40000}];//AbsoluteTiming
(*{0.086038,Null}*)

Table[lpkSolve[A,B],{40000}];//AbsoluteTiming

In executing In[10], I hear a beep sound and then the kernel crashes and there is no output. Why the kernel crashes at the second time of invoking the LibraryLink function? How to fix that?

Update

As Todd Gayley commented under halirutan's answer, this is a bug.

share|improve this question
    
The link to the source of the code is missing. –  Sjoerd C. de Vries Jul 25 '13 at 6:42
    
@b.gatessucks thanks for point it out, I meant to use lpkSolve instead of LinearSolve in In[10], but I see both cause the kernel to crash. I'll correct it so that it wouldn't cause confusion. –  xslittlegrass Jul 25 '13 at 6:49
    
@SjoerdC.deVries I've corrected that. Thanks for pointing out. –  xslittlegrass Jul 25 '13 at 6:55
    
@xslittlegrass, I don't know the answer to your question, but would you be interested in a top level interface to BLAS and/or LAPACK? Concerning the crash, I'd send that to the technical support. –  user21 Jul 25 '13 at 9:30
1  
@xslittlegrass For the sake of future visitors, could you reconsider taking any of the current answers as accepted answer? I explained in chat that what you experience is clearly a bug and cloning the tensor is not a solution neither is the documentation wrong. I talked to Todd Gayley, one of the persons in charge of LibraryLink at Wolfram. Please look at his comment under my answer. –  halirutan Aug 6 '13 at 1:13
show 6 more comments

2 Answers 2

An observation, rather than an answer to the why question.

The question suggests, that the problem might have something to do with the used lapack library. As it turns out, at least on my machine it doesn't. We can replace the C-code with the most simple routine which is taking an input tensor and returning exactly the same tensor as result.

MTensor ta = MArgument_getMTensor(Args[0]);
MArgument_setMTensor(Res,ta);
return LIBRARY_NO _ERROR;

Now, the underlying problem seems to be the very small numbers which were used as input. We can pin it down to 0. + 6.94908*^-13 I which seems to be the cause. Unfortunately, there is another important ingredient to crash the kernel: it's the usage of Table. When I try the same with Do or ConstantArray it works. Array crashed my Mathematica completely. Therefore, to make it more fun for you too, please ensure that you really have some unsaved notebooks open, when you try this.

So while this crashes the kernel

Table[testFunc[{0. - 6.94908*^-13 I}], {40000}];
Table[testFunc[{0. - 6.94908*^-13 I}], {40000}];

this does not, even if I evaluate it several times

Do[testFunc[{0. - 6.94908*^-13 I}], {40000}];
Do[testFunc[{0. - 6.94908*^-13 I}], {40000}];

the same is true for ConstantArray.

Test Code

MMASrc =     
  "# include <WolframLibrary.h>
   # include <WolframCompileLibrary.h>

   DLLEXPORT mint WolframLibrary_getVersion(){
     return WolframLibraryVersion;
   }

   DLLEXPORT int WolframLibrary_initialize(WolframLibraryData libData){
     return 0;
   }

   EXTERN_C DLLEXPORT int test(
     WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res){

     MTensor ta=MArgument_getMTensor(Args[0]);
     MArgument_setMTensor(Res,ta);
     return LIBRARY_NO_ERROR;
   }";

Needs["CCompilerDriver`"];
lib = CreateLibrary[MMASrc, "test", "Debug" -> True, "TargetDirectory" -> "/tmp"];
testFunc = LibraryFunctionLoad[lib, "test", {{Complex, 1}}, {Complex, 1}];
share|improve this answer
    
I feel your pain, but did you really have unsaved notebooks open when testing this? –  Sjoerd C. de Vries Jul 25 '13 at 11:16
    
@SjoerdC.deVries Fortunately, only the notebook with the toy-code above. Therefore, I only had to edit the OP's code again to get the reduced example. –  halirutan Jul 25 '13 at 13:34
    
In principle, you don't need to use FileNames to get the full path of the libs as "SystemIncludeDirectories" -> Automatic should handle that. –  rcollyer Jul 25 '13 at 13:37
    
@rcollyer Yes, thanks. I simplified the post. –  halirutan Jul 25 '13 at 13:40
6  
The crash in this sample code is clearly a bug. It should work with the default Automatic passing. And it cannot be made to work no matter what tensor-passing convention you specify. I've reported the bug, and if I get any further details I will report back. –  Todd Gayley Jul 30 '13 at 6:09
add comment

Have a look at the memory management write up for tensor. I suspect you are returning a pointer to an automatic tensor that is being cleaned up.

LibrayLink Interaction

Can you elaborate more on your thoughts, because I read the documentation exactly the opposite way: "For automatic passing, the MTensor is owned by the library function only when that call to the library is active. The function can return the MTensor whether it has been modified or not, but whatever happens, the MTensor cannot be used after the function has returned

Yes. I believe the docs are in error. If you clone an input tensor, or create a new one, and return it the kernel will not crash.

It is my belief that on automatic copy in a temporary pointer is created. On an automatic return a copy is not made but returned directly. So if you return an input tensor you are returning a tensor pointing to memory that is going to be destroyed. A crash is inevitable.

This works just fine ...

    MMASrc = 
"#include <WolframLibrary.h>
#include <WolframCompileLibrary.h>

DLLEXPORT mint WolframLibrary_getVersion(){
     return WolframLibraryVersion;
}

DLLEXPORT int WolframLibrary_initialize(WolframLibraryData libData){
     return 0;
}
EXTERN_C DLLEXPORT int test(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res){

MTensor T0;
T0=0;
MTensor ta=MArgument_getMTensor(Args[0]);
libData->MTensor_clone( ta, &T0);

MArgument_setMTensor(Res,T0);
     return LIBRARY_NO_ERROR;
}";

Needs["CCompilerDriver`"];
lib=CreateLibrary[MMASrc,"test","Debug"->True, "TargetDirectory"->"/tmp"];
testFunc=LibraryFunctionLoad[lib,"test",{{Complex,1}},{Complex,1}];
Table[testFunc[{0.-6.94908^-13I}],{40000}];
Table[testFunc[{0.-6.94908^-13I}],{40000}];
Table[testFunc[{0.-6.94908^-13I}],{40000}];
Table[testFunc[{0.-6.94908^-13I}],{40000}];
share|improve this answer
    
Can you elaborate more on your thoughts, because I read the documentation exactly the opposite way: "For automatic passing, the MTensor is owned by the library function only when that call to the library is active. The function can return the MTensor whether it has been modified or not, but whatever happens, the MTensor cannot be used after the function has returned." –  halirutan Jul 26 '13 at 1:01
    
See my added comments above. –  Eric M Jul 27 '13 at 15:56
    
OK, +1 for trying this out. So this suggests, that there are two copies needed, because automatic passing already creates a copy of input tensors and then you have to clone it again? And why does this happens with Table and the small numbers only? –  halirutan Jul 27 '13 at 21:49
    
Thanks for the answer. It works! You mentioned "automatic copy in a temporary pointer is created", what do you mean the temporary pointer, what's its relation to the tensor ta created in the library function? If it is due to the clean process, why it still crash when we use manual pass as testFunc=LibraryFunctionLoad[lib,"test",{{Complex,1,"Manual"}},{Complex,1}];? Since in the manual pass, it will "pass in a copy of the MTensor and do not automatically clean". –  xslittlegrass Jul 29 '13 at 1:58
    
And in order to save the time and memory for large matrixes, is it possible to use the disown features or shared passing etc instead of using the cloning process? –  xslittlegrass Jul 29 '13 at 2:06
show 1 more comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.