# Arpack puzzles: ask the maintainers?

**URL:** <https://discourse.julialang.org/t/arpack-puzzles-ask-the-maintainers/23082>\
**Category:** General Usage\
**Created:** [April 12, 2019, 2:49pm UTC](https://discourse.julialang.org/t/arpack-puzzles-ask-the-maintainers/23082 "2019-04-12T14:49:29Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [April 12, 2019, 2:49pm UTC](https://discourse.julialang.org/t/arpack-puzzles-ask-the-maintainers/23082/1 "2019-04-12T14:49:29Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [April 13, 2019, 5:48am UTC](https://discourse.julialang.org/t/arpack-puzzles-ask-the-maintainers/23082/2 "2019-04-13T05:48:43Z")

</div>

> [@PetrKryslUCSD](#):
>
> I obviously tried to email arpack@caam.rice.edu. That email address is dead.

AFAIK Arpack is considered abandonned:

> **[FS#28125 : \[arpack\] Switch upstream to arpack-ng](https://bugs.archlinux.org/task/28125)**
>
> Flyspray, a Bug Tracking System written in PHP.

and most distros have moved to

> **[GitHub - opencollab/arpack-ng: Collection of Fortran77 subroutines designed to...](https://github.com/opencollab/arpack-ng)**
>
> Collection of Fortran77 subroutines designed to solve large scale eigenvalue problems.

which appears to be actively maintained.

---

<div class="post-metadata">

**Author:** ![Ralph\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ralph_smith/32/10344_2.png) [@Ralph\_Smith](https://discourse.julialang.org/u/Ralph_Smith)\
**Post date:** [April 13, 2019, 2:21pm UTC](https://discourse.julialang.org/t/arpack-puzzles-ask-the-maintainers/23082/3 "2019-04-13T14:21:28Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [April 13, 2019, 2:37pm UTC](https://discourse.julialang.org/t/arpack-puzzles-ask-the-maintainers/23082/4 "2019-04-13T14:37:36Z")

</div>

Thanks, that is very helpful. I totally missed the transition to [arpack-ng](https://github.com/opencollab/arpack-ng).

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [April 13, 2019, 2:50pm UTC](https://discourse.julialang.org/t/arpack-puzzles-ask-the-maintainers/23082/5 "2019-04-13T14:50:17Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [April 13, 2019, 3:31pm UTC](https://discourse.julialang.org/t/arpack-puzzles-ask-the-maintainers/23082/6 "2019-04-13T15:31:14Z")

</div>

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:

```julia
 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!):

```julia
 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:

```julia
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]")

```

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [April 13, 2019, 4:53pm UTC](https://discourse.julialang.org/t/arpack-puzzles-ask-the-maintainers/23082/7 "2019-04-13T16:53:33Z")

</div>

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:

```julia
[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%.

---

<div class="post-metadata">

**Author:** ![Ralph\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ralph_smith/32/10344_2.png) [@Ralph\_Smith](https://discourse.julialang.org/u/Ralph_Smith)\
**Post date:** [April 13, 2019, 8:38pm UTC](https://discourse.julialang.org/t/arpack-puzzles-ask-the-maintainers/23082/8 "2019-04-13T20:38:12Z")

</div>

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`.

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [April 13, 2019, 8:52pm UTC](https://discourse.julialang.org/t/arpack-puzzles-ask-the-maintainers/23082/9 "2019-04-13T20:52:32Z")

</div>

> [@Ralph\_Smith](#):
>
> IIUC this problem is not difficult in a spectral sense. Thus only a few “iterations” (outer loops of the IRAM

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?

> [@Ralph\_Smith](#):
>
> P.S. I had to change your example since `H8cylindern` was not in the released `FinEtools` .

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

---

<div class="post-metadata">

**Author:** ![juthohaegeman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juthohaegeman/32/8620_2.png) [@juthohaegeman](https://discourse.julialang.org/u/juthohaegeman)\
**Post date:** [April 14, 2019, 11:02am UTC](https://discourse.julialang.org/t/arpack-puzzles-ask-the-maintainers/23082/10 "2019-04-14T11:02:14Z")

</div>

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](https://github.com/Jutho/KrylovKit.jl).

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [April 14, 2019, 1:54pm UTC](https://discourse.julialang.org/t/arpack-puzzles-ask-the-maintainers/23082/11 "2019-04-14T13:54:02Z")

</div>

Thanks, that is good to know about.