Why LibraryLink function crashes at a second time?

Question

Bug introduced in 9.0 or earlier and fixed in 10.3

---

Update

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

---

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?

Answer 1

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

Answer 2

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