* Bug #12137 fixed - eigs(A,B) returned incorrect result for sparse matrices. 95/9895/3
Guillaume Horel [Mon, 3 Dec 2012 11:07:51 +0000 (12:07 +0100)]
Change-Id: I73fcb0f8f3dbbd1955e687f018ddba79bc1e1ef7

scilab/CHANGES_5.4.X
scilab/modules/arnoldi/macros/eigs.sci
scilab/modules/arnoldi/tests/nonreg_tests/bug_12137.dia.ref [new file with mode: 0644]
scilab/modules/arnoldi/tests/nonreg_tests/bug_12137.tst [new file with mode: 0644]

index 5667c1a..a25b589 100644 (file)
@@ -219,6 +219,8 @@ Bug fixes
 
 * Bug #12135 fixed - stacksize('max') failed silently.
 
+* Bug #12137 fixed - eigs(A,B) returned incorrect result for sparse matrices.
+
 * Bug #12140 fixed - csvRead fails when comma is used as decimal mark.
 
 * Bug #12166 fixed - There was a bad label with drawaxis.
index ec4f5ea..c745747 100644 (file)
@@ -554,7 +554,15 @@ function [res_d, res_v] = speigs(A, %_B, nev, which, maxiter, tol, ncv, cholB, r
         if(~Breal)
             error(msprintf(gettext("%s: Impossible to use the Cholesky factorisation with complex sparse matrices.\n"), "eigs"));
         else
-            [R,P] = spchol(%_B);
+            if(issparse(%_B))
+                [R, P] = spchol(%_B);
+                perm = spget(P);
+                perm = perm(:,2);
+                iperm = spget(P');
+                iperm = iperm(:,2);
+            else
+                R = chol(%_B);
+            end
         end
         Rprime = R';
     end
@@ -628,10 +636,18 @@ function [res_d, res_v] = speigs(A, %_B, nev, which, maxiter, tol, ncv, cholB, r
 
         if(ido == -1 | ido == 1 | ido == 2)
             if(iparam(7) == 1)
-                if(matB == 0)
-                    workd(ipntr(2):ipntr(2)+nA-1) = A * workd(ipntr(1):ipntr(1)+nA-1);
+                if(ido==2)
+                    workd(ipntr(2):ipntr(2)+nA-1) = workd(ipntr(1):ipntr(1)+nA-1);
                 else
-                    workd(ipntr(2):ipntr(2)+nA-1) = inv(Rprime) * A * inv(R) * workd(ipntr(1):ipntr(1)+nA-1);
+                    if(matB == 0)
+                        workd(ipntr(2):ipntr(2)+nA-1) = A * workd(ipntr(1):ipntr(1)+nA-1);
+                    else
+                        if(issparse(%_B))
+                            workd(ipntr(2):ipntr(2)+nA-1) = R \ A(iperm,iperm) * (Rprime \ workd(ipntr(1):ipntr(1)+nA-1));
+                        else
+                            workd(ipntr(2):ipntr(2)+nA-1) = Rprime \ (A * (R \ workd(ipntr(1):ipntr(1)+nA-1)));
+                        end
+                    end
                 end
             elseif(iparam(7) == 3)
                 if(matB == 0)
@@ -722,6 +738,14 @@ function [res_d, res_v] = speigs(A, %_B, nev, which, maxiter, tol, ncv, cholB, r
             end
         end
     end
+    if rvec & iparam(7)==1 & matB<>0
+        if issparse(B) then
+             res_v = Rprime \ res_v;
+             res_v = res_v(perm,:);
+        else
+             res_v = R \ res_v;
+        end
+    end
 endfunction
 
 
diff --git a/scilab/modules/arnoldi/tests/nonreg_tests/bug_12137.dia.ref b/scilab/modules/arnoldi/tests/nonreg_tests/bug_12137.dia.ref
new file mode 100644 (file)
index 0000000..5a849b1
--- /dev/null
@@ -0,0 +1,34 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2012 - Scilab Enterprises - Sylvestre Ledru
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+// <-- CLI SHELL MODE -->
+// <-- Non-regression test for bug 12137 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=12137
+//
+// <-- Short Description -->
+//    eigs(A,B) returns incorrect result for sparse matrices
+// =============================================================================
+A=sparse(toeplitz([1 0 0 3 0 5 0 3 0 0 1]));
+B=[
+3.,  0.,  0.,  2.,  0.,  0.,  2.,  0.,  2.,  0.,  0. ;
+0.,  5.,  4.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0. ;
+0.,  4.,  5.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0. ;
+2.,  0.,  0.,  3.,  0.,  0.,  2.,  0.,  2.,  0.,  0. ;
+0.,  0.,  0.,  0. , 5.,  0.,  0.,  0.,  0.,  0.,  4. ;
+0.,  0.,  0.,  0.,  0.,  4.,  0.,  3.,  0.,  3.,  0. ;
+2.,  0.,  0.,  2.,  0.,  0.,  3.,  0.,  2.,  0.,  0. ;
+0.,  0.,  0.,  0.,  0.,  3.,  0.,  4.,  0.,  3.,  0. ;
+2.,  0.,  0.,  2.,  0.,  0.,  2.,  0.,  3.,  0.,  0. ;
+0.,  0.,  0.,  0.,  0.,  3.,  0.,  3.,  0.,  4.,  0. ;
+0.,  0.,  0.,  0.,  4.,  0.,  0.,  0.,  0.,  0.,  5.];
+B=sparse(B);
+C=eigs(A,B,10);
+[a b] = spec(full(A), full(B));
+C1=gsort(a./b, "g", "i");
+C1(1)=[]; // remove the small value
+assert_checkalmostequal(C1, C);
diff --git a/scilab/modules/arnoldi/tests/nonreg_tests/bug_12137.tst b/scilab/modules/arnoldi/tests/nonreg_tests/bug_12137.tst
new file mode 100644 (file)
index 0000000..4d593e1
--- /dev/null
@@ -0,0 +1,39 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2012 - Scilab Enterprises - Sylvestre Ledru
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+
+// <-- CLI SHELL MODE -->
+
+// <-- Non-regression test for bug 12137 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=12137
+//
+// <-- Short Description -->
+//    eigs(A,B) returns incorrect result for sparse matrices
+// =============================================================================
+
+A=sparse(toeplitz([1 0 0 3 0 5 0 3 0 0 1]));
+
+B=[
+3.,  0.,  0.,  2.,  0.,  0.,  2.,  0.,  2.,  0.,  0. ;
+0.,  5.,  4.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0. ;
+0.,  4.,  5.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0. ;
+2.,  0.,  0.,  3.,  0.,  0.,  2.,  0.,  2.,  0.,  0. ;
+0.,  0.,  0.,  0. , 5.,  0.,  0.,  0.,  0.,  0.,  4. ;
+0.,  0.,  0.,  0.,  0.,  4.,  0.,  3.,  0.,  3.,  0. ;
+2.,  0.,  0.,  2.,  0.,  0.,  3.,  0.,  2.,  0.,  0. ;
+0.,  0.,  0.,  0.,  0.,  3.,  0.,  4.,  0.,  3.,  0. ;
+2.,  0.,  0.,  2.,  0.,  0.,  2.,  0.,  3.,  0.,  0. ;
+0.,  0.,  0.,  0.,  0.,  3.,  0.,  3.,  0.,  4.,  0. ;
+0.,  0.,  0.,  0.,  4.,  0.,  0.,  0.,  0.,  0.,  5.];
+B=sparse(B);
+C=eigs(A,B,10);
+
+[a b] = spec(full(A), full(B));
+C1=gsort(a./b, "g", "i");
+C1(1)=[]; // remove the small value
+assert_checkalmostequal(C1, C);