------- Comment #39 from burnus at gcc dot gnu dot org 2010-08-22 21:02 ------- (In reply to comment #37) > PS: NORM2 is described as "careful calculation of Euclidean norm" in the BCS > slides and in the what's new in F2008 article. Currently, I use the trivial > brute-force method. Maybe something more careful should be done?
As Dominique points out - the algorithm can be made more robust by doing the calculation as tmp = max(abs(a)) NORM2(a) := tmp*sqrt(dot_product(a/tmp,a/tmp)) That helps a lot for "a" finite with (a^2 > huge(a)) [overflow] ;-) However, there is a method which only requires a single pass, cf. p. 38/39 in http://cpc.cs.qub.ac.uk/MRSN/higham.pdf. The algorithm by Sven Hammarling is also used in BLAS, cf. http://www.netlib.org/blas/snrm2.f -- http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33197
