Roland W Freund

Krylov-subspace methods for reduced-order modeling in circuit simulation

Keywords: Lanczos algorithm, Arnoldi process, Linear dynamical system, VLSI interconnect, Transfer function, Padà e approximation, Stability, Passivity, Positive real function

The simulation of electronic circuits involves the numerical solution of very large-scale, sparse, in general nonlinear, systems of di erential-algebraic equations. Often, the size of these systems can be reduced considerably by replacing the equations corresponding to linear subcircuits by approximate models of much smaller state-space dimension. In this paper, we describe the use of Krylov-subspace methods for generating such reduced-order models of linear subcircuits. Particular emphasis is on reduced-order modeling techniques that preserve the passivity of linear RLC subcircuits.

Introduction

Today's integrated electronic circuits are extremely complex, with up to tens of millions of devices. Prototyping of such circuits is no longer possible, and instead, computational methods are used to simulate and analyze the behavior of the electronic circuit at the design stage. This allows to correct the design before the circuit is actually fabricated in silicon.

The simulation of electronic circuits involves the numerical solution of very large-scale, sparse, in general nonlinear, systems of time-dependent di erential-algebraic equations (DAEs); see, e.g. [START_REF] Feldmann | Numerical Simulation of Electronic Circuits: State-of-the-Art Techniques and Challenges[END_REF][START_REF] Sangiovanni-Vincentelli | Circuit simulation[END_REF][START_REF] Vlach | Computer Methods for Circuit Analysis and Design[END_REF] and the references given there. These systems can be so large that time integration becomes ine cient or even prohibitive. On the other hand, electronic circuits often contain large linear subcircuits of passive components that contribute only linear equations to the system of DAEs describing the whole circuit. In particular, such linear subcircuits may result from extracted RLC models of the circuit's wiring, the so-called interconnect, models of the circuit's package, or models of wireless propagation channels. By replacing the equations corresponding to linear subcircuits by approximate models of much smaller state-space dimension, the size of the system of DAEs describing the whole circuit can be reduced signiÿcantly, so that time integration of the resulting system becomes feasible; see, e.g. [START_REF] Feldmann | Interconnect Extraction and Analysis in High-Frequency, Sub-Micron, Digital VLSI Design[END_REF][START_REF] Kim | Time-domain macromodels for VLSI interconnect analysis[END_REF][START_REF] Pileggi | Coping with RC(L) interconnect design headaches[END_REF][START_REF] Raghavan | AWESpice: A general tool for the accurate and e cient simulation of interconnect problems[END_REF] and the references given there. In recent years, there has been a lot of interest in generating suitable reduced-order models of linear subcircuits by means of Krylov-subspace methods, such as the Lanczos algorithm and the Arnoldi process. For a survey of these recent developments, we refer the reader to [START_REF] Freund | Reduced-order modeling techniques based on Krylov subspaces and their use in circuit simulation[END_REF].

In this paper, we describe the use of Krylov-subspace methods for generating reduced-order models of systems of linear DAEs, such as the ones arising in circuit simulation. Particular emphasis is on projection techniques that, when applied to a passive circuit, preserve the passivity of the circuit. We stress that the methods discussed in this paper are not restricted to systems of DAEs arising in circuit simulation and that they can be applied to general time-invariant linear dynamical systems. However, the development of these methods was mostly motivated by the need for reduced-order modeling in circuit simulation.

The remainder of the paper is organized as follows. In Section 2, we brie y review the systems of DAEs that arise in circuit simulation, and we describe how reduced-order models of linear subcircuits can be employed to reduce the dimension of these systems. In Section 3, we introduce our notion of block Krylov subspaces and review the construction of basis vectors via Lanczos and Arnoldi algorithms. In Section 4, we deÿne reduced-order models based on projection and describe their computation via Krylov-subspace methods. In Section 5, we discuss connections with Padà e and Padà e-type approximants. In Section 6, we establish results on the stability and passivity of reduced-order models obtained via projection. In Section 7, numerical results for two circuit examples are reported. Finally, in Section 8, we make some concluding remarks and mention a few open problems.

Throughout this article, we use boldface letters to denote vectors and matrices. Unless stated otherwise, vectors and matrices are allowed to have complex entries. As usual, M = [m jk ]; M T = [m kj ], and M H = M T = [m kj ] denote the complex conjugate, transpose, and the conjugate transpose, respectively, of the matrix M = [m jk ], and M ¿0 means that M is Hermitian positive semi-deÿnite. The vector norm x := √ x H x is always the Euclidean norm, and M := max x =1 Mx is the corresponding induced matrix norm. We use I n to denote the n × n identity matrix and 0 n×m to denote the n × m zero matrix; we will omit these indices whenever the actual dimensions of I and 0 are apparent from the context. The sets of real and complex numbers are denoted by R and C, respectively. For s ∈ C; Re(s) is the real part of s. Finally, C + := {s ∈ C | Re(s) ¿ 0} is the open right-half of the complex plane.

Circuit equations

In this section, we brie y describe the systems of DAEs that arise in circuit simulation and review how reduced-order modeling of linear subcircuits is employed in the numerical solution of such systems. For introductions to circuit simulation and overviews of typical simulation tasks, we refer the reader to [START_REF] Feldmann | Numerical Simulation of Electronic Circuits: State-of-the-Art Techniques and Challenges[END_REF][START_REF] Sangiovanni-Vincentelli | Circuit simulation[END_REF][START_REF] Vlach | Computer Methods for Circuit Analysis and Design[END_REF].

General circuit equations

Electronic circuits are usually modeled as networks whose branches correspond to the circuit elements and whose nodes correspond to the interconnections of the circuit elements; see, e.g. [START_REF] Feldmann | Numerical Simulation of Electronic Circuits: State-of-the-Art Techniques and Challenges[END_REF][START_REF] Sangiovanni-Vincentelli | Circuit simulation[END_REF][START_REF] Vlach | Computer Methods for Circuit Analysis and Design[END_REF]. Such networks are characterized by three types of equations: Kirchho 's current law (KCL), Kirchho 's voltage law (KVL), and branch constitutive relations (BCRs). The unknowns in these equations are the currents through the branches of the network, the voltage drops along the branches, and the voltages at the nodes of the network. The KCLs and KVLs are linear algebraic equations that only depend on the topology of the circuit. The KCLs state that, at each node N of the network, the currents owing in and out of N sum up to zero. The KVLs state that, for each closed loop L of the network, the voltage drops along L sum up to zero. The BCRs are equations that characterize the actual circuit elements. For example, the BCR of a linear resistor is Ohm's law. The BCRs are linear equations for simple devices, such as linear resistors, capacitors, and inductors, and they are nonlinear equations for more complex devices, such as diodes and transistors. Furthermore, in general, the BCRs involve ÿrst time derivatives of the unknowns, and thus they are ÿrst-order DAEs.

All the KCLs, KVLs, and BCRs characterizing a given circuit can be summarized as a system of ÿrst-order, in general nonlinear, DAEs of the form

f ( x; t) + d dt q( x; t) = 0; (1)
together with suitable initial conditions. Here, x = x(t) is the unknown vector of circuit variables at time t, the vector-valued function f ( x; t) represents the contributions of nonreactive elements such as resistors, sources, etc., and the vector-valued function (d=dt)q( x; t) represents the contributions of reactive elements such as capacitors and inductors. There are a number of established methods, such as sparse tableau, nodal formulation, and modiÿed nodal analysis, for writing down the system (1); see, e.g. [START_REF] Vlach | Computer Methods for Circuit Analysis and Design[END_REF]. The vector functions x; f ; q in (1), as well as their dimension, N , depend on the chosen formulation method. The most general method is sparse tableau, which consists of just listing all the KCLs, KVLs, and BCRs. The other formulation methods can be interpreted as starting from sparse tableau and eliminating some of the unknowns by using some of the KCL or KVL equations. For all the standard formulation methods, the dimension N is of the order of the number of devices in the circuit.

Linear subcircuits

Traditional circuit simulators are based on the numerical solution of the system of DAEs (1); see, e.g. [START_REF] Vlach | Computer Methods for Circuit Analysis and Design[END_REF]. However, the dimension of (1) can be so large that time integration of (1) is ine cient or even prohibitive. On the other hand, circuits often contain large linear subcircuits that can be well approximated by reduced-order models of much smaller dimension. By replacing the equations in (1) corresponding to such linear subcircuits by their respective reduced-order models, one obtains an approximate system of DAEs of much smaller dimension that can then be solved numerically by time integration. We now describe this process in more detail.

