Method of developing an analogical VLSI macro model in a global Arnoldi algorithm

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

1. Field of the Invention

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.

2. Description of Related Art

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.

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.

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.

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.

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.

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.

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.

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)⁻¹ of the original circuit, which is thus not practical to LSI.

In another prior art, the decomposed matrix (i_n−sA)⁻¹ 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).

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).

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

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.

For solution to the problems above, this invention provides a method of developing an analogical VLSI macro model in a global Arnoldi algorithm, comprising:

step 1 to input a net-shaped circuit;

step 2 to input a frequency expansion point;

step 3 to build up a state-space matrix for a circuit;

step 4 to generate a projection matrix by way of a global Arnoldi algorithm;

step 5 to determine a reduced model order by an iteration termination condition, and to execute a first model reduction; and

step 6 to build up a mathematics model for a perturbation system and to execute a second model reduction.

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,qV_g,q⁺−I_n)h_q+1,q^gV_g,q+1E_q^T(I_qs−s(H_g,qI_s)−sV_g,q⁺ΔV_g,q)∥E_r(s)∥_∞≦κ(S_qs)∥V_g,qV_g,q⁺−I_n∥₂|h_q+1,q^g|∥V_g,q⁺∥ ∥R∥₂ 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 ${\mu }_q=\frac{h_{q,q-1}^g}{h_{q+1,q}^g}<\varepsilon ,$ where ε is a permissible error that is small enough.

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\frac{d{x}_{\Delta }\left(t\right)}{dt}=x_{\Delta }\left(t\right)+Ru\left(t\right)\mathrm{and}y\left(t\right)=C{x}_{\Delta }\left(t\right)$

where Δ=h_q+1,q^gV_g,q+1E_q^TV_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.

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.

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

FIG. 1 is a flow chart of this invention.

FIG. 2 is a schematic view illustrating an embodiment of test on an interconnection circuit with 2 inputs and 2 outputs.

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.

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.

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.

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

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.

In a prior art, a linear, time-invariant, RLCG interconnect circuit in VLSI is generally represented in the following Modified Nodal Analysis (MNA) formula: $\begin{array}{cc}M\frac{dx\left(t\right)}{dt}+Nx\left(t\right)+Bu\left(t\right)=0\mathrm{and}y\left(t\right)=Cx\left(t\right).& \left(1\right)\end{array}$

where $M=\left[\begin{array}{cc}C& 0\\ 0& L\end{array}\right]$ comprises a capacitor and an inductor, and $N=\left[\begin{array}{cc}0& E\\ -E^T& R\end{array}\right]$ comprises a resistor and satisfies the incidence matrix E in Kirchhoff'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\left(t\right)=\left[\begin{array}{c}v\left(t\right)\\ i\left(t\right)\end{array}\right],$ while u(t) represents a function of system input signal. For A=−(N+s₀M)⁻¹M and R=(N+s₀M)⁻¹B, s₀ represents a frequency expansion point for selection, assuming that N+s₀M is a non-singular matrix. Then, equation (1) may be changed into $\begin{array}{cc}A\frac{dx\left(t\right)}{dt}=x\left(t\right)+Ru\left(t\right)\mathrm{and}y\left(t\right)=Cx\left(t\right)& \left(2\right)\end{array}$

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: $\begin{array}{cc}\hat{A}\frac{d\hat{x}\left(t\right)}{dt}=\hat{x}\left(t\right)+\hat{R}u\left(t\right)\mathrm{and}y\left(t\right)=\hat{C}\hat{x}\left(t\right)& \left(3\right)\end{array}$

where for  ε R^q×q and  ε R^q×q, q represents the order of reduced system. This invention is briefly described for s=k.

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−1V_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:Vqb, Hqb)
(1): /* Initialize */
Factor V0R=QR factorization of R
Vqb = [V1]
(2): /* Iteration */
for K = 1,2,···,q
Vk(0) = AVk−1
for j= 1,2,···,k
Hk−j,k−1b = Vk−jTVk(j−1)
Vk(j) = Vk(j−1) − vk−jHk−j,k−1b
end for
Factor (Vqb)k(Hqb)k,k−1 = QR factorization of Vk(k)
end for

