eigs routine in octave [closed]

I am using octave and observed a problem with the eigs-routine for non symmetric matrices. Using GNU octave version 3.8.1 the code below gives significant difference of eigenvalues although same commands are executed. The difference of eigenvalues is computed by subracting the eig- and eigs-eigenvalues.

% 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*(2*h-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');
if isreal(diag(LambdaEigs))
LambdaEig=sort(diag(LambdaEig));
lambdaEigs=sort(diag(LambdaEigs));
else
warning('eigenvalues of eigs are not real')
end
errmax=max(abs(lambdaEig-lambdaEigs));
disp(['eigs flag: ' num2str(fl) ', maximum error of eigenvalues: '...
num2str(errmax)])
comp=comp+1;
end

Do you have an explanation for this? GNU octave's manual reports that eigs is based on the ARPACK package. I guess ARPACK is widely used for large eigenvalue problems and tested well. But GNU octave's eigs-routine seems not to be reliable for this problem.

I would like to use the eigs routine for eigenvalue computations of larger general eigenvalue problems $A x=\lambda B x$ with $A\in\mathbb{R}^{2000\times2000}$ in GNU octave. Because of the matrix size Krylow subspace methods are reasonably to me and as far as I know eigs is a modified Arnoldi iteration.

Do you have any recommendations how to continue?

closed as off-topic by Brian Borchers, Christian Clason, Kirill, nicoguaro♦Apr 6 '18 at 23:48

This question appears to be off-topic. The users who voted to close gave this specific reason:

If this question can be reworded to fit the rules in the help center, please edit the question.

• What is the actual output you are getting? Possible duplicate of scicomp.stackexchange.com/questions/8073/… – Kirill Jan 20 '15 at 7:34
• Here is the output fl = 0 errmax = 1.5815e+04 fl = 0 errmax = 1.3520e-12 fl = 0 errmax = 2.2737e-12 fl = 0 errmax = 2.7285e-12 fl = 0 errmax = 1.8190e-12 fl = 0 errmax = 1.7906e-12 fl = 0 errmax = 1.5818e+04 fl = 0 errmax = 1.5805e+04 fl = 0 errmax = 2.2737e-12 fl = 0 errmax = 1.3074e-12 fl is the flag of eigs indicating convergence.I also saw the reference you gave. There it says that specifying the number of Lanczos vectors in opts is a solution. But it did not work out. – sebastian_g Jan 20 '15 at 8:29
• I haven't run the code, but one issue could be that you are calling sort on the eigenvalues, which might be complex in general since your boundary conditions create a nonsymmetric matrix. For complex numbers, sort defaults to ordering by nondecreasing modulus: there could be mismatches if there are eigenvalues with the same modulus (e.g., complex conjugate pairs). – Federico Poloni Jan 20 '15 at 15:28
• I experimented a bit with your code in octave 3.8.1. The fact that the results vary so widely from trial to trial is very disturbing. I modified your code slightly to print all 8 eigenvalues from both eig and eigs. Sometimes only one or two differ; sometimes all 8 are very different. I would encourage you to submit a bug report to the Octave team. (The MATLAB version of eigs doesn't exhibit this problem.) – Bill Greene Jan 20 '15 at 19:44
• @BillGreene I already submitted a bug report to the developers. @Kirill It is interesting that bad answers occur only if eigs is called with [V,D]=eigs(...). I will add this to my bug report. Thanks a lot. – sebastian_g Jan 21 '15 at 8:45