Let C l be a large linear subcircuit of a given circuit, and denote by C r the, in general nonlinear, remainder of the circuit. After a suitable reordering, the vector x of circuit variables in (1) can be partitioned as follows:

x =   xr y xl   : (2)
Here, xr and xl denote the circuit variables exclusive to C r and C l , respectively, and y represents the variables shared by C r and C l . Using the partitioning (2) and setting

x 0 := xr y and x := y xl ;

(3) the functions f and q in (1), after a suitable reordering of the equations in [START_REF] Aliaga | A Lanczos-type method for multiple starting vectors[END_REF], can be expressed as follows:

f ( x; t) = f 0 (x 0 ; t) 0 k×1 + 0 N -N ×1 Gx ; q( x; t) = q 0 (x 0 ; t) 0 k×1 + 0 N -N ×1 Cx : (4)
Here, f 0 and q 0 represent the contributions of resistive and reactive elements from the subcircuit C r , and the matrices G and C represent the contributions of resistive and reactive elements in the linear subcircuit C l . In (4), without loss of generality, we have assumed that the vector-valued functions f 0 and q 0 have the same number of components, that the zero vectors below f 0 and q 0 have the same length, k, and that the matrices G and C are square and of the same size, N × N ; this can always be achieved by padding f 0 ; q 0 ; G , and C with additional zeros, if necessary. Unless the subcircuit C l is completely decoupled form the remainder circuit C r , we have m := N -k ¿ 0. This means that, in (4), the last m components of the, in general nonlinear, functions f 0 and q 0 are connected with the ÿrst m components of the linear functions Gx and Cx. By introducing an additional m-dimensional vector, u = u(t), of circuit variables, these m connecting equations can be decoupled. Indeed, using (4), one readily veriÿes that the original system (1) is equivalent to the following system:

f 0 (x 0 ; t) + d dt q 0 (x 0 ; t) + 0 I m u = 0; (5)
C dx dt + Gx = I m 0 u: (6)
We remark that the additional variables u in ( 5) and ( 6) can be interpreted as interface signals between the subcircuits C r and C l . Let p denote the length of the vector y in the partitioning (2) of x, and set

B := I m 0 N -m×m
and L := I p 0 N -p×p :

Note that, by (3), the matrix L H = L T selects the subvector y from x, i.e.

y = L H x: (7)
Eqs. ( 6) and ( 7) constitute a linear dynamical system of the form

C dx dt = -Gx + Bu(t); y(t) = L H x(t): (8)
In [START_REF] Feldmann | E cient linear circuit analysis by Padà e approximation via the Lanczos process[END_REF], in general, C ; G ∈ C N ×N ; B ∈ C N ×m , and L ∈ C N ×p are given matrices, m and p denote the number of inputs and outputs, respectively, the components of the given vector-valued function u : [0; ∞) → C m are the inputs, and y : [0; ∞) → C p is the unknown function of outputs. The components of the unknown vector-valued function x : [0; ∞) → C N are the state variables, and N is the state-space dimension. In general, the matrices C and G in ( 8) are allowed to be singular. However, we assume that G +sC is a regular matrix pencil, i.e., G +sC is singular only for ÿnitely many values of s ∈ C. This condition is always satisÿed for linear dynamical systems [START_REF] Feldmann | E cient linear circuit analysis by Padà e approximation via the Lanczos process[END_REF] arising in circuit simulation.

A reduced-order model of ( 8) is a linear dynamical system of the same form as ( 8), but of smaller state-space dimension n ¡ N . More precisely, a reduced-order model of state-space dimension n is of the form

C n dz dt = -G n z + B n u(t); y(t) = L H n z(t): (9)
Here, C n ; G n ∈ C n×n ; B n ∈ C n×m , and L n ∈ C n×p are matrices that should be chosen such that the input-output mapping u(t) → y(t) of ( 9) somehow approximates the input-output mapping of the original system [START_REF] Feldmann | E cient linear circuit analysis by Padà e approximation via the Lanczos process[END_REF]; see Section 2.3 below. After a suitable reduced-order model ( 9) for systems ( 6) and ( 7) has been determined, the linear part (6) of the circuit equations is replaced by the ÿrst set of equations in [START_REF] Feldmann | Reduced-order modeling of large linear subcircuits via a block Lanczos algorithm[END_REF]. The result is a reduced-order system of DAEs that represents an approximation to the original system (1); see, e.g. [START_REF] Feldmann | Reduced-order modeling of large linear subcircuits via a block Lanczos algorithm[END_REF][START_REF] Raghavan | AWESpice: A general tool for the accurate and e cient simulation of interconnect problems[END_REF]. Provided that the size of C l dominates that of C r , the approximate system has a much smaller state-space dimension than (1), and thus time integration by means of standard circuit simulators becomes feasible.

Transfer functions

Next, we introduce the so-called transfer function, which describes the input-output behavior of a linear dynamical system (9) in frequency domain.

For vector-valued functions g(t); t ∈ [0; ∞), with g(0) = 0, we denote by

ĝ(s) = ∞ 0 g(t)e -st dt; s ∈ C; (10)
the (frequency-domain) Laplace transform of g. We remark that in [START_REF] Feldmann | Numerical Simulation of Electronic Circuits: State-of-the-Art Techniques and Challenges[END_REF], the purely imaginary values s = i!; !¿0, correspond to the frequency !; these are the physically meaningful values of the complex variable s. We now assume, for simplicity, zero initial conditions x(0) = 0 and u(0) = 0 in [START_REF] Feldmann | E cient linear circuit analysis by Padà e approximation via the Lanczos process[END_REF]. By applying [START_REF] Feldmann | Numerical Simulation of Electronic Circuits: State-of-the-Art Techniques and Challenges[END_REF] to the linear dynamical system (8), we obtain its frequency-domain formulation

sC x = -G x + B û(s); ŷ(s) = L H x(s): (11)
Eliminating x in [START_REF] Feldmann | Interconnect Extraction and Analysis in High-Frequency, Sub-Micron, Digital VLSI Design[END_REF] results in the frequency-domain input-output relation ŷ(s) = H (s) û(s), where H , the transfer function of (8), is given by

H (s) := L H (G + sC ) -1 B; s ∈ C: (12)
Note that H : C → (C ∪ {∞}) p×m is a matrix-valued rational function.

Similarly, the transfer function H n : C → (C ∪ {∞}) p×m of the reduced-order model ( 9) is given by

H n (s) := L H n (G n + sC n ) -1 B n ; s ∈ C: (13)
In terms of transfer functions, the problem of constructing a reduced-order model ( 9) of size n that approximates the input-output behavior of ( 8) can be stated as follows: Determine the matrices C n ; G n ; B n , and L n in ( 9) such that the reduced-order transfer function ( 13), H n , in some sense approximates the transfer function ( 12), H , of the original linear dynamical system [START_REF] Feldmann | E cient linear circuit analysis by Padà e approximation via the Lanczos process[END_REF].

For systems [START_REF] Feldmann | E cient linear circuit analysis by Padà e approximation via the Lanczos process[END_REF] of small-to-moderate state-space dimension N , there is a variety of techniques to construct reduced-order models such that, in some appropriate norm, H n approximates H over a whole range of values of s; see the references given in [START_REF] Freund | Reduced-order modeling techniques based on Krylov subspaces and their use in circuit simulation[END_REF]. However, these techniques are usually not applicable to large-scale systems [START_REF] Feldmann | E cient linear circuit analysis by Padà e approximation via the Lanczos process[END_REF], such as the ones arising in circuit simulation. In the latter case, the matrices C and G in ( 8) are large and sparse. Note that, in view of ( 12), the evaluation of H (s 0 ) at even a single point s 0 ∈ C requires the solution of systems of linear equations with the large sparse coe cient matrix G + s 0 C . Fortunately, the circuit matrices C and G are usually such that sparse Gaussian elimination can be employed to compute an LU factorization

G + s 0 C = P 1 L 0 U 0 P 2 ( 14
)
of the matrix G + s 0 C . In ( 14), P 1 and P 2 are permutation matrices that record pivoting for sparsity and numerical stability, L 0 is a lower triangular matrix, and U 0 is an upper triangular matrix. Pivoting for sparsity means that the original ordering of the rows and columns of G + s 0 C is changed so that potential ÿll-in in the factors L 0 and U 0 is reduced. For circuit matrices, typically very little ÿll-in occurs in L 0 and U 0 , although this cannot be guaranteed in general. Once the factorization ( 14) is computed, the solution of the linear systems needed to evaluate H (s 0 ) is obtained by sparse backsolves.