In the process of BA algorithm iteration, the relation is: AV_b,q=V_b,qH_b,q+V_b,qh_q,q−1^bE_q^T   (4)

where $E_q={\left[\begin{array}{cccc}{0}_s& \cdots & 0_s& I^{s\times s}\end{array}\right]}^T\in R^{\mathrm{qs}\times s}\mathrm{and}$$H_{b,q}=\left[\begin{array}{ccccc}{h}_{11}^b& h_{12}^b& \cdots & \cdots & h_{1q}^b\\ h_{21}^b& h_{22}^b& \cdots & \cdots & h_{2q}^b\\ 0& h_{32}^b& h_{33}^b& \cdots & ⋮\\ ⋮& ⋰& ⋰& ⋰& ⋮\\ 0& \cdots & 0& h_{q,q-1}^b& h_{qq}^b\end{array}\right]\in R^{qs\times qs}$

H_b,q is an upper Hessenberg matrix, and the relation between H_b,q and A is: H_b,q=V_b,q^TAV_b,q   (5)

where V_b,q^T V_b,q=I^qs×qs.

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^jV_b,1=V_b,qH_b,q^jE_q, ∀_j=0,1 . . . ,q−1.   (6) where for E₁└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,qE₁.

For a non-singular matrix, the system union is defined to X^(j)(s₀)=A^jR=A^jV_b,1G=V_b,qH_b,q^jE₁G. For j=0, . . . , q−1, X^(j)(s₀) relates to the matrix V_b,q. Thus, the reduced system may be defined to: Â=V_b,q^TAV_b,q=H_b,q, {circumflex over (R)}=V_b,q^TR, Ĉ=CV_b,q   (7)

For j=0, . . . ,q−1 the system union of reduced system may be corresponding to that of original system: $\begin{array}{c}{\hat{Y}}^{\left(j\right)}\left(s_o\right)=\hat{C}{\hat{A}}^j\hat{R}\\ =C{V}_{b,q}{V}_{b,q}^T{A}^j{V}_{b,q}{V}_{b,q}^{T}R\\ =C{A}^{j}R=Y^{\left(j\right)}\left(s_o\right)\end{array}$

where for the property of applied matrix operation, if V^TV=I and V^TV=I then V^TV=I.

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⁢:⁢ ⁢Vqg,Hqg)
(1): /* Initialize */
V0=R/RF
Vqg=[ ]⁢ ⁢and⁢ ⁢Vqg=[V1]
(2): /* Iteration */
for j = 1, 2, . . . , q
V~=AVjg
for j = 1, 2, . . . , j
hi,j=〈Vi,V~〉F
V~=V~-Hi,jg⁢Vig
end for
Hj+1,jg=V~F
Vj+1g=V~hj+1,j
Vqg=[VqgVj+1g]
end for

where for V_g,1 matrix, $V_{g,1}=\frac{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−1R} in the Krylov subspace K_q(A,R)=span{R, AR, . . . , A^q−1R}, 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, <V_g,i, V_g,j>_F=0   (8) For i=j, <V_g,i, V_g,j>_F=1   (9) where <.,.>_F is a Frobenius inner product, namely <A,B>_F=trace(A^TB)=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₁,A₂, . . . ,A_k, is linearly independent and the matrix of R^mn, vec(A₁),vec(A₂), . . . , 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)^Tvec(B)=trace(A^TB).   (11) (AB)(CD)=(AC×BD).   (12) Further, a Frobenius norm is defined to ∥A,B∥_F=√{square root over (|trace(A^TB)|)}=√{square root over (vec(A)^Tvec(B))} (proposed in year 1985 by P. Lancaster and M Tismenetsky, The Theory of Matrices: with Applications, Academic Press).

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)^TV_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ⁿ may be given in the BA algorithm, which is a differentia from the GA and BA algorithms.

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,qI_s)+h_q+1,q^gV_g,q+1E_q^T,   (13) where $H_{g,q}=\left[\begin{array}{ccccc}{h}_{11}^g& h_{12}^g& \dots & \dots & h_{1q}^g\\ h_{21}^g& h_{22}^g& \dots & \dots & h_{2q}^g\\ 0& h_{32}^g& h_{33}^g& \dots & ⋮\\ ⋮& ⋰& ⋰& ⋰& ⋮\\ 0& \cdots & 0& h_{q,q-1}^g& h_{\mathrm{qq}}^g\end{array}\right]\in {•}^{q\times q}.,$

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ⁿε R^n×n, and A and B Kronecker is made to be A{circle around (×)} B ε R^mn×mn, being defined below to: $A\otimes B=\left[\begin{array}{cccc}{a}_{11}B& a_{12}B& \dots & a_{1m}B\\ a_{21}B& a_{22}B& \dots & a_{2m}B\\ ⋮& ⋮& ⋮& ⋮\\ a_{m1}B& a_{m2}B& \dots & a_{mm}B\end{array}\right]={\left[a_{\mathrm{ij}}B\right]}_{i,j=1}^m$

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.

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₀)=A^jR and vectorization, vec(A^jR)=(I_sA^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).

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).

