5
$\begingroup$

I've the following code

#include "WolframLibrary.h"
#define bufferSize 60

EXTERN_C DLLEXPORT int getData(WolframLibraryData libData, mint Argc, MArgument*Args, MArgument Res)
{
    
    mint numOfData = MArgument_getInteger(Args[0]);
    
    mint  length[2]={numOfData,bufferSize}; 
    mint const *dimension=length;
    int rank =2, err;
    double **out;
    MTensor outTensor;

    err = libData->MTensor_new(MType_Real, rank, dimension, &outTensor);
    if (err) return err;
    
    out = libData->MTensor_getRealData(outTensor);

    for (int j=0;j<numOfData;j++)
       {
        for(int i=0;i<bufferSize;i++)
           {                 
             out[j][i] = i+1 + bufferSize*j; 
           } 
       }

    MArgument_setMTensor(Res, outTensor);

    return LIBRARY_NO_ERROR;
}

Problem is occurring in this line out = libData->MTensor_getRealData(outTensor);. This works file if it's a 1D double *out;, but I want to assign a 2D double **out;.

It's giving the following error

Needs["CCompilerDriver`"];
demoFile = FileNameJoin[{$CCompilerDirectory,"SystemFiles","CSource","mathNc.c"}];
lib = CreateLibrary[{demoFile}, "testDLL"];
datarcv = LibraryFunctionLoad[lib, "getData", {Integer}, {Real, 2}];

error: cannot convert 'mreal*' {aka 'double*'} to 'double**' in assignment

How may I assign a pointer to a rank 2 tensor/2D array?

$\endgroup$
1
  • $\begingroup$ "How may I assign a pointer to a rank 2 tensor/2D array?" Holding the shaft of the pointer, gently jab the tensor with the arrow end until the tensor "takes the point", so to speak. Works around 90% of the time. $\endgroup$ Feb 6, 2023 at 21:11

1 Answer 1

8
$\begingroup$

Don't even try that. This double indirection is not used in scientific programming because it makes it harder for the prefetcher to predict what you want. Thus it leads to higher latency and worse code execution speed. This is why the entries of MTensors are stored in a linear buffer. If out has bufferSize columns, then simply access the entry {i,j} with out[bufferSize * i + j].

You might argue that out[bufferSize * i + j] involves integer arithmetic. But modern processors typically ship with sufficiently many integer units to saturate the memory throughput. Moreover, the integer units are separate from the floating point units. So they do not compete with the latter; they can do their thing concurrently. Anyways, if you want to reduce the intensity of integer arithmetic, you can try things like this:

const mint n = numOfData * bufferSize;
for( mint row_begin = 0; row_begin < n; row_begin += bufferSize )
{
    for( mint i = 0; i < bufferSize; ++i )
    {
        out[row_begin + i] = row_begin + i + 1;
    }
}

This involves only additions. Not sure whether this improves anything. IMHO it obfuscates what one really wants to do. And telling from my experience, this is hardly any faster. I am not quite sure, but it may be that an integer addition costs exactly as much as an integer multiply-add. This depends on the hardware (and on whether the compiler is able to fuse the two operations into one). I have seen this idiom rarely for CPU code, but quite frequently for GPUs (whose cores are not that beefy as a CPU's).

Moreover, compiler developers are aware of these programming habits and thus they put some effort into optimizations like auto-vectorization. For most efficient vectorization, the data has to lie consecutively in memory (otherwise the copy-costs slow it down too much). Then short vectors of this data can be taken and processed as one. Problems arise at the beginning of a matrix row (due to alignment issues) and at the end of the row. The ends have to be peeled, i.e. processed in the "classical" way; this costs extra logic. So from this point of view, the loop might better be reformulated as follows:

const mint n = numOfData * bufferSize;

for( mint i = 0; i < n; ++i )
{
    out[i] = i + 1;
}

But maybe the compiler is clever enough to do it on its own (if bufferSize is a compile-time constant).

Anyways, the real bottleneck is always memory bandwidth and latency.

Edit:

If you have to deal with high-rank tensors, then you might want to use

libData->MTensor_getReal
libData->MTensor_setReal

to extract/modify a single entry. See here for details. I've got the hunch that this will be very slow because indices have to be stored in an array first and this array is read and processed for each query. Plus these routines might also involve bound checks which you might want to skip if you know what you are doing.

$\endgroup$
5
  • $\begingroup$ I've already written it for 1D similar to as you mentioned. I'm trying the 2D to keep rank and indices on par because it would be very difficult to keep track of the indices when the rank becomes 5 or 6. $\endgroup$
    – csk 7
    Feb 6, 2023 at 7:41
  • $\begingroup$ I see. Have a look at my edit. $\endgroup$ Feb 6, 2023 at 7:49
  • $\begingroup$ If I don't find any other simpler way I may need to use libData->MTensor_setReal directly on outTensor. I asked the chatGPT, it provided me with this function libData->MTensor_getRealDataArray() which doesn't exist! I don't know if it exists in 13.2. I'm using 12.3 $\endgroup$
    – csk 7
    Feb 6, 2023 at 8:34
  • 7
    $\begingroup$ " I asked the chatGPT" Don't do that. Really. $\endgroup$ Feb 6, 2023 at 8:48
  • $\begingroup$ You might be able to use experimental features of the Eigen C++ package to wrap a thin layer class around the pointer; the class provides the indexing. See eigen.tuxfamily.org/dox/unsupported/eigen_tensors.html $\endgroup$ Feb 6, 2023 at 8:50

Your Answer

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

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