Note that evaluating H (s 0 ) at several points s 0 would require the computation of a new factorization [START_REF] Freund | Passive reduced-order models for interconnect simulation and their computation via Krylov-subspace algorithms[END_REF] for each new point s 0 . Despite the limited ÿll-in for circuit matrices, the cost for factoring G + s 0 C is high enough that one tries to get away with computing a single factorization [START_REF] Freund | Passive reduced-order models for interconnect simulation and their computation via Krylov-subspace algorithms[END_REF]. This is the case for reduced-order models that are characterized by a matching of the leading terms in Taylor expansions of H and H n about a given expansion point s 0 . More precisely, such a reduced-order model of given size n is deÿned by

H n (s) = H (s) + O(s -s 0 ) q(n) : (15)
If q(n) in ( 15) is as large as possible, then H n is an nth matrix-Padà e approximant of H ; see, e.g. [START_REF] Baker | Padà e Approximants[END_REF]. In Section 5, we will also discuss certain matrix-Padà e-type approximants for which q(n) is not maximal.

Linear RLC subcircuits

In circuit simulation, an important special case is linear subcircuits that consist of only resistors, inductors, and capacitors. Such linear RLC subcircuits arise in the modeling of a circuit's interconnect and package; see, e.g. [START_REF] Freund | The SyMPVL algorithm and its applications to interconnect simulation[END_REF][START_REF] Freund | Reduced-order modeling of large linear passive multi-terminal circuits using matrix-Padà e approximation[END_REF][START_REF] Kim | Time-domain macromodels for VLSI interconnect analysis[END_REF][START_REF] Pileggi | Coping with RC(L) interconnect design headaches[END_REF].

The equations describing linear RLC subcircuits are of the form [START_REF] Feldmann | E cient linear circuit analysis by Padà e approximation via the Lanczos process[END_REF]. Furthermore, the equations can be formulated such that the matrices in (8) exhibit certain symmetries; see [START_REF] Freund | Reduced-order modeling of large passive linear circuits by means of the SyPVL algorithm[END_REF][START_REF] Freund | Reduced-order modeling of large linear passive multi-terminal circuits using matrix-Padà e approximation[END_REF]. More precisely, the N × N matrices G and C are real and symmetric, and have the following block structure:

G = G T = G 11 G 12 G T 12 0 and C = C T = C 11 0 0 -C 22 : (16)
Here, the submatrices G 11 ; C 11 ∈ R N1×N1 and C 22 ∈ R N2×N2 are symmetric positive semi-deÿnite, and N = N 1 + N 2 . Note that, except for the special case N 2 = 0, the matrices G and C are indeÿnite. The special case N 2 = 0 arises for RC subcircuits that contain only resistors and capacitors, but no inductors.

If the RLC subcircuit is viewed as an m-terminal component with m = p inputs and outputs, then the matrices B and L in ( 8) are identical and of the form

B = L = B 1 0 N2×m with B 1 ∈ R N1×m : (17)
For such an m-terminal RLC subcircuit, in view of ( 16) and ( 17), the transfer function [START_REF] Freund | Computation of matrix Padà e approximations of transfer functions via a Lanczos-type process[END_REF] reduces to

H (s) = B T (G + sC ) -1 B where G = G T ; C = C T : (18)
We call a transfer function H symmetric if it is of the form [START_REF] Guillemin | Synthesis of Passive Networks[END_REF] with real matrices G , C , and B. For symmetric transfer functions, we will always assume that the expansion point s 0 in ( 15) is chosen to be real:

s 0 ∈ R if H is symmetric: (19)
The condition [START_REF] Kerns | Preservation of passivity during RLC network reduction via split congruence transformations[END_REF] is necessary in order to generate passive reduced-order models of symmetric transfer functions. We will also use the following nonsymmetric formulation of [START_REF] Guillemin | Synthesis of Passive Networks[END_REF]. Let J be the block matrix

J = I N1 0 0 -I N2 : (20)
Note that, by [START_REF] Freund | Reduced-order modeling of large linear passive multi-terminal circuits using matrix-Padà e approximation[END_REF] and [START_REF] Kim | Time-domain macromodels for VLSI interconnect analysis[END_REF], we have B = JB. Using this relation, as well as [START_REF] Freund | The SyMPVL algorithm and its applications to interconnect simulation[END_REF], we can rewrite [START_REF] Guillemin | Synthesis of Passive Networks[END_REF] as follows:

H (s) = B T (JG + sJC ) -1 B; (21)
where

JG = G 11 G 12 -G T 12 0
and JC = C 11 0 0 C 22 :

In this formulation, the matrix JG is no longer symmetric, but now JG + (JG ) T ¿0 and JC ¿0:

Basis vectors for block Krylov subspaces

In this section, we introduce our notion of block Krylov subspaces for multiple starting vectors. We also review variants of the Arnoldi and Lanczos algorithms for generating basis vectors for block Krylov subspaces.

Reduction to one matrix

Let s 0 ∈ C be the expansion point that is to be used in the characterization [START_REF] Freund | Reduced-order modeling of large passive linear circuits by means of the SyPVL algorithm[END_REF] of the reduced-order transfer function H n . The only assumption on s 0 is that the matrix G + s 0 C be nonsingular; this guarantees that s 0 is not a pole of the original transfer function [START_REF] Freund | Computation of matrix Padà e approximations of transfer functions via a Lanczos-type process[END_REF], H .

An approximate transfer function H n satisfying (15) could be obtained by ÿrst explicitly computing the leading q(n) Taylor coe cients of the expansion of H about s 0 and then generating H n from these coe cients; see, e.g. [START_REF] Pillage | Asymptotic waveform evaluation for timing analysis[END_REF]. However, any approach based on explicitly computing the Taylor coe cients of H is inherently numerically unstable; see [START_REF] Feldmann | E cient linear circuit analysis by Padà e approximation via the Lanczos process[END_REF]. A much better alternative is to use block Krylov-subspace methods that obtain the same information as contained in the leading q(n) Taylor coe cients of H , but in a more stable manner.

Before block Krylov-subspace methods can be employed, the two matrices G and C in the deÿnition [START_REF] Freund | Computation of matrix Padà e approximations of transfer functions via a Lanczos-type process[END_REF] of H have to be reduced to a single matrix, denoted by A in the sequel. This can be done by rewriting [START_REF] Freund | Computation of matrix Padà e approximations of transfer functions via a Lanczos-type process[END_REF] as follows:

H (s) = L H (I + (s -s 0 )A) -1 R; (23)
where

A := (G + s 0 C ) -1 C and R := (G + s 0 C ) -1 B:
Although G and C are sparse matrices, in general, the matrix A in ( 23) is a dense matrix. However, block Krylov-subspace methods involve A only in the form of matrix-vector products AC and possibly A H w. To e ciently compute these products, one never needs to form A explicitly. Instead, one uses the sparse factorization ( 14) of G + s 0 C . Each matrix-vector product AC then requires one multiplication with the sparse matrix C and two backsolves with the sparse triangular matrices L 0 and U 0 from [START_REF] Freund | Passive reduced-order models for interconnect simulation and their computation via Krylov-subspace algorithms[END_REF]. Similarly, A H w requires one multiplication with C H and two backsolves with the sparse triangular matrices L H 0 and U H 0 .

Block Krylov subspaces

Next, we introduce block Krylov subspaces. The proper deÿnition of these subspaces is necessarily quite involved, and the reader may ask if block Krylov subspaces could not be avoided altogether by using standard Krylov subspaces induced by single vectors instead. For example, one can generate scalar approximations for all the p • m coe cient functions of the p × m-matrix-valued transfer function H via suitable basis vectors for m + p standard Krylov subspaces. However, the resulting approximation is not a matrix-Padà e approximant of H , and in fact, one can show that, in order to obtain an approximation of the same quality as the matrix-Padà e approximant, at least (m + p)=2 times more computational work is required compared to computing a matrix-Padà e approximant. Therefore, the use of block Krylov subspaces results in much more e cient reduced-order modeling techniques than those based on standard Krylov subspaces.

Let A ∈ C N ×N be a given N × N matrix and

R = [r 1 r 2 • • • r m ] ∈ C N ×m ( 24
)
be a given matrix of m right starting vectors, r 1 ; r 2 ; : : : ; r m . Before we introduce block Krylov subspaces induced by A and R, we brie y review the standard case m = 1 of a single starting vector r = r 1 . In this case, the usual nth Krylov subspace (induced by A and r) is given by K n (A; r) := span{r; Ar; A 2 r; : : : ; A n-1 r}:

Let n 0 be deÿned as the largest possible integer n such that in [START_REF] Pillage | Asymptotic waveform evaluation for timing analysis[END_REF], all the Krylov vectors, A j-1 r, 16j6n -1, are linearly independent. Note that n 0 6N . By the deÿnition of n 0 , the nth Krylov subspace [START_REF] Pillage | Asymptotic waveform evaluation for timing analysis[END_REF] has dimension n if 16n6n 0 and dimension n 0 if n ¿ n 0 . Moreover, K n (A; r) = K n0 (A; r) for all n ¿ n 0 . Thus, K n0 (A; r) is the largest possible Krylov subspace (induced by A and r), and we call the Krylov sequence r; Ar; A 2 r; : : : ; A n-1 r exhausted if n ¿ n 0 .

In the general case of m¿1 starting vectors [START_REF] Pileggi | Coping with RC(L) interconnect design headaches[END_REF], the situation is more complicated; we refer the reader to the discussion in [START_REF] Aliaga | A Lanczos-type method for multiple starting vectors[END_REF]. The main di culty is that in contrast to the case m = 1, linear independence of the columns in the block Krylov sequence, R; AR; A 2 R; : : : ; A j-1 R; : : : ;

is lost only gradually in general. More precisely, if the jth block, A j-1 R, contains a column that is linearly dependent on columns to its left in [START_REF] Raghavan | AWESpice: A general tool for the accurate and e cient simulation of interconnect problems[END_REF], then, in general, not all the columns of A j-1 R are linear dependent on columns to their left. Hence, the occurrence of a linear-dependent column does not necessarily imply that the block Krylov sequence R; AR; A 2 R; : : : ; A j-1 R is exhausted. As a result, in general, the construction of suitable basis vectors for the subspaces spanned by the columns of ( 26) needs to be continued even after a linear-dependent column has been found in [START_REF] Raghavan | AWESpice: A general tool for the accurate and e cient simulation of interconnect problems[END_REF]. However, a proper handling of such linear-dependent columns requires that the column itself and all its successive A-multiples need to be deleted. Formally, this can be done as follows. By scanning the columns of the matrices in ( 26) from left to right and deleting each column that is linearly dependent on earlier columns, we obtain the de ated block Krylov sequence

R 1 ; AR 2 ; A 2 R 3 ; : : : ; A jmax-1 R jmax : (27)
This process of deleting linearly dependent vectors is referred to as exact de ation in the following. In [START_REF] Rohrer | Passivity considerations in stability studies of numerical integration algorithms[END_REF], for each j =1; 2; : : : ; j max , R j is a submatrix of R j-1 , with R j = R j-1 if, and only if, de ation occurred within the jth Krylov block A j-1 R in [START_REF] Raghavan | AWESpice: A general tool for the accurate and e cient simulation of interconnect problems[END_REF]. Here, for j = 1, we set R 0 = R. Denoting by m j the number of columns of R j , we thus have

m¿m 1 ¿m 2 ¿ • • • ¿m jmax ¿1: (28)
By construction, the columns of the matrices (27) are linearly independent, and for each n, the subspace spanned by the ÿrst n of these columns is called the nth block Krylov subspace (induced by A and R). In the following, we denote the nth block Krylov subspace by K n (A; R). For later use, we remark that for

n = m 1 + m 2 + • • • + m j where 16j6j max ; (29)
the nth block Krylov subspace is given by

K n (A; R) = Colspan{R 1 ; AR 2 ; A 2 R 3 ; : : : ; A j-1 R j }: ( 30
)
For Lanczos-based reduced-order modeling techniques, we will also need the block Krylov subspaces induced by A H and a given matrix of p left starting vectors,

L = [l 1 l 2 • • • l p ] ∈ C N ×p : (31)
Applying the above construction to A H and L, the nth block Krylov subspace (induced by A H and L), K n (A H ; L), is deÿned as the subspace spanned by the ÿrst n columns of the de ated block Krylov sequence

L 1 ; A H L 2 ; (A H ) 2 L 3 ; : : : ; (A H ) kmax-1 L kmax : (32)
Denoting by p k the number of columns of L k , we have

p¿p 1 ¿p 2 ¿ • • • ¿p kmax ¿1: (33)
Note that for n = p 1 + p 2 + • • • + p k where 16k6k max ;

we have

K n (A H ; L) = Colspan{L 1 ; A H L 2 ; (A H ) 2 L 3 ; : : : ; (A H ) k-1 L k }:
We stress that in our construction of block Krylov subspaces, we only have used exact de ation of columns that are linearly dependent. In an actual algorithm for constructing basis vectors for K n (A; R) and K n (A H ; L) in ÿnite-precision arithmetic, one also needs to delete vectors that are in some sense "almost" linearly dependent on earlier vectors. We will refer to the deletion of such almost linearly dependent vectors as inexact de ation. In Sections 3.4 and 3.5 below, we describe how exact and inexact de ation can be built easily into Arnoldi-and Lanczos-type algorithms for multiple starting vectors. While inexact de ation is crucial in practice, concise statements of theoretical results are much simpler if only exact de ation is performed. Throughout this paper, theoretical results are thus formulated for exact de ation only.

For later use, we note the following invariance property of the block Krylov subspaces K n (A; R) induced by the matrices A and R deÿned in [START_REF] Odabasioglu | PRIMA: passive reduced-order interconnect macromodeling algorithm[END_REF].

Lemma 1. Let G ; C ; B be the matrix triplet used in the deÿnition of the matrices A and R in [START_REF] Odabasioglu | PRIMA: passive reduced-order interconnect macromodeling algorithm[END_REF]; and let J be any nonsingular matrix of the same size as A. Then the matrix triplets G ; C ; B and JG ; JC ; JB lead to the same nth block Krylov subspace K n (A; R).

Proof. By (23), we have

A = (G + s 0 C ) -1 C = (JG + s 0 JC ) -1 JC ; R = (G + s 0 C ) -1 B = (JG + s 0 JC ) -1 JB:
Thus both matrix triplets result in the same matrices A and R and the associated block Krylov subspaces are identical.

Basis vectors

The columns of the de ated block Krylov sequences ( 27) and (32), which are used in the above deÿnitions of K n (A; R) and K n (A H ; L), respectively, tend to be almost linearly dependent even for moderate values of n. Therefore, they should not be used in numerical computations. Instead, we construct other suitable basis vectors.

In the following,

C 1 ; C 2 ; : : : ; C n ∈ C N (34)
denotes a set of basis vectors for K n (A; R), i.e., span{C 1 ; C 2 ; : : : ;

C n } = K n (A; R):
The N × n matrix

V n := [C 1 C 2 • • • C n ] (35)
whose columns are the basis vectors (34) is called a basis matrix for K n (A; R). Similarly, w 1 ; w 2 ; : : : ;

w n ∈ C N (36)
denotes a set of basis vectors for K n (A H ; L), i.e., span{w 1 ; w 2 ; : : : ;

w n } = K n (A H ; L):
The N × n matrix

W n := [w 1 w 2 • • • w n ] ( 37
)
is called a basis matrix for K n (A H ; L). We stress that even though (34) and ( 36) are basis vectors of block Krylov subspaces, the algorithms discussed in this paper generate (34) and (36) in a vector-wise fashion, as opposed to the block-wise construction employed in more traditional block Krylov-subspace methods; see, e.g. the references given in [START_REF] Aliaga | A Lanczos-type method for multiple starting vectors[END_REF]. There are two main reasons why the vector-wise construction is preferable to a block-wise construction. First, it greatly simpliÿes both the detection of necessary de ation and the actual de ation itself. In fact, all that is required is checking if a suitable candidate vector for the next basis vector is the zero vector (for exact de ation) or close to the zero vector (for inexact de ation). Second, for Lanczos-type methods, which simultaneously construct basis vectors (34) and (36) for K n (A; R) and K n (A H ; L), respectively, only the vector-wise construction allows the treatment of the general case where the block sizes ( 28) and (33) of the de ated block Krylov sequences [START_REF] Rohrer | Passivity considerations in stability studies of numerical integration algorithms[END_REF] and (32) are not necessarily the same; for a detailed discussion, we refer the reader to [START_REF] Aliaga | A Lanczos-type method for multiple starting vectors[END_REF].

Arnoldi basis

The classical Arnoldi process [START_REF] Arnoldi | The principle of minimized iterations in the solution of the matrix eigenvalue problem[END_REF] generates orthonormal basis vectors for the sequence of Krylov subspaces K n (A; r), n¿1, induced by A and a single starting vector r. In this subsection, we state an Arnoldi-type algorithm that extends the classical algorithm to block-Krylov subspaces K n (A; R), n¿1.

