Filtering Frequencies in a Shift-and-invert Lanczos Algorithm for the Dynamic Analysis of Structures

The shift-and-invert Lanczos algorithm is a commonly used solution procedure to compute the eigenpairs of large, sparse eigenvalue problems that arise when approximating the elastic dynamic response of large structures under the influence of seismic forces. Not all eigenvectors are equally important to the response when the structure is exposed to a mass-dependent external force of the form g(t)Mb, where M is the mass matrix of the system and b the rigid body vector. Structural engineers select eigenvectors xj , j = 1, . . . , `, such that their cumulative mass participation, measured as ∑` j=1(x T j Mb) 2/(bTMb), is above a target threshold ξ. We show that when the starting vector for the unshifted Lanczos algorithm is the spatial distribution vector b, the Lanczos procedure can be used to provide an estimate of the cumulative mass participation. This allows us to identify intervals containing eigenvalues whose eigenvectors have a large contribution to the cumulative mass participation and filter out intervals containing eigenvalues whose eigenvectors have a negligible contribution. We use this information to devise a sequence of shifts σ1, . . . , σp for the shift-and-invert Lanczos algorithm as well as a stopping criterion for the iteration with shift σi so that the cumulative mass participation of the computed eigenvectors reaches the required level ξ. Numerical experiments on real engineering problems show that our approach computes up to 80% fewer eigenvectors and requires fewer shifts, on average, than the more general shifting strategy proposed by Ericsson and Ruhe (Math. Comp., 35 (1980)).


Introduction.
A structural dynamics problem consists of finding the response of a structure, for instance, a building or a bridge, given some dynamic loading.Such problems may be written in the form of a system of second order differential equations (1.1) M ü(t) + D u(t) + Ku(t) = f (t) that results from the finite element discretization of the equation of motion, together with some initial conditions.The mass matrix M ∈ R n×n is usually symmetric positive semidefinite (denoted by M ≥ 0), the stiffness matrix K ∈ R n×n is symmetric positive definite (K > 0), the damping matrix D is symmetric and often positive definite, u(t) is the displacement, and f (t) is the time-dependent external load on the given structure.Here, we concentrate on external forces of the form where g(t) is a scalar function and 0 = b ∈ R n is the spatial distribution vector (also called rigid body vector or spatial vector of loading patterns).External forces of the form (1.2) are particular to earthquake loading, where g(t) is the input earthquake acceleration.
Projection methods are usually employed to reduce the dimension n of the system (1.1).These methods consist of constructing a matrix X ∈ R n× of full rank and transforming (1.1) into the reduced system (1.3)X T M X v(t) + X T DX v(t) + X T KX v(t) = g(t)X T M b, which is then solved for v(t).An approximate solution to (1.1) is obtained as u(t) ≈ X v(t).In the mode superposition method, the columns x 1 , . . ., x of X are eigenvectors corresponding to finite eigenvalues of the associated generalized eigenvalue problem (GEP) (1.4) (K − λM )x = 0, where we assume generalized proportional damping, so that the damping matrix D in (1.1) is diagonalizable by the matrix of eigenvectors of (1.4).If the columns of X are M -orthonormal, i.e., x T i M x j = δ ij with δ ij the Kronecker delta, then the reduced system (1.3) can be rewritten as decoupled second order differential equations, (1.5) vj (t) + 2ζ j ω j vj (t) + ω 2 j v j (t) = g(t)x T j M b, j = 1, . . ., , where ω j = λ j , λ j > 0 is an eigenvalue of (1.4) with corresponding eigenvector x j and x T i Dx j = 2ζ j ω j δ ij for some ζ j ≥ 0 [16,Chap. 18].Not all eigenvectors are equally important to a given system of differential equations and it is clear from (1.5) that the response v j depends on both the frequencies ω j and the magnitude of x T j M b, which is called the mass participation factor of x j .They satisfy (1.6) (see Section 2.1).Within the context of seismic analysis and design, structural engineers aim to achieve 80 to 90% of mass participation in (1.6) using only a subset of the eigenvectors (see the justification in section 2.1).This leads to the following problem.
Problem 1.1 is easy to solve if we can compute all the eigenvectors of the GEP (1.4), but this is not feasible for problems of large dimensions.It is usually the eigenvectors corresponding to small eigenvalues that contribute the most to the total mass participation (see for example [2] or Section 2.1), so in previous work Problem 1.1 was relaxed to the following problem.
Problem 1.2.For a given proportion ξ ∈ (0, 1) and a spatial distribution vector b, find the M -orthonormal eigenvectors x i , i = 1, . . ., associated with the smallest eigenvalues of (1.4), where is the smallest integer such that The eigenvectors (modes of vibration) of a tall building structure.From left to right: the undeformed shape of the structure; its first eigenvector, x 1 corresponding to the smallest eigenvalue, which is predominantly a "sway mode" (i.e., large mass participation factor in the cartesian x-direction); and another eigenvector, x 20 , which is a "bouncing mode"(large participation along cartesian z-direction).Deformations are exaggerated.

