Arpack puzzles: ask the maintainers?

General Usage
1

I use Arpack in the shift-invert mode for the generalized eigenvalue problem. I try to set the maximum number of iterations and/or the tolerance, but both parameters appear to be ignored. Does anyone know who actually actively maintains Arpack (or does someone know what is going on inside the library with these parameters)?

Nota bene: I am not asking about Arpack.jl . I looked at that code and it appears to just pass these parameters along to the actual Arpack library.

EDIT: I forgot to mention that I obviously tried to email arpack@caam.rice.edu. That email address is dead.

1 Like
2

AFAIK Arpack is considered abandonned:

and most distros have moved to

which appears to be actively maintained.

1 Like
3

The Arpack.jl package is also based on arpack-ng.

@PetrKryslUCSD, before filing an issue at their site, I suggest providing more details of your problem here, since your report seems incongruous with others’ experience.
I am somewhat familiar with the Fortran code of ARPACK, and can confirm that those parameters are actually used. Further, when using Arpack.jl, I find that an inadequate value of maxiter leads to an ARPACKException(1), while increasing tol can reduce the iteration count.

4

Thanks, that is very helpful. I totally missed the transition to arpack-ng.

5

FWIW, I came across this incidentally too, and I was super-surprised to see actively maintained F77 code on Github. Given that I consider, functions longer than, say, 30 lines code smell, it is a very sobering experience.

1 Like
6

Here is a MishWE: It is a free-vibration problem for the solid cylinder. Code is included below.

With tol set to the machine-precision tolerance:

 40.832225 seconds (16.47 k allocations: 1.946 GiB, 0.78% gc time)                        
Natural frequencies: [0.0, 0.00018759, 0.000189618, 0.000225355, 0.000231742, 0.000251527, 2548.64, 2548.64, 2561.31, 4102.08, 4706.11, 4706.11, 5126.58, 6865.78, 6865.78, 6952.23, 6959.68, 6974.56, 6974.56, 7042.6, 7048.24, 7611.43, 7699.75, 7706.07, 7709.31, 8003.19, 8011.69, 8954.67, 8954.67, 9157.43, 9157.43, 9364.48, 9364.48, 9410.67, 9418.92, 9690.55, 9904.97, 10061.3, 10284.8, 10421.0, 10421.0, 10439.3, 10447.4, 10466.4, 10466.4, 10737.7, 10737.7, 10769.4, 10769.4, 11232.3, 11232.3, 11566.2, 11570.2, 11667.7, 11862.7, 11866.8, 11963.5, 11963.5, 12005.3, 12005.3, 12308.1, 12308.1, 12585.3, 12659.7, 12659.7, 12690.2, 12690.2, 12885.7, 12937.8, 13462.4, 13462.4, 13465.8, 13470.2, 13524.4, 13525.6, 13803.8, 13803.8, 14269.2, 14272.8, 14283.8, 14304.1, 14310.9, 14315.5, 14384.6, 14384.6, 14525.2, 14525.2, 14634.6, 14634.6, 14665.0, 14676.4, 14752.6, 14762.0, 14989.2, 14989.2, 15141.5, 15288.9, 15291.7, 15423.2, 15423.2, 15506.5, 15778.2, 15807.1, 15819.9, 15851.2, 15945.8, 15947.0, 16038.1, 16042.5, 16125.9, 16125.9, 16153.2, 16153.2, 16227.3, 16256.1, 16257.4, 16387.8, 16387.8, 16416.8, 16416.8, 16649.4, 16963.8, 16963.8, 16966.0, 17089.9, 17116.7, 17131.6, 17141.1, 17158.1, 17321.4, 17321.4, 17723.3, 17836.3, 17836.3, 17837.1, 17837.1, 17850.6, 17856.7, 17858.3, 17871.9, 17919.2, 17919.2, 17923.8, 17923.8, 17941.3, 17941.3, 17995.3, 18028.7, 18037.8, 18151.1] [Hz]

With tol set to 10^10 times the machine-precision tolerance (ten orders of magnitude bigger tolerance!):

 41.223566 seconds (16.47 k allocations: 1.946 GiB, 0.78% gc time)                        