Like the classical Arnoldi process, the algorithm constructs the basis vectors (34) to be orthonormal. In terms of the basis matrix (35), this orthonormality can be expressed as follows:

V H n V n = I : In addition to (34), the algorithm produces the so-called candidate vectors, Ĉn+1 ; Ĉn+2 ; : : : ; Ĉn+mc ;

(38)

for the next m c basis vectors C n+1 ; C n+2 ; : : : ; C n+mc . Here, m c = m c (n) is the number m of columns in the starting block [START_REF] Pileggi | Coping with RC(L) interconnect design headaches[END_REF], R, reduced by the number of exact and inexact de ations that have occurred so far. The candidate vectors (38) satisfy the orthogonality relation

V H n [ Ĉn+1 Ĉn+2 • • • Ĉn+mc ] = 0:
respectively. Here, m 1 is the integer given by ( 27) and [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF]. Moreover, in (40), all elements t i; l-m and t i; l that are not explicitly deÿned in Algorithm 1 are set to be zero. The compact statement of the recurrences used in Algorithm 1 is now as follows. For n¿1, we have

AV n = V n T n + [0 • • • 0 Ĉn+1 Ĉn+2 • • • Ĉn+mc ]: (41)
Furthermore, for n = m 1 , we have the relation

R = V m1 ; (42)
which re ects the fact that the ÿrst m 1 basis vectors are obtained by orthonormalizing the columns of the matrix (24), R. In ( 41) and (42), we assumed that only exact de ations are performed. If both exact and inexact de ations are performed, then an additional matrix term, say V de n , appears on the right-hand side of (41), respectively (42). The only non-zero columns of V de n are those non-zero vectors that satisÿed the de ation check (39). Since at any stage of Algorithm 1, at most m -m c = m -m c (n) vectors have been de ated, the additional matrix term is small in the sense that V de n 6dtol m -m c (n):

Lanczos basis

The classical Lanczos process [START_REF] Lanczos | An iteration method for the solution of the eigenvalue problem of linear di erential and integral operators[END_REF] generates two sequences of basis vectors (34) and (36) that span the Krylov subspaces K n (A; r) and K n (A H ; l), respectively, where r and l are single starting vectors. In [START_REF] Aliaga | A Lanczos-type method for multiple starting vectors[END_REF], a Lanczos-type method was presented that extends the classical algorithm to block-Krylov subspaces K n (A; R) and K n (A H ; L) with blocks R and L of multiple right and left starting vectors [START_REF] Pileggi | Coping with RC(L) interconnect design headaches[END_REF] and [START_REF] Wohlers | Lumped and Distributed Passive Networks[END_REF]. Such a Lanczos-type method is necessarily quite involved, and in order to keep this paper reasonably short, here we only state the governing equations that underlie the algorithm. For a complete listing of the actual algorithm, we refer the reader to [START_REF] Freund | Reduced-order modeling techniques based on Krylov subspaces and their use in circuit simulation[END_REF]Algorithm 9:2].

Like the classical Lanczos process, the extended algorithm constructs the basis vectors (34) and (36) to be biorthogonal. In terms of the associated basis matrices (35) and (37), the biorthogonality can be expressed as follows:

W H n V n = n := diag( 1 ; 2 ; : : : n ): (43)
Here, for simplicity, we have assumed that no look-ahead steps are necessary. This implies that n is a diagonal matrix, as stated in (43), and that all diagonal entries of n are nonzero. If look-ahead steps occur, then n is a block-diagonal matrix; see [START_REF] Aliaga | A Lanczos-type method for multiple starting vectors[END_REF] for further details. In addition to (34) and (36), the algorithm constructs the candidate vectors Ĉn+1 ; Ĉn+2 ; : : : ; Ĉn+mc and ŵn+1 ; ŵn+2 ; : : : ; ŵn+pc (44)

for the next m c basis vectors C n+1 ; C n+2 ; : : : ; C n+mc and the next p c basis vectors w n+1 ; w n+2 ; : : : ; w n+pc , respectively. Here, as in Section 3.4, m c = m c (n) is the number m of columns in the right starting block [START_REF] Pileggi | Coping with RC(L) interconnect design headaches[END_REF], R, reduced by the number of exact and inexact de ations of candidate vectors Ĉn that have occurred so far. Analogously, p c = p c (n) is the number p of columns in the left starting block [START_REF] Wohlers | Lumped and Distributed Passive Networks[END_REF], L, reduced by the number of exact and inexact de ations of candidate vectors ŵn that have occurred so far. Similar to (39) a de ation of the candidate vector Ĉn , respectively ŵn , is performed if, and only if, Ĉn 6dtol; respectively ŵn 6dtol;

where dtol¿0 is a suitably chosen de ation tolerance. If dtol = 0, then (45) only checks for exact de ation. The candidate vectors (44) are constructed to satisfy the following biorthogonality relations:

W H n [ Ĉn+1 Ĉn+2 • • • Ĉn+mc ] = 0; V H n [ ŵn+1 ŵn+2 • • • ŵn+pc ] = 0: (46)
The actual recurrences used to generate the basis vectors (34) and ( 36) and the candidate vectors (44) can be summarized compactly in matrix form. For simplicity, we again assume that only exact de ation is performed. Then, for n¿1, we have

AV n = V n T n + [0 • • • 0 Ĉn+1 Ĉn+2 • • • Ĉn+mc ]; A H W n = W n Tn + [0 • • • 0 ŵn+1 ŵn+2 • • • ŵn+pc ]: (47)
Moreover, for n = m 1 , respectively n = p 1 , we have the relations

V m1 = R; respectively W p1 Á = L; (48)
which summarize the recurrences for processing the starting blocks R and L. We note that the matrices T n and Tn in (47) essentially encode the same information. In fact, by pre-multiplying the ÿrst and second relation in (47) by W n and V n and by using ( 43) and ( 46), it follows that

W H n AV n = n T n = TH n n : (49)
In particular, (49) implies that TH n = n T n -1

n . Finally, we note that for symmetric transfer functions [START_REF] Guillemin | Synthesis of Passive Networks[END_REF], such as the ones describing RLC subcircuits, the Lanczos-type method sketched in this section can exploit the symmetry inherent in [START_REF] Guillemin | Synthesis of Passive Networks[END_REF]. Indeed, in this case, the Lanczos basis vectors (34) and (36) are connected as follows:

w n = n (G + s 0 C )C n for all n: (50)
Here, n = 0 are suitable normalization factors. In view of (50), only the vectors (34) need to be generated. This results in a symmetric variant of the Lanczos-type method that requires only half the computational work and storage of the general case; see [15 -17] for further details. For later use, we note that for symmetric transfer functions [START_REF] Guillemin | Synthesis of Passive Networks[END_REF], the coe cient matrices in (48) can be chosen to be identical:

= Á ∈ R m1×m : (51)
In fact, by [START_REF] Odabasioglu | PRIMA: passive reduced-order interconnect macromodeling algorithm[END_REF], [START_REF] Guillemin | Synthesis of Passive Networks[END_REF], and ( 19), (G + s 0 C )R = B = L and all these matrices are real. In view of (48) and (50), this implies (51).

Reduced-order models based on projection

In this section, we introduce two reduced-order models based on a one-sided projection onto K n (A; R), respectively a two-sided projection onto K n (A; R) and K n (A H ; L).

generated by SyMPVL is also a basis matrix for the nth block Krylov subspace associated with the nonsymmetric formulation [START_REF] Lanczos | An iteration method for the solution of the eigenvalue problem of linear di erential and integral operators[END_REF]. Hence we can also use V n to apply the one-sided projection of Section 4.1 to [START_REF] Lanczos | An iteration method for the solution of the eigenvalue problem of linear di erential and integral operators[END_REF]. By ( 21), (52), and (53), the resulting projected reduced-order transfer function is given by

H (1) n (s) := (V T n B) T (V T n (JG )V n + sV T n (JC )V n ) -1 (V T n B): (62)

Connections with Padà e-type and Padà e approximants

In this section, we show that the one-sided projection H (1) n is actually a matrix-Padà e-type approximant of H , and we review the matrix-Padà e property of H (2) n .

5.1. H (1) n is a matrix-Padà e-type approximant

