Abstract:
A new method for MIMO RLCG interconnects model order reduction technique using the global Arnoldi algorithm is proposed that is an extension of the standard Arnoldi algorithm for MIMO systems. Under this framework, the input matrix serves as a stacked vector form and the global Arnoldi algorithm will be the standard Arnoldi algorithm applied to a new matrix pair. This new matrix Krylov subspace from the Frobenius orthonormalization process is the union of system moments. By employing the congruence transformation with this matrix Krylov subspace, the one-sided projection method can be used to construct a reduced-order system. Connections of the reduced system and the original RLCG interconnect circuits are developed. The transfer matrix residual error of reduced system is derived analytically. This error information will be a guideline for the order selection scheme. Experimental results demonstrate the feasibility and the effectiveness of the proposed method.

Description:
BACKGROUND OF THE INVENTION  
       [0001]     1. Field of the Invention  
         [0002]     This invention relates to a method of developing an analogical VLSI macro model in a global Arnoldi algorithm and particularly to that of a fast and precise simplified design on an interconnect circuit model with multiple inputs and outputs.  
         [0003]     2. Description of Related Art  
         [0004]     From the development of CMOS process technology to that of Nano technology, the parasitic effect of an interconnection circuit cannot be ignored, such as IC Interconnect Analysis proposed by M. Celik, L. T. and A. Odabasioglu, Kluwer Academic Publishers in 2002, on the characteristics of interconnection circuit operation. The interconnect circuit is generally represented with mathematics models. Owing to the complexity of a circuit that is gradually going up, that of model order for precisely simulating the circuit also does; thus, a method of effectively reducing the model becomes an essential technology for interconnect circuit reduction and simulation, such as U.S. Pat. Nos. 6,810,506, 6,789,237, 6,687,658, 6,041,170, and U.S. Pat. No. 6,023,573, Krylov-Subspace Methods for Reduced-Order Modeling in Circuit Simulation, Journal of Computational and Applied Mathematics, Vol. 123, pp. 395-421, proposed by R. W. Freund in 2000; and On Projection Based Algorithm for Model Order Reduction of Interconnects, IEEE Trans. on Circuits and System-I: Fundamental Theory and Applications, Vol. 49, No. 11, pp. 1563-1585, proposed by J. M. Wang, C. C. Chu, Q. Yu and E. S. Khu in 2002.  
         [0005]     In the design of high-speed VLSI, an interconnect circuit modeling technology is highly concerned, and several methods have recently been proposed to solve the problems. A prior art is, for example, Asymptotic Waveform Evaluation (AWE), and Arnoldi algorithm, “On Projection-Based Algorithms for Model-Order Reduction of Interconnects,” IEEE Trans. on Circuit and Systems-I: Fundamental Theory and Applications, Vol. 49, No. 11, pp. 1563-1585, proposed by J. M. Wang, C. C. Chu, Q. Yu and E. S. Kuh in 2002. Besides, “Error Estimations of Arnoldi-Based Interconnect Model-Order Reductions,” IEICE Trans. Fundamentals, Vol. E88-A, No. 2, pp. 533-537 was proposed by C. C. Chu, H. J. Lee and W. S. Feng in 2005.  
         [0006]     PVL (Padé via Lanczos), “Efficient Linear Circuit Analysis by Padé Approximation via the Lanczos Process,” IEEE Trans. Comput-Aided Des. Integr. Circuits Syst., Vol. 15, No. 5, pp. 639-649, was proposed by P. Feldmann and R. W. Freund in 1995. Further, “Oblique Projection Methods for Large Scale Model Reduction,         ” SIAM J. Matrix Anal. Appl., Vol. 16, No. 2, pp. 602-627, was proposed by I. M. Jaimoukha and E. M. Kaseally in 1995. Next, “Krylov-subspace methods for reduced-order modeling in circuit simulation,” J. Comput. Appl. Math., vol. 123, pp. 395-421, was proposed by R. W. Freund in 2000.  
         [0007]     However, the prior arts mentioned above only deals with the Signal Input Signal Output (SISO) system; they have not yet dealt with the Multiple Input Multiple Output (MIMO) system.  
         [0008]     Additionally, among the prior arts, MPVL, Reduced-Order Modeling of Large Linear Subcircuits via a Block Lanczos Algorithm,” 32nd ACM/IEEE Design Automation Conference, pp. 474-479, was proposed by P. Feldmann and R. W. Freund.  
         [0009]     Block Arnoldi (BA) algorithm, “Krylov subspace techniques for reduced-order modeling of large-scale dynamical systems,i” Appl. Numer. Math., vol. 43, no. 1-2, pp. 9-44, was proposed by Z. Bai in 2002. Next, Krylov space methods on state-space control models,” Circuits Syst. Signal Process., vol. 13, no. 6, pp. 733-758, was proposed by D. L. Boley in 1994. RIMA: Passive Reduced-Order Interconnect Macromodeling Algorithm,i” IEEE Trans. on Computer-Aided Design of Integrated Circuits and Systems, Vol. 17, No. 8, pp. 645-654, was proposed by A. Odabasioglu, M. Celik and L. T. Pileggi in 1998.  
         [0010]     The MPVL and the block Arnoldi algorithm give the technology of MIMO system model reduction, but when the order of reduced system is higher, the value may not be stable yet.  
         [0011]     The relationship between the reduced circuit and the original circuit is then observed. A transfer function in a frequency domain is generally used to determine whether operation characteristics are consistent. Residual errors of two transfer functions may be regarded as a guideline of the model reduction algorithm to stop iteration process. In the prior arts, Bai used the PVL algorithm to derive a transfer function error E(s) between an original circuit and a reduced circuit. (“Error bound for reduced system model by Pade approximation via the Lanczos process,” IEEE Trans. On Computer-Aided Design of Integrated Circuits and Systems, Vol. 18 pp. 133-141, was proposed by Z. Bai, R. D. Slone, W. T. Smith and Q. Ye in 1999.) However, the expression of error involves the complicated operation of a decomposed matrix (i n −sA) −1  of the original circuit, which is thus not practical to LSI.  
         [0012]     In another prior art, the decomposed matrix (i n −sA) −1  is replaced with a decomposed matrix of reduced circuit, and the PRIMA algorithm is used to get the transfer function error (“Practical considerations for passive reduction of RLC circuits,” Proc. ICCAD, pp. 214-219, proposed by A. Odabasioglu, M. Celik, and L. T. Pileggi in 1999).  
         [0013]     Next, in another prior art, in order to avoid the complicated operation of transfer function error E(s), a residual error E r (s) is given to replace the technology of transfer function E(s) (“Krylov Projection Methods for Model Reduction,” Ph.D. thesis, University of Illinois at Urbana-Champaign, Urbana, Ill., proposed by E. J. Grimme in 1997).  
         [0014]     Consequently, because of the technical defects of described above, the applicant keeps on carving unflaggingly through wholehearted experience and research to develop the present invention, which can effectively improve the defects described above.  
       SUMMARY OF THE INVENTION  
       [0015]     This invention is to solve technical problems of prior arts, such as Asymptotic Waveform Evaluation, Arnoldi algorithm, PVL and so on, which only deal with a Signal Input Signal Output (SISO) system but not deal with a Multiple Input Multiple Output (MIMO) system. The MPVL and the block Arnoldi algorithm give the technology of MIMO system model reduction, but when the order of reduced system is higher, the value may not be stable yet.  
         [0016]     For solution to the problems above, this invention provides a method of developing an analogical VLSI macro model in a global Arnoldi algorithm, comprising:  
         [0017]     step 1 to input a net-shaped circuit;  
         [0018]     step 2 to input a frequency expansion point;  
         [0019]     step 3 to build up a state-space matrix for a circuit;  
         [0020]     step 4 to generate a projection matrix by way of a global Arnoldi algorithm;  
         [0021]     step 5 to determine a reduced model order by an iteration termination condition, and to execute a first model reduction; and  
         [0022]     step 6 to build up a mathematics model for a perturbation system and to execute a second model reduction.  
         [0023]     From the description above, under the iteration termination condition at step 5, the residual error for the first reduced system and the original system may be defined to be: 
 
 E   r ( s )= s ( V   g,q   V   g,q   +   −I   n ) h   q+1,q   g   V   g,q+1   E   q   T ( I   qs   −s ( H   g,q             I   s )− sV   g,q   +   ΔV   g,q )
 
