Fixing frequency-band mode: Adaptive Lanczos

[Bruno02468]

Some background

MYSTRAN, when solving eigen problems (SOL 103 or the eigen-step of SOL 105), supports either being told a number of modes ahead of time or receiving a frequency band. When working with a frequency band, it used to run an estimator based on Sylvester's law of inertia using the subroutine EST_NUM_EIGENS_BANDED.

The subroutine's name already indicates the first issue: it works with banded matrices, and ever since we got rid of BANDIT, anything that expects banded matrices can end up using massive amounts of time and memory: a deck with around 7k QUAD elements will require over 24 GB of memory for RFAC. That's one of the reasons why @victorkemp added SuperLU to the eigensolve step (RFAC is also present there).

Even if we brought back resequencing or adapted the subroutine to use sparse subroutines for the tridiagonal factorisation step, it's still inefficient. Bill knew this: he made it so that the eigen estimation subroutine never even runs when there are more DoF in the L-set than the PARAM EIGESTL, which has a default value of 5000. That's quite small, even for an L-set. If you increase EIGESTL, you're back to the problem RFAC with 24 GB and days of runtime.

All of that means that any eigen problem with over 5k DoF in the L-set, even if now quickly solvable thanks to @victorkemp's great work, will not run the eigen estimator, and thus default to 1 mode. So... frequency-band mode is broken for all but small models.

Why Lanczos?

Lanczos is the eigenvalue extraction method of choice for larger models and more modes, so it's the ideal starting point for a fix. Inverse-power is usually meant to extract just one (the fundamental) mode, and Givens (modified or not) is just not that powerful.

The existing structure

MYSTRAN's Lanczos procedure mainly lives in two modules:

• ARPACK_LANCZOS_EIG: contains the ARPACK Lanczos subroutines in F77, the most important being DSBAND, the "driver" routine.
• EIG_LANCZOS_ARPACK: called by LINK4, that's where eigen estimation happens, where DSBAND is called, etc.

ARPACK's DSBAND is where the actual eigensolve happens: the reverse communication loop, the factorisation of KMSM, etc. A key detail is that it requires a number of desired eigenvalues to be supplied from the start.

MYSTRAN's EIG_LANCZOS_ARPACK is responsible for calling DSBAND with the right parameters, including the number of eigenvalues (henceforth abbreviated to NEV).

Making it "adaptive"

The core idea

My idea was to make an adaptive version of MYSTRAN's Lanczos eigensolve procedure. The core of said idea was to start with some NEV (say, 10), ask Lanczos for that many, and filter for the ones inside the frequency band. Then double NEV, run Lanczos again, filter. Deciding when to stop is the hard part.

1. There is a hard limit on doublings (10).
2. If the search is one-sided (meaning the lower end is zero), stop when the in-range mode count stabilises.
3. If the search is two-sided, stop when modes outside both ends are found.

That can be achieved by just a new MYSTRAN routine that calls DSBAND more than once, handles reallocation as arrays grow up, the NEV-doubling loop, etcetera. That new routine was born as EIG_LANCZOS_ARPACK_ADAPTIVE.

Reusing the factorisation

One of the first things DSBAND has to do is factorise KMSM (K - σM). The "adaptive" subroutine we just saw calls DSBAND repeatedly with an increasing NEV. That means KMSM is factorised over and over again, wasting time.

So, to make this "adaptive Lanczos" efficient, a new DSBAND had to be put together. It's basically the same, except it takes a pre-done factorisation. Those subroutines live in the new module DSBAND_PREFAC.f F77 module. It's almost the same as DSBAND, and only the adaptive Lanczos procedure uses it.

Why two new modules?

I wanted to ensure this would not break the existing Lanczos procedure. So even though EIG_LANCZOS_ARPACK_ADAPTIVE could have bene made into a drop-in replacement for EIG_LANCZOS_ARPACK, I have modified LINK4 to only call the adaptive one for SOL 103, frequency-band problems.

Limitations

There are some drawbacks to this approach:

1. Inverse-power, Givens and modified-Givens are untouched by this work.
2. I did not cover SOL 105 (buckling) because it uses DSBAND's mode 2, and neither @victorkemp nor I could get it to work on mode 3 (even though it'd make more sense).
3. If modes are too clustered together, two-sided searches can stop early for small bandwidths, missing some modes near the ends. I have not observed this myself, but it's a possibility.

Giving ARPACK more juice

I have changed some ARPACK defaults. The PARAM ARP_TOL, the tolerance used with ARPACK used to have 1e-6 as the default, and I changed it to 1e-9. I also changed MXITERL's default from 50 to 300 -- that's max iterations for Lanczos. The reason for that is these defaults were written with early-2000s computers in mind. The 50-iteration default is particularly abysmal, 300 is alright (and not often reached in practice anyway).

I have also changed the Krylov subspace factor. Put simply, Lanczos requires you to specify the size of the Krylov subspace used in the solve, and it's usually a multiple of the number of requested eigenpairs. The minimum that makes mathematical sense is 1. MYSTRAN uses 2 (EIG_NCVFACL, a constant in MODEL_STUF). Nowadays, 3 is the standard recommendation, so I changed EIG_NCVFACL to 3. The performance impact was minor.

These changes have not majorly impacted results or performance for most use cases, but during development, bad defaults masked other issues (such as the MIN4 shell's issues with rigid-body modes in free-free decks because K6ROT partially grounds the model).

Final remarks

I have done some tests, and this looks sound. I'm not worried about the eigenpairs themselves (they are still calculated the same way). Still, I'd like to see your comments, especially on DSBAND_PREFAC -- I have little experience with Fortran 77.

I hope I got this one right so we can finally say we support frequency-band Lanczos. It's an important feature, so even though we discussed this extensively on Discord, test decks and feedback are more than welcome!