On Solving Maxwell Eigenvalue Problems for Accelerating Cavities

We investigate algorithms for computing steady state electromagnetic waves in cavities. The Maxwell equations for the strength of the electric ﬁeld are solved by a mixed method with quadratic ﬁnite edge elements for the ﬁeld values and corresponding node-based ﬁnite elements for the Lagrange multiplier. This approach avoids so-called spurious modes which are introduced if the divergence-free condition for the electric ﬁeld is not treated properly. To compute a few of the smallest positive eigenvalues and corresponding eigenmodes of the resulting large sparse matrix eigenvalue problems two algorithms habe been used, the implicitly restarted Lanczos and Jacobi-Davidson algorithm, both with shift-and-invert spectral transformation. Two-level hierarchical basis preconditioners have been employed for the iterative solution of the resulting systems of equations.


I. INTRODUCTION
The common way to produce the accelerating electromagnetic fields in cyclic accelerators is to excite standing waves in accelerating cavities.The mathematical model for these high frequency electromagnetic fields is the eigenvalue problem solving the Maxwell equations in a bounded volume [1].
Usually, the eigenfield corresponding to the fundamental mode of the cavity is used as the accelerating field.A few modes of higher order have to be analyzed as well because these modes can be excited due to higher harmonic components contained in the RF (radio frequency) power fed into the cavity and through interactions between the accelerated particles and the electromagnetic field.The RF engineer designing such an accelerating cavity therefore needs a tool to compute the fundamental and about ten to twenty of the following eigenfrequencies together with the corresponding electromagnetic eigenfields.The most interesting quantities besides the eigenfrequencies are local maxima of the eigenmodes as well as the fields on the surface that induce heat in the metallic boundary and determine the power loss from surface currents.
The finite element formulation for the solution of this † Electronic address: arbenz@inf.ethz.ch;URL: http://www.inf.ethz.ch/~arbenz/‡ Electronic address: geus@inf.ethz.ch;URL: http:/www.inf.ethz.ch/~geus/§ Electronic address: Stefan.Adam@psi.chproblem has been chosen for its high flexibility for the geometric modeling, in particular with strong variations of the scale of structural details.
For the solution of Maxwell eigenvalue problems Kikuchi [2] suggested a mixed formulation method based on Nédélec's edge elements [3,4] for the electric field and node-based elements for the Lagrange multiplier.In this way, the computed approximate eigenmodes are naturally split into (physically meaningless) curl-free and (in a discrete sense) divergence-free modes.Spurious modes cannot appear.In the engineering literature the advantages and disadvantages of Nédélec's edge elements vs. other finite element formulations are discussed intensely [5,6].
By numerical experiments reported in [7] we identified Sorensen's implicitly restarted Lanczos algorithm [8] and the Jacobi-Davidson algorithm [9] to be the algorithms best suited to solve the large sparse eigenvalue problems that arise with the above finite element discretizations.They are to be preferred to subspace iteration [10] or the block Lanczos algorithm [11].
In this paper we report on our continuing investigations on eigensolvers and their preconditioning.In contrast to [7] we conduct our experiments on a domain that approximates the new cavity design for the 590 MeV ring cyclotron at the Paul Scherrer Institute (PSI) in Villigen, Switzerland, see Fig. 1 on the right-hand side.
The outline of the paper is as follows.In section II we will give the mathematical formulation of the problem.In section III the Nédélec finite element approach for tetrahedral meshes is looked at from a linear algebra point of view.In sections IV to VI we detail how FIG.1: Cross section of two different designs for accelerating cavities.On the left the original design that is presently in operation; on the right the future design.The plane of these two cross sections is perpendicular to the accelerator mid plane and tangential to the orbit of the accelerated particles.Both cavities extend prismatically in the radial direction over approximately 3.3 meters.
we solve the large eigenvalue problems and in particular solve the systems of equations that are introduced by the spectral transformation.In section VII we compare the performance of several variants in solving the model problem.

