* Bug #9196 fixed - The threshold level for conditioning in backslash was too

                    small.

Change-Id: Ie0056d871047b63092542e73c19da0a609cb963b

diff --git a/scilab/CHANGES_5.4.X b/scilab/CHANGES_5.4.X
index 22ee4ed..1cf1360 100644 (file)
--- 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 (file)
--- 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 (file)
--- 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 (file)
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 (file)
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);