∥E r (s)∥ ∞ ≦κ(S qs )∥V g,q V g,q   + −I n ∥ 2 |h q+1,q   g |∥V g,q   + ∥ ∥R∥ 2 
 
 may be given when a norm is derived from E r (s); a value in the iteration process h q+1,q   g  may technically serve to evaluate the reduced model order; thus, an order q is determined to satisfy  
           μ   q     =              h     q   ,     q   -   1       g       h       q   +   1     ,   q     g            &lt;   ɛ       ,       
 
 where ε is a permissible error that is small enough. 
 
         [0024]     where for the second model reduction described at step 6, the perturbation system is added to serve as the perturbation of additive property for the transfer function H(s) of original circuit, and the transfer function H(s) of a corrective node analysis may be indicated below as:  
         A   ⁢       ⅆ       x   Δ     ⁡     (   t   )           ⅆ   t         =           x   Δ     ⁡     (   t   )       +     R   ⁢           ⁢     u   ⁡     (   t   )       ⁢           ⁢   and   ⁢           ⁢     y   ⁡     (   t   )           =     C   ⁢           ⁢       x   Δ     ⁡     (   t   )               
 
         [0025]     where Δ=h q+1,q   g V g,q+1 E q   T V g,q  and q are reduced model orders in a global Arnoldi algorithm, h q+1,q   g  and V g,q+1  may be given in the process of operation of the reduction system, V g,q   +  is a virtual inverse matrix of a projection matrix V g,q , and thus the transfer function H Δ (s) of perturbation system is equal to a transfer function Ĥ(s) of the system that is reduced.  
         [0026]     Compared with the prior arts for the effects, the global Arnoldi algorithm according to this invention may be regarded as an extension of conventional SISO Arnoldi algorithm. In the algorithm, Krylov subspace generated in Frobenius orthonormalization process should be used and is actually a transformation from a conventional method, in which an input matrix may be regarded as a stacked vector form, namely the union of original system moments. By employing the congruence transformation with the matrix Krylov subspace generated in the global Arnoldi algorithm according to this invention, the one-sided projection method can be used to construct a reduced model system; in comparison with the model reduction skill of the current block Arnoldi algorithm, it proves that the transfer functions of two reduced system are identical and the complicated calculation of global Arnoldi algorithm is easier than that of the conventional block Arnoldi algorithm. This invention provides a residual error relation for the reduced system and the original system and error formulae on which the order determination of reduced circuit is based. Further, in this invention, a math expression of the MIMO circuit perturbation system. It proves that the transfer matrix in a two-order reduced system is corresponding to the transfer function of the added perturbation matrix in the original system.  
         [0027]     However, in the description mentioned above, only the preferred embodiments according to this invention are provided without limit to this invention and the characteristics of this invention; all those skilled in the art without exception should include the equivalent changes and modifications as falling within the true scope and spirit of the present invention. 
     
    
     BRIEF DESCRIPTION OF THE DRAWINGS  
       [0028]      FIG. 1  is a flow chart of this invention.  
         [0029]      FIG. 2  is a schematic view illustrating an embodiment of test on an interconnection circuit with 2 inputs and 2 outputs.  
         [0030]      FIG. 3  is a schematic view illustrating the number of order of a parameter-determined reduced model in the process of algorithm iteration according to this invention.  
         [0031]      FIG. 4  is a schematic view illustrating an analysis on errors between the system union of three reduced models according to this invention and the system union of an original system.  
         [0032]      FIG. 5  is a schematic view illustrating the frequency response of a first reduced model and a reduced model of block Arnoldi algorithm according to this invention.  
         [0033]      FIG. 6  is a schematic view illustrating the frequency response of a second reduced model and an original system of added perturbation system according to this invention. 
     
    
     DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS  
       [0034]     Now, the present invention will be described more specifically with reference to the following embodiments. It is to be noted that the following descriptions of preferred embodiments of this invention are presented herein for purpose of illustration and description only; it is not intended to be exhaustive or to be limited to the precise form disclosed.  
         [0035]     In a prior art, a linear, time-invariant, RLCG interconnect circuit in VLSI is generally represented in the following Modified Nodal Analysis (MNA) formula:  
                 M   ⁢       ⅆ     x   ⁡     (   t   )           ⅆ   t         +     N   ⁢           ⁢     x   ⁡     (   t   )         +     B   ⁢           ⁢     u   ⁡     (   t   )           =       0   ⁢           ⁢   and   ⁢           ⁢     y   ⁡     (   t   )         =     C   ⁢           ⁢       x   ⁡     (   t   )       .                 (   1   )             
 
         [0036]     where  
       M   =     [         C       0           0       L         ]         
 
 comprises a capacitor and an inductor, and  
       N   =     [         0       E             -     E   T           R         ]         
 
 comprises a resistor and satisfies the incidence matrix E in Kirchhoff&#39;s voltage and current laws. M, N ε R n×n  and B ε R n×s  represents nodes of voltage sources applied, in which s is the number of voltage sources. C ε R k×n  represents a node that measures impulse response, in which k represents a single measured point. x(t) represents a system union function, comprising voltage union and current union, namely  
           x   ⁡     (   t   )       =     [           v   ⁡     (   t   )                 i   ⁡     (   t   )             ]       ,       
 
 while u(t) represents a function of system input signal. For A=−(N+s 0 M) −1 M and R=(N+s 0 M) −1 B, s 0  represents a frequency expansion point for selection, assuming that N+s 0 M is a non-singular matrix. Then, equation (1) may be changed into  
               A   ⁢       ⅆ     x   ⁡     (   t   )           ⅆ   t         =         x   ⁡     (   t   )       +     R   ⁢           ⁢     u   ⁡     (   t   )       ⁢           ⁢   and   ⁢           ⁢     y   ⁡     (   t   )           =     C   ⁢           ⁢     x   ⁡     (   t   )                   (   2   )             
 
         [0037]     The reduced model is used to find a smaller status space and provided for a designer to analyze the original system in limited hours through equivalent, original system. For the reduced system (MNA), equation (2) may be changed to:  
                 A   ^     ⁢       ⅆ       x   ^     ⁡     (   t   )           ⅆ   t         =           x   ^     ⁡     (   t   )       +       R   ^     ⁢           ⁢     u   ⁡     (   t   )       ⁢           ⁢   and   ⁢           ⁢     y   ⁡     (   t   )           =       C   ^     ⁢       x   ^     ⁡     (   t   )                   (   3   )             
 
         [0038]     where for Â ε R q×q  and Â ε R q×q , q represents the order of reduced system. This invention is briefly described for s=k.  
         [0039]     For the MIMO system, a standard Arnoldi algorithm must be corrected. The BA algorithm is a known method of solving MIMO and generates an orthonormalization base κ b (A,V b,0 )={V b,0  AV b,0  . . . A q−1 V b,0 } in the Krylov subspace by means of recursion, and generates a projection matrix V b,q =[V b,1  V b,2  . . . V b,q ]ε R n×qs  ε K b (A,V b,1 ), in which V b,1  is given from a matrix R through QR. The original system A in the MIMO system may be reduced to a smaller Upper Hessenberg Matrix. A pseudo code for the BA algorithm is shown below.  
                                                                                           Algorithm : (input :A ,R,q; output:V q   b , H q   b )           (1): /* Initialize */           Factor V 0 R=QR factorization of R           V q   b  = [V 1 ]           (2): /* Iteration */           for K = 1,2,···,q                V k   (0)  = AV k−1             for j= 1,2,···,k                H k−j,k−1   b  = V k−j   T V k   (j−1)             V k   (j)  = V k   (j−1)  − v k−j H k−j,k−1   b                  end for                Factor (V q   b ) k (H q   b ) k,k−1  = QR factorization of V k   (k)             end for                      
 
         [0040]     In the process of BA algorithm iteration, the relation is: 
 
 AV   b,q   =V   b,q   H   b,q   +V   b,q   h   q,q−1   b   E   q   T    (4) 
 
         [0041]     where  
         E   q     =         [           0   s         ⋯         0   s           I     s   ×   s             ]     T     ∈       R     qs   ×   s       ⁢           ⁢   and           
         H     b   ,   q       =       [           h   11   b           h   12   b         ⋯       ⋯         h     1   ⁢           ⁢   q     b               h   21   b           h   22   b         ⋯       ⋯         h     2   ⁢           ⁢   q     b             0         h   32   b           h   33   b         ⋯       ⋮           ⋮       ⋰       ⋰       ⋰       ⋮           0       ⋯       0         h     q   ,     q   -   1       b           h     q   ⁢           ⁢   q     b           ]     ∈     R     q   ⁢           ⁢   s   ×   q   ⁢           ⁢   s             
 
         [0042]     H b,q  is an upper Hessenberg matrix, and the relation between H b,q  and A is: 
 
H b,q =V b,q   T AV b,q    (5) 
 
         [0043]     where V b,q   T  V b,q =I qs×qs .  
         [0044]     If A j−1  is multiplied at the left side of equation (4) and A j−1  is multiplied at the right side of equation (4), then 
 
A j   V   b,1 =V b,q H b,q   j E q , ∀ j =0,1 . . . ,q−1.   (6) 
 
 where for E 1 └I s×s  0 s  . . . 0 s ┘ε R qs×s , V b,1  is an initial matrix in the Krylov subspace, and V b,1 =V b,q E 1 . 
 
         [0045]     For a non-singular matrix, the system union is defined to X (j) (s 0 )=A j R=A j V b,1 G=V b,q H b,q   j E 1 G. For j=0, . . . , q−1, X (j) (s 0 ) relates to the matrix V b,q . Thus, the reduced system may be defined to: 
 
Â=V b,q   T AV b,q =H b,q , {circumflex over (R)}=V b,q   T R, Ĉ=CV b,q    (7) 
 
         [0046]     For j=0, . . . ,q−1 the system union of reduced system may be corresponding to that of original system:  
                   Y   ^       (   j   )       ⁡     (     s   o     )       =       C   ^     ⁢       A   ^     j     ⁢     R   ^                   =     C   ⁢           ⁢     V     b   ,   q       ⁢     V     b   ,   q     T     ⁢     A   j     ⁢     V     b   ,   q       ⁢     V     b   ,   q     T     ⁢   R                 =       C   ⁢           ⁢     A   j     ⁢   R     =       Y     (   j   )       ⁡     (     s   o     )                   
 
         [0047]     where for the property of applied matrix operation, if V T V=I and V T V=I then V T V=I.  
         [0048]     What is disclosed in this invention is the method of developing the analogical VLSI macro model in the global Arnoldi algorithm, in which the global Arnoldi (GA) used in this invention is provided with a model reduction skill for the reduced system formed in the MIMO system. A pseudo code for the GA algorithm (proposed at 2005 by K. Jbilou, A. J. Riquet, “Projection Methods for Large Lyapunov Matrix Equations,” Linear Algebra and its Applications, to appear) is shown below:  
                                                                                                               Algorithm   ⁢     :      ⁢     (       input   ⁢     :      ⁢   A     ,   R   ,     q   ;     output   ⁢     :     ⁢           ⁢     V   q   g         ,     H   q   g       )                                 (1): /* Initialize */                                 V   0     =     R   /          R        F                                             V   q   g     =         [           ]     ⁢           ⁢   and   ⁢           ⁢     V   q   g       =     [     V   1     ]                                   (2): /* Iteration */           for j = 1, 2, . . . , q                                                  V   ~     =     AV   j   g                                 for j = 1, 2, . . . , j                                                  h     i   ,   j       =       〈       V   i     ,     V   ~       〉     F                                           V   ~     =       V   ~     -       H     i   ,   j     g     ⁢     V   i   g                                          end for                                 H       j   +   1     ,   j     g     =            V   ~          F                                           V     j   +   1     g     =       V   ~       h       j   +   1     ,   j                                                 V   q   g     =           [     V   q   g             V     j   +   1     g             ]                                    end for                      
 
         [0049]     where for V g,1  matrix,  
         V     g   ,   1       =       1          R        F       ⁢   R         
 
 is an initial matrix in the Krylov subspace. In the GA algorithm, from K q (A,R)=span{R, AR, . . . , A q−1 R} in the Krylov subspace K q (A,R)=span{R, AR, . . . , A q−1 R}, a Frobenius orthonormalization base V g,1  V g,2  . . . V g,q  is given, and the following properties are given: 
 