The equation (13) may be further reduced to: A^jV_g,1=V_g,q(H_g,q^je₁I_s)   (14)

For V_g,1=V_g,qE₁ and h_q+1,q^gV_g,q+1E_q^TE₁=0, E₁ is multiplied at the right side of equation (13) and E₁ is derived. Assuming that the equation (14) comes into existence for, in case of j=i+1, $\begin{array}{c}{A}^{i+1}{V}_{g,1}={\mathrm{AV}}_{g,q}\left(H_{g,q}{e}_1\otimes I_s\right)\\ =\left(V_{g,q}\left(H_{g,q}\otimes I_s\right)+h_{q+1,q}^g{V}_{g,q+1}{E}_q^T\right)\left(H_{i,q}{e}_1\otimes I_s\right)\\ =V_{g,q}\left(H_{i+1,q}{e}_1\otimes I_s\right)+h_{q+1,q}^g{V}_{g,q+1}{E}_q^T\left(H_{i,q}\otimes I_s\right)E_1\end{array}$

For all the matrices Z and q≠1, then E_q^TZE₁=0 exists. A^jV_g,1=V_g,q(H_g,q^je₁I_s), ∀_j32 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.

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,qAV_g,q, {circumflex over (R)}_g,1=V_g,qR, and Ĉ_g,1=CV_g,q,   (16)

where V_g,q⁺is a virtual inverse matrix of V_g,q and defined to V_g,q⁺=(V_g,q^TV_g,q)⁻¹V_g,q^T,

where the reduced system  may be further reduced to: $\begin{array}{c}{\hat{A}}_{g,1}=V_{g,q}^{+}{\mathrm{AV}}_{g,q}\\ =V_{g,q}^{+}\left(V_{g,q}\left(H_{g,q}\otimes I_s\right)+h_{q+1,q}^g{V}_{g,q+1}{E}_q^T\right)\\ =\left(H_{g,q}\otimes I_s\right)+V_{g,q}^{+}\Delta V_{g,q}\end{array}$

where Δ is defined to Δ=h_q+1,q^gV_g,q+1E_q^TV_g,q⁺. Also, {circumflex over (R)} may be reduced to: {circumflex over (R)}_g,1=V_g,q⁺R=V_g,q⁺∥R∥_FV_g,q(e₁I_s)=∥R∥_F (e₁I_s)

By means of the reduction technology developed in the GA algorithm, the property of union may be achieved.

In case of j=0,1, . . . , q−1, the output union Ŷ^(j)(s₀) of reduced system that is given in the global Arnoldi algorithm is equal to the Y^(j)(s₀) union output of original MNA system in equation (2), namely $\begin{array}{c}{\hat{Y}}^{\left(j\right)}\left(s_0\right)={\stackrel{̑}{C}}_{g,1}{\hat{A}}_{g,1}^j{\hat{R}}_{g,1}\\ ={\mathrm{CV}}_{g,q}{V}_{g,q}^{+}{A}^j{V}_{g,q}{V}_{g,q}^{+}R\\ ={\mathrm{CA}}^{j}R\\ =Y^{\left(j\right)}\left(s_0\right)\end{array}$

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.

