* Bug #12238 fixed - [d v] = eigs(A) was broken for sparse matrices.

test_run("arnoldi","bug_12238",["no_check_error_output"])
test_run("arnoldi",[],["no_check_error_output"])

Change-Id: I3cb71420c734e9e97cff3898906c0f0aa86fd41d

diff --git a/scilab/CHANGES_5.4.X b/scilab/CHANGES_5.4.X
index b8b0f06..c57b2ff 100644 (file)
--- a/scilab/CHANGES_5.4.X
+++ b/scilab/CHANGES_5.4.X
@@ -305,6 +305,8 @@ Bug fixes
 
 * Bug #12235 fixed - Matplot did not update on color_map change.
 
+* Bug #12238 fixed - [d v] = eigs(A) was broken for sparse matrices.
+
 * Bug #12247 fixed - Fix a typo in some error messages.
 
 

diff --git a/scilab/modules/arnoldi/macros/eigs.sci b/scilab/modules/arnoldi/macros/eigs.sci
index a8af185..b73456b 100644 (file)
--- a/scilab/modules/arnoldi/macros/eigs.sci
+++ b/scilab/modules/arnoldi/macros/eigs.sci
@@ -545,7 +545,7 @@ function [res_d, res_v] = speigs(A, %_B, nev, which, maxiter, tol, ncv, cholB, r
 
     if(cholB)
         if(or(triu(%_B) <> %_B))
-            error(msprintf(gettext("%s: Wrong type for input argument #%d: B must be symmetric or hermitian, definite, semi positive.\n"), "eigs", 2));
+            error(msprintf(gettext("%s: Wrong type for input argument #%d: if opts.cholB is true, B must be upper triangular.\n"), "eigs", 2));
         end
         if(issparse(%_B)) //sparse cholesky decomposition is reversed...
             Rprime = %_B;
@@ -605,7 +605,7 @@ function [res_d, res_v] = speigs(A, %_B, nev, which, maxiter, tol, ncv, cholB, r
             z = zeros(nA, nev);
         else
             lworkl = 3 * ncv * (ncv + 2);
-            v = zeros(nA, ncv);  
+            v = zeros(nA, ncv);
             workl = zeros(lworkl, 1);
             workd = zeros(3 * nA, 1);
             dr = zeros(nev+1, 1);
@@ -646,7 +646,7 @@ 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(ido==2)
+                if(ido == 2)
                     workd(ipntr(2):ipntr(2)+nA-1) = workd(ipntr(1):ipntr(1)+nA-1);
                 else
                     if(matB == 0)
@@ -704,8 +704,8 @@ function [res_d, res_v] = speigs(A, %_B, nev, which, maxiter, tol, ncv, cholB, r
             end
         end
     end
-    if(iparam(7)==3)
-        umf_ludel(Lup);
+    if(iparam(7) == 3)
+        umf_ludel(Lup);
     end
 
     if(Areal & Breal)
@@ -723,22 +723,40 @@ function [res_d, res_v] = speigs(A, %_B, nev, which, maxiter, tol, ncv, cholB, r
         else
             sigmar = real(sigma);
             sigmai = imag(sigma);
-            [dr, di, z, resid, v, iparam, ipntr, workd, workl, info_eupd] = dneupd(rvec, howmny, _select, dr, di, z, sigmar, sigmai, workev, bmat, nA, which, nev, tol, resid, ncv, v, iparam, ipntr, workd, workl, info_eupd);
+            computevec = rvec;
+            if iparam(7) == 3 & sigmai then
+                computevec = 1;
+            end
+            [dr, di, z, resid, v, iparam, ipntr, workd, workl, info_eupd] = dneupd(computevec, howmny, _select, dr, di, z, sigmar, sigmai, workev, bmat, nA, which, nev, tol, resid, ncv, v, iparam, ipntr, workd, workl, info_eupd);
             if(info_eupd <> 0)
                 error(msprintf(gettext("%s: Error with %s: info = %d.\n"), "eigs", "DNEUPD", info_eupd));
             else
-                res_d = complex(dr,di);
-                res_d(nev+1) = [];
+                if iparam(7) == 3 & sigmai then
+                    res_d = complex(zeros(nev + 1,1));
+                    i = 1;
+                    while i <= nev
+                        if(~di(i))
+                            res_d(i) = complex(z(:,i)'*A*z(:,i), 0);
+                            i = i + 1;
+                        else
+                            real_part = z(:,i)' * A * z(:,i) + z(:,i+1)' * A * z(:,i+1);
+                            imag_part = z(:,i)' * A * z(:,i+1) - z(:,i+1)' * A * z(:,i)
+                            res_d(i) = complex(real_part, imag_part);
+                            res_d(i+1) = complex(real_part, -imag_part);
+                            i = i + 2;
+                        end
+                    end
+                else
+                    res_d = complex(dr, di);
+                end
+                res_d = res_d(1:nev);
                 if(rvec)
-                    res_d = diag(res_d)
+                    index = find(di~=0);
+                    index = index(1:2:$);
                     res_v = z;
-                    c1 = 1:2:nev + 1;
-                    c2 = 2:2:nev + 1;
-                    if(modulo(nev + 1, 2) == 1)
-                        c1($) = [];
-                    end
-                    res_v(:,[c1, c2]) = [res_v(:,c1) + res_v(:,c2) * %i res_v(:,c1) - res_v(:,c2) * %i];
-                    res_v(:,$) = [];
+                    res_v(:,[index index+1]) = [complex(res_v(:,index),res_v(:,index+1)), complex(res_v(:,index),-res_v(:,index+1))];
+                    res_d = diag(res_d);
+                    res_v = res_v(:,1:nev);
                 end
             end
         end
@@ -755,7 +773,7 @@ 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(rvec & iparam(7) == 1 & matB)
         if issparse(%_B)
             res_v = umf_lusolve(Rprimefact, res_v);
             if(~cholB)
@@ -1134,7 +1152,6 @@ function [res_d, res_v] = feigs(A_fun, nA, %_B, nev, which, maxiter, tol, ncv, c
                 error(msprintf(gettext("%s: Error with %s: info = %d.\n"), "eigs", "DSEUPD", info));
             else
                 res_d = d;
-
                 if(rvec)
                     res_d = diag(res_d);
                     res_v = z;
@@ -1143,12 +1160,33 @@ function [res_d, res_v] = feigs(A_fun, nA, %_B, nev, which, maxiter, tol, ncv, c
         else
             sigmar = real(sigma);
             sigmai = imag(sigma);
-            [dr, di, z, resid, v, iparam, ipntr, workd, workl, info_eupd] = dneupd(rvec, howmny, _select, dr, di, z, sigmar, sigmai, workev, bmat, nA, which, nev, tol, resid, ncv, v, iparam, ipntr, workd, workl, info);
+            computevec = rvec;
+            if iparam(7) == 3 & sigmai then
+                computevec = 1;
+            end
+            [dr, di, z, resid, v, iparam, ipntr, workd, workl, info_eupd] = dneupd(computevec, howmny, _select, dr, di, z, sigmar, sigmai, workev, bmat, nA, which, nev, tol, resid, ncv, v, iparam, ipntr, workd, workl, info);
             if(info <> 0)
                 error(msprintf(gettext("%s: Error with %s: info = %d.\n"), "eigs", "DNEUPD", info));
             else
-                res_d = complex(dr,di);
-                res_d(nev+1) = [];
+                if iparam(7) == 3 & sigmai then
+                    res_d = complex(zeros(nev + 1,1));
+                    i = 1;
+                    while i <= nev
+                        if(~di(i))
+                            res_d(i) = complex(z(:,i)'*A*z(:,i), 0);
+                            i = i + 1;
+                        else
+                            real_part = z(:,i)' * A * z(:,i) + z(:,i+1)' * A * z(:,i+1);
+                            imag_part = z(:,i)' * A * z(:,i+1) - z(:,i+1)' * A * z(:,i)
+                            res_d(i) = complex(real_part, imag_part);
+                            res_d(i+1) = complex(real_part, -imag_part);
+                            i = i + 2;
+                        end
+                    end
+                else
+                    res_d = complex(dr,di);
+                end
+                res_d = res_d(1:nev);
                 if(rvec)
                     res_d = diag(res_d)
                     res_v = z;

diff --git a/scilab/modules/arnoldi/tests/nonreg_tests/bug_12238.dia.ref b/scilab/modules/arnoldi/tests/nonreg_tests/bug_12238.dia.ref
new file mode 100644 (file)
index 0000000..cdd3de4
--- /dev/null
+++ b/scilab/modules/arnoldi/tests/nonreg_tests/bug_12238.dia.ref
@@ -0,0 +1,19 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - Scilab Enterprises - Sylvestre Ledru
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+// <-- CLI SHELL MODE -->
+// <-- Non-regression test for bug 12238 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=12238
+//
+// <-- Short Description -->
+//    [d v] = eigs(A) is broken for sparse matrices
+// =============================================================================
+A = sparse(rand(10,10));
+[d v] = eigs(A);
+val=norm(A*v-v*d);
+assert_checkalmostequal(val,0,0,30*%eps);

diff --git a/scilab/modules/arnoldi/tests/nonreg_tests/bug_12238.tst b/scilab/modules/arnoldi/tests/nonreg_tests/bug_12238.tst
new file mode 100644 (file)
index 0000000..1c2a5c2
--- /dev/null
+++ b/scilab/modules/arnoldi/tests/nonreg_tests/bug_12238.tst
@@ -0,0 +1,23 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - Scilab Enterprises - Sylvestre Ledru
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+
+// <-- CLI SHELL MODE -->
+
+// <-- Non-regression test for bug 12238 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=12238
+//
+// <-- Short Description -->
+//    [d v] = eigs(A) is broken for sparse matrices
+// =============================================================================
+
+A = sparse(rand(10,10));
+[d v] = eigs(A);
+val=norm(A*v-v*d);
+
+assert_checkalmostequal(val,0,0,30*%eps);