Table 1
Mass participation factors for M -orthonormal eigenvectors x 1 and x 20 from Figure 1.Recovering sufficient mass participation can be a challenge depending on the geometry of the structure and orientation of b.The vector b represents the "rigid body" deformation of the structure in cartesian x-, y-, and z-directions.For instance a given structure may have a geometry that makes it easy to solve Problem 1.2 with only a handful of smallest eigenvectors (i.e., is small) when b is along the x-direction but might require a large value of when b is along z.This is graphically illustrated in Figure 1 which shows the first and twentieth eigenvectors, x 1 and x 20 , overlaid as displacements on the structure's geometry.As can be seen from the diagram in the middle, x 1 is a sway mode wherein large parts of the structure deform in the cartesian x-direction.By contrast, x 20 (right of middle) is a bouncing mode and therefore the structure deforms significantly in the z-direction.The mass participations of these eigenvectors are tabulated in Table 1 for b in the x-and z-directions.
A natural approach to solve Problem 1.2 is to apply the shift-and-invert Lanczos algorithm (see [1], [17] and references therein) to K − λM with or without shifts depending on how many eigenvectors will be needed to satisfy (1.8).Instead of using eigenvectors in (1.8), Wilson et al. [20] proposed to use the Ritz vectors resulting from a variant of the Arnoldi algorithm later known as the WYD algorithm [9], [19].The use of the Lanczos vectors was later suggested by Nour-Omid and Clough [13].
But as noted in [2], there is no analysis to support the use of Lanczos vectors and Ritz vectors (unless they have converged to eigenvectors of K − λM ).In particular, we show in section 2 that under suitable conditions on g(t) and the initial conditions for (1.1), the quantity (1 − ξ) 2 provides an upper bound on the relative error between the response u(t) to the ODE (1.1) and its approximation by X v(t) when the columns of X are eigenvectors chosen such that (1.7) or (1.8) holds.So we concentrate here on the computation of eigenvectors for Problems 1.1-1.2.
Given an initial shift σ 1 , a common strategy to determine the sequence of shifts σ 2 , σ 3 , . . . to be employed in a shift-and-invert Lanczos process for symmetric GEPs is to choose the new shift σ i , i > 1, such that the largest converged eigenvalue λ max is halfway between the old shift σ i−1 and the new shift σ i , namely σ i = 2λ max − σ i−1 [5].This way, if the eigenvalues are roughly evenly distributed across the spectrum then the same number of eigenvalues is expected to converge to the left of the new shift as the number of eigenvalues that have converged to the right of the old shift.The Lanczos iteration with shift σ i is then stopped when the smallest converged eigenvalue in that iteration coincides with the largest converged eigenvalue with shift σ i−1 .One issue with this shifting strategy for Problem 1.2 is that the shift-and-invert Lanczos algorithm may return too many eigenvectors with small normalized mass participation.As a result, the number of returned eigenpairs may be unnecessarily large, whereas eigenvectors with large normalized mass participations corresponding to eigenvalues that are not in the lower end of the spectrum may not be detected.
Our main contribution is the presentation of a new shifting strategy for the shiftand-invert Lanczos algorithm together with a stopping criterion for the iteration with shift σ i that are specifically designed to approximate the solution to Problem 1.1.For this, we use the theory of orthogonal polynomials to show that a few steps of the unshifted inverse Lanczos algorithm applied to K −1 M with starting vector b provides, at no additional cost, information about the location of the eigenvalues whose corresponding eigenvectors have non-negligible mass participation and also helps to identify intervals where eigenvalues have negligible mass participation.We use this information to devise a shifting strategy for the shift-and-invert Lanczos process so that condition (1.7) holds, albeit perhaps not for the smallest number of eigenvectors.This shifting strategy performs especially well in the cases where the eigenvectors with non-negligible mass participation correspond to larger eigenvalues while there are intervals of smaller eigenvalues whose corresponding eigenvectors have negligible contribution to the overall response of the structure.These are the cases that are the most problematic for available methods.
Numerical experiments performed on real structural engineering problems show an often large reduction in the number of eigenvectors computed to approximate the solution to Problem 1.1 using our new shifting strategy as opposed to the number of eigenvectors computed to approximate the solution to Problem 1.2 using the more general shifting strategy of Ericsson and Ruhe.The use of a new shift in the shiftand-invert Lanczos process has a cost since it leads to a new matrix factorization.For our set of test problems, our numerical experiments show that the number of shifts used with our new shifting strategy is, on average, smaller than with Ericsson and Ruhe's strategy.
Unlike previous attempts to solve Problems 1.1-1.2,we pay special attention to issues that can occur when the mass matrix M is singular.In this case, the Lanczos process proceeds with a quasi-inner product.In finite precision arithmetic, the computed Lanczos vectors and Ritz vectors have components in the null space of M and the magnitude of these unwanted components grows during the iterations.This can either delay or prevent convergence of the Ritz vectors [4].These issues are overcome with the use of an appropriate starting vector and implicit filtering [12].
In the next section we give some preliminary material that includes a justification of Problem 1.1 and a discussion of issues that arise when using a quasi-inner product in a shift-and-invert Lanczos process, as well as possible remedies.We describe and justify our new shifting strategy in section 3, and illustrate its performance on a number of real structural engineering problems in section 4.
2. Preliminaries.Following [2], we show in section 2.1 that, under some assumptions on the initial conditions and the input function g(t), choosing ξ close to 1 in (1.7) guarantees a small error between the exact response u(t) and its approximation as a linear combination of the eigenvectors x i k , k = 1, . . ., satisfying (1.7).
We recall in section 2.2 the shift-and-invert Lanczos process for the GEP (1.4) and discuss issues related to the use of a quasi-inner product defined by the symmetric semidefinite matrix M .
2.1.Upper bound on the response error.The n×n mass matrix M in (1.1) is often singular in applications.As a result, the GEP (1.4) has the eigendecomposition where I r is the r × r identity matrix with . ., λ n ) has real positive diagonal entries displaying the r finite eigenvalues as λ j , j = 1, . . ., r (the remaining n − r eigenvalues being at infinity), and X is a nonsingular matrix containing the corresponding eigenvectors x 1 , . . ., x n .Note that the eigenvectors x j , j = r + 1, . . ., n associated with the eigenvalues at infinity belong to the null space of M , i.e., (2.1) M x j = 0, j = r + 1, . . .n.
where v j (t) is the jth entry of the vector v(t).We assume generalized damping so that for some ζ j ≥ 0 and ω j = λ j , j = 1, . . ., n, and rewrite (1.1) as The system of uncoupled equations (2.3) can then be solved by direct integration, yielding where we let ω j = ω j (1 − ζ j ) 1/2 and assume for simplicity that u(0) = u(0) = 0 so that v(0) = v(0) = 0.As mentioned in the introduction, in practical applications, n is large and it is infeasible to compute all n eigenpairs of (1.4) so the solution u(t) in (2.2) is approximated instead by We use the M -quasi vector norm, to measure the relative error between the exact solution and its approximation, We rewrite v j (t) in (2.4) as v j (t) = h j (t)x T j M b.If the spatial distribution vector b has no components in the null space of M (this is usually the case), i.e., b = r j=1 b j x j , then which on using (2.1) is equivalent to (1.6).Now if we assume that h min ≤ |h j (t)| ≤ h max for t > 0, j = 1, . . ., r, and some positive scalars h min , h max , then it follows from (2.6) and (2.7) that where ξ is as in (1.7).Thus, under the above assumptions on the initial conditions u(0), u(0), and the functions h j , choosing a ξ close to 1 guarantees a small relative error between the exact solution u and its approximation ũ.We refer to section 4 for an illustration of the above analysis.Finally, we note that the factor 1/ω j in (2.4) suggests that the eigenvectors corresponding to smaller eigenvalues ω 2 j (i.e., lower frequencies) are more likely to have a larger contribution to the response u(t) in (2.2) than those corresponding to the higher frequencies.

