Computational Routines

  This chapter will discuss the implementation details of the main computational routines   of ARPACK. We first give an outline of the code structure. This shows how the Implicitly Restarted Arnoldi/Lanczos Method that was described in Algorithm 3 of Chapter 4 is modularized. Each of the basic steps described in the algorithm defines a subroutine module, i.e. a computational routine. The basic tasks and salient implementation details are presented here for each of these computational routines.

ARPACK relies heavily upon a number of basic operations and algorithms provided by the BLAS  and LAPACK . These have contributed greatly to the robustness, accuracy, and computational performance of ARPACK. The most important of these with respect to performance is the BLAS subroutine _gemv . For a fixed number nev of requested eigenvalues and a fixed length ncv Arnoldi basis, the computational cost   scales linearly with n, the order of the matrix. The rate of execution   (in FLOPS) for the IRA iteration is asymptotic to the rate of execution of _gemv.

In the outline of the implementation described in Figure 5.1, is a upper Hessenberg matrix, , and the residual vector is -orthogonal   to the columns of . For each j,

where OP and are defined with respect to one of the computational modes described in Chapter 3.

The integer k denotes the desired number of eigenvalue approximations and this may, at times, be greater than the users request for nev approximations. The integer will be the largest size factorization constructed. An eigenvalue of is a Ritz value of OP and is the corresponding Ritz vector when (A formal definition of Ritz value and vector is given in Chap. 4). The normalization is assumed throughout.

  

Figure 5.1: XYaupd - Implementation of the IRAM/IRLM in ARPACK

• Perform basic error checks, partition the internal workspace. Set $k={\tt nev}$ and $m={\tt ncv}$.

• Enter XYaup2. Generate a random initial vector ${\bf V}_m {\bf e}_1 ={\bf v}_1$ by calling Xgetv0, unless the initial vector is provided by the caller.

• Call XYaitr to compute the initial Arnoldi/Lanczos factorization
  ${\bf {\tt OP} V}_k = {\bf V}_k {\bf H}_k + {\bf f}_k {\bf e}^T_k$ of length $k = {\tt nev}.$

• For iter $= 1,\ldots,$ maxiter $+1$,
  1. Call XYaitr to extend the length $k$ Arnoldi/Lanczos factorization to a length $m$ factorization. Reverse communication is performed to compute matrix-vector products with $\bf {\tt OP}$ and possibly ${\bf B}$.

  2. Compute the eigenvalues of ${\bf H}_m$ and the associated error bounds. Call Xseigt for the symmetric eigenvalue problem or Xneigh otherwise.

  3. Call XYgets to partition the eigenvalues into two sets $\Omega_w$ and $\Omega_u$. The $k$ eigenvalues in the set $\Omega_w$ are the desired approximations while the remaining eigenvalues in the set $\Omega_u$ are to be used as shifts.

  4. Call XYconv to determine how many of the wanted Ritz values satisfy the convergence tolerance.

  5. Exit the loop if all $k$ eigenvalues in $\Omega_w$ satisfy the convergence criterion, or if iter > maxiter.

  6. Possibly increment $k.$ Determine $p=m-k$ shifts $\{\mu_j \}_{j=1}^p$. If the exact shift strategy is used, the eigenvalues of $\Omega_u$ are used as shifts, otherwise the $p$ shifts are provided through a reverse communication interface. If $p=0$, exit the loop.

  7. Implicitly restart the Arnoldi/Lanczos factorization by calling XYapps.
     1. Perform $p$ steps of the implicitly shifted QR algorithm on ${\bf H}_m$ with shifts $\{\mu_j \}_{j=1}^p$ to get ${\bf H}_m {\bf Q}_m = {\bf Q}_m {\bf H}_m^{+},$ where ${\bf H}_m^{+}$ and ${\bf Q}_m$ are upper Hessenberg and orthogonal matrices, respectively.

     2. Let ${\bf Q}_k$ denote the matrix consisting of the first $k$ columns of ${\bf Q}_m$ and ${\bf H}_k^{+}$ the leading principal submatrix of ${\bf H}_m^{+}$ of order $k.$
        Update ${\bf V}_k^+ = {\bf V}_m {\bf Q}_k$ to get the new length $k$ Arnoldi factorization ${\bf {\tt OP}} {\bf V}_m^+ = {\bf V}_m^{+} {\bf H}_k^{+} + {\bf f}_k^{+} {\bf e}^T_k.$ (See Algorithm 3, Chapter 4 ).

• End For
• Call XYeupd for computing Ritz and/or Schur vectors and for transforming the Ritz values of $\bf {\tt OP}$ if a spectral transformation was used.