Wed 25 Nov 2015 12:58:31 PM UTC, comment #9:
The original comment does not specified which version of ARPACK is involved. If 3.1.5, dneupd is buggy.
Marco
|
Wed 25 Nov 2015 12:19:44 AM UTC, comment #8:
Can this bug be closed as not an Octave bug, but an upstream bug in ARPACK's implementation of dneupd?
|
Fri 13 Feb 2015 11:03:52 AM UTC, comment #7:
Dear all,
there is a bug in ARPACK up to 3.1.5 related to dneupd. It was fixed in 3.2.0 and I'm proposing a better fix here https://github.com/opencollab/arpack-ng/pull/2
It would be nice if you can check that it really fix your problems.
Marco
|
Thu 29 Jan 2015 05:30:38 AM UTC, comment #6:
I have just thought of an ugly work-around for this problem. Before I create the patch, could someone (Jordi?) comment on whether or not it would be suitable?
It seems that dneupd returns the eigenvalues in decreasing order when it works correctly. However, when it fails the results are not monotonic (since they are the largest eigenvalues). We could check if the results are monotonic, and if not, repeat the calculation.
Yes, I said it is ugly... Thoughts?
|
Sat 24 Jan 2015 05:41:13 AM UTC, comment #5:
Thanks, Jordi, for pointing out that the initialization is random. However, that is not causing the problem here.
It seems that the call to the dnaupd function in ARPACK works reliably (and that is where the randomness in the algorithm is).
The problem is in the call to dneupd. Most of the time, it correctly returns the first k values of workl[iptr(5)] as expected. However sometimes it returns the last two (which are not valid eigenvalues) instead of the first two.
Sometimes it replaces three or four, depending on the value of k.
Note that the main difference between
D = eigs(...)
and
[V,D] = eigs(...)
is in the first parameter passed to dneupd.
bug #41305 also suggests that there is a (different) bug in dneupd, so I wouldn't rule out that possibility.
|
Wed 21 Jan 2015 02:33:19 PM UTC, comment #4:
As you can read [here](http://scicomp.stackexchange.com/questions/18688/eigs-routine-in-octave) , see the comments, the error does not occur if eigs is called with
D=eigs(...)
instead of
[V,D]=eigs(...).
Therefore, to me it seems like an octave bug not a problem with ARPACK. But I am not developer...
|
Tue 20 Jan 2015 09:09:28 PM UTC, comment #3:
> The function eigs should not return different results each time it is run.
The Arnoldi method involves some initial randomness. Getting slightly different results each time is normal:
https://en.wikipedia.org/wiki/Arnoldi_iteration#Krylov_subspaces_and_the_power_iteration
|
Tue 20 Jan 2015 06:14:53 PM UTC, comment #2:
I modified the code to ensure convergence (flag -- fl) and to check if the eigenvalues are real. If eigenvalues may be complex the sort command can cause the incorrect results.
% solves the Eigenvalue Problem
% -y''= lambda^2 y on the interval (a,b)
%
% with mixed boundary conditions
% -- dirichlet boundary condition y(a)=0
% -- robin boundary condition y'(b)=y(b)
clear all;
close all;
p=6;
%% number of points
n=2^p;
%% grid space
h=1/2^p;
%% numbers of eigenvalues to find
kmax=8;
%% assemble symmetric 2nd finite difference matrix
Lh=-gallery('tridiag',n,1,-2,1)./(h^2);
%% apply robin boundary in last line
Lh(n,:)=0.0;
Lh(n,n-1)=-2/h^2;
Lh(n,n)=-1/h^2(2h-2);
%% options for eigs
% opts.tol=1e-10;
opts.p=max(floor(n/2),2*kmax);
comp=0;
while comp<10 % do the same thing 10 times
%% use eig to get all eigenvalues
[~,LambdaEig]=eig(Lh);
if isreal(diag(LambdaEig))
LambdaEig=sort(diag(LambdaEig));
lambdaEig=LambdaEig(1:kmax);
else
warning('eigenvalues of eig are not real')
end
%% use eigs routine to get 1st kmax eigenvalues of smallest magnitude
[~,LambdaEigs,fl]=eigs(Lh,kmax,'sr',opts);
if isreal(diag(LambdaEigs))
LambdaEig=sort(diag(LambdaEig));
lambdaEigs=sort(diag(LambdaEigs));
else
warning('eigenvalues of eigs are not real')
end
%% compute difference of 1st kmax eigenvalues of smallest magnitude
errmax=max(abs(lambdaEig-lambdaEigs));
disp(['eigs flag: ' num2str(fl) ', maximum error of eigenvalues: '...
num2str(errmax)])
comp=comp+1;
end
|
Tue 20 Jan 2015 02:15:15 AM UTC, comment #1:
The function eigs should not return different results each time it is run. It is clearly returning incorrect results. However, this function uses a third-party library; we'll have to see if the bug is in that library or Octave's own code.
|
Wed 14 Jan 2015 06:06:42 PM UTC, original submission:
The code below demonstrates a large difference in the computed eigenvalues of the eig- and eigs-routine. eigs will find complex pairs whereas eig finds real ones.
It gets much better if you replace 'sm' by 'sr'. Execute the script a few times and observe that somehow large difference occure. This seems unreliable to me.
By the way, if you choose 'sr' and p=3, you will get the error
|