IEICE TRANS. FUNDAMENTALS, VOL.E88–A, NO.2 FEBRUARY 2005

513

PAPER

Special Section on Analog Circuit Techniques and Related Topics

A Reduction Technique for RLCG Interconnects Using Least
Squares Method
Junji KAWATA†a) , Yuichi TANJI†† , Yoshifumi NISHIO††† , Members, and Akio USHIDA† , Fellow

SUMMARY
In this paper, we propose a new algorithm for calculating
the exact poles of the admittance matrix of RLCG interconnects. Afterchoosing dominant poles and corresponding residues, each element of theexact admittance matrix is approximated by partial fraction. A procedure
to obtain the residues that guarantee the passivity is also provided, based on
experimental studies. In the procedure the residues are calculated by using
the least squares method so that the partial fraction matches each elementof the exact admittance matrix in the frequency-domain. From the partial
fraction representation, the asymptotic equivalent circuit models which canbe easily simulated with SPICE are synthesized. It is shown that an eﬃcientmodel-order reduction is possible for short-length interconnects.
key words: RLCG interconnects, Leverrier-Faddeeva algorithm, leastsquares method, passivity, asymptotic equivalent circuits, reduction

1.

Introduction

The analysis and design of high speed LSI chips are becoming more and more important, because the sub-circuits
coupled with interconnects embedded in the substrate sometimes cause the fault switching operations due to the signal
delays, crosstalks, reflections and so on [1]–[5]. The Elmore
resistance-capacitance (RC) delay metric is popular due toits simple closed-form expression, computation speed andfidelity with respect to the simulation [4]. The closed-form
combining with the delay and crosstalk is firstly presented
in the reference [5]. The modified algorithms are proposed
later for the improvement of the accuracy and the practicalapplications in the simulations [6]–[9].
AWE (asymptotic waveform evaluation method) [10] iswidely used as a reduction technique of the large scale of linear networks, the algorithm is based on a moment-matching
technique and Padé approximation. However, the method
sometimes become erroneous, if there exist the poles far
from the origin. To overcome the problem, CFH (complex
frequency hopping) [11] and multi-point Padé approximation [12] methods are proposed. In these reduction algorithms, the reduced circuits sometimes become unstable in
the time domain even if all the poles are located in the left
hand side of the complex plane. The ill-condition can be
Manuscript received June 20, 2004.
Manuscript revised September 15, 2004.
Final manuscript received October 20, 2004.
†
The authors are with the Faculty of Engineering, Tokushima
Bunri University, Sanuki-shi, 769-2193 Japan.
††
The author is with the Faculty of Engineering, Kagawa University, Takamatsu-shi, 761-0396 Japan.
†††
The author is with the Faculty of Engineering, Tokushima
University, Tokushima-shi, 770-8506 Japan.
a) E-mail: kawata@fe.bunri-u.ac.jp

overcome by PVL (Padé via Lanczos process) [13], and
PRIMA (passive reduced-order interconnect macromodeling algorithm) [2], [14]. In order to apply these algorithms
to the interconnects, we need two steps such that each interconnect is firstly modeled by a finite order system, andArnoldi-based congruence transformation is applied to thesystem to form its reduced order model.
In this paper, we consider LSIs such as ASIC or SoC(System on a Chip) are coupled with interconnects embedded in the substrate. In this case, the diﬀusion resistancecomponents of the interconnects are generally assumed tobe very large compared to those of PCBs [3] and the lengths
are very short. From the telegrapher’s equation of the interconnects, the admittance matrix can be derived from therelations at the near and far ends [1]. We propose here a
new computational algorithm for calculating the exact polesin the complex plane. Next, each of the elements of theadmittance matrix is approximated by partial fraction using
dominant poles around the origin, where the corresponding
residues are evaluated by the least squares method. A computational procedure to obtain the residues that guarantee
the passivity is also provided. Thus, the admittance matrix
are reduced as the partial fractions and are subsequently synthesized as asymptotic equivalent circuits.
We present the numerical methods for calculating theexact poles and corresponding residues in Sect. 2. Section
2.3 shows the asymptotic equivalent circuits of the interconnects. Illustrative examples are given in Sect. 3. First,
distribution of poles for some cases is investigated. Then,
using the residues estimated by the least squares method,
we calculate the frequency response curves and comparedthem with the exact ones. Here, it will be investigated that
how the passivity of partial fractions obtained from poleresidue pairs breaks and why inaccurate residues are obtained. Based on the considerations, we provide a procedure
to guarantee the passivity of the asymptotic circuits. Furthermore, SPICE simulations using the asymptotic equivalent circuit and the built-in lossy transmission line modelof SPICE are carried out and the transient waveforms arecompared. It is found that the transient responses obtainedby the macromodels are in good agreement with those of
SPICE model as the number of poles increases, and that theorder of partial fractions can be eﬃciently reduced for shortlength interconnects.

c 2005 The Institute of Electronics, Information and Communication Engineers
Copyright IEICE TRANS. FUNDAMENTALS, VOL.E88–A, NO.2 FEBRUARY 2005

514

2.

at the locations satisfying sinh γi (s) = 0, i = 1, 2, . . . , N.
Thus, we have the following theorem for calculation of theexact poles.

Calculation of Exact Poles of Interconnects

2.1 Admittance Matrix of Interconnects and Its Poles
Now, consider a uniform N coupled RLCG interconnects
described by the following telegrapher’s equations:dV(x, s)= −(R + sL)I(x, s) 



dx
(1)



