Pierre Hugonnet wrote:
> Dick Hendrickson wrote:
> >
> > >This goes totally against SELECTED_REAL_KIND spirit, whose
> > >intent is to request a given precision.
> > >
> > Not necessarily. Sometimes you need two different precisions. To
> accumulate
> > sums in high precision to reduce roundoff problems; or to check an
> > algorithm for accuracy.
>
> INTERFACE mysum
> MODULE PROCEDURE mysum_single, mysum_double
> END INTEFACE
>
> 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
>
> SUBROUTINE mysum_double(a)
> DOUBLE PRECISION :: a(:)
> REAL(selected_real_kind(precision(1.0d0)+1)) :: accu
This specifies quadruple precision, and would
not be portable.
> accu = 0.0
> DO i=1,size(a)
> accu = accu + a(i)
> END DO
> mysum_double = accu
>
> END SUBROUTINE
>
> Imagine that you are on a 64bit machine, where REAL is 64bits,
> DOUBLE PRECISION is 128bits "quad precision".
>
> in the main, I have two arrays:
>
> REAL(selected_real_kind(p=6)) :: a6(n)
> REAL(selected_real_kind(p=12)) :: a12(n)
>
> Both are in fact default REAL: mysum_single() will be called for both,
> but I actually don't need a DOUBLE PRECISION accu for a6, since
> a REAL accu has already 12 digits precision. The problem
> is that in mysum I don't know what was the requested precision
> of the actual argument.
>
> The solution is to write the two routines using explicit
> requested precisions:
>
> SUBROUTINE mysum_6(a)
> REAL(selected_real_kind(p=6)) :: a(:)
> REAL(selected_real_kind(p=12)) :: accu
>
> accu = 0.0
> DO i=1,size(a)
> accu = accu + a(i)
> END DO
> mysum = accu
>
> END SUBROUTINE
>
> SUBROUTINE mysum_12(a)
> REAL(selected_real_kind(p=12)) :: a(:)
> REAL(selected_real_kind(p=24)) :: accu
Again, this may get into trouble, as it's excessive
precision.
> accu = 0.0
> DO i=1,size(a)
> accu = accu + a(i)
> END DO
> mysum = accu
>
> END SUBROUTINE
>
> But now I can't define a generic name, since the kinds are
> identical.
We showed you how to resolve this before.
> Best regards
> | Pierre Hugonnet | mail....CGG |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|