4
$\begingroup$

I have a LibraryLink application and I have debugged it with VS. That is to say, the corresponding C code is totally right. However, when I test it in Mathematica, sometimes it gives the correct result, but at other times it also gives a result that contains the symbol Indeterminate with the same test case.

Has someone encountered this before?

enter image description here

Update

Here is a minimal example, the C code show as below:

src = "
#include <stdlib.h>
#include \"WolframLibrary.h\"

void matrix_band(const double *U, int n, double *K);

DLLEXPORT int MatrixBandValues(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res) {
    MTensor Ten_U = MArgument_getMTensor(Args[0]);
    MTensor Ten_K;
    mreal *U, *K;
    mint dims[1], n;
    int err;
    U = libData->MTensor_getRealData(Ten_U);
    n = libData->MTensor_getFlattenedLength(Ten_U) - 5;
    dims[0] = 4 * n - 2;
    err = libData->MTensor_new(MType_Real, 1, dims, &Ten_K);
    if (err) return err;
    K = libData->MTensor_getRealData(Ten_K);
    matrix_band(U, n, K);
    MArgument_setMTensor(Res, Ten_K);
    return LIBRARY_NO_ERROR;
}

//auxiliary function
void matrix_band(const double *U, int n, double *K) {
    double *d1, *d2, *d3;
    int i, idx;
    double *D, *E;
    double *K0, *K1, *K2;
    double *J0, *J1, *J2, *J3;
    //d1 ~ d3
    d1 = (double *) malloc((n + 1) * sizeof(double));
    d2 = (double *) malloc((n + 1) * sizeof(double));
    d3 = (double *) malloc((n + 1) * sizeof(double));
    //D E
    D = (double *) malloc((n + 3) * sizeof(double));
    E = (double *) malloc((n + 2) * sizeof(double));
    //K0 ~ K2
    K0 = (double *) malloc((n + 2) * sizeof(double));
    K1 = (double *) malloc((n + 1) * sizeof(double));
    K2 = (double *) malloc(n * sizeof(double));
    //J0 ~ J3
    J0 = (double *) malloc((n + 1) * sizeof(double));
    J1 = (double *) malloc((n + 1) * sizeof(double));
    J2 = (double *) malloc((n + 1) * sizeof(double));
    J3 = (double *) malloc((n + 1) * sizeof(double));
    //compute and store the difference of knots
    for (i = 3; i <= n; i++) d1[i] = U[i+1] - U[i];
    for (i = 2; i <= n; i++) d2[i] = U[i+2] - U[i];
    for (i = 1; i <= n; i++) d3[i] = U[i+3] - U[i];
    //compute D[i] <-- i = 2,...,n+2
    D[0] = D[1] = D[n+1] = D[n+2] = 0.0;
    D[2] = d1[3] / (3 * d2[2] * d2[2]);
    D[n] = d1[n] / (3 * d2[n] * d2[n]);
    for (i = 3; i <= n - 1; i++) D[i] = 1 / (3 * d2[i]);
    //compute E[i] <-- i = 2,...,n+1
    E[0] = E[1] = E[n] = E[n+1] = 0.0;
    for (i = 2; i <= n - 1; i++) E[i] = d1[i+1] / (6 * d2[i] * d2[i+1]);
    //K0[i], K1[i], K2[i]
    K0[0] = K0[n+1] = 0.0;
    for (i = 1; i <= n; i++) {
        K0[i] = 4 * (D[i] + D[i+1] - 2 * E[i]) / (d3[i] * d3[i]);
    }
    K1[0] = K1[n] = 0.0;
    for (i = 1; i <= n - 1; i++) {
        K1[i] = 4 * (E[i] + E[i+1] - D[i+1]) / (d3[i] * d3[i+1]);
    }
    K2[0] = K2[n-1] = 0.0;
    for (i = 1; i <= n - 1; i++) {
        K2[i] = -4 * E[i+1] / (d3[i] * d3[i+2]);
    }
    //J0 ~ J3
    for (i = 0; i <= n; i++)     J0[i] = 9 * (K0[i] + K0[i+1] - 2 * K1[i]);
    for (i = 0; i <= n - 1; i++) J1[i] = 9 * (K1[i] + K1[i+1] - K0[i+1] - K2[i]);
    for (i = 0; i <= n - 2; i++) J2[i] = 9 * (K2[i] + K2[i+1] - K1[i+1]);
    for (i = 0; i <= n - 3; i++) J3[i] = -9 * K2[i+1];
    //copy the J0~J3 to K
    for (i = 0; i <= n; i++) K[i] = J0[i];
    idx = n + 1;
    for (i = n + 1; i <= 2 * n; i++) K[i] = J1[i-idx];
    idx = 2 * n + 1;
    for (i = 2 * n + 1; i <= 3 * n - 1; i++) K[i] = J2[i-idx];
    idx = 3 * n;
    for (i = 3 * n; i <= 4 * n - 3; i++) K[i] = J3[i-idx];
    //free the pointers
    free(d1); free(d2); free(d3);
    free(D); free(E);
    free(K0); free(K1); free(K2);
    free(J0); free(J1); free(J2); free(J3);
}
"

