From 684f27a5101ca8b007d410db77f33cff9fff1ce4 Mon Sep 17 00:00:00 2001 From: Guillaume Horel Date: Fri, 25 Jan 2013 14:11:32 +0100 Subject: [PATCH] * 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 --- scilab/CHANGES_5.4.X | 2 + scilab/modules/arnoldi/macros/eigs.sci | 80 +++++++++++++++----- .../arnoldi/tests/nonreg_tests/bug_12238.dia.ref | 19 +++++ .../arnoldi/tests/nonreg_tests/bug_12238.tst | 23 ++++++ 4 files changed, 103 insertions(+), 21 deletions(-) create mode 100644 scilab/modules/arnoldi/tests/nonreg_tests/bug_12238.dia.ref create mode 100644 scilab/modules/arnoldi/tests/nonreg_tests/bug_12238.tst diff --git a/scilab/CHANGES_5.4.X b/scilab/CHANGES_5.4.X index b8b0f06..c57b2ff 100644 --- 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 --- 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 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 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); -- 1.7.9.5