Lett's look at Pierre's problem again:
Pierre Hugonnet wrote:
> ... my [Pierre's] problem:
>
> INTERFACE mysum
> MODULE PROCEDURE mysum_single, mysum_double
> END INTERFACE
>
> ...
>
> SUBROUTINE mysum_single(a)
> REAL :: a(:)
> DOUBLE PRECISION :: accu
>
> accu = 0.0
> DO i=1,size(a)
> accu = accu + a(i)
> END DO
> mysum_single = accu
>
> END SUBROUTINE
>
> I want to process data with 6 digits accuracy: in the main I declare
> REAL(selected_real_kind(6)) :: array(n)
>
> On a 32bits machine where REAL is 6 digits and DOUBLE 12 digits,
> array(:) is REAL:
> the routine sums array using a 12 digits accumulator to reduce
> round-off problems. This is exactly what I want.
>
> On a 64bits machine where REAL would be 15 digits and DOUBLE 30 digits,
> array(:) is still REAL:
> the routines sums array using a 30 digits accumulator. I will
> give a good result, but 30 digits is much more than I actually
> need, leading to wasted CPU time.
>
> The solution could be:
> REAL(selected_real_kind(p=12)) :: accu
>
> But now It doesn't work if I process data with 12 digits
> accuracy (REAL(selected_real_kind(12)) :: array(n)) :
> I really need a 24 digits accumulator.
But you have one. On your machine (64 bits), it's 30 digits.
You requested single precision (12 digits), but on a 64-bit machine
gives about 15 digits. That's all you can get.
To request 24 digits, you get 30 digits. That's all you can get.
It's what is provided by the hardware. There are no other options.
It's not a Fortran issue at all.
If you want additional precision (at least 24 digits)
for the intermediate computation, the required hardware
is 30 digits (i.e., double precision).
> I'm sorry, but I don't see any solution to this. I you
> have one, I'm very interested in...
Clearly, the request is satisfied by:
REAL(selected_real_kind(precision(1.0)+1) :: accu
Thus, on a 32-bit machine, accu has some 15 decimal digits,
allowing accu to be computed with the desired precision
of at least 12 decimal digits.
And on a 64-bit machine, accu has some 30 digits,
allowing accu to be computed with the desired precision
of at least 24 decimal digits.
And of course, whether the calling procedure specifies
the argument corresponding to dummy argument array "a" as
REAL (SELECTED_REAL_KIND (6) ) :: X
or as
REAL :: X,
the procedure MYSUM forms the intermediate sum in
double precision, guaranteeing at least 12 digits on a 32-bit machine.
On a 64-bit machine, the intermediate sum is formed
in double precision -- at least 24 decimal digits -- which
is exactly what you want.
Of course, if you want only 6 digit precision on any machine,
with 12-digit intermediate precision on any machine,
then it's only necessary to write the procedure MYSUM
with the declarations:
REAL, INTENT (IN) :: A( : )
REAL (KIND=SELECTED_REAL_KIND(12)) :: accu
and regardless of which machine (32-bit or 64-bit)
you get a minimum of 6 decimal digits for the dummy argument "a"
and a minimum of 12 decimal digits for local array "accu".
What's more, both variations are portable, and can be written
as a generic procedure for handling single and double
precision arguments.
> Best Regards
> | Pierre Hugonnet
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|