From 06949c1990a8f85aaaf7d6785019c092eb653090 Mon Sep 17 00:00:00 2001 From: =?utf8?q?Micha=C3=ABl=20Baudin?= Date: Wed, 6 Mar 2013 09:10:21 +0100 Subject: [PATCH] * Bug #9196 fixed - The threshold level for conditioning in backslash was too small. Change-Id: Ie0056d871047b63092542e73c19da0a609cb963b --- scilab/CHANGES_5.4.X | 3 + .../modules/linear_algebra/src/fortran/intdgesv3.f | 10 ++- .../modules/linear_algebra/src/fortran/intzgesv3.f | 7 +- .../tests/nonreg_tests/bug_9196.dia.ref | 69 ++++++++++++++++++ .../linear_algebra/tests/nonreg_tests/bug_9196.tst | 76 ++++++++++++++++++++ 5 files changed, 160 insertions(+), 5 deletions(-) create mode 100644 scilab/modules/linear_algebra/tests/nonreg_tests/bug_9196.dia.ref create mode 100644 scilab/modules/linear_algebra/tests/nonreg_tests/bug_9196.tst diff --git a/scilab/CHANGES_5.4.X b/scilab/CHANGES_5.4.X index 22ee4ed..1cf1360c 100644 --- a/scilab/CHANGES_5.4.X +++ b/scilab/CHANGES_5.4.X @@ -254,6 +254,9 @@ Bug fixes * Bug #9005 fixed - The bitset function did not have any tests. +* Bug #9196 fixed - The threshold level for conditioning in backslash was too + small. + * Bug #9305 fixed - In the optimisation help page, new chapter "Least Squares functions" created. diff --git a/scilab/modules/linear_algebra/src/fortran/intdgesv3.f b/scilab/modules/linear_algebra/src/fortran/intdgesv3.f index 669c72e..c71b1e6 100644 --- a/scilab/modules/linear_algebra/src/fortran/intdgesv3.f +++ b/scilab/modules/linear_algebra/src/fortran/intdgesv3.f @@ -1,5 +1,6 @@ c Scilab ( http://www.scilab.org/ ) - This file is part of Scilab c Copyright (C) INRIA +c Copyright (C) 2013 - Michael Baudin c$ c This file must be used under the terms of the CeCILL. c This source file is licensed as described in the file COPYING, which @@ -17,10 +18,12 @@ c a\b character fname*(*) double precision ANORM, EPS, RCOND double precision dlamch, dlange + double precision RCONDthresh integer vfinite external dlamch, dlange, vfinite intrinsic sqrt -c + +c minrhs=2 maxrhs=2 minlhs=1 @@ -74,6 +77,7 @@ c Check if A and B matrices contains Inf or NaN's if(.not.createvar(10,'d',1,LWORK,lDWORK)) return EPS = dlamch('eps') + RCOND_thresh=EPS*10 ANORM = dlange( '1', M, N, stk(lA), M, stk(lDWORK) ) c if(M.eq.N) then @@ -90,7 +94,7 @@ c SUBROUTINE DGETRF( N, N, A, LDA, IPIV, INFO ) $ istk(lIWORK), INFO ) c SUBROUTINE DGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, c $ IWORK, INFO ) - if(RCOND.gt.sqrt(EPS)) then + if(RCOND.gt.RCONDthresh) then call DGETRS( 'N', N, NRHS, stk(lAF), N, istk(lIPIV), $ stk(lB), N, INFO ) c SUBROUTINE DGETRS( TRANS, N, NRHS, A, LDA, IPIV, @@ -107,7 +111,7 @@ c SUBROUTINE DLACPY( UPLO, M, N, A, LDA, B, LDB ) c c M.ne.N or A singular c - RCOND = sqrt(EPS) + RCOND = RCONDthresh call DLACPY( 'F', M, NRHS, stk(lB), M, stk(lXB), max(M,N) ) c SUBROUTINE DLACPY( UPLO, M, N, A, LDA, B, LDB ) do 10 i = 1, N diff --git a/scilab/modules/linear_algebra/src/fortran/intzgesv3.f b/scilab/modules/linear_algebra/src/fortran/intzgesv3.f index 6ceb4ce..cb83a47 100644 --- a/scilab/modules/linear_algebra/src/fortran/intzgesv3.f +++ b/scilab/modules/linear_algebra/src/fortran/intzgesv3.f @@ -18,6 +18,7 @@ c a\b character fname*(*) double precision ANORM, EPS, RCOND double precision dlamch, zlange + double precision RCOND_thresh integer vfinite external dlamch, zlange, vfinite intrinsic sqrt @@ -26,6 +27,7 @@ c maxrhs=2 minlhs=1 maxlhs=1 + c if(.not.checkrhs(fname,minrhs,maxrhs)) return if(.not.checklhs(fname,minlhs,maxlhs)) return @@ -75,6 +77,7 @@ c Check if A and B matrices contains Inf or NaN's if(.not.createvar(10,'z',1,LWORK,lDWORK)) return EPS = dlamch('eps') + RCOND_thresh=EPS*10 ANORM = zlange( '1', M, N, zstk(lA), M, zstk(lDWORK) ) c if(M.eq.N) then @@ -92,7 +95,7 @@ c SUBROUTINE ZGETRF( N, N, A, LDA, IPIV, INFO ) $ zstk(lDWORK), stk(lRWORK), INFO ) c SUBROUTINE ZGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, c $ RWORK, INFO ) - if(RCOND.gt.sqrt(EPS)) then + if(RCOND.gt.RCOND_thresh) then call ZGETRS( 'N', N, NRHS, zstk(lAF), N, istk(lIPIV), $ zstk(lXB), N, INFO ) c SUBROUTINE ZGETRS( TRANS, N, NRHS, A, LDA, IPIV, @@ -111,7 +114,7 @@ c . ill conditioned problem c c M.ne.N or A singular c - RCOND = sqrt(EPS) + RCOND = RCOND_thresh call ZLACPY( 'F', M, NRHS, zstk(lB), M, zstk(lXB), max(M,N) ) c SUBROUTINE ZLACPY( UPLO, M, N, A, LDA, B, LDB ) do 10 i = 1, N diff --git a/scilab/modules/linear_algebra/tests/nonreg_tests/bug_9196.dia.ref b/scilab/modules/linear_algebra/tests/nonreg_tests/bug_9196.dia.ref new file mode 100644 index 0000000..ab9ef99 --- /dev/null +++ b/scilab/modules/linear_algebra/tests/nonreg_tests/bug_9196.dia.ref @@ -0,0 +1,69 @@ +// ============================================================================= +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) 2013 - Scilab Enterprises - Sylvestre Ledru +// Copyright (C) 2013 - Michaël Baudin +// +// This file is distributed under the same license as the Scilab package. +// ============================================================================= +// <-- CLI SHELL MODE --> +// <-- Non-regression test for bug 9196 --> +// +// <-- Bugzilla URL --> +// http://bugzilla.scilab.org/show_bug.cgi?id=9196 +// +// <-- Short Description --> +// The threshold level for conditioning in backslash is too small. +// ============================================================================= +n=9; +A = testmatrix("hilb",n); +b=ones(n,1); +xexpected=[7129/2520 + 4861/2520 + 42131/27720 + 35201/27720 + 395243/360360 + 348911/360360 + 62575/72072 + 113567/144144 + 1768477/2450448]; +x=A\b; +assert_checkalmostequal(x, xexpected); +Ac=complex(A,zeros(A)); +bc=complex(b,zeros(b)); +xc = Ac\bc; +xcexpected=complex(xexpected,zeros(xexpected)); +assert_checkalmostequal(xc, xexpected); +b=(1:n)'; +xexpected=[9; + 17819/2520; + 82609/13860; + 47959/9240; + 415567/90090; + 299737/72072; + 45533/12012; + 71761/20592; + 988277/306306]; +x=A\b; +assert_checkalmostequal(x, xexpected); +Ac=complex(A,zeros(A)); +bc=complex(b,zeros(b)); +xc = Ac\bc; +xcexpected=complex(xexpected,zeros(xexpected)); +assert_checkalmostequal(xc, xexpected); +b=[1;-1;1;-1;1;-1;1;-1;1]; +xexpected=[1879/2520; + 893/2520; + 6557/27720; + 4993/27720; + 52901/360360; + 44911/360360; + 39173/360360; + 69659/720720; + 1068047/12252240]; +x=A\b; +assert_checkalmostequal(x, xexpected); +Ac=complex(A,zeros(A)); +bc=complex(b,zeros(b)); +xc = Ac\bc; +xcexpected=complex(xexpected,zeros(xexpected)); +assert_checkalmostequal(xc, xexpected); diff --git a/scilab/modules/linear_algebra/tests/nonreg_tests/bug_9196.tst b/scilab/modules/linear_algebra/tests/nonreg_tests/bug_9196.tst new file mode 100644 index 0000000..7c91813 --- /dev/null +++ b/scilab/modules/linear_algebra/tests/nonreg_tests/bug_9196.tst @@ -0,0 +1,76 @@ +// ============================================================================= +// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab +// Copyright (C) 2013 - Scilab Enterprises - Sylvestre Ledru +// Copyright (C) 2013 - Michaël Baudin +// +// This file is distributed under the same license as the Scilab package. +// ============================================================================= + +// <-- CLI SHELL MODE --> + +// <-- Non-regression test for bug 9196 --> +// +// <-- Bugzilla URL --> +// http://bugzilla.scilab.org/show_bug.cgi?id=9196 +// +// <-- Short Description --> +// The threshold level for conditioning in backslash is too small. +// ============================================================================= + +n=9; +A = testmatrix("hilb",n); +b=ones(n,1); +xexpected=[7129/2520 + 4861/2520 + 42131/27720 + 35201/27720 + 395243/360360 + 348911/360360 + 62575/72072 + 113567/144144 + 1768477/2450448]; +x=A\b; +assert_checkalmostequal(x, xexpected); + + +Ac=complex(A,zeros(A)); +bc=complex(b,zeros(b)); +xc = Ac\bc; +xcexpected=complex(xexpected,zeros(xexpected)); +assert_checkalmostequal(xc, xexpected); + +b=(1:n)'; +xexpected=[9; + 17819/2520; + 82609/13860; + 47959/9240; + 415567/90090; + 299737/72072; + 45533/12012; + 71761/20592; + 988277/306306]; +x=A\b; +assert_checkalmostequal(x, xexpected); +Ac=complex(A,zeros(A)); +bc=complex(b,zeros(b)); +xc = Ac\bc; +xcexpected=complex(xexpected,zeros(xexpected)); +assert_checkalmostequal(xc, xexpected); + +b=[1;-1;1;-1;1;-1;1;-1;1]; +xexpected=[1879/2520; + 893/2520; + 6557/27720; + 4993/27720; + 52901/360360; + 44911/360360; + 39173/360360; + 69659/720720; + 1068047/12252240]; +x=A\b; +assert_checkalmostequal(x, xexpected); +Ac=complex(A,zeros(A)); +bc=complex(b,zeros(b)); +xc = Ac\bc; +xcexpected=complex(xexpected,zeros(xexpected)); +assert_checkalmostequal(xc, xexpected); -- 1.7.9.5