Thus, it is identical, proving that, if the system union is expressed as {X⁰(s₀),X¹(s₀), . . . X^q−1(s₀)}, 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: $\left[E_1{V}_{b,q}{E}_1{V}_{b,q}{H}_{b,q}{E}_1\dots V_{b,q}{H}_{b,q}^{q-1}\right]=V_{b,q}\left[E_1{H}_{b,q}{E}_1\dots {\left(H_{b,q}\right)}^{q-1}{E}_1\right]=V_{b,q}{\Re }^b$

H_b,q is an upper Hessenberg matrix, so is an triangle matrix on a block.

On the other hand, according to the equation (15), the system union given in the GA algorithm may be expressed as: $\left[V_{g,q}{E}_1{V}_{g,q}\left(H_{g,q}\otimes I_s\right)E_1\dots V_{g,q}\left(H_{g,g}^{q-1}\otimes I_s\right)E_1\right]=V_{g,q}\left[E_1\left(H_{g,q}\otimes I_s\right)E_1\dots \left({\left(H_{g,q}\right)}^{q-1}\otimes I_s\right)E_1\right]=V_{g,q}{\Re }^g$

where is also an upper triangle matrix, since $\mathrm{colsp}\left[X^0\left(s_0\right),X^1\left(s_0\right),\dots X^{q-1}\left(s_0\right)\right]=\mathrm{colsp}\left(V_{b,q}\right)=\mathrm{colsp}\left(V_{g,q}\right)·V_{b,q}=V_{g,q}{{\Re }^g\left({\Re }^b\right)}^{-1},{{\Re }^g\left({\Re }^b\right)}^{-1}$ is an triangle matrix and non-singular. Thus, it proves that the two transfer matrix are the same.

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)

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).

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.

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: $\begin{array}{cc}\begin{array}{c}{E}_r\left(s\right)=\left(I_n-\mathrm{sA}\right)V_{g,q}\hat{X}\left(s\right)-R\\ =\left(V_{g,q}-{\mathrm{sAV}}_{g,q}\right){\left(I_{\mathrm{qs}}-s\left(H_{g,q}\otimes I_s\right)\right)}^{-1}{V}_{g,q}^{+}R-R\\ =\left(\begin{array}{c}{V}_{g,q}\left(I_{\mathrm{qs}}-s\left(H_{g,q}\otimes I_s\right)\right)-\\ {\mathrm{sh}}_{q+1,q}^g{V}_{g,q+1}{E}_q^T\end{array}\right){\left(I_{\mathrm{qs}}-s\left(H_{g,q}\otimes I_s\right)\right)}^{-1}{\left(V_{g,q}\right)}^{+}R-R\end{array}& \left(18\right)\end{array}$

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^gV_g,q+1E_q^T(I_qs−s(H_g,qI_s))⁻¹(V_g,q)⁺R   (19)

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,qI_s)=S_qsΛ_qsS_qs⁻¹ is the property value resolution of (H_g,qI_s), the equation (19) being reduced to: $\begin{array}{cc}\begin{array}{c}{E}_r\left(s\right)=-{\mathrm{sh}}_{q+1,q}^g{V}_{g,q+1}{E}_q^T{S_{\mathrm{qs}}\left(I_{\mathrm{qs}}-s{\Lambda }_q\right)}^{-1}{S}_{\mathrm{qs}}^{-1}{V}_{g,q}^{+}R\\ ={\mathrm{sh}}_{q+1,q}^g{V}_{g,q+1}{E}_q^T{S}_{\mathrm{qs}}Z\left(s\right)S_{\mathrm{qs}}^{-1}{V}_{g,q}^{+}R\end{array}& \left(20\right)\end{array}$

