An `eye` function that accepts the default integer and another integer kind simultaneously

Help
25

Interoperability is one viable justification. But doing so for all procedures will lead to an exponentially large number of interfaces that no language could handle, including Fortran. In the end, one has to choose, activate, and build only a handful of interfaces to use, depending on the problem.

1 Like
26

On a side note, which version of MATLAB specifically does pass int64 by default? I have Fortran packages that can be called from MATLAB that pass the default integer (int32) without breaking the code. The only version of MATLAB that we have not tested yet, is MATLAB 2021b.

1 Like
27

Sorry, my claim was not precise.

Here is what happens in my MATLAB (2020a).

  1. MWSIZE is defined as INTEGER*8 (sorry, this expression is not standard-conforming but it is by MATLAB);
  2. mexCallMATLAB requires its first input to be INTEGER*4;
  3. mexFunction requires its first input to be INTEGER (this input even shares the same name with the first input of mexCallMATLAB, very conveniently).

We see a very nice, careful, and patient mixture of almost all the possibilities. Welllllll done, MathWorks!

We may praise or blame MathWorks for doing this but cannot stop any user from interoperating with MATLAB.

If we take care of only the default integer, then disaster will happen, including particularly segmentation faults, which did occur a few times on my side and cost me a few days to debug.

2 Likes
29

One option with respect to MATLAB would be to use the C Matrix API in combination with compilers that support Fortran-C interoperability. According to the MATLAB documentation the kinds are:

MATLAB C API C type Fortran kind
mwSize size_t c_size_t
mwIndex size_t c_size_t
mwSignedIndex ptrdiff_t c_ptrdiff_t

The downside of using the C route is you cannot rely upon on the MATLAB mex command to perform the compilation and linking of Fortran routines correctly. Instead you have to provide the compiled object files directly to the mex C driver. Edit: another peril is that size_t is an unsigned type in C, but is read as an signed type in Fortran. This means in Fortran you can only use have of the range of size_t (on a 64-bit platform, this can take values betweeen 0 and 18,446,744,073,709,551,615)

I don’t know what it would take to nudge Mathworks to improve standard Fortran interoperability, instead of relying upon extensions like INTEGER*8. I guess there is not much user demand. An example of calling Fortran through the C MEX API can be found in the thread below:

3 Likes
30

@zaikunzhang Here is a performance benchmark of the reshape implementation of Eye vs. the subroutine implementation getEye() and function implementation genEye() in my previous comment. It looks like your implementation via reshape can be up to 16 times slower than getEye() and up to ~ 5 times slower than genEye(). The compiler (gfortran), array size, and system spec likely play a role in the efficiency idscrepancy. This benchmark excludes timing for repeated allocations. Including repeated allocations would likely increase the performance discrepancy even more.

Regardless, I find the reshape implementation much more elegant and concise than the implementations used in getEye and genEye. The performance impact should not be a concern at all if the code is not to be called on the order of billions of times or more.

1 Like
31

Thank you for pointing this out.

Since the major concern of mine was to deal with integer(IK) and integer, I did not pay much attention to how the identity matrix should be generated. I simply copied the code from the person who suggested this implementation (it is not by me).

I did make a modification when copying the code: originally, eye will output an identity matrix of integer type. I am not sure how useful such a matrix would be. It is hard for me to conceive a situation where we truly need an integer identity matrix instead of a real one. So I changed it to real. In my project, it would be rather real(RP) with RP being a parameter representing the real precision.

BTW, checking eye in stdlib, I find that it also outputs a matrix of integer type. What is the motivation? Does it mean that the users are expected to use it in a fashion like real(eye(n))? Maybe there is a real version that I overlooked.

32

I did a search for eye in the stdlib repository now, but could not readily find anything relevant. My impression is that the stdlib is automatically generated, so the code is generally comprehensive, covering all possibilities, perhaps including integer eye matrices. I modified your reshape method to directly output real(real64) in the benchmarks that I posted above. Otherwise, outputting integer and then converting it to real will make it twice slower than the timing that is in the benchmark, because it has to loop over the entire matrix twice: setting values and then conversion to real.

1 Like
33

There’s only one eye in stdlib_linalg and it returns a matrix of integer(int8) type. It has the smallest memory footprint of all type kinds, and it’s automatically cast to the correct type on assignment or arithmetic operation (for example with reals or higher integer kinds). There is some computational cost associated with the extra cast. If the application is performance-critical, eye should be called once and its result stored for reuse in loops.

2 Likes
34

Thank you @milancurcic for the explanation.

1 Like