* Bug 16636 fixed: det(sparse) now actually implemented 91/21691/6
Samuel GOUGEON [Wed, 27 Jan 2021 22:05:08 +0000 (23:05 +0100)]
  http://bugzilla.scilab.org/16636

  det() page overhauled: http://bugzilla.scilab.org/attachment.cgi?id=5223

Change-Id: I4da698b36746b1989a788c683147e407144373d6

scilab/CHANGES.md
scilab/modules/linear_algebra/help/en_US/matrix/det.xml
scilab/modules/linear_algebra/help/fr_FR/matrix/det.xml
scilab/modules/linear_algebra/help/ja_JP/matrix/det.xml [deleted file]
scilab/modules/linear_algebra/help/pt_BR/matrix/det.xml [deleted file]
scilab/modules/linear_algebra/help/ru_RU/matrix/det.xml [new file with mode: 0644]
scilab/modules/linear_algebra/tests/unit_tests/det.tst
scilab/modules/overloading/macros/%sp_det.sci

index f654680..08399fc 100644 (file)
@@ -201,6 +201,7 @@ Feature changes and additions on 6.1.1
   - It can now sort any sparse 2D matrix, in all `g, r, c, lr, lc` methods, including sparse booleans and in multi-level mode. It was formerly limited to sparse real or complex vectors and only to the `g` mode.
   - Any hypermatrix can be sorted along a dimension > 2.
 * `unique` is enabled for any 2D sparse arrays, in simple, 'c' and 'r' modes.
+<<<<<<< HEAD
 * `%chars` constant added, to easily access to some selected sets of unicode symbols.
 * Lists are displayed in a more compact and comprehensive way.
 * `interp1` is upgraded:
@@ -226,6 +227,10 @@ Feature changes and additions on 6.1.1
 * `polyint` is introduced to compute polynomial antiderivatives.
 * Listbox uicontrol callback is now triggered by item click in single selection mode. For example, it allows successive execution of a demo in the demonstrations gui.
 
+=======
+* %chars constant added, to easily access to some selected sets of unicode symbols.
+* `det` is now actually extended to sparse matrices.
+>>>>>>> 608020c6bbd (* Bug 16636 fixed: det(sparse) now actually implemented)
 
 Help pages:
 -----------