Although the reduced-order transfer function (53), H (1) n , is deÿned via the simple one-sided projection (52), it satisÿes an approximation property of the form [START_REF] Freund | Reduced-order modeling of large passive linear circuits by means of the SyPVL algorithm[END_REF], where, however, q(n) is not maximal in general. This means that H (1) n is a matrix-Padà e-type approximant of H . For the special case of expansion point s 0 = 0 and a basis matrix V n generated by a simple block Arnoldi procedure without de ation, it was ÿrst observed in [START_REF] Odabasioglu | Provably passive RLC circuit reduction[END_REF][START_REF] Odabasioglu | PRIMA: passive reduced-order interconnect macromodeling algorithm[END_REF] that H (1) n satisÿes an approximation property [START_REF] Freund | Reduced-order modeling of large passive linear circuits by means of the SyPVL algorithm[END_REF]. Here, we extend this result to the most general case of arbitrary expansion points s 0 and arbitrary basis matrices V n for the properly deÿned block Krylov subspaces K n (A; R) that allow for necessary de ation of linearly dependent vectors. The only further assumption we need is that the matrix

G n + s 0 C n is nonsingular: (63)
This guarantees that s 0 is not a pole of H (1) n . Since, by (52),

G n + s 0 C n = V H n (G + s 0 C )V n ;
the condition (63) also ensures that the matrix G + s 0 C is nonsingular.

Expanding the transfer function H in (23) about s 0 , we get

H (s) = ∞ i=0 (-1) i M i (s -s 0 ) i ; where M i := L H A i R: (64)
On the other hand, expanding the reduced-order transfer function H (1) n in (52) about s 0 gives

H (1) n (s) = ∞ i=0 (-1) i M (1) i (s -s 0 ) i ; (65)
where

M (1) i := L H n ((G n + s 0 C n ) -1 C n ) i (G n + s 0 C n ) -1 B n :
We now show that for any n of the form (29), the ÿrst j terms in the expansions (64) and (65) are identical. To this end, we ÿrst establish the following proposition.

Theorem 7. Let n=m 1 +m 2 +• • •+m j ; where 16j6j max ; and let H (1) n be the reduced-order transfer function given by (52) and (53). Let s 0 ∈ C be an expansion point such that (63) is satisÿed. Then; H (1) n satisÿes H (1) n (s) = H (s) + O(s -s 0 ) j ;

i.e.; H (1) n is a matrix-Padà e-type approximant of the transfer function [START_REF] Freund | Computation of matrix Padà e approximations of transfer functions via a Lanczos-type process[END_REF]; H .

Proof. By (64) and (65), the assertion (76) is equivalent to M

i = M i for all i = 0; 1; : : : ; j -1:

By Proposition 6, the matrix F n deÿned in (66) satisÿes relations (67) and (68). Inserting (68) into (67) gives

M (1) i = L H F i n R = L H A i R = M i
for all i = 0; 1; : : : ; j -1; which is just the desired property (77).

Remark 8. By ( 28) and ( 29), we have j¿ n=m . Therefore, by Theorem 7, the Taylor expansions of H (1) n and H about s 0 match in at least the ÿrst n=m coe cient matrices. Remark 9. If p ¡ m, then a matrix-Padà e-type approximant that matches at least the ÿrst n=p Taylor coe cient matrices of H about s 0 can be obtained by performing the one-sided projection described in Section 4.1 onto K n (A H ; L), instead of K n (A; R).

H (2)

n is a matrix-Padà e approximant It turns out that, in general, the reduced-order transfer function H (2) n deÿned in (56) is even a better approximation to H than H (1) n . To properly state this result, we ÿrst deÿne the integers where the m j 's and p k 's are the integers given by ( 27), ( 28) and ( 32), (33), respectively. The main result of this section is then as follows.

Theorem 10. Let n min 6n6n max ; and let j = j(n) and k = k(n) be the maximal integers such that

m 1 + m 2 + • • • + m j 6n and p 1 + p 2 + • • • + p k 6n; (78)
respectively. Let s 0 ∈ C be an expansion point such that (63) is satisÿed; and let H (2) n be the reduced-order transfer function given by the two-sided projection (56). Then; H (2) n satisÿes H (2) n (s) = H (s) + O(s -s 0 ) q(n) ; where q

(n) = j(n) + k(n): (79)
Moreover; in (79); the exponent q(n) is as large as possible; and hence H (2) n is a matrix-Padà e approximant of the transfer function [START_REF] Freund | Computation of matrix Padà e approximations of transfer functions via a Lanczos-type process[END_REF]; H . internal terminals of the same pin and the neighboring pin no. 2, respectively. The plots show results with the projected model H (1) n and the Padà e model H (2) n , both of order n = 112, compared with an exact analysis.

In Fig. 3, we compare the relative error of the projected model H (1) 112 and the Padà e model H (2) 112 of the same size. Clearly, the Padà e model is more accurate. However, out of the 112 poles of H (2) 112 , 22 have positive real part, violating the passivity of the Padà e model. On the other hand, the projected model is passive.

An extracted RC circuit

This is an extracted RC circuit with about 4000 elements and m = 20 ports. The expansion point s 0 = 0 was used. Since the projected model and the Padà e model are identical for RC circuits, we only computed the Padà e model via SyMPVL. The point of this example is to illustrate the usefulness of the de ation procedure built into SyM-PVL. It turned out that sweeps through the ÿrst two Krylov blocks, R and AR, of the block Krylov sequence [START_REF] Raghavan | AWESpice: A general tool for the accurate and e cient simulation of interconnect problems[END_REF] were su cient to obtain a reduced-order model that matches the transfer function in the frequency range of interest. During the sweep through the second block, 6 almost linearly dependent vectors were discovered and de ated. As a result, the reduced-order model obtained with de ation is only of size n = 2m -6 = 34. When SyMPVL was rerun on this example, with de ation turned o , a reduced-order model of size n = 40 was needed to match the transfer function. In Fig. 4, we show the H 1; 11 component of the reduced-order model obtained with de ation and without de ation, compared to the exact transfer function. Clearly, de ation leads to a signiÿcantly smaller reduced-order model that is as accurate as the bigger one generated without de ation.

Concluding remarks

In the last few years, reduced-order modeling techniques based on Krylov subspaces have become indispensable tools for tackling the large linear subcircuits that arise in the simulation of electronic circuits. Much of this development was and continues to be driven by the emerging need to accurately simulate the interconnect of electronic circuits. Today, circuit interconnect is typically modeled as large linear passive subcircuits that are generated by automatic parasitics-extraction programs. Using reduced-order modeling techniques has become crucial in order to reduce these subcircuits to a size that is manageable for circuit simulators.

To guarantee stability of the overall simulation, it is crucial that passive subcircuits are approximated by passive reduced-order models. While reduced-order models based on projection are passive, they are -in terms of number of matched Taylor coe cients -only half as accurate as the corresponding, in general non-passive, Padà e models of the same size. It remains an open problem to describe and construct reduced-order models of a given size that are both passive and of maximal possible accuracy.

Finally, today's circuit simulation is based on the paradigm of lumped circuit elements, which leads to systems of DAEs. As circuit feature sizes continue to decrease and circuit speeds continue to increase, feature sizes are becoming comparable in size with signal wavelengths. As a result, at least parts of a circuit must be modeled as distributed elements, such as transmission lines. Including distributed elements in the simulation paradigm requires a fusion of traditional lumped circuit simulation and electromagnetic simulation. Electromagnetic simulation, however, involves systems of partial di erential equations (PDEs). Combining lumped circuit simulation with electromagnetic simulation will thus require e cient techniques for the solution of very large systems of DAEs coupled with PDEs. One of the challenges then is to develop reduced-order modeling techniques that allow to replace parts of such coupled systems with much smaller models. Research into and development of such techniques have hardly begun.

n

  min := max{m 1 ; p 1 } and n max := min

Fig. 1 .

 1 Fig. 1. Package: Pin no. 1 external to Pin no. 1 internal, exact, projected model, and Padà e model.

Fig. 2 .

 2 Fig. 2. Package: Pin no. 1 external to Pin no. 2 internal, exact, projected model, and Padà e model.

Fig. 3 .

 3 Fig. 3. Relative error of projected model and Padà e model.

Fig. 4 .

 4 Fig. 4. Impedance H1;11.

Due to the vector-wise construction of (34) and (38), detection of necessary de ation and the actual de ation becomes trivial. In fact, essentially the same proof as given in [START_REF] Aliaga | A Lanczos-type method for multiple starting vectors[END_REF] for the case of a Lanczos-type algorithm can be used to show that exact de ation at step n of the Arnoldi-type process occurs if, and only if, Ĉn = 0. Similarly, inexact de ation occurs if, and only if, Ĉn ≈ 0, but Ĉn = 0. Therefore, in the algorithm, one checks if Ĉn 6dtol;

