From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 12798 invoked by alias); 22 Aug 2010 21:02:59 -0000 Received: (qmail 10104 invoked by uid 48); 22 Aug 2010 21:02:45 -0000 Date: Sun, 22 Aug 2010 21:02:00 -0000 Message-ID: <20100822210245.10101.qmail@sourceware.org> X-Bugzilla-Reason: CC References: Subject: [Bug fortran/33197] Fortran 2008: math functions In-Reply-To: Reply-To: gcc-bugzilla@gcc.gnu.org To: gcc-bugs@gcc.gnu.org From: "burnus at gcc dot gnu dot org" Mailing-List: contact gcc-bugs-help@gcc.gnu.org; run by ezmlm Precedence: bulk List-Id: List-Archive: List-Post: List-Help: Sender: gcc-bugs-owner@gcc.gnu.org X-SW-Source: 2010-08/txt/msg01831.txt.bz2 ------- 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