172 CHAPTER 9. ARNOLDI AND LANCZOS ALGORITHMS

√

Reorthogonalization takes place in the j-th Lanczos step if maxk(ωj+1,k) > macheps. qj+1 is orthogonalized against all vectors qk with ωj+1,k > macheps3/4. In the following iteration step also qj+2 is orthogonalized against these vectors. In Table 9.2 the base10 logarithms of the values |wi,j |/macheps obtained with this procedure are listed where |wi,j | = |qi qj |, 1 ≤ j ≤ i ≤ 50 and macheps ≈ 2.2 · 10−16. In Table 9.3 the base-10 logarithms of the estimates |ωi,j |/macheps are given. The estimates are too high by (only) an order of magnitude. However, the procedure succeeds in that the resulting {qk} are semi-orthogonal.

9.7Block Lanczos

As we have seen, the Lanczos algorithm produces a sequence {qi} of orthonormal vectors. These Lanczos vectors build an orthonormal basis for the Krylov subspace Kj (x) = span{q1, . . . , qj } Rn. The restriction of A to Kj (x) is an unreduced tridiagonal matrix. However the Lanczos algorithm cannot detect the multiplicity of the eigenvalues it computes. This limitation prompted the development of the block version of the Lanczos process (Block Lanczos algorithm), which is capable of determining multiplicities of eigenvalues up to the block size.

The idea is not to start with a single vector q1 Rn but with a set of mutually orthogonal vectors which we take as the columns of the matrix Q1 Rn×p with the block size p > 1.

Associated with Q1 is the ‘big’ Krylov subspace

(We suppose, for simplicity, that Aj−1Q1 has rank p. Otherwise we would have to consider variable block sizes.)

The approach is similar to the scalar case with p = 1: Let Q1, . . . , Qj Rn×p be pairwise orthogonal block matrices (Qi Qk = O for i 6= k) with orthonormal columns (Qi Qi = Ip for all i ≤ j). Then, in the j-th iteration step, we obtain the matrix AQj and orthogonalize it against matrices Qi, i ≤ j. The columns of the matrices are obtained by means of the QR factorization or with the Gram–Schmidt orthonormalization process. We obtained the following:

Algorithm 9.5 Block Lanczos algorithm

1: Choose Q1 Fn×p such that Q1Q1 = Ip. Set j := 0 and Fn×p V := 0.

ˆ

This algorithm generates a block tridiagonal matrix Tj with the diagonal blocks Ai,

i ≤ j, the lower diagonal blocks Bi, i < j, and the Krylov basis [Q1, . . . , Qj ] of Kjp(Q1).

2:for j ≥ 0 do

3:if j > 0 then

4:V =: Qj+1Bj ; /* QR decomposition */

5:V := −Qj Bj ;

6:end if

7:j := j + 1;

8:Aj := Qj V ;

9:V := V − Qj Aj ;

10:Test for convergence (Ritz pairs, evaluation of error)

11:end for