Recently, I use the C to rewrite a Mathematica function NonzeroBasis[], please see here for fully descriptions.
Here, NonzeroBasis[i,p,u,U] calculates the values of $\{N_{i-p,p}(u),\cdots,N_{i-1,p}(u),N_{i,p}(u)\}$
The main algorithm as shown below:

LibraryLink Code
#include "WolframLibrary.h"
DLLEXPORT int NonzeroBasis(WolframLibraryData libData, mint Argc,
MArgument *Args, MArgument Res) {
/*define the argument-varible*/
mint i, p;
mreal u;
MTensor tensor_U;
mreal *U;
/*----------------------------*/
mint j, r;
mreal temp;
mreal *left, *right, *N;
MTensor tensor_N;
mint type = MType_Real;
mint dims[1];
mint rank = 1;
int err;
/*assign the argument-value to argument-varible*/
i = MArgument_getInteger(Args[0]);
p = MArgument_getInteger(Args[1]);
u = MArgument_getReal(Args[2]);
tensor_U = MArgument_getMTensor(Args[3]);
U = libData->MTensor_getRealData(tensor_U);
dims[0] = p + 1;
err = libData->MTensor_new(type, rank, dims, &tensor_N );
N = libData->MTensor_getRealData(tensor_N);
/*main implementation*/
N[0] = 1.0;
for(j = 1; j <= p; j++){
left[j] = u - U[i+1-j];
right[j] = U[i+j] - u;
saved = 0.0;
for(r = 0; r < j; r++){
temp = N[r]/(right[r+1] + left[j-r]);
N[r] = saved + right[r+1]*temp;
saved = left[j-r]*temp;
}
N[j] = saved;
}
/*return the result*/
MArgument_setReal(Res,tensor_N);
return LIBRARY_NO_ERROR;
}
However, when I compiled it via the following code in a Simplified Chinese operation system,
Needs["CCompilerDriver`"]
lib =
CreateLibrary[{"c:\\users\\Shutao TANG\\Desktop\\NonzeroBasis.c"},
"NonzeroBasis", "Debug" -> False]
(*load the .dll to Mathematica*)
NonzeroBasis =
LibraryFunctionLoad[lib, "NonzeroBasis",
{Integer, Integer, Real, {Real, 1}}, {Real, 1}]
it gives me this error(including garbled information)
Question
- How to debug the C code for LibraryLink? Because I cannot find useful information for debugging from
CreateLibrary::cmpper.
Update
Now I also implement the main algorithm in Mathematica. Note: the start index for the Mathematica is 1, while for C, the start index is 0.
mmaNonzeroBasis[i_, p_, u_, U_] :=
Module[{j, left, right, saved, r, temp, basis},
left = right = ConstantArray[0.0, {p + 1}];
basis = Prepend[ConstantArray[0.0, {p}], 1];
For[j = 1, j <= p, j++,
left[[j + 1]] = u - U[[i + 2 - j]];
right[[j + 1]] = U[[i + j + 1]] - u;
saved = 0.0;
For[r = 0, r < j, r++,
temp = basis[[r + 1]]/(right[[r + 2]] + left[[j - r + 1]]);
basis[[r + 1]] = saved + right[[r + 2]]*temp;
saved = left[[j - r + 1]]*temp
];
basis[[j + 1]] = saved;
];
basis
]
Thanks for Rolf Mertig's revision
add the variable declaration
mreal savedofsavedreplace the
MArgument_setReal()withMArgument_setMTensor()
knots = {0, 0, 0, 0, 0, 1, 2, 3, 4, 4, 5, 6, 7, 8, 9, 10, 10, 10, 10, 10};
mmaNonzeroBasis[9, 4, 4.5, knots]
NonzeroBasis[9, 4, 4.5, knots]
BSplineBasis[{4, knots}, #, 4.5] & /@ Range[5, 9]
(*{0.0104167,0.380208,0.545139,0.0616319,0.00260417}*)
Which is very strange. I don't know why? Could someone help me find the reason?