Compiling

Needs["CCompilerDriver`"]

lib = CreateLibrary[src, "matrix_band", "ShellCommandFunction" -> Print, "ShellOutputFunction" -> Print]

Test

matrixBandVals = LibraryFunctionLoad[lib, "MatrixBandValues", {{Real, 1}}, {Real, 1}]


U = {0., 0., 0., 0., 0.0362025, 0.0469845, 0.0627937, 0.066165, 
     0.0779439, 0.0930326, 0.108092, 0.11111, 0.123131, 0.138142, 
     0.153134, 0.168137, 0.183173, 0.188376, 0.198237, 0.200901, 
     0.213309, 0.228371, 0.243438, 0.246434, 0.258541, 0.273696, 
     0.288951, 0.304406, 0.320201, 0.336438, 0.353182, 0.370496, 
     0.388384, 0.406853, 0.425852, 0.445343, 0.465241, 0.485511, 
     0.506113, 0.526968, 0.531064, 0.547984, 0.569105, 0.590322, 
     0.611571, 0.632737, 0.637278, 0.653768, 0.674676, 0.678707, 
     0.695508, 0.71631, 0.737158, 0.758177, 0.779429, 0.800986, 
     0.822908, 0.82933, 0.845221, 0.867945, 0.89131, 0.91687, 0.923411, 
     0.942926, 1., 1., 1., 1.};

 matrixBandVals[U]

enter image description here

Update

Here is a minimal example, the C code show as below:

src = "
#include <stdlib.h>
#include \"WolframLibrary.h\"

void matrix_band(const double *U, int n, double *K);

DLLEXPORT int MatrixBandValues(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res) {
    MTensor Ten_U = MArgument_getMTensor(Args[0]);
    MTensor Ten_K;
    mreal *U, *K;
    mint dims[1], n;
    int err;
    U = libData->MTensor_getRealData(Ten_U);
    n = libData->MTensor_getFlattenedLength(Ten_U) - 5;
    dims[0] = 4 * n - 2;
    err = libData->MTensor_new(MType_Real, 1, dims, &Ten_K);
    if (err) return err;
    K = libData->MTensor_getRealData(Ten_K);
    matrix_band(U, n, K);
    MArgument_setMTensor(Res, Ten_K);
    return LIBRARY_NO_ERROR;
}

//auxiliary function
void matrix_band(const double *U, int n, double *K) {
    double *d1, *d2, *d3;
    int i, idx;
    double *D, *E;
    double *K0, *K1, *K2;
    double *J0, *J1, *J2, *J3;
    //d1 ~ d3
    d1 = (double *) malloc((n + 1) * sizeof(double));
    d2 = (double *) malloc((n + 1) * sizeof(double));
    d3 = (double *) malloc((n + 1) * sizeof(double));
    //D E
    D = (double *) malloc((n + 3) * sizeof(double));
    E = (double *) malloc((n + 2) * sizeof(double));
    //K0 ~ K2
    K0 = (double *) malloc((n + 2) * sizeof(double));
    K1 = (double *) malloc((n + 1) * sizeof(double));
    K2 = (double *) malloc(n * sizeof(double));
    //J0 ~ J3
    J0 = (double *) malloc((n + 1) * sizeof(double));
    J1 = (double *) malloc((n + 1) * sizeof(double));
    J2 = (double *) malloc((n + 1) * sizeof(double));
    J3 = (double *) malloc((n + 1) * sizeof(double));
    //compute and store the difference of knots
    for (i = 3; i <= n; i++) d1[i] = U[i+1] - U[i];
    for (i = 2; i <= n; i++) d2[i] = U[i+2] - U[i];
    for (i = 1; i <= n; i++) d3[i] = U[i+3] - U[i];
    //compute D[i] <-- i = 2,...,n+2
    D[0] = D[1] = D[n+1] = D[n+2] = 0.0;
    D[2] = d1[3] / (3 * d2[2] * d2[2]);
    D[n] = d1[n] / (3 * d2[n] * d2[n]);
    for (i = 3; i <= n - 1; i++) D[i] = 1 / (3 * d2[i]);
    //compute E[i] <-- i = 2,...,n+1
    E[0] = E[1] = E[n] = E[n+1] = 0.0;
    for (i = 2; i <= n - 1; i++) E[i] = d1[i+1] / (6 * d2[i] * d2[i+1]);
    //K0[i], K1[i], K2[i]
    K0[0] = K0[n+1] = 0.0;
    for (i = 1; i <= n; i++) {
        K0[i] = 4 * (D[i] + D[i+1] - 2 * E[i]) / (d3[i] * d3[i]);
    }
    K1[0] = K1[n] = 0.0;
    for (i = 1; i <= n - 1; i++) {
        K1[i] = 4 * (E[i] + E[i+1] - D[i+1]) / (d3[i] * d3[i+1]);
    }
    K2[0] = K2[n-1] = 0.0;
    for (i = 1; i <= n - 1; i++) {
        K2[i] = -4 * E[i+1] / (d3[i] * d3[i+2]);
    }
    //J0 ~ J3
    for (i = 0; i <= n; i++)     J0[i] = 9 * (K0[i] + K0[i+1] - 2 * K1[i]);
    for (i = 0; i <= n - 1; i++) J1[i] = 9 * (K1[i] + K1[i+1] - K0[i+1] - K2[i]);
    for (i = 0; i <= n - 2; i++) J2[i] = 9 * (K2[i] + K2[i+1] - K1[i+1]);
    for (i = 0; i <= n - 3; i++) J3[i] = -9 * K2[i+1];
    //copy the J0~J3 to K
    for (i = 0; i <= n; i++) K[i] = J0[i];
    idx = n + 1;
    for (i = n + 1; i <= 2 * n; i++) K[i] = J1[i-idx];
    idx = 2 * n + 1;
    for (i = 2 * n + 1; i <= 3 * n - 1; i++) K[i] = J2[i-idx];
    idx = 3 * n;
    for (i = 3 * n; i <= 4 * n - 3; i++) K[i] = J3[i-idx];
    //free the pointers
    free(d1); free(d2); free(d3);
    free(D); free(E);
    free(K0); free(K1); free(K2);
    free(J0); free(J1); free(J2); free(J3);
}
"

Compiling

Needs["CCompilerDriver`"]

lib = CreateLibrary[src, "matrix_band", "ShellCommandFunction" -> Print, "ShellOutputFunction" -> Print]

Test

matrixBandVals = LibraryFunctionLoad[lib, "MatrixBandValues", {{Real, 1}}, {Real, 1}]


U = {0., 0., 0., 0., 0.0362025, 0.0469845, 0.0627937, 0.066165, 
     0.0779439, 0.0930326, 0.108092, 0.11111, 0.123131, 0.138142, 
     0.153134, 0.168137, 0.183173, 0.188376, 0.198237, 0.200901, 
     0.213309, 0.228371, 0.243438, 0.246434, 0.258541, 0.273696, 
     0.288951, 0.304406, 0.320201, 0.336438, 0.353182, 0.370496, 
     0.388384, 0.406853, 0.425852, 0.445343, 0.465241, 0.485511, 
     0.506113, 0.526968, 0.531064, 0.547984, 0.569105, 0.590322, 
     0.611571, 0.632737, 0.637278, 0.653768, 0.674676, 0.678707, 
     0.695508, 0.71631, 0.737158, 0.758177, 0.779429, 0.800986, 
     0.822908, 0.82933, 0.845221, 0.867945, 0.89131, 0.91687, 0.923411, 
     0.942926, 1., 1., 1., 1.};

 matrixBandVals[U]

enter image description here