For i≠j; i,j=1,2, . . . ,q, &lt;V g,i , V g,j &gt; F =0   (8) 
 
For i=j, &lt;V g,i , V g,j &gt; F =1   (9) 
 
 where &lt;.,.&gt; F  is a Frobenius inner product, namely &lt;A,B&gt; F =trace(A T B)=vec(A) T  vec(B), in which if A=[A *1  A *2  A *3  . . . A *n ]ε R m×m  and A *j  ε R m , then vec(A)=[A *1   T  A *2   T  . . . A *n   T ] T  ε R mn  and it is called vectorization of A, which is a long vectoring formed by a row vector of stack A. When the matrix of R m×n , A 1 ,A 2 , . . . ,A k , is linearly independent and the matrix of R mn , vec(A 1 ),vec(A 2 ), . . . , vec(A k ), is linearly independent, the relation between the vectorization and Kronecker product may be found in a prior are (proposed in 1985 by P. Lancaster and M Tismenetsky, The Theory of Matrices: with Applications, Academic Press, pp. 410), conclusion is made below: 
 
vec(ABC)=(C T            A)vec(B).   (10) 
 
vec(A) T vec(B)=trace(A T B).   (11) 
 
(A         B)(C         D)=(AC×BD).   (12) 
 
 Further, a Frobenius norm is defined to ∥A,B∥ F =√{square root over (|trace(A T B)|)}=√{square root over (vec(A) T vec(B))} (proposed in year 1985 by P. Lancaster and M Tismenetsky, The Theory of Matrices: with Applications, Academic Press). 
 
         [0050]     From the GA algorithm, for i=1,2 . . . ,q of the matrix V g,j , the interdependent property between the column vectors does not impact the algorithm. In fact, when the matrix Krylov subspace is used, the Frobenius orthonormalization base {V g,1 , V g,2 , . . . , V g,q } may be given in the GA algorithm, namely, if i≠j, trace((V g,i ) T V g,j )=0 in the matrix Krylov subspace; each column of matrices given in the BA algorithm are orthonormal with each other, indicating that, in the space of real matrix, when the Frobenius orthonormalization base in the Krylov subspace is derived in the BA algorithm, the orthonormalization base in the block Krylov subspace of R n  may be given in the BA algorithm, which is a differentia from the GA and BA algorithms.  
         [0051]     In the GA algorithm, V g,q  is made to be a n×qs matrix and H g,q  is made to be q×q an upper Hessenberg matrix, satisfying the following relation: 
 
 AV   g,q   =V   g,q ( H   g,q             I   s )+ h   q+1,q   g   V   g,q+1   E   q   T ,   (13) 
 
 where  
           H     g   ,   q       =       [           h   11   g           h   12   g         …       …         h     1   ⁢   q     g               h   21   g           h   22   g         …       …         h     2   ⁢   q     g             0         h   32   g           h   33   g         …       ⋮           ⋮       ⋰       ⋰       ⋰       ⋮           0       ⋯       0         h     q   ,     q   -   1       g           h   qq   g           ]     ∈       •     q   ×   q       .         ,       
 
         [0052]     where           is the product of Kronecker, and the products of A=[a ij ] i,j=1   m ε R m×m , B=[b ij ] i,j=1   n ε R n×n , and A and B Kronecker is made to be A{circle around (×)} B ε R mn×mn , being defined below to:  
         A   ⊗   B     =       [             a   11     ⁢   B             a   12     ⁢   B         …           a     1   ⁢   m       ⁢   B                 a   21     ⁢   B             a   22     ⁢   B         …           a     2   ⁢   m       ⁢   B             ⋮       ⋮       ⋮       ⋮               a     m   ⁢           ⁢   1       ⁢   B             a     m   ⁢           ⁢   2       ⁢   B         …           a     m   ⁢           ⁢   m       ⁢   B           ]     =       [       a   ij     ⁢   B     ]       i   ,     j   =   1       m           
 
         [0053]     V g,q  is made to be a matrix, being defined to V g,q =[V g,1  V g,2  . . . V g,q ], where matrices V g,1 , V g,2 , . . . , V g,q  are given in the GA algorithm.  
         [0054]     Here, a vector value function is defined to close relate to the relevant matrices and Kronecker. For the relation between the system union matrix X (j) (s 0 )=A j R and vectorization, vec(A j R)=(I s           A j )vec(R) may be derived from equation (10), since all the matrices may be regarded as stacked vectors. The inner product is corrected at equation (11). The GA algorithm is quite similar to the standard Arnoldi algorithm, but the standard inner product is replaced by equation (11).  
         [0055]     In the GA algorithm, in case of h g+1,j   g =0, the algorithm stops until time j of iteration, but the same expanded subspaces may still be given. Opposite to the GA algorithm, breakdown more often occurs in the BA algorithm, which is discussed in many prior arts, and the method may be given the breakdown result better than the BA algorithm does. (proposed in year 2000 by Z. Bai, J. Demmel, J. Dongarra, A. Ruhe, H. van der Vorst, and editors, Templates for the Solution of Algebraic Eigenvalue Problems: A Practical Guide, SIAM, Philadelphia and proposed in year 1992 by Y. Saad, Numerical Methods for Large Eigenvalue Problems,” Manchester University Press).  
         [0056]     The equation (13) may be further reduced to: 
 
 A   j   V   g,1   =V   g,q ( H   g,q   j   e   1             I   s )   (14) 
 
         [0057]     For V g,1 =V g,q E 1  and h q+1,q   g V g,q+1 E q   T E 1 =0, E 1  is multiplied at the right side of equation (13) and E 1  is derived. Assuming that the equation (14) comes into existence for, in case of j=i+1,  
                 A     i   +   1       ⁢     V     g   ,   1         =       AV     g   ,   q       ⁡     (       H     g   ,   q       ⁢       e   1     ⊗     I   s         )                   =       (         V     g   ,   q       ⁡     (       H     g   ,   q       ⊗     I   s       )       +       h       q   +   1     ,   q     g     ⁢     V     g   ,     q   +   1         ⁢     E   q   T         )     ⁢     (       H     i   ,   q       ⁢       e   1     ⊗     I   s         )                   =         V     g   ,   q       ⁡     (       H       i   +   1     ,   q       ⁢       e   1     ⊗     I   s         )       +       h       q   +   1     ,   q     g     ⁢     V     g   ,     q   +   1         ⁢       E   q   T     ⁡     (       H     i   ,   q       ⊗     I   s       )       ⁢     E   1                   
 
         [0058]     For all the matrices Z and q≠1, then E q   T ZE 1 =0 exists. 
 
 A   j   V   g,1   =V   g,q ( H   g,q   j   e   1             I   s ), ∀ j 32 0,1 , . . . ,q− 1.   (15) 
 
 To create the reduced system, the relation between the system union and the Krylov subspace is developed. In this invention, two reduced systems is analyzed, and the two reduced systems maybe given by means of congruence transformation and union matching of order q may be achieved. 
 
         [0059]     In this invention, the two GA algorithms used to bring the reduced model is proposed. The first reduced system is defined to: 
 
