I’ll have to look at the FITS documentation to verify if this is true, but my experience with writing C-interop interfaces for large projects is that a lot of the code in the C side are support routines. My standard practice is to first look at the Fortran language manual and see what routines are called from Fortran and implement those first. Any routines not supported by the native Fortran interfaces but look to be useful (these were normally new functions added to the C side that have never found there way into the Fortran interfaces) are added next. I haven’t examined the FITS code and documentation in detail but if I have time I’ll take a look.
To illustrate just what you would face if you tried to replace the interfaces generated by cfortran.h with an equivalent F2003 C-interop interface, this is what cfortran.h generates for the FTMKKY function that generates a “properly formatted keyword record”
extern void ftmkky_( char * A1 , char * A2 , char * A3 , char * A4 , int * A5 , size_t C1 , size_t C2 , size_t C3 , size_t C4 ); extern void ftmkky_( char * A1 , char * A2 , char * A3 , char * A4 , int * A5 , size_t C1 , size_t C2 , size_t C3 , size_t C4 ) { char *B1=
# 341 "f77_wrap1.c" 3 4
((void *)0)
# 341 "f77_wrap1.c"
; char *B2=
# 341 "f77_wrap1.c" 3 4
((void *)0)
# 341 "f77_wrap1.c"
; char *B3=
# 341 "f77_wrap1.c" 3 4
((void *)0)
# 341 "f77_wrap1.c"
; char *B4=
# 341 "f77_wrap1.c" 3 4
((void *)0)
# 341 "f77_wrap1.c"
; ; ffmkky( ( !(C1<4||A1[0]||A1[1]||A1[2]||A1[3]) ) ? ((char*)0) : memchr(A1,'\0',C1) ? A1 : ((B1=(char*)malloc(( (C1>gMinStrLen) ? C1 : gMinStrLen )+1))[C1]='\0',memcpy(B1,A1,C1), kill_trailing(B1,' ')) , ( !(C2<4||A2[0]||A2[1]||A2[2]||A2[3]) ) ? ((char*)0) : memchr(A2,'\0',C2) ? A2 : ((B2=(char*)malloc(( (C2>gMinStrLen) ? C2 : gMinStrLen )+1))[C2]='\0',memcpy(B2,A2,C2), kill_trailing(B2,' ')) , ( !(C3<4||A3[0]||A3[1]||A3[2]||A3[3]) ) ? ((char*)0) : memchr(A3,'\0',C3) ? A3 : ((B3=(char*)malloc(( (C3>gMinStrLen) ? C3 : gMinStrLen )+1))[C3]='\0',memcpy(B3,A3,C3), kill_trailing(B3,' ')) , ((B4=(char*)malloc(( (C4>gMinStrLen) ? C4 : gMinStrLen )+1))[C4]='\0',memcpy(B4,A4,C4), kill_trailing(B4,' ')) , A5 ); if (B1) free(B1); if (B2) free(B2); if (B3) free(B3); if (B4) memcpy(A4,B4, (strlen(B4)<C4?strlen(B4):C4)), (C4>strlen(B4)?memset(A4+strlen(B4),' ', C4-strlen(B4)):0), free(B4); return ; }
This is just one function in one of the four Fortran wrapper C files. I think this illustrates why you are better off just rewritting everything in pure Fortran if possible. The level of effort to replace the cfortran.h generated interfaces with C-interop code is probably the same as just creating your own library in pure Fortran.
1 Like
It looks like LLVM has an optimization pass which can recognize byte swapping: llvm::recognizeBSwapOrBitReverseIdiom, but I don’t know if that’s what fired here.
If you compile the function with flang -O3 -c -mllvm -print-changed=quiet it shows the optimization passes and how they modify the intermediate representation. The byte swap instruction appears in the 15th or 16th pass that combines instructions:
*** IR Dump After InstCombinePass on le_to_be_ ***
define double @le_to_be_(ptr noalias nofree readonly captures(none) %0) local_unnamed_addr #0 {
%2 = load i64, ptr %0, align 8, !tbaa !4
%3 = call i64 @llvm.bswap.i64(i64 %2)
%4 = bitcast i64 %3 to double
ret double %4
}
The gfortran f951 Fortran frontend seems to be missing the registration of the BUILTIN_IN_BSWAP64 function: Code search results · GitHub
3 Likes
Ok, I experimented with this and the results are pretty interesting. Depending on how the fortran is written (inline code, do loop, array syntax), and the optimization options, compilers do and don’t reduce it down to that single instruction. The intel ifx compiler seems to do pretty well in general, some other compilers sometimes recognize the trick and sometimes don’t, and other compilers never see the trick. The arm64 architecture (called aarch64 here) also has 2-byte, 4-byte, and 8-byte swap instructions, but I could not get any of the fortran compilers to generate the 8-byte swap instruction for any of the code+option combinations that I tried. That was a little disappointing since that is a RISC architecture, and I would have expected the move-to-register, operate on register, move from register kind of operation to compile close to the hardware, but it doesn’t. I was also disappointed in the other LLVM compilers in general. I think they all should have the capability to recognize that reduction down to a single instruction as a common back-end optimization. Maybe I’m missing some critical compiler option? Also, I do not see the LFortran+arm64 option there – hopefully that can be added in the future, especially if the LFortran developers can get their compiler to work with MacOS Homebrew to make it easy to install, use, and upgrade. All in all though, this is an interesting exercise with the compiler explorer site.
2 Likes
I asked on the Compiler Explorer Discord, and they pointed out that using the --target=<value>/-target <value> flag should work with clang/flang. Currently it breaks the assembly syntax highlighting and something breaks when iso_fortran_env is used. After replacing the named constants real64 and int64 with the value 8 and adding --target=aarch64 I got this from flang-trunk (version 23):
le_to_be_:
ldr x8, [x0]
rev x8, x8
fmov d0, x8
ret
which uses the REV instruction.
Apologies if we have derailed the discussion from the FITS IO topic.
1 Like
Is that a feature or a bug? I am aware that the USE of some modules (e.g. ieee_features) can change the code generated by a compiler, even when nothing is referenced within the module. Does that principle apply also to the intrinsic iso_fortran_env module? If so, is that behavior affected by the only: clause, or is it triggered just by the USE statement?
FMU, Compiler Explorer is served from Amazon Web Services using x86-64 (virtual) machines. When you add the --target=aarch64 flag you are cross-compiling for Arm, but the (architecture-specific) flang Fortran runtime library (the libflang-rt archive) which includes the iso_fortran_env module is missing. The friendly folks at the CE Discord had the following to say:
“I think runtime libraries are somewhat less useful if not compiling for current arch”
“if its any consolation, we dont support clang libc++ for aarch64 either”
I guess it’s a “feature” that flang ships the iso_fortran_env.mod in the runtime library. But as long as your Fortran snippets don’t use any intrinsic modules, the -target flag should work if you just want to look at Arm 64-bit (or any other) assembly.
1 Like
I had a look at the file structure (using the original description) and found that it is fairly simple to read, at least if you limit the capabilities a bit. I used a sample file UITfuv2582gc.fits from the cited website. I have attached the program (very limited capabilities, but it can easily be extended 
). It merely demonstrates how you could read the data. Nothing beyond that.
readheader.f90 (3.7 KB)
2 Likes
Thanks for this example!
However, what can be the strategy to use the same function and the main program to work with different types of files? In IDL or python the type of an array can be changed on the fly, so data=readfits(filename.fits) will provide the type that corresponds to the FITS file. What will be the best practice us here?
Off-hand, since Fortran uses static types most of the time, you would probably want to convert the data to such a fixed type. 32-bits integers can hold all the data I have seen in the description. As for dimensionality, I would say support only the ones likely to turn up in practice. I have once used a five-dimensional array of data, but never more than that. And certainly not a 999-dimensional one.
A suggestion: use a routine to find out the dimensionality and branch on that, passing an array of the right dimension to the routine that retrieves the data. More sophisticated approaches are possible, but they require more text than I can type at the moment 
.
1 Like
No, they can not. BITPIX values of 64 (integers) and -64 (64-bit FP) are perfectly valid.
Also, although there is no native support for unsigned integers, FITS uses appropriate BSCALE (1.0) and BZERO keywords for images (TSCALEn, TZEROn for binary tables) to shift the unsigned value into the respective signed range. Yet, when read in, they should be shifted back and treated as unsigned.
Add image, ascii-table, binary-table extensions and several data compression algorithms (both lossless and not) to that and you get a pretty complicated world of data. Writing a full (or even partial) equivalent of CFITSIO in Fortran from scratch will be a long trip. I guess that creating modern interop interface to C might be a better choice. Going from the core functions up to more general ones.
Ah, then my information was incomplete . Thanks for this correction.