Natural frequencies: [0.0, 0.000187664, 0.000189626, 0.000225381, 0.000231761, 0.000251539, 2548.64, 2548.64, 2561.31, 4102.08, 4706.11, 4706.11, 5126.58, 6865.78, 6865.78, 6952.23, 6959.68, 6974.56, 6974.56, 7042.6, 7048.24, 7611.43, 7699.75, 7706.07, 7709.31, 8003.19, 8011.69, 8954.67, 8954.67, 9157.43, 9157.43, 9364.48, 9364.48, 9410.67, 9418.92, 9690.55, 9904.97, 10061.3, 10284.8, 10421.0, 10421.0, 10439.3, 10447.4, 10466.4, 10466.4, 10737.7, 10737.7, 10769.4, 10769.4, 11232.3, 11232.3, 11566.2, 11570.2, 11667.7, 11862.7, 11866.8, 11963.5, 11963.5, 12005.3, 12005.3, 12308.1, 12308.1, 12585.3, 12659.7, 12659.7, 12690.2, 12690.2, 12885.7, 12937.8, 13462.4, 13462.4, 13465.8, 13470.2, 13524.4, 13525.6, 13803.8, 13803.8, 14269.2, 14272.8, 14283.8, 14304.1, 14310.9, 14315.5, 14384.6, 14384.6, 14525.2, 14525.2, 14634.6, 14634.6, 14665.0, 14676.4, 14752.6, 14762.0, 14989.2, 14989.2, 15141.5, 15288.9, 15291.7, 15423.2, 15423.2, 15506.5, 15778.2, 15807.1, 15819.9, 15851.2, 15945.8, 15947.0, 16038.1, 16042.5, 16125.9, 16125.9, 16153.2, 16153.2, 16227.3, 16256.1, 16257.4, 16387.8, 16387.8, 16416.8, 16416.8, 16649.4, 16963.8, 16963.8, 16966.0, 17089.9, 17116.7, 17131.6, 17141.1, 17158.1, 17321.4, 17321.4, 17723.3, 17836.3, 17836.3, 17837.1, 17837.1, 17850.6, 17856.7, 17858.3, 17871.9, 17919.2, 17919.2, 17923.8, 17923.8, 17941.3, 17941.3, 17995.3, 18028.7, 18037.8, 18151.1] [Hz]

The code:

using FinEtools
using Arpack

E = 70000*phun("MPa")::FFlt;
nu = 0.33::FFlt;
rho = 2700*phun("KG/M^3")::FFlt;
radius = 0.5*phun("ft"); 
leng = 2*phun("ft"); 
omegashift = (2*pi*100) ^ 2; # to resolve rigid body modes
nr = 10
nL = 40
tol = 1e10 * eps(1.0)   
neigvs = 150

fens, fes = H8cylindern(radius, leng, nr, nL)

geom = NodalField(fens.xyz)
u = NodalField(zeros(size(fens.xyz,1),3))
applyebc!(u)
numberdofs!(u)

MR = DeforModelRed3D
material = MatDeforElastIso(MR, rho, E, nu, 0.0)
femm = FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(3,2)), material)
K = stiffness(femm, geom, u)
K .= 0.5 * (K .+ transpose(K))
femm = FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(3,3)), material)
M = mass(femm, geom, u)
M .= 0.5 * (M .+ transpose(M))

@time evals,evecs,nev,nconv = eigs(K + omegashift*M, M; nev=neigvs, which=:SM, tol = tol)
evals .-=  omegashift;
fs = real(sqrt.(complex(evals)))/(2*pi)

println("Natural frequencies: $fs [Hz]")
7

A sweep through tolerances from 10^0 * eps(1.0) to 10^20 * eps(1.0) (twenty orders of magnitude), yields this set of times of the eigenvalue problem solution with Arpack:

[41.7089, 42.2171, 43.1736, 43.1846, 42.8208, 43.3437, 43.0354, 43.3041, 43.468, 43.1582, 43.0319, 43.1361, 43.1094, 35.9275, 36.3748, 29.1554, 30.0258, 29.2839, 29.317, 43.2936, 29.4419]

So from around 10^12 * eps(1.0) the solution TIME is affected to some degree by setting the tolerance.

Concerning the ACCURACY: Until the tolerance of 10^14 * eps(1.0) is used, the natural frequencies are identical to six decimal digits. With the tolerance of 2*10^14 * eps(1.0) 28 natural frequencies change, with difference up to 4%.

8

IIUC this problem is not difficult in a spectral sense. Thus only a few “iterations” (outer loops of the IRAM scheme) are needed for full convergence. Using a huge tolerance can reduce the number of iterations, as you saw. Each “iteration” involves many solves with the factored sparse matrix. Unfortunately such a solve is currently a single-threaded operation, with (AFAICT) poor memory locality, so performance is sad.

I’d be glad to learn if Krylov-Schur methods can reduce the number of solves needed in shift-invert schemes, or if there are convenient alternatives to the inefficiencies of CHOLMOD. (Perhaps @juthohaegeman knows?)

P.S. I had to change your example since H8cylindern was not in the released FinEtools.

9

I think you are right about the difficulty of the spectral problem, and if I understand the code of the Fortran library correctly, it will only exit once the IRAM iteration completed. I think that was your point about why the relaxed tolerance doesn’t result in much reduced time?

My apologies, I did not remember that this was a new function, available only on master.

10

Is your generalized eigenvalue problem hermitian/symmetric with one of the two matrices positive definite (I assume so since you talk about natural frequencies and list real eigenvalues). Then there is another way of solving it, and you might want to look at geneigsolve in KrylovKit.jl.

2 Likes
11

Thanks, that is good to know about.