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

Thank you @FortranFan for the detailed solution and explanation!

Does this work if on the user’s platform real_kind has only two elements rather than three? Or does the standard ensure that size(real_kind) == 3?

A minor side note is that real_kinds and integer_kinds are F2008 but I hope to have an implementation compatible with F2003.

Thank you @NormanKirkby. Yes, this would completely resolve the problem provided that I am not so lazy and careless to forget to include the kind argument almost every time.

I totally agree. The problem is, which kinds should be included in “different kinds of integers”? Is there a way to figure out which kinds are supported on the user’s platform before the compilation? Thank you.

No, the code will fail to compile if one of int8, int16, int32, or
int64 is not a valid kind type. You would need to look at the
available kinds by inspecting integer_kinds from iso_fortran_env.
In the places that I have used iso_fortran_env, you would have
something like use mytypes, only : ik1 => XXX, ik2 => YYY where
XXX and YYY have been data mined from iso_fortran_env.

2 Likes

Thank you @kargl for the explanation. Is there a way to do this data mining automatically so that I include it in my code and the users do not even need to know about it?

It would not be exactly ideal if the users have to modify the code manually, since they may not be capable of doing so.

This has been discussed previously here:

2 Likes

Thank you @certik for pointing this out. According to this nice blog by @everythingfunctional , the answer is “no, there is no automatic way”. The “most automatic” way is to write a program that outputs integer_kinds, and then write some code that acts according to the output. Do I understand it correctly or do I overlook something?

1 Like

Well, you want to support back to F2003, which means you’ll need
write a code to look for the available kinds. Try this

program foo
  integer j, k, n
  integer, parameter   :: i(40) = [(selected_int_kind(j), j= 1, 40)]
  write(*,'(A)')  'module mytype'
  write(*,'(A)')  '   implicit none'
  n = 1
  k = i(1)
  write(*, '(A,I0,A,I0)') &
  &   '   integer, parameter :: ik', n, ' = ', k
  n = n + 1
  do j = 1, 40
     if (i(j) == k) cycle
     if (i(j) == -1) exit
     k = i(j)
     write(*, '(A,I0,A,I0)') &
     &   '   integer, parameter :: ik', n, ' = ', k
     n = n + 1
  end do
  write(*,'(A)') 'end module mytype'
end program foo

gfortran on FreeBSD produces

module mytype
  implicit none
  integer, parameter :: ik1 = 1
  integer, parameter :: ik2 = 2
  integer, parameter :: ik3 = 4
  integer, parameter :: ik4 = 8
  integer, parameter :: ik5 = 16
end module mytype
5 Likes

Thank everyone for the valuable input!

The best solution I have seen up to now is the following (it is not by me).

module eye_mod
use, intrinsic :: iso_fortran_env, only : INTEGER_KINDS
implicit none

integer, parameter :: IK = INTEGER_KINDS(1)  ! Only an example. It can be whatever kind that is valid
integer, parameter, private :: ID = kind(1)
integer, parameter, private :: IOTHER = merge(ID, minval(INTEGER_KINDS, INTEGER_KINDS /= IK), ID /= IK)

interface eye
    module procedure eye_ik, eye_iother
end interface eye

contains

function eye_ik(n)
implicit none
integer(IK), intent(in) :: n
real :: eye_ik(n, n)
integer(IK) :: j, k
eye_ik = reshape([1.0, ([(0.0, k=1, n)], 1.0, j=1, n - 1)], shape(eye_ik))
end function

function eye_iother(n)
implicit none
integer(IOTHER), intent(in) :: n
real :: eye_iother(n, n)
integer(IOTHER) :: j, k
eye_iother = reshape([1.0, ([(0.0, k=1, n)], 1.0, j=1, n - 1)], shape(eye_iother))
end function

end module eye_mod

program testeye
use eye_mod
implicit none
print *, eye(3)
print *, eye(3_IK)
end program testeye

It works provided that the following assumptions hold.

  1. The user’s platform supports use, intrinsic :: iso_fortran_env, only : INTEGER_KINDS. This is F2008 rather than F2003 as I wanted, but it is OK since there is no better solution yet.

  2. On the user’s platform, size(INTEGER_KINDS) > 1. This is not guaranteed by the standard, but not too stringent.

  3. The user’s platform allows us to use ‘minval’ for initialization. Absoft 21.0 does not allow this. Does the standard say anything about it?

I have tested the code by the following compilers: gfortran (pass), ifort (pass), ifx (pass), nagfor (pass), nvfortran (pass), sunf95 (pass), flang (pass), AOCC flang (fail), absoft Fortran 21.0 (fail), g95 (fail).

1 Like

Here is a standalone generic subroutine and function implementation of the Eye matrix for all kinds of complex and real types (but only with default input integer kind in the function interface. Any non-default integer has really no practical storage or performance advantage from my perspective). The interface and usage doc is the same as getEye() subroutine and genEye function I posted above.

The test can be compiled and run with

gfortran -cpp -DRK64_ENABLED Constants_mod.f90 MatrixDiag_mod.f90 MatrixDiag_mod@Routines_smod.f90 main.f90 -o main.exe
./main.exe

To enable or disable interfaces, one has to set or remove the corresponding preprocessor flags in the compile command,

-DRK128_ENABLED 
-DRK64_ENABLED 
-DRK32_ENABLED 
-DCK128_ENABLED 
-DCK64_ENABLED 
-DCK32_ENABLED 

like -DRK64_ENABLED that is set for the example.

1 Like

I agree, but I do not think we can decide which integer kind will be used by the user.

For example, the latest MATLAB uses 64-bit integer as the default in MEX files even if it is not the default integer. I do not know whether this is reasonable, but it is a fact.

1 Like

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

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

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

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

@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

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.

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

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

Thank you @milancurcic for the explanation.

1 Like