search

opencollab / arpack-ng

Collection of Fortran77 subroutines designed to solve large scale eigenvalue problems.
Other
284 stars 124 forks source link

not able to compute eigenvalues for generalized non-hermitian eigenvalue problem. #149

Open manavbhatia opened 6 years ago

manavbhatia commented 6 years ago

Hi,

I am having trouble computing eigenvalues with arpack for the following eigenvalue problem: 

   M1 x = lambda M0 x

The matrices M1 and M0 are 178x178 and are attached with this email. They are from a finite element discretization of a beam.

While LAPACK’s dggev is able to compute the eigenvalues for my matrices without issues, I am unable to get arpack to do this form me. I am using “dnaupd" and it is returning with "info=-8". I am also getting the same error when attempting to use "eigs" command in octave. 

I have struggled with this for a few days and am not sure how to proceed. I have made sure that the driver routine calling “dnaupd" is fine, and the same error from octave gives me some assurance on the reproducibility. 

I would greatly appreciate if you could provide some guidance on how to identify the issue here.

M0.txt M1.txt

caliarim commented 6 years ago

Which eigenvalues are you trying to compute? Which mode?

manavbhatia commented 6 years ago

I have tried both smallest_magnitude and largest_magnitude with no success. I have played around with number of modes from 2 to 20.

caliarim commented 6 years ago

I'm away from my computer. Is it possible that M0 is not SPD? Can you compute chol (M0) in octave?

manavbhatia commented 6 years ago

chol(m0) and chol(m1) returns without errors on octave. However, m0 is symmetric, while m1 is not.

caliarim commented 6 years ago

Ok, I will give a look when back from my holidays.

manavbhatia commented 6 years ago

Here is a related thread on octave concerning the same issue: http://octave.1599824.n4.nabble.com/Issues-with-eigs-td4689314.html#a4689320 .

While different people are able to converge to somewhat different eigenvalues on different releases of octave, I did not get convergence using octave 4.4.1 on Mac OS with arpack installed through macports. I am linking to the same arpack in my user code.

caliarim commented 6 years ago

Which version of arpack are you using?

manavbhatia commented 6 years ago

v 3.6.2 (installed from macports on Mac 10.13.6).

fghoussen commented 6 years ago

@manavbhatia : did you try with slepc ? slepc + lapack should be OK (as you tried directly with lapack). if slepc + arpack is OK, it means the way you use arpack may be buggy (would be likely eigs use arpack in matlab / octave)... Let's be honest, arpack has never been user-friendly and the doc can be misleading !...

caliarim commented 6 years ago

In the documentation of dnaup2 it is written that info = -8 (lapack error) should never happen. So, yes, either there is a problem with lapack libraries (and this affects octave, too) or there is a problem with 3.6.2 (since other users with octave do not see the problem, but maybe they use a lower version).

manavbhatia commented 6 years ago

@fghoussen Yes, I have successfully solved this with slepc using its Krylov Schur solver. I have not tried to use its arpack interface (need to recompile slepc for that). My experience with arpack is through my own interface and through octave.

manavbhatia commented 6 years ago

In order to overcome the user-friendly challenge, I would recommend putting together a wrapper to ARPACK in Eigen (C++). This would be a header-only file that someone can use to access the functionality in ARPACK. I have one, which is what I am using to interface to arpack from my C++ code. I will be happy to share, if you think it will be worthwhile.

caliarim commented 6 years ago

@manavbhatia: Am I right that you are using lapack <= 3.5.0?

manavbhatia commented 6 years ago

@caliarim Since this is on a Mac, the library is linking to the system lapack provided by the Accelerate framework: 