dI(x, s)= −(G + sC)V(x, s)dx
where, V(x, s) ∈ RN and I(x, s) ∈ RN are the voltage andcurrent vectors at x, and R ∈ RN×N , L ∈ RN×N ,C ∈ RN×N
and G ∈ RN×N are the resistance, inductance, capacitance and conductance matrices per unit length, respectively.
These are non-negative definite symmetric matrices. From(1), we have
dV2 (x, s)= (R + sL)(G + sC)V(x, s) 



dx
(2)


2

dI (x, s)= (G + sC)(R + sL)I(x, s) dx
Let us introduce the transform matrices Pv (s) and Pc (s) inorder to transform (2) into diagonal forms. Thus, we havediag[γi (s)2 ] = Pv (s)−1 (R + sL)(G + sC)Pv (s)(3)
diag[γi (s)2 ] = Pc (s)−1 (G + sC)(R + sL)Pc (s)When a complex frequency s is given, the eigenvalues γi (s)2and the corresponding eigenvectors can be uniquely obtained. Because Pv (s) and Pc (s) have the eigenvectors as thecolumns, these matrices can be also uniquely determined forthe complex frequency s. Pc (s) satisfies the following relations:
Pc (s) = (R + sL)−1 Pv (s)Γ(s)PTc (s) = Pv (s)−1

(4)
(5)

for Γ(s) = diag[γi (s)]. Then, the input and output relations
at the near and far ends are described by the admittance matrix as follows [1]:

 


I(0, s)Y11 (s) Y12 (s)V(0, s)=(6)
I(d, s)Y21 (s) Y22 (s) V(d, s)
where d is the length of the interconnects, and {V(0, s),V(d, s)} and {I(0, s), I(d, s)} denote the terminal voltage andcurrent vectors of the interconnects, respectively. The submatrices Yi j (s) ∈ RN×N of the admittance matrix are givenby:
Y11 (s) = Y22 (s)



−1


= Pc (s)diag[coth γi (s)d]Pv (s)Y12 (s) = Y21 (s) 
(7)






1


Pv (s)−1 = −Pc (s)diagsinh γi (s)dObserve that all the poles of admittance matrix can be foundTheorem 1: The locations of poles satisfying the relations(7) are found by solving the following equation:nπ 2 
I = 0, n = 1, 2, . . .
(8)
(R + sL)(G + sC) +dProof: In the case of n  0, we have from (7) that the poles
satisfy the following relations:|Pv (s)diag[sinh γi (s)d]Pc (s)−1 | = 0|Pv (s)diag[tanh γi (s)d]Pc (s)−1 | = 0

(9)
(10)

Since the transform matrices Pv (s) and Pc (s) are nonsingular
for the nonzero eigenvalues, the poles satisfying the above
two relations are given byγi (s)d = jnπ, i = 1, 2, . . . , N, n = 1, 2, . . .
Therefore, the characteristic equation from the telegrapher’sequation (1) has to satisfy the relation (8).
Q.E.D.
Corollary 1.1: The poles at n = 0 satisfy the following
relation:|R + sL| = 0(11)

Proof: We have γi = 0 at n = 0. Thus, using the relationof (4), we have the following relation:
1
lim Pc (s)diag
Pv (s)−1γi →0sinh γi (s)d1
Pv (s)−1= (R + sL)−1 Pv (s)diag[γi (s)]diag[γi (s)d]1= (R + sL)−1(12)d
Q.E.D.
2.2 Numerical Methods for Calculating Exact Poles
Generally, it is not easy to calculate the poles using the relation (8). Hence, we apply Leverrier-Faddeeva algorithm[16] in order to obtain the characteristic equation |sI−A| = 0
for the relation (8). Setting (8) into


nπ 2
−1
2
I+ s LC(13)
A(s) = −(RC+LG) RG+dwe have|sI−A(s)| = α0 (s)+ sα1 (s)+· · ·+ sN−1 αN−1 (s)+ sN(14)where αk (s) is found by the following recurrence formulae:BN−1 = I,⇒ αN−1 (s) = −tr(A(s))





BN−2 (s) = A(s)BN−1 + αN−1 (s)I,




1


⇒ αN−2 (s) = − tr(BN−2 (s)A(s)) 
(15)

2




............................................





B0 (s) = A(s)B1 (s) + α1 (s)I,





1



⇒ α0 (s) = − tr(B0 (s)A(s))
NKAWATA et al.: A REDUCTION TECHNIQUE FOR RLCG INTERCONNECTS USING LEAST SQUARES METHOD

515

Note that the maximum degree of (14) eventually becomes
2N, because αk (s) is polynomial with s. This algebraic
equation can be numerically solved by the use of Bairstow
method and so on. In this way, we can calculate the exact
poles of the admittance matrix (7).
We choose dominant poles located around the origin ofcomplex plane, because the poles have large influence on
the transient responses. Thus, an arbitrary (i, j) element of
the submatrices Y11 (s), Y12 (s), Y21 (s), and Y22 (s) of (7) isapproximated by the partial fractions that is composed ofthe dominant poles and the corresponding residues, in thefollowing forms† :Fig. 1 Equivalent circuit models of each term in the partial fractions (16)and (17).

Ŷ11 (s) = Ŷ22 (s) =
N

N

M

k01,i
b11,i,n s + b10,i,n+2s − p0,i i=1 n=1 s + a1,i,n s + a0,i,n
i=1

(16)

Ŷ12 (s) = Ŷ21 (s) =
N

N

M

k02,i
b21,i,n s + b20,i,n+2s − p0,i i=1 n=1 s + a1,i,n s + a0,i,n
i=1

(17)

