Underlying algorithm for HYPOT and NORM2 functions

[avinashs]

Is there a reference for the underlying algorithms for computing the Euclidian norms in HYPOT and NORM2 functions in IVF per the Fortran 2008 standard? I am particularly interested in the way overflow and underflow are handled. Also, whether they are the same as the PYTHAG and DNRM2 routines previously distributed as part of BLAS and SLATEC. The function reference in the IVF documentation has no details on the actual implementation.

[jimdempseyatthecove]

As to the internal algorithm I cannot say, I do suspect the results will vary dependent upon the /Qimf... option(s) selected. In particular the precision of the square root.

Jim Dempsey

[Steve_Lionel]

The Fortran standard does not go into that level of detail (and pretty much washes its hands of any out-of-range behavior.) Intel does not document algorithms or implementation of its math library (few vendors do.) It's very unlikely to match exactly some external library.

[avinashs]

Thanks for the responses. Having replaced a legacy code that uses a BLAS function to calculate sqrt(x*x + y*y) without overflow or destructive underflow with the HYPOT function in Fortran 2008, I am experiencing a loss in numerical performance on problems that have ill-conditioned matrices when using the HYPOT function. With regard to the response by @jimdempseyatthecove ... it is interesting to note that the legacy BLAS function implements the HYPOT function without actually calling a SQRT function. The internal HYPOT function results do not exactly match the result of the BLAS function, which suggests they are not the same and hence not an exact replacement. Hence, I will use these internal functions with due caution for production codes.

It may also be noted that the IVF documentation for HYPOT states that "the result has a value equal to a processor-dependent approximation to the Euclidean distance, without undue overflow or underflow."

[jimdempseyatthecove]

Note, if the internal code of HYPOT, as well as your version, is using a programmed algorithm, and the results are less precise as with older compiler code, the cause is likely due to the newer compile using AVX... registers as opposed to using the FPU instruction set.

[Steve_Lionel]

The Intel math library does a lot of alternate code paths depending on the CPU type, no matter what the compiler does. This can also affect the results. I will put in another plug for Improving Numerical Reproducibility in C/C++/Fortran

[avinashs]

Thanks once again, @Steve Lionel (Ret.) and @jimdempseyatthecove. Steve's article is very interesting and could explain some anamolies that we are seeing in Fortran codes written for nonlinear optimization. This is an area of numerical analysis where iterative methods are employed. An improved solution XNEW is obtained from a current solution XC by adding a correction vector DX. As convergence is approached, DX must become very small (of the order of 10^(-8) in double precision) before the iteration is terminated. This is where we begin to see the anamolies. Using HYPOT, a particular case took 182 iterations to declare convergence whereas using a library PYTHAG routine took 61 iterations.

[jimdempseyatthecove]

Note, it is not the magnitude of dX alone that matters, but rather the magnitude differences between X and dX. Try to make use of the epsilon of X in determining convergence. E.g. within n bits of the epsilon of X...

Jim Dempsey