2

I am writing a mex gateway for a piece of Fortran code.

In the Fortran code, for portability, the floating-point variables are declared as 

REAL(kind(0.0D0)) :: x, y, etc

(BTW, I am aware that there are better ways to do it, as discussed at Fortran: integer*4 vs integer(4) vs integer(kind=4), What does "real*8" mean?, and https://software.intel.com/en-us/blogs/2017/03/27/doctor-fortran-in-it-takes-all-kinds )

However, it seems to me that mex supports only REAL*8 and REAL*4, the former being Double, the latter being Single. I got this impression from the following functions/subroutines:

mxIsDouble, mxIsSingle, mxCopyPtrToReal8, mxCopyReal8ToPtr, mxCopyPtrToReal4, mxCopyReal4ToPtr

My questions are as follows.

  1. Is it true that mex supports only REAL*8 and REAL*4?

  2. Does it improve the portability of the mex gateway if I declare double-precision floating-point variables as

    REAL(kind(0.0D0)) :: x, y, etc

or even

integer, parameter :: dp = selected_real_kind(15, 307)
real(kind=dp) :: x, y, etc

Or should I simply declare 

REAL*8 :: x, y, etc

  1. Are REAL*8 and/or REAL*4 supported on all platforms? If no, does this mean that MATLAB mex is intrinsically unportable?

  2. What is the best way to specify the kind of floating-point variables in mex gateways for Fortran code?

The following code is an example. See the declaration of x, y, and xs. 

#include "fintrf.h"

      subroutine mexFunction(nlhs, plhs, nrhs, prhs)
C     y = square (x) 
C     x: a floating point scalar
C     y: x^2

      implicit none

C     mexFunction arguments 
      integer, intent(in) :: nlhs, nrhs 
      mwPointer, intent(in) :: prhs(nrhs)
      mwPointer, intent(inout) :: plhs(nlhs)

C     function declarations:
      mwPointer, external :: mxCreateDoubleScalar, mxGetPr
      mwSize, external :: mxGetM, mxGetN
      integer*4, external :: mxIsDouble, mxIsSingle 

C     variables
      mwSize, parameter :: mwOne = 1
      integer, parameter :: dKind = kind(0.0D0) 
      integer, parameter :: sKind = kind(0.0) 
      real(kind=dKind) :: x, y ! Does this improve the portablity?
      real(kind=sKind) :: xs ! Does this improve the portablity?

C     validate number of arguments
      if (nrhs .ne. 1) then
         call mexErrMsgIdAndTxt ('mex:nInput', '1 input required.')
      endif
      if (nlhs .gt. 1) then
         call mexErrMsgIdAndTxt ('mex:nOutput', 'At most 1 output.')
      endif

C     validate input
      if (mxIsDouble(prhs(1)) .ne. 1 .and. mxIsSingle(prhs(1)) .ne. 1) 
      ! What if the input is a floating point number but neither Double nor Single? 
     + then 
          call mexErrMsgIdAndTxt ('mex:Input', 'Input a real number.')
      endif
      if (mxGetM(prhs(1)) .ne. 1 .or. mxGetN(prhs(1)) .ne. 1) then
          call mexErrMsgIdAndTxt ('mex:Input', 'Input a scalar.')
      endif

C     read input
      if (mxIsDouble(prhs(1)) .eq. 1) then
         call mxCopyPtrToReal8(mxGetPr(prhs(1)), x, mwOne) 
      else
         call mxCopyPtrToReal4(mxGetPr(prhs(1)), xs, mwOne)
         x = real(xs, dKind) 
         ! What if the input is a floating point number but neither REAL*8 nor REAL*4
      endif

C     do the calculation
      y = x**2

C     write output
      plhs(1) = mxCreateDoubleScalar(y)

      return
      end subroutine mexFunction

The code runs correctly. Yet I am not sure whether it is portable.

Nuno
  • 256
  • 1
  • 11
  • Maybe also `real*4 single; real(kind(single)) x`. Just as unportable, but limiting the number of parts to worry about. – francescalus Jul 28 '19 at 20:56
  • 1
    "does this mean that MATLAB mex is intrinsically unportable?". MATLAB is available only for 3 OSes, and [the documentation lists exactly which version of which compilers you can use to make MEX-files](https://www.mathworks.com/support/requirements/supported-compilers.html). So yes, MEX is intrinsically unportable, in that you cannot compile MEX-files on other platforms or with other compilers. – Cris Luengo Jul 29 '19 at 04:59
  • Thank you, @francescalus and Cris Luengo, for the informative comments! – Nuno Jul 29 '19 at 05:26

2 Answers2

3

REAL*4 and REAL*8 are non-standard and non-portable. REAL(KIND(0.0D0) gets you DOUBLE PRECISION on every platform, as this is required by the Fortran standard.

I can't speak to MEX gateways, but you should avoid obvious non-standard features.

A popular choice is to define a module that declares named (PARAMETER) constants for the kinds in use. For example:

module kinds integer, parameter :: SP = KIND(0.0) integer, parameter :: DP = KIND(0.0D0) end module kinds

Then you can use SP and DP as kind values. If you ever need to change these, just edit the module.

Steve Lionel
  • 6,972
  • 18
  • 31
  • Thank you, Steve. May I ask two more questions? Q1. Is REAL(DP) equivalent to REAL*8 and REAL(SP) equivalent to REAL*4 on ALL platforms? Q2. How should I declare INTEGER*4 in a similar way? I ask these questions because mex explicitly requires REAL*8, REAL*4, and INTEGER*4. Thank you very much! – Nuno Jul 29 '19 at 01:45
  • 1
    I would say that whoever wrote those "requirements" needs to be sent in for reeducation. Seriously, you can interpret "REAL*4" as "single precision" and "REAL*8" as "double precision". On compilers that support those extensions, that's what they mean. As for INTEGER*4, you could use KIND(0), though I might actually recommend use of the standard intrinsics SELECTED_INT_KIND and SELECTED_REAL_KIND, with appropriate values for arguments. – Steve Lionel Jul 29 '19 at 17:14
  • 1
    Thanks again, @Steve Lionel. I totally agree with the reeducation --- these nonstandard things cause me a lot a headache. Sorry to bother you with a final question about the "appropriate values for arguments": Is it correct that SELECTED_INT_KIND(7), SELECTED_INT_KIND(8), SELECTED_INT_KIND(9) will all give us INTEGER*4, SELECTED_INT_KIND(6, 37) will give us REAL*4, and SELECTED_INT_KIND(15,307) will give us REAL*8? I read these values from elsewhere, yet it would be great to hear opinions from a true expert. Thank you very much! – Nuno Jul 30 '19 at 00:34
  • 1
    I assume you meant to say SELECTED_REAL_KIND for the real types. Those are the values I would use (9 for INT), though I know of a platform (in an unusual configuration) where (15,307) would not get what you want. An alternative (though I don't like it) is the INT32, REAL32 and REAL64 constants from intrinsic module ISO_FORTRAN_ENV. (I have some philosophical objections to how these are often used, but they work in your case.) – Steve Lionel Aug 01 '19 at 06:35
  • If I understand correctly, `KIND` was introduced in Fortran 90. I'm pretty sure that most of the Fortran code in MATLAB is from before then. When the `REAL*8` code was written, it was the correct way of doing so. I would even bet that, besides Cleve himself, there's no Fortran programmers working at the MathWorks, and there haven't been any in a loooong time. – Cris Luengo Sep 11 '19 at 16:47
  • 1
    `REAL*8` was never standard. The pre-F90 way of expressing that was `DOUBLE PRECISION` – Steve Lionel Sep 12 '19 at 00:29
2

Currently, it makes no difference whether you define variables as REAL*8/REAL*4 or REAL(REAL64)/REAL(REAL32). In the future MathWorks may come around and rewrite their functions to use portable variable declarations, but in my opinion this is unlikely for many reasons.

If you look in the fintrf.h file (included in every Fortran MEX gateway source file), you'll see that all of the MEX-specific procedures are defined with "asterisk notation," e.g. # define MWPOINTER INTEGER*8. So even if you define all of your variables with kinds from iso_fortran_env or selected_real_kind, any time you use a MathWorks variable type you're still using "asterisk notation" types, unless you go through that header file and redefine every symbol using your chosen kind specification.

zaen
  • 326
  • 3
  • 14