(39

where dtol¿0 is a suitably chosen de ation tolerance. If (39) is satisÿed, then Ĉn is de ated by deleting Ĉn , shifting the indices of all remaining candidate vectors by -1, and setting m c = m c -1.

If this results in m c = 0, then the block-Krylov subspace is exhausted and the algorithm is stopped.

Otherwise, the de ation procedure is repeated until a vector Ĉn with Ĉn ¿ dtol is found. This vector is then turned into C n by normalizing it to have Euclidean norm 1.

A complete statement of the resulting Arnoldi-type algorithm is as follows.

Algorithm 1 (Construction of Arnoldi basis for K n (A; R).).

(0) Set Ĉi = r i for i = 1; 2; : : : ; m. Set m c = m. For n = 1; 2; : : : ; do:

(1) Compute Ĉn and check if the de ation criterion (39) is fulÿlled.

If yes; Ĉn is de ated by doing the following: Set m c = m c -1. If m c = 0; set n = n -1 and stop. Set Ĉi = Ĉi+1 for i = n; n + 1; : : : ; n + m c -1.

Return to

Step (1). (2) Set t n; n-mc = Ĉn and C n = Ĉn =t n; n-mc .

(3) Compute Ĉn+mc = AC n . (4) For i = 1; 2; : : : ; n do:

Remark 2. If dtol = 0 in (39), then Algorithm 1 performs only exact de ation.

Remark 3. Other block-Arnoldi algorithms (all without de ation though) can be found in [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF]Section 6:12].

After n passes through the main loop, Algorithm 1 has constructed the ÿrst n basis vectors (34) and the candidate vectors (38) for the next m c basis vectors. In terms of the basis matrix (35), V n , the recurrences used to generate all these vectors can be written compactly in matrix format. To this end, we collect the recurrence coe cients computed during the ÿrst n = m 1 and n¿1 passes through the main loop of Algorithm 1 in the matrices := [t i; l-m ] i=1; 2; :::; m1 l=1; 2; :::; m and T n := [t i; l ] i=1; 2; :::; n l=1; 2; :::; n ;

(40)

Let V n ∈ C N ×n be a given matrix with full column rank n. By simply projecting the matrices G , C , B, and L in the original linear dynamical system (8) onto the subspace spanned by the columns of V n , we obtain a reduced-order model [START_REF] Feldmann | Reduced-order modeling of large linear subcircuits via a block Lanczos algorithm[END_REF] with matrices G n , C n , B n , and L n given by

In particular, we now assume that V n is a basis matrix (35) for K n (A; R). In this case, the reduced-order model deÿned by ( 9) and (35) represents a (one-sided) projection of the original system (8) onto the nth block-Krylov subspace K n (A; R). In the sequel, we denote the associated reduced-order transfer function by

where the superscript (1) indicates the one-sided projection. In Section 5.1 below, we show that H (1) n is a certain Padà e-type approximant of the original transfer function H .

The following proposition shows that H (1) n is independent of the actual choice of the basis matrix V n for K n (A; R). Proposition 4. The reduced-order transfer function H (1) n given by (52) and (53) does not depend on the particular choice of the basis matrix (35); V n ; for the nth block Krylov subspace K n (A; R).

Proof. Let V n be the basis matrix for K n (A; R) that is used in (52). Let Ṽn be any other basis matrix for K n (A; R). In analogy to (52) and (53), Ṽn induces the reduced-order transfer function

where

and Ṽn are basis matrices for the same subspace, there exists a nonsingular n × n matrix M such that Ṽn = V n M . Using this relation, we obtain from (54) and (52) that

Inserting (55) into (54), we get

n (s): Thus, the reduced-order transfer functions H (1) n and H n are identical.

Two-sided projection onto K

Let V n and W n be any two basis matrices of K n (A; R) and K n (A H ; L), and consider the representation [START_REF] Odabasioglu | PRIMA: passive reduced-order interconnect macromodeling algorithm[END_REF] of the transfer function H of [START_REF] Feldmann | E cient linear circuit analysis by Padà e approximation via the Lanczos process[END_REF]. By projecting the matrices in [START_REF] Odabasioglu | PRIMA: passive reduced-order interconnect macromodeling algorithm[END_REF] from the right and left onto the columns of V n and W n , respectively, we obtain the reduced-order transfer function

In analogy to Proposition 4, we have the following result.

Proposition 5. The reduced-order transfer function H (2) n given by (56) does not depend on the particular choice of the basis matrices V n and W n for the nth block Krylov subspaces K n (A; R) and K n (A H ; L).

Proof. Analogous to the proof of Proposition 4.

Computation via Krylov-subspace algorithms

In view of Proposition 5, the reduced-order transfer function (56), H (2) n , can be generated from any two basis matrices V n and W n . However, there is one distinguished choice of V n and W n that eliminates the need to explicitly compute the projections in (56). This choice is the Lanczos basis described in Section 3.5.

Indeed, let n¿max{m 1 ; p 1 }, and assume that the Lanczos-type algorithm is run for n steps. Let n be the diagonal matrix deÿned in (43), and let T n , , and Á be the matrices of recurrence coe cients given by ( 47) and (48). Then, from (43) and (48), it readily follows that

where

Furthermore, by multiplying the ÿrst relation in (47) from the left by W H n and using the ÿrst relation in (46), as well as (43), we get

Inserting ( 57) and ( 59) into (56), it readily follows that

The MPVL (matrix-Padà e via Lanczos) algorithm, which was ÿrst proposed in [START_REF] Feldmann | Reduced-order modeling of large linear subcircuits via a block Lanczos algorithm[END_REF], is a numerical procedure for computing H (2) n via the formula (60). For symmetric transfer functions [START_REF] Guillemin | Synthesis of Passive Networks[END_REF], by (51) and (58), the reduced-order transfer function (60) is also symmetric:

where -1

are real symmetric matrices. The SyMPVL algorithm [START_REF] Freund | The SyMPVL algorithm and its applications to interconnect simulation[END_REF][START_REF] Freund | Reduced-order modeling of large linear passive multi-terminal circuits using matrix-Padà e approximation[END_REF] is a special symmetric variant of the general MPVL algorithm that computes symmetric reduced-order transfer functions (61).

Furthermore, recall from Section 2.4 that RLC subcircuits are described by special symmetric transfer functions [START_REF] Guillemin | Synthesis of Passive Networks[END_REF] with matrices G , C , and B of the form ( 16) and [START_REF] Freund | Reduced-order modeling of large linear passive multi-terminal circuits using matrix-Padà e approximation[END_REF]. In this case, as we will discuss in Section 6, the reduced-order transfer function (60) in general does not preserve the passivity of the RLC subcircuit. However, one can easily augment the SyMPVL algorithm to generate a second projected reduced-order model that, by Corollary 14 below, is always passive. To this end, let J be the matrix deÿned in [START_REF] Kim | Time-domain macromodels for VLSI interconnect analysis[END_REF], and consider the nonsymmetric formulation [START_REF] Lanczos | An iteration method for the solution of the eigenvalue problem of linear di erential and integral operators[END_REF] of the symmetric transfer function [START_REF] Guillemin | Synthesis of Passive Networks[END_REF]. Note that by Lemma 1, both formulations [START_REF] Guillemin | Synthesis of Passive Networks[END_REF] and [START_REF] Lanczos | An iteration method for the solution of the eigenvalue problem of linear di erential and integral operators[END_REF] result in the same nth block Krylov subspace K n (A; R). In particular, the Lanczos basis matrix V n Proposition 6. Let n = m 1 + m 2 + • • • + m j ; where 16j6j max . Then; the matrix

satisÿes the following relations:

for all i = 0; 1; : : : ;

(67)

R for all i = 0; 1; : : : ; j -1: (68)

Proof. By ( 29) and ( 30), for each i = 1; 2; : : : ; j, the columns of the matrix A i R are all contained in K n (A; R). Since V n is a basis matrix for K n (A; R), for each i = 1; 2; : : : ; j, there exists an n × m matrix E i such that

We now prove (67). From ( 23) and (69) (for i = 1), we get

Multiplying (70) from the left by V H n and using (52), it follows that

Moreover, from (66), we obtain the relation

C n , which, by induction on i, implies that

for all i = 0; 1; : : : :

Note that, by (52), L H n = L H V n . Using this relation, as well as (65), (71), (72), and (69) (for i = 1), it follows that, for all i = 0; 1; : : : ;

Next, we prove (68) using induction on i. For i = 0, (68) is trivially satisÿed. Now assume that (68) is true for some 06i ¡ j -1. We show that (68) then also holds true for i + 1, i.e.,

Using ( 23), (68), and (69) (with i replaced by i + 1), we get

Multiplying (74) from the left by

