* Bug 16369 fixed: sparse right divisions / restored 43/21443/5
Samuel GOUGEON [Sun, 15 Mar 2020 02:20:57 +0000 (03:20 +0100)]
  http://bugzilla.scilab.org/16369

  After p = a * b:   a = p / b
  * restored for p or/and b sparses, and b and p real
  * added for b square and complex

  After p = a * b:   b = a \ p
  * added for a square and complex

  Unit tests added

Change-Id: I0738f4e542f6c527ab74d1debfc0b5570cb9d8b2

scilab/CHANGES.md
scilab/modules/overloading/macros/%s_r_sp.sci [new file with mode: 0644]
scilab/modules/overloading/macros/%sp_l_sp.sci
scilab/modules/overloading/macros/%sp_r_s.sci [new file with mode: 0644]
scilab/modules/overloading/macros/%sp_r_sp.sci [new file with mode: 0644]
scilab/modules/sparse/tests/unit_tests/sparse_divide.tst [new file with mode: 0644]

index 96085d1..0c16a53 100644 (file)
@@ -270,6 +270,7 @@ Bug Fixes
 * [#16350](https://bugzilla.scilab.org/16350): in if/while conditions, the empty sparse boolean was considered as TRUE.
 * [#16365](https://bugzilla.scilab.org/16365): `median(m,"r")` and `median(m,"c")` yielded wrong results (6.1.0 regression)
 * [#16366](https://bugzilla.scilab.org/16366): `plot([0 1], ":")` plotted a dash-dotted curve instead of a dotted one.
+* [#16369](https://bugzilla.scilab.org/16369): Right divisions / involving one or two sparse numerical matrices were no longer supported.
 * [#16397](https://bugzilla.scilab.org/16397): display of long (real) column vectors was slow (regression).
 * [#16399](https://bugzilla.scilab.org/16399): `mtlb_zeros([])` was crashing Scilab.
 * [#16401](https://bugzilla.scilab.org/16401): global `external_object_java` class was crashing Scilab.
diff --git a/scilab/modules/overloading/macros/%s_r_sp.sci b/scilab/modules/overloading/macros/%s_r_sp.sci
new file mode 100644 (file)
index 0000000..6045801
--- /dev/null
@@ -0,0 +1,18 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 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.
+
+function x = %s_r_sp(a, b)
+    // a / b , a dense, b sparse
+
+    x = a / full(b)
+
+endfunction
index 9cb3ad1..fe31068 100644 (file)
@@ -1,7 +1,7 @@
 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
 // Copyright (C) INRIA
-//
 // Copyright (C) 2012 - 2016 - Scilab Enterprises
+// Copyright (C) 2020 - Samuel GOUGEON
 //
 // 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.
 // For more information, see the COPYING file which you should have received
 // along with this program.
 
-function x = %sp_l_sp(a,b)
-    // a\b , a sparse b sparse
+function b = %sp_l_sp(a, p)
+    // a sparse, p sparse, b sparse
+    // b = a \ p   such that p = a * b
 
-    [ma,na]=size(a)
-    [mb,nb]=size(b)
-    if ma<>mb then
+    [ma, na] = size(a)
+    [mp, np] = size(p)
+    if ma <> mp then
         msg = _("%s: Arguments #%d and #%d: Same numbers of rows expected.\n")
         error(msprintf(msg, "%sp_l_sp", 1, 2))
     end
-    if ma<>na then
-        b=a'*b;a=a'*a
+    if issquare(a)
+        if det(a) <> 0 then
+            b = inv(a) * p
+            return
+        end
+    else
+        p = a.' * p
+        a = a.' * a
     end
 
-    [h,rk]=lufact(a)
-    if rk<min(ma,na) then warning("deficient rank: rank = "+string(rk)),end
-    x=[]
-    for k=1:nb
-        x=[x,sparse(lusolve(h,full(b(:,k))))]
+    [h,rk] = lufact(a)
+    if rk < min(ma,na) then
+        warning(msprintf(_("sparse \ sparse: Deficient rank: rank = %d"),rk))
+    end
+    b = []
+    for k = 1:np
+        b = [b, sparse(lusolve(h,full(p(:,k))))]
     end
     ludel(h)
 endfunction
diff --git a/scilab/modules/overloading/macros/%sp_r_s.sci b/scilab/modules/overloading/macros/%sp_r_s.sci
new file mode 100644 (file)
index 0000000..833e947
--- /dev/null
@@ -0,0 +1,18 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 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.
+
+function x = %sp_r_s(a, b)
+    // a / b , a sparse, b dense
+
+    x = full(a) / b
+
+endfunction
diff --git a/scilab/modules/overloading/macros/%sp_r_sp.sci b/scilab/modules/overloading/macros/%sp_r_sp.sci
new file mode 100644 (file)
index 0000000..ef44524
--- /dev/null
@@ -0,0 +1,43 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA
+// Copyright (C) 2012 - 2016 - Scilab Enterprises
+// Copyright (C) 2020 - Samuel GOUGEON
+//
+// 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.
+
+function a = %sp_r_sp(p, b)
+    // a = p / b, such that p = a * b
+    // p sparse, b sparse
+    // not-square complex b not supported
+
+    [rp, cp] = size(p)
+    [rb, cb] = size(b)
+    if cp <> cb then
+        msg = _("%s: Arguments #%d and #%d: Same numbers of columns expected.\n")
+        error(msprintf(msg, "%sp_r_sp", 1, 2))
+    end
+    if issquare(b)
+        if det(b) <> 0 then
+            a = p * inv(b)
+            return
+        end
+    else
+        p = p * b.'
+        b = b * b.'   // Makes b square
+    end
+
+    [h,rk] = lufact(b.')
+    if rk < min(rb,cb) then
+        warning(msprintf(_("sparse / sparse: Deficient rank: rank = %d"),rk))
+    end
+    a = []
+    for k = 1:rp
+        a = [a ; sparse(lusolve(h,full(p(k,:)).').')]
+    end
+    ludel(h)
+endfunction
diff --git a/scilab/modules/sparse/tests/unit_tests/sparse_divide.tst b/scilab/modules/sparse/tests/unit_tests/sparse_divide.tst
new file mode 100644 (file)
index 0000000..d5b9227
--- /dev/null
@@ -0,0 +1,111 @@
+
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2020 - Samuel GOUGEON
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+
+// <-- CLI SHELL MODE -->
+// <-- NO CHECK REF -->
+
+A = list(rand(4,5), rand(4,5) + rand(4,5)*%i, ..
+         sprand(4,5,0.5), sprand(4,5,0.5) + sprand(4,5,0.5)*%i);
+B = list(rand(5,6), rand(5,6)+rand(5,6)*%i, ..
+         sprand(5,6,0.5), sprand(5,6,0.5) + sprand(5,6,0.5)*%i);
+
+// ----------------
+// RIGHT DIVISION /
+// ----------------
+// p = a * b   => a = p / b
+// b rectangular
+// -------------
+for a = A
+    isspa = issparse(a)
+    for b = B
+        isspb = issparse(b)
+        // complex sparse / complex sparse is not supported
+        if ~isspa & ~isspb | (isspa & isspb & (~isreal(a) | ~isreal(b)))
+            continue
+        end
+        p = a * b;
+        a2 = p / b;
+        if ~issparse(p) | ~isspb
+            assert_checkfalse(issparse(a2));
+        else
+            assert_checktrue(issparse(a2));
+        end
+        assert_checkequal(size(a2), size(a));
+        assert_checkalmostequal(clean(a2*b,%eps,%eps), p);
+    end
+end
+// Square b
+// --------
+SB = list(rand(5,5), rand(5,5)+rand(5,5)*%i, ..
+         sprand(5,5,0.5), sprand(5,5,0.5) + sprand(5,5,0.5)*%i);
+for a = A
+    isspa = issparse(a)
+    for b = SB
+        isspb = issparse(b)
+        if ~isspa & ~isspb
+            continue
+        end
+        p = a * b;
+        a2 = p / b;
+        if ~issparse(p) | ~isspb
+            assert_checkfalse(issparse(a2));
+        else
+            assert_checktrue(issparse(a2));
+        end
+        assert_checkequal(size(a2), size(a));
+        assert_checkalmostequal(clean(a2*b,%eps,%eps), p);
+    end
+end
+
+// ---------------
+// LEFT DIVISION \
+// ---------------
+// p = a * b   => b = a \ p
+// a rectangular
+// -------------
+for a = A
+    isspa = issparse(a)
+    for b = B
+        isspb = issparse(b)
+        // complex sparse \ complex sparse is not supported
+        if ~isspa & ~isspb | (isspa & isspb & (~isreal(a) | ~isreal(b)))
+            continue
+        end
+        p = a * b;
+        b2 = a \ p;
+        if ~isspa | ~issparse(p)
+            assert_checkfalse(issparse(b2));
+        else
+            assert_checktrue(issparse(b2));
+        end
+        assert_checkequal(size(b2), size(b));
+        assert_checkalmostequal(clean(a*b2,%eps,%eps), p);
+    end
+end
+// Square a
+// --------
+// A = SB
+for a = SB
+    isspa = issparse(a)
+    for b = B
+        isspb = issparse(b)
+        if ~isspa & ~isspb
+            continue
+        end
+        p = a * b;
+        b2 = a \ p;
+        if ~isspa | ~issparse(p)
+            assert_checkfalse(issparse(b2));
+        else
+            assert_checktrue(issparse(b2));
+        end
+        assert_checkequal(size(b2), size(b));
+        assert_checkalmostequal(clean(a*b2,%eps,%eps), p);
+    end
+end
+