# Defining multiple real precisions in a program

**URL:** <https://fortran-lang.discourse.group/t/defining-multiple-real-precisions-in-a-program/2550>\
**Category:** Uncategorized\
**Created:** [January 9, 2022, 3:59pm UTC](https://fortran-lang.discourse.group/t/defining-multiple-real-precisions-in-a-program/2550 "2022-01-09T15:59:11Z")\
**Posts on this page:** 16\
**Page:** 1

<div class="post-metadata">

**Author:** ![Beliavsky](https://avatars.discourse-cdn.com/v4/letter/b/ba8739/32.png) [@Beliavsky](https://fortran-lang.discourse.group/u/Beliavsky)\
**Post date:** [January 9, 2022, 3:59pm UTC](https://fortran-lang.discourse.group/t/defining-multiple-real-precisions-in-a-program/2550/1 "2022-01-09T15:59:11Z")

</div>

To define multiple real precisions in a program, which are relative to each other so that some calculations can be done in a higher precision than others, I came up with

```auto
program kinds
implicit none
integer, parameter :: qp_extra_p = 10 ! set to 1 to get qp = 10 for gfortran
integer, parameter :: sp = kind(1.0) ! could be sp = kind(1.0d0)
integer, parameter :: dp = selected_real_kind(p=1+precision(1.0_sp),r=1+exponent(1.0_sp))
integer, parameter :: qp = selected_real_kind(p=qp_extra_p+precision(1.0_dp),r=1+exponent(1.0_dp))
real(kind=sp) :: xsingle
real(kind=dp) :: xdouble
real(kind=qp) :: xquad ! will not compile if sp = kind(1.0d0)
print*," kind",sp,dp,qp
print*,"precision",precision(xsingle),precision(xdouble),precision(xquad)
print*," range",range(xsingle),range(xdouble),range(xquad)
end program kinds

```

Output with gfortran, Intel Fortran, and g95:

```auto
      kind 4 8 16
 precision 6 15 33
     range 37 307 4931

```

Are there better ways to do it?

---

<div class="post-metadata">

**Author:** ![Beliavsky](https://avatars.discourse-cdn.com/v4/letter/b/ba8739/32.png) [@Beliavsky](https://fortran-lang.discourse.group/u/Beliavsky)\
**Post date:** [January 9, 2022, 5:51pm UTC](https://fortran-lang.discourse.group/t/defining-multiple-real-precisions-in-a-program/2550/2 "2022-01-09T17:51:06Z")

</div>

Using the method above of setting kinds, I have written a version of [norm2](https://gcc.gnu.org/onlinedocs/gfortran/NORM2.html) that appears to handle overflow better than Intel Fortran does.

```auto
module norm_mod
implicit none
private
public :: wp, my_norm2, my_norm2_extra
integer, parameter :: wp = kind(1.0)
integer, parameter :: wp2 = selected_real_kind(p=10+precision(1.0_wp),r=1+exponent(1.0_wp))
contains
pure function my_norm2(x) result(y)
real(kind=wp), intent(in) :: x(:)
real(kind=wp) :: y
y = sqrt(sum(x**2))
end function my_norm2
!
pure function my_norm2_extra(x) result(y)
! use extra precision to calculate L2 norm
real(kind=wp), intent(in) :: x(:)
real(kind=wp) :: y
y = sqrt(sum(real(x,kind=wp2)**2))
end function my_norm2_extra
end module norm_mod
!
program test_norm2
use norm_mod, only: wp, my_norm2, my_norm2_extra
implicit none
integer, parameter :: n = 1000, nhuge = 100
real(kind=wp) :: huge_vec(n)
huge_vec = 0.0_wp
huge_vec(:nhuge) = huge(0.0_wp)/nhuge
print*," true",huge(0.0_wp)/sqrt(real(nhuge,kind=wp))
print*,"intrinsic",norm2(huge_vec)
print*," naive",my_norm2(huge_vec)
print*," extra",my_norm2_extra(huge_vec)
end program test_norm2

```

Intel Fortran output:

```auto
      true 3.4028235E+37
 intrinsic Infinity
     naive Infinity
     extra 3.4028235E+37

```

gfortran output:

```auto
      true 3.40282347E+37
 intrinsic 3.40282347E+37
     naive Infinity
     extra 3.40282347E+37

```

Similar results are obtained if the program uses

`wp = kind(1.0d0)`

---

<div class="post-metadata">

**Author:** ![Arjen](https://avatars.discourse-cdn.com/v4/letter/a/b9bd4f/32.png) [@Arjen](https://fortran-lang.discourse.group/u/Arjen)\
**Post date:** [January 10, 2022, 8:12am UTC](https://fortran-lang.discourse.group/t/defining-multiple-real-precisions-in-a-program/2550/3 "2022-01-10T08:12:06Z")

</div>

You might want to guard against a compilation failure on qp by:

```auto
integer, parameter :: qptry = selected_real_kind( ... )
integer, parameter :: qp = merge( qptry, qp, qptry > 0 )

```

Also: you increase both the precision and the exponent range. I can imagine that it is not always desirable - you might want to increase precision and leave the range as is. Though probably there are only few systems where these two parameters are disconnected.

---

<div class="post-metadata">

**Author:** ![Beliavsky](https://avatars.discourse-cdn.com/v4/letter/b/ba8739/32.png) [@Beliavsky](https://fortran-lang.discourse.group/u/Beliavsky)\
**Post date:** [January 10, 2022, 12:36pm UTC](https://fortran-lang.discourse.group/t/defining-multiple-real-precisions-in-a-program/2550/4 "2022-01-10T12:36:55Z")

</div>

Thanks for the guard idea. I think it needs to modified to handle the case that `qptry` is -1.

```auto
implicit none
integer, parameter :: dp = selected_real_kind(15, 307)
integer, parameter :: qptry = selected_real_kind(133, 4931) ! -1 for non-existent KIND
integer, parameter :: qp = merge(qptry, dp, qptry > 0) ! use qptry if available, otherwise dp
print*,precision(1.0_qp) ! 15
end

```

---

<div class="post-metadata">

**Author:** ![Arjen](https://avatars.discourse-cdn.com/v4/letter/a/b9bd4f/32.png) [@Arjen](https://fortran-lang.discourse.group/u/Arjen)\
**Post date:** [January 10, 2022, 1:08pm UTC](https://fortran-lang.discourse.group/t/defining-multiple-real-precisions-in-a-program/2550/5 "2022-01-10T13:08:47Z")

</div>

Yes, of course - I did not check the code properly. This is indeed what I had intended.

---

<div class="post-metadata">

**Author:** ![Jcollins](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/jcollins/32/540_2.png) [@Jcollins](https://fortran-lang.discourse.group/u/Jcollins)\
**Post date:** [January 25, 2022, 7:10pm UTC](https://fortran-lang.discourse.group/t/defining-multiple-real-precisions-in-a-program/2550/6 "2022-01-25T19:10:04Z")

</div>

A few years ago we looked at a number format where the precision is traded off against the exponent . So, for example, if the exponent fits in 4 bits we have more bits for precision. The advantages are:  
i The representation is more even in log space;  
ii. The data size can be kept down for interprocessor communication.  
But  
SELECTED\_REAL\_KIND gets more interesting.

Has this been investigated elsewhere?

John

---

<div class="post-metadata">

**Author:** ![Arjen](https://avatars.discourse-cdn.com/v4/letter/a/b9bd4f/32.png) [@Arjen](https://fortran-lang.discourse.group/u/Arjen)\
**Post date:** [January 26, 2022, 7:37am UTC](https://fortran-lang.discourse.group/t/defining-multiple-real-precisions-in-a-program/2550/7 "2022-01-26T07:37:44Z")

</div>

I have never seen anything of the sort ;). Something that is of considerable interest in the data science community is half-precision reals, because they have little to gain with six decimals but everything with less memory per number. That would, however, still fit in a more or less ordered sequence of number formats - the precision is reduced along with the range. Your example would upset that.

I have seen a Convex machine in years long gone by where there were two number formats. The main difference being a different interpretation of the exponent encoding and a slightly faster processing with the native number format (that was not IEEE-compliant). But that was before Fortran 90 was a serious option. You had to select the format via a compiler option.

Ads it happens there is a thread on the gfortran mailing list that concerns the handling of signalling NaNs. A large number of platforms are mentioned there each with its own ideas about what a number should be (\*). Perhaps one of them has something akin to the number format you mention.

(\*) Is it my too poor mastering of the English language that I have to use the word “number” in at least two different meanings?

---

<div class="post-metadata">

**Author:** ![JohnCampbell](https://avatars.discourse-cdn.com/v4/letter/j/5daacb/32.png) [@JohnCampbell](https://fortran-lang.discourse.group/u/JohnCampbell)\
**Post date:** [January 26, 2022, 8:18am UTC](https://fortran-lang.discourse.group/t/defining-multiple-real-precisions-in-a-program/2550/8 "2022-01-26T08:18:09Z")

</div>

> [@Arjen](#):
>
> Also: you increase both the precision and the exponent range. I can imagine that it is not always desirable - you might want to increase precision and leave the range as is. Though probably there are only few systems where these two parameters are disconnected.

I presume these “few systems” would not be IEEE-compliant.  
Do any of these who do support modern Fortran actually exist ?

---

<div class="post-metadata">

**Author:** ![Arjen](https://avatars.discourse-cdn.com/v4/letter/a/b9bd4f/32.png) [@Arjen](https://fortran-lang.discourse.group/u/Arjen)\
**Post date:** [January 26, 2022, 8:34am UTC](https://fortran-lang.discourse.group/t/defining-multiple-real-precisions-in-a-program/2550/9 "2022-01-26T08:34:52Z")

</div>

I honestly would not know. I guess most of the ones I knew about or heard about have crumbled to dust, but the GCC compiler suite, hence gfortran, is tested for a large range of devices.

---

<div class="post-metadata">

**Author:** ![Jcollins](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/jcollins/32/540_2.png) [@Jcollins](https://fortran-lang.discourse.group/u/Jcollins)\
**Post date:** [January 26, 2022, 9:26am UTC](https://fortran-lang.discourse.group/t/defining-multiple-real-precisions-in-a-program/2550/10 "2022-01-26T09:26:01Z")

</div>

I think Itanium VMS Fortran can run in D-format. The real and double precision numbers have the same exponent format (and therefore range). Double precision just has an additional 32 bits of precision at the low precision end. A side effect is that you can read a single precision value from the address of a double precision variable. This was common practice and makes migration to IEEE exciting.

Gould-SEL numbers have the same characteristic but the exponent is HEX. But the Fortran is extended FORTRAN 77.

---

<div class="post-metadata">

**Author:** ![nncarlson](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/nncarlson/32/1175_2.png) [@nncarlson](https://fortran-lang.discourse.group/u/nncarlson)\
**Post date:** [January 26, 2022, 7:42pm UTC](https://fortran-lang.discourse.group/t/defining-multiple-real-precisions-in-a-program/2550/11 "2022-01-26T19:42:27Z")

</div>

With respect to variable mantissa and exponent ranges, I’ll point out the unum number format promoted in recent years by John Gustafson. The [wikipedia page](https://en.wikipedia.org/wiki/Unum_(number_format)) is a decent place to start.

---

<div class="post-metadata">

**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\
**Post date:** [January 26, 2022, 8:44pm UTC](https://fortran-lang.discourse.group/t/defining-multiple-real-precisions-in-a-program/2550/12 "2022-01-26T20:44:46Z")

</div>

> [@Jcollins](#):
>
> A few years ago we looked at a number format where the precision is traded off against the exponent . So, for example, if the exponent fits in 4 bits we have more bits for precision.

In my field of work, there was a publication recently which investigated something along these lines but for 16-bit formats:

> **[On the accuracy and performance of the lattice Boltzmann method with 64-bit,...](https://arxiv.org/abs/2112.08926)**
>
> Fluid dynamics simulations with the lattice Boltzmann method (LBM) are very memory-intensive. Alongside reduction in memory footprint, significant performance benefits can be achieved by using FP32 (single) precision compared to FP64 (double)...

---

<div class="post-metadata">

**Author:** ![JohnCampbell](https://avatars.discourse-cdn.com/v4/letter/j/5daacb/32.png) [@JohnCampbell](https://fortran-lang.discourse.group/u/JohnCampbell)\
**Post date:** [January 27, 2022, 6:25am UTC](https://fortran-lang.discourse.group/t/defining-multiple-real-precisions-in-a-program/2550/13 "2022-01-27T06:25:32Z")

</div>

I am surprised that you are investigating a 16-bit real to look for an improved performance.

I have no idea of the hardware you have available, but improved efficiency of vector instructions must be preferable to software emulation of smaller memory real calculations, given the ratios are only 50% change of memory.

I would have thought that utilising avx-256 or avx-512 could be more effective, combined with targeting the L1 cache, which is so important for efficient AVX. Surely a software real(2) would not achieve the vector real(4) performance.  
The alternative is to utilise a processor with larger cache and increased memory bandwidth. I have been trying to understand the “black art” AVX inefficiency for a few years now.  
My latest attempt with Ryzen 5900X: more cache and faster memory was only moderately successful, but I am hoping that newer hardware and DDR5 memory might show an improvement.  
For me, multi thread computation has extra cores, but they share the same limited memory addressing capacity/bandwidth. Am I now arguing for smaller reals ?

---

<div class="post-metadata">

**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\
**Post date:** [January 27, 2022, 10:15am UTC](https://fortran-lang.discourse.group/t/defining-multiple-real-precisions-in-a-program/2550/14 "2022-01-27T10:15:33Z")

</div>

I am not involved directly, just following what’s going on in the field.

The reasons why they are looking at 16-bit reals are several:

- the application is extremely memory intensive; for each cell of a 3D spatial grid, you need to store 19 or even 27 DDF’s (discrete distribution functions), plus any additional scalar fields (temperature, pressure, etc.)
- given how memory intensive it is, the available GPU memory quickly becomes a restriction on the domain size (or grid size); this also puts a cap on the Reynolds number you can simulate (without resorting to turbulence models)
- the algorithm is bandwith-limited; a large share of time is spent just shifting the DDF’s in memory (think particles moving along a grid)
- It turns out with some tricks, and assuming a stable simulation, the DDF’s will fall between -2 and 2, mostly aggregating around 0. Parts of the mantissa and exponent are therefore completely unused.

As far I can understand, the application was written in OpenCL, for hardware that had support for 16-bit floats ([cl\_khr\_fp16](https://www.khronos.org/registry/OpenCL/sdk/2.2/docs/man/html/cl_khr_fp16.html)), and not in Fortran.

---

<div class="post-metadata">

**Author:** ![lkedward](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/lkedward/32/72_2.png) [@lkedward](https://fortran-lang.discourse.group/u/lkedward)\
**Post date:** [January 27, 2022, 10:58am UTC](https://fortran-lang.discourse.group/t/defining-multiple-real-precisions-in-a-program/2550/15 "2022-01-27T10:58:34Z")

</div>

> [@JohnCampbell](#):
>
> I am surprised that you are investigating a 16-bit real to look for an improved performance.
> 
> I have no idea of the hardware you have available…

As Ivan points out, any way to exploit half precision can give potential benefits on modern GPUs which now seem to prioritise FP16 performance for machine learning applications. In fact any way to exploit lower precision on GPUs is beneficial due to the financial cost now required for full double precision hardware support.

I did some [simple investigations](https://www.researchgate.net/publication/348239254_Optimisation_of_a_Finite-Volume_Test-bench_Code_for_Highly_Parallel_Architectures) into mixed-precision solvers for CFD on GPUs and got a 2x speedup using FP32 by alleviating the effective memory bottleneck.

Some further reading on the topic:

- [Investigating Half Precision Arithmetic to Accelerate Dense Linear System Solvers (utk.edu)](https://www.icl.utk.edu/files/publications/2017/icl-utk-985-2017.pdf)

- [Mixed-precision AMG as linear equation solver for definite systems - ScienceDirect](https://www.sciencedirect.com/science/article/pii/S1877050910000219)

- [Full article: Performance and accuracy of hardware-oriented native-, emulated- and mixed-precision solvers in FEM simulations (tandfonline.com)](https://www.tandfonline.com/doi/full/10.1080/17445760601122076)

---

<div class="post-metadata">

**Author:** ![oscardssmith](https://avatars.discourse-cdn.com/v4/letter/o/b9e5f3/32.png) [@oscardssmith](https://fortran-lang.discourse.group/u/oscardssmith)\
**Post date:** [January 27, 2022, 1:51pm UTC](https://fortran-lang.discourse.group/t/defining-multiple-real-precisions-in-a-program/2550/16 "2022-01-27T13:51:48Z")

</div>

One place where 16 bit can be really useful is if you are doing iterative algorithms. Running the first few iterations in a low precision type can get you a good approximation that you can then refine in higher precision.
