From patchwork Sun Aug 22 09:55:40 2010 Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit X-Patchwork-Submitter: Tobias Burnus X-Patchwork-Id: 62368 Return-Path: X-Original-To: incoming@patchwork.ozlabs.org Delivered-To: patchwork-incoming@bilbo.ozlabs.org Received: from sourceware.org (server1.sourceware.org [209.132.180.131]) by ozlabs.org (Postfix) with SMTP id 915F4B70A3 for ; Sun, 22 Aug 2010 19:55:58 +1000 (EST) Received: (qmail 22816 invoked by alias); 22 Aug 2010 09:55:54 -0000 Received: (qmail 22794 invoked by uid 22791); 22 Aug 2010 09:55:51 -0000 X-SWARE-Spam-Status: No, hits=-1.9 required=5.0 tests=AWL, BAYES_00, RCVD_IN_DNSWL_NONE X-Spam-Check-By: sourceware.org Received: from mx02.qsc.de (HELO mx02.qsc.de) (213.148.130.14) by sourceware.org (qpsmtpd/0.43rc1) with ESMTP; Sun, 22 Aug 2010 09:55:44 +0000 Received: from [192.168.178.22] (port-92-204-8-164.dynamic.qsc.de [92.204.8.164]) by mx02.qsc.de (Postfix) with ESMTP id 849741DAA4; Sun, 22 Aug 2010 11:55:41 +0200 (CEST) Message-ID: <4C70F41C.8050801@net-b.de> Date: Sun, 22 Aug 2010 11:55:40 +0200 From: Tobias Burnus User-Agent: Mozilla/5.0 (X11; U; Linux x86_64; de; rv:1.9.2.7) Gecko/20100714 SUSE/3.1.1 Thunderbird/3.1.1 MIME-Version: 1.0 To: gcc patches , gfortran Subject: [Patch, Fortran, committed] PR 45367/36158 - Increase numeric tolerance in bessel_{6, 7}.f90 Mailing-List: contact gcc-patches-help@gcc.gnu.org; run by ezmlm Precedence: bulk List-Id: List-Unsubscribe: List-Archive: List-Post: List-Help: Sender: gcc-patches-owner@gcc.gnu.org Delivered-To: mailing list gcc-patches@gcc.gnu.org As HJ and Dominique reported, for Linux/ia32 and powerpc-apple-darwin9, the tolerance values were too tight for the comparison of the recurrence algorithm with the library version. I increased now the tolerance based on Dominique's values and cross checked using -m32 on x86-64-linux. I hope it will now work on most systems, but we still might need to tweak the values (or use xfail) for some architectures as the initial values are slightly different and as the architectures might do different rounding or take different code paths. Thanks to Dominique for testing and adjusting the values - and to HJ for keeping an eye on regressions. Committed as Rev. 163454/163455. * * * Note of warning to users: Be careful when applying the recurrence algorithm (i.e. when using BESSEL_JN/YN(N1, N2, X): The results can be quite different! For instance, using BESSEL_JN(0, 100, 475.78) the deviation is < 6*epsilon(0.0) - except for the following N, where it is (much) larger: N=63: < 9*epsilon(0.0) N=91: < 13*epsilon(0.0) N=49: < 326*epsilon(0.0) Note in particular that checking the last item in the recurrence algorithm - for JN that's lowest, here: N=0, does not help: For the example above the difference to the value returned by the library for JN is just 0.5*epsilon while for N=49 the difference is much larger! (The reason is that for N=49 the result is -0.8E-04 while the other values are of the order 0.xE-01 and 0.xE-02.) . Tobias Index: gcc/testsuite/ChangeLog =================================================================== --- gcc/testsuite/ChangeLog (Revision 163454) +++ gcc/testsuite/ChangeLog (Arbeitskopie) @@ -1,4 +1,9 @@ 2010-08-22 Tobias Burnus + + PR fortran/36158 + * gfortran.dg/bessel_7.f90: Disable accidently enabled debug output. + +2010-08-22 Tobias Burnus Dominique d'Humieres PR fortran/45367 Index: gcc/testsuite/gfortran.dg/bessel_7.f90 =================================================================== --- gcc/testsuite/gfortran.dg/bessel_7.f90 (Revision 163454) +++ gcc/testsuite/gfortran.dg/bessel_7.f90 (Arbeitskopie) @@ -32,13 +32,13 @@ rec = BESSEL_YN(0, Nmax, X) lib = [ (BESSEL_YN(i, X), i=0,Nmax) ] -print *, 'YN for X = ', X, ' -- Epsilon = ',epsilon(x) +!print *, 'YN for X = ', X, ' -- Epsilon = ',epsilon(x) do i = 0, Nmax - print '(i2,2e17.9,e12.2,f14.10,2l3)', i, rec(i), lib(i), & - rec(i)-lib(i), ((rec(i)-lib(i))/rec(i))/epsilon(x), & - i > nit .or. rec(i) == lib(i) & - .or. abs((rec(i)-lib(i))/rec(i)) < myeps2, & - rec(i) == lib(i) .or. abs((rec(i)-lib(i))/rec(i)) < myeps +! print '(i2,2e17.9,e12.2,f14.10,2l3)', i, rec(i), lib(i), & +! rec(i)-lib(i), ((rec(i)-lib(i))/rec(i))/epsilon(x), & +! i > nit .or. rec(i) == lib(i) & +! .or. abs((rec(i)-lib(i))/rec(i)) < myeps2, & +! rec(i) == lib(i) .or. abs((rec(i)-lib(i))/rec(i)) < myeps if (.not. (i > nit .or. rec(i) == lib(i) & .or. abs((rec(i)-lib(i))/rec(i)) < myeps2)) & call abort ()