Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Bug? MATLAB MEX changes the kind of the default logical

When interfacing a piece of Fortran 2003 (or above) code with MATLAB by MEX, I am surprised to find that MEX changes the kind of the default logical. This is fatal, because a piece of perfectly compilable Fortran code may fail to be mexified due to a mismatch of types, which did happen in my project.

Here is a minimal working example.

Name the following code as "test_kind.F", compile it by mex test_kind.F in MATLAB, and then run test_kind in MATLAB. This will produce a plain text file named fort.99, which contains two numbers "4" and then "8" as the result of the WRITE instructions.

! test_kind.F
! Tested by MATLAB 9.8.0.1323502 (R2020a) with GNU Fortran (Ubuntu 9.3.0-17ubuntu1~20.04) 9.3.0

#include "fintrf.h"

      subroutine mexFunction(nlhs, plhs, nrhs, prhs)

      use ieee_arithmetic, only : ieee_is_nan
      implicit none
      mwPointer plhs(*), prhs(*)
      integer nlhs, nrhs

      write(99, *) kind(ieee_is_nan(1.0))  ! This prints a number in fort.99
      write(99, *) kind(.false.)  ! A benchmark, which should print the same number in fort.99
      close(99)

      end subroutine mexFunction

I thought the two printed numbers should always be equal to each other, although the specific value depends on the compiler and needs not to be 4 or 8. (As stressed by Doctor Fortran @SteveLionel, the Fortran standard imposes no relation between these kind numbers and the number of bytes used to represent the data. See Steve's blog Doctor Fortran in “It Takes All KINDs” for a nice elaboration on this topic.)

[Update: The above speculation about kind(ieee_is_nan(1.0)=kind(.false.) turns out wrong even though it is what the Fortran 2018 standard suggests. With certain compiler options, kind(ieee_is_nan(1.0) and kind(.false.) can actually differ from each other and hence violate the specifications in the Fortran standard. See the answer by @francescalus and my summary at the end of this question.]

[Update: in the original question, I flipped "4" and "8", so that I thought kind(ieee_is_nan(1.0)) was changed; indeed, it is kind(.false.) that is altered from 4 to 8, while kind(ieee_is_nan(1.0)) remains 4 always. So everything is fully explained by the nice answer from @francescalus, thank you!]

For comparison, here is the same code without interfaced with MATLAB, which prints "4" and "4" on the screen when compiled with gfortran. With the nagfor compiler, the numbers remain equal although become 3.

! test_kind.f90
! Tested by GNU Fortran (Ubuntu 9.3.0-17ubuntu1~20.04) 9.3.0

program test_kind

use ieee_arithmetic, only: ieee_is_nan
implicit none

write(*, *) kind(ieee_is_nan(1.0))  ! This prints a number on STDOUT (screen)
write(*, *) kind(.false.)   ! A benchmark, which should print the same number on STDOUT (screen)

end program test_kind

For your reference, below is the section on ieee_is_nan in the Fortran 2018 standard. It specifies that the ieee_is_nan returns a "default logical", which I guess should be identical to the type of the intrinsic .true. or .false. constant --- or did I get it wrongly?

17.11.13 IEEE_IS_NAN (X)
1 Description. Whether a value is an IEEE NaN.
2 Class. Elemental function.
3 Argument. X shall be of type real.
4 Restriction. IEEE_IS_NAN (X) shall not be invoked if IEEE_SUPPORT_NAN (X) has the value false.
5 Result Characteristics. Default logical.
6 Result Value. The result has the value true if the value of X is an IEEE NaN; otherwise, it has the value false.

It seems extraordinary to me that MEX can change the type of the default logical without taking care of ieee_is_nan. Maybe there exists an option of MEX that can rectify this behavior, but why should it be the default in the first place?

I tried the same code on more machines with other versions of MATLAB and Fortran compiler. The results are the same.

  1. MATLAB 9.7.0.1319299 (R2019b) Update 5 with GNU Fortran (GCC) 8.3.1 20191121 (Red Hat 8.3.1-5):

    kind(ieee_is_nan(1.0))= 4, kind(.false.) = 8

    The same compiler without interfacing with MATLAB:

    kind(ieee_is_nan(1.0)) = 4 = kind(.false.)

  2. MATLAB 9.5.0.1049112 (R2018b) Update 3 with GNU Fortran (Ubuntu 9.3.0-17ubuntu1~20.04) 9.3.0:

    kind(ieee_is_nan(1.0)) = 4, kind(.false.) = 8

    The same compiler without interfacing with MATLAB:

    kind(ieee_is_nan(1.0)) = 4 = kind(.false.)

  3. (On Windows 10) MATLAB 9.5.0.944444 (R2018b) with Intel(R) Visual Fortran Intel (R) 64 Compiler, Version 19.1.1.216 Build 20200306

    kind(ieee_is_nan(1.0)) = 4, kind(.false.) = 8

    The same compiler without interfacing with MATLAB:

    kind(ieee_is_nan(1.0)) = 4 = kind(.false.)


A summary after accepting @francescalus' answer

  1. It turns out that the mismatch between kind(.false.) and kind(ieee_is_nan(1.0)) comes from the gfortran option -fdefault-integer-8, which is adopted by MEX as the default. This option enforces gfortran to use 64-bit integer and 64-bit logical as the default kind, but does not change the returned kind of ieee_is_nan, even though the Fortran standard specifies that ieee_is_nan should return the default logical kind. This may be because ieee_is_nan is not really an intrinsic procedure but only a procedure from the intrinsic module ieee_arithmetic.

  2. Note that ifort (version ifort (IFORT) 2021.2.0 20210228) and nagfor (NAG Fortran Compiler Release 7.0 (Yurakucho) Build 7036) behave also in the way described above, the corresponding option being -i8 for both of them. So the compiler vendors agree that it is fine to break the consistency of some intrinsic modules when certain options are enforced. This is surprising to me. Fortunately, flang (under clang 7.1.0) follows the Fortran standard even with -fdefault-integer-8 imposed, keeping kind(is_ieee_nan(1.0)) == kind(.false.) --- so it is not a mission impossible.

  3. The NAG compiler nagfor warns about ieee_is_nan when -i8 is adopted, pointing out that they are incompatible; gfortran and ifort, however, keep absolutely silent even if you invoke them with -Wall -Wexta and -warn all respectively. This is even more surprising to me.

  4. Given these facts, I decide not to use ieee_is_nan but to implement my own is_nan and stay vigilant about (or away from) all the procedures provided by intrinsic modules. Otherwise, my package will fail to compile due to a mismatch of types if the users choose to enforce the 64-bit integer as the default (without thinking about the logical kind; why should they?). More seriously, MATLAB has already made such a choice for all its users without telling them.

  5. I am happy to see that my question leads to interesting discussions. Since some compilers seem to be in need of improvements regarding the problems discovered here (you may have different opinions), I made a post on Fortran Discourse about this topic. I hope that it will draw more attention from the community and someone will patch at least the open-source compilers like gfortran.

Thank you very much for any comments or criticism.

like image 573
Nuno Avatar asked Sep 05 '21 03:09

Nuno


2 Answers

By default MEX compiles with the gfortran option -fdefault-integer-8. The way gfortran handles this results in what you see.

Consider the non-MEX program

  use, intrinsic :: ieee_arithmetic, only : ieee_is_nan, ieee_support_nan
  implicit none

  if (.not.ieee_support_nan(1.)) error stop "Lack of support"
  print*, KIND(ieee_is_nan(1.)), KIND(.TRUE.)

end

compiled with/without -fdefault-integer-8.

This option makes default integer and logical variables 8 bytes wide (and have kind parameter 8).

This all makes sense. However, it seems that gfortran doesn't change the function result kind parameter for functions in the intrinsic module ieee_arithmetic with this option (this is reasonable in some way: it's just taking what the "pre-compiled" module files say, which is "return is logical(4)".

(If things seem confusing when you see that it does seem to use the correct kinds in another intrinsic module like iso_c_binding, note that this second intrinsic module isn't supported by a supplied module file. You will see the same behaviour for IEEE modules with the OpenMP intrinsic module which is also provided as a file.)

As a workaround, you can use LOGICAL(ieee_is_nan(1.)) to make things match again.


Compiler options like -fdefault-integer-8 are not toys. They should be used only in the rarest of rare circumstances. With gfortran, if you use this option by itself you are telling the compiler "Go ahead and do whatever you like: I don't care whether or not you treat my program as the Fortran standard specifies." In this case, this is a choice that MATLAB has made for you.

In particular, if you use a compiler option to change kind parameters of things (such as default kinds, or the kind numbering scheme as seen in the NAG compiler for example), you must compile all constituent parts of your program with that option to be sure of consistency. As seen in this case, that's not always easy, making these horribly dangerous options to have as defaults. Even intrinsic modules may be included in the list of things which need to be considered for consistency.

like image 179
francescalus Avatar answered Dec 03 '22 23:12

francescalus


(Edit - I misread the statement. See comment below.)

Certainly, these two numbers depend on the compiler --- they are 3 for the nagfor compiler, but I thought they should always be equal to each other.

NO! There is no implicit correspondence between the kind numbers of one type and the kind numbers of another type! Do not be fooled by the common use of byte sizes for kind numbers - this is just a convention some compilers adopted to look more familiar to programmers used to the integer*4 extension.

It would be perfectly valid for a compiler to have kinds 14, 27 and 830 for LOGICAL and 2,9 and 13 for REAL. The only things you can count on are that:

  • Default integer, default real and default logical each occupy one "numeric storage unit"
  • Double precision and default complex occupy two numeric storage units

(See Fortran 2018 19.5.3.2 Storage Sequence)

Of course, as already pointed out, if you use compiler options to change the default kind for one type but not another, you break this rule.

For further reading, see my blog post, Doctor Fortran in “It Takes All KINDs”

like image 38
Steve Lionel Avatar answered Dec 04 '22 00:12

Steve Lionel