Improve this question
$\endgroup$
5
  • $\begingroup$ I am afraid that we will need more details to help out here. Perhaps some code, to start. Also, can you isolate an input case that repeatably generates that Indeterminate output? $\endgroup$
    – MarcoB
    Feb 6, 2017 at 6:14
  • $\begingroup$ Regarding the update, have you written some code to verify that your function does not produce NaNs? $\endgroup$
    – Szabolcs
    Feb 8, 2017 at 10:01
  • $\begingroup$ In Visual Studio, I debug the C function with printf, all the values are right.@Szabolcs $\endgroup$
    – mma
    Feb 8, 2017 at 10:09
  • $\begingroup$ @mma Are you sure you didn't make a mistake? See my update to the answer. $\endgroup$
    – Szabolcs
    Feb 8, 2017 at 10:35
  • $\begingroup$ Very strange because I think the corresponding theory is right. Thanks for your help and I will check them in detail. @Szabolcs $\endgroup$
    – mma
    Feb 8, 2017 at 11:08

1 Answer 1

Reset to default
6
$\begingroup$

You will see an Indeterminate in a LibraryLink array (i.e. a packed array) if it contains a NaN value ("Not a Number").

You should be extremely cautious with these values, as they aren't properly supported by Mathematica.

By experience, I found that when I return a NaN or infinity through LibraryLink or MathLink, I get a packed array whose element prints as either Indeterminate or as NaN`. I do not understand what determines the difference between the two, though I noticed that using LibraryLink passing tends to result in Indeterminate and using MathLink passing tends to result in NaN`.

If I get NaN` or Inf` as a result, all bets are off about its behaviour. Just do not use it, because it triggers weird and inconsistent behaviour. There is no reliable way to work with these values in Mathematica.

If I get Indeterminate or Infinity, then I found that unpacking the array (Developer`FromPackedArray) produces a proper Indeterminate or Infinity symbol. However, if you leave these values in the packed array, and try to do operations on the packed array, it may again behave in inconsistent ways. Note that if you leave this values in the packed array, they aren't actual symbols (regardless or how they are printed on the screen).

Can we check if a packed array contains any Indeterminates? The only reliable way to do this in Mathematica is to first unpack the array. If the array is large, this can be slow and consume a lot of memory. Thus I recommend checking on the C side using isnan().

Warning: Do not trust any results from packed arrays containing such values. You may experiment to see if some operations work. When you find that they do, you may be tempted to use them. I did this in the past and I was bitten hard: the behaviour changed in a new version of Mathematica and my program started returning bad results.

Only use these arrays after unpacking, and only use them if they contain Indeterminate or Infinity. Never try to use values that print as NaN` or Inf` in Mathematica.

Update:

Your C code is indeed producing NaN values. The problem is not with transferring data from C to Mathematica. The problem is that the values you compute in C already contains NaNs.

I verified this by adding the following to MatrixBandValues:

for (i=0; i < 4*n-2; ++i)
    if (isnan(K[i]))
        printf("Element %d is NaN.\n", i);
fflush(stdout);

and running Mathematica in a terminal so I can see the printf output. Flushing is important, otherwise you may not see the output until Mathematica exits. It tells me:

Element 126 is NaN.
Element 188 is NaN.
Improve this answer
$\endgroup$
5
  • $\begingroup$ Thanks for your detailed explanation. My librarylink function returns the result that contained Inderteminate, then the corrponding result of Internal`PartitionRagged@res consists of NaN. $\endgroup$
    – mma
    Feb 6, 2017 at 10:14
  • $\begingroup$ According to my debug process in Visual Studio, I discovered that my C code is right. However, the librarylink function return the result that contained Indeterninate. It is very strange. So I would like to know is there a method to resolve it. $\endgroup$
    – mma
    Feb 6, 2017 at 10:17
  • 1
    $\begingroup$ @mma So you are saying that on the C side you have no NaN values in the array, but they suddenly appear when the array is transferred to Mathematica? What sort of evidence do you have for this? Did you loop through the array in C and test it with isnan()? Can you provide a minimal example? Also, I do not understand how PartitionRagged comes in here. Are you looking directly at the return value of the LibraryLink function? Or did you apply some functions to it before looking at it? $\endgroup$
    – Szabolcs
    Feb 7, 2017 at 7:59
  • 1
    $\begingroup$ If you find it hard to construct a minimal example, consider writing the array to a binary file and uploading that, so I can experiment with these bad values and I have them exactly as they are stored in memory. Do this in C (just fwrite the array), not in Mathematica. $\endgroup$
    – Szabolcs
    Feb 7, 2017 at 8:08
  • $\begingroup$ Please see my update @Szabolcs $\endgroup$
    – mma
    Feb 8, 2017 at 9:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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