@@ -428,6 +433,7 @@ Bug Fixes
 * [#16629](https://bugzilla.scilab.org/16629): `interp1`'s documentation did not tell the spline edges conditions ; extrapolation modes were poorly explained. ; the description of the result's size was completely wrong ; x as an option was not documented. A wrong extrapolation value could silently return a wrong result. There was some dead code like `if varargin(5)==%nan`. A bugged error message yielded its own error. When x is implicit, the argument index in error messages could be wrong. `periodic` and `edgevalue` extrapolation modes were not available. `linear` extrapolation was not available for splines. When `xp` is an hypermatrix with `size(xp,1)==1`, the size of the result was irregular/wrong.
 * [#16631](https://bugzilla.scilab.org/16631): read-only handle properties were reported as unknown when trying to set them.
 * [#16632](https://bugzilla.scilab.org/16632): Scilab did not start with unsupported locale on macOS.
+* [#16636](https://bugzilla.scilab.org/16636): `det(sparse)` most often yielded `Nan`. `[e,m]=det(sparse)` was not actually implemented.
 * [#16644](https://bugzilla.scilab.org/16644): `input("message:")` yielded a wrong error message about `mprintf` in case of non-interpretable input.
 * [#16654](https://bugzilla.scilab.org/16654): `interp` was leaking memory.
 
index a032d96..c6df862 100644 (file)
@@ -2,8 +2,8 @@
 <!--
  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
  * Copyright (C) 2008 - INRIA
- *
  * Copyright (C) 2012 - 2016 - Scilab Enterprises
+ * Copyright (C) 2021 - Samuel GOUGEON - Le Mans Université
  *
  * This file is hereby licensed under the terms of the GNU GPL v2.0,
  * pursuant to article 5.3.4 of the CeCILL v.2.1.
  * along with this program.
  *
  -->
-<refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:svg="http://www.w3.org/2000/svg" xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org" xml:lang="en" xml:id="det">
+<refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink"
+          xmlns:svg="http://www.w3.org/2000/svg" xmlns:mml="http://www.w3.org/1998/Math/MathML"
+          xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org"
+          xml:lang="en" xml:id="det">
     <refnamediv>
         <refname>det</refname>
-        <refpurpose>determinant</refpurpose>
+        <refpurpose>determinant of a square matrix</refpurpose>
     </refnamediv>
     <refsynopsisdiv>
         <title>Syntax</title>
-        <synopsis>det(X)
-            [e,m]=det(X)
+        <synopsis>
+            d = det(X)
+            [e,m] = det(X)
         </synopsis>
     </refsynopsisdiv>
     <refsection>
             <varlistentry>
                 <term>X</term>
                 <listitem>
-                    <para>real or complex square matrix, polynomial or rational matrix.</para>
+                    square matrix of real or complex numbers, polynomials, or rationals.
+                    Sparse-encoded matrices accepted.
+                    <para/>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>d</term>
+                <listitem>
+                    Scalar of the <varname>X</varname>'s type: the determinant of
+                    <varname>X</varname>. If <varname>X</varname> is sparse-encoded,
+                    <varname>d</varname> is dense.
+                    <para/>
                 </listitem>
             </varlistentry>
             <varlistentry>
                 <term>m</term>
                 <listitem>
-                    <para>real or complex number, the determinant base 10 mantissae</para>
+                    real or complex number: the determinant base 10 mantissa, with
+                    <literal>abs(m) ∈ [1,10)</literal>. Not supported for <varname>X</varname>
+                    polynomial or rational.
+                    <para/>
                 </listitem>
             </varlistentry>
             <varlistentry>
                 <term>e</term>
                 <listitem>
-                    <para>integer, the determinant base 10 exponent</para>
+                    integer: the determinant base 10 exponent, such that
+                    <literal>d = m * 10<superscript>e</superscript></literal>.
+                    Not supported for <varname>X</varname> polynomial or rational.
+                    <para/>
                 </listitem>
             </varlistentry>
         </variablelist>
     <refsection>
         <title>Description</title>
         <para>
-            <literal>det(X)</literal> ( <literal>m*10^e</literal> is the determinant of the square matrix <literal>X</literal>.
+            <emphasis role="bold">d = det(X)</emphasis> yields the determinant of the matrix
+            <varname>X</varname>.
+        </para>
+        <para>
+            For a polynomial or rational matrix, <literal>d=det(X)</literal> uses <literal>determ(..)</literal>
+            whose algorithm is based on the FFT.
+            <literal>d=detr(X)</literal> can be alternatively used, based on the Leverrier algorithm.
+            Both methods yield equivalent results. For rational matrices, turning off <code>simp_mode(%f)</code>
+            might be required to get identical results.
         </para>
         <para>
-            For polynomial matrix <literal>det(X)</literal> is equivalent to <literal>determ(X)</literal>.
+            <emphasis role="bold">[e, m] = det(X)</emphasis> can be used only for a matrix of numbers.
+            This syntax allows to overcome computation's underflow or overflow, when <literal>abs(d)</literal>
+            is smaller than
+            <literal>number_properties("tiny")</literal> ≈ 2.23 10<superscript>-308</superscript> or
+            bigger than <literal>number_properties("huge")</literal> ≈ 1.80 10<superscript>308</superscript>.
         </para>
         <para>
-            For rational matrices <literal>det(X)</literal> is equivalent to <literal>detr(X)</literal>.
+            For denses matrices, <literal>det(..)</literal> is based on the Lapack routines
+            DGETRF for real matrices and  ZGETRF for the complex case.
         </para>
         <para>
-            <important>
-                The <literal>det</literal> and <literal>detr</literal> functions don't use the same algorithm.
-                For a rational fraction, <literal>det(X)</literal> is overloaded by <literal>%r_det(X)</literal> which is based on the <literal>determ</literal> function.
-                <literal>detr()</literal> uses the Leverrier method.
-            </important>
-            <warning>
-                Sometimes the <literal>det</literal> and <literal>detr</literal> functions may return different values for rational matrices.
-                In such cases you should set rational simplification mode off by using <code>simp_mode(%f)</code> to get the same result.
-            </warning>
+            For sparse matrices, the determinant is obtained from LU factorization thanks to the umfpack library.
         </para>
     </refsection>
     <refsection>
-        <title>References</title>
+        <title>Examples</title>
+        <programlisting role="example"><![CDATA[
+A = rand(3,3)*5;
+det(A)
+[e, m] = det(A)
+
+// Matrix of complex numbers:
+// A = grand(3,3,"uin",0,10) + grand(3,3,"uin",0,10)*%i
+A = [3+%i, 9+%i*3, 9+%i ; 8+%i*8, 4+%i*3, 7+%i*7 ; 4, 6+%i*2, 6+%i*9]
+det(A)
+[e, m] = det(A)
+abs(m)  // in [1, 10)
+     ]]></programlisting>
+        <screen><![CDATA[
+--> A = rand(3,3)*5;
+--> det(A)
+ ans  =
+  -10.805163
+
+--> [e, m] = det(A)
+ e  =
+   1.
+ m  =
+  -1.0805163
+
+--> // Matrix of complex numbers:
+--> A = [3+%i, 9+%i*3, 9+%i ; 8+%i*8, 4+%i*3, 7+%i*7 ; 4, 6+%i*2, 6+%i*9]
+ A  =
+   3. + i     9. + 3.i   9. + i
+   8. + 8.i   4. + 3.i   7. + 7.i
+   4. + 0.i   6. + 2.i   6. + 9.i
+
+--> det(A)
+ ans  =
+   745. - 225.i
+
+--> [e, m] = det(A)
+ e  =
+   2.
+ m  =
+   7.45 - 2.25i
+
+--> abs(m)  // in [1, 10)
+ ans  =
+   7.7823518
+]]></screen>
+        <para/>
         <para>
-            det computations are based on the Lapack routines
-            DGETRF for  real matrices and  ZGETRF for the complex case.
+            Very big or small determinants: underflow and overflow handling:
         </para>
+        <programlisting role="example"><![CDATA[
+// Very big determinant:
+n = 1000;
+A = rand(n, n);
+det(A)
+[e, m] = det(A)
+
+// Very small determinant (of a sparse-encoded matrix):
+A = (triu(sprand(n,n,1)) + diag(rand(1,n)))/1.5;
+det(A)
+prod(diag(A))
+[e, m] = det(A)
+A = A/2;
+det(A)
+[e, m] = det(A)
+     ]]></programlisting>
+        <screen><![CDATA[
+--> // Very big determinant:
+--> A = rand(n, n);
+--> det(A)
+ ans  =
+  -Inf
+
+--> [e, m] = det(A)  // -3.1199e743
+ e  =
+   743.
+ m  =
+  -3.1198687
+
+--> // Very small determinant (of a sparse-encoded matrix):
+--> n = 1000;
+--> A = (triu(sprand(n,n,1)) + diag(rand(1,n)))/1.5;
+--> det(A)
+ ans  =
+   5.21D-236
+
+--> prod(diag(A))
+ ans  =
+   5.21D-236
+
+--> [e, m] = det(A)
+ e  =
+  -236.
+ m  =
+   5.2119757
+
+--> A = A/2;
+--> det(A)
+ ans  =
+   0.
+
+--> [e, m] = det(A)
+ e  =
+  -537.
+ m  =
+   4.8641473
+]]></screen>
+        <para/>
         <para>
-            Concerning sparse matrices, the determinant is obtained from LU factorization of umfpack library.
+            Determinant of a polynomial matrix:
         </para>
-    </refsection>
-    <refsection>
-        <title>Examples</title>
         <programlisting role="example"><![CDATA[
-x=poly(0,'x');
-det([x,1+x;2-x,x^2])
-w=ssrand(2,2,4);roots(det(systmat(w))),trzeros(w)   //zeros of linear system
-A=rand(3,3);
-det(A), prod(spec(A))
- ]]></programlisting>
+s = %s;
+det([s, 1+s ; 2-s, s^2])
+
+w = ssrand(2,2,4);
+roots(det(systmat(w))),trzeros(w)   //zeros of linear system
+     ]]></programlisting>
+        <screen><![CDATA[
+--> det([s, 1+s ; 2-s, s^2])
+ ans  =
+  -2 -s +s² +s³
+
+--> w = ssrand(2,2,4);
+--> roots(det(systmat(w))),trzeros(w)
+ ans  =
+  -3.1907522 + 0.i
+   2.3596502 + 0.i
+
+ ans  =
+   2.3596502 + 0.i
+  -3.1907522 + 0.i
+]]></screen>
     </refsection>
     <refsection role="see also">
         <title>See also</title>
@@ -104,4 +245,15 @@ det(A), prod(spec(A))
             </member>
         </simplelist>
     </refsection>
+    <refsection role="history">
+        <title>History</title>
+        <revhistory>
+            <revision>
+                <revnumber>6.1.1</revnumber>
+                <revdescription>
+                    [e,m]=det(X) syntax extended to sparse matrices.
+                </revdescription>
+            </revision>
+        </revhistory>
+    </refsection>
 </refentry>
index c598296..8220cf2 100644 (file)
@@ -2,8 +2,8 @@
 <!--
  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
  * Copyright (C) 2008 - INRIA
- *
  * Copyright (C) 2012 - 2016 - Scilab Enterprises
+ * Copyright (C) 2021 - Samuel GOUGEON - Le Mans Université
  *
  * This file is hereby licensed under the terms of the GNU GPL v2.0,
  * pursuant to article 5.3.4 of the CeCILL v.2.1.
  * along with this program.
  *
  -->
-<refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:svg="http://www.w3.org/2000/svg" xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org" xml:lang="fr" xml:id="det">
+<refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink"
+          xmlns:svg="http://www.w3.org/2000/svg" xmlns:mml="http://www.w3.org/1998/Math/MathML"
+          xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org"
+          xml:lang="fr" xml:id="det">
     <refnamediv>
         <refname>det </refname>
         <refpurpose>déterminant  </refpurpose>
     </refnamediv>
     <refsynopsisdiv>
         <title>Séquence d'appel</title>
-        <synopsis>det(X)
-            [e,m]=det(X)
+        <synopsis>
+            d = det(X)
+            [e,m] = det(X)
         </synopsis>
     </refsynopsisdiv>
     <refsection>
         <title>Paramètres</title>
         <variablelist>
             <varlistentry>
-                <term>X  </term>
+                <term>X</term>
+                <listitem>
+                    matrice carrée de nombres ou polynômes ou fractions rationnelles (à coefficients)
+                    réels ou complexes. Matrices creuses acceptées.
+                    <para/>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>d</term>
                 <listitem>
-                    <para>matrice carrée réelle ou complexe (creuse ou pleine), polynomiale ou rationnelle
-                    </para>
+                    Scalaire du type de <varname>X</varname> : le déterminant de la matrice
+                    <varname>X</varname>. Si <varname>X</varname> est creuse,
+                    <varname>d</varname> est toujours dense.
+                    <para/>
                 </listitem>
             </varlistentry>
             <varlistentry>
-                <term>m  </term>
+                <term>m</term>
                 <listitem>
-                    <para>nombre réel ou complexe, mantisse du déterminant en base 10
-                    </para>
+                    nombre réel ou complexe : mantisse du déterminant en base 10, telle que
+                    <literal>abs(m) ∈ [1,10)</literal>. Argument non supporté lorsque <varname>X</varname>
+                    est polynômiale ou en fractions rationnelles.
+                    <para/>
                 </listitem>
             </varlistentry>
             <varlistentry>
-                <term>e  </term>
+                <term>e</term>
                 <listitem>
-                    <para>entier, exposant du déterminant en base 10
-                    </para>
+                    entier : la puissance de 10 du déterminant, telle que
+                    <literal>d = m * 10<superscript>e</superscript></literal>.
+                    Argument non supporté lorsque <varname>X</varname> est polynômiale ou
+                    en fractions rationnelles.
+                    <para/>
                 </listitem>
             </varlistentry>
         </variablelist>
     <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>.
+            <emphasis role="bold">d = det(X)</emphasis> calcule et donne le déterminant de la matrice
+            <varname>X</varname>.
+        </para>
+        <para>
+            Pour une matrice de polynômes ou de fractions rationelles, <literal>d=det(X)</literal>
+            utilise <literal>determ(..)</literal> dont l'algorithme est basé sur la transformée de
+            Fourier.
+            <literal>d=detr(X)</literal> peut être utilisée de manière alternative, utilisant
+            l'algorithme de Leverrier. Les deux méthodes produisent des résultats équivalents.
+            Pour une matrice de fractions rationnelles, neutraliser la simplification avec
+            <code>simp_mode(%f)</code> peut être requis pour obtenir deux résultats identiques.
+        </para>
+        <para>
+            <emphasis role="bold">[e, m] = det(X)</emphasis> peut être utilisé uniquement pour des
+            matrices de nombres. Cette syntaxe permet de remédier aux débordements numériques vers 0
+            ou vers l'infini, lorsque <literal>abs(d)</literal> est plus petit que
+            <literal>number_properties("tiny")</literal> ≈ 2.23 10<superscript>-308</superscript> ou
+            plus grand que <literal>number_properties("huge")</literal> ≈ 1.80 10<superscript>308</superscript>.
         </para>
         <para>
-            Pour les matrices polynomiales <literal>det(X)</literal> est équivalent à <literal>determ(X)</literal>.
+            Pour les matrices numériques denses, <literal>det(..)</literal> est basée sur les routines LAPACK
+            DGETRF pour les matrices réelles et  ZGETRF pour les matrices de nombres complexes.
         </para>
         <para>
-            Pour les matrices rationnelles <literal>det(X)</literal> est équivalent à <literal>detr(X)</literal>.
+            Pour les matrices numériques creuses, le déterminant est calculé d'après la factorisation
+            LU  de <varname>X</varname>, réalisée à l'aide de la bibliothèque UMFPACK.
         </para>
     </refsection>
     <refsection>
         <title>Exemples</title>
         <programlisting role="example"><![CDATA[
-x=poly(0,'x');
-det([x,1+x;2-x,x^2])
-w=ssrand(2,2,4);roots(det(systmat(w))),trzeros(w)   // zéros du système linéaire
-A=rand(3,3);
-det(A), prod(spec(A))
- ]]></programlisting>
+A = rand(3,3)*5;
+det(A)
+[e, m] = det(A)
+
+// Matrice de nombres complexes :
+// A = grand(3,3,"uin",0,10) + grand(3,3,"uin",0,10)*%i
+A = [3+%i, 9+%i*3, 9+%i ; 8+%i*8, 4+%i*3, 7+%i*7 ; 4, 6+%i*2, 6+%i*9]
+det(A)
+[e, m] = det(A)
+abs(m)  // dans [1, 10)
+     ]]></programlisting>
+        <screen><![CDATA[
+--> A = rand(3,3)*5;
+--> det(A)
+ ans  =
+  -10.805163
+
+--> [e, m] = det(A)
+ e  =
+   1.
+ m  =
+  -1.0805163
+
+--> // Matrice de nombres complexes :
+--> A = [3+%i,9+%i*3,9+%i;8+%i*8,4+%i*3,7+%i*7;4,6+%i*2,6+%i*9]
+ A  =
+   3. + i     9. + 3.i   9. + i
+   8. + 8.i   4. + 3.i   7. + 7.i
+   4. + 0.i   6. + 2.i   6. + 9.i
+
+--> det(A)
+ ans  =
+   745. - 225.i
+
+--> [e, m] = det(A)
+ e  =
+   2.
+ m  =
+   7.45 - 2.25i
+
+--> abs(m)  // dans [1, 10)
+ ans  =
+   7.7823518
+]]></screen>
+        <para/>
+        <para>
+            Déterminants très grands ou très petits : gestion des débordements numériques :
+        </para>
+        <programlisting role="example"><![CDATA[
+// Très grand déterminant :
+n = 1000;
+A = rand(n, n);
+det(A)
+[e, m] = det(A)
+
+// Très petit déterminant (d'une matrice encodée creuse) :
+A = (triu(sprand(n,n,1)) + diag(rand(1,n)))/1.5;
+det(A)
+prod(diag(A))
+[e, m] = det(A)
+A = A/2;
+det(A)
+[e, m] = det(A)
+     ]]></programlisting>
+        <screen><![CDATA[
+--> // Très grand déterminant :
+--> A = rand(n, n);
+--> det(A)
+ ans  =
+  -Inf
+
+--> [e, m] = det(A)  // -3.1199e743
+ e  =
+   743.
+ m  =
+  -3.1198687
+
+--> // Très petit déterminant (d'une matrice encodée creuse) :
+--> n = 1000;
+--> A = (triu(sprand(n,n,1)) + diag(rand(1,n)))/1.5;
+--> det(A)
+ ans  =
+   5.21D-236
+
+--> prod(diag(A))
+ ans  =
+   5.21D-236
+
+--> [e, m] = det(A)
+ e  =
+  -236.
+ m  =
+   5.2119757
+
+--> A = A/2;
+--> det(A)
+ ans  =
+   0.
+
+--> [e, m] = det(A)
+ e  =
+  -537.
+ m  =
+   4.8641473
+]]></screen>
+        <para/>
+        <para>
+            Déterminant de matrices polynômiales :
+        </para>
+        <programlisting role="example"><![CDATA[
+s = %s;
+det([s, 1+s ; 2-s, s^2])
+
+w = ssrand(2,2,4);
+roots(det(systmat(w))),trzeros(w)   //zeros of linear system
+     ]]></programlisting>
+        <screen><![CDATA[
+--> det([s, 1+s ; 2-s, s^2])
+ ans  =
+  -2 -s +s² +s³
+
+--> w = ssrand(2,2,4);
+--> roots(det(systmat(w))),trzeros(w)
+ ans  =
+  -3.1907522 + 0.i
+   2.3596502 + 0.i
+
+ ans  =
+   2.3596502 + 0.i
+  -3.1907522 + 0.i
+]]></screen>
     </refsection>
     <refsection role="see also">
         <title>Voir aussi</title>
@@ -81,6 +244,9 @@ det(A), prod(spec(A))
             <member>
                 <link linkend="determ">determ</link>
             </member>
+            <member>
+                <link linkend="simp_mode">simp_mode</link>
+            </member>
         </simplelist>
     </refsection>
     <refsection>
@@ -94,4 +260,15 @@ det(A), prod(spec(A))
             à partir de la décomposition LU de la librairie umfpack.
         </para>
     </refsection>
+    <refsection role="history">
+        <title>Historique</title>
+        <revhistory>
+            <revision>
+                <revnumber>6.1.1</revnumber>
+                <revdescription>
+                    [e,m]=det(X) syntax extended to sparse matrices.
+                </revdescription>
+            </revision>
+        </revhistory>
+    </refsection>
 </refentry>
diff --git a/scilab/modules/linear_algebra/help/ja_JP/matrix/det.xml b/scilab/modules/linear_algebra/help/ja_JP/matrix/det.xml
deleted file mode 100644 (file)
index f651842..0000000
+++ /dev/null
@@ -1,215 +0,0 @@
-<?xml version="1.0" encoding="UTF-8"?>
-
-<!--
- * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
- * Copyright (C) 2008 - INRIA
- *
- * Copyright (C) 2012 - 2016 - Scilab Enterprises
- *
- * This file is hereby licensed under the terms of the GNU GPL v2.0,
- * pursuant to article 5.3.4 of the CeCILL v.2.1.
- * This file was originally licensed under the terms of the CeCILL v2.1,
- * and continues to be available under such terms.
- * For more information, see the COPYING file which you should have received
- * along with this program.
- *
- -->
-
-<refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:svg="http://www.w3.org/2000/svg" xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org" xml:lang="ja" xml:id="det">
-
-    <refnamediv>
-
-        <refname>det</refname>
-
-        <refpurpose>行列式</refpurpose>
-
-    </refnamediv>
-
-    <refsynopsisdiv>
-
-        <title>呼出し手順</title>
-
-        <synopsis>det(X)
-
-            [e,m]=det(X)
-
-        </synopsis>
-
-    </refsynopsisdiv>
-
-    <refsection>
-
-        <title>引数</title>
-
-        <variablelist>
-
-            <varlistentry>
-
-                <term>X</term>
-
-                <listitem>
-
-                    <para>実数または複素正方行列, 多項式または有理行列.</para>
-
-                </listitem>
-
-            </varlistentry>
-
-            <varlistentry>
-
-                <term>m</term>
-
-                <listitem>
-
-                    <para>実数または複素数, 行列式の 10 を基底とする仮数</para>
-
-                </listitem>
-
-            </varlistentry>
-
-            <varlistentry>
-
-                <term>e</term>
-
-                <listitem>
-
-                    <para>整数, 行列式の 10 を基底とする指数</para>
-
-                </listitem>
-
-            </varlistentry>
-
-        </variablelist>
-
-    </refsection>
-
-    <refsection>
-
-        <title>説明</title>
-
-        <para>
-
-            <literal>det(X)</literal> (<literal>m*10^e</literal>)は,
-
-            正方行列<literal>X</literal>の行列式です.
-
-        </para>
-
-        <para>
-
-            多項式行列の場合,<literal>det(X)</literal> は
-
-            <literal>determ(X)</literal>と等しくなります.
-
-        </para>
-
-        <para>
-
-            有理数行列の場合, <literal>det(X)</literal> は
-
-            <literal>detr(X)</literal>と等しくなります.
-
-        </para>
-
-        <para>
-
-            <important>
-
-                <literal>det</literal> および <literal>detr</literal> 関数は
-
-                同じアルゴリズムを使用しません.
-
-                有理数関数の場合, <literal>det(X)</literal> は
-
-                <literal>determ</literal>関数に基づく
-
-                <literal>%r_det(X)</literal> でオーバーロードされます.
-
-                <literal>detr()</literal> は, Leverrier法を使用します.
-
-            </important>
-
-            <warning>
-
-                時々,
-
-                <literal>det</literal> および <literal>detr</literal> 関数は
-
-                有理数関数と異なる値を返す可能性があります.
-
-                このような場合,同じ結果を得るために,
-
-                有理数は<code>simp_mode(%f)</code>を使用することにより
-
-                有理数を簡単化するモードを無効にする必要があります.
-
-            </warning>
-
-        </para>
-
-    </refsection>
-
-    <refsection role="see also">
-
-        <title>参照</title>
-
-        <para>
-
-            det の計算は Lapack ルーチン DGETRF (実数行列の場合) および
-
-            ZGETRF (複素数の場合)に基づいています.
-
-        </para>
-
-        <para>
-
-            疎行列の場合, 行列式は umfpack ライブラリのLU分解により得られます.
-
-        </para>
-
-    </refsection>
-
-    <refsection>
-
-        <title>例</title>
-
-        <programlisting role="example"><![CDATA[
-x=poly(0,'x');
-det([x,1+x;2-x,x^2])
-w=ssrand(2,2,4);roots(det(systmat(w))),trzeros(w)   //線形システムのゼロ
-A=rand(3,3);
-det(A), prod(spec(A))
- ]]></programlisting>
-
-    </refsection>
-
-    <refsection>
-
-        <title>参照</title>
-
-        <simplelist type="inline">
-
-            <member>
-
-                <link linkend="detr">detr</link>
-
-            </member>
-
-            <member>
-
-                <link linkend="determ">determ</link>
-
-            </member>
-
-            <member>
-
-                <link linkend="simp_mode">simp_mode</link>
-
-            </member>
-
-        </simplelist>
-
-    </refsection>
-
-</refentry>
-
diff --git a/scilab/modules/linear_algebra/help/pt_BR/matrix/det.xml b/scilab/modules/linear_algebra/help/pt_BR/matrix/det.xml
deleted file mode 100644 (file)
index cea4772..0000000
+++ /dev/null
@@ -1,97 +0,0 @@
-<?xml version="1.0" encoding="UTF-8"?>
-<!--
- * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
- * Copyright (C) 2008 - INRIA
- *
- * Copyright (C) 2012 - 2016 - Scilab Enterprises
- *
- * This file is hereby licensed under the terms of the GNU GPL v2.0,
- * pursuant to article 5.3.4 of the CeCILL v.2.1.
- * This file was originally licensed under the terms of the CeCILL v2.1,
- * and continues to be available under such terms.
- * For more information, see the COPYING file which you should have received
- * along with this program.
- *
- -->
-<refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:svg="http://www.w3.org/2000/svg" xmlns:ns3="http://www.w3.org/1999/xhtml" xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org" xml:id="det" xml:lang="pt">
-    <refnamediv>
-        <refname>det</refname>
-        <refpurpose>determinante</refpurpose>
-    </refnamediv>
-    <refsynopsisdiv>
-        <title>Seqüência de Chamamento</title>
-        <synopsis>det(X)
-            [e,m]=det(X)
-        </synopsis>
-    </refsynopsisdiv>
-    <refsection>
-        <title>Parâmetros</title>
-        <variablelist>
-            <varlistentry>
-                <term>X</term>
-                <listitem>
-                    <para>matriz quadrada de reais ou complexos, matriz de polinômios ou
-                        de razões de polinômios
-                    </para>
-                </listitem>
-            </varlistentry>
-            <varlistentry>
-                <term>m</term>
-                <listitem>
-                    <para>número real ou complexo, a mantissa de base 10 do
-                        determinante
-                    </para>
-                </listitem>
-            </varlistentry>
-            <varlistentry>
-                <term>e</term>
-                <listitem>
-                    <para>inteiro, o expoente de base 10 do determinante</para>
-                </listitem>
-            </varlistentry>
-        </variablelist>
-    </refsection>
-    <refsection>
-        <title>Descrição</title>
-        <para>
-            <literal>det(X)</literal> ( <literal>m*10^e</literal> é o
-            determinante da matriz quadrada <literal>X)</literal>.
-        </para>
-        <para>
-            Para uma matriz de polinômios, <literal>det(X)</literal> é
-            equivalente a <literal>determ(X)</literal>.
-        </para>
-        <para>
-            Para matrizes de razões de polinômios <literal>det(X)</literal> é
-            equivalente a <literal>detr(X)</literal>.
-        </para>
-    </refsection>
-    <refsection>
-        <title>Referências</title>
-        <para>As computações da função det são baseadas nas rotinas do LAPACK
-            DGETRF para matrizes de reais e ZGETRF para o caso de matrizes de
-            complexos.
-        </para>
-    </refsection>
-    <refsection>
-        <title>Exemplos</title>
-        <programlisting role="example"><![CDATA[
-x=poly(0,'x');
-det([x,1+x;2-x,x^2])
-w=ssrand(2,2,4);roots(det(systmat(w))),trzeros(w)   //zeros do sistema linear
-A=rand(3,3);
-det(A), prod(spec(A))
- ]]></programlisting>
-    </refsection>
-    <refsection>
-        <title> Ver Também</title>
-        <simplelist type="inline">
-            <member>
-                <link linkend="detr">detr</link>
-            </member>
-            <member>
-                <link linkend="determ">determ</link>
-            </member>
-        </simplelist>
-    </refsection>
-</refentry>
diff --git a/scilab/modules/linear_algebra/help/ru_RU/matrix/det.xml b/scilab/modules/linear_algebra/help/ru_RU/matrix/det.xml
new file mode 100644 (file)
index 0000000..fd6fc26
--- /dev/null
@@ -0,0 +1,264 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<!--
+ * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+ * Copyright (C) 2008 - INRIA
+ * Copyright (C) 2012 - 2016 - Scilab Enterprises
+ * Copyright (C) 2021 - Samuel GOUGEON - Le Mans Université
+ *
+ * This file is hereby licensed under the terms of the GNU GPL v2.0,
+ * pursuant to article 5.3.4 of the CeCILL v.2.1.
+ * This file was originally licensed under the terms of the CeCILL v2.1,
+ * and continues to be available under such terms.
+ * For more information, see the COPYING file which you should have received
+ * along with this program.
+ *
+ -->
+<refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink"
+          xmlns:svg="http://www.w3.org/2000/svg" xmlns:mml="http://www.w3.org/1998/Math/MathML"
+          xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org"
+          xml:lang="ru" xml:id="det">
+    <refnamediv>
+        <refname>det</refname>
+        <refpurpose>определитель квадратной матрицы</refpurpose>
+    </refnamediv>
+    <refsynopsisdiv>
+        <title>Синтаксис</title>
+        <synopsis>
+            d = det(X)
+            [e,m] = det(X)
+        </synopsis>
+    </refsynopsisdiv>
+    <refsection>
+        <title>Аргументы</title>
+        <variablelist>
+            <varlistentry>
+                <term>X</term>
+                <listitem>
+                    квадратная матрица вещественных или комплексных чисел, полиномов или
+                    дробно-рациональных выражений. Принимаются разрежённо-кодированные
+                    матрицы.
+                    <para/>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>d</term>
+                <listitem>
+                    Скаляр типа <varname>X</varname>: определитель <varname>X</varname>.
+                    Если <varname>X</varname> является разрежённо-кодированной, то
+                    <varname>d</varname> является плотной.
+                    <para/>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>m</term>
+                <listitem>
+                    вещественное или комплексное число: определитель мантиссы по основанию 10
+                    <literal>abs(m) ∈ [1,10)</literal>. Не поддерживается для полиномиального
+                    или дробно-рационального <varname>X</varname>.
+                    <para/>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>e</term>
+                <listitem>
+                    целое: определитель экспоненты по основанию 10, такой, что
+                    <literal>d = m * 10<superscript>e</superscript></literal>.
+                    Не поддерживается для полиномиального или дробно-рационального
+                    <varname>X</varname>.
+                    <para/>
+                </listitem>
+            </varlistentry>
+        </variablelist>
+    </refsection>
+    <refsection>
+        <title>Описание</title>
+        <para>
+            <emphasis role="bold">d = det(X)</emphasis> выдаёт определитель матрицы
+            <varname>X</varname>.
+        </para>
+        <para>
+            Для полиномиальной или дробно-рациональной матрицы <literal>d=det(X)</literal>
+            использует <literal>determ(..)</literal>, чей алгоритм основан на БПФ.
+            Альтернативно может использоваться <literal>d=detr(X)</literal>, основанный
+            на алгоритме Леверье. Оба метода выдают одинаковые результаты. Для
+            дробно-рациональных матриц может потребоваться отключение <code>simp_mode(%f)</code>
+            для получения идентичных результатов.
+        </para>
+        <para>
+            <emphasis role="bold">[e, m] = det(X)</emphasis> может использоваться только
+            для матрицы чисел. Этот синтаксис позволяет преодолеть переполнение снизу
+            или сверху, когда <literal>abs(d)</literal> меньше, чем
+            <literal>number_properties("tiny")</literal> ≈ 2,23 10<superscript>-308</superscript>
+            либо больше, чем
+            <literal>number_properties("huge")</literal> ≈ 1,80 10<superscript>308</superscript>.
+        </para>
+        <para>
+            Для плотных матриц <literal>det(..)</literal> основан на процедуре из Lapack
+            DGETRF для вещественных матриц и ZGETRF для комплексных.
+        </para>
+        <para>
+            Для разрежённых матриц определитель получается из LU-разложения, благодаря
+            библиотеке umfpack.
+        </para>
+    </refsection>
+    <refsection>
+        <title>Примеры</title>
+        <programlisting role="example"><![CDATA[
+A = rand(3,3)*5;
+det(A)
+[e, m] = det(A)
+
+// Матрица комплексных чисел:
+// A = grand(3,3,"uin",0,10) + grand(3,3,"uin",0,10)*%i
+A = [3+%i, 9+%i*3, 9+%i ; 8+%i*8, 4+%i*3, 7+%i*7 ; 4, 6+%i*2, 6+%i*9]
+det(A)
+[e, m] = det(A)
+abs(m)  // in [1, 10)
+     ]]></programlisting>
+        <screen><![CDATA[
+--> A = rand(3,3)*5;
+--> det(A)
+ ans  =
+  -10.805163
+
+--> [e, m] = det(A)
+ e  =
+   1.
+ m  =
+  -1.0805163
+
+--> // Matrix of complex numbers:
+--> A = [3+%i, 9+%i*3, 9+%i ; 8+%i*8, 4+%i*3, 7+%i*7 ; 4, 6+%i*2, 6+%i*9]
+ A  =
+   3. + i     9. + 3.i   9. + i
+   8. + 8.i   4. + 3.i   7. + 7.i
+   4. + 0.i   6. + 2.i   6. + 9.i
+
+--> det(A)
+ ans  =
+   745. - 225.i
+
+--> [e, m] = det(A)
+ e  =
+   2.
+ m  =
+   7.45 - 2.25i
+
+--> abs(m)  // in [1, 10)
+ ans  =
+   7.7823518
+]]></screen>
+        <para/>
+        <para>
+            Очень большой или маленький определители: обработка переполнения снизу и сверху:
+        </para>
+        <programlisting role="example"><![CDATA[
+// Очень большой определитель:
+n = 1000;
+A = rand(n, n);
+det(A)
+[e, m] = det(A)
+
+// Очень маленький определитель (разрежённой кодированной матрицы):
+A = (triu(sprand(n,n,1)) + diag(rand(1,n)))/1.5;
+det(A)
+prod(diag(A))
+[e, m] = det(A)
+A = A/2;
+det(A)
+[e, m] = det(A)
+     ]]></programlisting>
+        <screen><![CDATA[
+--> // Very big determinant:
+--> A = rand(n, n);
+--> det(A)
+ ans  =
+  -Inf
+
+--> [e, m] = det(A)  // -3.1199e743
+ e  =
+   743.
+ m  =
+  -3.1198687
+
+--> // Very small determinant (of a sparse-encoded matrix):
+--> n = 1000;
+--> A = (triu(sprand(n,n,1)) + diag(rand(1,n)))/1.5;
+--> det(A)
+ ans  =
+   5.21D-236
+
+--> prod(diag(A))
+ ans  =
+   5.21D-236
+
+--> [e, m] = det(A)
+ e  =
+  -236.
+ m  =
+   5.2119757
+
+--> A = A/2;
+--> det(A)
+ ans  =
+   0.
+
+--> [e, m] = det(A)
+ e  =
+  -537.
+ m  =
+   4.8641473
+]]></screen>
+        <para/>
+        <para>
+            Определитель полиномиальной матрицы:
+        </para>
+        <programlisting role="example"><![CDATA[
+s = %s;
+det([s, 1+s ; 2-s, s^2])
+
+w = ssrand(2,2,4);
+roots(det(systmat(w))),trzeros(w)   //нули линейной системы
+     ]]></programlisting>
+        <screen><![CDATA[
+--> det([s, 1+s ; 2-s, s^2])
+ ans  =
+  -2 -s +s² +s³
+
+--> w = ssrand(2,2,4);
+--> roots(det(systmat(w))),trzeros(w)
+ ans  =
+  -3.1907522 + 0.i
+   2.3596502 + 0.i
+
+ ans  =
+   2.3596502 + 0.i
+  -3.1907522 + 0.i
+]]></screen>
+    </refsection>
+    <refsection role="see also">
+        <title>Смотрите также</title>
+        <simplelist type="inline">
+            <member>
+                <link linkend="detr">detr</link>
+            </member>
+            <member>
+                <link linkend="determ">determ</link>
+            </member>
+            <member>
+                <link linkend="simp_mode">simp_mode</link>
+            </member>
+        </simplelist>
+    </refsection>
+    <refsection role="history">
+        <title>История</title>
+        <revhistory>
+            <revision>
+                <revnumber>6.1.1</revnumber>
+                <revdescription>
+                    Синтаксис [e,m]=det(X) расширен до разрежённых матриц.
+                </revdescription>
+            </revision>
+        </revhistory>
+    </refsection>
+</refentry>
index 0b7b78f..1230371 100644 (file)
@@ -1,12 +1,13 @@
-// <-- CLI SHELL MODE -->
-// <-- NO CHECK REF -->
 // =============================================================================
 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
 // Copyright (C) ????-2008 - INRIA Michael Baudin
 // Copyright (C) 2015 - Scilab Enterprises - John Gliksberg
+// Copyright (C) 2021 - Le Mans Université - Samuel GOUGEON
 //
 //  This file is distributed under the same license as the Scilab package.
 // =============================================================================
+// <-- CLI SHELL MODE -->
+// <-- NO CHECK REF -->
 
 //==========================================================================
 //==============================   det        ==============================
@@ -116,5 +117,79 @@ b
 b
 b
 b];
-assert_checkalmostequal(det(A), 0);
+assert_checkalmostequal(det(A), 0, %eps, norm(A)*%eps);
+
+a = sprand(1000,1000,0.1); a(:,1:2) = 1;
+assert_checkequal(det(a), 0);
 
+// Underflow and overflow management for sparse matrices. [e, m] = det(s)  syntax
+// -------------------------------------------------------------------------------
+// http://bugzilla.scilab.org/16636 :
+// With a real matrix
+// ..................
+// Underflow:
+n = 75;
+while n < 3000
+    s = sparse(triu(rand(n,n)));
+    ref = prod(diag(s));
+    if n < 1000
+        assert_checkalmostequal(det(s), ref, 10*n*%eps);
+    else
+        assert_checkequal(det(s), 0);
+    end
+    [e, m] = det(s);
+    ref = sum(log10(full(diag(s))));
+    assert_checkequal(e, int(ref)-1);
+    assert_checkalmostequal(m, 10^(ref-int(ref)+1), 10*n*%eps);
+    n = n*2;
+end
+// Overflow:
+n = 75;
+while n < 3000
+    s = sparse(triu(rand(n,n)))*6;
+    ref = prod(diag(s));
+    if n < 1000
+        assert_checkalmostequal(det(s), ref, 10*n*%eps);
+    else
+        assert_checkequal(det(s), %inf);
+    end
+    [e, m] = det(s);
+    ref = sum(log10(full(diag(s))));
+    assert_checkequal(e, int(ref));
+    assert_checkalmostequal(m, 10^(ref-int(ref)), 10*n*%eps);
+    n = n*2;
+end
+// With a complex matrix
+// .....................
+// Underflow:
+n = 75;
+while n < 3000
+    s = sparse(triu(complex(rand(n,n),rand(n,n))));
+    ref = prod(diag(s));
+    if n < 2000
+        assert_checkalmostequal(det(s), ref, 100*n*%eps);
+    else
+        assert_checkequal(det(s), 0*%i);
+    end
+    [e, m] = det(s);
+    [eref, mref] = det(full(s));
+    assert_checkequal(e, eref);
+    assert_checkalmostequal(m, mref, 100*n*%eps);
+    n = n*2;
+end
+// Overflow:
+n = 75;
+while n < 3000
+    s = sparse(triu(complex(rand(n,n),rand(n,n))))*6;
+    ref = prod(diag(s));
+    if n < 600
+        assert_checkalmostequal(det(s), ref, 100*n*%eps);
+    else
+        assert_checkequal(abs(real(det(s))), %inf);
+    end
+    [e, m] = det(s);
+    [eref, mref] = det(full(s));
+    assert_checkequal(e, eref);
+    assert_checkalmostequal(m, mref, 200*n*%eps);
+    n = n*2;
+end
index 113fb98..a5da027 100644 (file)
@@ -1,6 +1,7 @@
 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
 // Copyright (C) 2013 - Scilab Enterprises - Charlotte HECQUET
 // Copyright (C) 2012 - 2016 - Scilab Enterprises
+// Copyright (C) 2021 - Samuel GOUGEON - Le Mans Université
 //
 // This file is hereby licensed under the terms of the GNU GPL v2.0,
 // pursuant to article 5.3.4 of the CeCILL v.2.1.
@@ -18,19 +19,46 @@ function [res1, res2] = %sp_det(A)
         end
         return
     end
-    hand = umf_lufact(A); //umfpack is used for complex sparse matrix
+    if nnz(A)==0 then
+        [res1, res2] = (0,0)
+        return
+    end
+
+    wstatus = warning("query"); warning("off")  // to avoid the "is singular" message
+    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;
+    umf_ludel(hand);
+    warning(wstatus)
+
+    dU = clean(diag(U), 0, %eps);
+    r = clean(r, 0, %eps);
+    if isreal(A) then
+        p = r .* dU;
+        tmp = sum(log10(full(abs(p))))
+        res1 = int(tmp)
+        if res1 < 0
+            res1 = res1 - 1
+        end
+        res2 = 10^(tmp-res1) * prod(sign(p))
+    else
+        tmp = sum(log([r full(dU)]))
+        phase = imag(tmp)
+        e10 = real(tmp)/log(10) + log10(abs(cos(phase)))
+        res1 = int(e10)
+        res2 = 10^(real(tmp)/log(10)-res1) * exp(%i*phase)
         while abs(res2) >= 10
-            if abs(res2) < 1 then
-                break;
-            end
-            res2 = res2 / 10;
-            res1 = res1 + 1;
+            res2 = res2 / 10
+            res1 = res1 + 1
+        end
+        while abs(res2)<1 & abs(res2)<>0
+            res2 = res2 * 10
+            res1 = res1 - 1
         end
     end
-    umf_ludel(hand);
+    if res1 == -%inf
+        [res1, res2] = (0,0)
+    end
+    if lhs == 1 then
+        res1 = res2 * 10^res1
+    end
 endfunction