Dhcp-90-164:build manav$ otool -L /opt/local/lib/libarpack.dylib 
/opt/local/lib/libarpack.dylib:
    /opt/local/lib/libarpack.2.dylib (compatibility version 3.0.0, current version 3.0.0)
    /opt/local/lib/libvecLibFort.dylib (compatibility version 0.0.0, current version 0.0.0)
    /opt/local/lib/libgcc/libgfortran.4.dylib (compatibility version 5.0.0, current version 5.0.0)
    /usr/lib/libSystem.B.dylib (compatibility version 1.0.0, current version 1252.50.4)
    /opt/local/lib/libgcc/libgcc_s.1.dylib (compatibility version 1.0.0, current version 1.0.0)
    /opt/local/lib/libgcc/libquadmath.0.dylib (compatibility version 1.0.0, current version 1.0.0)
    /System/Library/Frameworks/Accelerate.framework/Versions/A/Accelerate (compatibility version 1.0.0, current version 4.0.0)
caliarim commented 6 years ago

So, 3.0.0? Anyway, I can reproduce the problem with lapack 3.5.0. The problem is in dlahqr.f. It was fixed in 3.6.0 http://www.netlib.org/lapack/errata_from_350_to_360.html#_strong_span_class_green_bug125_span_strong_xlahqr_v3_5_does_not_converge_on_some_matrices_while_xlahqr_v2_converges

manavbhatia commented 6 years ago

Interesting. So, I should install a more recent lapack and link to that instead. I will give this a shot and report back.

caliarim commented 6 years ago

Did you use a plain fortran implementation when you discovered the bug? If yes, you can download dlahqr.f from lapack 3.8.0 and add it to the compilation. I did

gfortran bug149.f dlahqr.f -larpack -llapack -lblas

where bug149.f is my own implementation of your problem.

manavbhatia commented 6 years ago

My application code is in C++, which uses Eigen. I link to arpack provided by macport (3.6.2), which links to lapack provided by Accelerate framework.

I did also try compiling 3.6.2 from source but linked that to the system lapack from Accelerate framework, with the same issue.

I could compile dlahqr with gfortran and then provide to my linker before -llapack, which would pick the dlahqr routine from this new file instead of lapack. But that would be a temporary solution.

Since I am now looking to build lapack from source, do you think it would matter if I install ATLAS or lapack? Not sure what version of lapack is provided by ATLAS.

caliarim commented 6 years ago

ATLAS provides only few lapack routines. And, unfortunately, I'm not sure ATLAS is still developed: the last email from its main author is from August 2017. I recently switched to openblas.

manavbhatia commented 6 years ago

I built arpack v 3.6.2 from source using lapack v 3.8.0 and forced octave to link to this new dyld. I am now able to get the correct values. 

octave:1> m0 = load('M0.txt');
octave:2> m1 = load('M1.txt');
octave:3> eigs(m1, m0)
ans =

        -0.00002 - 2672301.02581i
        -0.00002 + 2672301.02581i
        -0.00017 + 2639902.96203i
        -0.00017 - 2639902.96203i
         0.00006 - 2587669.31754i
         0.00006 + 2587669.31754i
manavbhatia commented 6 years ago

Following are the steps that I took to plug my arpack into the octave installed with macports.

arpack installed with octave (from macports) links to system lapack from Accelerate framework. This is a fork from an older version of lapack and has some bugs that were fixed in a v 3.6 of lapack.

So, we need to change the linking of octave from the macports arpack to one built from source and linked against new lapack:

  1. install octave (this will install all dependencies of octave including libVecFort and arpack). 
    sudo port install octave
  2. install lapack 
    sudo port install lapack
  3. uninstall arpack and vecLibFort. Answer yes when macports warns about breaking dependencies. 
    sudo port uninstall arpack vecLibFort
  4. clone arpack-ng. In CMakeLists.txt modify arpack PROPERTIES VERSION and SOVERSION from 2 to 3 (this is because octave complains about needing version 3)
  5. build and install arpack from source and link to new lapack in macports (not the system lapack)
  6. create a link /opt/local/lib/libvecFort.dylib pointing to /opt/local/lib/lapack/liblapack.dylib
  7. create a link /opt/local/lib/libarpack.2.dylib to the new libarpack.3.0.0.dylib installed in step 5

octave should now work with the new arpack that links to the new lapack.

manavbhatia commented 6 years ago

On a related note, is the macport build of arpack handled by someone in this group?