II. STATEMENT OF THE PROBLEM
Without changing the basic structure of the problem one can assume that the metallic surfaces are perfectly conducting and that the inside volume of the cavity Ω is all in vacuum.The electromagnetic field in the cavity is described by the Maxwell equations [1].After separation of time and space and after elimination of the magnetic field intensity the differential equations are obtained for the electric field intensity e. Equations (1) have solutions only for certain values ω called eigenfrequencies of the electromagnetic oscillations of the resonator.These solutions are called eigenmodes.
The eigenvalue problem (1) can be solved analytically for a few particular domains, e.g. when Ω is a rectangular box or a right circular cylinder [1, p.355].We intend to solve (1) in more complicated domains Ω by the finite element method.To that end we have to write the eigenvalue problem in variational form.There are a number of ways to do so.The mixed formulation proposed by Kikuchi [2] Find (λ, e, p) ∈ R × H 0 (curl; Ω) × H 1 0 (Ω) such that e = 0 and (a) (curl e, curl Ψ) is useful as it leads to a mixed finite element formulation that is free of spurious modes [5,6].Here, (u, v) = Ω u(x)•v(x) dx denotes the usual inner product in L 2 (Ω) 3 , the space of 3-vector valued square integrable functions on Ω.The well-known Sobolev space [12, p.16] 3 , q = 0 on Γ} and the Hilbert space [13,14] It is easy to see, that p = 0 for all solutions of (2), see [2], and that there is a one-to-one correspondence between the solutions (λ, e) of ( 1) and the solutions (λ, e, 0) of ( 2).
In the next section we investigate the finite element discretization of H 0 (curl; Ω) and H 1 0 (Ω) by Nédélec edge elements and Lagrange node elements, respectively.If also H 0 (curl; Ω) is approximated by Lagrange elements the eigenvalue zero will not be confined in the neighborhood of zero and can cause spectral pollution.

III. FINITE ELEMENT DISCRETIZATIONS
In this section we discuss finite element approximations of the mixed formulation (2).We assume that Ω is an open, bounded domain in R 3 with a polyhedral boundary Γ.The domain is triangulated by tetrahedrons.
We replace H 0 (curl; Ω) and H 1 0 (Ω) by finite dimensional subspaces V h and W h , respectively, to obtain The space W h ⊂ H 1 0 (Ω) is the well-known Lagrange (or node-based) finite element space [12], where P k (Ω e ) is the space of polynomials in Ω e of degree ≤ k.In this paper we exclusivly deal with W h .We work with a hierarchical basis and write The natural basis of W (1) h is formed by the piecewise linear polynomials that are 1 in one node (vertex) k and 0 in all others.Thus, h equals the number of interior vertices of the triangulation.The basis of W h consists of those piecewise quadratic polynomials that are 1 in the midpoint of one edge and 0 in the midpoints of all other edges and in all vertices.So, h is equal to the number of interior edges of the triangulation.
The space V h is made up of the H(curl; Ω)-conforming finite elements proposed by Nédélec [4,14,15].The basis of V (1) h consists of piecewise linear vector valued functions with tangential component equal 1 on one interior edge and equal zero on all ther edges.Therefore, h consists of two sets of functions.The first set is made up of piecewise quadratic polynomials that have a linearly varying tangential component along one interior edge and are equal to zero on all other edges.Each of these basis functions is curl-free and equals the gradient of the basis functions of W (2) h corresponding to the same edge!The second set of basis functions of V h is made up of piecewise quadratic polynomials that have their degrees of freedom on the element faces.There are two basis functions per face.They are defined by the two zeroth moments of the tangential components on the face.So, h is equal to the number of interior edges plus twice the number of interior faces of the triangulation.
h .Then (3) defines the matrix eigenvalue problem where A and M are n-by-n and C is n-by-m.We partition the matrices A, M , and C according to the hierarchical bases, h , respectively.Formally, spaces consisting of common Lagrange finite elements could be chosen for V h and W h in (3).However, Nédélec's edge elements have the crucial advantage over Lagrange finite elements that [14, In a way similar to the continuous case [16], we can decompose any -in a discrete sense -into a divergence-free and a curl-free part.We write this decomposition as For We therefore cannot expect monotonic convergence of the eigenvalues.Equation ( 8) corresponds to the 'non-pollution condition' for the conforming finite element approximation of the model problem by Gruber and Rappaz [17, p.20].Therefore, the infinitely multiple eigenvalue zero of curl curl is not spread over the real axis but remains confined in the origin.
In order that or, equivalently, which is the second equation in (6).Because of ( 8) we can write grad ϕ l = n j=1 η jl Φ j .Therefore, In a similar way one obtains H is the system matrix that is obtained when solving the Poisson equation with the Lagrange finite elements is the vector form of the Helmholtz decomposition (9).
Notice that Y is very sparse.We have already mentioned that the gradient of a basis function h is an element of the first set of basis functions of V (2) h .So, m 1 rows of Y have a single entry 1.The gradient of the piecewise linear basis function corresponding to vertex k, say, is a linear combination (with coefficients ±1) of the basis functions of V (1) h whose corresponding edge has vertex k as one of its endpoints.Equation (10) means that C T x = 0 is equivalent to requiring x to be M -orthogonal to the eigenspace N (A) = R(Y ) corresponding to the eigenvalue 0. Thus, the solutions of ( 6) are precisely the eigenpairs of corresponding to the positive eigenvalues.We could therefore compute the desired eigenpairs (λ j , x j ) of ( 6) by means of ( 11) alone.The computed eigenvectors corresponding to positive eigenvalues would automatically satisfy the constraint C T x j = 0.This is actually done if the linear systems of the form (A−σM )x = M y that appear in the eigensolver can be solved by direct methods.Numerical experiments showed however that the highdimensional zero eigenspace has a negative effect on the convergence rates if preconditioned iterative solvers have to be used.