Shift-and-invert Lanczos process with semi-definite inner product.
Applying the Lanczos algorithm with an M -quasi inner product to compute approximate eigenpairs of the definite pencil K − λM with M ≥ 0 was first suggested by Scott [18].Given a shift σ near the eigenvalues of interest and a starting vector w, this Lanczos procedure constructs a matrix Q k whose columns, called the Lanczos vectors, form an M -orthonormal basis for the kth order Krylov subspace w/ w M , . . ., q k and q 0 = 0, the next Lanczos vector q k+1 is obtained from the three-term recurrence (2.8) where The three-term recurrence (2.8) can be rewritten in matrix form as where Since, by construction, q k+1 is M -orthogonal to the columns of Now, for an eigenpair (θ is called a Ritz pair for K − λM and is considered as an approximate eigenpair of K − λM if the scaled residual (which approximates the backward error for (λ j , x j ) [10, Thm.2.1]) (2.12) is below a given tolerance.
An algorithm implementing the three-term recurrence (2.8) is provided in Algorithm 2.1.This is essentially the algorithm provided in [14].Breakdown occurs when β k+1 ≤ 0 in step 10.The number of eigenpairs returned will depend of the convergence criterion in step 3. We discuss the latter in section 3.3.Some more comments follow.
(a) Choice of starting vector and implicit filtering.When M is singular, Nour-Omid et al [14, Sec.2.2] recommend choosing a starting vector w in the range of (K − σM For an eigenvector x i of K − λM with finite eigenvalue λ i , we have that (2.13) If the starting vector w for the shift-and-invert Lanczos procedure is in the range of (K − σM ) −1 M then so is the first Lanczos vector q 1 and the subsequent Lanczos vectors q 2 , . . ., q k if operations are performed in exact arithmetic.As a result, the Ritz vectors lie in range((K − σM ) −1 M ).Now if Algorithm 2.1 Shift-and-invert Lanczos algorithm This algorithm takes as input two n×n symmetric matrices M ≥ 0 and K > 0, a shift σ ≥ 0 such that K − σM is nonsingular, and a starting vector 0 = w ∈ range((K − σM ) −1 M ).It returns eigenpairs (converged Ritz pairs) (λ ij , x ij ), j = 1, . . ., .
12 Check for convergence of the Ritz pairs (λ j , x j ) = ( 1 ) then the Ritz vectors may have unwanted components in the null space of M , which in turn slows down their convergence.Note that in finite precision arithmetic, even if w ∈ range((K − σM ) −1 M ), rounding errors prevent the computed Lanczos vectors from staying in the range of (K − σM ) −1 M .The unwanted components in null(M ) are mainly introduced when solving the linear systems with K−σM .Since this operation is performed at each iteration, the accumulation can be rapid.The latter is set off at the beginning when the starting vector w is constructed as the product of (K − σM ) −1 M with some y ∈ R n , since this operation itself can introduce a null space component.The starting vector can be put in the right space explicitly by forming a particular projection matrix [4, Section 2.3], or solving the linear system to higher precision, however both methods come at a substantial cost.Instead, as suggested by Meerbergen [12], we apply an implicit filter that alters the starting vector implicitly, producing Lanczos vectors lying in range((K − σM ) −1 M ).This approach is briefly discussed below.For q = q R + q N ∈ R n with q R ∈ range(M ) and q N ∈ null(M ), we have that q M = q R M , i.e., the components in the null space of M are undetectable by the M -norm.So q N 2 can be arbitrarily large even when q M = 1.
Hence, if at step k of the Lanczos process we have for some tolerance tol 1 (we choose tol = 10 4 ) then we apply implicit filtering.For this let (2.15) , and if we let then it is not difficult to show that the matrix T k−1 obtained from T k−1 by deleting its last row is tridiagonal and satisfies the Lanczos recurrence where The implicit filter implicitly premultiplies the Lanczos vectors by K −1 M , thereby removing any components in null(M ).Since the dimension of the projection space is reduced by one, we continue the Lanczos algorithm by forming the next Lanczos vector qk+1 from qk as well as updating T k in step 11 with T k−1 .Note that the implicit filtering described above is essentially one step of the unshifted QR algorithm applied to T k .A Lanczos vector q k with a large 2-norm means that q k has large components in the null space of M .As a result, the smallest eigenvalue of T k becomes close to zero and the unshifted QR step pushes that very small eigenvalue to the bottom of the tridiagonal matrix.Indeed, if T k had a zero eigenvalue then after one step of unshifted QR, the last row of T k would be zero.The construction of T k−1 corresponds to a deflation of the eigenpair of T k with smallest eigenvalue, thereby removing the Ritz pair in (2.11) with a large (possibly infinite) eigenvalue.(b) Testing for convergence of eigenpairs.As explained earlier, we consider the Ritz pair (λ j , x j ) := (1/θ , where (θ j ) is an eigenpair of T k , to have converged if its backward error η(λ j , x j ) in (2.12) is below a given tolerance.We suggest to only compute η(λ j , x j ) when the Ritz pair (λ j , x j ) is likely to have converged.This can be checked as follows.On using (2.9), we have that , where we used Cauchy-Schwarz for the first inequality.Without loss of generality, we can assume that All the quantities in the above approximate upper bound are readily available during the Lanczos steps including q k+1 2 (see point (a) and (2.14)).Hence, we suggest to only compute η(λ j , x j ) when the upper bound in (2.17) is below, say, 10 × tol.

Shifting strategies for an approximate solution to problems 1.1-1.2.
A considerable number of eigenvectors may be required to solve Problem 1.2 for problems where there are large normalized mass participation factors corresponding to eigenvalues λ j that lie further away from the lower end of the spectrum.In this case, the Lanczos algorithm (Algorithm 2.1 with σ = 0 as we are aiming for the small eigenvalues first) becomes increasingly slow and memory intensive.So if condition (1.8) is not satisfied by the converged eigenvectors after a given number, say k max , steps of Algorithm 2.1 with zero shift, we then restart Algorithm 2.1 with a sequence of shifts.As mentioned in the introduction, Ericsson and Ruhe [5] propose to use the following sequence: (3.2) where λ max is the largest converged eigenvalue.Algorithm 2.1 with shift σ i is stopped when the smallest converged eigenvalue coincides with the largest converged eigenvalue with the shift σ i−1 , or once condition (1.8) is satisfied, in which case we terminate the computation.The shift-and-invert Lanczos algorithm with this shifting strategy, which we refer to as SIL, approximates the solution to Problem 1.2.An issue with Problem 1.2 is that its solution may include many eigenvectors with very small or negligible normalized mass participation factors.Also, the shifting strategy (3.2) is not as effective when K − λM has clustered eigenvalues, as is the case for the real structural engineering problems we consider in section 4. In what follows, we show that Algorithm 2.1 with σ = 0 and starting vector w equal to the spatial distribution vector b provides at almost no additional cost information about where the eigenvalues associated with eigenvectors of non negligible normalized mass participation lie, while at the same time identifies the parts of the spectrum that do not contribute much to the total mass participation.We use this information to devise a new shifting strategy.
3.1.Estimating the cumulative mass participation.Let q 1 , q 2 , . . .be the Lanczos vectors generated by the three-term recurrence (2.8) with σ = 0 and 0 = q 1 ∈ range(K −1 M ).Each Lanczos vector q i can be written as where p i is called the ith Lanczos polynomial.It is well known from the theory of orthogonal polynomials [6, Chap.2], [8,Chap. 4] that these polynomials are orthogonal with respect to the inner product defined in terms of the Riemann-Stieltjes integral where a ≤ µ min , and b ≥ µ max , with µ min and µ max the smallest and the largest eigenvalues of K −1 M , respectively (i.e., µ min = 1/λ max and µ max = 1/λ min where λ max and λ min are the largest and smallest eigenvalues of K − λM ).The distribution function φ(µ) is a step function with jumps at the eigenvalues µ i of K −1 M , and is given by where is the Heaviside function and the φ i 's are the coefficients of the first Lanczos vector q 1 when expressed in the M -orthonormal basis x 1 , . . ., x r for range(K −1 M ), namely ) (and φ i = 0 for i = r + 1, . . ., n).
It turns out that the step function φ(µ) coincides with the cumulative mass participation sum when the starting vector for the Lanczos algorithm with zero shift is the spatial distribution vector b in (1.2) with b ∈ range(K −1 M ), as we now show.
where m(x i ) is the normalized mass participation of the eigenvector x i in (3.1).
Proof.The first Lanczos vector is given by q 1 = b/ b M and, for i = 1, . . ., r, For i = r + 1, . . ., n, x i ∈ null(M ) so that M x i = 0 and hence m(x i ) = 0 = φ 2 i .If the nonzero eigenvalues µ j , j = 1, . . ., r, of K −1 M are ordered by decreasing values then it follows from (3.5) and Proposition 3.1 that i ), i = 1, . . ., k, be the eigenpairs of the tridiagonal matrix T k resulting from k steps of the unshifted Lanczos algorithm (Algorithm 2.1 with σ = 0) and ordered such that θ k .The Lanczos polynomials p i in (3.3) are not only orthogonal with respect to the inner product (3.4),they are also orthogonal with respect to the inner product induced by the step function where τ k,i is the first entry of the eigenvector s (see [6,Chap. 2], [8]).As a result, τ k (µ) and the distribution function φ(µ) in (3.5) have the same modified moments up to degree 2k − 1, namely In turn, by the following theorem due to Karlin and Shapley [11,Thm. 22.2], it follows that τ k (µ) serves as a good approximation to φ(µ).Thus if τ k and φ do not coincide, the vertical and horizontal steps of τ k will intersect φ exactly 2k −1 times.This theory has been used, for instance, in estimating eigenvalue distribution [6], and in constructing polynomial preconditioners [7].In our case, we use the step function τ k , obtained after k steps of the Lanczos algorithm, as an approximation to the cumulative mass participation sum, i.e., This is illustrated in Figure 2(a)-(b) for a real structural engineering problem called chilled problem.On plot (a), the step functions τ k and φ appear to coincide at least to the eye but plot (b), which is a closeup of plot (a) around λ = 7.3 × 10 3 shows that the step functions indeed intersect.We will return to Figure 2 at the end of section 3.2.

3.2.
A new shifting strategy.The step function τ k (µ) in (3.8) is readily available after k steps of Algorithm 2.1 with σ = 0 and starting vector w = b.Intuitively, the eigenvalues corresponding to the eigenvectors with largest normalized mass participation should lie under tall and narrow steps of φ(λ), whereas short and wide steps would indicate intervals of eigenvalues corresponding to eigenvectors of negligible mass participation (see Figure 2 for illustration).To make this formal, denote by  Step functions φ(µ) and τ k (µ) for the chilled problem in the x-direction.The small yellow stars in (a) correspond to the smallest eigenvalues λ i and their normalized mass participation factors m(x i ) such that i m(x i ) ≥ 0.9.The green and red circles correspond to the pairs (λ i , m(x i )) obtained by MASIL and pMASIL, respectively, and such that i m(x i ) ≥ 0.9.The shading shows the intervals selected by the new shifting strategy.

Ψ(c, d) the total mass participation of eigenvectors whose corresponding eigenvalues lie in the interval
where λ j = 1/µ j is the jth largest eigenvalue of K − λM .Since by Theorem 3.2 the function φ(µ) − τ k (µ) has exactly 2k − 1 sign changes , we have that for 1 We can now make use of the bounds in (3.11) to identify the union of intervals of smallest total length over which we are guaranteed to satisfy (3.12) for some ξ ∈ (0, 1).Since we require that the total mass participation over that union of intervals is at least ξ, we must look a the lower bounds of (3.11).Let us denote by γ ij the following ratios , and small γ ij indicate low relative mass participation.Since the lower bound for Ψ(1/θ ]. Thus to simplify the following discussion, define γ i to be Suppose that after k steps of the Lanczos algorithm the first Ritz pairs (λ ν , x ν ), ν = 1, . . ., of K − λM have converged.Let us denote by ξ the sum of their mass participation factors, that is, To construct a union of intervals of smallest total length over which we are guaranteed to attain the remaining mass participation ξ − ξ , we pick the s largest γ i , say The wanted union of intervals is then . By merging the neighbouring and the overlapping intervals we can construct a set of s ≤ s disjoint intervals jν ], ν = 1, . . ., s .
We then choose the shifts σ ν to be the midpoints of those intervals, namely, We end the shift-and-invert Lanczos iteration with shift σ ν whenever the sum of the mass participation of the converged Ritz vectors with eigenvalues in [1/θ jν ] attains the minimum in (3.11), i.e., jν −1 =iν +1 τ 2 k, .Remark 3.3.The inequality (3.13) may not hold if even s = k − + 1 although in practice, it is usually satisfied after a small number of steps k.If (3.13) does not hold then we can increase k, or reduce ξ.
As an illustration, let us look at Figure 2(a).The jumps of τ k (µ) correspond to the points (1/θ There are 34 more jumps outside of plot (a) corresponding to those θ The plot shows as small yellow stars the normalized mass participation factors m(x i ) of the eigenvectors x i corresponding to the smallest eigenvalues λ i , computed by the unshifted Lanczos algorithm and thus solving Problem 1.2.Although not visible on plot (a), the eigenvalues of the chilled problem are clustered: for example, the smallest eigenvalue is around 507, there are 47 eigenvalues in the interval [535,538] and 25 in [570,581].The first shaded region from the left corresponds to the interval [1/θ the converged Ritz values computed by the unshifted Lanczos algorithm.The remaining shaded regions correspond to the s = 9 disjoint intervals in (3.14).They define nine shifts for the chilled problem (there is a small interval just after the first interval that corresponds to the shift σ = 0).The non shaded areas correspond to intervals containing eigenvalues associated with eigenvectors of negligible mass participation.These intervals are ignored by our approach.

3.
3. An algorithm for the approximate solution to Problem 1.1.Given two n×n matrices M ≥ 0 and K > 0, a spatial distribution vector b ∈ range(K −1 M ), a proportion ξ ∈ (0, 1), and a maximum number of iterations k max our algorithm for the approximate solution to Problem 1.1 goes through the following steps.step 1 Call Algorithm 2.1 with σ = 0, w = b, and the implicit filtering turned off.
Stop the Lanczos iterations at step k when either (a) the converged Ritz vectors x j , j = 1, . . ., 0 are such that 0 j=1 m(x j ) ≥ ξ, or (b) k has reached a number of k max iterations (or larger if k j=1 τ 2 k,j ≥ ξ is not satisfied).If (a) holds then return the converged Ritz vectors x j , j = 1, . . ., 0 as an approximate solution to Problem 1.1.End the algorithm.If (b) holds then • if q k 2 ≤ tol q 1 2 then save the converged Ritz vectors x j , j = 1, . . ., 0 , let ξ 0 = 0 j=1 m(x j ) and proceed to step 2 with the converged and unconverged Ritz pairs (1/θ where w is the sum of the unconverged Ritz vectors from the previous step.Apply implicit filtering when q k 2 > tol q 1 2 .Stop the Lanczos iterations when either (i) the converged Ritz vectors x j , j = p+1, . . ., p+ ν with p = 0 + ν−1 i=1 i , are such that Stop the for loop when (ii) holds.Return the converged Ritz vectors x j , j = 1, . . ., p with p = 0 + ν i=1 i as an approximate solution to Problem 1.1.End the algorithm.Implicit filtering cannot be used in the first call to Algorithm 2.1 in step 1 since this would alter the starting vector and we could not apply our shifting strategy.In step 1(b) the condition k j=1 τ 2 k,j ≥ ξ is usually satisfied after a small number of steps k (see remark 3.3).The choice of the maximum number of iterations k max is important.If k max is too small then the approximation τ k (λ) of φ(λ) is too rough and leads to large intervals and shifts that are not close enough to eigenvalues with eigenvectors that have large mass participation.On the other hand, a too large k max can lead to unnecessary computations, in particular of converged eigenvectors with negligible mass participation but can also result in too many shifts being identified by our shifting strategy.This last point is illustrated in Figure 3 for a real engineering problem.The plot shows that the number of shifts (or intervals in (3.14)) determined by our shifting strategy increases as the number k of unshifted Lanczos steps increases.In practice, we found that taking k max = 200 is a reasonable choice.For tol we choose tol = 10 4 , and we consider that the eigenpair (λ j , x j ) has converged if η(λ j , x j ) ≤ nu, where u is the machine precision and η(λ j , x j ) is defined in (2.12).

Numerical experiments.
For our numerical experiments we used matrices M and K, and spatial distribution vectors b provided by Arup Group Limited that were constructed by the finite element software package Oasys GSA [15] from models Table 2 List of test problems together with their size n and the number kmax of unshifted Lanczos steps for SIL and the new method MASIL.The last four columns list the number of eigenvectors needed to satisfy the MPF condition (3.12) with ξ = 0.9 for SIL and MASIL, and their purged versions pSIL and pMASIL.The number of shifts used by SIL and MASIL is provided inside brackets.We compare the following approaches to solve our problem:

Problem
• SIL: shift-and-invert Lanczos algorithm with the Ericsson and Ruhe shifting strategy (3.2), • MASIL: mass accumulating shift-and-invert Lanczos algorithm with the new shifting strategy described in section (3.3).When no shifts are used, SIL and MASIL are identical, except for their starting vectors: SIL uses w = K −1 M b as suggested in [20] and MASIL uses w = b.
In practice, the cumulative sum of the mass participation factors of the eigenvectors returned by SIL and MASIL is always slightly larger than the wanted proportion ξ.So we have the possibility to remove from the list of computed eigenvectors those with smallest mass participation so that the cumulative sum of the mass participation factors of the remaining eigenvectors is exactly ξ or just above.We refer to this small modification of SIL and MASIL as pSIL and pMASIL, respectively, where the "p" stands for extra purging step.An illustration of the purging step can be seen in Figure 2(a), where the pairs (λ i , m(x i )) from MASIL are shown as small green circles and those kept by pMASIL are shown as large red circles.
In Table 2, we compare the number of computed eigenvectors required to satisfy the mass participation condition (3.12) with ξ = 0.9.Some directions for the problems are excluded from the table.This is either because shifts were not employed (e.g.ccnb in the y-or the z-direction, where the 90% mass participation was reached in fewer than k max steps and no shift), or too many eigenvectors were required to achieve 90% mass participation, exceeding time and memory constraints (e.g., for the chilled problem in z-direction and ξ = 0.75 0.9, SIL returns 6959 eigenvectors whereas MASIL returns 4201 eigenvectors).
The new shifting strategy allows us to exclude intervals containing eigenvalues whose eigenvectors have a negligible mass participation and shift in the middle of intervals containing eigenvalues with eigenvectors of large mass participation.As a result, the number of eigenvectors returned by MASIL can be much smaller than that returned by SIL.This is, for example, the case for the local modes problem in the z-direction and the chilled problem in the x-direction.The number of shifts used by SIL or MASIL depends on the problem and the shifting strategy employed.Notice the large number of shifts employed by SIL for the local modes problem in the zdirection and for the TT problem in the z-direction.For MASIL, the number of shifts needed is known before hand so if this number is too high then there is always the possibility to increase the value of k max or to see if there is a k < k max that leads to larger search intervals but fewer of them.We reported results for two different values of k max for the chilled problem in the y-direction.For k max = 100, SIL uses only one shift but returns 107 eigenvectors whereas MASIL uses 4 shifts but returns only 14 eigenvectors.No shifts are needed if we increase k max to 200.Now comparing the number of eigenvectors returned by pSIL to that of SIL shows that a large proportion of the eigenvectors computed by SIL have a negligible mass participation and can be purged away.The last column in Table 2 shows that MASIL still returns many eigenvectors with negligible mass participation that can be removed while still maintaining the condition (3.12) but the reduction is not as drastic as for SIL.Note that the number of eigenvectors returned by pMASIL is lower or equal than those returned by pSIL.
Let us now look at the quality of the approximation ũ(t) = p j=1 v j (t)x j to the response u(t) of (1.1) when the p eigenvectors x j are computed by SIL, MASIL, pSIL or pMASIL.For the time dependent external load f (t) = g(t)M b we use for the loading function g as suggested in [3, Fig. 2].We employ 2% damping which consists of setting ζ j = 0.02 in (2.4) for all j.For the relative error between the exact solution u(t) and its approximation ũ(t), the analysis in section 2.1 tells us to expect (4.1) We do not have access to u(t).Therefore as a reference we use the solution obtained from SIL with ξ = 0.99.The relative errors for the approximate responses ũ(t) returned by SIL, pSIL, MASIL and pMASIL with ξ = 0.9, are shown in Figure 4 for the local modes problem in the z-direction in (a) and for the chilled problem in the x-direction in (b).The relative errors agree with (4.1).Not surprisingly, the relative errors for the pSIL and pMASIL solutions are close to √ 1 − ξ ≈ 0.3.We note that for the chilled problem, the MASIL solution with 14 eigenvectors is almost as good as the SIL solution which requires 101 eigenvectors, whereas for the local modes problem pMASIL algorithm outperforms pSIL algorithm with significantly fewer eigenvectors.
Finally, we plot in Figure 5(a)-(b) the ith entry of the reference solution u(t) and its approximations obtained by SIL, pSIL, MASIL, and pMASIL for the ccnb and TT problems.We chose the index i for which the reference solution has the largest amplitude.For the TT problem in the y-direction in plot (a), all solutions agree with the reference solution even the pMASIL solution which uses significantly fewer  Relative error between the reference solution u(t) and its approximation ũ(t) obtained from SIL, MASIL, pSIL, and pMASIL with a proportion of ξ = 90% for the total mass participation.
eigenvectors.Although (4.1) holds for all the approximate solutions, the SIL and pSIL solutions do not agree as well with the reference solution for the ccnb problem in the x-direction in plot (b), whereas the MASIL and pMASIL solutions do, thereby suggesting a better selection of eigenvectors representing the solution.

Conclusion.
We have shown that if the Lanczos process is applied to K −1 M with starting vector equal to the spatial distribution vector b, then the Lanczos polynomials are orthogonal with respect to the inner product induced by a step function φ that coincides with the cumulative mass participation sum, that is, the quantity we are interested in to solve Problems 1.1-1.2.The Lanczos polynomials are also orthogonal with respect to an inner product induced by a step function τ k which, unlike φ, is readily available at step k of the Lanczos process.The step function τ k

Fig. 1 .
Fig.1.The eigenvectors (modes of vibration) of a tall building structure.From left to right: the undeformed shape of the structure; its first eigenvector, x 1 corresponding to the smallest eigenvalue, which is predominantly a "sway mode" (i.e., large mass participation factor in the cartesian x-direction); and another eigenvector, x 20 , which is a "bouncing mode"(large participation along cartesian z-direction).Deformations are exaggerated.

Eigenvector x j
Mass participation factor (x T j M b) 2 /(b T M b) b along cartesian x b along cartesian z x

Proposition 3 . 1 .
Let 0 = b be the spatial distribution vector in (1.2) and assume that b ∈ range(K −1 M ).If w = b in Algorithm 2.1 with σ = 0, then for the step function φ(λ) in (3.5) we have that
chilled problem in the x-direction, n = 93445, kmax = 200.Closeup the above stair graph.

Fig. 2 .
Fig. 2.Step functions φ(µ) and τ k (µ) for the chilled problem in the x-direction.The small yellow stars in (a) correspond to the smallest eigenvalues λ i and their normalized mass participation factors m(x i ) such that i m(x i ) ≥ 0.9.The green and red circles correspond to the pairs (λ i , m(x i )) obtained by MASIL and pMASIL, respectively, and such that i m(x i ) ≥ 0.9.The shading shows the intervals selected by the new shifting strategy.

k for step 2 .Fig. 3 .
Fig. 3. Number of intervals s in (3.14) as a function of the number k of unshifted Lanczos steps for the local modes problem in the z-direction.
local modes problem in the z-direction, kmax = 200.chilled problem in the x-direction, kmax = 200.

Fig. 4 .
Fig. 4. Relative error between the reference solution u(t) and its approximation ũ(t) obtained from SIL, MASIL, pSIL, and pMASIL with a proportion of ξ = 90% for the total mass participation.
These are listed in Table2.Our numerical experiments are performed with MATLAB.As mentioned in the introduction, b is always associated with a direction, either x, y, or z.So for a given problem, we will have different solutions depending on the chosen direction since they correspond to different spatial distribution vectors.