where N is the number of conductors and M shows the number of complex conjugate pairs for each γi (s).
It has been shown in [18] that the absolute values ofresidues are the same, but their signs change alternately withn; i.e.,k02,i = −k01,ib21,i,n = (−1)n+1 b11,i,n 
(18)



b
= (−1)n+1 b
20,i,n10,i,nTherefore, there is no need to evaluate the residues for (17).
Equations (16) and (17) are divided into two parts depending on whether n is odd or even, where the first term in theequations corresponds to the case of n = 0 and is included in
the even part. Since all elements of the approximated submatrices Ŷ11 , Ŷ12 , Ŷ21 , and Ŷ22 are separated into the odd
and even parts, we can express the submatrices with the two
new submatrices Ya and Yb as follows.
 Ŷ11 (s) Ŷ12 (s)Ya (s)+Yb (s) Ya (s)−Yb (s)=Ya (s)−Yb (s) Ya (s)+Yb (s)Ŷ21 (s) Ŷ22 (s)(19)
where Ya and Yb are corresponding to the odd and evenparts, respectively. The asymptotic equivalent circuit models to be synthesized must satisfy the relations (16)–(19).
2.3 Asymptotic Equivalent Circuit of Interconnects
In this subsection, we derive the asymptotic equivalent circuit models that satisfy the relations (16)–(19). Each termof partial fractions (16) and (17) is synthesized by the circuit models shown in Figs. 1(a) and (b), where the circuit
parameters are given by
Li,0 =

p0,i1, Ri,0 =
k0,i
k0,i

(20)

Fig. 2 Equivalent circuit models of the matrix (19) for single-conductor
interconnect. (a) Ya and (b) Yb .








b1,i,n








Ci,n =2
a0,i,n b1,i,n + (b0,i,n − a1,i,n b1,i,n )b0,i,n 




2

b1,i,n b0,i,nGi,n =
2
a0,i,n b1,i,n + (b0,i,n − a1,i,n b1,i,n )b0,i,n 
Li,n =

1, Ri,n =a1,i,n b1,i,n − b0,i,nb21,i,n3
b1,i,n(21)

Each element of the admittance matrix (19) is modeled byparallel connection of the sub-circuits shown in Fig. 1. For
example, in the case of single-conductor interconnect, wehave the circuits shown in Fig. 2 for Ya (s) and Yb (s), where
the subscript i is omitted for simplicity. In the case of singleconductor interconnect, we can construct two equivalentcircuit models in the diﬀerent forms shown in Figs. 3 and4(a), where Ir (0, s) and Ir (d, s) in Fig. 4(a) are the voltagecontrolled current sources given by Figs. 4(b) and 4(c). Observe that the circuit of Fig. 3 does not contain any controlled source, so that the circuit is passive if all of thecircuit parameters have positive values. Also, it should be
noted that the circuit shown in Fig. 4(a) is passive if thecircuit shown in Fig. 3 is passive. The macromodels can†
To avoid confusion about suﬃx notation, the suﬃx (i, j) isomitted in (16)–(17).

IEICE TRANS. FUNDAMENTALS, VOL.E88–A, NO.2 FEBRUARY 2005

516

Fig. 3 Equivalent circuit model for single-conductor interconnect usinga bridge circuit.

Fig. 6Equivalent circuit model to multi-conductor interconnect.

2.4 Numerical Procedure for Calculating Residues
The residues satisfying (16) can be calculated by the leastsquares method. The frequency response curve of each element of the admittance matrix (7) must be matched to the exact one, namely Y(sk ) = Ŷ(sk ) at sampled frequency points
sk , k = 1, . . . , K. Thus we have
AX = B(22)

where X is unknown vector containing the residues, andB is column vector consisted of exact sampled frequencydomain data:
Fig. 4 Equivalent circuit models for single-conductor interconnect. The
subsystems Ya and Yb are shown in Fig. 2.

X = k0,1 · · · k0,N b1,1,1 b0,1,1 · · · b1,N,M b0,N,M

T(23)
B = Y(s1 ) Y(s2 ) · · · Y(sK )

T

(24)

and each element of A is given as:
Ak, j =

1,sk − p0, j

Ak,N+2 j−1 =Ak,N+2 j =j = 1, . . . , N

(25)


sk



2

sk + a1, j sk + a0, j 


(26)

1

s2k + a1, j sk + a0, j

,







j = (i − 1)M + n, i = 1, . . . , N, n = 1, . . . , M
Fig. 5Equivalent circuit model to two-conductors interconnect.

be easily implemented on any kind of SPICE-like simulators, since the asymptotic circuits are composed of typicalcircuit elements only. The model of two conductors interconnect is realized by the use of 4 current transformers asshown in Fig. 5. The asymptotic equivalent circuit for multiconductor interconnect is shown in Fig. 6, which is easilymodeled with SPICE.

The least squares solution X can be found by using QR factorization so that ||B − AX||2 is minimized.
In order to guarantee the passivity of the macromodels, the following trial-and-error approach for calculating
the residues is provided. In Sect. 3, it will be explained indetail how the following procedure ensures the passivity ofthe macromodels.
Step 1) Determine the highest frequency fmax considering
the frequency spectrum of input pulse signals.
Step 2) Choose all poles with the imaginary part less thanor equal to 2π fmax and further several extra poles beyond the frequency. Let the number of the former andthe latter be M and α, respectively. Thus the value ofKAWATA et al.: A REDUCTION TECHNIQUE FOR RLCG INTERCONNECTS USING LEAST SQUARES METHOD

517

M in (23) and (26) is replaced with Mext = M + α. The
appropriate initial value of α is about 5 ∼ 10.
Step 3) Perform the least squares fitting and subsequently
truncate the extra α pole-residue pairs beyond 2π fmax .
Step 4) If all of the residues (23) and the resulting circuit
parameters (20) and (21) are positive, finish the computation; otherwise increase the number of extra poles
α and return to Step 3.
It should be noted that the proposed method enforces
the passivity of the interconnect macromodels. All the values of the passive elements included in the macromodels
are constrained to be positive. Hence, the macromodels areguaranteed to be passive.
3.

Illustrative Examples

3.1 Distribution of PolesFirst, we investigate distribution of poles in the complex
plane by the computational method described in the subsection 2.2, because this gives us very important information tochoose dominant poles. The parameters of interconnects are
set as follows:
Ri,i = r [Ω/mm], other elements of R are 0.
Li,i = 10 [nH/mm], Li,i+1 = Li−1,i = 1 [nH/mm],
Li,i+2 = Li−2,i = 0.1 [nH/mm], other elements of L are 0.
Ci,i = 4 [pF/mm], Ci,i+1 = Ci−1,i = −0.21 [pF/mm],
Ci,i+2 = Ci−2,i = −0.01 [pF/mm], other elements of C are 0.

Gi,i = 0.5 [mS/mm], Gi,i+1 = Gi−1,i = −0.05 [mS/mm],
Gi,i+2 = Gi−2,i = −0.01 [mS/mm], other elements of G are 0.
and the values of N, r and d (resp., the number of conductors, the resistance component and the length of interconnect) are changed in the examples.
Figure 7 shows distribution of poles for interconnectswith diﬀerent values of N. The poles with negative imaginary part are not shown in these figures. Observe that they
consist of N real poles and complex conjugate poles for eachn in (8). As shown in Fig. 7 the imaginary parts of the complex conjugate poles increase monotonically as n in (8). We
assume an input pulse waveform which has large spectrum
in the lower frequency range. Thus model accuracy in thelower frequency range, which is closely related to poles located near the origin, is very important for transient simulations. Since poles for large n are far from the origin, they do
not play important roles for the transient responses.
Next, either the length d or the resistance component
r of 4-conductors interconnect is changed and the distributions are examined. The results are shown in Fig. 8. Compared Fig. 7(c) with Fig. 8(a), it is obvious that the imaginary parts of the poles for short-length interconnect become
quickly large for n. Further, from Figs. 7(c) and 8(b), we
found that the real parts of poles become larger, as the resistance components of the interconnect increase.
Since some of poles keep away from the origin do notplay important roles for the transient response, it is expectedthat for either short-length or high-resistive interconnects
the number of dominant poles will become small, i.e., ef-(a) N = 1(b) N = 2(c) N = 4
Fig. 7 Distribution of poles for N-conductor interconnects. The poles with negative imaginary partare not shown in this figure. (r = 0.5 [Ω/mm], d = 5.0 [mm])(a) r = 0.5 [Ω/mm], d = 1.0 [mm](b) r = 25.0 [Ω/mm], d = 5.0 [mm]Fig. 8 Distribution of poles for 4-conductor interconnect for the cases that (a) the length d becomesshort and (b) the resistance component r becomes large.

IEICE TRANS. FUNDAMENTALS, VOL.E88–A, NO.2 FEBRUARY 2005

518(a) |Y11 ( jω)|(b) |Y12 ( jω)|Fig. 9 Frequency response curves for single conductor interconnect.
5.0 [mm], fmax = 6.0 [GHz], M = 10)

fective reduction may be possible.
3.2 Calculation of Residues by Least Squares Method
We mention here the reason why the residues calculation
procedure described in subsection 2.4 is required, throughthe following numerical results. Throughout our simulations, we set the real part of the complex frequency sk aszero (i.e., sk = jωk = j2π fk ) in (24) and (25), and the sampling interval is taken to be 0.01 [GHz] equally.
Case 1: Let the highest frequency and the number of poles
be fmax = 6 [GHz] and M = 10, respectively. In this case thefrequency corresponding to the pole p10 (the suﬃx i = 1 isomitted) with largest imaginary part is about 5 [GHz]. The
frequency response curve obtained by substituting the calculated poles and residues to the partial fraction (16) or (17)and the exact curve obtained from (7) are shown in Fig. 9.
Here, the above-mentioned procedure of residues calculation is not used, and |Ŷ12 ( jω)| is calculated using the relation (18). The approximated curves by the partial fractions
have relatively good agreement with the exact one up to 11th peak around 5 [GHz], because we took a real pole and
10 of complex conjugate poles (i.e., M = 10). From thefigure, it seems that the least squares fitting is successfully
performed. However the calculated residues contain negative values as listed in the left-half of Table 1, which yields
negative circuit parameters in (21). Thus, it is clear that thepassivity is not satisfied in this case.
Case 2: For the case that the resistance component of interconnect becomes large, the least squares fitting fails toapproximate the frequency-domain response accurately, as
shown in Fig. 10. In this example, |Ŷ12 ( jω)| is also calculated using the relation (18). From these two examples, we
can conclude that the pole selection in this manner is notsuitable for our purpose.
Case 3: Next we tried to choose some extra poles beyond
the highest frequency fmax . After choosing 6 of extra poles
beyond fmax = 6 [GHz] (i.e., M = 18), the residues arecalculated by performing the least squares fitting, which
are listed in the right-half of Table 1 and the frequency response curves are shown in Fig. 11. As expected, |Ŷ11 ( jω)|

Table 1
3.
n
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