IV. SOLVING THE MATRIX EIGENVALUE PROBLEM
The classic method for computing a few selected eigenpairs of Ax = λM x is subspace iteration [10,18].Its linear convergence rate is however too low to be competitive with modern algorithms.The Lanczos algorithm in turn converges exponentially [10,11].Its main drawback is that it requires memory space proportional to the number of iteration steps that have to be made until convergence.These huge memory requirements can cause the Lanczos algorithm to stop prematurely before convergence has taken place.These well-known facts have been verified for the rectangular cavity [7].The algorithms that successfully solved all problem sizes were the implicitly restarted Lanczos algorithm [8,19,20], and the Jacobi-Davidson algorithm [9,21,22].Like Lanczos these algorithms converge superlinearly but like subspace iteration their memory requirements are bounded by the order of the problem n times a small multiple of the desired eigenvalues p.
In this study we use the implicitly restarted Lanczos algorithm and the Jacobi-Davidson algorithm as described in [7] to solve the matrix eigenvalue problem (6) that we write in the form Note that M −1 A is symmetric (self-adjoint) with respect to the inner product x, y = x T M y.Admissible vectors x are orthogonal to the nullspace The lowest part of its spectrum is desired.

A Spectral transformations
When computing a few, say p, eigenvalues of (12) closest to a number τ it is essential to apply a spec-tral transformation to get a reasonable speed of convergence [11,23,24,25].With the shift-and-invert spectral transformation (12) becomes The shift σ is chosen close to τ .The spectral transformation leaves the eigenvectors unchanged.The eigenvalues of ( 12) close to the shift σ become the largest absolute of (13).In addition they are relatively well-separated which improves the speed of convergence of the eigensolvers [10].The cost of the improved convergence rate is solving a system of equations of the form in each iteration step of the eigensolver.We assume that A and M are so big that they cannot be factored such that (14) has to be solved iteratively.In the subsequent sections we are going to discuss ways how to do this and, in particular, how to precondition (14).

B Implicitly Restarted Lanczos algorithm (IRL)
To overcome the memory consumption of the Lanczos algorithm it has to be restarted in some way.Sorensen [8] suggested an elegant way to do so.He showed how to extract a Krylov subspace with given spectral properties from the actual trial space (also Krylov subspace).These ideas that apply to the Arnoldi algorithm for nonsymmetric eigenvalue problems as well are implemented in the Fortran subroutine package ARPACK [26].
The Lanczos iteration process with shift-and-invert spectral transformation is executed until j = p+k, where k is some positive integer.We used k = p as suggested in [26].Complete reorthogonalization to guarantee stability is feasible as by assumption p + k is not big.If j = p + k then the dimension of the the trial space is reduced to p with a procedure related to the implicit shifted QR algorithm.The shifts are chosen to be the k eigenvalues of the tridiagonal matrix produced by the Lanczos iteration process that are furthest away from the target τ .As the p dimensional subspace left after the above procedure is again a Krylov subspace the Lanczos iteration process can be restarted.
In IRL, the system of equations ( 14) has to be solved to high accuracy as otherwise the Lanczos relation among the basis vectors does not hold true.
Besides the storage for the matrices A and M , the memory requirements of IRL are essentially the space needed to store 2n(p + k) floating point variables.

