Try module for scaling constants

[weslleyspereira]

This PR adds a module for scaling constants proposed by https://doi.org/10.1145/3061665. It was modified following the discussion in #496.

• Since the codes were written in the Fortran90 standard, I removed

# Tell CMake that our Fortran sources are written in fixed format.
set(CMAKE_Fortran_FORMAT FIXED)

from the CMakeLists.txt. The compiler now determines the Fortran layout by the file extension.

[codecov]

Codecov Report

Merging #501 (5b73654) into master (8bc7b0a) will not change coverage.
The diff coverage is n/a.

@@           Coverage Diff           @@
##           master     #501   +/-   ##
=======================================
  Coverage   83.33%   83.33%           
=======================================
  Files        1820     1820           
  Lines      170857   170857           
=======================================
  Hits       142384   142384           
  Misses      28473    28473           

Impacted Files	Coverage Δ
SRC/dlasq2.f	76.92% <0.00%> (ø)
SRC/slasq2.f	76.61% <0.00%> (ø)

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 8bc7b0a...5b73654. Read the comment docs.

[zerothi]

I made a PR against your branch. Also, at some point it may be nice to have an explicit quad precision, for the record I put it here:

   !  Standard constants for 
   integer, parameter :: qp = selected_real_kind(p=33)

   real(qp), parameter :: qzero = 0.0_qp
   real(qp), parameter :: qhalf = 0.5_qp
   real(qp), parameter :: qone = 1.0_qp
   real(qp), parameter :: qtwo = 2.0_qp
   real(qp), parameter :: qthree = 3.0_qp
   real(qp), parameter :: qfour = 4.0_qp
   real(qp), parameter :: qeight = 8.0_qp
   real(qp), parameter :: qten = 10.0_qp
   complex(qp), parameter :: uzero = ( 0.0_qp, 0.0_qp )
   complex(qp), parameter :: uhalf = ( 0.5_qp, 0.0_qp )
   complex(qp), parameter :: uone = ( 1.0_qp, 0.0_qp )
   character*1, parameter :: qprefix = 'Q'
   character*1, parameter :: uprefix = 'U'

!  Scaling constants
   real(qp), parameter :: qulp = epsilon(0._qp)
   real(qp), parameter :: qeps = qulp * 0.5_qp
   real(qp), parameter :: qsafmin = real(radix(0._qp),qp)**max( &
      minexponent(0._qp)-1, &
      1-maxexponent(0._qp) &
   )
   real(qp), parameter :: qsafmax = qone / qsafmin
   real(qp), parameter :: qsmlnum = qsafmin / qulp
   real(qp), parameter :: qbignum = qsafmax * qulp
   real(qp), parameter :: qrtmin = sqrt(qsmlnum)
   real(qp), parameter :: qrtmax = sqrt(qbignum)

!  Blue's scaling constants
   real(qp), parameter :: qtsml = real(radix(0._qp), qp)**ceiling( &
       (minexponent(0._qp) - 1) * 0.5_qp)
   real(qp), parameter :: qtbig = real(radix(0._qp), qp)**floor( &
       (maxexponent(0._qp) - digits(0._qp) + 1) * 0.5_qp)
!  ssml = 1/s, where s was defined in https://doi.org/10.1145/355769.355771
   real(qp), parameter :: qssml = real(radix(0._qp), qp)**( - floor( &
       (minexponent(0._qp) - 1) * 0.5_qp))
!  ssml = 1/S, where S was defined in https://doi.org/10.1145/355769.355771
   real(qp), parameter :: qsbig = real(radix(0._qp), qp)**( - ceiling( &
       (maxexponent(0._qp) - digits(0._qp) + 1) * 0.5_qp))

however, note that compilers/architectures which doesn't support this will error out at compile time. So probably not the best to put in... :(
It at least requires some build magic.

[weslleyspereira]

I made a PR against your branch.

I don't see your PR. In any case, I shall mention (and you've probably already noticed) that I used

integer, parameter :: sp = kind(1.e0)

and

integer, parameter :: dp = kind(1.d0)

instead of the selected_real_kind(p=...) as you suggested. The reason is to keep real(sp) = real and real(dp) = double precision regardless the use of flags
-freal-4-real-8
-freal-8-real-16
-fdefault-real-8
This was based on the results of your code snippet #496 (comment) in my computer:

           4
    single  4  4  4  4
    double  8  8  8  8  8
   1.19209290E-07   2.2204460492503131E-016

-freal-4-real-8
           8
    single  8  8  8  8
    double  8  8  8  8  8
   2.2204460492503131E-016   2.2204460492503131E-016

-freal-8-real-16
           4
    single  4  4  4  4
    double 16 16 16 16 16
   1.19209290E-07   1.92592994438723585305597794258492732E-0034

-freal-8-real-16 -freal-4-real-8
           8
    single  8 16 16 16
    double 16 16 16 16 16
   2.2204460492503131E-016   1.92592994438723585305597794258492732E-0034

-fdefault-real-8
           8
    single  8  4  4  4
    double 16  8  8 16  8
   2.2204460492503131E-016   1.92592994438723585305597794258492732E-0034

Notice that the flags we mentioned do not work with Lapack.

[zerothi]

See here weslleyspereira#1

[zerothi]

Looks good to me!

[zerothi]

For the record I files a bug report with gcc about the inconsistent handling of the -freal-X-real-Y handling.

From 10.3 onwards (fix merged!) one should be able to compile without promoting everything to quad precision, i.e. having a double-quad library. See here.

EDIT: bug fixed in gcc!