bugGNU Octave - Bugs: bug #31479, Crash & bugs in eigs

 
 

bug #31479: Crash & bugs in eigs

Submitted by:  Otto Seiskari <oseiskar>
Submitted on:  Wed 27 Oct 2010 04:29:08 PM UTC  
 
Category: NoneSeverity: 3 - Normal
Priority: 5 - NormalItem Group: Segfault, Bus Error, etc.
Status: FixedAssigned to: None
Originator Name: Open/Closed: Closed
Release: 3.2.4Operating System: GNU/Linux

Add a New Comment(Rich Markup)
   

You are not logged in

Please log in, so followups can be emailed to you.

 

(Jump to the original submission Jump to the original submission)

Sat 29 Jan 2011 08:12:38 PM UTC, comment #18:

I went and benchmarked eigs for a sparse 1000x1000 matrix. Everything looks fine.

The newly embedded version uses my CFLAGS and FFLAGS

The default Debian package compiled with -O3. I ran 100 tests on the same matrix.

External ARPACK
1.1105 average run time in seconds

Internal ARPACK
1.1204 average run time in seconds

Rik <rik5>
Project Administrator
Sat 29 Jan 2011 07:39:29 PM UTC, comment #17:

Yes, this bug has been addressed and I will close the report.

The dneupd -8 error is linked to the initial guess for the starting solution and to the number of Lanzcos basis vectors used. By supplying a static starting point for eigs through the OPTS structure I can't get it to fail with several 100K executions. Similarly, If I adjust the number of basis vectors using OPTS.p to 3 I can't get it to fail. The basis vector issue won't generally be a problem on a truly sparse matrix where the dimension N is much greater than the number of eigenvalues requested.

Rik <rik5>
Project Administrator
Sat 29 Jan 2011 05:22:34 PM UTC, comment #16:

I checked again, and running

In my copy of Octave 3.3.90, which is built with ARPACK from the Debian package, I see a crash after a few calls:

With a copy of Octave built from the current sources (ARPACK sources moved to Octave), I don't see the crash, but I do sometimes (every few thousand calls?) see

Perhaps this is just because the algorithm is iterative and picks a starting point using some random numbers?

I don't see either of the above errors in 10,000 calls if I set

so this seems somehow dependent on the small size of the test matrix in the original bug report, which Rik has noted is not the type of matrix that eigs was really designed for.

So it would be great if someone wants to try to debug the "error -8" problem, but I don't think we should worry about this for the release. The more important thing for me is that the crash is gone and I'm no longer seeing complaints from valgrind about uninitialized values. Given that, I'd like to close this report. OK?

John W. Eaton <jwe>
Project Administrator
Sat 29 Jan 2011 03:16:36 AM UTC, comment #15:

I started with the sources from the Debian package, but now I see that I missed the patches from the debian/patches directory.

Yes, it looks like David's patch is there.

But I thought I was still seeing the the problem with the Debian ARPACK package.

In any case, I would still like to fix the problem for all users of Octave, and given that there is little hope for getting the patch applied in all distributions or upstream, I think it still makes sense to include the ARPACK sources with Octave.

John W. Eaton <jwe>
Project Administrator
Fri 28 Jan 2011 11:24:23 PM UTC, comment #14:

I think that ARPACK in Debian is built with David's patch. I'm not sure because the links I gave in the Debian bug report do not work anymore (which is strange, because I'm really surce I copy/pasted them back then).

The bug report against ARPACK in Debian:
http://bugs.debian.org/572935

Thomas Weber <tweber>
Fri 28 Jan 2011 07:34:51 PM UTC, comment #13:

I just checked with Ubuntu's, i.e., Debian's packaging of ARPACK. As expected, they use slightly different CFLAGS and FFLAGS than is perhaps normal. From the packaging files:

The license file included with the Debian distribution is copyright 2001, while the one from the ARPACK website is copyright 2008 and slightly shorter. Thus, I'd propose going with the shorter, newer version.

Rik <rik5>
Project Administrator
Fri 28 Jan 2011 07:33:55 PM UTC, comment #12:

I've checked in the changes. Unfortunately, just after checking them in, I noticed that I'm sometimes getting this:

This is quite odd, as I thought I was using the same sources for ARPACK as the debian package, which is what I was using before. I'll try to debug later.

The comments in dneupd say that error -8 is

Hmm...

John W. Eaton <jwe>
Project Administrator
Fri 28 Jan 2011 07:28:35 PM UTC, comment #11:

I also started a discussion about this topic on the maintainers list, which is archived starting here:

https://mailman.cae.wisc.edu/pipermail/octave-maintainers/2011-January/022538.html

Rik, in the debian package, it looks like they compile with -O3 but no other special options. I doubt there would be too much of a difference between that and the default -O2 when using gfortran, but I haven't tried to do any testing. ARPACK already uses the BLAS and LAPACK for linear algebra, so what is used for that will probably have more of an effect than the compiler options for ARPACK itself.

In any case, the response so far has been to include ARPACK in Octave, so I'll go ahead with that.

Note that the patch from David does not solve the problem of sometimes getting different results. But as Rik has pointed out, "The ARPACK library is an iterative method and it's results are absolutely going to depend on where in the solution space the search begins".

Marco, do you always get the exact same value when you use eigs(A,1), or is it sometimes a complex value? If it is sometimes complex, then I think you really should be getting a complex pair.

John W. Eaton <jwe>
Project Administrator
Fri 28 Jan 2011 06:44:18 PM UTC, comment #10:

Right, I forgot about the advantages of Open Source; it hadn't occurred to me to just grab the code and fix it. From a process perspective, I think this is the best option as I see little hope of getting the original developers to update their code. The license for ARPACK is at http://www.caam.rice.edu/software/ARPACK/. Does the top, BSD Software License, incorporate terms by reference other than the one's listed on this page? If not, the license isn't onerous and it should go in before the 3.4 release.

After it goes in, we might want to take a look at a few of the distro-provided ARPACK packages to see if the package maintainers are throwing in any special flags for the compilation. It might also be worth benchmarking the current situation. It would be a shame to introduce an Octave-compiled version of the library which results in lower performance than a distribution installed version.

Rik <rik5>
Project Administrator
Thu 27 Jan 2011 08:21:46 AM UTC, comment #9:

I have a series of patches ready to commit that will

  • add ARPACK to libcruft
  • fix a declaration problem with ETIME in an ARPACK function
  • apply David's patch
  • fix an uninitialized memory problem I discovered with valgrind while doing some testing (see the attached changeset)

Does anyone object to making these changes?

If not, then I will check these in, and also add a note to the NEWS file and update install.txi and anywhere else that arpack is mentioned.

(file #22523)

John W. Eaton <jwe>
Project Administrator
Thu 27 Jan 2011 06:09:41 AM UTC, comment #8:

I think we should do something about this problem before releasing 3.4.

Since ARPACK does not seem to be maintained now, and getting everyone who packages ARPACK to fix this bug for us is probably not going to happen either, maybe we should just put ARPACK into the Octave sources, same as we do with some other numerical routines like lsode, dassl, daspk, etc.

Then we can just patch it with David's fix that he posted more than 4 years ago here:

https://mailman.cae.wisc.edu/pipermail/octave-maintainers/2006-November/005005.html

David and Rik, please let me know whether you object to adding ARPACK to Octave. If you think it is OK, I'll make the changes.

John W. Eaton <jwe>
Project Administrator
Fri 03 Dec 2010 10:42:04 PM UTC, comment #7:

I had no idea that there was such a long history to this problem; this gives me no confidence that the problem will be fixed in ARPACK since they have already ignored it for 4 years.

You and I did come independently to the same conclusion as to the source of the problem and that is good. As part of this most recent investigation, jwe extracted a self-contained example which is attached to this bug report (.cc files, .f files, and Makefile). Would this be enough evidence to convince either the ARPACK maintainers or the Linux distro packagers?

My opinion now is that we should apply for getting ARPACK fixed, but separately solve the issue in Octave. It seems unnecessary to force users to wait another 4 years for this issue to be resolved, when we can suppress it ourselves. I went back and using the problem matrix from the first bug report in 2006 proved that my solution #3 is insufficient. We really will need to scale all the matrices to the size of p. I re-ran a test of solution #2 and it works for both the 2006 and 2010 problem matrices. According to your earlier work, this change would only need to be made for the NonSymmetric cases in eigs-base.cc, correct?

Rik <rik5>
Project Administrator
Fri 03 Dec 2010 09:27:13 PM UTC, comment #6:

My original bug report is here

http://www-old.cae.wisc.edu/pipermail/octave-maintainers/2006-November/001230.html

D.

David Bateman <dbateman>
Project Member
Fri 03 Dec 2010 09:20:42 PM UTC, comment #5:

I hope you read my mails

http://octave.1599824.n4.nabble.com/eigs-crash-on-a-particular-spare-matrix-td1673964.html#a1673970

and numerous others about this problem which has dated back to 2006? Making the change in ARPACK seems impossible (though maybe in PARPACK it might be possible). Getting a distribution to patch ARPACK themselves seems to best route to me, and for this we need a self-contained ARPACK example demonstrating the problem something I tried to do a few times, but never got around to finalizing. I'd prefer this solution to patching Octave to workaround a bug in ARPACK

D.

David Bateman <dbateman>
Project Member
Fri 03 Dec 2010 06:32:58 PM UTC, comment #4:

I spent some time over the Thanksgiving holiday looking into this issue. I have 3 possible solutions, but each one has some disadvantages.

The problem is in the Fortran code for dneupd.f. The C++ code passes in three arrays which are sized based on the number of requested eigenvalues (k) + 1. The extra 1 accounts for the fact that sometimes a complex conjugate pair will be found and so precisely one extra eigenvalue will be returned.

Within the dneupd code the variable in question is nconv, which is the number of converged Ritz values. In general, it is equal to k, but not always. The call which modifies it is a call to dtrsen. The documentation for dtrsen says that the modified value for k is in the range 0 <= k <= p (the number of Lanzcos basis vectors). For the example in question, p = 3, and occasionally it happens that nconv reaches this value whereas k + 1 = 2 and a segmentation violation results.

It appears, from my reading of the dneupd code, that the authors may not have realized that nconv was both an input and output from dtrsen. By replacing nconv with numcnv (an existing temporary variable containing the number of converged Ritz values) the program no longer segfaults. This is a one-line fix, but there are several problems with it. First, it is a fix in a code base we don't control and getting someone to look at this issue, understand it, and make the change may take quite some time. Second, although I reviewed the code and it produces the correct eigenvalues, I can't say that this change wouldn't have other effects on the algorithm. (I really don't think it does, but who knows?).

If solution 1, fixing ARPACK, is not the right course then we can handle it on the Octave-side by re-dimensioning everything to be large enough to avoid memory overwrites. The three arrays to consider are dr, di, and z. These are all double (8 byte) arrays. dr and di are of size (k + 1) while z is n * (k + 1). These would now have to be sized based on p which defaults in Octave code to 2 * k. The extra memory requirements are (k - 1) for both dr and di, but n * (k-1) * 8 bytes for z. When used with true sparse matrices n (# of rows in matrix) might be very large indeed (A brief look at the ARPACK web page showed examples where n was 90,000 and 2.2e6). So, while this solution is code correct, it suffers from excess memory usage. To be fair, if you're running the extremely large simulations your machine or cluster probably has a lot of memory. Even the large simulation I cited would need only an additional 176 MB of memory when calculating 10 eigenvalues.

Solution 3 is a modification of the previous solution, re-dimension the arrays, but only to size k + 2 rather than to size p. I have not found a situation where nconv exceeds k + 2. Admittedly, I have tried only a few matrices. If I understood the entire ARPACK algorithm better one might be able to assert that nconv will never exceed k + 2. It seems possible that there is another +1 limit on the change in k. So, the first call to dnaupd increases k to k+1. The second call in the code to dneupd then increases k+1 to k+2, but not more than that.

Any votes or words of advice for which solution to implement?

Rik <rik5>
Project Administrator
Tue 02 Nov 2010 10:33:24 AM UTC, comment #3:

Just my 2 cents: the call eigs(A,1) at least returns always a single eigenvalue (I tried 500 times). So, the problem with different number of eigenvalues is in the extraction procedure for the eigenvectors.

Marco

Marco Caliari <caliari>
Project Member
Fri 29 Oct 2010 01:44:02 AM UTC, comment #2:

The ARPACK library is an iterative method and it's results are absolutely going to depend on where in the solution space the search begins. There may be an opportunity to initialize the starting point based on some characteristic of the initial matrix which yields better convergence, but I'm not sure what that would be.

As an aside, you want to use the right tool for the job and eigs is not the right tool for a small matrix. You can just use eig and it will be more computationally efficient and deliver all of the eigenvalues, not just one as you were requesting.

Finally, I can confirm the a memory leak and crash in the eigs code. Attached is a small script which runs eigs repeatedly and always causes a crash on my system.

According to the backtrace from run-octave, the problem is in the destructor for the Array class.

(file #21855)

Rik <rik5>
Project Administrator
Thu 28 Oct 2010 02:39:56 AM UTC, comment #1:

I'm copying David Bateman on this bug report because he wrote the eigs function and the interface to arpack for Octave. Perhaps he will have some insight.

I can reproduce this problem on a Debian system with the current development sources for Octave.

I don't have a complete fix for the problem yet, but maybe the following info will help someone else to solve the problems.

I'm attaching a file, arpack-bug.tar.gz that is a pared-down version of the function that Octave uses to do this particular eigs call. In the arpack-bug.cc file there are a couple of comments about array dimensions that are different in Octave from what they are in arpack-bug.cc. When I use the values from Octave, I see warnings from valgrind (see the valgrind-output-before-changes file in the tarball) and those warnings disappear after making the changes.

I don't know what causes arpack to sometimes find an extra eigenvalue. I thought it might be the random value of resid that is used in the initial call to dnaupd, but even when that value is fixed as in my arpack-bug.cc file, it still sometimes finds a complex pair instead of a single real value. And valgrind is not telling me anything about uninitialized values.

OK, I see now that if info is zero on the first call to dnaupd, it will use its own random number generator to initialize the resid vector. So that explains why I was getting apparently random results when I thought I was using the same input values on each call. And I can confirm that by setting info to 1 on the first call instead of 0, and setting the same initial values for resid on each initial, arpack will return the same values each time. Unfortunately, it still may converge to a complex conjugate pair instead of a single real value, depending on the values selected for resid. So I still don't see how to fix the problem and ensure that arpack will converge to the same correct point each time.

(file #21848)

John W. Eaton <jwe>
Project Administrator
Wed 27 Oct 2010 04:29:08 PM UTC, original submission:

I am trying to calculate an eigenpair of a small matrix (and check convergence)
The result is (randomly) one of:

  1. The correct answer
  2. An additional eigenpair I did not ask for
  3. Crash

The following problem is probably also related

Otto Seiskari <oseiskar>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach File(s):
   
   
Comment:
   

Attached Files
file #22523:  diffs.txt added by jwe (2KiB - text/plain)
file #21855:  force_eigs_crash.m added by rik5 (99B - application/octet-stream)
file #21848:  arpack-bug.tar.gz added by jwe (52KiB - application/x-gzip)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -unavailable- added by tweber (Posted a comment)
  • -unavailable- added by dbateman (Posted a comment)
  • -unavailable- added by caliari (Posted a comment)
  • -unavailable- added by rik5 (Updated the item)
  • -unavailable- added by jwe (Updated the item)
  • -unavailable- added by jwe
  • -unavailable- added by oseiskar (Submitted the item)
  • -unavailable- added by oseiskar
  •  

    Do you think this task is very important?
    If so, you can click here to add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.

     

    Please enter the title of George Orwell's famous dystopian book (it's a date):

     

     

    Follow 8 latest changes.

    Date Changed By Updated Field Previous Value => Replaced By
    Sat 29 Jan 2011 07:39:29 PM UTCrik5StatusConfirmed=>Fixed
      Open/ClosedOpen=>Closed
    Thu 27 Jan 2011 08:21:46 AM UTCjweAttached File-=>Added diffs.txt, #22523
    Fri 29 Oct 2010 01:44:02 AM UTCrik5Attached File-=>Added force_eigs_crash.m, #21855
    Thu 28 Oct 2010 02:39:56 AM UTCjweAttached File-=>Added arpack-bug.tar.gz, #21848
      StatusNone=>Confirmed
      Carbon-Copy-=>Added david bateman
    Wed 27 Oct 2010 04:29:08 PM UTCoseiskarCarbon-Copy-=>Added -unavailable-

    Back to the top


    Powered by Savane 3.1-cleanup1