C Jacobi-Davidson algorithm (JD)
The Jacobi-Davidson algorithm has been introduced by Sleijpen and van der Vorst [9] for solving arbitrary generalized eigenvalue problems.We have adapted their algorithm JDQZ to the generalized eigenvalue problem with symmetric A and symmetric positive definite M [7].
The Jacobi-Davidson algorithm is not a Krylov subspace method.The trial space spanned by the Morthogonal columns of V j ∈ R n×j used in JD is extended by means of the solution t of the correction equation where Qk = [q 1 , . . ., q k , q] and Wk = M Qk .Here, ( λ, q) is a Ritz pair in the actual trial space and q 1 , . . ., q k are the already computed eigenvectors.Initially, k = 0. Notice, that (I− Wk QT k )r = r such that ( 15) is consistent if A − λM is invertible.The solution t of ( 15) is Morthogonalized against the columns of V j and appended to V j to obtain V j+1 .
If j = j max the dimension of the trial space is reduced to j = j min ≥ p by choosing the j min most desirable Ritz vectors from R(V jmax ).
If ( 15) is solved exactly, then one step of the algorithm turns out to be one step of the Rayleigh quotient iteration [9] which converges cubicly for symmetric eigenvalue problems [10].If (15) is solved only approximately this high convergence rate gets lost.There is a trade-off between speed of convergence and the amount of work one is willing to spend for solving (15).
If equation ( 15) is solved approximately by a Krylov subspace method a good preconditioner has to be employed to get a satisfactory convergence rates.Notice that the matrix in (15) maps on the orthogonal complement of Qk instead of Wk .Therefore, the preconditioner also has the duty to map the solutions in the correct subspace.In a similar way as Fokkema et al. [22] we use a preconditioner of the form where the symmetric matrix K is a good and easily invertible approximation of A − λM .It is straightforward to show that the inverse of this mapping is given by where Ỹk = K −1 Wk [16].The preconditioned (15) equation then becomes Memory space is needed for the blocks Qk , Wk , Ỹk , and V .Thus, besides the storage for the matrices A − σM and M , (3p + j max )n memory locations are needed.