Â g,1 =V g,q AV g,q , {circumflex over (R)} g,1 =V g,q R, and Ĉ g,1 =CV g,q ,   (16) 
 
         [0060]     where V g,q   + is a virtual inverse matrix of V g,q  and defined to V g,q   + =(V g,q   T V g,q ) −1 V g,q   T ,  
         [0061]     where the reduced system Â may be further reduced to:  
                 A   ^       g   ,   1       =       V     g   ,   q     +     ⁢     AV     g   ,   q                     =       V     g   ,   q     +     ⁡     (         V     g   ,   q       ⁡     (       H     g   ,   q       ⊗     I   s       )       +       h       q   +   1     ,   q     g     ⁢     V     g   ,     q   +   1         ⁢     E   q   T         )                   =       (       H     g   ,   q       ⊗     I   s       )     +       V     g   ,   q     +     ⁢   Δ   ⁢           ⁢     V     g   ,   q                     
 
         [0062]     where Δ is defined to Δ=h q+1,q   g V g,q+1 E q   T V g,q   + . Also, {circumflex over (R)} may be reduced to: 
 
 {circumflex over (R)}   g,1   =V   g,q   +   R=V   g,q   +   ∥R∥   F   V   g,q ( e   1             I   s )=∥ R∥   F  ( e   1   I   s ) 
 
         [0063]     By means of the reduction technology developed in the GA algorithm, the property of union may be achieved.  
         [0064]     In case of j=0,1, . . . , q−1, the output union Ŷ (j) (s 0 ) of reduced system that is given in the global Arnoldi algorithm is equal to the Y (j) (s 0 ) union output of original MNA system in equation (2), namely  
                   Y   ^       (   j   )       ⁡     (     s   0     )       =         C   ̑       g   ,   1       ⁢       A   ^       g   ,   1     j     ⁢       R   ^       g   ,   1                     =       CV     g   ,   q       ⁢     V     g   ,   q     +     ⁢     A   j     ⁢     V     g   ,   q       ⁢     V     g   ,   q     +     ⁢   R                 =       CA   j     ⁢   R                 =       Y     (   j   )       ⁡     (     s   0     )                 
 
         [0065]     Provided in a prior art (proposed in year 2004 by K. Gallivan, A. Vandendorpe and P. Van Dooren,           §Sylvester Equations and Projection-based. Model Reduction,         ”  J. Computational and Applied Mathematics,  Vol. 162, pp. 213-229), the reduction transformation matrix is defined to Ĥ(s), and Ĥ(s), Ĥ(s) and Ĥ(s). If projection V and Z is replaced with other matrices {circumflex over (V)}=VR and {circumflex over (Z)}=ZL in a space similarly expanded and R and L may be oppositely transformable, the projection transformation matrix is not changed.  
         [0066]     Thus, it is identical, proving that, if the system union is expressed as {X 0 (s 0 ),X 1 (s 0 ), . . . X q−1 (s 0 )}, the transformation matrix of reduced system is given in the BA algorithm and that is given in the GA algorithm. According to the equation (9), the system union given in the BA algorithm may be expressed as:  
         [       E   1     ⁢     V     b   ,   q       ⁢     E   1     ⁢     V     b   ,   q       ⁢     H     b   ,   q       ⁢     E   1     ⁢           ⁢   …   ⁢           ⁢     V     b   ,   q       ⁢     H     b   ,   q       q   -   1         ]     =         V     b   ,   q       ⁡     [       E   1     ⁢     H     b   ,   q       ⁢     E   1     ⁢           ⁢   …   ⁢           ⁢       (     H     b   ,   q       )       q   -   1       ⁢     E   1       ]       =       V     b   ,   q       ⁢     ℜ   b             
 
         [0067]     H b,q  is an upper Hessenberg matrix, so           is an triangle matrix on a block.  
         [0068]     On the other hand, according to the equation (15), the system union given in the GA algorithm may be expressed as:  
         [       V     g   ,   q       ⁢     E   1     ⁢       V     g   ,   q       ⁡     (       H     g   ,   q       ⊗     I   s       )       ⁢     E   1     ⁢           ⁢   …   ⁢           ⁢       V     g   ,   q       ⁡     (       H     g   ,   g       q   -   1       ⊗     I   s       )       ⁢     E   1       ]     =         V     g   ,   q       ⁡     [         E   1     ⁡     (       H     g   ,   q       ⊗     I   s       )       ⁢     E   1     ⁢           ⁢   …   ⁢           ⁢     (         (     H     g   ,   q       )       q   -   1       ⊗     I   s       )     ⁢     E   1       ]       =       V     g   ,   q       ⁢     ℜ   g             
 
         [0069]     where           is also an upper triangle matrix, since  
           colsp   ⁡     [         X   0     ⁡     (     s   0     )       ,       X   1     ⁡     (     s   0     )       ,     …   ⁢           ⁢       X     q   -   1       ⁡     (     s   0     )           ]       =       colsp   ⁡     (     V     b   ,   q       )       =         colsp   ⁡     (     V     g   ,   q       )       ·     V     b   ,   q         =       V     g   ,   q       ⁢         ℜ   g     ⁡     (     ℜ   b     )         -   1               ,         ℜ   g     ⁡     (     ℜ   b     )         -   1           
 
 is an triangle matrix and non-singular. Thus, it proves that the two transfer matrix are the same. 
 
         [0070]     An transfer matrix error between the original MNA and reduced systems is not easily analyzed and given, so for the conventional SISO system, the prior art here gives the difference of the two systems from the conception of a residual error (proposed in year 2005 by C. C. Chu, H. J. Lee and W. S. Feng, “Error Estimations of Arnoldi-Based Interconnect Model-Order Reductions,”  IEICE Trans. Fundamentals,  Vol. E88-A, No. 2, pp. 533-537). From the definition of MIMO system in this invention, the residual error E r (s) is 
 
 E   r ( s )=( I   n   −sA ) {tilde over (X)} ( s )− R    (17) 
 
         [0071]     where {tilde over (X)}(s) is an approximate solution to {tilde over (X)}(s). If {tilde over (X)}(s)=X(s), then {tilde over (X)}(s)=X(s). In the BA algorithm or the GA algorithm, the approximate solution {tilde over (X)}(s) must belong to the Krylov subspace, namely {tilde over (X)}(s)=V g,q {tilde over (X)}(s) or V b,q {circumflex over (X)}(s).  
         [0072]     If the GA algorithm executes q times of iteration, a Frobenius orthonormalization matrix V g,q  and a corresponding upper Hessenberg matrix H g,q  may be given. {tilde over (X)}(s) is made to be an approximate solution to X(s), while {tilde over (X)}(s) is made to be that after q times of the iteration operation in the Arnoldi algorithm, namely X(s)=V g,q {circumflex over (X)}(s), in which E r (s) is a residual error. The calculation of residual error E r (s) is shown as follows.  
         [0073]     Because of {tilde over (X)}(s) ε κ q (A,R), {tilde over (X)}(s) may be expressed as the linear combination of column vector V g,q , namely {tilde over (X)}(s)=V g,q {circumflex over (X)}(s). The residual error may be given from the following operation:  
                       E   r     ⁡     (   s   )       =         (       I   n     -   sA     )     ⁢     V     g   ,   q       ⁢       X   ^     ⁡     (   s   )         -   R                 =         (       V     g   ,   q       -     sAV     g   ,   q         )     ⁢       (       I   qs     -     s   ⁡     (       H     g   ,   q       ⊗     I   s       )         )       -   1       ⁢     V     g   ,   q     +     ⁢   R     -   R                 =         (               V     g   ,   q       ⁡     (       I   qs     -     s   ⁡     (       H     g   ,   q       ⊗     I   s       )         )       -                 sh       q   +   1     ,   q     g     ⁢     V     g   ,     q   +   1         ⁢     E   q   T             )     ⁢       (       I   qs     -     s   ⁡     (       H     g   ,   q       ⊗     I   s       )         )       -   1       ⁢       (     V     g   ,   q       )     +     ⁢   R     -   R                   (   18   )             
 
         [0074]     R belongs to the expanded subspace of V g,q , so V g,q (V g,q ) + R=R proves. Through the simple algebra operation, the above equation may be changed into: 
 
 E   r ( s )=− sh   q+1,q   g   V   g,q+1   E   q   T ( I   qs   −s ( H   g,q             I   s )) −1 ( V   g,q ) +   R    (19) 
 
         [0075]     An error range may be estimated from the following means. Assuming that all the properties of H g,q  are very simple and (H g,q           I s )=S qs Λ qs   S   qs   −1  is the property value resolution of (H g,q           I s ), the equation (19) being reduced to:  
                       E   r     ⁡     (   s   )       =       -     sh       q   +   1     ,   q     g       ⁢     V     g   ,     q   +   1         ⁢     E   q   T     ⁢         S   qs     ⁡     (       I   qs     -     s   ⁢           ⁢     Λ   q         )         -   1       ⁢     S   qs     -   1       ⁢     V     g   ,   q     +     ⁢   R                 =       sh       q   +   1     ,   q     g     ⁢     V     g   ,     q   +   1         ⁢     E   q   T     ⁢     S   qs     ⁢     Z   ⁡     (   s   )       ⁢     S   qs     -   1       ⁢     V     g   ,   q     +     ⁢   R                   (   20   )             
 
         [0076]     where  
           Z   ⁡     (   s   )       =         diag   ⁡     [     s     1   -     s   ⁢           ⁢     λ   i           ]       qs     .       ,       
 
 since Z(s) is a high-pass matrix,  
                Z   ⁡     (   s   )            ∞     =       min     i   =   1     qs     ⁢       1                ⁢     λ   i            .           
 
 After a norm L ∞  is employed from the two sides of equation (20), the following equation is given: 
 
