Patent Publication Number: US-6041170-A

Title: Apparatus and method for analyzing passive circuits using reduced-order modeling of large linear subcircuits

Description:
FIELD OF THE INVENTION 
     The invention pertains to a method and apparatus for analyzing electrical circuits. More particularly, the invention pertains to a method and apparatus for analyzing linear subcircuits of a larger circuit by means of a Lanczos type procedure to reduce original linear circuit matrices to a pair of banded smaller matrices. 
     BACKGROUND OF THE INVENTION 
     Electrical circuits typically are comprised of linear circuit elements, such as resistors, inductors, and capacitors as well as non-linear circuit elements, such as transistors and power sources. In the process of designing circuits, particularly large scale integrated circuits, it is normal practice to mathematically model the circuit. Particularly, the outputs of the circuit are modeled as a function of the inputs to the circuit. The mathematical model is used to determine various response characteristics of the circuit. 
     In our prior patent applications, application Ser. No. 08/269,230, filed Jun. 30, 1994, now U.S. Pat. No. 5,537,329 and application Ser. No. 08/489,270, filed Jun. 9, 1995, now U.S. Pat. No. 5,689,685, both of which are incorporated herein by reference, we disclosed new methods and apparatus for modeling linear subcircuits within a larger electrical circuit. 
     An electrical circuit, such as an integrated circuit, frequently contains large portions thereof comprised entirely of linear circuit elements. For instance, a wire or other connector usually can be modeled as a network of resistances and capacitances. Accordingly, the entire interconnect structure between a last non-linear circuit element in an integrated circuit and the end of a package pin can be mathematically modeled as a linear subcircuit. 
     The modeling of the larger overall circuit can be extremely complex in large scale integrated (LSI) circuits and very large scale integrated (VLSI) circuits. Our aforementioned patent application Ser. No. 08/489,270 discloses a technique for modeling the linear subcircuit portions within the larger circuit by a greatly reduced matrix compared to the original matrix equation for the linear subcircuit. 
     More specifically, each connection point of a linear subcircuit to elements outside of the bounds of the linear subcircuit is modeled as both an input and an output in a matrix formula for the linear subcircuit. Accordingly, a linear subcircuit having p electrical connection points to circuitry outside of the bounds of the linear subcircuit is modeled using a p by p matrix function. 
     Referring to our aforementioned patent applications, in U.S. patent application Ser. No. 08/269,230, we disclosed a Pade via Lanczos method for determining a scalar Pade approximant of a transfer function of a linear circuit for use in determining, with substantial accuracy, the frequency response behavior of the linear circuit. 
     In U.S. patent application Ser. No. 08/489,270, we disclosed an extension to the previous application in which the Pade via Lanczos method was extended to matrix transfer functions. More specifically, U.S. patent application Ser. No. 08/489,270 discloses a numerically stable method to determine a matrix Pade approximation of the frequency response of a linear circuit or subcircuit with little numerical degradation, and providing an acceptable computational cost per order of approximation. The method reduces the very large matrices used to represent the p by p matrix transfer function to a pair of much smaller matrices so that the resulting approximate matrix transfer function has basically the same characteristics as the original matrix transfer function. Using the method and apparatus disclosed therein, the original matrix can be reduced to any particular size desired. The accuracy of the approximation depends on the order of reduction in size of the matrices. Generally, the matrices can be reduced to a very large extent without significant decrease in accuracy. 
     One potential application for the reduction methods disclosed in the aforementioned patent applications is use on one or more linear subcircuits in a larger circuit. In such an application, the reduced equations can be inserted into the equation for the larger overall circuit, thus greatly reducing the complexity of the mathematical model for the overall non-linear circuit. While such application of these techniques can greatly simplify analysis of the larger non-linear circuit, the reductions still may remain very complex and be very time consuming in themselves. Accordingly, it is desireable to further reduce the complexity of the procedure for generating reduced-order circuit models to whatever extent possible in order to minimize the required processing power and/or time. Also, there is an ever-present drive to increase the accuracy of the models. 
     Therefore, it is an object of the present invention to provide an improved method and apparatus for generating a reduced order matrix approximant of the transfer function associated with a passive circuit. 
     SUMMARY OF THE INVENTION 
     The present invention is a method and apparatus for analyzing linear circuits using reduced-order mathematical matrix approximations thereof generated by a matrix Pade via Lanczos type procedure. In particular, it is a specialized adaptation of the previously known matrix Pade via Lanczos method. The present invention is particularly adapted for use with passive, linear circuits. In particular, the restriction to passive, linear circuits guarantees that the circuit matrices will be symmetric. When the additional condition is met that the circuit (or subcircuit) contains only two types of linear circuit elements, (i.e., resistor and inductor, resistor and capacitor, or inductor and capacitor), the matrices are guaranteed to be not only symmetric, but also positive definite, which guarantees stability and passivity of the approximation. 
     The present invention is a method and apparatus for the computation of the symmetric multi-port transfer function of an RLC circuit. It employs a novel symmetric block-Lanczos procedure to reduce original circuit matrices to a pair of banded symmetric reduced-order matrices, herein termed SyMPVL (for Symmetric Matrix Pade via Lanczos type procedure). When the circuit comprises only two of the three RLC components, the matrices are also positive definite, guaranteeing stability and passivity. The reduced transfer function represents a matrix-Pade approximation of the original multi-port matrix transfer function. 
    
    
     BRIEF DESCRIPTION OF THE DRAWINGS 
     FIG. 1 is a block diagram of the components of an exemplary apparatus in accordance with the present invention. 
     FIG. 2 is a flow chart illustrating an exemplary set of steps in accordance with the present invention employing a symmetric matrix Pade approximant of the frequency response of a passive, linear circuit. 
     FIG. 3 is a flow chart illustrating the details of step 34 of FIG. 2. 
     FIG. 4 is a partial circuit diagram of an exemplary circuit that can be analyzed using the method and apparatus of the present invention. 
     FIG. 5 is a diagram of the magnitude of the transfer function of the circuit illustrated in FIG. 4. 
     FIG. 6 is a diagram of an exemplary integrated circuit package that can be analyzed using the method and apparatus of the present invention. 
     FIG. 7 is a diagram of the magnitude of the transfer function with respect to a first pin of the exemplary integrated circuit package of FIG. 6. 
     FIG. 8 is a diagram of the magnitude of the transfer function with respect to a second pin of an exemplary integrated circuit package. 
    
    
     DETAILED DESCRIPTION OF THE INVENTION 
     Referring to the drawings, in which like reference numerals identify similar or identical elements, as shown in FIG. 1, the present disclosure describes an apparatus and method implementing a circuit analyzer 10, which includes a processing unit 12 having a processor 14, memory 16, and stored programs 18, including a matrix transfer function program 20; an input device 22; and an output device 24. 
     The present invention is an improvement over the method and apparatus disclosed in aforementioned application Ser. No. 08/489,270, which can cut the processing load by approximately one-half while still generating a highly accurate model. The new method and apparatus disclosed herein is particularly adapted for use in connection with purely passive circuits or subcircuits, i.e., circuits which can be modeled with symmetric matrices. Further, the method and apparatus guarantees stability and passivity of the model if the circuit is an RL, RC, or LC circuit. 
     In an exemplary embodiment, the processing unit 12 is a SPARC workstation available from SUN MICROSYSTEMS, INC. having about 32 MB associated RAM memory and a hard or fixed drive memory 15. The processor 14 operates using the UNIX operating system to run application software 18 implementing the disclosed circuit analyzer 10 and method. 
     It should be understood that the above described processor is merely exemplary and that any computer or processing system having sufficient processing ability and memory may be employed. 
     A description of the circuit to be analyzed is provided from a circuit characteristic data source 26 through an input device 22 to the processor 14. The data source 26 may be any source for providing the circuit characteristic information to the computer, such as a floppy disk containing a circuit characteristic file or a keyboard on which is typed the circuit characteristic data. The input device 22 may be a keyboard, a floppy disk drive or other data reading device. 
     The processor 14 receives commands and the circuit characteristic data from the circuit characteristic data source 26 through the input device 22. The received circuit characteristic data may be stored in memory 16 for further processing to generate a frequency response signal representing the approximated or predicted frequency response of the circuit. The generated frequency response signal of a linear subcircuit may be a sub-goal in the overall analysis of a larger, non-linear circuit. 
     The generated frequency response signal may be used by the processor 14 to generate a graphic representation for plotting the approximated frequency response in volts, amps, or other units as a function of frequency. The graphic representation is sent to an output device 24 such as a display for displaying the predicted frequency response of the circuit under analysis. Alternately, the output device 24 may include specialized graphics programs to convert the generated frequency response signal to a displayed graphic. In additional embodiments, the generated frequency response signal may include determined poles, zeros, and/or residues listed in a file for output as columns or tables of text by the output device 24 which may be a display or a hard copy printer. 
     The circuit analyzer 10 performs the application programs and subroutines described hereinbelow in conjunction with FIGS. 2-3, which are implemented from compiled source code in the C++ and/or FORTRAN programming languages. 
     Referring to FIG. 2, the circuit analyzer 10 includes software for generating a frequency response signal of a circuit. In step 28, the program for generating the frequency response of the circuits is started. In step 30, the circuit characteristic data is received and stored in memory 16 for use with the program. In step 32, the program generates two circuit characteristic matrices G and C that mathematically describe the circuit. These matrices are square matrices of the same size, the size typically being of the order of the number of circuit elements in the circuit or subcircuit. 
     In step 34, the matrices are reduced in size, remaining symmetric, by use of a look-ahead block Lanczos procedure. The results of this reduction are the matrices T n  and Δ n . The look-ahead block Lanczos procedure carried out in step 34 also produces the matrix, ρ, which contains input/output information. 
     In step 36, the matrices of ρ, T n  and Δ n  are used to generate a matrix Pade approximant, Z n  (s), of a transfer function of the circuit, where s is complex frequency. In step 38, the matrix Pade approximant is applied to generate a graphic representation of the circuit response predicted by the approximant. The representation typically may be a graph of the signal (e.g., voltage or current magnitude) at a particular terminal of the circuit plotted against frequency, such as illustrated in FIGS. 5, 7 and 8. 
     FIG. 3 shows the details of step 34 for reducing the characteristic matrices G and C. The subroutine is started in step 40. In step 42, the Lanczos vectors are initialized to initial iteration values. In step 43, a value of zero is assigned to variable n. Since the steps to follow will be iterated for values of n from 1 to n max , n is incremented by one in step 44. In step 46, a normalization factor t n , n-p is determined as an element of the banded matrix T n . 
     In step 48, it is determined if the normalization factor is substantially near zero. If so, it is then determined, in step 49, whether p is equal to 1, where p started out as the number of connection points to outside circuitry in the circuit under analysis, but which is being decremented every time the normalization factor is determined to be substantially near zero, as explained more fully below. If p has reached one, the iterations are stopped in step 50. Otherwise, p is decremented by 1 in step 52 and control flows back to step 46 for further iterations through steps 46 to 52 until p is equal to 1. If the normalization factor is not substantially near zero, control flows to steps 54 through 60. 
     In step 54, the Lanczos vector for the present value of n as well as corresponding entries in the matrix ρ n  are determined. In step 56, the values of auxiliary vectors are determined, and, in step 58, the value of the auxiliary vectors are updated. In step 60, it is determined if the last iteration, i.e., for n=n max , has been performed. If not, the process returns to loop through steps 44 through 60 until either n has reached n max  or p has reached one. When the iterations are completed by virtue of either n having reached n max  or p having reached one, flow proceeds to step 62 where the banded matrix T n  is formed from the iterated values of the variables and vectors and is returned to step 36 in FIG. 2 along with the matrices ρ n  and ρ n . 
     Returning to FIG. 2, in step 36, the generated banded matrices T n  and Δ n , together with the matrix ρ generated in step 34 are used to determine a matrix Pade approximant of the frequency response of the circuit in a manner comparable to the method described in aforementioned U.S. application Ser. No. 08/489,270. For example, eigenvalues and eigenvectors of T n  may be generated to determine the poles and residues of the frequency response of the circuit, and T n  may be used to calculate the zeros of the frequency response of the circuit. 
     As mentioned above, U.S. application Ser. No. 08/489,270, U.S. Pat. No. 5,689,685, discloses a method and apparatus termed MPVL (Matrix Pade via Lanczos-type process), for the accurate and efficient computation of reduced-order models of large linear circuits. MPVL generates a reduced-order model, the transfer function of which represents the matrix-Pade approximation of the original circuit&#39;s matrix transfer function. MPVL is a general procedure, in the sense that it can stably compute reduced-order models for any linear system and for different number of inputs and outputs up to the desired level of accuracy. In contrast, RLC circuits can be characterized in terms of square and symmetric matrix transfer functions. Moreover, all the matrices resulting from formulating the RLC circuit equations are symmetric. Direct application of the MPVL procedure to RLC circuits does not take advantage of this special structure and cannot guarantee in general the preservation of stability and passivity of the approximation. 
     The preservation of the stability and passivity properties of the circuits being modeled is a significant requirement for reduced-order models. RLC circuits, being composed entirely of passive components, are always stable and passive. This property should not be violated by the reduced-order model. 
     The invention is a method and apparatus that computes symmetric matrix-Pade approximations of the circuit matrix transfer function. It employs a novel symmetric block-Lanczos procedure to reduce original circuit matrices to a pair of banded symmetric reduced-order matrices, herein termed SyMPVL (for Symmetric Matrix Pade via Lanczos type procedure). 
     I. Symmetric Formulation of Circuit Equations 
     This section formulates the equations that describe general RLC circuits, and discusses the special cases of RC, RL, and LC circuits. The matrices associated with the circuit equations are always symmetric and, in the special cases, are also positive semi-definite. These properties are important for SyMPVL and for stability and passivity proofs. 
     A. General RLC Circuits 
     This subsection discloses the details of step 32 of FIG. 2. The connectivity of a circuit can be captured using a directional graph. The nodes of the graph correspond to the nodes of the circuit, and the edges of the graph correspond to each of the circuit elements. An arbitrary direction is assigned to graph edges, so one can distinguish between the source and destination nodes. The adjacency matrix, A, of the directional graph describes the connectivity of a circuit. Each row of the matrix corresponds to a graph edge and, therefore, to a circuit element. Each column of the matrix corresponds to a graph or circuit node. The column corresponding to the datum (ground) node of the circuit is omitted in order to remove redundancy. By convention, a row of the connectivity matrix will contain +1 in the column corresponding to the source node, -1 in the column corresponding to the destination node, and 0 everywhere else. It is easy to see that Kirchhoff&#39;s laws, which depend only on connectivity can be expressed using the adjacency matrix 
     
         KCL: A.sup.T i.sub.b =0 
    
     
         KCL: Av.sub.n =v.sub.b                                     (Eqs. 1) 
    
     where i b  and v b  are the vectors of branch currents and voltages, respectively, and v n  is the vector of the non-datum node voltages. 
     We are interested in analyzing RLC circuits and, for simplicity, we assume that the circuit is excited only by current sources. In this case, the adjacency matrix and the branch current and voltage vectors can be partitioned according to circuit-element types as follows: ##EQU1## Here, the subscripts i, g, c, and l stand for branches containing current sources, resistors, capacitors, and inductors, respectively. 
     The set of circuit equations is completed by adding the so-called branch constitutive relationships (BCR&#39;s), which describe the physical behavior of the circuit elements. In the case of RLC circuits, the BCR&#39;s are as follows: ##EQU2## Here, I l  (t) is the vector of current-source values, ζ and  are appropriately-sized diagonal matrices whose diagonal entries are the conductance and capacitance values of each element. It is clear that these values are positive for any physical circuit. The matrix L also is diagonal in the absence of inductive coupling. Inductive coupling introduces off-diagonal terms in the inductance matrix, but L remains symmetric and positive definite. 
     The modified nodal formulation (MNA) of the circuit equations is obtained by combining the Kirchhoff equations (Eqs. 1) with the BCRs (Eq. 3), and eliminating as many current unknowns as possible. For the case of RLC circuits, only inductor currents need to be left as unknowns. The resulting MNA equations are as follows: ##EQU3## Writing (Eq. 4) in matrix form, we obtain ##EQU4## Note that the matrices G and C are symmetric. They are however indefinite in general. 
     We are interested in determining the network functions of the RLC block viewed as a p-terminal component. Since we allowed only the current sources in our formulation, it is natural to determine the Z-parameters. By applying the Laplace transform and assuming zero initial conditions, we obtain 
     
         (G+sC)X=BI.sub.s (s), V.sub.i =B.sup.T X.                  (Eq. 6) 
    
     Here, X, I s  (s), and V i  represent the Laplace transforms of the unknown vector x, the excitation current I l  (t), and the vector of voltages across the excitation sources, respectively. Eliminating X in (Eq. 6) gives ##EQU5## While for general RLC circuits the matrices G and C are indefinite, in some important special cases the formulation can be modified such that the matrices become positive semi-definite. In these special cases, important passivity and stability results can be proved. 
     B. Special Cases of RC, RL and LC Circuits 
     In the case of RC circuits, the matrices L, A l , and the vector i l  are empty. The Laplace transform of the MNA equations (Eq. 4) thus reduces to ##EQU6## Clearly, the matrices C and G are symmetric positive semi-definite. 
     In the case of RL circuits, the matrices  and A c  are empty. The Laplace transform of the MNA equations (Eq. 4) now becomes 
     
         A.sub.g.sup.T ζA.sub.g V.sub.n +A.sub.l.sup.T I.sub.l =A.sub.i.sup.T I.sub.s (s), A.sub.l V.sub.n -LsI.sub.l =0.               (Eq. 9) 
    
     By multiplying the first equation in (Eq. 9) by s and substituting sI l  from the second, we obtain ##EQU7## where the matrices C and G are again symmetric positive semi-definite. 
     In the case of LC circuits, the matrices ζ and A g  are empty. The Laplace transform of the MNA equations (Eq. 4) reduced to 
     
         A.sub.c.sup.T A.sub.c sV.sub.n +A.sub.l.sup.T I.sub.l =A.sub.i.sup.T I.sub.s (s), A.sub.l V.sub.n -LsI.sub.l =0.               (Eq. 11) 
    
     By multiplying the first equation in (Eq. 11) by s and substituting sI L  from the second, we obtain ##EQU8## Again, the matrices C and G are symmetric positive semi-definite and a change of variables σ=s 2  brings (Eq. 12) to the usual form. 
     II. Pade Approximation of Matrix Transfer Functions 
     Recall from (Eq. 7) that the transfer function of a p-port is of the form 
     
         Z(s)=B.sup.T (G+sC).sup.-1 B, sε,                  (Eq. 13) 
    
     where G and C are real symmetric N×N matrices, and B is a real N×p matrix. Note that Z is a matrix-valued rational function: 
     
         Z: ({∞}).sup.p×p.                              (Eq. 14) 
    
     Moreover, each pole s of Z is also an eigenvalue of the matrix pencil G+λC, i.e., each pole s of Z satisfies the equation det (G+sC)=0. 
     A. Review of the Case p=1 
     For the moment, assume that p=1. In this case, Z is a scalar-valued rational function, and for each n, we can define a scalar-valued Pade approximant Z n  to Z as follows. A function of the form ##EQU9## where φ n-1  and ψ n  are polynomials of degree at most n-1 and n, respectively, is called an nth Pade approximant to Z if 
     
         Z(s)=Z.sub.n (s)+θ(s.sup.q(n)),                      (Eq. 16) 
    
     where q(n) is as large as possible. The condition (Eq. 16) means that the Taylor expansions of Z and Z n  about s=0 agree in as many leading Taylor coefficients (the so-called moments) as possible. In general, q(n)=2n. 
     The standard approach to computing Z n  is based on explicit moment generation. First, one computes the leading q(n) Taylor coefficients of Z, and from these, one then generates the coefficients of the polynomials φ n-1  and ψ n  in (Eq. 15). This standard approach to computing Z n  is employed in the asymptotic waveform evaluation (AWE) technique described in L. T. Pillage and R. A. Rohrer, &#34;Asymptotic Waveform Evaluation for Timing Analysis,&#34; IEEE Trans. Computer-Aided Design, vol. 9, pp. 352-366, April 1990; V. Raghavan, R. A. Rohrer, L. T. Pillage, J. Y. Lee, J. E. Bracken and M. M. Alaybeyi, &#34;AWE-inspired,&#34; in Proc. IEEE Custom Integrated Circuits Conference, May 1993. However, computing Pade approximants using explicit moment computations is inherently numerically unstable, and indeed, in practice, this approach can be used for very moderate values of n, such as n≦10. These numerical instabilities can be avoided by exploiting the Lanczos-Pade connection and generating the Pade approximant Z n  via the Lanczos procedure. The resulting procedures for stably computing Pade approximants of transfer functions are PVL for general circuits and its special variant SyPVL for RLC circuits. In SyPVL, the Pade approximant Z n  is computed using the formula 
     
         Z.sub.n (s)=ρ.sub.n.sup.T (A.sub.n +sB.sub.n).sup.-1 ρ.sub.n,(Eq. 17) 
    
     where the n×n matrices A n , B n  and the vector ρ n  are generated from a symmetric version of the Lanczos process. 
     B. The General Case p≧1 
     Returning to the general case of p-ports with p≧1, one approach to obtaining approximations of Z is to compute scalar Pade approximants for each of the p 2  entries of Z by means of p 2  runs of PVL. However, a much more efficient approach is to use the concept of matrix-Pade approximation [Eq. 3] that generates a matrix-valued approximation Z n  for all entries of Z in one run. Moreover, the reduced-order model generated by Z n  is much smaller than the reduced order model obtained from p 2  individual PVL runs. 
     Matrix-Pade approximants can be represented by means of a pair of numerator and denominator matrix polynomials whose coefficients can be obtained via explicit moment computations. However, this approach suffers from the same instability mentioned above for the case p=1. Instead, we use an appropriate extension of the formula (Eq. 17) to the general matrix case p≧1 in our definition of matrix-Pade approximants. More precisely, we say that a matrix-valued function 
     
         Z.sub.n :({∞}).sup.p×p.                        (Eq. 18) 
    
     is an nth matrix-Pade approximant to Z if Z n  is of the form 
     
         Z.sub.n (s)=ρ.sub.n.sup.T (A.sub.n +sB.sub.n).sup.-1 ρ.sub.n,(Eq. 19) 
    
     where ρ n  ε n ×p and A n , B n  ε n ×n 
     and 
     
         Z(s)=Z.sub.n (s)+θ(s.sup.q(n))                       (Eq. 20) 
    
     with maximal q(n). 
     In general, we have q(n)≧2.left brkt-bot.n/p.right brkt-bot. if, and only if, so-called deflation occurs due to certain linear dependencies. 
     III. A Symmetric Lanczos-Type Procedure 
     In this section, the symmetric Lanczos-type procedure is formulated that generates appropriate matrices A n , B n , and ρ n  such that the function Z n  defined in (Eq. 19) is indeed an nth matrix-Pade approximant to Z. This section corresponds to steps 40-60 of FIG. 3, which, in turn, are the details of step 34 of FIG. 2. 
     In order to compute matrix-Pade approximants to matrix-valued transfer functions (Eq. 13), we need a Lanczos-type procedure that can handle multiple starting vectors, namely the p columns of the Matrix B in (Eq. 13). We used a special symmetric variant of the procedure described in R. W. Freund, &#34;An Extension of the Lanczos-Pade Connection to the Matrix Case,&#34; Numerical Analysis Manuscript, Bell Laboratories, Murray Hill, N.J., in preparation, that exploits the symmetry of the matrices G and C in (Eq. 13). Lanczos-type procedures for multiple starting vectors are necessarily quite involved for two reasons. First, in the course of the procedure, some of the generated vectors can become linearly dependent and thus need to be deflated. Second, so-called look-ahead techniques are required to avoid potential breakdowns due to division by quantities that cannot be excluded to be zero. SyMPVL is the only existing Lanczos-type process that has both deflation and look-ahead built in. 
     Next, we give a precise statement of the symmetric Lanczos-type procedure used to compute matrix-Pade approximants to Z. This procedure could be directly stated in terms of G, C, and B, but then we would have to either store additional vectors or abandon the modified Gram-Schmidt orthogonalization used within the procedure. Instead, we assume that G is factored as 
     
         G=MJM.sup.T,                                               (Eq. 21) 
    
     where M, Jε N ×N, 
     and J is a &#34;simple&#34; matrix, such as, e.g., a diagonal matrix. For example, if M is real symmetric positive definite, then we can choose J=I (the N×N identity matrix). A factorization (Eq. 21) can be computed via a suitable version of the Bunch-Parlett-Kaufman procedure if G is indefinite, or a version of the Cholesky procedure if G is symmetric positive definite. We now state the Lanczos-type procedure in terms of M, M T , J, C, and B. This corresponds to rewriting the transfer function (Eq. 13) as 
     
         Z(s)=(M.sup.-1 B).sup.T (J+sM.sup.-1 CM.sup.-T).sup.-1 (M.sup.-1 B),(Eq. 22) 
    
     and applying the Lanczos procedure to the matrix 
     
         F=M.sup.-1 CM.sup.-T                                       (Eq. 23) 
    
     with starting vectors M -1  B, and constructing the Lanczos vectors v 1 , v 2 , . . . , to be J-orthogonal. More precisely, if no look-ahead steps are necessary, the Lanczos vectors satisfy ##EQU10## If look-ahead steps do occur, then the Lanczos vectors are only cluster-wise J-orthogonal, instead of the vector-wise J-orthogonality (Eq. 24). 
     We note that after n steps of the Lanczos-type procedure below, the Lanczos vectors 
     
         v.sub.1, v.sub.2, . . . , v.sub.n                          (Eq. 25) 
    
     have been constructed. In addition, at the end of each nth step, there are p c  =p c  (n) &#34;auxiliary&#34; vectors, 
     
         v.sub.n+1, v.sub.n+2, . . . , v.sub.p.sbsb.c,              (Eq. 26) 
    
     available. These vectors will be turned into Lanczos vectors or deflated in successive iterations. Here, p c  denotes the current block size. Initially, p c  =p, and then within the procedure p c  is reduced by one every time a deflation occurs. 
     Procedure 1 (Symmetric Lanczos-type method with deflation and look-ahead.) 
     INPUT: 
     Matrices G=G T  =MJM T , C=C T  ε N ×N ; 
     A block B=[b 1  b 2  . . . b p  ] ε N ×p of p starting vectors (the vectors b i ). 
     OUTPUT: 
     An p l  ×p matrix ρ; 
     (p l  =p-(# of deflations during the first p iterations).) 
     The non-zero entries of the n×n matrices T n  and Δ n . 
     (Here, n is the value of the iteration counter at termination.) 
     0) For i=1, 2, . . . , p, set 
     
         v.sub.i =J.sup.-1 M.sup.-1 b.sub.i. 
    
     Set p c  =p. 
     (p c  is the current block size.) 
     Set Π v  =.O slashed.. 
     (Π v  records indices of vectors that must be preserved due to inexact deflation. If kεΠ v , then the newly constructed v vectors need to be J-orthogonalized against all vectors v i  with iε k .) 
     Set γ=1, .sub.γ =Π, and V.sup.(γ) =Π. 
     γ is the number of the currently constructed look-ahead cluster, and is the set of indices of the already constructed Lanczos vectors in the γth cluster V.sup.(γ).) 
     For n=1, 2, . . . , do (Build nth Lanczos vector v n .): 
     1) (Perform deflation (if necessary) and obtain v n .) 
     1a) (Decide if ν n  should be deflated.) 
     (i) Set μ=n-p c . 
     (ii) If the current γth cluster contains at least one vector, i.e., ≠.O slashed., then orthogonalize ν n  against the vectors v i  in that cluster: 
     For all iε.sub.γ in ascending order), set 
     (iii) (Deflate if ∥ν n  ∥ is &#34;small&#34;.) ##EQU11## If ∥ν n  ∥&gt;dtol, then continue with step 1b). Otherwise, deflate ν n  by doing the following: 
     (iv) If p c  =1, then stop. 
     (There are no more Krylov vectors.) 
     (v) If μ&gt;0 and the deflated vector ν n  is non-zero, then set ∥ v  =∥ v  ∪{γ(μ)}, and save the vectors v i , iε.sub.γ (μ). 
     (vi) Set p c  =p c  -1. 
     (The current block size is reduced by one.) 
     (vii) For i=n, n+1, . . . , n+p c  -1, set ν i  =ν i+1 . 
     (The auxiliary vector ν n  is deflated. The indices of the remaining auxiliary vectors are reduced by one.) 
     (viii) Repeat step 1a) 
     1b) (Normalize ν n  to obtain the nth Lanczos vector v n .) 
     Set ##EQU12## where t n ,n-p.sbsb.c is a suitable scaling factor, e.g., 
     
         t.sub.n,n-p.sbsb.c =∥v.sub.n ∥.sub.2. 
    
     If n≦p c , then set ρ n ,n-p.sbsb.c +m  =t n ,n-p.sbsb.c. 
     1c) (Update cluster information.) 
     If .sub.γ =.O slashed., set γ v  =γ(max{1,n-p c  }). 
     Set V.sup.(γ) =[V.sup.(γ) v n  ], .sub.γ =.sub.γ ∪{n}, and γ(n)=γ. 
     (γ(n) is the index of the cluster containing the nth Lanczos vector v n .) 
     2) (Compute Δ.sup.(γ) and check for end of look-ahead cluster.) 
     2a) Form 
     
         Δ.sup.(γ) =(V.sup.(γ)).sup.T JV.sup.(γ). 
    
     2b) (Decide if the current cluster is complete.) If the matrix Δ.sup.(γ) is singular (or in some sense &#34;close&#34; to singular), then continue with step 3). Otherwise, the current cluster is complete and the following &#34;end-of-cluster&#34; updates are performed: 
     2c) (J-orthogonalize the vectors ν i , n+1≦i≦n+p c  -1, against the vectors v j  in the current cluster.) 
     For i=n+1, n+2, . . . , n+p c  -1, set 
     
         [t.sub.j,i-p.sbsb.c ].sub.jε.sbsb.γ =(Δ.sup.(γ)).sup.-1 (V.sup.(γ)).sup.T J.sub.νi, 
    
     
         v.sub.i =v.sub.i -ν.sup.(γ) [t.sub.j,i-p.sbsb.c ].sub.jε.sbsb.γ, 
    
     and if i≦p c , set [ρ j ,i-p.sbsb.c +m  ] j ε.sbsb.γ =[t j ,i-p.sbsb.c ] j ε.sbsb.γ. 
     2d) Increment the cluster counter γ=γ+1; 
     Set γ(n+1)=γ, =.O slashed., and V.sup.(γ) =.O slashed.. 
     (The current cluster is complete, and so the vector V n+1  constructed in the next iteration will start a new cluster.) 
     3) (Advance the block Krylov subspace and obtain new vector ν n+p .sbsb.c.) 
     3a) Obtain 
     
         V=J.sup.-1 M.sup.-1 CM.sup.-T v.sub.n 
    
     by first solving the linear system 
     
         M.sup.T t=v.sub.n 
    
     for t, then solving the linear system 
     
         Mv=Ct 
    
     for v, and finally setting v=J -1  v. 
     If the current γth cluster contains at least one vector, i.e., C.sub.γ ≠.O slashed., we may add arbitrary combinations of the vectors in the γth cluster to v: 
     Set 
     
         v=v+V.sup.(γ) [t.sub.i,μ ].sub.iε.sbsb.γ, 
    
      with arbitrary t i ,μ ε. 
     3b) (J-orthogonalize v against previous Lanczos vectors.) 
     (i) (J-orthogonalize v against clusters V.sup.(k), γ v  ≦γ-2.) 
     Set ##EQU13## (δ is the last column of Δ.sup.(γ-1).) 
     For all k=γ v , γ v  +1, . . . , γ-2, set 
     
         [t.sub.j,n ].sub.jε.sbsb.k =(Δ.sup.(k)).sup.-1 ([t.sub.i,j ].sub.iε.sbsb.γ-1,.sub.jε.sbsb.k).sup.T δ, v=v-V.sup.(k) [t.sub.j,n ].sub.jε.sbsb.k. 
    
     (ii) (Perform J-orthogonalization of v due to inexact deflation.) 
     For all k ε∥ v  with k&lt;γ v  (in ascending order), set 
     
         [t.sub.j,n ].sub.jε.sbsb.k =(Δ.sup.(k)).sup.-1 (V.sup.(k)).sup.T Jv, v=v-V.sup.(k) [t.sub.j,n ].sub.jε.sbsb.k. 
    
     (iii) (J-orthogonalize v against cluster V.sup.(γ-1).) 
     Set 
     
         [t.sub.j,n ].sub.jε.sbsb.γ-1 =(Δ.sup.(γ-1)).sup.-1 (V.sup.(γ-1)).sup.T Jv, v.sub.n+p.sbsb.c =v-V.sup.(γ-1) [t.sub.j,n ].sub.jε.sbsb.γ-1. 
    
     To get the nth matrix-Pade approximant Z n  to Z, we need the following quantities from Procedure 1: 
     
         ρ,=[ρ.sub.i,j ].sub.1≦i≦p.sbsb.1.sub., 1≦j≦p, 
    
     
         T.sub.n =[t.sub.i,j ].sub.1≦i≦n, 1≦j≦n, 
    
     
         Δ.sub.n =diag(Δ.sup.(1), Δ.sup.(2), . . . , Δ.sup.(γ)).                                   (Eq. 27) 
    
     The values, , T n  and Δ n  provided from step 34 (FIG. 2) are needed in step 36 for generating the matrix Pade approximant Z n  (s) . 
     Here, for ρ and T n  we use the convention that all the entries ρ ij  and t ij  that are not explicitly defined in Procedure 1 are set to zero. Then ##EQU14## is the nth matrix-Pade approximant to Z generated in step 36 (FIG. 2). For later use, we remark that the matrix T n  can be viewed as a certain oblique projection of the matrix F defined in (Eq. 23) onto the space spanned by the first n Lanczos vectors. More precisely, assuming for simplicity that only zero vectors ν n  =0 are deflated in Procedure 1, we have 
     
         V.sub.n.sup.T JFV.sub.n =Δ.sub.n T.sub.n,            (Eq. 29) 
    
     where 
     
         V.sub.n =[V.sub.1 V.sub.2 . . . V.sub.n ].                 (Eq. 30) 
    
     IV. Stability and Passivity 
     For RC, RL and LC circuits, the method and apparatus generates a stable and passive model, as discussed more fully below. For RLC circuits, Pade-based reduced-order models are in general not stable and not passive. However, if the order n is large enough so that the reduced-order models are sufficiently accurate, then the reduced-order models defined by Z n  are almost stable and passive, and can in fact be made stable and passive by suitable post-processing of Z n . 
     If the reduced-order model is used only for purposes of simulating the linear circuit or subcircuit itself (e.g., to validate the desired frequency response of a circuit design), then the model not being stable and/or passive is not a significant issue. However, if the reduced-order model is used in a larger simulated circuit or intended to be synthesized as an actual circuit, then it should be modified by suitable post-processing to insure stability and passivity. 
     The above-noted post-processing has not been fully been developed at this time and does not form a part of the present invention. However, by way of example, to achieve stability, a pole/residue representation of the matrix-Pade approximant Z n , may be computed. This can be done by computing an eigen decomposition of the matrix T n  in Equation 28. If the reduced-order model is unstable, then at least one of the poles in the pole/residue representation lies in the right half of the complex plane. A first approach to make the model stable is to simply drop the terms in the pole/residue representation of Z n  that correspond to poles in the right half of the complex plane. Let us denote the resulting function by Z n . If the original reduced-order model generated by SyMPVL was accurate enough, then it can be shown that the dropped terms are always small. Still, the modified function Z n  is no longer a matrix-Pade approximant, since dropping the unstable terms will change the leading Taylor coefficients of Z n  and thus the condition of matching leading Taylor coefficients of Z, which defines a matrix-Pade approximant, is no longer satisfied. A second approach is to make up for the dropped terms by modifying the retained stable poles and residues in the pole/residue representation of Z n  in such a way that as many as possible leading Taylor co-efficients of Z n  remain unchanged. The resulting stable reduced-order model is then at least a matrix-Pade-type approximant in the sense that it matches a number of leading Taylor co-efficients of Z, although this number is smaller than the number of co-efficients matched by a true matrix-Pade approximant. In principal, similar approaches can be employed to turn a non-passive reduced-order model into a passive one. 
     A. Stability of RL, RC and LC Circuits 
     Recall from the discussion in Section (I)(B) that the matrices G and C in the formula (Eq. 13) of the matrix transfer function Z are symmetric positive semi-definite for RC and RL circuits. Furthermore, for LC circuits, Z is given by (Eq. 13) with s replaced by s 2 , where G and C are again symmetric positive semi-definite. 
     Since G is symmetric positive semi-definite, it follows that J=I in (Eq. 21). As a result, the Lanczos vectors generated by Procedure 1 are actually orthogonal, and thus Δ n  =I n  (the n×n identity matrix) for all n. Therefore, equation (28) reduces to 
     
         Z.sub.n (s)=ρ.sub.n.sup.T (I.sub.n +sT.sub.n).sup.-1 ρ.sub.n.(Eq. 31) 
    
     Furthermore, the relation (Eq. 29) now becomes 
     
         T.sub.n =V.sub.n.sup.T FV.sub.n.                           (Eq. 32) 
    
     By (Eq. 23), the matrix F has the same eigenvalues as the symmetric positive semi-definite matrix C, and thus all eigenvalues of F are non-negative. Together with (Eq. 32), it follows that T n  is symmetric positive semi-definite. 
     By (Eq. 31), all poles of Z n  are of the form ##EQU15## Since T n  is symmetric positive semi-definite, λ≧0, and thus all poles of Z n  are non-positive. Moreover, it can be shown that a possible pole s=0 is simple. Altogether, this proves that the reduced-order models defined by Z n  are stable for RC and RL circuits. The case of LC circuits can be handled similarly, using the formula (Eq. 12) and the fact that the transformation of s→s 2  maps the purely imaginary poles of an LC circuit into points on the negative real axis. 
     B. Passivity of RL, RC and LC Circuits 
     The reduced-order models defined by Z n  also are passive. Again, we only treat the case of RC and RL circuits. The case of LC circuits can be handled similarly after the transformation s→s 2  has been employed. 
     It is well known that the reduced-order model defined by the matrix-valued rational matrix Z n  given by (Eq. 31) is passive if, and only if, the following three conditions are satisfied: 
     (i) Z n  (s) has no poles in {sε|Re s&gt;0} (the right half of the complex plane); 
     (ii) Z n  (S)=T n  (s) for all sε; 
     (iii) For all sε with Re s&gt;0, 
     
         Re(x.sup.H Z.sub.n (s)x)≧0                          (Eq. 34) 
    
     for all xε p . 
     It thus remains to verify (i)-(iii). Condition (i) is satisfied in view of the stability of T n . Condition (ii) follows immediately from (Eq. 31) and the fact that T n  and ρ n  are real matrices. Finally, with respect to condition (iii), let s be any complex number with Re s&gt;0. Since T n  is symmetric positive semi-definite, we have 
     
         Re(y.sup.H (I.sub.n +sT.sub.n)y)=∥y∥.sub.2.sup.2 +(Re s)y.sup.H T.sub.n y≧0                              (Eq. 35) 
    
     for all yε n . 
     For any given xε p , we set 
     
         y=(I.sub.n +sT.sub.n)ρ.sub.n x.                        (Eq. 36) 
    
     Then yε p , and inserting (Eq. 36) into (Eq. 35) gives 
     
         0≦Re(y.sup.H (I.sub.n +sT.sub.n)y)=Re(x.sup.H ρ.sub.n (I.sub.n +sT.sub.n).sup.-1 ρ.sub.n x)=Re(x.sup.H Z.sub.n (s)x).(Eq. 37) 
    
     Hence the reduced-order model given by Z n  is passive. 
     V. EXAMPLES 
     In this section, numerical results are applied to two example circuits. 
     A. The PEEC Circuit 
     FIG. 4 shows an LC circuit 64 of interest. The details of the circuit 64 are not shown and are not necessary for this high level discussion other than the facts that the circuit consists of only inductors, capacitors, and inductive couplings and is driven by current I x  from a finite impedance source comprising voltage source E and resistor R c . We are interested in computing the response of this circuit; in this case, the current I o  flowing through one of the inductors, L, in the circuit. 
     The equations of the LC two-port are obtained by writing the KCL equations and differentiating them with respect to time gives 
     
         A.sub.l.sup.T sI.sub.l +A.sub.c.sup.T sI.sub.c +asI.sub.x =0.(Eq. 38) 
    
     We substitute the inductor and capacitor equations 
     
         sI.sub.l =L.sup.-1 V.sub.l =L.sup.-1 A.sub.l V 
    
     and 
     
         I.sub.c =sV.sub.c =A.sub.c sV                              (Eq. 39) 
    
     into the differentiated KCL equation (Eq. 38) and obtain the nodal circuit equations 
     
         A.sub.l.sup.T L.sup.-1 A.sub.l V+A.sub.c.sup.T A.sub.c s.sup.2 V+asI.sub.x =0.                                                       (Eq. 40) 
    
     The output of interest is selected among the inductor currents, using the column vector b. This results in the relation ##EQU16## 
     Setting 
     
         G=A.sub.l.sup.T L.sup.-1 A.sub.l 
    
     and 
     
         C=A.sub.c.sup.T A.sub.c, 
    
     and simplifying s, the output of interest can be expressed as ##EQU17## The input impedance of the two-port is defined as the voltage V x  when a unit valued current is applied: ##EQU18## The quantities needed to compute the response of the network are obtained from the two-by-two transfer matrix of the LC two-port 
     
         Z(s)=B.sup.T (G+s.sup.2 C).sup.-1 B,                       (Eq. 44) 
    
     where B denotes the two-column matrix [a l]. Observe that 
     
         Z.sub.in =-sh.sub.11 (s) 
    
     and 
     
         α=-h.sub.21 (s).                                     (Eq. 45) 
    
     SyMPVL is applied to compute the matrix-Pade approximation of this two-by-two transfer function. In this case, the matrix G is singular, or in electrical terms, there is no DC path to the reference voltage (ground) from every node. The transfer function computation requires therefore a frequency shift ##EQU19## In view of (Eq. 28), the approximation will have the form 
     
         Z.sub.n (s)=ρ.sub.n.sup.T (I.sub.n +σT.sub.n).sup.-1 ρ.sub.n.(Eq. 47) 
    
     An approximation of order 50 was needed to obtain a good match of the function as illustrated in FIG. 5. Since the change of variables included the square of s, we have, in fact, a matrix-Pade approximation of order n=50, and we match 2.left brkt-bot.n/p.right brkt-bot.=2.left brkt-bot.50/2.right brkt-bot.=50 matrix moments. Running the procedure six more iterations results in a perfect match. Note that in this example the matrices C and G are real, symmetric and positive definite, and therefore, the approximation always preserves all stability and passivity properties. 
     B. A Package Model 
     As a further example, we analyze a 64-pin package model used for an RF integrated circuit. The integrated circuit is illustrated in FIG. 6 and comprises a die 70, external pins 74 and wirebonds 71 from the die connector pads 73 to the package external pins 74. The package is represented in the figure by dashed line 72. Only eight of the package pins 74, pins 1-4 and 61-64, carry signals, the rest being either unused or necessary for supply voltages. The goal is to characterize the package as a 16-port component (8 exterior and 8 interior terminals) and apply the SyMPVL model reduction on it. The reduced-order model can then replace the large package subcircuit in a circuit simulator resulting in considerable savings in simulation time. 
     The package model is described by approximately 4000 circuit elements, resistors, capacitors, inductors, and inductive couplings. The size of the nodal circuit matrices for the package is about 2000. 
     FIG. 7 shows the voltage-to-voltage transfer function between the external terminal 74 of pin no. 1 and the internal terminal 73 of the same pin. FIG. 8 shows the voltage-to-voltage transfer function between the external terminal 74 of pin no. 1 and the internal terminal 73 of the neighboring pin no. 2. The voltage-to-voltage transfer function between external terminal of pin no. 1 and internal terminal of pin no. 2 illustrated in FIG. 8 is a parasitic interference signal. The plots show results with a reduced model of order 48, 64, and 80 compared with an exact analysis. The results show that reduced-order models can indeed replace the full package subcircuit with little loss of accuracy. The reduction level depends on the desired accuracy. The most accurate of the models that gives an almost perfect match of the frequency response requires only 80 state variables compared to 2000 for the full subcircuit. Since the cost of nonlinear circuit simulation is superlinear in the number of state variables, the computational savings can be significant. 
     As noted above, the application of post-processing techniques to the approximant Z n  (S) can be used to assure stability and passivity. 
     VI. Conclusions 
     We have introduced the present method and apparatus as a tool to compute reduced-order models (based on the matrix-Pade approximation) for large passive linear circuit multi-ports. It is more efficient than the more general MPVL and has the property of generating guaranteed stable and passive models for the special cases of RC, RL and LC circuits. The invention is based on a novel symmetric block-Lanczos procedure. 
     While the reduced-order models generated in accordance with the present invention cannot be guaranteed to be stable and passive for general RLC circuits, they can be made stable and passive using the suitable &#34;post-processing&#34; techniques. 
     Having thus described a few particular embodiments of the invention, various alterations, modifications, and improvements will readily occur to those skilled in the art. Such alterations, modifications and improvements as are made obvious by this disclosure are intended to be part of this description though not expressly stated herein, and are intended to be within the spirit and scope of the invention. Accordingly, the foregoing description is by way of example only, and not limiting. The invention is limited only as defined in the following claims and equivalents thereto.