* Bug #10180 fixed - det was not defined for sparse matrices. 58/11658/5
Charlotte HECQUET [Thu, 6 Jun 2013 08:21:48 +0000 (10:21 +0200)]
Change-Id: I75fd51d884f3c54aca4751be87d5eb48dfafc454

scilab/CHANGES_5.5.X
scilab/modules/linear_algebra/help/en_US/matrix/det.xml
scilab/modules/linear_algebra/help/fr_FR/matrix/det.xml
scilab/modules/linear_algebra/tests/unit_tests/det.dia.ref
scilab/modules/linear_algebra/tests/unit_tests/det.tst
scilab/modules/overloading/macros/%sp_det.sci [new file with mode: 0644]
scilab/modules/overloading/tests/nonreg_tests/bug_10180.dia.ref [new file with mode: 0644]
scilab/modules/overloading/tests/nonreg_tests/bug_10180.tst [new file with mode: 0644]

index 970f035..e2b6690 100644 (file)
@@ -275,6 +275,8 @@ Bug fixes
 
 * Bug #10146 fixed - In SciNotes, 'help on keyword' was at the end of the popup menu.
 
+* Bug #10180 fixed - det was not defined for sparse matrices.
+
 * Bug #10213 fixed - sci2exp can be impacted by format which has been mentionned in its help page.
 
 * Bug #10226 fixed - When a // <empty session> line is delete, all sessions
index 0edd884..97053f0 100644 (file)
@@ -27,7 +27,7 @@
             <varlistentry>
                 <term>X</term>
                 <listitem>
-                    <para>real or complex square matrix, polynomial or rational matrix.</para>
+                    <para>real or complex square matrix (full or sparse), polynomial or rational matrix</para>
                 </listitem>
             </varlistentry>
             <varlistentry>
@@ -47,7 +47,7 @@
     <refsection>
         <title>Description</title>
         <para>