∥ E ( s )∥ ∞ ≦κ( S   qs )| h   q+1,q   g |∥( V   q   g ) +   ∥∥R∥   2    (21) 
 
         [0077]     where κ(·) is the status number of matrix. From the above error estimation, only κ(S qs ), V qs   +  R, and h q+1,q   g  are included. Compared with the representation of error in prior art (proposed in year 1999 by A. Odabasioglu, M. Celik and L. T. Pileggi, practical Considerations for Passive Reduction of RLC Circuits,         ”  Proc. ICCAD , pp. 214-219), the suggested formula is involved in less calculation. Since it wastes time in κ(S qs ) calculation, only h q+1,q  may be considered among candidate systems. Instead of an absolute value, a relative value  
         μ   q     =            h     q   ,     q   -   1       g       h       q   +   1     ,   q     g                
 
 is used to serve the basis of iteration process termination. If μ q  is quite small, then the original system and the reduced system are almost equal. 
 
         [0078]     From the second reduced system according to this invention, in comparison of the given system given in the BA algorithm (10) with that given in the GA algorithm (18), it is apparent that the reduced system Â is similar; except the term, related to Δ, added in the GA algorithm, the column vector of matrix V g,q  does not need interdependent orthonormalization. To keep the simple formula of matrix Â in the BA algorithm, the second reduced system is here defined to: 
 
 Â   g,2   =V   g,q   + ( A −Δ) V   g,q    (22) 
 
         [0079]     where Δ=h q+1,q   g V g,q+1 E q   T V g,q . From the equation (12), it is known that the reduced Â may be further reduced to: 
 