If so, I would recommend creating a variant to use lapack installed via macports. Currently the only two variants that exist are accelerate and atlas. Even with atlas, I suspect that the routines not provided by atlas will come from system lapack.

fghoussen commented 6 years ago

@manavbhatia : so all is OK now if I got it well ? What about adding version check for lapack (> 3.6.1 ?) in cmake/autotools ?

personal opinion: since the day I "discovered" arpack (!), I feel a simple utility is missing to run/test arpack on user matrix [matrix market format] just to see/check eigen val/vec (provided fortran examples are hard-coded => new users can't "really" use arpack => need matlab or digging the doc/examples during weeks to get dummy code to do simple things...). I wrote such a utility at the time: I relied on PETSc (av, mv, inversion when necessary, ...): I could take an hour to open source this (as EXAMPLE/CUSTOM for instance ?). The major "drawback" would be that it needs petsc, but, this could be packaged in such a way that petsc is involved only if some kind of "EXAMPLE_CUSTOM" option is set. Not sure how this could sound to you (good/bad)...

manavbhatia commented 6 years ago

I support the addition of such a utility. I think it will be great, so that users don't have to write and debug their own reverse-communication code.

At my end, I have put together a wrapper in Eigen/C++ which links to the fortran version of arpack (I add the "_" after the function name). This could be added as a contribution to both arpack and eigen. Currently I can only find this: https://eigen.tuxfamily.org/dox/unsupported/group__ArpackSupport__Module.html

I have put together my wrappers for both lapack dggev and arpack dnaupd. I am happy to make my wrapper available. It will need some modifications to add the template parameters so that it can be used with arbitrary types of matrix/vector and linear solver types.

What you recommend will be great if people need to read in their matrices and compute eigenvalues. But typically eigenvalue computations need to fit within a context of a large computation. So, wrappers in octave/Eigen would be greatly useful. Slepc also provides a wrapper to arpack.

fghoussen commented 6 years ago

I support the addition of such a utility. I think it will be great, so that users don't have to write and debug their own reverse-communication code.

:+1: ... But I guess everybody (including me !) will agree on the fact that adding a dependency on petsc is not great but mandatory. I had a quick look at that: my old utility is still running fine (would just need to be adapted a bit to be open-sourced - no such a big deal). But the "problem" is more about the dependency:

... Not sure what would be the best. For me the idea is good, but, the way to do it in a satisfactory way is no clear. Advices and / or points of view of you guys are welcome !

manavbhatia commented 6 years ago

petsc vs mumps: I would vote for petsc, which can leverage mumps as a third party solver through command-line options. Also, on some systems (I just checked ubuntu) mumps is a default dependency of petsc. If you are compiling petsc from source, then it can download and install mumps for you automatically.

petsc vs eigen: I would vote for both, since they serve two different user bases. Eigen is a good matlab replacement in the c++ world with a very natural API for matrix/vector manipulation. All my quick and dirty code happens here (and I would assume for others too). It also has support for sparse matrices and solvers. What I have missed so far is the ability to compute eigenvalues (its native eigensolvers don't always do the job), hence my wrapper class for arpack. Petsc is awesome, but uses a more complex API, especially since it was built with MPI and HPC in mind.

So, I would vote for both eigen and petsc.

Also, note that slepc, which builds on petsc, has a wrapper for arpack. So, you may not have to really build one for petsc, assuming the users can install and use slepc instead of petsc. Downside is that there is one more dependency.

caliarim commented 6 years ago

I feel a simple utility is missing to run/test arpack on user matrix [matrix market format] just to see/check eigen val/vec

I think we should distinguish run and see from test and check. I ran arpack during my phd thesis, in plain fortran. Now I use it from time to time through octave. I do not know what it is necessary today to "work" with arpack. But, in order to test/check it, I think arpack+lapack is ok. All the bug reports I worked on, were nice enough to use small matrices. Most of times the matrices were even in full format (like in this report). I used the drivers in EXAMPLES and modified them. In fact, they are hardcoded with special matrices (usually tridiagonal). In my opinion, in order to test/check arpack, it would be enough to generalise those drivers in order to read user matrices (full or coordinate format). And use the proper lapack routines for general matrices. This is the opinion of an old guy which feels confortable with fortran77...

manavbhatia commented 6 years ago

Just a F77 utility to work with user matrices, like that suggested by @caliarim, would be very helpful. This is the minimum that could be met without a need for other dependencies.

This would have helped in debugging my arpack wrapper. I found myself looking through the source-code for octave and slepc to figure out what I might be doing wrong.

fghoussen commented 6 years ago

So, I would vote for both eigen and petsc.

petsc: cmake file is not provided, pc file is provided (if configured) but not shipped by .deb package => impossible to integrate this in arpack autotools/cmake... (+ CI KO !)

eigen: may be a solution (provide both cmake and pc files + shipped by .deb package) that make sense.

I think arpack+lapack is ok.

Yes, but, for "small enough" matrices only... Not for large sparse matrices I had to deal with (more and more of them with, for instance, the "Big Data" oriented stuffs and friends)

Just a F77 utility to work with user matrices, like that suggested by @caliarim, would be very helpful.

No way ! :D Reading ASCII file is a nigthmare in fortran (I guess that's why examples provided are hard coded !)

it would be enough to generalise those drivers in order to read user matrices

Not sure to understand. Huge job. Many examples. Different operations to apply depending on mode, ido value... No simple way to do that (?)

should distinguish run and see from test and check.

The utility I designed (based on petsc) at the time was more meant to "play" (change option: shift, invert, shift+invert, large/small mag, ...) => check/compare ev, iterations, time to solution => it's very instructive play (time to solution) depending on matrices (not sure matlab/octave can really enable that as you use them as some kind of black boxes... And don't know what options/modes have been used and why/how - actually, I was using slepc and had huge slow down... So I get back to arpack to understand how it works... and how to set up slepc !)

fghoussen commented 6 years ago

If I find some free time, I can try to prototype something...

caliarim commented 6 years ago

No way ! :D Reading ASCII file is a nigthmare in fortran (I guess that's why examples provided are hard > coded !)

open(4,file='M0.txt')
do 5 i = 1,n
read(4,*) (m0(i,j),j=1,n)
5    continue
close(4)
fghoussen commented 6 years ago

@caliarim : you forgot to handle comments %, empty lines, hints (coordinates, matrix, symmetric, ..) in the first line, handling optional parameters (nnz)... Fortran is not the "good" solution for reading txt files ! :D Try to write an xml file someday...

caliarim commented 6 years ago

@fghoussen: we are speaking of different things. I am speaking of debugging. You see a problem using arpack from {octave,eigen,petsc,slepc,whatever_you_want}. You restrict the problem to a small example. Now you want to understand whether the problem is in arpack itself or in the previous set of interfaces. So, you need to run plain F77 arpack onto your small example.

fghoussen commented 6 years ago

@caliarim : oh yeah, I see ! I was more thinking of a utility with switches you could play with (no debug). Quite busy these days, but, I'll try to prototype something (based on what I did) if I can save some time...

fghoussen commented 6 years ago

Found bugs trying to (begin to !) backport my old utility in current code base...

I've just PR'ed a few commits. If more experienced guys than me (@caliarim ?) could double check #153, #154, #155 would be a good thing...

Also I need #155 to be merged if I want to have a chance to get this done !

fghoussen commented 6 years ago

@sylvestre : BTW, would you be OK to add a (conditional) dependency on eigen3 (if configure --enable-mm for instance) to get such kind of utility ? If yes, I can try to finish this next week.

sylvestre commented 6 years ago

why not but what for?

fghoussen commented 6 years ago

This is where I had to stop for now: check #157 . Any help would be really appreciated ! :D Do NOT merge #157 !

fghoussen commented 5 years ago

OK, at that point the arpackmm utility enables to use/try all build-in eigen solvers. Could be possible to add support for external solvers like superlu or umfpack but this would make packaging a lot more difficult (and add arpack-ng dependencies on these packages): I won't go into this. Also could be possible to extract some kind of arpackSolver.hpp from arpackmm, but, as eigen provides already it's own eigen solvers, not sure this is worth.