V. SOLVING THE CONSTRAINED SYSTEM OF EQUATION I: THE STRAIGHTFORWARD APPROACH
We investigate two ways how to solve (14).In the straightforward approach of this section we apply the conjugate gradient method forcing the iterates into R(C) ⊥ = R(Y ) ⊥ M .In the second approach of section VI we treat ( 14) as an augmented system of order n + m.
We exploit that is symmetric positive definite and apply the preconditioned conjugate gradient (cg) method to solve M −1 (A− σM )x = y.Applying the preconditioner M −1 amounts to approximately solving Notice that both M and A − σM map R(Y ) ⊥ M one-toone onto R(Y ) ⊥ .Two problems arise if ( 19) is solved as it stands.First, (A − σM ) as a map of the whole n-space is indefinite and is bound to slow down convergence [27].Second, if (19) is not solved to high accuracy the approximate solution will fail to satisfy the constraints in (14).By choosing the preconditioner the cg iterates stay in the correct space.In the Jacobi-Davidson algorithm the approximate solution t of ( 15) is additionally forced to be orthogonal to the Ritz vector q and to the previously computed eigenvectors.
Here, we are mainly interested in two-level preconditioners induced by hierarchical finite element basis.Proceeding as Bank [28] we approximate A, M , and C partitioned as in (7) by the block diagonal A 1 , M 1 , and C 1 have similar properties as A, M , and Initialization: (i) r 0 := M y − (A − σM )x 0 , (C T y = 0 and, in general, (iii) p 0 := z 0 , j := 0.
Notice that the M -inner products have disappeared.Operations with M occur in the computation of the residuals and implicitely in the M -orthogonal projector.In the algorithm in Fig. 2 the various iterates satisfy Y T r j = C T z j = C T p j = C T x j = 0 for all j, such that the approximation x j is in the correct space whenever we stop the iteration.Step (vii) in Fig. 2 consumes most of the execution time.Two system of equations have to be solved (1) with A 1 − σ 1 M 1 and ( 2) with H inside the projector.In fact, we use a direct solver for the (1, 1) block A 11 − σ 1 M 11 and just execute one step of the Jacobi or symmetric Gauss-Seidel iteration for approximating the inverse of A 22 − σ 1 M 22 , see Bank [28].For solving Poisson's equation with H we use a direct method as well.

VI. SOLVING THE CONSTRAINED SYSTEM OF EQUATION II: THE AUGMENTED SYSTEM APPROACH
An alternative way to write ( 14) is as an augmented linear system, If A − σM is nonsingular then the block LU factorization of the 2-by-2 block matrix in (23) immediately shows that x = 0. Therefore, the first component x of the solution x of ( 23) satisfies (A−σM )x = M y by the first equation and the constraints by the second equation.
The symmetric indefinite matrix in equation ( 23) has order n+m.It has n positive and m negative eigenvalues.Solving (23) instead of ( 19) has the advantage that there are no (explicit) constraints to be taken into account.There are however more unknowns and quite an accurate solution is needed to satisfy the constraints.
We considered preconditioners of the form when solving (23) by the preconditioned conjugate gradient method [16].A 1 and M 1 have been defined in (20).It turned out that this preconditioner leads to an algorithm that is very similar to the one of the straightforward approach.The difference is that the projector We therefore tried out preconditioners of the form We can easily transcribe the algorithm of Fig. 2 into a solver for (23).Each vector in the algorithm has now two components, the first being in R n the second in R m , see Fig. 3.If x 0 is chosen such that C T x 0 = 0 then traversing stepwise through the algorithm in Fig. 3 one sees that Initialization: Iteration: until convergence ( r j < tol r 0 ) do FIG. 3: Preconditioned conjugate gradient method for the augmented linear system provided that σ = σ 1 which we assume henceforth.To verify (26) we observe that initially C T y = C T x 0 = 0 such that Y T r 0 = 0 and r 0 = 0. Let us assume that (26) holds for all indices up to j.From step (v) in Fig. 3 we see that x j+1 = 0. From step (vi) we have The second component of the solution of the equation in step (vii) of the algorithm is Here, we have used that K 1 = A 1 − σM 1 is nonsingular and that K 1 Y = −σC 1 .p j+1 = 0 follows from step (ix) in Fig. 3.
The relations in (26) simplify the algorithm considerably.Only the r j have a nonzero second component.It is used in step (vii) of the algorithm to compute z j , in fact, z j as z j = 0 as we have just seen.From (27) we get So, actually, r j is never used and thus does not have to be computed at all.It can be recovered at any time as r j = (−1/σ)Y T r j .Taking all the above simplifications into account the algorithm in Fig. 3 with preconditioner (25) becomes Algorithm 4. Initialization: p 0 := z 0 , j := 0.
FIG. 4: Preconditioned conjugate gradient method for the augmented linear system after simplifications Notice that the algorithm in Fig. 4 differs from the one in Fig. 2 in that the projectors I −Y H −1 C T are omitted.In fact, in Fig. 4 plain pcg is applied to the consistent but indefinite system of equations (A−σM )x = M y.It is not yet clear to us if this algorithm is stable at all.Certainly, the constraints C T x j+1 = 0 are not satisfied before convergence.If we solve systems as in the Jacobi-Davidson algorithm with only low accuracy then we enforce the constraints explicitely at the end of the iteration.

VII. NUMERICAL EXPERIMENTS
We performed a few numerical experiments computing the 10 lowest eigenvalues and corresponding eigenvectors of a model of the future design of the copper cavity, cf.Fig. 1.We used three different discretizations, a coarse discretization with only n = 1880 degrees of freedom, a fine discretization with 27'088 and an even finer discretization with 45'040 degrees of freedom.The corresponding values for n 1 and m 1 , cf. (7), are given in Tables I and II.We combined the implicitely restarted Lanczos algorithm (IRL) and Jacobi-Davidson algorithm (JD) with the straightforward ('projection') or modified augmented approach for solving (14).With IRL we solved the positive definite system of equations in the straightforward approach with the preconditioned conjugate gradient method.With JD the shifts vary dynamically and can get inside the convex hull of the spectrum.That's why we used MINRES as our system solver [29].
The system of equations in the augmented system where the preconditioner is indefinite were solved with a special variant of QMR [30].In either case, we set σ 1 = σ to be a value slightly below the lowest positive eigenvalue.When solving (14), we employed no, diagonal or hierarchical basis preconditioning.The last was combined with two different treatments of the A 22 block.Either we performed one step of damped Jacobi (damping factor 2/3) or one step of symmetric Gauss-Seidel.The results for IRL are summarized in Tab.I those for JD in Tab.II.The number it int provides the average number of inner iterations, i.e., the steps until convergence of the linear solver.As systems are solved to much higher accuracy with IRL (to establish the Lanczos three-term recurrence) the iteration count it int is much bigger in Tab.I than in Tab.II.
In fact, with JD, the convergence tolerance of the inner iteration is initially very loose and becomes more stringent as the outer iteration converges [7].In contrast, the number of outer iterations it out , i.e., the number of systems that are solved is smaller with IRL.In both IRL and JD we let the dimension of the trial space vary between 10 and 25. t tot gives the overall solution time, including grid handling, matrix assembly, eigensolver, etc., while t eig is the time spent for solving the eigensystem only.Finally, t 11 gives the time that is spent for solving (forward/backward substitution) the upper diagonal (1, 1) block by a direct solver.t eig includes t 11 .These two numbers are the times of interest.Our discussion is based on them.Notice that in the JD timings the overall solution time for the augmented approach does not contain the building of matrix C and H while it is still included in the IRL times.
The numerical experiments were run in sequential mode on a Sun Enterprise 3500 with six 336 Mhz Ul-traSparc processors and 3 GB main memory running the Solaris 2.6 operating system.All times are user time.
In all examples, the Jacobi-Davidson algorithm is the faster eigensolver than the implicitly restarted Lanczos algorithm.IRL executes many more steps in the inner iteration as the Lanczos vectors have to be computed to quite a high precision [26].
Higher sophistication in the preconditioner decreases not only the iteration counts but also the execution times.This becomes more prominent as the problem size increases.This is in contrast to the experiments that we conducted with node elements in the box shaped cavity [31].The analysis of the hierarchical basis preconditioner states that it int should not increase with the problem size.Here it even decreases substantially.
When comparing the straightforward with the (modified) augmented system approach, IRL and JD behave quite differently.With IRL we observe a big gain of more that 25% when employing the augmented system approach combined with damped Jacobi smoothing.The gain is almost negligible if Jacobi is replaced by symmetric Gauss-Seidel.The average number of inner iterations differs only slightly in the two approaches.With Jacobi-Davidson there is a minor gain with damped Jacobi smoothing and a considerable loss with symmetric Gauss-Seidel smoothing.The average number of inner iterations is twice as high for the augmented system approach.The reduced cost per iteration by saving the projector is compensated by a higher number of iterations.
It is worth noting that the time t 11 for solving the (1,1) block takes a considerable fraction of the solution of the eigenvalue problem.t 11 is about 20% of t eig in the medium sized problems and about 25% in the large problems.With the augmented system of equation this portion is even above 50%.In this case, the direct solver appears to not yet having been implemented optimally.In any case, the present two-level preconditioner will not suffice for larger problem sizes.A multi-level method will be necessary for solving the system of equations corresponding to the (1,1) block A 11 − σM 11 .

VIII. CONCLUSIONS
We have investigated algorithms for solving the eigenvalue problem occurring in the design of resonant cavities.The governing Maxwell equations for the electric field intensity have been discretized by quadratic Nédélec edge elements.The finite element spaces had dimensions up to 45'040.The resulting symmetric matrix eigenvalue problems have been solved by two methods, the Jacobi-Davidson (JD) and the implicitly restarted Lanczos (IRL) algorithm.In both methods, a shift-and-invert spectral transformation has been applied.The such introduced symmetric indefinite systems of equations have been solved iteratively by a Krylov subspace method.The most effective preconditioners with respect to time and iteration numbers were hierarchical basis preconditioners.The average number of iterations to solve these systems not only stayed constant but decreased with increasing problem size.
The JD algorithm clearly outperformed IRL.The execution times of the latter are about twice as long due to the high accuracy requirement of the construction of a correct three-term recurrence.
The symmetric indefinite systems of equations have been solved in two ways.First, by projecting in the subspace where the system is positive definite, second, by augmenting the system, i.e., by introducing constraints and Lagrange multipliers.With the first we obtained very satisfactory results.The latter gave good results with IRL but lagged clearly behind with JD.The augmented systems approach is not yet completely understood with regard to its convergence behavior and stability.Further investigations in this direction are required.

TABLE I :
Implicitly Restarted Lanczos (ARPACK) results for the copper cavity.Times are in seconds.