(r = 0.5 [Ω/mm], d =Residues obtained by the least squares method for Cases 1 andCase 1b1,nb0,n3.815e−02
3.979e−03
3.815e−02
1.939e−03
3.816e−02−1.568e−03
3.818e−02−6.727e−03
3.820e−02−1.385e−02
3.822e−02−2.349e−02
3.825e−02−3.664e−02
3.828e−02−5.537e−02
3.829e−02−8.508e−02
3.818e−02−1.483e−01

Case 3b1,nb0,n4.000e−02
5.002e−03
4.000e−02
5.008e−03
4.000e−02
5.017e−03
4.000e−02
5.016e−03
4.000e−02
5.004e−03
4.000e−02
4.979e−03
4.000e−02
4.940e−03
4.000e−02
4.960e−03
4.000e−02
5.049e−03
4.000e−02
5.127e−03
4.000e−02
4.623e−03
3.990e−02
1.383e−02
4.149e−02
3.551e+00
1.704e+00−1.005e+02−1.902e+01
6.750e+026.859e+01−1.764e+03−9.761e+01
1.972e+03
4.815e+01−7.910e+02

matches the exact one within the whole frequency range(0 ≤ f ≤ fmax ) of interest, as shown in Fig. 11(a). Unfortunately, as shown in Fig. 11(b), |Ŷ12 ( jω)| obtained from(17) and (18) is entirely diﬀerent from the exact one. Furthermore, the values of |Ŷ11 ( jω)| at certain frequency range
beyond the highest matched frequency ( fmax = 6 [GHz]) become very large and diﬀer completely from the exact one,as shown in Fig. 11(c). Note that the only diﬀerence ofFigs. 11(a) and (c) is the scales of axes.
From the right-half of Table 1, the residues for n >12, which correspond to the extra poles, contain negative or
large values and consequently do not guarantee the passivity.
Thus, it is concluded that disagreement of |Ŷ12 ( jω)| and hugevalues of |Ŷ11 ( jω)| at frequency range beyond fmax are dueto the extra pole-residue pairs.
Case 4: In this example, to remove the pole-residue pairs
beyond the highest frequency is considered. In the same
conditions as Case 3, 6 extra pole-residue pairs for n > 12
are truncated and the frequency response curves are computed. The results are illustrated in Fig. 12. The eﬀect oftruncation begins to appear at the last half of the frequency
range; that is, the accuracy of partial fractions (16) and (17)at higher frequency range is degraded. However, at lower

KAWATA et al.: A REDUCTION TECHNIQUE FOR RLCG INTERCONNECTS USING LEAST SQUARES METHOD

519(a) |Y11 ( jω)|(b) |Y12 ( jω)|Fig. 10 Frequency response curves for single conductor interconnect. (r = 25.0 [Ω/mm], d =5.0 [mm], fmax = 6.0 [GHz], M = 10)(a) |Y11 ( jω)|(b) |Y12 ( jω)|(a) |Y11 ( jω)|(b) |Y12 ( jω)|Fig. 12 Frequency response curves for the case that 6 poles beyond fmax
are truncated. (r = 25.0 [Ω/mm], d = 5.0 [mm], fmax = 6.0 [GHz], M =12, α = 6)

ber of extra poles, even if negative residues and/or circuit
elements appeared around fmax .

(c) |Y11 ( jω)|Fig. 11 Frequency response curves for the case that extra poles beyond
fmax are used. (r = 25.0 [Ω/mm], d = 5.0 [mm], fmax = 6.0 [GHz], M =18)frequency range both curves remain matched. Furthermore
the calculated response |Ŷ12 ( jω)| using the relation (18) is
also matched to the exact one at lower frequency range,
unlike the Case 3. Note that the partial fractions (16) and(17) (or the resulting asymptotic equivalent circuits) become
more precise as the value of fmax increases, because matchedfrequency range becomes wide accordingly. Further it has
been confirmed that positive residues and/or circuit elementsof the macromodels can be obtained by increasing the num-

In order to investigate the eﬀect of the truncation andhow to determine the number of extra poles α, we examinedroot mean square (rms) error between the exact and the approximated frequency responses in 0 ≤ f ≤ fmax , with different line parameters, highest matched frequency fmax andnumber of extra poles α. The results are shown in Table 2,
where the most right column indicates whether all passive
elements have positive values or not when the extra poleresidue pairs are truncated. Note that the rms error of Y11
increases due to truncating the extra poles, while the one
for Y12 always decreases. This coincides with the results
shown in Figs. 11 and 12. When the same line parametersand diﬀerent α are taken, the diﬀerences between the rms
errors are small. Therefore, a small α is preferable fromcomputational eﬃciency point of view. However, this does
not necessarily obtain the passive macromodels as shownin Table 2. Roughly speaking, α depends on the number

IEICE TRANS. FUNDAMENTALS, VOL.E88–A, NO.2 FEBRUARY 2005

520
Table 2 Root mean square error between the approximated and exact responses, where “truncated” and “not truncated” meanthat the terms corresponding to the extra poles are truncated or not, and “True” and “False” in the column of passivity indicatethat the obtained circuit is passive or not, respectively.
line
parameters

fmax [GHz]
(M)
5.0
(10)

R = 25.0 [Ω/mm],
L = 100.0 [nH/mm],
C = 0.4 [pF/mm],
G = 0.5 [mS/mm],
d = 5.0 [mm]

10.0
(20)

15.0
(30)

5.0
(10)
R = 0.5 [Ω/mm],
L = 1.0 [nH/mm],
C = 40.0 [pF/mm],
G = 0.5 [mS/mm],
d = 5.0 [mm]

10.0
(20)

15.0
(30)

5.0
(2)
R = 25.0 [Ω/mm],
L = 10.0 [nH/mm],
C = 4.0 [pF/mm],
G = 0.5 [mS/mm],
d = 1.0 [mm]

10.0
(4)

15.0
(6)

5.0
(2)
R = 25.0 [Ω/mm],
L = 1.0 [nH/mm],
C = 40.0 [pF/mm],
G = 0.5 [mS/mm],
d = 1.0 [mm]

10.0
(4)

15.0
(6)

α
3
5
10
20
3
5
10
20
3
5
10
20
3
5
10
20
3
5
10
20
3
5
10
20
3
5
10
20
3
5
10
20
3
5
10
20

rms error for Y11
not truncatedtruncated
2.911e−06
9.570e−04
5.968e−08
9.552e−04
4.420e−119.560e−04
3.305e−119.560e−04
9.816e−06
1.025e−03
4.938e−07
1.021e−03
8.583e−11
1.025e−03
1.129e−10
1.025e−03
1.767e−05
1.052e−03
1.397e−06
1.044e−03
8.882e−101.056e−03
5.192e−11
1.055e−03
3.699e−04
9.635e−02
4.611e−06
9.575e−02
1.568e−08
9.574e−02
6.777e−08
9.574e−02
1.591e−03
1.042e−01
7.420e−05
1.028e−01
6.906e−08
1.027e−01
6.899e−08
1.027e−01
3.031e−03
1.079e−012.410e−04
1.059e−01
1.301e−07
1.057e−01
5.662e−08
1.057e−01
3.271e−07
6.790e−037.379e−106.789e−03
1.965e−106.789e−037.529e+05
6.789e−033.736e−06
8.188e−03
1.032e−08
8.180e−03
4.771e−108.180e−03
9.821e−108.180e−03
1.283e−05
8.874e−038.669e−08
8.854e−03
6.770e−108.854e−03
1.590e−108.854e−03

rms error for Y12
not truncatedtruncated
1.053e−02
1.090e−04
7.438e−02
1.061e−04
3.847e+02
1.052e−04
1.177e+05
1.051e−04
1.945e−02
1.135e−04
2.205e−01
8.557e−05
7.327e+01
7.906e−05
2.945e+05
7.907e−05
2.623e−02
1.435e−04
3.887e−01
8.303e−05
2.090e+02
6.601e−05
2.298e+05
6.607e−05
8.652e−01
9.866e−03
6.237e+00
1.066e−02
1.023e+05
1.068e−02
1.850e+08
1.068e−02
1.565e+00
6.394e−03
1.810e+01
7.926e−038.539e+03
8.040e−03
3.270e+08
8.040e−03
2.085e+00
4.963e−03
3.150e+016.459e−03
1.821e+04
6.719e−03
1.566e+08
6.719e−03
1.513e−02
1.618e−03
4.215e−02
1.619e−03
1.083e+06
1.619e−03
1.596e+06
1.619e−03
3.310e−02
1.419e−03
1.368e−011.429e−03
5.440e+05
1.429e−03
9.321e+05
1.429e−03
5.124e−02
1.245e−032.655e−01
1.270e−032.420e+04
1.270e−03
4.734e+06
1.270e−03

3
5
10
20
3
5
10
20
3
5
10
20

1.911e−06
1.072e−09
1.005e−132.510e+02
1.493e−05
4.768e−08
6.061e−12
3.881e−12
4.096e−05
3.089e−07
5.741e−125.656e−129.103e−01
3.702e+00
9.716e+02
1.898e+02
1.519e+00
7.608e+00
1.021e+05
6.500e+04
2.115e+00
1.279e+01
1.006e+045.911e+04

of poles which have the imaginary part less than or equal
to 2π fmax . With increase of the number of poles, M and α
should be taken largely in order to obtain the passive macromodels. However, α is also dependent on the line parameters. Therefore, it is diﬃcult to determine a suitable valueof α in advance, and the proposed procedure provided in
Subsect. 2.4 is necessary.
3.3 Transient Simulations
Now, let us calculate the transient responses for a simple

6.813e−02
6.473e−02
6.470e−02
6.470e−02
7.983e−02
7.866e−02
7.843e−02
7.843e−028.544e−02
8.611e−028.539e−028.539e−021.220e−02
1.424e−02
1.424e−02
1.424e−02
1.249e−02
1.211e−021.230e−021.230e−02
1.821e−02
1.012e−02
1.086e−02
1.086e−02

Passivity
False
True
True
True
False
False
True
True
False
False
True
True
False
True
True
True
False
True
True
True
False
False
True
True
True
True
True
True
True
True
True
True
True
True
True
True
False
True
True
True
False
True
True
True
False
False
True
True

linear circuit with a single-conductor interconnect shown in
Fig. 13. The SPICE simulations using the asymptotic equivalent circuit of Fig. 4(a) are carried out with varying thehighest matched frequency fmax (accordingly M is changed),
and compared with the simulated results of built-in SPICE
model. Figure 15 shows the transient response waveforms
for r = 0.5 [Ω/mm] and d = 5.0 [mm]. All the waveforms have relatively good agreement, although ringing occurs at rising and falling edges of waveforms. The ringing can be suppressed with increasing M, but remarkable

KAWATA et al.: A REDUCTION TECHNIQUE FOR RLCG INTERCONNECTS USING LEAST SQUARES METHOD

521

Fig. 13 Simple linear circuit with single conductor interconnect. (RS =10[Ω], RL = 10[Ω])Fig. 14 Full-adder circuit. The parameters of all interconnects T 1 ∼ T 4
are R = 500 [Ω/mm], L = 10 [µH/mm], C = 0.01 [pF/mm], G =0.5 [µS/mm], d = 5 [mm].

changes are not observed even for M = 100 (in this casefmax = 50 [GHz]).
Figure 16 is the transient waveforms for the case thatthe resistance component r of interconnect was increased.
In this case, the poles are far from the imaginary axis (seeFig. 8 (b)), but this does not aﬀect the transient waveforms;
that is, the number of poles M can not be reduced for thecase of interconnect with high resistance component.
Next, after setting the length d of interconnect to be
short, the transient simulation is carried out. In this case, thepoles go far away from the real axis with n of (8), which isillustrated in Fig. 8(a). The transient waveforms for M = 10
(in this case fmax = 25 [GHz]) are approximately coincidewith the SPICE results, as shown in Fig. 17. Thus, for thecase of short-length interconnect, eﬀective reduction is possible.
Finally, as an example of circuit with nonlinear element, a full-adder circuit shown in Fig. 14 is simulated and(a) Near end(b) Far end
Fig. 15 Transient response of the circuit shown in Fig. 13. (r = 0.5 [Ω/mm], d = 5.0 [mm])(a) Near end(b) Far end
Fig. 16 Transient response of the circuit shown in Fig. 13. (r = 25.0 [Ω/mm], d = 5.0 [mm])(a) Near end(b) Far end
Fig. 17 Transient response of the circuit shown in Fig. 13. (r = 25.0 [Ω/mm], d = 1.0 [mm])IEICE TRANS. FUNDAMENTALS, VOL.E88–A, NO.2 FEBRUARY 2005

522(a) vS
Fig. 18(b) vC
Transient response for the full-adder circuit shown in Fig. 14.

the results are illustrated in Fig. 18. Replacing the interconnect with the corresponding asymptotic equivalent circuit,
where M = 20, eﬀective simulation using SPICE is possible
for nonlinear cases.
4.

Conclusion

In this paper, we have proposed an algorithm for calculating
the exact poles, the corresponding residues are estimated by
the least squares method. A procedure for calculation of theresidues that guarantee the passivity has also provided. Using the poles and residues, each element of the admittancematrix of interconnects is approximated by the partial fraction. Though the frequency response curves calculated using the partial fractions becomes inaccurate around the highest matched frequency fmax , the partial fractions and the resulting asymptotic equivalent circuit models are always passive. Further, it has been shown, from the investigationsof pole distributions and the transient simulations, that theeﬀective model-order reduction is possible for the case of
short-length interconnects. Thus our method is eﬃciently
applied for VLSIs.
References[1] J.A. Brandao Faria, Multiconductor Transmission-Line Structures:
Modal Analysis Techniques, John Wiley & Sons, 1993.
[2] M. Celik, L. Pileggi, and A. Odabasioglu, IC Interconnect Analysis,
Kluwer Academic Pub., 2002.
[3] J.L. Wyatt, Jr., “Signal propagation delay in RC models for interconnect,” Circuit Analysis, Simulation and Design, ed. A.E. Ruehli,
pp.254–291, Elsevier Science Pub. B.V., North-Holland, 1987.
[4] W.C. Elmore, “The transient response of damped linear networks
with particular regard to wideband amplifiers,” J. Appl. Phys.,vol.19, pp.55–63, Jan. 1948.
[5] T. Sakurai, “Closed-form expressions for interconnection delay, coupling, and crosstalk in VLSI’s,” IEEE Trans. Electron Devices,
vol.40, no.1, pp.118–124, 1993.
[6] S.-Y. Kim, N. Gopal, and L.T. Pillage, “Time-domain macromodels
for VLSI interconnect analysis,” IEEE Trans. Comput.-Aided Des.
Integr. Circuits Syst., vol.13, no.10, pp.1257–1270, 1994.
[7] R. Gupta, B. Tutuianu, and L.T. Pileggi, “The Elmore delay as abound for RC trees with generalized input signals,” IEEE Trans.
Comput.-Aided Des. Integr. Circuits Syst., vol.16, no.1, pp.95–104,
1997.
[8] C.J. Alpert, A. Devgan, and C.V. Kashap, “RC delay metrics for performance optimization,” IEEE Trans. Comput.-Aided Des. Integr.
Circuits Syst., vol.20, no.5, pp.571–582, 2001.

[9] T. Kimura and M. Okumura, “An eﬃcient reduction method of a
substrate RC network model,” IEICE Trans. Fundamentals, vol.E84A, no.3, pp.698–704, March 2001.
[10] E. Chiprout and M.S. Nakhla, Asymptotic Waveform Evaluationand Moment Matching for Interconnect Analysis, Kluwer Academic
Pub., 1994.
[11] E. Chiprout and M.S. Nakhla, “Analysis of interconnect networks
using complex frequency hopping (CFH),” IEEE Trans. Comput.Aided Des. Integr. Circuits Syst., vol.14, no.2, pp.186–200, 1995.
[12] Q. Yu, J. Meiling, and E.S. Kuh, “Passive multipoint moment matching model order reduction algorithm on multiport distributed interconnect networks,” IEEE Trans. Circuits Syst. I, vol.46, no.1,
pp.140–160, 1999.
[13] P. Feldmann and R.W. Freund, “Eﬃcient linear analysis by Padé approximation via the Lanczos process,” IEEE Trans. Comput.-Aided
Des. Integr. Circuits Syst., vol.14, no.5, pp.639–649, 1995.
[14] A. Odabasioglu, M. Celik, and L.T. Pileggi, “PRIMA: Passive
reduced-order interconnect macromodeling algorithm,” IEEE Trans.
Comput.-Aided Des. Integr. Circuits Syst., vol.17, no.8, pp.645–
654, 1998.
[15] A. Ushida, K. Urabe, Y. Yamagami, and Y. Nishio, “Asymptotic
equivalent circuits of interconnects based on complex frequency
method,” Proc. ECCTD’01, vol.3, pp.29–32, 2001.
[16] D.K. Faddeeva and V.N. Faddeeva, Computational Methods of Linear Algebra, W.H. Freeman and Co., 1963.
[17] J. Kawata, Y. Tanji, Y. Nishio, and A. Ushida, “Asymptotic equivalent circuits of interconnects and the passivity,” IEICE Technical
Report, NLP2003-134, 2003.
[18] Y. Yamagami, Y. Tanji, A. Hattori, Y. Nishio, and A. Ushida, “A
reduction technique of large scale RCG interconnects in the complex
frequency domain,” Int. J. Circuit Theory Appl., (to appear).
[19] A. Yonemoto, T. Hisakado, and K. Okumura, “An improvementof convergence of FFT-based numerical inversion of Laplace transforms,” Proc. ISCAS’02, vol.V, pp.769–772, 2002.
[20] W.T. Beyene and J.E. Schutt-Aine, “Eﬃcient transient simulationof high-speed interconnects characterized by sampled data,” IEEE
Trans. Compon. Packag. Manuf. Technol. B, vol.21, no.1, pp.105–
114, Feb. 1998.
[21] M. Celik and A.C. Cangellaris, “Eﬃcient transient simulation of
lossy packaging interconnects using moment-matching techniques,”
IEEE Trans. Compon. Packag. Manuf. Technol. B, vol.19, no.1,
pp.64–73, Feb. 1996.
[22] K.L. Choi, N. Na, and M. Swaminathan, “Characterization of embedded passives using macromodels in LTTC technology,” IEEE
Trans. Compon. Packag. Manuf. Technol. B, vol.21, no.3, pp.258–
268, Aug. 1998.
[23] T. Watanabe and H. Asai, “Synthesis of time-domain models forinterconnects having 3-D structure based on FDTD method,” IEEE
Trans. Circuits Syst. II, vol.47, no.4, pp.302–305, April 2000.
[24] Y. Tanji and M. Tanaka, “Hierarchical least-squares algorithm for
macromodeling high-speed interconnects characterized by sampled
data,” IEICE Trans. Fundamentals, vol.E83-A, no.9, pp.1833–1843,

KAWATA et al.: A REDUCTION TECHNIQUE FOR RLCG INTERCONNECTS USING LEAST SQUARES METHOD

523

Sept. 2000.
[25] M. Suzuki, H. Miyashita, A. Kamo, T. Watanabe, and H. Asai, “A
synthesis of time-domain interconnect models by MIMO type of selective orthogonal least-square method,” IEEE Trans. Microw. Theory Tech., vol.49, no.10, pp.1708–1714, Oct. 2001.

Junji Kawata
received the B.E., M.E.,
and Ph.D. degrees from Tokushima University,
Tokushima, Japan, in 1990, 1992 and 2001, respectively. Since 1992, he has been with the
Department of Mechanical and Electronic Engineering at Tokushima Bunri University, where
he is currently a Lecturer. His research interests
are in circuit simulation and chaotic communication. Dr. Kawata is a member of the IEEE.

Yuichi Tanji
received the B.E., M.E.,
and Ph.D. degrees from Tokushima University,
Tokushima, Japan, in 1993, 1995 and 1998, respectively. After graduation, he held Research
Associate positions at Sophia University and
Kagawa University, and is currently an Associate Professor of the Department of Reliabilitybased Information Systems Engineering at Kagawa University. His research interests arein circuit simulation, artificial neural networks,
and image and signal processing. Dr. Tanji is amember of the IEEE.

Yoshifumi Nishio
received the B.E., M.E.,
and Ph.D. degrees in electrical engineering from
Keio University, Yokohama Japan, in 1988,
1990, and 1993, respectively. In 1993, he joined
the Department of Electrical and Electronic Engineering at Tokushima University, Tokushima
Japan, where he is currently an Associate Professor. From May 2000 he spent a year in the
Laboratory of Nonlinear Systems (LANOS) at
the Swiss Federal Institute of Technology Lausanne (EPFL) as a Visiting Professor. His research interests include analysis and application of chaos in electrical circuits, analysis of synchronization in coupled oscillators circuits, development of analyzing methods for nonlinear circuits and theory and application of cellular neural networks. He is currently the Chair of the IEEE CAS
Technical Committee on Nonlinear Circuits and Systems (NCAS) and theSecretary of the IEICE Research Society of Nonlinear Theory and its Applications (NOLTA). He is serving as an Associate Editor for the IEEE
Transactions on Circuits and Systems - Part I, for the area of Theory and
Systems. He is a member of the IEEE and the RISP.

Akio Ushidareceived the B.E. and M.E.
degrees in electrical engineering from Tokushima University in 1961 and 1966, respectively, and the Ph.D. degree in electrical engineering from University of Osaka Prefecture in
1974. He was an associate professor from 1973
to 1980 at Tokushima University. From 1980 to
2003, he was a Professor in the Department ofElectrical Engineering at Tokushima University.
Since 2003, he has been with the Department of
Mechanical and Electronic Engineering, Tokushima Bunri University. From 1974 to 1975 he spent one year as a visiting
scholar at the Department of Electrical Engineering and Computer Sciencesat the University of California, Berkeley. His current research interests include numerical methods and computer-aided analysis of nonlinear system.
Dr. Ushida is a member of the IEEE.