Using (66) and (69) (with i replaced by i + 1), we obtain from (75) that

which is just the desired relation (73).

Proof. In [START_REF] Freund | Computation of matrix Padà e approximations of transfer functions via a Lanczos-type process[END_REF], we studied the reduced-order transfer function H (2) n given by (60), where n , Á n , n , and T n are the matrices generated by n steps of the Lanczos-type method sketched in Section 3.5. In particular, in [12, Theorem 1], we showed that H (2) n satisÿes the properties listed in Theorem 10 above. Recall from Section 4.2 that the reduced-order transfer functions deÿned in (56) via a two-sided projection and the one given by (60) in terms of the Lanczos-type method are identical. Therefore, the assertions in Theorem 10 follow from [12, Theorem 1].

Remark 11. In view of ( 28), (33), and (78), we have j(n)¿ n=m and k(n)¿ n=p . Therefore, by Theorem 10, the Taylor expansions of H (2) n and H about s 0 match in at least the ÿrst n=m + n=p coe cient matrices.

Passivity and stability

As we discussed in Section 2, in circuit simulation, reduced-order modeling is mostly applied to large passive linear subcircuits, such as RLC networks consisting of only resistors, inductors, and capacitors. When reduced-order models of such subcircuits are used within a simulation of the whole circuit, stability of the overall simulation can only be guaranteed if the reduced-order models preserve the passivity of the original subcircuits; see, e.g. [START_REF] Chirlian | Integrated and Active Network Analysis and Synthesis[END_REF][START_REF] Rohrer | Passivity considerations in stability studies of numerical integration algorithms[END_REF]. Unfortunately, except for special cases such as RC subcircuits consisting of only resistors and capacitors, the Padà e model given by H (2) n is not passive in general; see, e.g. [START_REF] Bai | How to make theoretically passive reduced-order models passive in practice[END_REF][START_REF] Elfadel | Zeros and passivity of Arnoldi-reduced-order models for interconnect networks[END_REF][START_REF] Freund | Passive reduced-order models for interconnect simulation and their computation via Krylov-subspace algorithms[END_REF][START_REF] Freund | Reduced-order modeling of large passive linear circuits by means of the SyPVL algorithm[END_REF][START_REF] Kerns | Preservation of passivity during RLC network reduction via split congruence transformations[END_REF]. In this section, we derive a su cient criterion for passivity of general transfer functions, and then apply the criterion to establish passivity of the particular projected model given by (62).

Roughly speaking, a (linear or nonlinear) dynamical system is passive if it does not generate energy. The concept was ÿrst used in circuit theory; see, e.g. [START_REF] Anderson | Network Analysis and Synthesis[END_REF][START_REF] Guillemin | Synthesis of Passive Networks[END_REF]. For example, networks consisting of only resistors, inductors, and capacitors are passive.

The following well-known theorem (see, e.g. [START_REF] Anderson | Network Analysis and Synthesis[END_REF][START_REF] Wohlers | Lumped and Distributed Passive Networks[END_REF]) relates the passivity of the linear dynamical system [START_REF] Feldmann | E cient linear circuit analysis by Padà e approximation via the Lanczos process[END_REF] to the positive realness of its transfer function. Here and in the sequel, we assume that m = p in (8).

Theorem A. The linear dynamical system (8) is passive if; and only if; the associated transfer function [START_REF] Freund | Computation of matrix Padà e approximations of transfer functions via a Lanczos-type process[END_REF]; H , is positive real.

The deÿnition of a positive real matrix-valued function is as follows; see, e.g. [START_REF] Anderson | Network Analysis and Synthesis[END_REF].

Recall that a function H : C → (C ∪ {∞}) m×m is stable if H has no poles in C + and if all possible purely imaginary poles of H are simple. It is well known that any positive real function is necessary stable.

Next, we prove the following su cient condition for positive realness.

Theorem 13. Let G ; C ∈ R N ×N ; and B ∈ R N ×m . Assume that G + G T ¿0; C = C T ¿0; and that G + sC is a regular matrix pencil. Then; the transfer function

is positive real.

Proof. We need to show that H satisÿes the conditions (i) -(iii) given in Deÿnition 12. Condition (ii) follows directly from the fact that the matrices G , C , and B in (80) are real. Next, we verify condition (iii). Let s=s 1 +is 2 ∈ C + be arbitrary, but ÿxed. Here i := √ -1 denotes the purely imaginary unit. Note that

where

, and s 1 = Re s ¿ 0. This guarantees that S¿0 and K = -K H . These properties imply that y H Sy¿0 and Re(y H Ky) = 0 for all y ∈ C N . Therefore, by (81), we have Re(y H (G + sC ) H y) = y H Sy¿0 for all y ∈ C N :

(82)

We now assume that s ∈ C + is such that the matrix G + sC is nonsingular. Furthermore, let x ∈ C m be arbitrary, and set

Then, by (80) and (83), we have

Combining ( 82) and (84), it follows that Re(x H H (s)x) = Re(y H (G + sC ) H y)¿0 for all x ∈ C m (85) and for all s ∈ C + for which G + sC is nonsingular. Now let ŝ ∈ C + be such that G + ŝC is singular. Note that there are at most N such "singular" points ŝ, since G + sC is assumed to be a regular matrix pencil. Therefore, each ŝ is an isolated point in C + , i.e., there exists an = ( ŝ) ¿ 0 and a (punctured) neighborhood

of ŝ such that D ⊂ C + and the matrix G + sC is nonsingular for all s ∈ D . Thus (85) holds true for all s ∈ D . If ŝ is not a pole of the rational function H , then H ( ŝ) = lim s→ ŝ H (s) is a ÿnite m × m matrix. In this case, by taking limits s → ŝ in (85), it follows that (85) also holds true for s = ŝ. Now suppose that ŝ is a pole of H . Then at least one of the components h jk (s) of H (s) = [h jk (s)] 16j; k6m has a pole at ŝ. Such an h jk (s) maps D onto a suitable neighborhood of ∞ in the complex plane and, in particular, attains large negative numbers in D . By selecting a suitable component h jk of H and an associated vector x ∈ C N , it is thus possible to ÿnd a point s ∈ D such that Re(x H H (s)x) ¡ 0. However, this is a contradiction to (85), and therefore, ŝ cannot be a pole of H . This concludes the proof of (iii). It remains to verify condition (i). By (80), if ŝ is a pole of H , then the matrix G + ŝC is necessarily singular, i.e., ŝ is a singular point. However, we have just shown that there are no such singular points ŝ ∈ C + . Consequently, H cannot have poles in C + .

The matrix function H satisÿes all three conditions (i) -(iii), and hence H is positive real.

Finally, we apply Theorem 13 to the reduced-order transfer function (62).

Corollary 14. Let H be the transfer function given by [START_REF] Lanczos | An iteration method for the solution of the eigenvalue problem of linear di erential and integral operators[END_REF] with matrices that satisfy [START_REF] Odabasioglu | Provably passive RLC circuit reduction[END_REF]. Let V n ∈ R N ×n have rank n and assume that the matrix pencil

is regular. Then; the reduced-order transfer function

is positive real; and thus the reduced-order model given by (87) is passive.

Proof. By ( 22) and ( 86), it follows that G n + G T n ¿0 and C n = C T n ¿0. The transfer function (87), H (1) n , is thus positive real by Theorem 13.

Numerical examples

In this section, we present two circuit examples.

A package model

The ÿrst example arises in the analysis of a 64-pin package model used for an RF integrated circuit. Only eight of the package pins carry signals, the rest being either unused or carrying supply voltages. The package is characterized as a passive linear dynamical system with m = p = 16 inputs and outputs, representing eight exterior and eight interior terminals. The package model is described by approximately 4000 circuit elements, resistors, capacitors, inductors, and inductive couplings, resulting in a linear dynamical system with a state-space dimension of about 2000.

In [START_REF] Freund | The SyMPVL algorithm and its applications to interconnect simulation[END_REF], SyMPVL was used to compute a Padà e-based reduced-order model (61) of the package, and it was found that a model H (2) n of order n = 80 is su cient to match the transfer-function components of interest. However, the model H (2) n has a few poles in the right half of the complex plane, and therefore, it is not passive.

In order to obtain a passive reduced-order model, we ran SyMPVL again on the package example, and this time, also generated the projected reduced-order model H (1) n given by (62). The expansion point s 0 =5 × 10 9 was used. Recall that H (1) n is only a Padà e-type approximant and thus less accurate than the Padà e approximant H (2) n . Therefore, one now has to go to order n = 112 to obtain a projected reduced-order model H (1) n that matches the transfer-function components of interest.