where $Z\left(s\right)={\mathrm{diag}\left[\frac{s}{1-s{\lambda }_i}\right]}^{\mathrm{qs}}.,$ since Z(s) is a high-pass matrix, ${Z\left(s\right)}_{\infty }={\min }_{i=1}^{\mathrm{qs}}\frac{1}{{\lambda }_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∥₂   (21)

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 ${\mu }_q=\frac{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.

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)

where Δ=h_q+1,q^gV_g,q+1E_q^TV_g,q. From the equation (12), it is known that the reduced  may be further reduced to: Â_g,2=H_qI_s

where in case of i=0,1, . . . ,q−1, the union matching property still exists. $\begin{array}{cc}\begin{array}{c}{\stackrel{̑}{C}}_{g,2}{\stackrel{̑}{A}}_{g,2}^i{\stackrel{̑}{R}}_{g,2}={{\mathrm{CV}}_{g,q}\left(H_{g,q}\otimes I_s\right)}^i{R}_F\left(e_1\otimes I_s\right)\\ ={\mathrm{CV}}_{g,q}\left(H_{g,q}^i{e}_1\otimes I_s\right){R}_F\\ ={\mathrm{CA}}^i{V}_{g,1}{R}_F\\ ={\mathrm{CA}}^{i}R\end{array}& \left(23\right)\end{array}$

Thus, the transfer matrix of reduced system may be changed into: Ĥ_g,2(s₀+σ)=Ĉ_g,2(I_q−σÂ_g,2)⁻¹{circumflex over (R)}_g,2   (24)

Next, an equation of the status of a perturbation system that is defined in this invention is: $\begin{array}{cc}A\frac{d{x}_{\Delta }\left(t\right)}{dt}=x_{\Delta }\left(t\right)+\mathrm{Ru}\left(t\right)\mathrm{and}y\left(t\right)={\mathrm{Cx}}_{\Delta }\left(t\right)& \left(25\right)\end{array}$

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₀+σ)=H_Δ(s₀+σ).

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)⁻¹=(I−σ(A−Δ))⁻¹V_g,q

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−σÂ)⁻¹∥R∥_F(e₁I_s)=C(I_n−σ(A−Δ))⁻¹V_g,q∥R∥_F(e₁I_s)   (26)

Since Ĉ_g,2=CV_g,q, Ĉ_g,2=CV_g,q, and Ĉ_g,2=CVg,q, the equation (26) may be changed into: Ĉ_g,2(I_qs−σÂ_g,2)⁻¹{circumflex over (R)}_g,2=C(I_n−σ(A−Δ))⁻¹R

The transfer matrix H_Δ(S₀+σ) of perturbation system defined in the equation (25) is actually equal to the transfer matrix Ĥ_g,2(s₀+σ) of reduced system defined in the equation (22), namely Ĥ_g,2(s₀+σ)=H_Δ(s₀+σ).

Δ 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,{\mu }_q=\frac{h_{q,q-1}^g}{h_{q+1,q}^g}$ may serves as the basis of GA iteration process termination.

As shown in FIG. 1, a flow chart is used in this invention to describe the whole process of model reduction.

At step 101, a nodal analysis equation for an original circuit is inputted to create a circuit model equation (1).

At step 102, A=−(N+s₀M)⁻¹M and R=(N+s₀M)⁻¹B are set.

At step 103, an orthonormalization base K(A, R, q)=colsp[R, AR, . . . , A^q−1R]=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.

At step 104, a first method of reduction is used to find a reduced matrix at step 106.

At step 105, a second reduced system is provided to find the reduced system provided at step 106 in a method of recursion.

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.

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.

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.

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.

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.

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.

Claims

1. Method of developing an analogical VLSI macro model in a global Arnoldi algorithm, comprising: step 1 to input a net-shaped circuit; step 2 to input a frequency expansion point; step 3 to build up a state-space matrix for a circuit; step 4 to generate a projection matrix by way of a global Arnoldi algorithm; step 5 to determine a reduced model order by an iteration termination condition, and to execute a first model reduction; and step 6 to build up a mathematics model for a perturbation system and to execute a second model reduction.

2. The method of developing an analogical VLSI macro model in a global Arnoldi algorithm according to claim 1, wherein from the iteration termination condition described at step 5, a residual error for a first reduced system and an original system may be defined to:

3. The method of developing an analogical VLSI macro model in a global Arnoldi algorithm according to claim 1, wherein 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 nodal analysis may be indicated below as: