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?
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]
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]
Indeterminate
output? $\endgroup$printf
, all the values are right.@Szabolcs $\endgroup$