Â g,2 =H q           I s  
 
         [0080]     where in case of i=0,1, . . . ,q−1, the union matching property still exists.  
                         C   ̑       g   ,   2       ⁢       A   ̑       g   ,   2     i     ⁢       R   ̑       g   ,   2         =           CV     g   ,   q       ⁡     (       H     g   ,   q       ⊗     I   s       )       i     ⁢          R        F     ⁢     (       e   1     ⊗     I   s       )                   =         CV     g   ,   q       ⁡     (       H     g   ,   q     i     ⁢       e   1     ⊗     I   s         )       ⁢          R        F                   =       CA   i     ⁢     V     g   ,   1       ⁢          R        F                   =       CA   i     ⁢   R                   (   23   )             
 
         [0081]     Thus, the transfer matrix of reduced system may be changed into: 
 
 Ĥ   g,2 ( s   0 +σ)= Ĉ   g,2 ( I   q   −σÂ   g,2 ) −1   {circumflex over (R)}   g,2    (24) 
 
         [0082]     Next, an equation of the status of a perturbation system that is defined in this invention is:  
               A   ⁢       ⅆ       x   Δ     ⁡     (   t   )           ⅆ   t         =           x   Δ     ⁡     (   t   )       +       Ru   ⁡     (   t   )       ⁢           ⁢   and   ⁢           ⁢     y   ⁡     (   t   )           =       Cx   Δ     ⁡     (   t   )                 (   25   )             
 
         [0083]     The transfer matrix of perturbation system is made to be H Δ (s) The transfer matrix of reduced system defined in the equation (22) is actually equal to the transfer matrix of original MNA provided with perturbation system defined in the equation (25), namely Ĥ g,2 (s 0 +σ)=H Δ (s 0 +σ).  
         [0084]     Since Â g,2 =V g,q   + (A−Δ)V g,q  may be changed into (A−Δ)V g,q =V g,q Â g,2 , −σ is multiplied and then plus V q  at the two sides of equation, thereby the equation being changed into: 
 
 V   g,q ( I   qs   −σÂ   g,2 ) −1 =( I −σ( A −Δ)) −1   V   g,q  
 
         [0085]     Finally, C is multiplied at the left side of equation and C is multiplied at the right side of equation, and then the result is: 
 
 CV   g,q ( I   qs   −σÂ)   −1   ∥R∥   F ( e   1             I   s )= C ( I   n −σ( A −Δ)) −1   V   g,q   ∥R∥   F ( e   1           I s )   (26) 
 
         [0086]     Since Ĉ g,2 =CV g,q , Ĉ g,2 =CV g,q , and Ĉ g,2=CV   g,q , the equation (26) may be changed into: 
 
 Ĉ   g,2 ( I   qs   −σÂ   g,2 ) −1   {circumflex over (R)}   g,2   =C ( I   n −σ( A −Δ)) −1   R  
 
         [0087]     The transfer matrix H Δ (S 0 +σ) of perturbation system defined in the equation (25) is actually equal to the transfer matrix Ĥ g,2 (s 0 +σ) of reduced system defined in the equation (22), namely Ĥ g,2 (s 0 +σ)=H Δ (s 0 +σ).  
         [0088]     Δ may be regarded as s of the rank of original matrix, and such perturbation shows the reason of the dynamic action of reduced system that is quite approximate to that of original system MNA. It meanwhile reflects the result of some restriction on the reduced system when circuits are interconnected. Since Δ comprises  
         h       q   +   1     ,   q     g     ,       μ   q     =            h     q   ,     q   -   1       g       h       q   +   1     ,   q     g                  
 
 may serves as the basis of GA iteration process termination. 
 
         [0089]     As shown in  FIG. 1 , a flow chart is used in this invention to describe the whole process of model reduction.  
         [0090]     At step  101 , a nodal analysis equation for an original circuit is inputted to create a circuit model equation (1).  
         [0091]     At step  102 , A=−(N+s 0 M) −1 M and R=(N+s 0 M) −1 B are set.  
         [0092]     At step  103 , an orthonormalization base K(A, R, q)=colsp[R, AR, . . . , A q−1 R]=colsp[V r,q ] in the Krylov subspace is found, and in the GA algorithm, a projection matrix V g,q  may be given; next, the projection is used in this invention to create a reduced system.  
         [0093]     At step  104 , a first method of reduction is used to find a reduced matrix at step  106 .  
         [0094]     At step  105 , a second reduced system is provided to find the reduced system provided at step  106  in a method of recursion.  
         [0095]     At step  107 , a residual error of the reduced system is deduced, and transfer functions of the reduced systems in a global Arnoldi algorithm and in a regional Arnoldi algorithm may be verified to be identical to each other.  
         [0096]     At step  108 , a math model for a perturbation system is deduced, and the transfer function of reduced matrix and that of original system additionally provided with the perturbation system is verified to be corresponding to each other.  
         [0097]     Again, a simple embodiment is used for test to verify the accuracy of algorithm according to this invention. In  FIG. 2 , an RLC circuit with 12 lines is provided. Line parameters are a resistor of 1.0 Ω/cm, a capacitor of 5.0 pF/cm, an inductor of 1.5 nH/cm, a drive resistor of 3Ω, and a load capacitor of 1.0 pF, respectively. Each line is 30 mm long and divided into 30 mm pieces. Thus, for the dimension of MNA matrix, n=1198 In the embodiment, a frequency range lower than {0,15 GH z } is used and the expanded frequency of a reduced system is {0,15 GH z }. When the GA algorithm starts, the values h q+1,q  and h q+1,q  are recorded.  
         [0098]     In  FIG. 5  showing the concluded results of simulation, it is observed that μ q  becomes small when 12 times of iteration reach, and thus 12 is the recommended order setting of the reduced system, thereby the order of reduced system being qs=24. In the method of model reduction provided in this invention, the relevant error of system union before q order(s) for the original system compared with the three reduced systems Ĥ g,1 (s), Ĥ g,2 (s), and Ĥ b,q (s) in  FIG. 4 , in which the system union of first reduced model in the Global Arnoldi algorithm is apparent, is accurate, compared with the reduced model generated in the conventional Block Arnoldi algorithm. However, the system union of the three reduced systems are identical to the original system. Also, H ij (s) indicates the impact of an input source i on a receiving end of output J.  FIG. 5  shows transfer matrices of original system H(s) and two reduced systems H(s) and H(s). In the figure, it is observed that the transfer matrix of reduced system and the transfer matrix of original system are corresponding at a frequency expansion point, and that the frequency response of two reduced systems are completely corresponding to each other, and thus frequency response curves fully overlap in the whole simulation frequency domain. Further,  FIG. 6  shows the frequency response H(s) of original system, the frequency response Ĥ g,2 (s) of second reduced system in the Global Arnoldi algorithm, and the frequency response H Δ (s) of original system additionally provided with the perturbation system, in which the frequency response Ĥ g,2 (s) and H Δ (s) is corresponding to the original system near the frequency expansion point; also, from the result of frequency response simulation, it is verified that the two are corresponding and the overlap of response curves proves that high accuracy achieves.  
         [0099]     This invention gives a method of reducing the model in the MIMO interconnect circuit system, in which the operation complexity of simulation and analysis on the interconnect circuit may be reduced in the global Arnoldi algorithm. From the two reduced model system according to this invention, it is verified that the system union of the preceding q orders is completely corresponding to the original system. It is verified that the output transfer function of the first reduced system is completely corresponding to the reduced system in the block Arnoldi algorithm. Besides, this invention provides the residual error that may serve as a reference material for determination of the reduced model order. Next, from the second reduced system according to this invention, the math model of perturbation system may be derived, in which it is verified that the output transfer function of second reduced system is corresponding to the transfer function of original system additionally provided with the perturbation system in height.  
         [0100]     While the invention has been described in terms of what is presently considered to be the most practical and preferred embodiments, it is to be understood that the invention needs not be limited to the disclosed embodiment. On the contrary, it is intended to cover various modifications and similar arrangements included within the spirit and scope of the appended claims which are to be accorded with the broadest interpretation so as to encompass all such modifications and similar structures.