-            <literal>det(X)</literal> ( <literal>m*10^e</literal> is the determinant of the square matrix <literal>X</literal>.
+            <literal>det(X)</literal> ( <literal>m*10^e</literal> ) is the determinant of the square matrix <literal>X</literal>.
         </para>
         <para>
             For polynomial matrix <literal>det(X)</literal> is equivalent to <literal>determ(X)</literal>.
@@ -62,6 +62,9 @@
             det computations are based on the Lapack routines
             DGETRF for  real matrices and  ZGETRF for the complex case.
         </para>
+        <para>
+            Concerning sparse matrices, the determinant is obtained from LU factorisation of umfpack library.
+        </para>
     </refsection>
     <refsection>
         <title>Examples</title>
index d99e66e..7b26e0b 100644 (file)
@@ -27,7 +27,7 @@
             <varlistentry>
                 <term>X  </term>
                 <listitem>
-                    <para>matrice réelle, complexe, polynomiale, rationnelle
+                    <para>matrice carrée réelle ou complexe (creuse ou pleine), polynomiale ou rationnelle
                     </para>
                 </listitem>
             </varlistentry>
@@ -50,7 +50,7 @@
     <refsection>
         <title>Description</title>
         <para>
-            <literal>det(X)</literal> ( <literal>m*10^e</literal> est le déterminant de la matrice carrée <literal>X</literal>.
+            <literal>det(X)</literal> ( <literal>m*10^e</literal> ) est le déterminant de la matrice carrée <literal>X</literal>.
         </para>
         <para>
             Pour les matrices polynomiales <literal>det(X)</literal> est équivalent à <literal>determ(X)</literal>.
@@ -86,5 +86,9 @@ det(A), prod(spec(A))
             Le calcul du determinant est basé sur les routines Lapack :
             DGETRF pour les matrices réelles et  ZGETRF pour le cas complexe.
         </para>
+        <para>
+            Concernant le cas des matrices creuses, le calcul du déterminant est effectué
+            à partir de la décomposition LU de la librairie umfpack.
+        </para>
     </refsection>
 </refentry>
index c98bd71..3a218b7 100644 (file)
 //
 //  This file is distributed under the same license as the Scilab package.
 // =============================================================================
-function r=Err(x),r=norm(x,1),endfunction
-rand('normal')
 //==========================================================================
-//==============================   det        ==============================
+//==============================   det        ============================== 
 //==========================================================================
-if execstr('det([1 2;3 4;5 6])','errcatch')==0 then bugmes();quit;end
-//Small dimension
+//Small dimension 
 //Real
 A=[1 1; 1 2];
-if Err(det(A)-1)>10*%eps then bugmes();quit;end
+assert_checkalmostequal(det(A), 1);
 [e,m]=det(A);
-if e<>0 |Err(m-1)>10*%eps then bugmes();quit;end
+assert_checkalmostequal(e, 0);
+assert_checkalmostequal(m, 1);
 //Complex
 A=A+%i;
-if Err(det(A)-1-%i)>10*%eps then bugmes();quit;end
+assert_checkalmostequal(real(det(A)), 1);
+assert_checkalmostequal(imag(det(A)), 1);
 [e,m]=det(A);
-if e<>0 |Err(m-1-%i)>10*%eps then bugmes();quit;end
+assert_checkalmostequal(e, 0);
+assert_checkalmostequal(real(m), 1);
+assert_checkalmostequal(imag(m), 1);
+//Sparse
+A=[1 1; 1 2];
+A=sparse(A);
+assert_checkalmostequal(det(A), 1);
+[e,m]=det(A)
+ m  =
+    1.  
+ e  =
+    0.  
+assert_checkalmostequal(e, 0);
+assert_checkalmostequal(m, 1);
+//Polynomials
+A=[1+%s 1; 1 2+%s];
+assert_checkequal(det(A), 1+3*%s+%s*%s);
+//Rationals
+A=[1+%s 1/%s; 1 2+%s];
+assert_checkequal(numer(det(A)), -1+2*%s+3*%s^2+%s^3);
+assert_checkequal(denom(det(A)), %s);
+//Sparse complex
+A=[1 1; 1 2];
+A=A+%i;
+A=sparse(A);
+assert_checkalmostequal(real(det(A)), 1);
+assert_checkalmostequal(imag(det(A)), 1);
+[e,m]=det(A);
+assert_checkalmostequal(e, 0);
+assert_checkalmostequal(real(m), 1);
+assert_checkalmostequal(imag(m), 1);
 //Large dimension
 //Real
 v=rand(1,21);
 A=rand(21,21); A=(triu(A,1)+diag(v))*(tril(A,-1)+diag(ones(1,21)));
-if Err(det(A)-prod(v))>400000*%eps then bugmes();quit;end
+assert_checktrue(abs(det(A) - prod(v)) < 1D-7);
 [e,m]=det(A);
-if Err(m*(10^e)-prod(v))>400000*%eps then bugmes();quit;end
+assert_checktrue(abs(m*(10^e) - prod(v)) < 1D-7);
 //Complex
 v=(v+rand(v)*%i)/2;
 A=rand(21,21); A=(triu(A,1)+diag(v))*(tril(A,-1)+diag(ones(1,21)));
-if Err(det(A)-prod(v))>10000*%eps then bugmes();quit;end
+assert_checktrue(abs(det(A) - prod(v)) < 1D-7);
+[e,m]=det(A);
+assert_checktrue(abs(m*(10^e) - prod(v)) < 1D-7);
+//Sparse
+v=rand(1,21);
+v=sparse(v);
+A=rand(21,21); A=(triu(A,1)+diag(v))*(tril(A,-1)+diag(ones(1,21)));
+assert_checktrue(abs(det(A) - prod(v)) < 1D-7);
 [e,m]=det(A);
-if Err(m*(10^e)-prod(v))>10000*%eps then bugmes();quit;end
+assert_checktrue(abs(m*(10^e) - prod(v)) < 1d-7);
+//Polynomials
+v=rand(1,21)
+ v  =
+         column 1 to 5
+    0.5866105    0.0477090    0.9416931    0.9613204    0.5583350  
+         column  6 to 10
+    0.5700629    0.3169258    0.9932628    0.8074780    0.8554797  
+         column 11 to 15
+    0.5031461    0.0963323    0.7058098    0.8630577    0.0076185  
+         column 16 to 20
+    0.8048951    0.5963762    0.1176836    0.8010095    0.5132340  
+         column 21
+    0.201091  
+v=v+%s;
+A=rand(21,21); A=(triu(A,1)+diag(v))*(tril(A,-1)+diag(ones(1,21)));
+assert_checktrue(abs(coeff(det(A)-prod(v))) < 1D-7);
+//Rationals
+v=rand(1,21);
+v=v/%s;
+A=rand(21,21); A=(triu(A,1)+diag(v))*(tril(A,-1)+diag(ones(1,21)));
+assert_checktrue(abs(coeff(numer(det(A))-numer(prod(v)))) < 1D-7);
+//Sparse complex
+v=rand(1,21);
+v=(v+rand(v)*%i)/2;
+v=sparse(v);
+A=rand(21,21);
+A=(triu(A,1)+diag(v))*(tril(A,-1)+diag(ones(1,21)));
+A=sparse(A);
+assert_checktrue(abs(det(A) - prod(v)) < 1D-7);
+[e,m]=det(A);
+assert_checktrue(abs(m*(10^e) - prod(v)) < 1d-7);
+//Error messages
+A=[1 1; 1 2];
+errmsg1 = msprintf(_("Wrong type for first argument: Square matrix expected.\n"));
+assert_checkerror("det([1,2;3,4;5,6])", errmsg1, 20);
+errmsg2 = msprintf(_("%s: Wrong number of input argument(s): %d expected.\n"), "det", 1);
+assert_checkerror("det(A,1)", errmsg2, 77);
index c8f9c39..4b02407 100644 (file)
@@ -4,35 +4,95 @@
 //
 //  This file is distributed under the same license as the Scilab package.
 // =============================================================================
-function r=Err(x),r=norm(x,1),endfunction
-rand('normal')
 
 //==========================================================================
 //==============================   det        ============================== 
 //==========================================================================
-if execstr('det([1 2;3 4;5 6])','errcatch')==0 then pause,end
+
 //Small dimension 
 //Real
 A=[1 1; 1 2];
-if Err(det(A)-1)>10*%eps then pause,end
+assert_checkalmostequal(det(A), 1);
 [e,m]=det(A);
-if e<>0 |Err(m-1)>10*%eps then pause,end
+assert_checkalmostequal(e, 0);
+assert_checkalmostequal(m, 1);
 //Complex
 A=A+%i;
-if Err(det(A)-1-%i)>10*%eps then pause,end
+assert_checkalmostequal(real(det(A)), 1);
+assert_checkalmostequal(imag(det(A)), 1);
+[e,m]=det(A);
+assert_checkalmostequal(e, 0);
+assert_checkalmostequal(real(m), 1);
+assert_checkalmostequal(imag(m), 1);
+//Sparse
+A=[1 1; 1 2];
+A=sparse(A);
+assert_checkalmostequal(det(A), 1);
+[e,m]=det(A)
+assert_checkalmostequal(e, 0);
+assert_checkalmostequal(m, 1);
+//Polynomials
+A=[1+%s 1; 1 2+%s];
+assert_checkequal(det(A), 1+3*%s+%s*%s);
+//Rationals
+A=[1+%s 1/%s; 1 2+%s];
+assert_checkequal(numer(det(A)), -1+2*%s+3*%s^2+%s^3);
+assert_checkequal(denom(det(A)), %s);
+//Sparse complex
+A=[1 1; 1 2];
+A=A+%i;
+A=sparse(A);
+assert_checkalmostequal(real(det(A)), 1);
+assert_checkalmostequal(imag(det(A)), 1);
 [e,m]=det(A);
-if e<>0 |Err(m-1-%i)>10*%eps then pause,end
+assert_checkalmostequal(e, 0);
+assert_checkalmostequal(real(m), 1);
+assert_checkalmostequal(imag(m), 1);
+
 //Large dimension
 //Real
 v=rand(1,21);
 A=rand(21,21); A=(triu(A,1)+diag(v))*(tril(A,-1)+diag(ones(1,21)));
-if Err(det(A)-prod(v))>400000*%eps then pause,end
+assert_checktrue(abs(det(A) - prod(v)) < 1D-7);
 [e,m]=det(A);
-if Err(m*(10^e)-prod(v))>400000*%eps then pause,end
+assert_checktrue(abs(m*(10^e) - prod(v)) < 1D-7);
 //Complex
 v=(v+rand(v)*%i)/2;
 A=rand(21,21); A=(triu(A,1)+diag(v))*(tril(A,-1)+diag(ones(1,21)));
-if Err(det(A)-prod(v))>10000*%eps then pause,end
+assert_checktrue(abs(det(A) - prod(v)) < 1D-7);
+[e,m]=det(A);
+assert_checktrue(abs(m*(10^e) - prod(v)) < 1D-7);
+//Sparse
+v=rand(1,21);
+v=sparse(v);
+A=rand(21,21); A=(triu(A,1)+diag(v))*(tril(A,-1)+diag(ones(1,21)));
+assert_checktrue(abs(det(A) - prod(v)) < 1D-7);
+[e,m]=det(A);
+assert_checktrue(abs(m*(10^e) - prod(v)) < 1d-7);
+//Polynomials
+v=rand(1,21)
+v=v+%s;
+A=rand(21,21); A=(triu(A,1)+diag(v))*(tril(A,-1)+diag(ones(1,21)));
+assert_checktrue(abs(coeff(det(A)-prod(v))) < 1D-7);
+//Rationals
+v=rand(1,21);
+v=v/%s;
+A=rand(21,21); A=(triu(A,1)+diag(v))*(tril(A,-1)+diag(ones(1,21)));
+assert_checktrue(abs(coeff(numer(det(A))-numer(prod(v)))) < 1D-7);
+//Sparse complex
+v=rand(1,21);
+v=(v+rand(v)*%i)/2;
+v=sparse(v);
+A=rand(21,21);
+A=(triu(A,1)+diag(v))*(tril(A,-1)+diag(ones(1,21)));
+A=sparse(A);
+assert_checktrue(abs(det(A) - prod(v)) < 1D-7);
 [e,m]=det(A);
-if Err(m*(10^e)-prod(v))>10000*%eps then pause,end
+assert_checktrue(abs(m*(10^e) - prod(v)) < 1d-7);
 
+//Error messages
+A=[1 1; 1 2];
+errmsg1 = msprintf(_("Wrong type for first argument: Square matrix expected.\n"));
+assert_checkerror("det([1,2;3,4;5,6])", errmsg1, 20);
+errmsg2 = msprintf(_("%s: Wrong number of input argument(s): %d expected.\n"), "det", 1);
+assert_checkerror("det(A,1)", errmsg2, 77);
diff --git a/scilab/modules/overloading/macros/%sp_det.sci b/scilab/modules/overloading/macros/%sp_det.sci
new file mode 100644 (file)
index 0000000..372334f
--- /dev/null
@@ -0,0 +1,27 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - Scilab Enterprises - Charlotte HECQUET
+// 
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at    
+// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
+
+function [res1, res2]=%sp_det(A)
+    [lhs, rhs]=argn(0);
+    hand = umf_lufact(A); //umfpack is used for complex sparse matrix
+    [L,U,P,Q,r] = umf_luget(hand);
+    res1=prod(r)*prod(diag(U));
+    res2=res1;
+    if (lhs == 2) then
+        res1=0;
+        while abs(res2) >= 10
+            if abs(res2) < 1 then
+                break;
+            end
+            res2 = res2 / 10;
+            res1 = res1 + 1;
+        end
+    end
+    umf_ludel(hand);
+endfunction
diff --git a/scilab/modules/overloading/tests/nonreg_tests/bug_10180.dia.ref b/scilab/modules/overloading/tests/nonreg_tests/bug_10180.dia.ref
new file mode 100644 (file)
index 0000000..f84ca35
--- /dev/null
@@ -0,0 +1,19 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - Scilab Enterprises - Charlotte HECQUET
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+// <-- Non-regression test for bug 10180 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=10180
+//
+// <-- Short Description -->
+// det is not defined for sparse matrix
+A=[1,2,3;3,4,5;0,0,5];
+sp=sparse(A);
+assert_checkequal(det(sp), det(A));
+[e,m]=det(sp);
+assert_checkequal(e, 1);
+assert_checkequal(m, -1);
diff --git a/scilab/modules/overloading/tests/nonreg_tests/bug_10180.tst b/scilab/modules/overloading/tests/nonreg_tests/bug_10180.tst
new file mode 100644 (file)
index 0000000..d0b7a5f
--- /dev/null
@@ -0,0 +1,21 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - Scilab Enterprises - Charlotte HECQUET
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+
+// <-- Non-regression test for bug 10180 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=10180
+//
+// <-- Short Description -->
+// det is not defined for sparse matrix
+
+A=[1,2,3;3,4,5;0,0,5];
+sp=sparse(A);
+assert_checkequal(det(sp), det(A));
+[e,m]=det(sp);
+assert_checkequal(e, 1);
+assert_checkequal(m, -1);